Stability of Expansive Soil Slopes under Wetting–Drying Cycles Based on the Discrete Element Method

: The swelling-shrinkage behavior of expansive soils under climate changes will cause the crack development, which can be destructive of expansive soil slopes. This study investigated the effect of drying/wetting cycles on the stability of an expansive soil slope using the discrete element method (DEM), in consideration of the crack development induced by climate changes. The strength reduction method was adopted in the DEM calculations, which was coupled with the unsaturated seepage analysis given by the finite element method. The slope stability and the failure model of the slope after different times of wetting–drying cycles were analyzed, and the results were compared with those calculated by the limit equilibrium method and the finite element method. The results indicated that the failure pattern of the expansive soil slope was strongly influenced by the wetting–drying cycles. A shallow sliding surface of the expansive soil slope occurred after several wetting–drying cycles. Similarly, the safety factor of the expansive soil slope decreased gradually with the wetting–drying cycles. Considering the cracks’ evolution inside the expansive soil slope from the drying/wetting cycles, a shallower sliding surface with a smaller safety factor was obtained from the strength reduction method of the DEM, in comparison with the two conventional methods of the Limit equilibrium method and finite element method. Therefore, cracks play an essential role in the expansive soil slope stability. The strength reduction method of the DEM, which considers the cracks’ evolution during drying/wetting cycles, is more reliable.


Introduction
Expansive soil is a problematic soil that swells when it absorbs water and shrinks when it loses water, leading to cracks.This behavior often triggers geological disasters and engineering problems, causing significant economic losses and even threatening people's lives [1,2].The engineering properties of expansive soil are very sensitive to external climate and environmental changes.Expansive soil slopes are exposed to climate changes in the long term.The wetting-drying cycles induced by climate changes not only deteriorate the engineering properties of expansive soil but also exacerbate crack formation due to expansion and shrinkage of the soil [3,4].Cracks break the integrity of slope soil, significantly reduce its shear strength, and provide pathways for rainwater infiltration [5,6].Therefore, it is essential to evaluate the long-term stability of expansive soil slopes under climate changes [7][8][9].
In recent decades, scholars worldwide have conducted comprehensive studies on the development of cracks in expansive soil slopes under the influence of wetting-drying cycles.Various studies have explored the failure mechanisms of these slopes using methods such as centrifugal or shrinkage model tests [10][11][12][13].For instance, Cai et al. [14] conducted a centrifugal model test to examine crack development in expansive soil slopes under wetting-drying cycles.The results indicated that crack extension increased with each cycle, correlating to the expansive strength of the foundation soil.However, this method has a significant size effect and struggles to accurately simulate the impact of rainfall on expansive soil slopes.To address these limitations, some researchers have constructed full-sized slope models aimed at analyzing crack patterns, formation, and mitigation in expansive soil slopes subjected to alternating drying and wetting conditions [15][16][17].The process of rainwater infiltration and its influence on crack formation is a key focus of these studies.Zhang et al. [18] conducted a comprehensive full-scale model test under conditions of rainfall and evaporation, providing detailed insights into the relationship between crack parameters and the number of wetting-drying cycles.This research emphasized the role of matrix suction in crack formation and concluded that composite cracks significantly contribute to the destabilization of slope.The aforementioned work has laid a solid foundation for further research.
At present, numerical simulation methods are widely applied for studying the stability of expansive soil slopes.Continuum methods, such as the finite element method (FEM) [19,20], finite difference method (FDM) [21], material point method (MPM) [22][23][24], and numerical manifold method (NMM) [25], have demonstrated their effectiveness in analyzing slope stability.Zhou et al. [26] observed that wetting-drying cycles significantly reduce the shear strength and safety factor of expansive soil slopes by altering soil structure and pore water dynamics.However, the finite element analysis method is not suitable for dealing with the problems of crack simulations in expansive soil under wetting/drying cycles.The discrete element method has been accepted by many researchers and used for slope stability studies [27][28][29].Ma et al. [30] implemented moisture migration in soil through thermal coupling and calibrated the swelling-shrinkage coefficients through laboratory tests to simulate the cracking and shrinking of expansive soil using discrete element simulation.Chen et al. [31] investigated the development of cracks in expansive soil slopes and their impact on slope stability under the wetting-drying cycle through large-scale indoor slope tests and numerical simulations.However, the simulations only examined the effect of different initial numbers of cracks on slope stability and could not reveal the pattern of crack development.Current studies do not consider the evolutionary nature of crack formation and its cumulative effect on the stability of expansive soil slopes under periodic environmental conditions.This dynamic process is crucial for understanding all the factors contributing to slope failure and remains an area for further exploration.Therefore, this study investigates the stability of expansive soil slopes with emphasis on the evolution of cracks induced by the wetting/drying cycles through the discrete element method.
This study proposed a discrete element model for analyzing the stability of the expansive soil slopes.This model coupled the unsaturated seepage simulation with the discrete element strength subtraction method to calculate the safety factors of the expansive soil slopes during wetting-drying cycles.To validate the proposed model, a discrete element simulation was conducted using the high-performance discrete element software MatDEM version 1.32 [32][33][34][35][36]. Furthermore, this discrete element method was also validated through a comparison with the limit equilibrium method and finite element strength reduction method.

