Review of Methods to Solve Desiccation Cracks in Clayey Soils

: This paper reviews numerical methods used to simulate desiccation cracks in clayey soils. It examines ﬁve numerical approaches: Finite Element (FEM), Lattice Boltzmann (LBM), Discrete Element (DEM), Cellular Automaton (CAM), and Phase Field (PFM) Methods. The paper presents a simpliﬁed description of the methods, including their basic numerical formulations. Several factors such as the multiphase nature of soils, heterogeneity, nonlinearities, coupling, scales of analysis, and computational aspects are discussed. The review highlights the characteristics, strengths, and limitations of each method. FEM shows a good capacity to deal with the thermo-hydromechanical behavior of clays when drying that complement well with the ability of DEM to deal with particle interactions as well as LBM, PFM, and CAM to deal with complex crack patterns. The article concludes by reviewing the integration of multiple numerical methods to enhance the simulation of desiccation cracks in clayey soils and proposing what is the best option to continue improving the study of this problem.

Simulating desiccation cracks in clayey soils is a complex task due to several reasons. Firstly, clayey soils are multiphase media composed of soil particles, and pores that contain only air, only water, or water and air, depending on if the condition is dry, saturated, or unsaturated. Secondly, clayey soils exhibit coupled nonlinear THM behavior. As moisture content decreases, the soil undergoes significant volume changes due to the suction generated in the soil matrix, resulting in shrinkage first and cracking when the soil strength is reached. This nonlinear behavior requires advanced constitutive models and numerical techniques to accurately capture the soil's response to environmental contour conditions [32]. The coupling of multiple physical processes, including fluid flow, heat transfer, and deformation is accounted for by incorporating temperature-dependent and moisture-dependent properties. Additionally, constitutive relationships are employed to couple moisture content and mechanical behavior. Thermal expansion coefficients and thermal conductivity should be considered as functions of mechanical strain/stress to account for the coupling between temperature and deformation. Thirdly, the behavior of desiccation cracks is influenced by various factors such as soil composition, mineralogy, pore structure, and initial moisture content. The inherent variability and uncertainties associated with these factors make it difficult to predict crack formation and propagation accurately [33][34][35]. Fourthly, simulating desiccation cracks often requires considering different scenarios, from Geotechnics 2023, 3 809 laboratory specimens to field-scale applications. Accounting for scale effects and capturing the heterogeneity of soil properties is crucial for realistic simulations [36]. Finally, simulating desiccation cracks in clayey soils requires computationally intensive simulations due to the need to solve complex nonlinear coupled equations and handle large deformation and long-term drying processes.
This problem involves several interconnected physical processes in a portion of a multiphase soil system that is in mechanical equilibrium. The shrinkage that takes place as moisture is lost is governed by the principle of balance, known as Richards' equation, and the generalized Darcy's law. At the same time, there will be a distribution of temperature governed by the heat balance equation and Fourier's law. Mechanical deformation arises from shrinkage, exhibiting elastic and plastic behavior described by the THM constitutive model [30].
In the pores of the soil, there are physical processes that occur during the desiccation and cracking. The air dissolution in water is governed by Henry's law, and the diffusion of air in water is governed by Fick's law, with water molecules moving from areas of high moisture content to low moisture content. Additionally, heat transfer occurs under effects such as Soret's thermal diffusion of water vapors in the air because of pressure gradients produced by temperature gradients [37][38][39], vapor effusion, and Stefan's flow [40]. Fortunately, in many cases, they can be neglected and continue capturing the main mechanisms that govern the desiccation and cracking processes. Thus, they are not necessarily needed in a formulation and a numerical model.
At the beginning of the desiccation, the soil is a saturated fluid slurry, but with time, the condition turns to compacted unsaturated soil. For this reason, the degree of saturation must be included in the formulation and simulation by using, for example, a simplified Van Genuchten's formula [41] or more complex variations that include the effect of temperature. To accurately model the mechanical behavior of clayey soils during desiccation, THM constitutive equations are necessary. These equations define the relationship between stress and strain and capture the material's response when suction and temperature change. The simplest stress-strain relationship is the generalized isotropic linear elastic model commonly known as Hooke's law, which characterizes stress-strain behavior in the linear elastic range. Even if the soil's mechanical behavior is considered elastic, the equation must include the effect of the temperature and suction to couple the thermal, hydraulic, and mechanical processes. More sophisticated models are nonlinear, viscoelastic, or plasticity, such as bilinear, Mohr-Coulomb models, and many others.
To simulate the initiation and propagation of cracks, Griffith's criterion (tensile strength controls the initiation of the cracks) or linear elastic fracture mechanics (LEFM) equations are commonly used due to their simplicity. LEFM principles and its rules determine the critical crack length and assess crack propagation. Incorporating LEFM principles allows for the analysis of crack formation and growth during desiccation. The numerical simulation of crack propagation is in particular a very challenging problem in the context of the Finite Element Method (FEM). For this reason, several approaches to simulate the cracking process have been proposed apart from the FEM. Lattice Boltzmann (LBM), Discrete Element (DEM), Cellular Automaton (CAM), and Phase Field (PFM) methods are all alternatives to the FEM that can be used to effectively simulate the desiccation cracks in clayey soils and cracks in other problems.
These numerical techniques enable the solution of complex equations and the simulation of desiccation cracks in clayey soils. Boundary conditions, initial conditions, and appropriate numerical algorithms also play a crucial role in accurately capturing the behavior of desiccation cracks. Figure 1 shows a laboratory test made on a cylindrical sample to study the problem of desiccation cracks in clayey soils under controlled conditions. A whole cycle of drying, wetting, flooding, drying, and cracking demonstrated that flooding produces more cracks and wetting modifies suction profiles. Even when this problem is usually studied as a desiccation problem, the first semi-cycle in Figure 1, wetting and flooding are part of the problem and significantly affect the cracking process. Today, the research community is working mainly on semi-cycles of desiccation. To fully understand this problem, the whole cycle must be understood.
Geotechnics 2023, 3, FOR PEER REVIEW 3 flooding are part of the problem and significantly affect the cracking process. Today, the research community is working mainly on semi-cycles of desiccation. To fully understand this problem, the whole cycle must be understood. In Section 2, physical and non-physical methods are reviewed, then, in Section 3, the integration of these methods is reviewed and commented on. Finally, in Section 4, the combination of FEM and CAM is presented as a promising alternative to tackle desiccation cracks in clayey soils.

