1. Introduction
A large amount of statistical data show that rainfall, especially heavy rain, is a pivotal cause of geological disasters such as landslides [
1,
2]. Therefore, slope instability analysis under rainfall conditions is of great significance for construction safety and geological disaster protection. Rainwater soaking reduces soil cohesion, matrix suction, and shear strength, making it more susceptible to damage. Moreover, with rainfall infiltration, there will be erosion and migration of fine particles in granular soils, leading to an increase in porosity and corresponding permeability coefficient, which in turn leads to a further increase in pore water pressure. All of these are the causes of slope instability. Therefore, to accurately simulate the possible instability of soil slopes under rainfall infiltration conditions, it is necessary to reflect the physical coupling of the seepage, erosion, and deformation involved.
Rainfall-induced slope instability is a typical soil–pore fluid coupling problem with strong nonlinearity. The finite element model (FEM) is still the mainstream numerical method for the seepage–stress coupling problem [
3,
4,
5,
6]. For the cofferdam of tidal flats behind Changxing Island, Li established its seepage–stress coupling FEM considering the soil strength degradation with the wave cyclic loading, and explored the seepage–stress coupling properties of the cofferdam on a soft clay foundation under a storm surge attack [
3]. Xu established a seepage/stress-damage coupling FEM of a deep excavated canal in the Xichuan Section, and explored its long-term performance, including settlement and damage [
4]. Zhou used a seepage–stress coupling FEM model to study the rainfall-induced loess landslide including the hazard’s occurrence and evolution [
5]. Liu et al. [
6] established a finite element numerical model for slope under rainfall load, which realized the unidirectional coupling of the hydro-mechanical behavior of the unsaturated soil slope; then, they studied the influence of rainfall history and saturated permeability coefficient on the internal stress and pore pressure of the slope. The finite element method together with strength reductions is usually used to evaluate slope stability. The slope is judged to be unstable when the slope experiences phenomena such as sudden changes in displacement or non-convergence of the calculation [
7]. However, the instability of slopes under rainfall infiltration shows a significant nonlocal landslide morphology, that is, an inconspicuous sliding band, while the results of finite element plastic or damage model simulation show a narrow sliding band morphology, which is not consistent with reality.
When faced with discontinuous media containing crack extension, holes, etc., classical continuum mechanics will fail. To solve this problem, researchers have proposed different calculation methods such as the extended finite element (XFEM) [
8,
9] and phase field [
10,
11]. The XFEM avoids the mesh meshing continuously to fit the crack extension by introducing the enriched shape functions reflecting the displacement jump at the crack and the stress singularity at the crack tip into the standard FEM displacement model. In the XFEM, the crack topology is represented implicitly by level sets, which enables the cracks to propagate completely independent of the fixed mesh. However, the large number of enriched nodal degrees of freedom, and complex crack propagation processing involving the determination of the extension direction of crack and iteratively converging the propagation length, greatly increased the computational time. For the phase field method, it is based on the variational approach of fracture mechanics and introduces phase field scalars characterizing the material damage state to simulate discontinuities in a smeared form in narrow banded regions. However, for problems involving multiple physical fields, the fracture phase field method will involve more fields, greatly increasing the difficulty of solving the problem. Different from other methods, Silling proposed peridynamics based on the concept of a nonlocal interaction. It describes material deformation by solving spatial integral equations instead of partial differential equations [
12]. It is usually solved using the meshless particle method, which not only has good applicability in simulating damage evolution and crack propagation but also has the potential to calculate problems with spatial characteristic scales or significant nonlocal effects. Because of the above advantages, many scholars have applied peridynamics to the deformation and failure analysis of rock and soil with nonlocal characteristics. Jabakhanji [
13] established a peridynamics model for coupling unsaturated soil deformation with transient water–gas flow, and verified the simulation accuracy by comparing it with experimental results. Gu et al. [
14] proposed a fluid–structure coupling scheme based on peridynamics to simulate the interaction between soil and pore fluid in the liquefaction analysis of saturated granular soil, and used this method to analyze the liquefaction process of liquefiable soil on the horizontal and inclined ground under seismic action. Zhou et al. [
15] established a one-dimensional saturated soil frost heave model based on peridynamics and accelerated the calculation of convolution in the model through a fast Fourier transform. This method greatly improved the calculation efficiency while ensuring the calculation accuracy. Menon and Song [
16] established an updated Lagrangian peridynamics model to simulate large deformations of unsaturated soils under drainage conditions. The model introduced the concept of “bond-related secondary neighborhood”, which greatly enhanced the stability and accuracy of the soil skeleton under large deformations. The effectiveness and robustness of the method were verified through numerical examples. Song and Silling [
17] used the energy balance of unsaturated porous media to derive the effective stress state and suction state of the porous media skeleton that is conjugated to the nonlocal deformation state. They established the multiphase constitutive correspondence principle between classical unsaturated seepage mechanics and peridynamics through energy equivalence, providing a new way to introduce various constitutive relations into peridynamics. Song and Hossein [
18] combined the Cosserat continuum theory with peridynamics and introduced the length scale related to the microstructure as a nonlocal parameter in the model to simulate the shear band bifurcation and cracking phenomena in dry porous media. Menon and Song [
19] proposed a nonlocal fluid model for state-based peridynamics, introduced the classical local constitutive model of the solid skeleton and the generalized unsaturated Darcy law into peridynamics, and simulated the shear band of unsaturated soil. Liu et al. [
20] proposed an improved bond-based peridynamics moisture–mechanical coupling model to simulate the shrinkage cracking of soil and performed a detailed simulation of the moisture migration and deformation of the soil ring and soil strip. Ren et al. [
21] proposed a peridynamics-smoothed particle hydrodynamics coupling strategy to simulate the damage of buried explosions to soil. Liu et al. [
22] used fully coupled hygro-mechanical ordinary state-based peridynamics to establish a model of soil strip desiccation deformation and curling, successfully capturing the entire curling process and exploring the influencing factors of the liquid limit, additional evaporation surface, and thickness on the curling performance. Artificial viscosity and virtual particle methods were used in the simulation to improve the calculation accuracy and eliminate the numerical instability caused by shock wave propagation. Sedighi et al. [
23] proposed a nonlocal clay erosion formula based on peridynamics, integrating clay swelling, particle separation, and transport into one model. The final simulation results were in good agreement with the experimental data.
In slope instability simulation, even if the slope enters a critical failure state, peridynamics can continue the calculation, which is of great significance for analyzing the failure process of the slope. Lai et al. [
24] applied peridynamics to slope stability analysis and used non-ordinary state-based peridynamics combined with the Drucker–Prager constitutive model to model a failure process of a two-dimensional soil slope. The results were consistent with the finite element method and the dynamic sliding process of slope failure was given. Zhang and Zhang [
25] used ordinary state-based peridynamics combined with the Drucker–Prager yield criterion to successfully simulate the localized deformation of the slope and locate the critical failure surface of the slope. This method successfully avoided the phenomenon of non-physical plastic property changes under extreme non-uniform deformation. Wang et al. [
26] combined non-ordinary state-based peridynamics with random field models to study the influence of spatial variability of soil parameters on slope strength. They also compared the effects of pulse and non-pulse earthquakes on the landslide process. Zhou et al. [
27] used peridynamics to research the relationship between shear bandwidth and horizon size, and combined non-ordinary state-based peridynamics with the strength reduction method to conduct safety analysis on intact slopes and slopes with defects, respectively, providing guidance for real engineering. Although peridynamics has many advantages over the traditional finite element method, since it considers nonlocal effects, it requires more computing resources than traditional methods during discretization, which greatly prolongs the calculation time. It also limits the application of peridynamics in practical problems. To solve this problem, many scholars have conducted research. Yang and Liu [
28] established a multiscale method based on the combination of a boundary element method and peridynamics. By introducing boundary element nodes at the same position as the PD material points on the interface, a shared node coupling model was established to give full play to the advantages of both methods. While accurately simulating object deformation and crack propagation, higher computational efficiency was achieved. Ni et al. [
29] used PD and FEM coupling to simulate the propagation of hydraulic fractures in saturated porous media. Liu et al. [
30] proposed an element-based coupling method of peridynamics and the finite element method. Sun and Fish [
31] used non-ordinary state-based peridynamics coupled with the finite element method combined with Boit theory to simulate crack propagation in saturated porous media. Jin et al. [
32] established a coupling model of non-ordinary state-based peridynamics and the finite element method, connecting the two methods by embedding peridynamics material points in the interface elements, and proposed two coupling force distribution schemes: volume coupling and contact coupling. Although some researchers have incorporated peridynamics into the slope stability analysis, the integration of erosion effects with peridynamics remains an area underexplored. Consequently, the simulation of slope instability under the influence of rainfall could benefit from further refinement.
Fine particle erosion is one of the important factors affecting slope instability. The basic mathematical model of erosion is mainly developed based on experimental simulation and porous continuum mechanics. Dahaghi et al. [
33] established an expression for the permeability during the erosion process based on the mass conservation equation for particle transport in porous media, derived an analytical model for particle movement and deposition, and considered the effects of the inflow fluid rate and the concentration of particles it carries on the erosion process. Yang et al. [
34] divided the soil into four components: a stable skeleton composed of coarse particles, fine particles that can be eroded, liquefied particles that have been eroded, and seepage liquid. They then established a mass conservation equation to describe the migration of liquefied particles. Golay et al. [
35] studied the erosion between fluid and soil particle surfaces at the pore scale and described the evolution of the water–soil interface using a level set function. Chang [
36] conducted experiments to study the corrosion process under complex conditions, analyzed the influence of hydraulic gradient and stress state on erosion, and derived the control equation of the critical failure hydraulic gradient. Sterpi [
37] experimentally studied the erosion of fine particles in soil samples by controlled seepage and established an FEM model to analyze the effect of fine particle loss on the stress–strain distribution in the soil. Zhang et al. [
38] established an unsaturated seepage and erosion coupled model through a study of the interaction between erosion and infiltration and analyzed the impact of erosion on soil slope stability.
Considering the above statement, this paper attempts to propose a PD-FEM coupling model and a modified erosion constitutive model to simulate and analyze the slope instability caused by heavy rainfall, and explore the influence of rainfall duration, rainfall intensity, erosion amount, and initial saturated permeability coefficient on slope safety. Specifically,
Section 2 gives the basic formula of non-ordinary state-based peridynamics and establishes a numerical model of the slope.
Section 3 proposes a coupling scheme of non-ordinary state-based peridynamics and finite element method, and explores the effects of rainfall duration and rainfall intensity on slope stability.
Section 4 further considers the erosion of soil, proposes a modified erosion constitutive model, and explores the effects of erosion and initial saturated permeability on slope stability. Finally,
Section 5 summarizes the main conclusions of this paper.
3. Slope Stability Analysis under Rainfall Conditions with NOSBPD-FEM Modeling
3.1. Fluid–Structure Interaction Solution
Unsaturated soil is composed of three phases: soil particles, water, and air. Its physical properties will change significantly when affected by the complex external environment. Among them, water content has the strongest impact on unsaturated soil. During rainfall, rainwater falls on the slope surface and is absorbed by the surface soil particles under the action of matrix suction. When the absorbed water reaches an extreme value, the microcracks in the soil will become infiltration channels for water, resulting in unsaturated seepage, which will increase the volume moisture content of the lower soil layer and cause the pore water pressure
of the soil to increase. Since the pores in the soil are assumed to be always connected to the outside atmosphere, the pore air pressure
will not change significantly, resulting in a continuous decrease in the matrix suction
, which in turn affects physical parameters such as soil permeability coefficient, cohesion, and internal friction angle [
40].
According to the Terzaghi principle of soil mechanics, saturated soil can be regarded as consisting of two phases: soil particle skeleton and pore water. When the soil pores are filled with water, the pore water pressure will bear part of the stress, while the deformation of the soil is mainly controlled by the effective stress of the soil skeleton, and the pore water pressure does not contribute to the strength and deformation of the soil. The situation of unsaturated soil is more complicated. It is composed of three phases: soil particle skeleton, pore water, and gas in the pores. The pore gas pressure also plays a role in sharing some of the stress. To describe the effective stress of unsaturated soil, Zhao et al. [
41] proposed a specific expression for the effective stress of unsaturated soil based on the work-based effective stress principle of unsaturated soil proposed by Dean and Houlsby [
42]:
where
is the effective stress;
is the total stress;
is the saturation; and
is the second-order unit tensor.
To consider the effect of seepage on the original stress state of the soil, effective stress needs to be used to correct the peridynamic equation of motion. Decompose the Cauchy stress in Equation (6) into two parts: spherical stress and deviatoric stress, and obtain the Cauchy spherical stress term reflecting the pore water pressure, and obtain the corrected Cauchy stress. Use the corrected Cauchy stress to continue calculating the force density state.
Soil is a loose structure with many pores inside. When rainwater from outside flows in and displaces the original air in the pores, the density of the soil will increase, increasing the sliding force of the slope and harming an adverse effect on the stability of the slope. The updated expression of soil density is
where
is the dry soil density;
is the current volume moisture content of the soil;
is the density of water.
During rainfall, the infiltration water flow will also drag the soil particles, generating a seepage volume force
in the same direction as the seepage velocity, and its magnitude is proportional to the water head gradient.
where
is the soil specific weight;
is the total water head, which can be obtained by solving the seepage control equation [
43]
where
and
are the permeability coefficients in the direction x and the direction y respectively;
is the volume water content.
This paper selects the Van-Genuchten water–soil characteristic curve equation and the unsaturated permeability coefficient equation to describe the seepage characteristics of the soil. The Van-Genuchten water–soil characteristic curve expression is [
44]
where
is the saturated volume moisture content;
is the residual volume moisture content;
is the pore gas pressure;
is a parameter related to the inverse of the air intake value;
is a parameter related to the slope of the water–soil characteristic curve;
is a parameter related to the residual state.
The unsaturated permeability coefficient equation is
where
is the saturated permeability coefficient.
After substituting the seepage volume force from Equation (19), it can be added to the peridynamic equation of motion as an external body force density term.
Soil deformation can also affect seepage. After calculating the volume strain of the soil through peridynamics, the porosity after deformation can be updated according to soil mechanics as
where
is the initial porosity of the soil and
is the volume strain. The change in porosity leads to a decrease in the cross-sectional area available for water flow, and the saturated permeability of the soil decreases. The relationship between the saturated permeability and porosity is as follows:
where
is the initial saturated permeability coefficient. The change in saturated permeability coefficient will cause the change in soil permeability function, thus affecting the seepage distribution of slope.
In summary, after establishing the fluid domain and solid domain models, the pore water pressure distribution is obtained according to the initial and boundary conditions of the seepage, and the corresponding seepage volume force and density change are calculated. Considering the effect of the pore water pressure on the solid domain that causes the effective stress to decrease, the slope displacement and stress under the influence of seepage are calculated. The volume strain generated in the solid domain will cause the porosity to change, changing the permeability coefficient of the fluid domain, thereby affecting the pore water pressure distribution of the seepage field in the next time step, thus achieving the coupling between the solid domain and the fluid domain. Compared with other researchers on fluid–structure coupling [
45,
46,
47], although the coupling model in this paper adopts a nonlocal method, its purpose is to capture the nonlocal characteristics of the slope slip zone and does not consider the influence of crack damage. At the same time, the soil uses a perfect elastoplastic model, and the influence of parameters such as soil saturation is not reflected in the constitutive model, which greatly simplifies the soil model.
In the coupling process, since peridynamics uses uniform meshfree particle discretization, while the finite element method uses four-node quadrilateral isoparametric element discretization, the data storage points of the two do not correspond one to one. Usually, the distribution of peridynamics material points is denser than that of FEM nodes. Therefore, the element number and local coordinates of each material point within the element, as well as the number of material points near each node, should be calculated in advance. When the coupled calculation is started, the required data for the peridynamics material points can be obtained from the element node data using shape function interpolation. Since the distribution of material points is denser, each FEM node can be considered to be located in the area formed by the nearest four PD material points. If this area is regarded as an element, the data required by the node can still be obtained by interpolating the data of nearby material points and shape functions. Thus, a data exchange model between the deformation field and the seepage field is established.
In terms of program architecture, since the peridynamics part uses explicit iteration to solve the motion equations, the maximum time step that satisfies the iterative convergence is much smaller than the timescale of the rainfall process. It is unrealistic to calculate the entire rainfall process with the peridynamics time step. Therefore, when writing a fluid–solid coupling program, the time step taken by the finite element method is used as the real time. The peridynamics only imports the pore water pressure and other data and calculates them until convergence after the FEM module calculates one or several time steps. The program flowchart
Figure 5 is as follows:
3.2. Transient Flow in Unsaturated Soil Pillar
To verify the accuracy of the finite element method for simulating unsaturated and unstable seepage problems, the soil column model [
48] shown in
Figure 6 was established. The soil column height
H = 1.0 m, width
, saturated permeability coefficient
. The parameters of the VG model are
,
,
a = 1.0 m
−1, and
n = 1.53. The initial pressure head of the soil column is
.
The soil column FEM model uses a four-node quadrilateral element, which is evenly divided into 500 elements, with a time step of , and a calculation time of 468,000 time steps. The AB and CD sides of the soil column are constant head boundaries, with head values of 0 m and −8 m respectively, and the AD and BC sides are impermeable boundary conditions.
Figure 7 shows the comparison results of the FEM solution and Warrick’s analytical solution [
49] at different instants. The results show that as the seepage progresses, the upper part of the soil column gradually becomes saturated, and the pressure head below the saturated zone decreases sharply, while the pressure head below the soil column still maintains the initial state and is not affected by the seepage. The simulation results of the finite element can correspond well with the analytical solution, verifying the accuracy of the model.
3.3. Numerical Validation of Fluid–Structure Interaction Scheme
To verify the feasibility and accuracy of the above fluid–solid coupling model, a soil column model with a height of 1m and a width of was established. Its mass density is , the permeability coefficient is , and the compression modulus is 3.4 Mpa. The sides and bottom of the soil column are impermeable, and displacement constraints in the x direction are applied to the sides of the column, and displacement constraints in the x and y directions are applied to the bottom of the column. Water seepage and soil compression only occur in the vertical direction. The top of the soil column is drained on one side and is subjected to a uniformly distributed downward load of 100 kPa.
Reference [
50] gives an analytical solution for the pore water pressure at any depth z of a soil column at different times t.
where
is the depth of soil column,
m is a positive odd number,
is the consolidation coefficient of soil,
is the density of water,
is the dimensionless time factor,
is the compression modulus,
is the uniformly distributed load applied on the top.
Figure 8 shows a comparison between the analytical solution and the numerical solution of the pore water pressure distribution at
t = 100 s. The pore water pressure in the lower section of the soil column is linearly related to the depth, while the pore water pressure in the upper section decreases slowly. It can be seen from the images that the numerical solution is in good agreement with the analytical solution, verifying the accuracy of the coupling model.
Figure 9 is a comparison of the Abaqus finite element solution and the homemade finite element solution of the vertical displacement distribution at
t = 100 s. Since the pore water pressure bears part of the external load, the vertical displacement of the soil column regarded as an elastic body is not linearly distributed along the height, and the maximum settlement value of the soil column is reduced. It can be seen from the image that the development trends of the two curves are the same, which further verifies the accuracy of the coupling model.
3.4. Effect of Rainfall Duration on Slope Stability
Select the slope model in the example in
Section 2 for further simulations. The two side boundaries AB and CD of the slope constrain displacement in the x direction; the bottom BC of the slope constrains displacement in the x and y directions; AB and BC are impermeable boundaries. The constant head boundary conditions are applied on the CD and DE sides, where the pressure head on the CD side is linearly distributed
and the pressure head on the DE side is always 0. Flow boundary conditions are applied on the AF and FE edges to simulate rainfall infiltration. Heavy rainfall often forms runoff and water accumulation on the slope surface. Since the depth of water accumulation is difficult to determine, this paper assumes that there is no surface water accumulation during the rainfall process, and the rainwater that does not infiltrate will be drained away immediately. The initial saturation distribution of the slope is shown in
Figure 10. The groundwater level is 3.25 m high and parallel to the bottom of the slope.
If the flow boundary condition is simply applied, pressure infiltration will occur on the slope surface, which is unreasonable if there is no surface water accumulation. Therefore, this paper reduces the flow rate by adding a proportional coefficient to avoid pressure infiltration, and the reduced flow rate is regarded as discharge from the slope. The initial value of the proportional coefficient is 1. When the pressure head of the flow node is greater than zero, let , be the single reduction size, and it can be taken as 0.5 for the first reduction. If the pressure head is still greater than zero after the reduction, it will continue to be reduced with the current reduction size until the pressure head is less than zero. At this time, it is determined that the true proportional coefficient exists in the interval composed of the coefficients before and after the reduction. The average value of the coefficients before and after the reduction is used as the search result for this time and substituted into the trial pressure head. If the accuracy requirement is met, it will be used as the flow boundary condition to continue the seepage calculation; otherwise, the reduction coefficient will be halved and the binary search will continue.
To observe the effect of rainfall duration on slope stability, the pressure head distribution cloud diagrams before rainfall, 6 h, 12 h, and 24 h of rainfall are given when the rainfall intensity is 5 mm/h, namely
Figure 11. When there is no rainfall, the pressure head is distributed linearly with height, decreasing from 3 m to −10 m from bottom to top. As the rainfall progresses, the top and surface of the slope quickly reach saturation 6 h after rainfall, the pressure head becomes zero, a transient saturated zone is formed, and a wetting front is generated at the junction with the unsaturated zone. The maximum negative pressure head appears below the wetting front. Since the rainfall boundary at the foot of the slope is close to the groundwater level, rainwater can quickly flow into the groundwater, causing the water level below it to rise slightly; while at the top of the slope, the seepage distance is long and there is a moist front, and there is still a large unsaturated area below, so the corresponding groundwater level does not change significantly.
When the rainfall reaches 12 h, the wetting front descends, the slope matrix suction decreases, the groundwater level rises further, and the maximum pressure head appears on the right side of the slope bottom. After 24 h of rainfall, the saturated zone above the moist front further developed, wrapping the unsaturated zone in the middle of the slope. The groundwater level increased significantly and gradually developed into the interior of the slope with the continuous influence of seepage.
The equivalent plastic strain distribution before rainfall, 6 h of rainfall, 12 h of rainfall, and 24 h of rainfall is shown in
Figure 12. The equivalent plastic strain is the largest at the slope toe. As rainfall continues, a transient saturated zone is formed at the top and below the slope, weakening the total stress of the soil. At the same time, the dragging effect of the volume force generated by rainwater seepage on soil particles will also harm the slope stability. Since the density of rainwater is much greater than that of air, the weight of soil in the transient saturated zone increases significantly and the sliding force increases. Under the influence of the above factors, the high plastic strain zone at the toe of the slope continues to develop forward, the plastic strain zone gradually penetrates and expands, and eventually leads to instability of the slope along the slip zone.
3.5. Effect of Rainfall Intensity on Slope Stability
To observe the effect of rainfall intensity on slope stability, the slope pressure head and equivalent plastic strain after 24 h of rainfall with intensities of 1 mm/h, 2.5 mm/h, 5 mm/h, and 7.5 mm/h were compared. The pressure head distribution is shown in
Figure 13. As the rainfall intensity increases, the transient saturated zone will develop deeper into the slope, the pressure head at the wetting front will increase significantly, and the matrix suction will dissipate faster. When the rain intensity is 1 mm/h, there is no significant change in the pressure head at the foot of the slope, and when the rainfall intensity increases to 5 mm/h, the pressure head at the slope foot increases significantly, weakening the total stress of the soil, making the slope foot more susceptible to damage. There is no obvious change in the pressure head for rainfall intensities of 7.5 mm/h and 5 mm/h. This is because the rainfall surface permeability has reached its limit. Although the rainfall intensity of 7.5 mm/h is higher, the flow rate higher than the surface permeability does not participate in the seepage process.
The equivalent plastic strain change is shown in
Figure 14. As rainfall intensity increases, the equivalent plastic strain zone gradually penetrates to the top of the slope. However, no matter how the rainfall intensity changes, the plastic strain at the bottom of the slope is always much larger than that at the top of the slope. On the one hand, this is because the bottom of the slope is subjected to a large gravity load from the soil above, and there is a geometric mutation at the toe of the slope that leads to stress concentration. On the other hand, since the ED side is a constant head boundary condition and the EF side connected to it is a flow boundary condition, part of the rainwater will flow out from the ED side through the toe of the slope as soon as it infiltrates into the slope, resulting in a large head gradient here. The resulting penetration force drags the soil particles, which also increases the plastic strain. At the same time, the groundwater level makes the pressure head at the bottom of the slope much greater than that at the top of the slope, further weakening the effective stress of the soil and harming the stability of the slope.
The equivalent plastic strain at the slope foot changes with time under different rainfall intensities, as shown in
Figure 15. The growth rate of equivalent plastic strain was high before 5 h, and then it gradually slowed down and tended to be stable. The reason is that 5 h ago, a part of the soil near the slope toe was still unsaturated. As the rainfall progressed, the weight of this part of the soil increased, and the pore water pressure increased, which jointly affected the development of equivalent plastic strain. When all the nearby soils reached saturation, the weight tended to be stable. At this time, only the pore water pressure increased with time, resulting in a slowdown in the growth rate.