A State of the Art: Numerical Methods to Simulate Desiccation Cracks in Clayey Soils

1. Introduction

After more than a century of research from experimental [1,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55], theoretical [2,49,62,63,64], and numerical points of view [3,56,57,58,59,60,61], desiccation cracks in clayey soils are still an open research field due to their complexity. It has two very different components, desiccation, which is the loss of water due to water evaporation, and cracking, a failure produced when reaching the strength of the soil. The first is a thermo-hydromechanical (THM) problem [4] and the second is a fracture mechanics (FM) problem [5]. This topic is important because when clayey soil desiccates and cracks, its properties change becoming more permeable and less strong against loads.

Simulating desiccation cracks in clayey soils is a complex task due to several reasons.

Firstly, clayey soils are multiphase mediums 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 [4].

Secondly, clayey soils exhibit coupled nonlinear THM behavior when drying then is a Multiphysics problem. 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 [6]. 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 [4].

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 [7,15,20].

Fourthly, simulating desiccation cracks often requires considering different scenarios, from laboratory specimens to field-scale applications. Accounting for scale effects and capturing the heterogeneity of soil properties is crucial for realistic simulations [8].

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.

In section 4.1 the equations for the THM problem are presented. This problem involves several interconnected physical processes in a portion of a multiphase soil system that is in mechanical equilibrium, equation (1). The shrinkage that takes place as moisture is lost is governed by the principle of balance, known as Richards' equation (2), and the generalized Darcy’s Law (10). At the same time, there will be a distribution of temperature governed by the heat balance equation (3) and Fourier’s law (12). Mechanical deformation arises from shrinkage, exhibiting elastic and plastic behavior described by the THM constitutive model, equation (7).

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 [4], and the diffusion of air in water is governed by Fick’s law [4], 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 [9,10,11], vapor effusion, and Stefan’s flow [12]. Fortunately, in many cases, they can be neglected and continue capturing the main mechanisms that govern the desiccation and cracking processes. So, they are not necessarily needed in a formulation and a numerical model [15].

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 [13] 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 [5]. LEFM principles and its rules determine the critical crack length and assess crack propagation [6]. 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 FEM [6]. For this reason, several approaches to simulate the cracking process have been proposed apart from the FEM. LBM, DEM, CAM, and PFM are all alternatives to the FEM that have been used to effectively simulate the desiccation cracks in clayey soils and cracks in other problems [16,17,18,19,22,23,24,25,26,29].

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 [15]. 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.

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 reviewed in detail as a promising alternative.

2. Methods to simulate desiccation cracks in clayey soils

In the attempt of simulating desiccation cracks in soils, researchers have used physical-based 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 main characteristics, strengths, and limitations. The strengths of these methods are presented in Table 1.

Physical-based Models

2.1. The Finite Element Method (FEM)

FEM and its variant XFEM [15,21] have been extensively used for simulating desiccation cracks in clayey soils at a macroscale level and can be applied too 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 (Section 4).

FEM has been successful in capturing the complex behavior of desiccation cracks since it is able to 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 [15].

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.

2.2. Lattice Boltzmann Method (LBM)

LBM was used to study clayey soils undergoing desiccation at a mesoscale level [24] and it was introduced to simulate fluid flow for the first time in 1986 [39].

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 [16,22,23]. LBM can solve physical equations of balance and equilibrium. LBM is based on the Boltzmann equation that is derived from the principles of conservation of mass, momentum, and energy, hence, a physical-based model.

The limitations of the LBM are that 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.

2.3. 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 [17,24,25,29,34].

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.

PFM is commonly employed to simulate phase transitions and evolving interfaces in various physical systems, such as solidification, solid-state transformations, and fluid dynamics, in combination with LBM [35].

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.

Non-physical-based Models

2.4. 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 [18,26]. 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.

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 particle-particle 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.

2.5. Cellular Automaton Method (CAM)

CAM (Cellular Automaton or Cellular Automata) has been employed to simulate the growth and pattern formation of desiccation cracks in clayey soils and is able to work from micro to macro scale levels [19,32,33]. 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. So, one limitation of CAM is the challenge of accurately representing the complex behavior of clayey soils within the simplified cellular automaton 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.

Nothing stops scientists from improving CAM to use as rules the classic physical-based equations, but this is something that is yet an open research area.

3. Integration of methods to improve simulations and analysis

All the methods in the previous section share limitations 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 [34]. 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 FEM with DEM [27] or LBM [28] has been investigated. This hybrid approach allows 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.