Basic Discrete Element Model
In MatDEM, a soil block model is created using randomly compacted spherical elements.These particles interact through normal and tangential connections based on the linear elastic particle contact model, as depicted in Figure 1.The normal force F n and normal deformation X n are transmitted through the normal connections, which are governed by Equation (1): Water 2024, 16, 861 3 of 20 where F n = normal force; K n = normal stiffness; X n = the normal relative displacement; and X b = the breaking displacement, respectively.When X n is less than X b , the tension or pressure can exist on the connection until X n exceeds the X b , which means the normal connection breaks.After the breakage of the normal connection, only compression force works.
The shear force is transmitted through the tangential connection mechanism. ( where = tangential force; = tangential stiffness; and = the tangential displacement, respectively.
The failure criterion for a connection in the tangential direction is based on the Mohr-Coulomb criterion:

−
(3) where = the maximum shear force; = initial shear resistance; and = interparticle friction coefficient, respectively.When exceeds , the tangential connection breaks and only sliding friction exists between particles.In MatDEM, five macro parameters are used, including Young's modulus E, Poisson's ratio V, tensile strength , compressive strength , and internal friction coefficient .The intergranular mesoscopic parameters, such as Normal stiffness , tangential stiffness , normal fracture displacement , initial shear resistance , and particle friction coefficient , can be calculated using the derivation formulae in reference [37]: The shear force F s is transmitted through the tangential connection mechanism.
where F s = tangential force; K s = tangential stiffness; and X s = the tangential displacement, respectively.
The failure criterion for a connection in the tangential direction is based on the Mohr-Coulomb criterion: where F smax = the maximum shear force; F s0 = initial shear resistance; and µ p = interparticle friction coefficient, respectively.When F s exceeds F smax , the tangential connection breaks and only sliding friction exists between particles.
In MatDEM, five macro parameters are used, including Young's modulus E, Poisson's ratio V, tensile strength T u , compressive strength C u , and internal friction coefficient µ i .The intergranular mesoscopic parameters, such as Normal stiffness K n , tangential stiffness K s , normal fracture displacement X b , initial shear resistance F s0 , and particle friction coefficient µ p , can be calculated using the derivation formulae in reference [37]: where d = particle unit diameter.The closest packing, which was used in the above macroscopical transformation formula, was not the same as the input macroscopical parameters.Therefore, it is necessary to calibrate the microscopical parameters again in practice.
The damping setting method in MatDEM differs from that of PFC.In MatDEM, global damping is utilized to dissipate kinetic energy, as shown in Equation ( 9): Water 2024, 16, 861 4 of 20 where F v = damping force; η = damping coefficient; and v = velocity, respectively.The damping force is proportional to the particle velocity, but it is inversely proportional to the direction of the velocity.In MatDEM, the optimal damping coefficient is calculated using the semi-empirical formula derived by Liu Chun et al. [33] with Equation (10).This allows the kinetic energy of the system to be dissipated rapidly, leading to improved computational efficiency and enabling quasi-static simulation.
where m = particle mass; k = particle stiffness coefficient; D = particle diameter; and V = the volume of the model, respectively.
In the DEM, an essential assumption is that the velocity and acceleration of each element remain constant within very small time steps.This fundamental concept forms the basis for employing an iterative algorithm, which adheres to Newton's law of motion.Through this algorithm, the resultant force-comprising connection forces, damping forces, and gravity forces in all directions-is calculated.These forces enable the dynamic simulation of discrete elements within the system.

Discrete Element Simulation Method of Crack Evolution
Previous researchers have investigated the mechanism of crack evolutions in expansive soil [31,38,39].Generally, the tensile stress develops as the water content of the soil decreases, leading to increased suction and volume reduction.Cracks occur when the tensile stress within the soil exceeds its tensile strength and lateral pressure.Water infiltration into expansive soil can cause cracks to shrink or close completely.
In the discrete element method (DEM), each particle unit is represented as a combination of soil particles and pores, as shown in Figure 2.This approach enables the monitoring of changes in particle unit volume, which reveals the swelling and shrinkage processes of expansive soil.For two-dimensional particle elements, the equation is provided in Equation (11).R = αR 0 (11) where R = particle radius and R 0 = initial particle radius; α is the expansion/contraction coefficient, which can be calculated according to Equations (12) or (13): α = e + 1 e 0 + 1 (13) where ε v = the volumetric strain from the initial water content to the target water content; e = the pore ratios corresponding to the target water content; and e 0 = the pore ratios corresponding to the initial water content, respectively.During discrete element simulation, α was calibrated by the results of laboratory swelling contraction tests.Based on the above method, the crack development process of the expansive soil slope is simulated.Firstly, the unsaturated seepage of the expansive soil slope under atmospheric action is calculated by the finite element method, and the change in water content of the slope during the evaporation-rainfall process is obtained.Secondly, based on the finite element calculation results, the discrete element method is used to simulate the soil swelling and drying shrinkage and crack development process of the expansive soil slope.The specific steps are as follows: (1) The finite element method was used to simulate the unsaturated seepage analysis of expansive soil slope soil under the action of evaporation-rainfall, and the water content of each part of the soil at different times was obtained by taking evaporation for 14 days and then rainfall for 1 day as a wetting-drying cycle.The two-dimensional seepage partial differential governing equation is as follows: where H = the total head; k x = and the permeability coefficients in x direction; k y = the permeability coefficients in y direction; Q = the boundary flow; m w = the slope of soil water characteristic curve; γ w = the bulk density of water; and t = time, respectively.
Water 2024, 16, x FOR PEER REVIEW 5 of 21 Based on the above method, the crack development process of the expansive soil slope is simulated.Firstly, the unsaturated seepage of the expansive soil slope under atmospheric action is calculated by the finite element method, and the change in water content of the slope during the evaporation-rainfall process is obtained.Secondly, based on the finite element calculation results, the discrete element method is used to simulate the soil swelling and drying shrinkage and crack development process of the expansive soil slope.The specific steps are as follows: (1) The finite element method was used to simulate the unsaturated seepage analysis of expansive soil slope soil under the action of evaporation-rainfall, and the water content of each part of the soil at different times was obtained by taking evaporation for 14 days and then rainfall for 1 day as a wetting-drying cycle.The two-dimensional seepage partial differential governing equation is as follows: where H = the total head; = and the permeability coefficients in x direction; = the permeability coefficients in y direction; Q = the boundary flow; = the slope of soil water characteristic curve; = the bulk density of water; and t = time, respectively.The calculation model is shown in Figure A1.Considering the intense exchange of water between the surface soil and the atmosphere, the grid elements within 1m of the surface are refined, and a total of 1780 units are obtained.The initial water content of the slope soil is 20%.The soil-water characteristic curve (SWCC) and saturated permeability coefficient of the soil are shown in Figures A2 and A3.According to the measured soilwater characteristic curve, the pore water pressure in the slope is set to −387.9 kPa, without considering the effect of groundwater.The boundary conditions are set as follows: the left, right, and bottom of the slope are impermeable boundaries (total flow Q = 0), and the surface of the slope (green line part) is set as the surface atmospheric interaction boundary.
(2) The same slope model is established in the discrete element software MatDEM 2.0, and the moisture content parameters are set for the particles.According to the coordinates of the particle center point, the moisture content at the corresponding position of the finite element calculation result is obtained as the average moisture content of the particle unit.In the process of particle expansion and contraction, it is assumed that each particle remains unchanged at the corresponding position in the finite element model.
(3) Assuming that the expansion-shrinkage coefficient of particles in each stage remains unchanged, the expansion-shrinkage coefficient of each particle in this stage is cal- The calculation model is shown in Figure A1.Considering the intense exchange of water between the surface soil and the atmosphere, the grid elements within 1 m of the surface are refined, and a total of 1780 units are obtained.The initial water content of the slope soil is 20%.The soil-water characteristic curve (SWCC) and saturated permeability coefficient of the soil are shown in Figures A2 and A3.According to the measured soilwater characteristic curve, the pore water pressure in the slope is set to −387.9 kPa, without considering the effect of groundwater.The boundary conditions are set as follows: the left, right, and bottom of the slope are impermeable boundaries (total flow Q = 0), and the surface of the slope (green line part) is set as the surface atmospheric interaction boundary.
(2) The same slope model is established in the discrete element software MatDEM 2.0, and the moisture content parameters are set for the particles.According to the coordinates of the particle center point, the moisture content at the corresponding position of the finite element calculation result is obtained as the average moisture content of the particle unit.In the process of particle expansion and contraction, it is assumed that each particle remains unchanged at the corresponding position in the finite element model.
(3) Assuming that the expansion-shrinkage coefficient of particles in each stage remains unchanged, the expansion-shrinkage coefficient of each particle in this stage is calculated based on laboratory tests.Taking the adjacent time as a stage, according to the change in particle moisture content in each stage, the discrete element method is used to simulate the evolution process of soil cracks.

Principle and Method of Discrete Element Slope Stability Calculation
The model applied the discrete element strength reduction method to analyze the stability of slopes made of expansive soil, considering the development of cracks during wetting-drying cycles.Previous researchers [40,41] have highlighted the crucial role of connection strength and particle friction coefficients in determining the stability of discrete element slopes.
Consequently, the normal breaking force b f (normal breaking displacement x b ), initial shear resistance f s0 , and particle friction coefficient µ p in the linear elastic model are systematically reduced using a consistent reduction coefficient.For this purpose, critical values are defined for normal fracture force ( b) f cr , critical initial shear resistance ( f s0 ) cr , Water 2024, 16, 861 6 of 20 and critical friction coefficient µ p cr .When the conditions b f ≤ b f cr , f s0 ≤ ( f s0 ) cr , and µ p ≤ µ p cr are met, it indicates that the slope displacement surpasses the critical displacement, signifying slope instability.Conversely, when b f > b f cr , f s0 > ( f s0 ) cr , µ p > µ p cr , the slope displacement remains below the critical displacement, indicating slope stability.
The assessment of slope stability relies on the relationship between displacement at the characteristic position and critical displacement.Therefore, the reduction factor for slope instability serves as the safety factor of the slope with Equation (15).
where FOS = the safety factor; b f cr = the critical normal fracture force of particles; ( f s0 ) cr = the critical initial shear resistance of particles; and µ p cr = the critical friction coefficient of particles, respectively.Based on the secondary development functionality of MatDEM, the function of automatic calculation of slope safety factor is programmed.The slope safety factor can be calculated automatically by inputting the calculation range of the safety factor.In the discrete element slope model, five monitoring points have been strategically positioned to track the displacement of various sections of the slope during iterative calculations.The layout positions of these monitoring points are depicted in Figure 3a.
crete element slopes.Consequently, the normal breaking force (normal breaking displacement ), initial shear resistance , and particle friction coefficient in the linear elastic model are systematically reduced using a consistent reduction coefficient.For this purpose, critical values are defined for normal fracture force ( , critical initial shear resistance , and critical friction coefficient .When the conditions , , and are met, it indicates that the slope displacement surpasses the critical displacement, signifying slope instability.Conversely, when , , , the slope displacement remains below the critical displacement, indicating slope stability. The assessment of slope stability relies on the relationship between displacement at the characteristic position and critical displacement.Therefore, the reduction factor for slope instability serves as the safety factor of the slope with Equation (15).

Initial Model Establishment
The average particle radius was set as 1 cm with a particle size distribution coefficient of 0.2, weighing both the accuracy and efficiency of calculations.This coefficient implies that the particle radius follows a normal distribution ranging from 0.836 cm to 1.204 cm.The discrete element slope accumulation model, depicted in Figure 3b, consists of a total of 294,183 particles.

Initial Model Establishment
The average particle radius was set as 1 cm with a particle size distribution coefficient of 0.2, weighing both the accuracy and efficiency of calculations.This coefficient implies that the particle radius follows a normal distribution ranging from 0.836 cm to 1.204 cm.The discrete element slope accumulation model, depicted in Figure 3b, consists of a total of 294,183 particles.

Material Setting and Parameter Determination
The particle mesoscopic parameters of the discrete element slope model were calibrated using the biaxial compression test, taking into account the physical and mechanical properties and the macroscopic parameters of Nanyang expansive soil, as shown in Tables 1 and 2, respectively, after subjecting it to wetting-drying cycles.The mineral composition of Nanyang expansive soil is illustrated in Figure A4 and Tables A1 and A2, and the particle grading curve is shown in Figure A5.To enhance the analysis, the following assumptions were made.(1) The elastic modulus, Poisson's ratio, and tensile strength of expansive soil are invariant, irrespective of variations in water content and the number of wetting-drying cycles.
(2) The normal stiffness k n , the tangential stiffness k s , and the normal fracture displacement x b (the normal fracture force F n ) were calculated using the macro and micro transformation equation.Nevertheless, the initial shear resistance and particle friction coefficient necessitate calibration.
The calibration process yielded particles of the same size and distribution in the biaxial compression test as in the slope model.The biaxial compression model had dimensions of 1 m in width and 2 m in height, consisting of a total of 5202 particle units.The confining pressures applied during the biaxial compression test matched the laboratory test conditions, set as 50, 100, 150, and 200 kPa, respectively, in accordance with the laboratory test conditions.The mesoscopic parameters of the particle obtained from the final calibration are presented in Table 3.

The Effects of Crack on the Stability of Expansive Soil Slopes
The discrete element model presented in this study examines the impact of cracks on the stability of expansive soil slopes, focused on two crucial factors: (1) Slope model with cracks: The discrete element model accurately simulated the development of cracks during wetting-drying cycles.During the rainfall infiltration process, certain cracks may partially close as soil particles expand, while others remain open.
(2) Strength parameters for soil with cracks: In order to consider the effect of cracks on the simulation of unsaturated seepage, slopes made of expansive soil are divided into two zones: crack zone and no-crack zone.This division is based on different depths of cracks, as shown in Figure 4. Varying permeability coefficients are established within these zones.Additionally, the particle characteristics of the soil in the crack zone are established by analyzing the shear strength properties of samples that have undergone various cycles of drying and wetting.In contrast, the strength characteristics of the soil in the no-crack zone are not influenced by the process of drying and wetting.
ious cycles of drying and wetting.In contrast, the strength characteristics of the soil in the no-crack zone are not influenced by the process of drying and wetting.
It is worth noting that the depth of the crack area is determined based on the average maximum crack depth observed from the expansive soil slopes throughout various wetting-drying cycles, as presented in Table 4.

Failure Model
The displacement diagrams (Figures 5 -10) show that the displacement changes in the expansive soil slopes over various wetting-drying cycles.Figure 5 shows the displacement diagram without any wetting-drying cycle.After 25,000 time steps, the slope remains stable, but a noticeable crack zone appears at the base of the slope.At 50,000 time steps, the sliding surface forms in the initial section, causing a portion of the soil to detach from the slope.Furthermore, at 75,000 time steps, the sliding distance of the first section exceeds 1 m, which suggests the progression of the landslide into a deeper stage.It is worth noting that the depth of the crack area is determined based on the average maximum crack depth observed from the expansive soil slopes throughout various wettingdrying cycles, as presented in Table 4.

Failure Model
The displacement diagrams (Figures 5 -10) show that the displacement changes in the expansive soil slopes over various wetting-drying cycles.Figure 5 shows the displacement diagram without any wetting-drying cycle.After 25,000 time steps, the slope remains stable, but a noticeable crack zone appears at the base of the slope.At 50,000 time steps, the sliding surface forms in the initial section, causing a portion of the soil to detach from the slope.Furthermore, at 75,000 time steps, the sliding distance of the first section exceeds 1 m, which suggests the progression of the landslide into a deeper stage.
The behavior of expansive soil slopes after one and two cycles of wetting and drying, is depicted in Figures 6 and 7, respectively.At 25,000 time steps, there is no clear evidence of any slippage in these slopes.Nevertheless, cracks emerge and a sliding surface forms at the higher section of the slope, progressively expanding in a circular manner within the slope.After 50,000 time steps, the sliding crack plane achieves complete connectivity, which signifies a profound circular slide.
In the early stage of the wetting-drying cycle, although the shear strength of expansive soil decreases rapidly, the number of slope cracks and the depth of crack development is still small, resulting in a similar pattern to that without any wetting-drying cycles.Due to the expansion and contraction cracks caused by atmospheric action in all parts of the slope after the wetting-drying cycle, the slope soil mass will continue to develop along the existing cracks at the top of the slope under the action of gravity, and the slope sliding surface will first develop into the crack at the top of the slope, and finally form a circular arc overall sliding.
The expansive soil slopes after three to five wetting-drying cycles exhibit similar patterns of shallow failure, as shown in Figures 8-10.At 25,000 time steps, shallow cracks contribute to the formation of a sliding surface.By 50,000 time steps, this sliding crack plane becomes fully connected, resulting in a shallow landslide.The behavior of expansive soil slopes after one and two cycles of wetting and drying, is depicted in Figures 6 and 7, respectively.At 25,000 time steps, there is no clear evidence of any slippage in these slopes.Nevertheless, cracks emerge and a sliding surface forms at the higher section of the slope, progressively expanding in a circular manner within the slope.After 50,000 time steps, the sliding crack plane achieves complete connectivity, which signifies a profound circular slide.
In the early stage of the wetting-drying cycle, although the shear strength of expansive soil decreases rapidly, the number of slope cracks and the depth of crack development is still small, resulting in a similar pattern to that without any wetting-drying cycles.Due to the expansion and contraction cracks caused by atmospheric action in all parts of the slope after the wetting-drying cycle, the slope soil mass will continue to develop along the existing cracks at the top of the slope under the action of gravity, and the slope sliding surface will first develop into the crack at the top of the slope, and finally form a circular arc overall sliding.The expansive soil slopes after three to five wetting-drying cycles exhibit similar patterns of shallow failure, as shown in Figures 8-10.At 25,000 time steps, shallow cracks contribute to the formation of a sliding surface.By 50,000 time steps, this sliding crack plane becomes fully connected, resulting in a shallow landslide.
During the later stage of the wetting-drying cycle (3-5 times), the shear strength of the expansive soil decreases, increasing the number and depth of cracks in the slope.The lower soil strength and increased expansion and contraction cracks on the surface of the slope are the main causes of shallow failure.There are differences in the morphology of shallow failure on expansive soil slopes after 3-5 wetting-drying cycles: (1) There are variations in the starting position of the sliding surface.( 2) As the number of wetting-drying cycles increases, the depth of the shallow sliding surface decreases gradually, and the volume of the sliding body also decreases.This is because in the later stage of wetting-drying cycle, the difference in soil strength parameters is not large, and the increase in the number of surface expansion and contraction cracks is the main reason for the above phenomenon.As the number of wetting-drying cycles increases, the number of cracks on the slope surface also increases.Most of the newly formed cracks are shallow in depth.The increased number of cracks divides the surface layer of the slope into smaller soil blocks, further reducing the integrity of the soil.The slope sliding will initially occur in the area with more surface cracks, i.e., weak parts.Therefore, the depth of the slope sliding surface gradually decreases.During the later stage of the wetting-drying cycle (3-5 times), the shear strength of the expansive soil decreases, increasing the number and depth of cracks in the slope.The lower soil strength and increased expansion and contraction cracks on the surface of the slope are the main causes of shallow failure.There are differences in the morphology of shallow failure on expansive soil slopes after 3-5 wetting-drying cycles: (1) There are variations in the starting position of the sliding surface.(2) As the number of wetting-drying cycles increases, the depth of the shallow sliding surface decreases gradually, and the volume of the sliding body also decreases.This is because in the later stage of wetting-drying cycle, the difference in soil strength parameters is not large, and the increase in the number of surface expansion and contraction cracks is the main reason for the above phenomenon.As the number of wetting-drying cycles increases, the number of cracks on the slope surface also increases.Most of the newly formed cracks are shallow in depth.The increased number of cracks divides the surface layer of the slope into smaller soil blocks, further reducing the integrity of the soil.The slope sliding will initially occur in the area with more surface cracks, i.e., weak parts.Therefore, the depth of the slope sliding surface gradually decreases.

Slope Stability
The analysis shows that the displacements at the measuring points, evaluated after 100,000 iterations, accurately indicate the instability of the slope in relation to the corresponding reduction factor.Hence, the number of iterations following each reduction is set as 100,000.Figure 11 illustrates the displacement of individual measuring points on the slope of the expansive soil at different stages throughout the wetting-drying cycle.The data indicate a significant increase in in the displacement at each measurement point when the strength reduction coefficient approaches the critical value.The safety factor of the slope can be calculated using these displacement data.
Table 5 provides a comprehensive overview of the safety factors associated with expansive soil slopes under various wetting-drying cycles.Initially, the safety factor of the expansive soil slopes without any wetting-drying cycle is about 3.30.After the first wetting-drying cycle, it decreases to 3.22.The safety factor decreases to 2.91 after the second cycle, and subsequently drops to 2.23 after three cycles.These findings suggest that the safety factor decreases at an increasing rate as the number of wetting-drying cycles increases.After undergoing four to five cycles, the decrease in the safety factor becomes less significant, with respective reductions of 0.23 and 0.38.During the initial wetting-drying process, a noticeable decrease in the strength of the expansive soil is seen.Simultaneously, a large quantity of cracks and crack connections take place, which leads to a faster deterioration in the slope's safety factor.However, after three to five cycles of wetting and drying, the shear strength of the soil begins to reach a stable value, with a decreasing rate of crack development.Consequently, the safety factor of the slope decreases at a lower rate.During the initial wetting-drying process, a noticeable decrease in the strength of the expansive soil is seen.Simultaneously, a large quantity of cracks and crack connections take place, which leads to a faster deterioration in the slope's safety factor.However, after three to five cycles of wetting and drying, the shear strength of the soil begins to reach a stable value, with a decreasing rate of crack development.Consequently, the safety factor of the slope decreases at a lower rate.The sliding patterns of expansive soil slopes during wetting-drying cycles were analyzed using three different methods: the limit equilibrium method, the finite element strength reduction method, and the discrete element method (DEM).The analyses pre-  The sliding patterns of expansive soil slopes during wetting-drying cycles were analyzed using three different methods: the limit equilibrium method, the finite element strength reduction method, and the discrete element method (DEM).The analyses presented in Figures 12-15 offer valuable insights into the behavior of soil slopes under varying conditions.
During the initial 0 to 1 wetting-drying cycles, all three methods indicated a uniform occurrence of profound circular sliding during the initial wetting-drying cycles.The small number of wetting-drying cycles has resulted in a limited number and depth of cracks in the DEM.As a result, the integrity of the shallow soil of the slope remains high.The presence of these cracks has little effect on the safety factor or sliding surface of the slope.After three wetting-drying cycles, the DEM simulations revealed an obvious occurrence of cracks in the slope, which had a substantial impact on the properties of the sliding surface.More precisely, the DEM showed a relatively shallow sliding pattern, which differed from the deep circular sliding pattern that was still evident in the limit equilibrium and finite element strength reduction methods.
After five wetting-drying cycles, shallow landslides occurred in all three methods.This highlights the pronounced effect of multiple wetting-drying cycles on the stability and failure mechanisms of expansive soil slopes.

The Slope Safety Factor
The safety factors of expansive soil slopes, as obtained from different methods, are presented in Figure 16 and summarized in Table 6.The DEM yields a smaller safety factor, decreasing from 3.30 to 1.68.trends are observed with the limit equilibrium method and the finite element strength reduction method, starting with initial values of 3.385 and 3.475, respectively.With an increase in the number of wetting-drying cycles, the safety coefficient gradually decreases with minimal variations.Specifically, the results of the limit equilibrium method decrease by 0.459, from 3.385 to 2.926, while the results of the finite element method decrease by 0.459, from 3.475 to 2.978.
Water 2024, 16, x FOR PEER REVIEW After five wetting-drying cycles, shallow landslides occurred in all three m This highlights the pronounced effect of multiple wetting-drying cycles on the and failure mechanisms of expansive soil slopes.

The Slope Safety Factor
The safety factors of expansive soil slopes, as obtained from different meth presented in Figure 16 and summarized in Table 6.The DEM yields a smaller safet decreasing from 3.30 to 1.68.Similar trends are observed with the limit equ method and the finite element strength reduction method, starting with initial v 3.385 and 3.475, respectively.With an increase in the number of wetting-drying cy safety coefficient gradually decreases with minimal variations.Specifically, the re the limit equilibrium method decrease by 0.459, from 3.385 to 2.926, while the re the finite element method decrease by 0.459, from 3.475 to 2.978.In the case of the expansive soil slopes without any wetting-drying cycles, th safety factors obtained by the DEM were similar to those obtained by the limit equi method and finite element method.Note that the model parameters in the DEM brated through biaxial tests.The macro and micro transformation formula in the m training function is specifically applicable to the compact regular packing model roidal elements.Randomly packed particles typically exhibit macro mechanical ties that are lower than the values predicted by the other two methods.Addition calculation results of the DEM are still comparable to those of the other method  In the case of the expansive soil slopes without any wetting-drying cycles, the slope safety factors obtained by the DEM were similar to those obtained by the limit equilibrium method and finite element method.Note that the model parameters in the DEM are calibrated through biaxial tests.The macro and micro transformation formula in the material training function is specifically applicable to the compact regular packing model of spheroidal elements.Randomly packed particles typically exhibit macro mechanical properties that are lower than the values predicted by the other two methods.Additionally, the calculation results of the DEM are still comparable to those of the other methods, with only a negligible discrepancy.This validates the effectiveness of the discrete element strength reduction technique in computing slope safety factors.
Comparatively, in the case of the expansive soil slopes undergoing several cycles of wetting-drying cycles, similar results were obtained by the limit equilibrium method and finite element strength reduction method: the safety factor slightly decreased as the increasing number of wetting-drying cycles did.However, the reduction in soil strength was only considered in the crack zone.Conversely, the DEM took into account both the reduction in soil strength and the evolution of cracks in the soil slope during the drying and wetting cycles.Therefore, the safety factor determined by the DEM is considerably lower than that of the other two methods.The results indicate that when dealing with the stability analysis of expansive soil slopes, it is essential to study the impact of swelling and shrinking cracks on slope stability.

Conclusions
This study investigated the crack evolution in the expansive soil slopes during the wetting-drying cycles by the discrete element method.Coupled with the strength discount method, the stability of expansive soil slopes under climate changes was computed in the discrete element model.The following results can be obtained: 1.
The slope stability can be obtained by the discrete element strength reduction method.
Compared with the calculation results of the limit equilibrium method and finite element method, a much lower stability was computed in the DEM considering the crack evolutions in the expansive soil slopes, and the effectiveness of the discrete element strength reduction method is verified.2.
The safety factor of an expansive soil slope, calculated by the discrete element strength reduction method, gradually decreases with an increase in the number of wettingdrying cycles.The decrease rate of the slope safety factor increases gradually when the number of wetting-drying cycles is between 0 and 3.However, during three to five wetting-drying cycles, the reduction rate of the slope safety factor gradually decreases and tends to stabilize.The trend in variation in the safety factor, as calculated by the DEM, differs significantly from that of the limit equilibrium method and the finite element method.As the number of wetting-drying cycles increases, the calculation results of the discrete element method are much smaller than those of the other two methods.

3.
Following 0-2 wetting-drying cycles, the discrete element slope reaches a critical state, presenting a deep arc-shaped sliding surface.The sliding surface initially develops at the foot of the slope without wetting-drying cycles, and after 1-2 wetting-drying cycles, it begins with the expansion and contraction cracks at the top of the slope.After 3-5 wetting-drying cycles, the discrete element slope reaches a critical state of shallow sliding.As the number of wetting-drying cycles increases, the volume of shallow sliding soil gradually decreases, and the depth of the sliding surface decreases accordingly.There is a significant discrepancy between the sliding surface calculated using the DEM and the results obtained from the limit equilibrium method and the finite element method.

4.
The study put emphasis on the evolution of cracks in the stability analysis of expansive soil slopes.Important factors were considered: not only the reduction in soil strength but also the presence, number, depth, and distribution of cracks.These factors greatly influence the stability of expansive soil slopes.

Figure 2 .
Figure 2. Discrete element particles simulate the swelling and shrinkage of expansive soil in the wetting-drying cycle.

Figure 2 .
Figure 2. Discrete element particles simulate the swelling and shrinkage of expansive soil in the wetting-drying cycle.

Figure 3 .
Figure 3. (a) Measuring points layout of the discrete element slope model; (b) discrete element slope model.

Figure 3 .
Figure 3. (a) Measuring points layout of the discrete element slope model; (b) discrete element slope model.

Figure 4 .
Figure 4. Division of crack area of expansive soil slopes.

Figure 4 .
Figure 4. Division of crack area of expansive soil slopes.

Figure 5 .
Figure 5. Displacement diagrams of expansive soil slopes without wetting-drying cycle: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 6 .Figure 7 .
Figure 6.Displacement diagrams of expansive soil slopes after one cycle of wetting-drying: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 7 .
Figure 7. Displacement diagrams of expansive soil slopes after two cycles of wetting-drying: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 8 .
Figure 8. Displacement diagrams of expansive soil slopes after three cycles of wetting-drying: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 8 .
Figure 8. Displacement diagrams of expansive soil slopes after three cycles of wetting-drying: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 9 .Figure 10 .
Figure 9. Displacement diagrams of expansive soil slopes after four cycles of wetting-drying: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 10 .
Figure 10.Displacement diagrams of expansive soil slopes after five cycles of wetting-drying: (a) at 25,000 time steps, (b) at 50,000 time steps, (c) at 75,000 time steps, (d) at 100,000 time steps.

Figure 11 .
Figure 11.Relationship between displacement and reduction coefficient of slope measuring points after different times of wetting-drying cycles: (a) without wetting-drying cycle, (b) the first cycle of wetting-drying, (c) the second cycle of wetting-drying, (d) the third cycle of wetting-drying, (e) the fourth cycle of wetting-drying, (f) the fifth cycle of wetting-drying.

Figure 11 .
Figure 11.Relationship between displacement and reduction coefficient of slope measuring points after different times of wetting-drying cycles: (a) without wetting-drying cycle, (b) the first cycle of wetting-drying, (c) the second cycle of wetting-drying, (d) the third cycle of wetting-drying, (e) the fourth cycle of wetting-drying, (f) the fifth cycle of wetting-drying.

4. 3 .
Comparison between the Limit Equilibrium Method and the Finite Element Method 4.3.1.The Slope Failure Pattern

Figure 13 .
Figure 13.Position of sliding surface of expansive soil slopes after one wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 14 .
Figure 14.Position of sliding surface of expansive soil slopes after three wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 13 .
Figure 13.Position of sliding surface of expansive soil slopes after one wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 14 .
Figure 14.Position of sliding surface of expansive soil slopes after three wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 13 .
Figure 13.Position of sliding surface of expansive soil slopes after one wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 14 .
Figure 14.Position of sliding surface of expansive soil slopes after three wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 15 .
Figure 15.Position of sliding surface of expansive soil slopes after five wetting-drying cycles: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 13 .
Figure 13.Position of sliding surface of expansive soil slopes after one wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 14 .
Figure 14.Position of sliding surface of expansive soil slopes after three wetting-drying cycle: (a) DEM, (b) limit equilibrium method, (c) finite element method.

Figure 16 .
Figure 16.Variation trend of safety factors of expansive soil slopes after different times of drying cycles.

Figure 16 .
Figure 16.Variation trend of safety factors of expansive soil slopes after different times of wettingdrying cycles.
is montmorillonite, It is illite, Kao is kaolinite, C is chlorite, I/S is illite-montmo layer, C/S is green-montmorillonite mixed layer.

Figure A5 .
Figure A5.Grain size distribution curve of Nanyang expansive soil.

Figure A4 .
Figure A4.X-ray diffraction results of Nanyang expansive soil.
Figure A4.X-ray diffraction results of Nanyang expansive soil.
is montmorillonite, It is illite, Kao is kaolinite, C is chlorite, I/S is illite-montmo layer, C/S is green-montmorillonite mixed layer.

Figure A5 .
Figure A5.Grain size distribution curve of Nanyang expansive soil.

Figure A5 .
Figure A5.Grain size distribution curve of Nanyang expansive soil.

Table 1 .
Basic physical properties of Nanyang expansive soil.

Table 2 .
Macroscopic mechanical parameters of Nanyang expansive soil.

Table 3 .
Granular meso-parameters of slope model.

Table 4 .
Depth of crack area in the model.

Table 4 .
Depth of crack area in the model.

Table 5 .
Calculation results of the safety factor of expansive soil slope in DEM.

Table 5 .
Calculation results of the safety factor of expansive soil slope in DEM.

Table 6 .
The values of safety factor of expansive soil slopes.

Table 6 .
The values of safety factor of expansive soil slopes.

Table A1 .
Mineral Composition of Nanyang expansive Soil.

Table A2 .
Clay minerals of Nanyang expansive soil.

Table A1 .
Mineral Composition of Nanyang expansive Soil.

Table A2 .
Clay minerals of Nanyang expansive soil.

Table A1 .
Mineral Composition of Nanyang expansive Soil.

Table A2 .
Clay minerals of Nanyang expansive soil.