Methods to Simulate Desiccation Cracks in Clayey Soils
In the attempt of simulating desiccation cracks in soils, researchers have used physicalbased approaches and non-physical-based approaches. In this section, five of the most effective methods to resolve this problem are commented on in terms of the main characteristics, strengths, and limitations. The strengths of these methods are presented in Table 1. Table 1. Methods that effectively tackle challenges when simulating desiccation cracks in clayey soils, and the scale levels they work into. The methods are classified into non-physical-based (nPb) and physical-based (Pb) methods into columns for every scale level. In Section 2, physical and non-physical methods are reviewed, then, in Section 3, the integration of these methods is reviewed and commented on. Finally, in Section 4, the combination of FEM and CAM is presented as a promising alternative to tackle desiccation cracks in clayey soils.

Methods to Simulate Desiccation Cracks in Clayey Soils
In the attempt of simulating desiccation cracks in soils, researchers have used physicalbased approaches and non-physical-based approaches. In this section, five of the most effective methods to resolve this problem are commented on in terms of the main characteristics, strengths, and limitations. The strengths of these methods are presented in Table 1.

The Finite Element Method (FEM)
FEM [23,24,26,32,34] and its variant Extended Finite Element Method (XFEM) [42] have been extensively used for simulating desiccation cracks in clayey soils at a macroscale level and can also be applied at micro-and mesoscale levels. Researchers have employed FEM to study moisture diffusion, shrinkage, and cracking behavior during drying.
FEM resolves classic transient continuum mechanics' partial differential equations that describe the phenomenon. FEM has been successful in capturing the complex behavior of desiccation cracks since it can deal with complex geometries and heterogeneity. It can map the distribution of stress in the soil mass locating the areas of concentration of stresses that produce cracks. FEM deals well with coupling the thermo-hydromechanical physical processes that the problem includes. The accuracy of this method relies on the appropriate implementation of constitutive models and boundary conditions.
During the desiccation process, the three phases of the soil interact in general thermally, hydraulically, and mechanically. Once the contour conditions in suction, temperature, and displacements are set, and if the soil-structure and soil-atmosphere interactions are neglected, the main equations that define the THM problem of desiccation in clayey soils are the governing equations and constitutive equations from continuum mechanics, respectively. Table 1. Methods that effectively tackle challenges when simulating desiccation cracks in clayey soils, and the scale levels they work into. The methods are classified into non-physical-based (nPb) and physical-based (Pb) methods into columns for every scale level. The limitations FEM has shown are mesh dependency, making it challenging to capture intricate crack patterns. Additionally, it has shown difficulties in accurately predicting crack propagation without explicit crack geometry modeling. FEM and the other methods presented here neglect soil-structure and soil-atmosphere interaction (Figure 2), which can significantly influence crack formation and behavior. Finally, calibrating constitutive models to accurately represent soil behavior is always a challenge when using FEM and any other numerical approximation. These limitations drive the need for specialized techniques and alternative numerical methods to overcome these challenges and improve the accuracy of simulations.  [42] have been extensively used for simulating desiccation cracks in clayey soils at a macroscale level and can also be applied at micro-and mesoscale levels. Researchers have employed FEM to study moisture diffusion, shrinkage, and cracking behavior during drying.
FEM resolves classic transient continuum mechanics' partial differential equations that describe the phenomenon. FEM has been successful in capturing the complex behavior of desiccation cracks since it can deal with complex geometries and heterogeneity. It can map the distribution of stress in the soil mass locating the areas of concentration of stresses that produce cracks. FEM deals well with coupling the thermo-hydromechanical physical processes that the problem includes. The accuracy of this method relies on the appropriate implementation of constitutive models and boundary conditions.
The limitations FEM has shown are mesh dependency, making it challenging to capture intricate crack patterns. Additionally, it has shown difficulties in accurately predicting crack propagation without explicit crack geometry modeling. FEM and the other methods presented here neglect soil-structure and soil-atmosphere interaction (Figure 2), which can significantly influence crack formation and behavior. Finally, calibrating constitutive models to accurately represent soil behavior is always a challenge when using FEM and any other numerical approximation. These limitations drive the need for specialized techniques and alternative numerical methods to overcome these challenges and improve the accuracy of simulations. During the desiccation process, the three phases of the soil interact in general thermally, hydraulically, and mechanically. Once the contour conditions in suction, temperature, and displacements are set, and if the soil-structure and soil-atmosphere interactions are neglected, the main equations that define the THM problem of desiccation in clayey soils are the governing equations and constitutive equations from continuum mechanics, respectively.
In the next section, the formulation for the desiccation cracks problem using FEM is presented. All the details of this formulation, including numerical examples, can be found In the next section, the formulation for the desiccation cracks problem using FEM is presented. All the details of this formulation, including numerical examples, can be found in [32].

Equilibrium Equation (Cauchy Equation of Motion)
If no dynamic effects are considered, the equilibrium equation of the soil matrix is as follows: ∇·σ Equation (1) is an elliptic partial differential equation where σ is the total stress tensor, and ρ is the average density of the multiphase medium (soil, water, or air). The vector g is the gravity vector.