Another approach involves combining FEM with PFM [30] to simulate crack propagation. The FEM provides an accurate representation of soil deformation, while the PFM handles the evolution and propagation of cracks. This combination enables the simulation of complex crack patterns and the prediction of crack paths without explicitly tracking them.

Additionally, researchers have explored the integration of different methods using a multi-scale approach [31]. 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 macro-scale. 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 macro scale levels being computationally efficient for large-scale simulation. 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. For this reason, in the next section, the fundamental equations that combine FEM and CAM are presented as a promising approach to tackling the desiccation cracks in clayey soils.

4. Finite Element and Cellular Automaton Method (FEM-CAM)

The problem of desiccation cracks in clayey soils is a Multiphysics and multiphase problem (soil matrix + water and air in the pores) that can be resolved using FEM for the THM process and CAM for the cracking problem by developing a FEM-CAM method.

4.1. Govern and constitutive equations for desiccation in clayey soil problems.

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 interaction 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.

Governing equations

4.1.1. Equilibrium equation (Cauchy Equation of Motion)

If no dynamic effects are considered, the equilibrium equation of the soil matrix is as follows.

$$\nabla ·\mathsf{\sigma }+\rho \mathbf{g}=\mathbf{0}$$

(1)

Equation (1) is an elliptic partial differential equation where, $\mathsf{\sigma }$, is the total stress tensor, and $\rho$ is the average density of the multiphase medium (soil, water, and air). The vector $\mathbf{g}$ is the gravity vector.

4.1.2. 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:

$$\nabla ·\left({\rho }^w\mathbf{q}\right)+\frac{\partial }{\partial t}\left({\rho }^{w}n{S}_r\right)=0$$

(2)

