Abstract
An ABAQUS UMAT subroutine was used for the secondary development of the established coupled hydromechanical constitutive model of unsaturated soil considering the effect of the microscopic pore structure. Combined with Euler’s backward implicit integration algorithm, a numerical program was established for simulating the proposed model. The developed numerical program was used to simulate the rainfall infiltration process of an actual slope engineering example, and the effects of rainfall intensity and rainfall duration on the pore pressure, fluid velocity, and displacement of the unsaturated soil slope were analyzed. The results show that the developed numerical program can reasonably analyze the changes in the seepage field and displacement field of unsaturated soil slopes under rainfall infiltration.
1. Introduction
Rainfall infiltration causes complex changes in the hydromechanical characteristics of unsaturated soil slopes. Rainfall is the main cause of landslide instability [1–3], which has an important influence on the safety and stability of unsaturated soil slopes. During rainfall infiltration, the water content and bulk density of the slope soil will increase, and the matric suction will decrease, causing a decrease in the shear strength of the soil. In addition, rainfall infiltration is a continuous dynamic process, and slope instability usually occurs during continuous rainfall or a short time after the rain stops [4]. Therefore, it is of great significance to comprehensively and accurately analyze the hydromechanical characteristics of unsaturated soil slopes under rainfall infiltration.
At present, many scholars use the seepage-deformation coupling theory to study the slope rainfall infiltration problem [5–14]. Davies et al. [5] established a coupled hydromechanical model for slope stability analysis, revealing the importance of pore pressure. Wang et al. [6] studied the coupling characteristics of seepage and deformation of slopes under rainfall, revealing the interaction characteristics between slope deformation and matric suction. Wang and Li [7] found that seepage in unsaturated soil is very sensitive to the seepage-deformation coupling effect and related boundary conditions. Lu and Godt [8] proposed a definition of the effective stress of unsaturated soil considering the effect of stable seepage and applied the effective stress theory to the analysis of infinite slopes under rainfall. Cai and Ugai [9] combined the finite element strength reduction method with saturated-unsaturated seepage analysis to analyze the influence of different soil properties on the slope stability under rainfall. The rainfall infiltration of an unsaturated soil slope is a complicated hydromechanical coupling problem; in such cases, it is necessary to use advanced coupled hydromechanical constitutive models [15–18] for more accurate and reasonable analysis. Most of the existing research results are based on the effective stress principle of unsaturated soil, and ideal elastic-plastic theory is often used to calculate the deformation of unsaturated soil without considering the influence of saturation variation on its material properties and mechanical properties. Seepage theory mostly adopts the calculation method of uncoupled seepage and deformation, which fails to consider the hydromechanical coupling effect of unsaturated soil, so it cannot truly reflect the change in the seepage field under rainfall conditions. In addition, the microscopic pore structure of unsaturated soil has a significant effect on its hydromechanical coupling characteristics [19, 20], and the influence of the microscopic pore structure during rainfall infiltration cannot be ignored. However, most of the existing coupled hydromechanical models for unsaturated soil do not consider the effect of microscopic pore structure.
In this paper, the coupled hydromechanical constitutive model of unsaturated soil considering the influence of the microscopic pore structure (the modified Glasgow coupled model, MGCM) established by the authors in [21] is summarized first. An ABAQUS UMAT subroutine is used for the secondary development of the established model. By incorporating Euler’s backward implicit integration algorithm, a numerical program that can analyze the hydromechanical coupling characteristics of unsaturated soil slopes under rainfall infiltration is established. The established numerical program is used for the numerical simulation of an actual slope project under rainfall infiltration. The seepage field and displacement field of the unsaturated soil slope under rainfall are given. The influences of rainfall intensity and rainfall duration on the pore pressure, fluid velocity, and displacement of the slope are analyzed. The results show that the numerical program developed in this paper can reasonably describe the hydromechanical coupling characteristics of unsaturated soil slopes under rainfall infiltration.
2. MGCM Overview
2.1. Constitutive Variables
The effective stress [22] of the micropores (Bishop parameter replaced with the effective degree of saturation [23]) and the modified suction are used as the two constitutive variables:where is the mean net stress, is the total mean stress, is the major principal stress, is the intermediate principal stress, is the minor principal stress, is the effective degree of saturation, and are the air pressure and the pore pressure of unsaturated soil, respectively, is the matric suction, and is the porosity.
The incremental expressions of the constitutive variables can be obtained as follows according to equation (1):where volumetric strain increment and effective degree of saturation increment are the conjugate strain increments corresponding to the effective stress and modified suction considering the micropores, respectively.
2.2. Yield Equations
Yield curve equations for the LC, SI, and SD curves are given as follows:where , , and are hardening parameters, defining the current locations of the three curves.
The yield equation in the stress plane is defined by the modified Cambridge Clay model [24]:where is the deviator stress and is the slope of the critical state line.
Without considering the influence of the Lode Angle, the LC yield surface is extruded to the principal stress space and changed towhere is the second invariant of the deviator stress tensor. Thus, in the modified stress space, LC, SI, and SD yield surfaces constitute a closed stress space, which form the yield surfaces of the MGCM model.
2.3. Hardening Law
The hardening equations of the solid phase and liquid phase are expressed as follows:where ν is ratio volume, and are coupling parameters, is the slope of the saturated normal compression line, and is the slope of the boundary line of the SWCC. is the slope of the saturated unloading line in the plane. is the slope of the scanning line of the SWCC. is the plastic volumetric strain increment, and is the plastic effective degree of saturation increment.
2.4. Flow Rule
During yielding, one or more plastic mechanisms are activated. In the MGCM model, associated flow rule is assumed in LC, SD, and SI yield curves. Yielding on only the LC curve produces plastic volumetric strains but is assumed to not cause plastic changes in the effective degree of saturation.
Similarly, yielding on the SI or SD curves produces plastic changes in the effective degree of saturation but is assumed to not induce plastic volumetric strain. Therefore, the associated flow rules applied to the SD or SI curve are expressed as follows:
2.5. Stress-Strain Relation
The elastic strain increment equations in isotropic compressive stress state are as follows:and the plastic strain increment relationship of unsaturated soil in isotropic compressive stress state is expressed as follows:
In real stress space, the stress and strain tensors are given by the following vectors:and the matrix form of the stress-strain relationship is as follows:where is the stress increment vector, is the 7×7 generalized elastic matrix, is the elastic strain increment vector, is the plastic strain increment vector, is the total strain increment vector, and is the 7×7 elastic-plastic matrix. The form of depends on the specific plastic mechanisms that are active during loading. In the MGCM model, there are four plastic active mechanisms. The LC plastic mechanism alone, either the SD or SI mechanism, or both the LC and either the SD or SI plastic mechanisms are active simultaneously.
To derive the plastic multipliers of the MGCM model, the consistency conditions of the yield functions of LC and SI/D are as follows:
According to flow rules (7) and (8) and the consistency conditions (13) and (14), for all the possible elastoplastic cases can be acquired in Appendix A.
The developed numerical program uses Euler’s backward implicit integration algorithm for numerical integration.
2.6. Permeability Coefficient Function
In this paper, the influence of rainfall on slope instability is mainly discussed, so the water retention characteristics and permeability of soil need to be quantified. The water retention capacity of soil has been included in the coupled hydromechanical constitutive model of unsaturated soil considering the microscopic pore structure; that is, the influence of deformation on the SWCC is considered in the model. Therefore, it is not necessary to define SWCC in the software when using ABAQUS for numerical simulation. The permeability coefficient of unsaturated soil is related to its matric suction (volumetric moisture content) and void ratio. In this paper, the permeability coefficient function refers to the calculation formula in [25], which only defines the relationship between the permeability coefficient and matric suction, without considering the influence of deformation on seepage. The permeability coefficient function is given as follows:where is the permeability coefficient of unsaturated soil, is the permeability coefficient of saturated soil, and , , and are the model parameters. The relationship curve corresponding to equation (15) is shown in Figure 1.
3. Numerical Analysis of an Unsaturated Soil Slope under Rainfall Infiltration
3.1. Project Overview
A slope engineering project in the suburbs of a South China city is taken as the research object. There are five layers of soil within this slope. From bottom to top, the first layer in the slope is Quaternary fully renewed soil with poor engineering geological properties. It is characterized by moderate compressibility, moderate-to-low strength, and low permeability. The second layer is the Quaternary Pleistocene Shu loess with moderate-to-low compressibility, moderate-to-high strength, low permeability, and weak expansion properties. The third layer is the Quaternary-lower Pleistocene gravel layer with moderate-to-high strength and high permeability. The fourth layer is Quaternary-lower Pleistocene clay with low compressibility, high strength, and low permeability. The fifth layer is upper Tertiary medium-grained sandy gravel with high permeability. The groundwater in the slope is mainly phreatic water supplied by atmospheric precipitation infiltration and concentrates in the gravel layer. The bedrock is composed of silty mudstone and basalt with tight fractures, low permeability, and an abundance of water, including some bedrock fissure water. The topography of the slope fluctuates significantly, the water content in the slope is nonuniformly distributed, and the location of the groundwater level is closely related to the terrain elevation and rainfall. The stable water level of the phreatic water is approximately 8.20–17.24 m, which is not affected by precipitation. A schematic diagram of the slope model is shown in Figure 2.
3.2. Slope Model of Unsaturated Soil
According to slope engineering conditions, on the basis of properly optimizing the slope structure and simplifying the soil layer distribution, a schematic diagram of the ABAQUS slope calculation model is shown in Figure 3. Points A, B, C, and D are the selected points in the slope. The horizontal and vertical displacements of the bottom boundary and the horizontal displacements of the left and right boundaries of the slope are fixed, and the other boundaries are free. It is assumed that the groundwater level is located at the top of the sandy gravel layer; that is, the water level on the left side is at height of 10 m and that on the right side is at height of 15 m. The soil below the groundwater level is saturated, and the soil above the groundwater level is unsaturated. The initial pore pressure distribution inside the slope is calculated as a linear distribution of , where is the distance from the groundwater level. The top boundary of the model is set as the rainfall boundary, and the influence of the slope angle on the infiltration intensity is considered. The related parameter values of each layer of soil are shown in Table 1.
(a)
(b)
The simulation process is divided into two steps. The first step is the initial stress balance step, in which the gravity load, initial stress field, initial void ratio, and initial pore pressure are set. The second step is the rainfall analysis step, in which the infiltration intensity is set on the slope surface according to the designed rainfall scheme.
3.3. Rainfall Scheme
The city is located in a humid and semihumid monsoon climate zone, with an annual average temperature of 13.6°C. The minimum annual rainfall is 348.0 ml (1978), the maximum annual rainfall is 2015.2 ml (2009), and the average annual rainfall is 682.5 ml. According to the distribution of rainfall in 2017 (see Figure 4), the precipitation from June to October accounts for approximately 60% of the precipitation that occurs throughout the year. The precipitation time is relatively concentrated every year, and heavy rain often occurs. Therefore, the influences of rainfall intensity and rainfall duration on the seepage field and displacement field of the slope are mainly considered in this paper.
Referring to the meteorological conditions of this area and the rainfall classification standard issued by the meteorological department, it is determined that the rainfall intensities considered for analysis are 40 mm/d (heavy rain), 80 mm/d (torrential rain), and 160 mm/d (downpour), and the rainfall durations considered are 12 h and 24 h.
3.4. Seepage Field Simulation of Unsaturated Soil Slope under Rainfall Infiltration
3.4.1. Influence of Rainfall Intensity on the Slope Seepage Field
Figure 5 shows the distributions of pore pressure inside the slope at the end of 12 h continuous rainfall events under different rainfall intensities. With the increase in rainfall intensity, the rainwater infiltration depth and the minimum pore pressure of the slope increase. The initial minimum pore pressure is −304.6 kPa; the minimum pore pressure is −205.1 kPa when the rainfall intensity is 40 mm/d, the minimum pore pressure is −141.2 kPa when the rainfall intensity is 80 mm/d, and the minimum pore pressure is −98.9 kPa when the rainfall intensity is 160 mm/d; that is, with the increase in rainfall intensity, the rate of increase in the minimum pore pressure continuously decreases. This trend forms because when the rainfall intensity is 40 mm/d, the rainfall intensity is much smaller than the saturation permeability coefficient of the silty clay on the top of the slope, so the rainwater can quickly infiltrate into the slope, and its impact range is only the top of the slope. When the rainfall intensity is 80 mm/d, the rainfall intensity is slightly greater than the permeability coefficient of the silty clay at the top of the slope, and its influence range partially extends to the clay layer above the slope. Due to the low permeability of the clay layer, the rate of increase in the minimum pore pressure will decrease. When the rainfall intensity is 160 mm/d, the influence range of rainwater infiltration mostly covers the clay layer above the slope, so the rate of increase in the minimum pore pressure value is lower. Because the groundwater level inside the slope is low and the soil layers above the groundwater level are clay and silty clay, respectively, the permeability is poor, so the increase in the groundwater level caused by the rainfall infiltration is not obvious.
(a)
(b)
(c)
Figure 6 shows the distributions of pore pressure inside the slope at the end of the 12 h continuous rainfall events under different rainfall intensities. As the rainfall intensity increases, the maximum velocity inside the slope continues to increase, which is caused by part of the rainwater infiltrating into the slope, resulting in a wider fluid flow channel. In addition, with the increase in rainfall intensity, the range of flow velocity change continuously extends from the left boundary to the right one. This is because the initial pore pressure and saturation of the soil on the left side of the slope are high, resulting in a large gradient. As the influence range of rainwater infiltration increases, the difference in pore pressure between the left and right sides of the slope continues to decrease, so the infiltration rate gradually becomes uniformly distributed. Because the gravelly sand layer present in the middle of the slope has a high permeability and is completely in the saturated area, the maximum value of the velocity is always located in the gravel layer. The presence of the gravelly sand layer causes the infiltration of rainwater to be quickly discharged, and it is necessary to pay special attention to the surface water accumulation phenomenon caused by groundwater overflow in the surrounding low-lying areas.
(a)
(b)
(c)
Figures 7 and 8 show the simulation results of pore pressure and fluid velocity with rainfall duration at point A under different rainfall intensities. The variations in pore pressure and fluid velocity at points B, C, and D with rainfall duration are similar to those at point A. Clearly, with the development of rainfall duration, the pore pressure and fluid velocity both increase; additionally, the pore pressure and fluid velocity both increase with the rainfall intensity.
3.4.2. Influence of Rainfall Duration on the Slope Seepage Field
Figure 9 shows the distribution of pore pressure inside the slope at the end of the 24 h continuous rainfall event under an 80 mm/d rainfall intensity. A comparison of Figures 9 and 5(b) shows that, with the same rainfall intensity, the infiltration depth of the rainwater will continue to increase with the extension of the rainfall duration, but, due to the presence of the middle clay layer, the infiltration depth will first increase slowly and then increase rapidly, and the rate of increase in the pore pressure exhibits the same trend. In addition, a comparison of Figures 9 and 5(c) shows that the minimum pore pressure is −134.2 kPa when the rainfall intensity is 80 mm/d for 24 h, and the minimum pore pressure is −98.9 kPa when the rainfall intensity is 160 mm/d for 12 h. Under these two rainfall conditions, the slope experienced the same amount of rainfall, but the trends of the increase in pore pressure are different, indicating that the corresponding rainfall infiltration depths are different. Therefore, for the rainfall infiltration depth of this slope, under the same amount of rainfall, the impact of rainfall intensity is greater than the impact of rainfall duration.
Figure 10 shows the distribution of the fluid velocity inside the slope at the end of the continuous rainfall event under an 80 mm/d rainfall intensity for 24 h. A comparison of Figures 10 and 6(b) shows that, with the increase in rainfall duration, the maximum fluid velocity inside the slope will continue to increase, but its rate of increase will gradually decrease. Similar to the increase in rainfall intensity, with the increase in rainfall duration, the range of the change in fluid velocity continuously extends from the left boundary to the right boundary.
Figures 11 and 12 show the simulation results of pore pressure and fluid velocity with rainfall duration at points A, B, C, and D under an 80 mm/d rainfall intensity for 24 h. At the beginning of the rainfall event, the pore pressure at point A increases the fastest, followed by those at points B and C, while that at point D is the lowest. With the extension of the rainfall duration, the rate of increase in the pore pressure at point D hardly changes, while the rates of increase in the pore pressure at points A, B, and C all decrease, especially at point A. Figure 11 also indicates that the rates of increase in the pore pressure at points A, B, C, and D will eventually become equal. After 24 h of rainfall, the pore pressure values at points A, B, C, and D are −133.62 kPa, −114.67 kPa, −114.99 kPa, and −94.86 kPa, while the fluid velocity values at points A, B, C, and D are 1.9E-3, 3.1E-3, 3.1E-3, and 3.3E-3, respectively, so it can also be proven that the velocity and pore pressure at any point in the slope are proportional during rainfall. The pore pressure and fluid velocity of points A, B, C, and D generally increase as the height of the selected point decreases.
3.5. Displacement Field Simulation of Unsaturated Soil Slope under Rainfall Infiltration
3.5.1. Influence of Rainfall Intensity on the Slope Displacement Field
Figure 13 shows the distributions of vertical strain inside the slope at the end of the 12 h continuous rainfall events under different rainfall intensities. As the rainfall intensity increases, the vertical strain inside the slope changes due to rain infiltration. This is because, within the influence range of the rain infiltration, the matric suction inside the slopes decreases, resulting in the constant change of yield surface of the MGCM, which causes part of the soil to yield and undergo changes in strain.
(a)
(b)
(c)
Figure 14 shows the distributions of horizontal displacement of the slope at the end of the 12 h continuous rainfall events under different rainfall intensities. The developed numerical program can accurately reflect the hydromechanical coupling effect inside the unsaturated slope during rainfall, and the trend of the change in the horizontal displacement is consistent with that obtained by the limit equilibrium analysis method. In addition, during the rainfall event, the area of change in the slope horizontal displacement generally forms an irregular arc passing through the slope foot. With the increase in rainfall intensity, this area continues to expand, and the maximum horizontal displacement continues to increase.
(a)
(b)
(c)
3.5.2. Influence of Rainfall Duration on the Slope Displacement Field
Figure 15 shows the distribution of the vertical strain inside the slope at the end of the continuous rainfall event under an 80 mm/d rainfall intensity for 24 h. A comparison of Figures 15 and 13(b) indicates that, under the same rainfall intensity, with the increase in rainfall duration, the vertical strain inside the slope continues to increase, but, due to the presence of the middle clay layer, the earlier vertical strain develops quickly, and the later vertical strain develops slowly.
Figure 16 shows the distribution of the horizontal displacement inside the slope at the end of the continuous rainfall event under an 80 mm/d rainfall intensity for 24 h. A comparison of Figure 16 with Figure 14(b) indicates that as the rainfall duration increases, the horizontal displacement generated by the slope continues to increase, but the rate of increase gradually decreases. In addition, a comparison of Figures 16 and 14(c) shows that the maximum horizontal displacement is 6.955E-3 m when the rainfall intensity is 80 mm/d for 24 h and that the maximum horizontal displacement is 8.65E-3 m when the rainfall intensity is 160 mm/d for 12 h. Under these two rainfall conditions, the slope experienced the same amount of rainfall, but the horizontal displacement results vary, indicating that the corresponding areas where the water-force coupling effect occurs are not the same. Therefore, regarding the stability of this slope, under the same amount of rainfall, the impact of rainfall intensity is greater than the impact of rainfall duration.
4. Conclusion
In this paper, the coupled unsaturated soil-water-force model considering the influence of the microscopic pore structure previously established by the author is imported into an ABAQUS and a UMAT subroutine for secondary development. The developed numerical program is used for the numerical simulation of the seepage field and displacement field for an actual slope engineering example during rainfall infiltration, and the following conclusions are obtained:(1)With the increase in rainfall intensity or increase in rainfall duration, the pore pressure, fluid velocity, vertical strain, and horizontal displacement of the slope all increase.(2)The existence of the clay layer in the middle of the slope has a certain influence on the changes in the seepage field and displacement field in the slope.(3)When the amount of rainfall is the same, the impact of rainfall intensity on the internal seepage field and displacement field of the slope is greater than the impact of rainfall duration. Therefore, when analyzing the stability of a slope, special attention should be paid to short-term heavy rainfall.
Appendix
Matrix ep
When yielding is taking place on the LC surface alone, the matrix of ep iswhere
The expressions of and are
When yielding on only SD or SI (SI/D) is taking place,where
and the expressions of and are
When LC and SI/D plastic mechanisms act simultaneously,where
The expressions of , , , , and areand nT = (1, 1, 1, 0, 0, 0).
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
All authors approved the manuscript.
Acknowledgments
This research was financially supported by the National Natural Science Foundation of China (51722802, 51678041, and U1834206), Beijing Natural Science Foundation (8202038), and Beijing Nova Program (Z181100006218005).