Abstract
Seepage flow is one of the primary factors that trigger slopes and landslides’ failure. In this study, the slope stability under saturated or unsaturated conditions is analyzed. The influence of a complex stress state on the slope stability with the saturated or unsaturated seepage flow is proposed in this paper. Firstly, an elastoplastic constituted model for the soil under a complex stress state is established and as a user subroutine of the finite element method code of FLAC. Secondly, the 2D and 3D problems of slope stability influenced by the saturated or unsaturated seepage flow are analyzed via the finite difference method with the influence of the complex stress state. Finally, the influence of the intermediate principal stress and the saturated or unsaturated seepage flow on the slope stability is analyzed in this study.
1. Introduction
Soil is a typical polyporous and multiphase material on our planet. Therefore, soils usually contain soil particles, pore water, and pore air. The soil can be divided into dry, saturated, and unsaturated states via the pore filling composition. The seepage flow of pore fluid will influence the stability of the soil significantly. Soil slope stability is one of the key and classical problems in soil mechanics. The variation of the water content in soil (due to underground water, seepage flow, or rainfall) is a vital influence factor for soil slope safety and stability. In recent years, soil slope stability, i.e., hill slope, embankment, and cutting high slope, is influenced by the water content variation due to climatic and anthropogenic factors [1–3]. Many true triaxial tests have verified that the intermediate principal stress has specific influences on mechanical behavior and strength characteristics of soil, and this problem was summarized by Ma et al. [4–6].
The factor of safety (FOS) is an essential parameter to determine the slope or landslide stability, and several methods have been suggested for FOS calculation, including analytical and numerical methods. The limit equilibrium method (LEM) is a general analytical method for slope stability analysis and FOS calculation, e.g., the circular slip surface method or the general methods of slices [7–11]. The limit analysis method is a semianalytic method for the problem of slope stability analysis. The stability of three-dimensional undrained slopes was analyzed by Li et al. using the numerical limit analysis method [12]. The FOS of a slope also can be calculated via the numerical method with the strength reduction method (SRM), e.g., finite element method (FEM) or finite difference method (FDM). SRM is a method that the original shear strength parameters (e.g., cohesion and friction angle ) of soil mass are divided by the strength reduction factor (SRF) to bring the slope to the point of failure. The shear strength parameters of the soil mass are decreased gradually by SRF until the slope becomes unstable, and the value of SRF when the failure is initiated is equal to FOS for the slope [13–16]. The numerical method has several advantages over the analytical method (e.g., LEM) for slope stability analysis, e.g., finding the critical failure surface automatically. The relationship between the LEM and SRM in slope stability analysis was discussed by Griffiths and Lane [17]. A model test for soil slope failure due to the seepage flow is proposed by Jia et al. [3]. The soil slope failure caused by rainfall and infiltration is analyzed by Chen et al. [1].
In this paper, an explicit finite difference elastoplastic model is developed based on the twin shear strength theory (TSST). TSST is proposed by Yu [18, 19] that takes the effect of intermediate principal stress on the failure of materials into account. The model is written in C++ and compiled as a dynamic link library file that can be loaded into an explicit finite difference code FLAC. According to the model test and field investigation, the stability analyses of an embankment slope and a high loess slope are performed with FLAC, where the complex stress state, groundwater level variation, and seepage flow are taken into account. The problems of fractured wells in a reservoir of hydrogen and carbon based on a dual-porosity medium model were discussed by Xue et al. and Liu et al., respectively [20, 21].
2. Theory and Method
2.1. Twin Shear Strength Theory
The principal stress form of TSST can be written as [18, 19] where is the ratio of tensile and compression yield strength, is the tensile yield strength, and is the coefficient which reflects the influence of intermediate principal stress. The three principal stresses are denoted as , and the direction of principal stress is positive in tension and negative in compression. TSST can be written as the form of the strength parameters and , which are widely used in geotechnical engineering; it follows that
2.2. Explicit Finite Difference Form
From the incremental elastoplastic theory, the total strain can be written as
The elastic portion follows increment generalized Hooke’s law: where , , is bulk modulus, and is shear modulus. The tensor form of Eq. (4) is given:
From the plastic flow rule, the plastic incremental deformation part can be determined: where is plastic potential and is the plastic multiplier.
In order to deduce the expression of , the critical state between elastic and plastic is assumed: an increment stress loaded to a structure when , and the state of the structure has no change (). From Eq. (5), the component of increment stress can be written as
Substituting Eq. (6) into Eq. (7) and simplifying where is the corresponding stress of plastic strain. can also be written as where , and then, follows:
Substituting Eq. (8) into Eq. (10) and then can be deduced:
In order to process the material nonlinearity, the stress and is defined as follows [22]: where is the total strain of the elastic stage, is a stress tensor in the elastic stage, and is a stress tensor in the plastic stage. can be written as
Substituting Eq. (14) into Eq. (11) and then can be deduced:
From Eqs. (6), (12), and (13), the expression can be deduced:
Eq. (2) can also be written as where . The nonassociated flow is employed, the dilation angle substitutes the internal friction angle in Eq. (17) and as the equation of the plastic potential, where , and the item of can be determined:
Substituting Eq. (18) into Eq. (16), the component of yields where the expression of is as follows:
Eqs. (18) and (19) are the formats of Lagrangian explicit finite difference for the twin shear elastoplastic model. The model is written in the language of C++ and compiled as a file of dynamic link library that can be loaded in FLAC code [22]. Ma et al. propose a numerical study of gravel soil ground’s dynamic compaction using the explicit discrete element method [23].
2.3. Slope Stability due to Saturated Seepage Flow
Due to groundwater level variations, the FOS of the slope is calculated according to the saturated seepage flow theory without the influence of unsaturation on slope stability during seepage flow. FOS values for a 2D slope (with the plane strain condition) are calculated by SRM using TSST. The FOS of the soil slope is calculated via SRM, i.e., the same factor SRF is applied to both of the shear strength parameters ( and ). The reduced strength parameters ( and ) are given by where and are the original cohesion and internal friction angle of the slope soil, and are the cohesion and internal friction angle of soil after the reduction, and SRF (strength reduction factor) is the intensity reduction coefficient. The strength parameter of slope soil is divided by a reduction factor SRF step by step until the slope failure, and then, the value of SRF equals FOS of the slope. The values of SRF at the critical failure state of the slope are taken as FOS of the slope [17].
Analysis and calculation parameters for slope stability with seepage flow are shown in Table 1. Figure 1 shows the relationship between FOS and groundwater level calculated by saturated seepage. It can be seen from Figure 1 that under saturated seepage, the slope safety coefficient calculated according to Mohr-Coulomb theory is lower than 1.0 when the groundwater level is raised. The values of FOS are greater than 1.0 when the effect of intermediate principal stress is considered. The model test results in the literature [3] showed that the slope did not fail in the water injection stage (the groundwater level is raised), so the calculation results without taking the effect of intermediate principal stress into account do not agree with the actual situation. Figure 2 shows the maximum shear strain contours of the slope calculated according to TSST () when the slope water level drops to under saturated seepage condition. Figure 2 shows that the nonassociated flow () calculation results show that there are multiple sliding surfaces of the slope. However, the results obtained by associated flow () show only one sliding surface.
(a) Nonassociated flow
(b) Associated flow
2.4. Slope Stability due to Unsaturated Seepage Flow
The effect stress can be formulated as follows using Bishop’s unsaturated effect stress theory (compression is negative) [24]: where is the total stress, is the saturation of water, is the air saturation, and Pw and Pa are the water pressure and air pressure. TSST can establish the formula of unsaturated shearing strength, and TSST can be written by the twin shear stress ( and ) and normal stress ( and ) as follows [18, 19]: where and is the strength parameters and can be determined by cohesion and friction angle (, ); TSST also can be formulated by and as follows: where the normal stress and under the unsaturated condition is as follows:
The shear strength formula of unsaturated soil can be established via substituting Eq. (25) into Eq. (24), and the shear strength formula of unsaturated soil can be obtained as follows:
The principal shear stress (, , and ) and normal stress (, , and ) can be expressed as follows:
The shear strength formula of unsaturated soil also can be formulated via substituting Eq. (27) into Eq. (26):
The variation of suction of the unsaturated soil can be described by the empirical formula recommended by van Genuchten [25–27]: where is the volume water content of the soil, is the residual volume water content of the soil, and is the measured parameters. The relationship between saturation and volume water content is , where is soil’s porosity. The shear strength formula of unsaturated soil based on TSST and the empirical formula of SWCC van Genuchten were used to analyze the stability analysis of unsaturated seepage slope in the slope model test in literature [3]. The lsqcurvefit function in Matlab’s software is used to fit the data of the soil-water characteristic curve (SWCC) tested by literature [3] and the empirical formula proposed by van Genuchten. Figure 3 is the fitting result of SWCC data of the slope model test soil and the empirical formula of van Genuchten. It can be seen that the hygroscopic and dehumidification curves of slope soil have an obvious hysteretic relationship. The matric suction of the soil’s unsaturated state in the calculation process is controlled by the empirical formula of van Genuchten, in which the stage of water level rise and fall is calculated by the fitting results of the empirical formula of van Genuchten in SWCC in Figure 3. According to the strength formula of unsaturated soil with TSST derived in the previous part, the slope’s safety factor is calculated by the strength reduction method. The numerical calculation model and boundary conditions are consistent with the saturated seepage calculation. Figure 4 shows the slope to water injection after reaching the top, saturated and unsaturated seepage flow theory to calculate the safety factor of slope, and the relationship between the theory of parameter ; it can be seen that the value of FOS yielded by Mohr-Coulomb theory is less than 1.0 both under saturated and unsaturated conditions and does not agree with the results of the slope model test in literature [3]. The value of FOS yielded by TSST () is greater than 1.0. Therefore, the intermediate principal stress influence should be considered in slope stability analysis due to the seepage flow. Figure 5 shows the maximum shear strain contours of slope under the unsaturated seepage when the water level dropped to yield by Mohr-Coulomb and TSST. The calculation result yield by TSST () in the associated flow () shows that the multiple sliding surfaces have occurred. However, the calculation result yield by Mohr-Coulomb theory or associated flow () shows that the slope only has one sliding surface.
(a) TSST
(b) TSST
(c) TSST
(d) TSST
2.5. Three-Dimensional Analysis of High Slope Stability due to Saturated Seepage Flow
The high loess slope is situated along the northern boundary of a loess plateau close to Jingyang, Shaanxi province, China. The distance from the high loess slope to the downtown of Xi’an is about 25 km. The Jing river flows beside the northern boundary of the loess plateau (see Figure 6). Due to the Jing river’s erosion, many high loess slopes came into being along the northern loess plateau boundary. The northern loess plateau boundary is about 30 km long, the elevation of 30-90 m above the Jing river, and the slope angle is about 45-80°.
The landslide mass’ length and width are about 400 and 410 m, respectively, and the slope angle is about 45-55°. Seven cross-sections are measured for the landslide mass in a longitudinal direction (see Figure 7). The continuous line represents the edge of the slope crest before the landslide event, and the dashed line represents the current edge of the slope crest, as shown in Figure 7. The cross-section sketch is drawn according to the topography data, as shown in Figure 8. The groundwater of the loess plateau is stored in a phreatic aquifer. In 1976, the groundwater level was kept equal to the Jing river (elevation is 380 m). However, large-scale irrigation and raining above the loess plateau induced a significant rise in the groundwater level. In 1992, the groundwater level’s depth under the loess plateau was 37 m (elevation is 425 m). Since 1976, the landslides in this area have occurred over 40 times, which enormously influenced the local area.
A series of laboratory tests, including the particle-size analysis test and consolidated-undrained triaxial tests, were conducted. The particle-size analysis test shows that the percent of particles with diameters greater than 0.075 mm is 27.8%. The percent of particles with diameters less than 0.005 mm is greater than 10%. The plasticity index ranges from 3.0 to 10.0. Therefore, the soil can be classified as clayey silt. The cohesion and internal friction angle of soil samples obtained from the consolidated-undrained triaxial test under different confining pressures are 3.56 kPa and 17.9°, respectively. Table 2 shows the physical properties of the soil samples from the high loess slope.
2.5.1. Influence of Intermediate Principal Stress and Groundwater Level
The previous terrain of the high loess slope before the landslide is created based on terrain from field investigation, and the three-dimensional model of high loess slope before the landslide event is established. The cross-sections of A-A and G-G are selected as the boundaries of two sides of the slope model, and the dimension and meshing of the model are shown in Figure 9. The dashed line represents the groundwater level, and is the depth of groundwater under the loess plateau. Node O of the slope top on the cross-section of D-D is selected. The displacement of node O in the direction of the -axis is recorded during the calculation process to analyze the slope stability. The calculation parameters are as follows: Young’s modulus MPa, Poisson’s ratio , and dilation angle . The calculations are performed under the groundwater level in 1976 ( m) and 1992 ( m) to compare the slope stability under different groundwater levels, respectively. The relationship between the displacement at node O in the -direction and the unbalanced force rate with the different magnitude of is shown in Figure 10. The unbalanced force rate can be defined as the maximum unbalanced force magnitude for all nodes divided by the average applied force magnitude for all the nodes. The equilibrium in the calculation is closely related to the unbalanced force rate. By default, the ratio limit for convergence is .
Figure 10 indicated that the slope changes towards instability with the decrease of the parameter , and the slope failure occurs when the magnitude of falls below a specific value. Because the high loess slope was stable before 1976, the stability analysis taking the intermediate principal stress’ influence into account can reflect the actual situation more accurately. The factor of safety for the high loess slope is calculated via the strength reduction method. The strength parameters of soil mass are decreased gradually by a factor until the slope becomes unstable. The factor when the failure is initiated is the factor of safety for the slope. The relationship between FOS values and the magnitude of under the groundwater level in 1976 and 1992 is shown in Figure 11. The same results can be observed from Figure 11, the values of FOS decreased with the decreasing value of the parameter , and the values of FOS also decrease with the rising level of the groundwater. The magnitudes of under the groundwater levels in 1976 and 1992 are equal to 0.3 and 0.5, respectively, when FOS values equal to 1.0. Figure 12 shows the pore-water pressure and displacement vector of slope in m and . As shown in Figure 12, the high loess slope becomes unstable when m and .
The stability analysis of the high loess slope during the rise in the groundwater level caused by raining and irrigation is performed with different depths (-50 m) of the groundwater level. Figure 13 shows the relationship between the displacement of node O in the -direction and unbalanced force rate under different depths of the groundwater level when . The results suggest that the displacement of node O in the -direction increases significantly as the groundwater rises gradually. The slope failure occurred until the depth of the groundwater level is close to 35 m.
2.5.2. Influence of Seepage Flow
According to statistics, landslides of the loess area induced by the distinct change of the groundwater level occurred frequently during the shot-time raining in the rainy season. The rainfall in the district of Xi’an fluctuates seasonally, and the most rainfall is yielded during the months of July to September, and the time of raining is very short. The mean monthly rainfall of July to September is about 12.5-16.4% of the mean annual rainfall. For example, the total rainfall of 2003 is 898.7 mm, and the total rainfall of July in 2003 is close to 140 mm, which is about 15% of the total rainfall in 2003. The high loess slope on the southern bank of the Jing river slided in July of 2003. Thus, it is necessary to analyze the influence of the distinct rise of the groundwater level and seepage flow which are caused by raining and irrigation on the stability of the high loess slope.
Parameters of soil are the same as the calculation of no seepage flow, and the soil permeability is 0.1 m/d. The calculations of seepage flow and stress field are parallel performed, and the stress field of the slope is calculated by the mechanical strains caused by pore-water pressure changes. Based on the calculation results of groundwater depth m, the calculation of seepage flow is performed with four scenarios: the rise of groundwater level m, 35 m, 40 m, and 45 m under the loess plateau.
Comparing with the numerical results of no seepage, the results of seepage flow calculation show that the displacement of node O in the -direction increases significantly when the unbalanced force rate is close to the target of convergence and the slope failure occurs suddenly when reaches a certain value. However, the results of no seepage calculation show that the slope reaches a static state when the unbalanced force rate decreased. On the other hand, the -displacement at node O from the seepage calculation is less than the result of no seepage calculation. This phenomenon suggests that the deformation of the loess slope is too limited to absorb the overload of seepage flow completely; thus, the sudden failure of the loess slope may occur more easily during the process of seepage flow.
3. Conclusion
The form of the explicit finite difference of the twin shear elastoplastic model is established and used to analyze the slope stability with saturated and unsaturated seepage flow. The following conclusions can be stated. (1)The intermediate principal stress effect of soil strength greatly influences FOS of saturated and unsaturated slopes under the seepage flow condition. The results of slope stability analysis under saturated and unsaturated seepage flow indicate that the slope FOS calculated by Mohr-Coulomb strength theory is less than 1.0, which is not consistent with the results of the slope model test(2)The saturated and unsaturated seepage calculation results show that the slope will have multiple sliding surfaces when the associated flow is used (). The comparison of slope stability analysis results under saturated and unsaturated seepage conditions shows that the shear dilatancy (flow rule) of the unsaturated soil has little influence on the slope stability, and it is mainly related to the unsaturated seepage mechanical behavior of the soil. In the case of saturated seepage, the soil’s dilatancy (flow rule) significantly influences the slope stability under the condition of seepage, so the soil’s dilatancy should be considered in the slope stability analysis under the condition of saturation(3)The loess area’s arid climate often causes a significant rise in the groundwater level due to the large-scale irrigation and seasonal rainfall. According to the numerical results, it is necessary to reduce the reclamation and irrigation of farmland along the loess plateau’s boundary. The stability of the high loess slope can be enhanced
Data Availability
The FEM calculation and the graph 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.
Acknowledgments
This research described in this paper was funded by the National Natural Science Foundation of China (No. 51979225, 51308456).