In equation (2), ${\rho }^w$ is the water density, $\mathbf{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.

4.1.3. 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.

$$\nabla ·\left({\mathrm{K}}^{\theta }\nabla \theta \right)+{\mathbf{q}}^{\theta }-{\rho }^s\mathrm{c}\frac{d\theta }{dt}=0$$

(3)

In equation (3), ${\rho }^s$ is the density of the soil, $c$ is the specific heat capacity, $\theta$ is the temperature, ${\mathbf{q}}^{\theta }$ is the heat transfer rate, and $K^{\theta }$ is the thermal conductivity (could be scalar for isotropic permeability or tensorial for anisotropic permeability.

Constitutive equations

4.1.4. Stress-strain thermos-mechanical constitutive law

For the most general case of unsaturated soils, the effective stress tensor (${\mathsf{\sigma }}^{\mathrm{\text{'}}}$) is:

$${\mathsf{\sigma }}^{\mathrm{\text{'}}}=\mathsf{\sigma }-u_a\mathbf{1}+\chi \left(u_a-u_w\right)\mathbf{1}$$

(4)

In equation (4) $\mathit{\sigma }$ is the total stress tensor. The air and water pressure are respectively $u_a$ and $u_w$, $\chi$ is a parameter that depends on the degree of saturation, the stress history and the soil’s fabric and $\mathbf{1}\equiv {\delta }_{ij}$, is the identity tensor.

In this formulation, the matrix suction and the net mean stress define the effective stress tensor ${\mathsf{\sigma }}^{\mathrm{\text{'}}}$ through equation (4).

The net stress ${\mathsf{\sigma }}^{net}$ and the suction $s$ are:

$$1.{\mathsf{\sigma }}^{net}=\mathsf{\sigma }-u_a\mathbf{1}$$

(5)

$$2.s=u_a-u_w$$

(6)

The general strain-stress relation must be written in differential form, because of the nonlinearity of the material behavior.

$$d\mathsf{\sigma }=\mathbf{D}(\mathit{\epsilon },\theta ,s)d\mathsf{\epsilon }$$

(7)

For the most general THM case, $\mathbf{D}$ is a tangent matrix in function of the strain, $\mathit{\epsilon }$, temperature, $\theta$, and suction, $s$. Equation (7) establishes the coupling between temperature, suction, and mechanical effects.

The deformations are calculated by 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:

$$d\mathsf{\epsilon }=d{\mathsf{\epsilon }}^{net}+d{\mathsf{\epsilon }}^s+d{\mathsf{\epsilon }}^{\theta }=\mathbf{C}\left(\mathbb{K},\mathbb{G}\right)d{\mathsf{\sigma }}^{net}+\mathbf{h}\left({\mathbb{K}}^s\right)ds+\mathbf{t}\left({\mathbb{K}}^{\theta }\right)d\theta$$

(8)

In equation (8), the parameter $\mathbb{K}$ is the volumetric modulus and $\mathbb{G}$ is the shear modulus of the soil matrix; the parameter, ${\mathbb{K}}^s$ is the volumetric modulus due to changes in suction; the parameter, ${\mathbb{K}}^{\theta }$ is the volumetric modulus due to changes in temperature. These parameters must be established depending on the constitutive model chosen (linear elastic, non-linear elastic, viscoelastic, plastic, etc.) the factor $\mathbf{C}$ is a 4th order compliance tensor and $\mathbf{h},\mathbf{t}$ are 2nd order tensors.

The net stress increments can be obtained from (8):

$$d{\mathsf{\sigma }}^{net}={\mathbf{C}}^{-1}\left(d\mathsf{\epsilon }-\mathbf{h}\left({\mathbb{K}}^s\right)ds-\mathbf{t}\left({\mathbb{K}}^{\theta }\right)d\theta \right)=\mathbf{D}\left(d\mathsf{\epsilon }-\mathbf{h}\left({\mathbb{K}}^s\right)ds-\mathbf{t}\left({\mathbb{K}}^{\theta }\right)d\theta \right)$$

(9)

In equation (9), $\mathbf{D}={\mathbf{C}}^{-1}$, is the tangent stiffness tensor.

4.1.5. Generalized Darcy’s law for unsaturated soils and permeability tensor.

The generalized Darcy’s law for unsaturated soils is:

$$\mathbf{q}=-\mathbf{K}\left(S_r\right)·\left(\nabla s-\mathbf{g}{\rho }^w\right)$$

(10)

In equation (10) $\mathbf{q}$ is the velocity of Darcy vector; $\nabla s$ is the gradient of the suction; $\mathbf{K}\left(S_r,n,\theta \right)$ is a permeability tensor that changes with water saturation degree, $S_r$, porosity, $n$, and temperature, $\theta$; $\mathbf{g}$ is the gravity vector and ${\rho }^w$ is the water density. The permeability tensor, $\mathbf{K}$, can be isotropic or anisotropic.

4.1.5.1. 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$.

$$S_r={\left[1+{\left(\frac{s}{P_0·f_n}\right)}^{\frac{1}{1-\lambda }}\right]}^{-\lambda }{f}_n=\mathrm{e}\mathrm{x}\mathrm{p}\left[-\eta \left(n-n_0\right)\right]$$

(11)

Where, $\lambda$, 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 by means of a parameter, $\eta$. For non-deformable soils, $f_n=1$, because porosity is constant.

4.1.6. Fourier’s law

This is the constitutive law for the thermal problem.

$${\mathbf{q}}^{\theta }=-{\mathrm{K}}^{\theta }\nabla \mathsf{\theta }$$

(12)

In equation (12), ${\mathbf{q}}^{\theta }$ is the heat transfer rate and ${\mathrm{K}}^{\theta }$ is the thermal conductivity.

4.2. 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.

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.

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 [36,37,38].

Finally, a direct, iterative, or embedded strategy is necessary to couple FEM with CAM.

4.3. 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. 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 equations (13) in what is known as a u-p formulation [6]. This assumption simplifies the problem considerably since equations (3) and (12) are not needed and equations (7), (8) are (9) are less complex.

$$\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{H}\end{array}\right]\left[\begin{array}{c}\overline{\mathbf{u}}\\ \overline{\mathbf{p}}\end{array}\right]+\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{T}}& {\mathbf{Q}}_{\mathrm{T}}\\ \mathbf{P}& \mathbf{S}\end{array}\right]\left[\begin{array}{c}\frac{d\overline{\mathbf{u}}}{dt}\\ \frac{d\overline{\mathbf{p}}}{dt}\end{array}\right]=\left[\begin{array}{c}\frac{d{\mathbf{f}}^u}{dt}\\ {\mathbf{f}}^p\end{array}\right]$$

(13)

In equation (13), $\mathit{u}$, are the displacements and $\mathit{p}$, is the suction (or negative pore water pressure is the air pressure is considered zero).

The factors in equation (13) are:

Permeability Matrix:

$$\mathbf{H}=\underset{\mathsf{\Omega }}{\overset{}{\int }}{\left(\nabla {\mathbf{N}}_{\mathrm{p}}\right)}^{\mathrm{T}}\mathbf{K}\left(S_r\right)\nabla {\mathbf{N}}_{\mathrm{p}}d\mathsf{\Omega }$$

(14)

Stiffness Matrix:

$${\mathbf{K}}_{\mathrm{T}}=\underset{\mathsf{\Omega }}{\overset{}{\int }}{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\mathbf{B}d\mathsf{\Omega }$$

(15)

Coupling Matrix:

$${\mathbf{Q}}_{\mathrm{T}}=\underset{\mathsf{\Omega }}{\overset{}{\int }}\frac{1}{{3\mathbb{K}}_t^s}{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\mathbf{m}{\mathbf{N}}_{\mathrm{p}}d\mathsf{\Omega }$$

(16)

Coupling Matrix:

$$\mathbf{P}=\underset{\mathsf{\Omega }}{\overset{}{\int }}{\left({\mathbf{N}}_{\mathrm{p}}\right)}^{\mathrm{T}}{S}_r{\mathbf{m}}^{\mathrm{T}}\mathbf{B}d\mathsf{\Omega }$$

(17)

Compressibility Matrix:

$$\mathbf{S}=\underset{\mathsf{\Omega }}{\overset{}{\int }}{\left({\mathbf{N}}_{\mathrm{p}}\right)}^{\mathrm{T}}n\frac{\partial S_r}{\partial u_w}{\mathbf{N}}_{\mathrm{p}}d\mathsf{\Omega }+\underset{\mathsf{\Omega }}{\overset{}{\int }}{\left({\mathbf{N}}_{\mathrm{p}}\right)}^{\mathrm{T}}\frac{n{S}_r}{K^w}{\mathbf{N}}_{\mathrm{p}}d\mathsf{\Omega }$$

(18)

Vector of Nodal Forces:

$$\frac{\partial {\mathbf{f}}^u}{\partial t}=\underset{\mathsf{\Omega }}{\overset{}{\int }}{\mathbf{N}}_u\rho \frac{\partial \mathbf{g}}{\partial t}d\mathsf{\Omega }+\underset{\mathsf{\Omega }}{\overset{}{\int }}{\mathbf{N}}_u\frac{\partial \overline{\mathbf{t}}}{\partial t}d\mathsf{\Omega }$$

(19)

Vector of Nodal Flow:

$${\mathbf{f}}^p=\underset{\mathsf{\Omega }}{\overset{}{\int }}{\rho }^w{\left(\nabla {\mathbf{N}}_{\mathrm{p}}\right)}^{\mathrm{T}}\mathbf{K}\left(S_r\right)\mathbf{g}d\mathsf{\Omega }-\underset{\mathsf{\Gamma }}{\overset{}{\int }}{\left({\mathbf{N}}_{\mathrm{p}}\right)}^{\mathrm{T}}{q}^{w}d\mathsf{\Gamma }$$

(20)

5. 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 on this paper include the Finite Element (FEM), Lattice Boltzmann (LBM), Phase Field (PFM), Discrete Element (DEM), and Cellular Automaton (CAM) Methods have each their strengths and limitations when applied to this problem.

FEM captures very well and with consistency the THM desiccation process at macro scale level and can be applied to micro and mesoscale levels also. Being based on the continuum mechanics equations make the method reliable but limited when dealing with complex crack patterns due to the need of complex remeshing that are computational highly demanding. LBM is a good method to simulate the flow of the initial slurry of clayey soil at mesoscale level but not so good for deformation in connection with the complex THM behaviour of clayey soils when becomes a compacted unsaturated medium. 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. DEM is good at simulating the behavior of particles interacting at microscale level but not so good capturing the THM nature of the soil behavior and is not a physical-based model. CAM is good for simulating complex crack patterns at micro, meso and macroscale levels but limited to deal with the THM desiccation process.

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 computational highly demanding, highlight the need for further research.

Since the desiccation cracks in clayey soils is a problem well divided in two coupled processes, desiccation, and cracking, and considering that FEM simulate well the desiccation and other methods simulate well the cracking process, the combination of them 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 large-scale crack pattern simulations.

The future holds promising opportunities for the development and refinement of these combined approaches. Advancements in combining different methods and using complementary techniques will likely lead to more sophisticated simulations and a deeper understanding of desiccation cracks in clayey soils.

Each of these methods will continue evolving and improving. All of them have the potential to tackle the desiccation cracks in clayey soils and other challenging problems in isolation. Since the technology on computers continues improving, some of the limitations of these methods, like the computational requirements, will be reduced.

All these efforts will ultimately contribute to improved engineering practices, risk assessments, and mitigation strategies in various fields such as geotechnical engineering, soil mechanics, and environmental sciences since desiccation cracks in clayey soil is a topic that define clayey soils failure.