Balance Equation (Continuity Equation also Known as Richards' Equation)
Equation (2) is a parabolic partial differential equation that represents the balance of water in the pores of the soil. In an unsaturated porous medium (the general case that includes the saturated case, S r = 1, and the dry case, S r = 0), the water mass balance equation is written as follows: In Equation (2), ρ w is the water density, q is Darcy's velocity vector, t is time, n is the porosity of the soil, and S r is the degree of saturation of water in the soil pores.

Conservation of Energy Equation (First Law of Thermodynamics)
If thermal effects are considered, the first law of thermodynamics establishes the need for the heat transfer equation in the soil. Equation (3) is a parabolic partial differential equation: In Equation (3), ρ s is the density of the soil, c is the specific heat capacity, θ is the temperature, q θ is the heat transfer rate, and K θ is the thermal conductivity (could be scalar for isotropic permeability or tensorial for anisotropic permeability).

Stress-Strain Thermos-Mechanical Constitutive Law
For the most general case of unsaturated soils, the effective stress tensor (σ ' ) is: In Equation (4), σ is the total stress tensor. The air and water pressure are, respectively, u a and u w , χ is a parameter that depends on the degree of saturation, the stress history and the soil's fabric and 1 ≡ δ ij are the identity tensor.
In this formulation, the matrix suction and the net mean stress define the effective stress tensor σ ' through Equation (4): The net stress σ net and the suction s are: The general strain-stress relation must be written in differential form, because of the nonlinearity of the material behavior.
For the most general THM case, D is a tangent matrix in function of the strain, ε, temperature, θ, and suction, s. Equation (7) establishes the coupling between temperature, suction, and mechanical effects.
The deformations are calculated by the addition of a component due to the net stress plus a component due to the suction plus a component due to temperature. Equation (8) considers, then, the additive deformation hypothesis: In Equation (8), the parameter K is the volumetric modulus and G is the shear modulus of the soil matrix; the parameter K s is the volumetric modulus due to changes in suction; the parameter K θ is the volumetric modulus due to temperature changes. These parameters must be established depending on the constitutive model chosen (linear elastic, non-linear elastic, viscoelastic, plastic, etc.); the factor C is a fourth-order compliance tensor and h, t are second-order tensors.
The net stress increments can be obtained from (8): In Equation (9), D = C −1 is the tangent stiffness tensor.

Generalized Darcy's Law for Unsaturated Soils and Permeability Tensor
The generalized Darcy's law for unsaturated soils is: In Equation (10), q is the velocity of Darcy's vector; ∇s is the gradient of the suction; K(S r , n, θ) is a permeability tensor that changes with water saturation degree, S r , porosity, n, and temperature, θ; g is the gravity vector and ρ w is the water density. The permeability tensor, K, can be isotropic or anisotropic.

Water Retention Curve
The van Genuchten function [13] is usually adopted in this formulation to relate changes between the degree of saturation and the suction, s: where, λ is a material parameter and P 0 is the air entry value for the initial porosity n 0 , adopted as the reference value. Function f n considers the changes of porosity during desiccation and its effect in the water retention curve using a parameter, η. For nondeformable soils, f n = 1, because porosity is constant.

Fourier's Law
This is the constitutive law for the thermal problem.
In Equation (12), q θ is the heat transfer rate and K θ is the thermal conductivity.

Hydro-Mechanical Formulation to Resolve Desiccation Cracks in Clayey Soils
As it was stated in previous sections, desiccation cracks in clayey soils are in general a THM problem. However, since the experimental programs are made under controlled conditions, it is usually chosen to work under isothermal conditions to simplify the study of the process [32,34]. Under these assumptions, the numerical simulations can be HM since the temperature will remain constant during the whole process. Under HM conditions, the problem is resolved by FEM resulting in the system of Equation (13) in what is known as a u-p formulation. This assumption simplifies the problem considerably since Equations (3) and (12) are not needed and Equations (7)-(9) are less complex.
In Equation (13), u are the displacements and p is the suction. The factors in Equation (13) are: Stiffness Matrix : Coupling Matrix : Vector of Nodal Flow :

Lattice Boltzmann Method (LBM)
LBM was introduced to simulate fluid flow for the first time in 1986 [43]. It was used to study fluid flow and fluid-solid interaction in [44,45]. LBM was used for dealing with several other problems like melting and solidification [46], gas transport into rocks [47], hydraulics fracturing [48], multiphase flow through micro-and macropores [49], fracture and flow [50], and finally, drying [51] and desiccation cracks [52].
LBM enables the consideration of pore-scale processes during drying. The fluid domain is discretized into a lattice structure, with each lattice node representing a small volume or pore. Instead of solving the governing equations at a continuum level, LBM simulates the behavior of individual fluid particles, represented by lattice cells or lattice Boltzmann particles, that move and interact within the lattice. By explicitly representing the individual fluid particles and their interactions, LBM allows for the consideration of various pore-scale processes during dryings, such as capillary effects, evaporation, fluid-solid interactions, and convective flows.
LBM employs a simplified kinetic model to describe the motion of fluid particles. These particles propagate along discrete lattice directions and undergo collisions with neighboring particles, leading to the redistribution of mass, momentum, and energy. The particle interactions at the pore scale directly influence the macroscopic behavior of the fluid.
LBM provides a means to couple pore-scale simulations with larger-scale models, such as FEM, to capture the interactions between the microscopic and macroscopic phenomena. This enables the integration of pore-scale information into continuum-based simulations and improves the accuracy of predictions at larger scales.
LBM provides insights into the fundamental mechanisms of crack formation, but its computational cost can be relatively high due to the need for fine spatial resolution. LBM can solve physical equations of balance and equilibrium. LBM is based on the Boltzmann equation that is derived from the principles of the conservation of mass, momentum, and energy, hence, a physical-based model.
The limitations of the LBM are that the accurate representation of soil behavior in LBM requires proper material characterization, which can be challenging due to the complexity of clayey soil slurry behavior during desiccation. Modeling crack propagation may need additional techniques, as the inherent lattice structure of LBM may not directly capture the process. Additionally, LBM primarily focuses on fluid flow, potentially overlooking some soil mechanics aspects when the soil acquires consistency, such as the mechanical behavior of the soil matrix, which influences crack initiation and propagation.

Formulation of LBM
LBM can be considered as a variant of the Finite Difference Method [53] since it generates a grid (lattice) where the fluid flow is simulated. However, its roots are statistical rather than deterministic [54].
The method calculates the probability of fictitious particles, that compose the material being modeled, being located at a particular point on a lattice as a function of time. It is executed by performing two algorithms. Algorithm 1 models the streaming of particles traveling from one lattice node to connected neighbors. Algorithm 2 mimics their collision ( Figure 3). The response of the particles to collision and the lattice used determine the model behavior.
phenomena. This enables the integration of pore-scale information into continuum-based simulations and improves the accuracy of predictions at larger scales.
LBM provides insights into the fundamental mechanisms of crack formation, but its computational cost can be relatively high due to the need for fine spatial resolution. LBM can solve physical equations of balance and equilibrium. LBM is based on the Boltzmann equation that is derived from the principles of the conservation of mass, momentum, and energy, hence, a physical-based model.
The limitations of the LBM are that the accurate representation of soil behavior in LBM requires proper material characterization, which can be challenging due to the complexity of clayey soil slurry behavior during desiccation. Modeling crack propagation may need additional techniques, as the inherent lattice structure of LBM may not directly capture the process. Additionally, LBM primarily focuses on fluid flow, potentially overlooking some soil mechanics aspects when the soil acquires consistency, such as the mechanical behavior of the soil matrix, which influences crack initiation and propagation.

Formulation of LBM
LBM can be considered as a variant of the Finite Difference Method [53] since it generates a grid (lattice) where the fluid flow is simulated. However, its roots are statistical rather than deterministic [54].
The method calculates the probability of fictitious particles, that compose the material being modeled, being located at a particular point on a lattice as a function of time. It is executed by performing two algorithms. Algorithm 1 models the streaming of particles traveling from one lattice node to connected neighbors. Algorithm 2 mimics their collision ( Figure 3). The response of the particles to collision and the lattice used determine the model behavior.

A Multi-Material Multiple-Relaxation-Time LBM
LBM models the statistical mechanics of many particles under collision [49,50]. In the formulation presented in [50], the objective of the method is to simulate the evolution of a single-particle-distribution function of statistical mechanics, ( ⃗, ⃗ , ). Where, ⃗ is the position of the particle, ⃗ is the velocity of the particle, and is time. The function, , is the probability of finding a particle at the position, ⃗, velocity, ⃗ , at time, .
In this method, the mass and momentum conservation equations for an isotropic material can be written as Equations (21) and (22):

A Multi-Material Multiple-Relaxation-Time LBM
LBM models the statistical mechanics of many particles under collision [49,50]. In the formulation presented in [50], the objective of the method is to simulate the evolution of a single-particle-distribution function of statistical mechanics, f p is the position of the particle, → ξ is the velocity of the particle, and t is time. The function, f , is the probability of finding a particle at the position, → p , velocity, → ξ , at time, t. In this method, the mass and momentum conservation equations for an isotropic material can be written as Equations (21) and (22): where, ρ is the ghost density; → j is the momentum in the three directions of the space; J B and J D are the bulk and deviatoric components of the momentum tensor; P is the mean stress and τ is the deviatoric stress. ν and ζ bulk and shear coefficients of viscosity.
In three dimensions, the stress is defined as in Equations (23) and (24): In LBM, the connection between the probability distribution function, f , of particles and the conservation equations (mass and momentum) is established through the streaming and collision steps.
Streaming step: After collision, the particles' f are streamed to neighboring lattice nodes according to their velocities. This step ensures the propagation of information through the fluid domain. By streaming particles, the macroscopic properties, such as density and velocity, are advected throughout the simulation domain.
Collision step: In this step, particles collide with each other locally, leading to the redistribution of particle velocities. The collision process obeys conservation laws, such as mass and momentum conservation. During collision, the equilibrium of f is determined based on local fluid properties, like density and velocity, at a given lattice point. The equilibrium of f represents the f that the system would tend toward if no collisions occurred.
The key point is that the macroscopic properties of the fluid, such as density and velocity, are related to the moments (summations) of the f . The zeroth moment corresponds to the particle density (mass conservation), and the first moment corresponds to the momentum density (momentum conservation). By tracking and updating the function f through the collision and streaming steps, the LBM indirectly solves the conservation equations.

Phase Field Method (PFM)
PFM appeared in 1992 and it was used as an effective tool for simulating desiccation cracks in clayey soils from micro-to mesoscale levels. PFM describes the evolution of a system with multiple phases and has been applied to represent the degree of saturation or water content in the soil [55][56][57][58].
PFM provides a continuous representation of crack formation and propagation, enabling the study of complex crack patterns. However, the computational cost associated with PFM is high. PFM is a mathematical framework that can be used to solve physical equations of balance and equilibrium, so, it is a physical-based method.
The limitations of PFM are the difficulty in accurately calibrating the model parameters to represent the specific behavior of clayey soils. The constitutive relations and material properties used in the PFM may need to be carefully tuned to capture the unique characteristics of clayey soils, such as their complex moisture retention and swelling-shrinkage behavior. Additionally, the PFM tends to smooth out crack features due to its diffuse interface representation of cracks, potentially overlooking small-scale crack details. The cracks change the contour conditions of the problem; since the method treats the cracks with continuous functions, the method cannot update the contour conditions.

General Phase Field for Fracture Model
In this section, the formulation for the desiccation problem from [55] is resumed to their basic elements. The necessary variational approach for cracks was proposed initially by the authors of [59], the phase field implementation is due to [60] and an overview can be found in [61]. In this method, crack surfaces are represented by a surface density function in terms of an auxiliary phase field. Such an approach is known to be thermodynamically consistent and recent works have illustrated its potential to be predictive for a wide range of fracture problems. Phase field for fracture models have succeeded in capturing complex crack patterns, including branching and coalescence in two and three dimensions [62][63][64][65][66].
In Figure 4, the solid Ω is a three-dimensional continuum where the main unknown variables are the displacement field, u, and the phase field, d. The contour conditions are the usual Dirichlet boundary ∂Ω u and Neumann boundary ∂Ω t conditions. function in terms of an auxiliary phase field. Such an approach is known to be thermodynamically consistent and recent works have illustrated its potential to be predictive for a wide range of fracture problems. Phase field for fracture models have succeeded in capturing complex crack patterns, including branching and coalescence in two and three dimensions [62][63][64][65][66].
In Figure 4, the solid Ω is a three-dimensional continuum where the main unknown variables are the displacement field, , and the phase field, . The contour conditions are the usual Dirichlet boundary Ω and Neumann boundary Ω conditions.  The formulation is presented for infinitesimal deformations in Equation (25): The total energy in the system is equal to the summation of the external and internal energy as shown in Equation (26): In a mechanical process, the internal energy is given by Equation (27): where the elastic energy density is given by Equation (28): where and are the Lame's elastic constants. Once a crack is produced in the crack set Γ, the total energy can be expressed as in Equation (29) |Γ is the energy associated with the fracture at Γ, which is in function of the fracture toughness as expressed in Equation (30): The method requires an explicit surface tracking algorithm and a method to handle the discontinuities. The surface integral is approximated by a volume integral in Equation (31): ̃| Ω is the approximation to the energy released during the cracking process and is the crack-density function. The formulation is presented for infinitesimal deformations in Equation (25): The total energy in the system is equal to the summation of the external and internal energy as shown in Equation (26): In a mechanical process, the internal energy is given by Equation (27): where the elastic energy density is given by Equation (28): where λ e and µ e are the Lame's elastic constants. Once a crack is produced in the crack set Γ, the total energy can be expressed as in Equation (29): ψ f racture|Γ is the energy associated with the fracture at Γ, which is in function of the fracture toughness ζ c as expressed in Equation (30): The method requires an explicit surface tracking algorithm and a method to handle the discontinuities. The surface integral is approximated by a volume integral in Equation (31): ψ internal|Ω is the approximation to the energy released during the cracking process and γ is the crack-density function.
For the definition of the crack-density function, γ, Allen-Cahn approximation to the crack surface is used, assuming a smooth transition from d = 1 along the crack set to d = 0 away from the crack set.
In Equation (32), w(d) ∈ C 2 ([0, 1]) is the local dissipation function and l is a regularization parameter carrying units of length.
The system can be resolved by minimizing the total energy, Equation (33):

Govern and Constitutive Equations
The continuum mechanics problem is expressed by the equation of equilibrium and the contour conditions, Equations (34)-(37): The constitutive laws are, Equations (38) and (39): whereσ is the degraded stress tensor, ξ is the thermodynamic conjugate to ∇d, κ = 2l 2 is the interfacial coefficient, and M = ζ c c 0 l is often referred to the mobility in keeping with general Allen-Cahn phase field models.

Discrete Element Method (DEM)
DEM has been used to simulate the behavior of granular materials and clayey soils. DEM considers individual particles and their interactions, enabling the simulation of crack formation during drying. DEM has proven effective in capturing the deformation and interaction between soil particles during drying, but its applicability to large-scale problems can be limited due to high computational costs [67][68][69][70].
While the DEM can accurately simulate the behavior of granular materials, it does not solve the macroscopic equations of balance and equilibrium. Instead, it focuses on capturing the microscale interactions between individual particles and their resulting collective behavior.
The limitations of the DEM method are that, firstly, DEM requires a substantial number of discrete particles to accurately represent the soil structure, making it computationally demanding for large-scale simulations.
Additionally, the accurate characterization of material properties, such as particleparticle interactions, contact forces, and soil-water interaction, can be challenging in clayey soils. The calibration of DEM parameters specific to clayey soils is often complex and time-consuming. Furthermore, DEM struggles to capture intricate crack patterns and accurately predict crack propagation due to its discrete nature. The method may also overlook important factors, such as the influence of soil matrix behavior and complex moisture redistribution phenomena, which are crucial in desiccation crack simulations.
Thus, while DEM offers valuable insights into the microscale behavior of soil particles, its limitations need careful consideration and validation when simulating desiccation cracks in clayey soils.

DEM Overview
In DEM, the soil is discretized by rigid particles that interact. For example, in [67], the contact behavior is modeled by a system consisting of springs, dashpots, and sliders ( Figure 5). The motion of each particle i is governed by the following force and moment balance in Equations (40) and (41): .. → . → ω i are the translational and rotational accelerations of a particle i, respectively; m i is the total mass of the particle i; → g is the acceleration of gravity; I i is the moment of inertia of the particle i; finally, the force ( → F i ) and moment ( → M i ) acting on particle I can be calculated as follows: Geotechnics 2023, 3, FOR PEER REVIEW 12 factors, such as the influence of soil matrix behavior and complex moisture redistribution phenomena, which are crucial in desiccation crack simulations. Thus, while DEM offers valuable insights into the microscale behavior of soil particles, its limitations need careful consideration and validation when simulating desiccation cracks in clayey soils.

DEM Overview
In DEM, the soil is discretized by rigid particles that interact. For example, in [67], the contact behavior is modeled by a system consisting of springs, dashpots, and sliders ( Figure 5). The motion of each particle i is governed by the following force and moment balance in Equations (40) and (41): where ⃗⃗⃗⃗̈ and ⃗⃗⃗⃗⃗ are the translational and rotational accelerations of a particle i, respectively; is the total mass of the particle i; ⃗ is the acceleration of gravity; is the moment of inertia of the particle i; finally, the force ( ⃗⃗⃗ ) and moment ( ⃗⃗⃗⃗⃗ ) acting on particle I can be calculated as follows: In Equations (42) and (43), ⃗⃗⃗⃗ and ⃗⃗⃗⃗ are the normal and shear contact forces at contact c; ⃗⃗⃗⃗⃗ and ⃗⃗⃗⃗⃗ are the normal and shear dashpot forces or damping forces; N is the total number of particle i contacts; ⃗⃗⃗ is the vector from the center of particle i to contact c; and is the mass corresponding with the particle i. By having particles that represent the soil and particles that represent water, this model [67] can reproduce the desiccation process of a thin layer of clay capturing the characteristic curling behavior during drying. are the normal and shear dashpot forces or damping forces; N is the total number of particle i contacts; → l c is the vector from the center of particle i to contact c; and m i is the mass corresponding with the particle i.
By having particles that represent the soil and particles that represent water, this model [67] can reproduce the desiccation process of a thin layer of clay capturing the characteristic curling behavior during drying.

Cellular Automaton Method (CAM)
CAM (singular: Cellular Automaton or plural: Cellular Automata) has been employed to simulate the growth and pattern formation of cracks in rocks/clays and is able to work from micro-to macroscale levels [71][72][73][74]. This method uses a grid of cells with different states to represent the crack initiation, propagation, and interaction. CAM offers a simplified representation of crack evolution and is computationally efficient, allowing for the simulation of large-scale crack patterns. The method can simulate heterogeneity in the soil; however, it may lack accuracy in capturing the mechanical behavior of the soil. CAM represents the soil as a 2D or 3D grid of cells and uses rules to simulate the drying process and the resulting stress distribution. The model calculates the evolution of the system over time based on these rules and the initial conditions specified.
This method is what is called a mechanistic method that does not necessarily use physical-based equations to establish its rules. Thus, one limitation of CAM is the challenge of accurately representing the complex behavior of clayey soils within the simplified cellular automata framework. CAM relies on predefined rules and assumptions, which may not fully capture the intricacies of crack formation and propagation in clayey soils. Additionally, the model's grid-based nature may result in limited spatial resolution, potentially overlooking fine-scale details of crack patterns.
These limitations underscore the need for cautious interpretation and validation of results when applying CAM to simulate desiccation cracks in clayey soils, as well as the potential for combining CAM with other methods to address these shortcomings.

Formulation of CAM
As an example of formulation, the rock discontinuous automaton (RDCA) from [71] is resumed in this section. RDCA includes: (1) Strong discontinuity; (2) Tracking a moving crack; (3) Cellular automaton updating rule; (4) Yield criteria and crack growth laws.
In Figure 6, this method generates a grid of cells to discretize the continuum. In this example, the grids are squares but the grids can have other shapes (triangles, hexagons, etc.). The elements of the grid are of three types: normal elements (white), elements intersected by a crack (yellow), and elements containing a crack tip (light blue).

Cellular Automaton Method (CAM)
CAM (singular: Cellular Automaton or plural: Cellular Automata) has been employed to simulate the growth and pattern formation of cracks in rocks/clays and is able to work from micro-to macroscale levels [71][72][73][74]. This method uses a grid of cells with different states to represent the crack initiation, propagation, and interaction. CAM offers a simplified representation of crack evolution and is computationally efficient, allowing for the simulation of large-scale crack patterns. The method can simulate heterogeneity in the soil; however, it may lack accuracy in capturing the mechanical behavior of the soil. CAM represents the soil as a 2D or 3D grid of cells and uses rules to simulate the drying process and the resulting stress distribution. The model calculates the evolution of the system over time based on these rules and the initial conditions specified.
This method is what is called a mechanistic method that does not necessarily use physical-based equations to establish its rules. Thus, one limitation of CAM is the challenge of accurately representing the complex behavior of clayey soils within the simplified cellular automata framework. CAM relies on predefined rules and assumptions, which may not fully capture the intricacies of crack formation and propagation in clayey soils. Additionally, the model's grid-based nature may result in limited spatial resolution, potentially overlooking fine-scale details of crack patterns.
These limitations underscore the need for cautious interpretation and validation of results when applying CAM to simulate desiccation cracks in clayey soils, as well as the potential for combining CAM with other methods to address these shortcomings.

Formulation of CAM
As an example of formulation, the rock discontinuous automaton (RDCA) from [71] is resumed in this section. RDCA includes: (1) Strong discontinuity; (2) Tracking a moving crack; (3) Cellular automaton updating rule; (4) Yield criteria and crack growth laws.
In Figure 6, this method generates a grid of cells to discretize the continuum. In this example, the grids are squares but the grids can have other shapes (triangles, hexagons, etc.). The elements of the grid are of three types: normal elements (white), elements intersected by a crack (yellow), and elements containing a crack tip (light blue).

Strong discontinuity
To capture the discontinuous displacement field produced by cracks, Equation (44) was selected [75]:

Strong discontinuity
To capture the discontinuous displacement field produced by cracks, Equation (44) was selected [75]: where N represents the total nodal number of the element; N dis is the total nodal number on an element intersected completely by a crack; N asy1 and N asy2 are the sets of nodes associated with crack tips 1 and 2 in their influence domain, respectively; N i , N j and N k are the shape functions of the associated node; u i is the nodal displacements (standard degrees of freedom); and a j , b α1 k and b α2 k are the vectors of the additional degrees of nodal freedom for modeling crack faces and the two crack tips, respectively.
H(ξ(x)) is the Heaviside enrichment function to model the strong discontinuity. It is defined as in Equation (45): where ξ is the value of the level set function.

Tracking a moving crack
To track the moving crack, the level set method is used considering the domain Ω is divided into two non-overlapping subdomains, Ω 1 and Ω 2 , that share an interphase Γ as seen in Figure 7a. The level set function φ(x) is defined [76] as in Equation (46): Geotechnics 2023, 3, FOR PEER REVIEW 14 where represents the total nodal number of the element; is the total nodal number on an element intersected completely by a crack; 1 and 2 are the sets of nodes associated with crack tips 1 and 2 in their influence domain, respectively; , and are the shape functions of the associated node; is the nodal displacements (standard degrees of freedom); and , 1 and 2 are the vectors of the additional degrees of nodal freedom for modeling crack faces and the two crack tips, respectively.
( ( )) is the Heaviside enrichment function to model the strong discontinuity. It is defined as in Equation (45): where is the value of the level set function.

Tracking a moving crack
To track the moving crack, the level set method is used considering the domain Ω is divided into two non-overlapping subdomains, Ω 1 and Ω 2 , that share an interphase Γ as seen in Figure 7a. The level set function ( ) is defined [76] as in Equation (46): An option for the level set function can be simply defined in terms of the signed distance function: where d is the normal distance from a point x to the interphase Γ.

Cellular automaton updating rule
In terms of the Cellular Automaton updating rule for continuous problems, they can be elastic, elasto-plastic, elasto-visco plastic, etc. [77]. Discontinuities introduce more complexity. In this method, a cell has nodes , related cell elements and neighbor cell nodes. For every cell, the equilibrium equation is An option for the level set function can be simply defined in terms of the signed distance function: where d is the normal distance from a point x to the interphase Γ.

Cellular automaton updating rule
In terms of the Cellular Automaton updating rule for continuous problems, they can be elastic, elasto-plastic, elasto-visco plastic, etc. [77]. Discontinuities introduce more complexity. In this method, a cell has nodes N i , related cell elements and neighbor cell nodes. For every cell, the equilibrium equation is In Equation (48), K ij is the local stiffness on the cell node, ∆u j is the incremental degree of freedom, and ∆F i is the nodal force. The dimension of K ij , ∆u j , and ∆F i are dependent on the cell type, i.e., a normal cell or a discontinuous cell. When a cell is intersected by a crack, the cell should be enriched by an enrichment function. The increment nodal force can be obtained from the equation In Equation (49), B j im is the cell element stiffness. For different element types, such as normal elements, elements with intersecting cracks or elements with crack tips, the cell element stiffness differs.

Yield criteria and crack growth laws
Different yield criteria can be implemented in this method, Mohr-Coulomb, Drucker-Prager, etc., so, this is suitable for rocks and soils.
For the crack growth, several laws can be used, including maximum circumferential stress, maximum tangential stress, minimum strain energy criteria, etc.
For example, the maximum circumferential stress criterion will be In Equation (50), θ c is measured with respect to a local polar coordinate system with its origin at the crack tip and is aligned with the direction of the existing crack.

Governing Equations and Local Stiffness Matrix
The problem is solved in the domain Ω and the boudary Γ = Γ t ∪ Γ u ∪ Γ c as can be seen in Figure 7b. ∇·σ In Equations (51)- (56), σ is the Cauchy stress tensor; u is the displacement; b is the body force per unit volume; n is the unit outward normal; t and u are the prescribed traction and displacement, respectively; ε is the strain tensor; and C is the Hooke tensor.

Integration of Methods to Improve Simulations and Analysis
All the methods in the previous section share limitations that encompass difficulties in accurately characterizing material properties and behavior, representing complex interactions between soil particles, cracks, and fluid flow, and addressing computational demands, particularly for large-scale simulations. Challenges also arise in capturing intricate crack patterns, accurately predicting crack propagation, and incorporating the mechanical behavior of the soil matrix. Furthermore, simulating soil-structure and soil-atmosphere interactions, moisture redistribution, and the microstructure of clayey soils can be challenging. The limitations underscore the need for careful consideration of method selection, calibration of constitutive models, and the exploration of techniques to overcome these challenges and enhance the accuracy and reliability of desiccation crack simulations in clayey soils.
Researchers have explored the combination of multiple numerical methods to simulate the process of desiccation cracks in clayey soils. By combining different methods, they aim to leverage the strengths of each approach to overcome individual limitations and improve the overall accuracy and reliability of the simulations.
For example, the coupling of Scaled Boundary Finite Element Method (SBFEM) with adaptive Phase Field Modeling (PFM) was published in 2019 [78]. This method can model materials that show quasi-brittle fracture properties and reduces the computational effort required by the phase field model. The FEM provides an accurate representation of soil deformation and handles the crack zone by densifying the mesh in the area, while the PFM handles the evolution and propagation of cracks. For hydraulic fracture propagation in rocks, a hydromechanical coupled LBM-DEM model was published in 2018 [79]. This is an interesting combination of physical-based and non-physical based models.
The FEM was coupled with CAM to simulate residual stress of dual-phase steel subjected to laser treatments in a work in 2021 [80]. CAM simulated the microstructure orientation that was after the input for a thermos-mechanical finite element model. This model was able to capture the anisotropy of the material, the dual-phase ferrite-martensite, the temperature sensitivity, and the strain rate influence in the thermos-mechanical behavior. In 2007, three techniques were used to study wave propagation problems. LBM, FEM, and CAM [81] were used in different parts of the geometry of the problem to take advantage of the power of every method.
Hybrid approaches allow for the simultaneous modeling of soil deformation and crack propagation, considering the discrete behavior of soil particles or the fluid flow within the soil matrix. This combination enables capturing both the macroscale behavior of the soil structure and the microscale interactions between particles or fluid.
Additionally, researchers have explored the integration of different methods using a multi-scale double threshold segmentation algorithm [82], methods like free-mesh smoothed particle hydrodynamics method (SPH) [83] and particle discretization scheme finite element method (PDS-FEM) [52]. This involves coupling methods such as FEM, DEM, or LBM at different length scales to capture the behavior of the soil from the microscale to the macroscale. This allows for a more comprehensive understanding of desiccation crack formation and evolution.
While the combination of multiple methods shows promise, it is still an active area of research, and the specific combinations and approaches vary depending on the research objectives and available computational resources.
FEM and CAM are the only methods able to work from micro-to macroscale levels being computationally efficient for large-scale simulation. In the opinion of the author of this paper, FEM is the best method to simulate the desiccation process taking into consideration the complexities of the soil behavior and is a physical-based method.

Finite Element and Cellular Automaton Method (FEM-CAM)
After the study of the state of the art, the author arrives to the conclusion that since the problem of desiccation cracks in clayey soils is a multiphysics and multiphase problem (soil matrix + water and air in the pores), it can be resolved using FEM for the THM process and CAM for the cracking problem by developing a FEM-CAM method. The method proposed in this section (for future research) is based on the one proposed by the authors of [80] to simulate residual stress by combining FEM and CAM.

Integration of FEM with CAM to Simulate Desiccation Cracks in Clayey Soils
CAM models the soil as a grid of cells, and each cell can be in different states representing soil moisture content, stress, or cracking. The mathematical formulation of the desiccation crack problem using CAM involves expressing the evolution of moisture content, suction, temperature, and stress, in the soil domain that can be calculated by a THM or HM model resolved by FEM.
In Figure 8, a simplified flow chart shows the coupling between FEM and CAM to solve the desiccation cracks problem. To start the resolution of the problem, the soil domain must be defined by including the geometry, contour conditions, loads, and material properties. FEM is used to discretize the soil domain and then this mesh is used to define the CAM grid and then define the initial conditions. As it was shown in the previous examples of the review, CAM first produces a grid representation of the soil and can establish heterogeneity.
Geotechnics 2023, 3, FOR PEER REVIEW 17 content, suction, temperature, and stress, in the soil domain that can be calculated by a THM or HM model resolved by FEM.
In Figure 8, a simplified flow chart shows the coupling between FEM and CAM to solve the desiccation cracks problem. To start the resolution of the problem, the soil domain must be defined by including the geometry, contour conditions, loads, and material properties. FEM is used to discretize the soil domain and then this mesh is used to define the CAM grid and then define the initial conditions. As it was shown in the previous examples of the review, CAM first produces a grid representation of the soil and can establish heterogeneity. A cracking criterion is defined based on stress thresholds or stress gradients. The criterion determines when a cell transitions from an intact state to a cracked state. For example, a simple criterion could be the crack initiate when the tensile strength is reached in any cell. Once a cell transitions to a cracked state, crack propagation rules determine how cracks propagate to neighboring cells. This can be based on stress redistribution or local crack propagation rules. The direction and extent of crack propagation can be influenced by factors such as stress concentration, crack coalescence, and crack branching.
One the CAM model is updated with the new cracks, the moisture content, suction, stress, temperature, and cracking states are updated at each time step using FEM. The A cracking criterion is defined based on stress thresholds or stress gradients. The criterion determines when a cell transitions from an intact state to a cracked state. For example, a simple criterion could be the crack initiate when the tensile strength is reached in any cell. Once a cell transitions to a cracked state, crack propagation rules determine how cracks propagate to neighboring cells. This can be based on stress redistribution or local crack propagation rules. The direction and extent of crack propagation can be influenced by factors such as stress concentration, crack coalescence, and crack branching. One the CAM model is updated with the new cracks, the moisture content, suction, stress, temperature, and cracking states are updated at each time step using FEM. The specific equations and constitutive relationships used will depend on the chosen modeling approach and the characteristics of the soil being studied.
Finally, a direct, iterative, or embedded strategy is necessary to couple FEM with CAM. The degree of convergence of every time step will depend on the accuracy needed to resolve specific problems. The final step in the method will be the postprocessing of the results in terms of stresses, strains, humidity, temperature, etc., and should show the final configuration of cracks in the soil.

Conclusions
Modeling and simulating desiccation cracks in clayey include cycles of drying, wetting, flooding, and re-drying ( Figure 1) that will complicate the formulation and implementation of numerical models in the future even more.
The five methods reviewed in this paper include the Finite Element (FEM), Lattice Boltzmann (LBM), Phase Field (PFM), Discrete Element (DEM), and Cellular Automaton (CAM) Methods, whereby each one has their strengths and limitations when applied to this problem.
FEM captures very well and with consistency the THM desiccation process at a macroscale level and can be applied to micro-and mesoscale levels also. Being based on the continuum mechanics equations makes the method reliable but limited when dealing with complex crack patterns due to the need for complex remeshing that are computationally highly demanding. For example, Ref. [32] presents the simulation of the initiation and propagation of only one crack in 2D, showing a good agreement with laboratory results in terms of water loss and shrinkage. LBM is a good method to simulate the flow of the initial slurry of clayey soil at a mesoscale level but not so well for deformation in connection with the complex THM behavior of clayey soils when it becomes a compacted unsaturated medium. However, Ref. [50] demonstrated the ability of simulating multiple cracks in 2D. PFM is good to simulate micro-and mesoscale continuous cracks, but not so good to capture the THM nature of the process in the soil mass. Again, Ref. [55] has shown the ability to simulate multiple cracks in 2D. DEM is good at simulating the behavior of particles interacting at a microscale level but not so good at capturing the THM nature of the soil behavior and is not a physical-based model. However, Ref. [67] demonstrated that this method can simulate the curling process when clays dry. CAM is good for simulating complex crack patterns at micro-, meso-and macroscale levels, but limited in dealing with the THM desiccation process, which is what makes it ideal to work in combination with FEM. In [71], the ability of this method to simulate cracks in concrete, which can be extrapolated to clays, is demonstrated.
The limitations of these methods working separately, such as introducing heterogeneity, difficulties in accurately characterizing material properties, capturing intricate crack patterns, considering complex soil-fluid, soil-structure, soil-atmosphere interactions, and being computationally highly demanding, highlight the need for further research.
Since the desiccation cracks in clayey soils is a problem well divided into two coupled processes, desiccation and cracking, and considering that FEM simulates the desiccation well, and other methods simulate the cracking process well, their combination seems to be the best option. Then, to address these limitations and improve the simulations, researchers started to explore the combination of different methods. By integrating the strengths of multiple approaches, such as coupling FEM with DEM or LBM, or combining FEM with PFM, it becomes possible to overcome individual limitations and obtain more accurate and reliable simulations. Furthermore, a multi-scale approach, integrating methods at different length scales, allows for a comprehensive understanding of desiccation crack formation and evolution.
The combination of FEM and CAM, presented in Section 4, is a promising alternative since FEM and CAM can capture the process from micro-to macroscale levels and CAM is particularly efficient when computationally intense calculations are needed for largescale crack pattern simulations. In addition, since they both have the same governing equations and can share the mesh, it seems that they can efficiently work together to resolve desiccation cracks in clayey soils. It is of course difficult to know which method or combination of methods will be the best option to tackle the desiccation crack problem in clayey soils in the coming years, but the aim of this paper is to identify the most promising alternatives based on the research performed up until today.