# Slope Stability Analysis under Complex Stress State with Saturated and Unsaturated Seepage Flow

Authors: Zongyuan Ma; Caihui Zhu; Xiaoliang Yao; Faning Dang
Journal: Geofluids
 (2021)

---

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


---

## Body

## 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 c 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]
(1)f=σ1−α1+bbσ2+σ3−σt,whenσ2≤σ1+ασ31+α,f′=11+bσ1+bσ2−ασ3−σt,whenσ2≥σ1+ασ31+α,where α is the ratio of tensile and compression yield strength, σt is the tensile yield strength, and b is the coefficient which reflects the influence of intermediate principal stress. The three principal stresses are denoted as σ1≥σ2≥σ3, 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 c and φ, which are widely used in geotechnical engineering; it follows that
(2)f=cosφ1+sinφ1+bbσ2+σ3−σ1+2C⋅cosφ1+sinφ,whenσ2≤1+sinφ2σ1+1−sinφ2σ3,f′=cosφ1+sinφ⋅σ3−11+bσ1+bσ2+2C⋅cosφ1+sinφ,whenσ2≥1+sinφ2σ1+1−sinφ2σ3.


### 2.2. Explicit Finite Difference Form

From the incremental elastoplastic theory, the total strain can be written as(3)Δεi=Δεie+Δεip.The elastic portion follows increment generalized Hooke’s law:(4)Δσ1=α1Δε1e+α2Δε2e+Δε3e,Δσ2=α1Δε2e+α2Δε1e+Δε3e,Δσ3=α1Δε3e+α2Δε1e+Δε2e,where α1=K+4G/3, α2=K−2G/3, K is bulk modulus, and G is shear modulus. The tensor form of Eq. (4) is given:
(5)Δσi=SiΔε1e,Δε2e,Δε3e=SiΔεne,i=1,2,3.From the plastic flow rule, the plastic incremental deformation part can be determined:(6)Δεip=λ∂g∂σi,where g 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 Δσn=σ1,σ2,σ3 loaded to a structure when fσn=0, and the state of the structure has no change (fσn+Δσn=0). From Eq. (5), the component of increment stress Δσn can be written as
(7)Δσi=SiΔεn−Δεnp.Substituting Eq. (6) into Eq. (7) and simplifying
(8)Δσi=SiΔεn−λ∂g∂σn=SiΔεn−λ⋅Si∂g∂σn,where λ·Si∂g/∂σn is the corresponding stress of plastic strain. fσn+Δσn can also be written as
(9)fσn+Δσn=fσn+f∗σn,where fσn=0,f∗σn=fΔσn−f0, and then, fσn+Δσn follows:
(10)fσn+Δσn=f∗Δσn.Substituting Eq. (8) into Eq. (10) and then fσn+Δσn can be deduced:
(11)fσn+Δσn=f∗SnΔεn−λ⋅f∗Sn∂g∂σn=0.In order to process the material nonlinearity, the stressσiN and σiI is defined as follows [22]:
(12)σiN=σi+Δσi=σi+SiΔεn−SiΔεnp,(13)σiI=σi+SiΔεn,where Δεn is the total strain of the elastic stage, σiI is a stress tensor in the elastic stage, and σiN is a stress tensor in the plastic stage. fσnI can be written as
(14)fσnI=fσn+SnΔεn=f∗SnΔεn.Substituting Eq. (14) into Eq. (11) and then λ can be deduced:
(15)λ=fσnIf∗Sn∂g/∂σn.From Eqs. (6), (12), and (13), the expression σiN can be deduced:
(16)σiN=σiI−λ⋅Si∂g∂σn.Eq. (2) can also be written as
(17)f=11+bNφbσ2+σ3−σ1+2cNφ,whenσ2≤1+sinφ2σ1+1−sinφ2σ3,f′=σ3Nφ−11+bσ1+bσ2+2cNφ,whenσ2≥1+sinφ2σ1+1−sinφ2σ3,where Nφ=1+sinφ/1−sinφ. 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 Nψ=1+sinψ/1−sinψ, and the item of λ·Si∂g/∂σn can be determined:
(18)S1λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα2Nψ−α1,S2λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λb⋅α11+bNψ+α211+bNψ−1,S3λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα11+bNψ+α2b1+bNψ−10,whenσ2≤1+sinϕ2σ1+1−sinϕ2σ3,S1λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα21Nψ−b1+b−α11+b,S2λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα21Nψ−11+b−α1⋅b1+b,S3λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα11Nψ−α2,whenσ2≥1+sinϕ2σ1+1−sinϕ2σ3.Substituting Eq. (18) into Eq. (16), the component of σiN yields
(19)σ1N=σ1I−λα21Nψ−α1,σ2N=σ2I−λα1b1+bNψ+α211+bNψ−1,σ3N=σ3I−λα111+bNψ+α2b1+bNψ−1,whenσ2≤1+sinϕ2σ1+1−sinϕ2σ3,σ1N=σ1I−λα21Nψ−b1+b−α11+b,σ2N=σ2I−λα21Nψ−11+b−α1⋅b1+b,σ3N=σ31I−λα11Nψ−α2,whenσ2≥1+sinϕ2σ1+1−sinϕ2σ3,where the expression of λ is as follows:
(20)λ=fσ1I,σ2I,σ3IfS1∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S2∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S3∂g/∂σ1,∂g/∂σ2,∂g/∂σ3−f0,whenσ2≤1+sinφ2σ1+1−sinφ2σ3,λ=f′σ1I,σ2I,σ3If′fS1∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S2∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S3∂g/∂σ1,∂g/∂σ2,∂g/∂σ3−f′0,whenσ2≥1+sinφ2σ1+1−sinφ2σ3.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 (c and φ). The reduced strength parameters (cf and φf) are given by
(21)cf=cSRF,φf=arctantanφSRF,where c and φ are the original cohesion and internal friction angle of the slope soil, cf and φf 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 Table1. 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 (b=1.0) when the slope water level drops to L/H=1.0 under saturated seepage condition. Figure 2 shows that the nonassociated flow (ψ=0) calculation results show that there are multiple sliding surfaces of the slope. However, the results obtained by associated flow (ψ=φ) show only one sliding surface.Table 1
Analysis and calculation parameters of slope stability under seepage flow condition.
Saturated seepage flowUnsaturated seepage flowPermeabilityk (m/d)Porosityn (%)Residual volume water contentθrSaturated volume water contentθsαm0.6500.0812 sorption0.4155 sorption0.0573 sorption4.0974 sorption0.1028 desorption0.4465 desorption0.0325 desorption7.4845 desorptionFigure 1
Relation between slope FOS and water levelL/H during the falling of water level.Figure 2
Maximum shear strain counters obtained by TSSTb=1.0 when the water level drops to L/H=1.0.
(a)
Nonassociated flowψ=0(b)
Associated flowψ=φ


### 2.4. Slope Stability due to Unsaturated Seepage Flow

The effect stressσb can be formulated as follows using Bishop’s unsaturated effect stress theory (compression is negative) [24]:
(22)σb=σ−SwPw+SaPa,where σ is the total stress, Sw is the saturation of water, Sa=1−Sw 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 (τ12 and τ23) and normal stress (σ12 and σ23) as follows [18, 19]:
(23)τ13+bτ12=C−βσ13b+bσ12b,whenτ12+βσ12b≥τ23+βσ23b,τ13+bτ23=C−βσ13b+bσ23b,whenτ12+βσ12b≤τ23+βσ23b,where β and C is the strength parameters and can be determined by cohesion c and friction angle φ (β=sinφ, C=2ccosφ); TSST also can be formulated by c and φ as follows:
(24)τ13+bτ12=2ccosφ−σ13b+bσ12bsinφ,whenτ12+σ12bsinφ≥τ23+σ23bsinφ,τ13+bτ23=2ccosφ−σ13b+bσ23bsinφ,whenτ12+σ12bsinφ≤τ23+σ23bsinφ,where the normal stress σ12 and σ23 under the unsaturated condition is as follows:
(25)σ13b=σ13−SwPw+SaPa=σ13−Pa−SwPa+Pw,σ12b=σ12−SwPw+SaPa=σ12−Pa−SwPa+Pw,σ23b=σ23−SwPw+SaPa=σ23−Pa−SwPa+Pw.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:
(26)τ13+bτ12=2ccosφ−σ13−Pa+bσ12−Pasinφ+1+bSwPa−Pwsinφ,whenτ12+σ12−Pasinφ≥τ23+σ23−Pasinφ,τ13+bτ23=2ccosφ−σ13−Pa+bσ23−Pasinφ+1+bSwPa−Pwsinφ,whenτ12+σ12−Pasinφ≤τ23+σ23−Pasinφ.The principal shear stress (τ13, τ12, and τ23) and normal stress (σ13, σ12, and σ23) can be expressed as follows:
(27)τ13=σ1−σ32,τ12=σ1−σ22,τ23=σ2−σ32,σ13=σ1+σ32,σ12=σ1+σ22,σ23=σ2+σ32.The shear strength formula of unsaturated soil also can be formulated via substituting Eq. (27) into Eq. (26):
(28)1+bσ1−bσ2−σ3=4ccosφ−1+bσ1+bσ2+σ3−21+bPasinφ+21+bSwPa−Pwsinφ,whenσ1−σ2+σ1+σ2−2Pasinφ≥σ2−σ3+σ2+σ3−2Pasinφ,σ1+bσ2−1+bσ3=4ccosφ−σ1+bσ2+1+bσ3−21+bPasinφ+21+bSwPa−Pwsinφ,whenσ1−σ2+σ1+σ2−2Pasinφ≤σ2−σ3+σ2+σ3−2Pasinφ.The variation of suctionPa−Pw of the unsaturated soil can be described by the empirical formula recommended by van Genuchten [25–27]:
(29)θw=θr+θs−θr1+αPa−Pwm1−1/m,where θw is the volume water content of the soil, R is the residual volume water content of the soil, and m is the measured parameters. The relationship between saturation Sw and volume water content θw is θw=Swn, where n 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 b; 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 (b=1.0) 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 L/H=1.0 yield by Mohr-Coulomb and TSST. The calculation result yield by TSST (b=1.0) in the associated flow (ψ=0) 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.Figure 3
Soil-water characteristic curves.Figure 4
Relation between FOS and parameterb when the underground water level rises to the slope crest.Figure 5
Maximum shear strain contours of slope (unsaturated seepage) when the water level drops toL/H=1.0.
(a)
TSSTb=0.0,ψ=0(b)
TSSTb=1.0,ψ=0(c)
TSSTb=0.0,ψ=φ(d)
TSSTb=1.0,ψ=φ


### 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 Figure6). 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°.Figure 6
Location of the high loess slope and terrain of Jingyang loess plateau.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 Figure7). 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.Figure 7
General view of landslide and cross-sections.Figure 8
Sketch of the cross-section of D-D′.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 indexIp 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.Table 2
Physical properties of soil for high loess slope.
Densityρ (kg/m3)Specific gravitydsWater contentwn (%)Void ratioePlasticity indexIp15402.7217.81.088.6


#### 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 h 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 x-axis is recorded during the calculation process to analyze the slope stability. The calculation parameters are as follows: Young’s modulus E=5.0 MPa, Poisson’s ratio ν=0.33, and dilation angle ψ=10°. The calculations are performed under the groundwater level in 1976 (h=80 m) and 1992 (h=37 m) to compare the slope stability under different groundwater levels, respectively. The relationship between the displacement at node O in the x-direction and the unbalanced force rate with the different magnitude of b 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 1.0×10−5.Figure 9
Meshing model of slope before landslide event.Figure 10
Relationship betweenx-displacement of node O and unbalanced force rate with different b.Figure10 indicated that the slope changes towards instability with the decrease of the parameter b, and the slope failure occurs when the magnitude of b 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 b 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 b, and the values of FOS also decrease with the rising level of the groundwater. The magnitudes of b 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 h=37 m and b=0.5. As shown in Figure 12, the high loess slope becomes unstable when h=37 m and b=0.5.Figure 11
Relationship between FOS andb under the groundwater levels in 1976 and 1992.Figure 12
Pore-water pressure and displacement vector of slope inh=37 m and b=0.5.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 (h=35-50 m) of the groundwater level. Figure 13 shows the relationship between the displacement of node O in the x-direction and unbalanced force rate under different depths of the groundwater level when b=0.5. The results suggest that the displacement of node O in the x-direction increases significantly as the groundwater rises gradually. The slope failure occurred until the depth of the groundwater level is close to 35 m.Figure 13
Relationship between displacement of node O inx-direction and unbalanced force rate (b=0.5).


#### 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 depthh=80 m, the calculation of seepage flow is performed with four scenarios: the rise of groundwater level Δh=30 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 thex-direction increases significantly when the unbalanced force rate is close to the target of convergence and the slope failure occurs suddenly when Δh 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 x-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.


## 2.1. Twin Shear Strength Theory

The principal stress form of TSST can be written as [18, 19]
(1)f=σ1−α1+bbσ2+σ3−σt,whenσ2≤σ1+ασ31+α,f′=11+bσ1+bσ2−ασ3−σt,whenσ2≥σ1+ασ31+α,where α is the ratio of tensile and compression yield strength, σt is the tensile yield strength, and b is the coefficient which reflects the influence of intermediate principal stress. The three principal stresses are denoted as σ1≥σ2≥σ3, 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 c and φ, which are widely used in geotechnical engineering; it follows that
(2)f=cosφ1+sinφ1+bbσ2+σ3−σ1+2C⋅cosφ1+sinφ,whenσ2≤1+sinφ2σ1+1−sinφ2σ3,f′=cosφ1+sinφ⋅σ3−11+bσ1+bσ2+2C⋅cosφ1+sinφ,whenσ2≥1+sinφ2σ1+1−sinφ2σ3.


## 2.2. Explicit Finite Difference Form

From the incremental elastoplastic theory, the total strain can be written as(3)Δεi=Δεie+Δεip.The elastic portion follows increment generalized Hooke’s law:(4)Δσ1=α1Δε1e+α2Δε2e+Δε3e,Δσ2=α1Δε2e+α2Δε1e+Δε3e,Δσ3=α1Δε3e+α2Δε1e+Δε2e,where α1=K+4G/3, α2=K−2G/3, K is bulk modulus, and G is shear modulus. The tensor form of Eq. (4) is given:
(5)Δσi=SiΔε1e,Δε2e,Δε3e=SiΔεne,i=1,2,3.From the plastic flow rule, the plastic incremental deformation part can be determined:(6)Δεip=λ∂g∂σi,where g 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 Δσn=σ1,σ2,σ3 loaded to a structure when fσn=0, and the state of the structure has no change (fσn+Δσn=0). From Eq. (5), the component of increment stress Δσn can be written as
(7)Δσi=SiΔεn−Δεnp.Substituting Eq. (6) into Eq. (7) and simplifying
(8)Δσi=SiΔεn−λ∂g∂σn=SiΔεn−λ⋅Si∂g∂σn,where λ·Si∂g/∂σn is the corresponding stress of plastic strain. fσn+Δσn can also be written as
(9)fσn+Δσn=fσn+f∗σn,where fσn=0,f∗σn=fΔσn−f0, and then, fσn+Δσn follows:
(10)fσn+Δσn=f∗Δσn.Substituting Eq. (8) into Eq. (10) and then fσn+Δσn can be deduced:
(11)fσn+Δσn=f∗SnΔεn−λ⋅f∗Sn∂g∂σn=0.In order to process the material nonlinearity, the stressσiN and σiI is defined as follows [22]:
(12)σiN=σi+Δσi=σi+SiΔεn−SiΔεnp,(13)σiI=σi+SiΔεn,where Δεn is the total strain of the elastic stage, σiI is a stress tensor in the elastic stage, and σiN is a stress tensor in the plastic stage. fσnI can be written as
(14)fσnI=fσn+SnΔεn=f∗SnΔεn.Substituting Eq. (14) into Eq. (11) and then λ can be deduced:
(15)λ=fσnIf∗Sn∂g/∂σn.From Eqs. (6), (12), and (13), the expression σiN can be deduced:
(16)σiN=σiI−λ⋅Si∂g∂σn.Eq. (2) can also be written as
(17)f=11+bNφbσ2+σ3−σ1+2cNφ,whenσ2≤1+sinφ2σ1+1−sinφ2σ3,f′=σ3Nφ−11+bσ1+bσ2+2cNφ,whenσ2≥1+sinφ2σ1+1−sinφ2σ3,where Nφ=1+sinφ/1−sinφ. 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 Nψ=1+sinψ/1−sinψ, and the item of λ·Si∂g/∂σn can be determined:
(18)S1λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα2Nψ−α1,S2λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λb⋅α11+bNψ+α211+bNψ−1,S3λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα11+bNψ+α2b1+bNψ−10,whenσ2≤1+sinϕ2σ1+1−sinϕ2σ3,S1λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα21Nψ−b1+b−α11+b,S2λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα21Nψ−11+b−α1⋅b1+b,S3λ∂g∂σ1,λ∂g∂σ2,λ∂g∂σ3=λα11Nψ−α2,whenσ2≥1+sinϕ2σ1+1−sinϕ2σ3.Substituting Eq. (18) into Eq. (16), the component of σiN yields
(19)σ1N=σ1I−λα21Nψ−α1,σ2N=σ2I−λα1b1+bNψ+α211+bNψ−1,σ3N=σ3I−λα111+bNψ+α2b1+bNψ−1,whenσ2≤1+sinϕ2σ1+1−sinϕ2σ3,σ1N=σ1I−λα21Nψ−b1+b−α11+b,σ2N=σ2I−λα21Nψ−11+b−α1⋅b1+b,σ3N=σ31I−λα11Nψ−α2,whenσ2≥1+sinϕ2σ1+1−sinϕ2σ3,where the expression of λ is as follows:
(20)λ=fσ1I,σ2I,σ3IfS1∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S2∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S3∂g/∂σ1,∂g/∂σ2,∂g/∂σ3−f0,whenσ2≤1+sinφ2σ1+1−sinφ2σ3,λ=f′σ1I,σ2I,σ3If′fS1∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S2∂g/∂σ1,∂g/∂σ2,∂g/∂σ3,S3∂g/∂σ1,∂g/∂σ2,∂g/∂σ3−f′0,whenσ2≥1+sinφ2σ1+1−sinφ2σ3.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 (c and φ). The reduced strength parameters (cf and φf) are given by
(21)cf=cSRF,φf=arctantanφSRF,where c and φ are the original cohesion and internal friction angle of the slope soil, cf and φf 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 Table1. 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 (b=1.0) when the slope water level drops to L/H=1.0 under saturated seepage condition. Figure 2 shows that the nonassociated flow (ψ=0) calculation results show that there are multiple sliding surfaces of the slope. However, the results obtained by associated flow (ψ=φ) show only one sliding surface.Table 1
Analysis and calculation parameters of slope stability under seepage flow condition.
Saturated seepage flowUnsaturated seepage flowPermeabilityk (m/d)Porosityn (%)Residual volume water contentθrSaturated volume water contentθsαm0.6500.0812 sorption0.4155 sorption0.0573 sorption4.0974 sorption0.1028 desorption0.4465 desorption0.0325 desorption7.4845 desorptionFigure 1
Relation between slope FOS and water levelL/H during the falling of water level.Figure 2
Maximum shear strain counters obtained by TSSTb=1.0 when the water level drops to L/H=1.0.
(a)
Nonassociated flowψ=0(b)
Associated flowψ=φ


## 2.4. Slope Stability due to Unsaturated Seepage Flow

The effect stressσb can be formulated as follows using Bishop’s unsaturated effect stress theory (compression is negative) [24]:
(22)σb=σ−SwPw+SaPa,where σ is the total stress, Sw is the saturation of water, Sa=1−Sw 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 (τ12 and τ23) and normal stress (σ12 and σ23) as follows [18, 19]:
(23)τ13+bτ12=C−βσ13b+bσ12b,whenτ12+βσ12b≥τ23+βσ23b,τ13+bτ23=C−βσ13b+bσ23b,whenτ12+βσ12b≤τ23+βσ23b,where β and C is the strength parameters and can be determined by cohesion c and friction angle φ (β=sinφ, C=2ccosφ); TSST also can be formulated by c and φ as follows:
(24)τ13+bτ12=2ccosφ−σ13b+bσ12bsinφ,whenτ12+σ12bsinφ≥τ23+σ23bsinφ,τ13+bτ23=2ccosφ−σ13b+bσ23bsinφ,whenτ12+σ12bsinφ≤τ23+σ23bsinφ,where the normal stress σ12 and σ23 under the unsaturated condition is as follows:
(25)σ13b=σ13−SwPw+SaPa=σ13−Pa−SwPa+Pw,σ12b=σ12−SwPw+SaPa=σ12−Pa−SwPa+Pw,σ23b=σ23−SwPw+SaPa=σ23−Pa−SwPa+Pw.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:
(26)τ13+bτ12=2ccosφ−σ13−Pa+bσ12−Pasinφ+1+bSwPa−Pwsinφ,whenτ12+σ12−Pasinφ≥τ23+σ23−Pasinφ,τ13+bτ23=2ccosφ−σ13−Pa+bσ23−Pasinφ+1+bSwPa−Pwsinφ,whenτ12+σ12−Pasinφ≤τ23+σ23−Pasinφ.The principal shear stress (τ13, τ12, and τ23) and normal stress (σ13, σ12, and σ23) can be expressed as follows:
(27)τ13=σ1−σ32,τ12=σ1−σ22,τ23=σ2−σ32,σ13=σ1+σ32,σ12=σ1+σ22,σ23=σ2+σ32.The shear strength formula of unsaturated soil also can be formulated via substituting Eq. (27) into Eq. (26):
(28)1+bσ1−bσ2−σ3=4ccosφ−1+bσ1+bσ2+σ3−21+bPasinφ+21+bSwPa−Pwsinφ,whenσ1−σ2+σ1+σ2−2Pasinφ≥σ2−σ3+σ2+σ3−2Pasinφ,σ1+bσ2−1+bσ3=4ccosφ−σ1+bσ2+1+bσ3−21+bPasinφ+21+bSwPa−Pwsinφ,whenσ1−σ2+σ1+σ2−2Pasinφ≤σ2−σ3+σ2+σ3−2Pasinφ.The variation of suctionPa−Pw of the unsaturated soil can be described by the empirical formula recommended by van Genuchten [25–27]:
(29)θw=θr+θs−θr1+αPa−Pwm1−1/m,where θw is the volume water content of the soil, R is the residual volume water content of the soil, and m is the measured parameters. The relationship between saturation Sw and volume water content θw is θw=Swn, where n 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 b; 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 (b=1.0) 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 L/H=1.0 yield by Mohr-Coulomb and TSST. The calculation result yield by TSST (b=1.0) in the associated flow (ψ=0) 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.Figure 3
Soil-water characteristic curves.Figure 4
Relation between FOS and parameterb when the underground water level rises to the slope crest.Figure 5
Maximum shear strain contours of slope (unsaturated seepage) when the water level drops toL/H=1.0.
(a)
TSSTb=0.0,ψ=0(b)
TSSTb=1.0,ψ=0(c)
TSSTb=0.0,ψ=φ(d)
TSSTb=1.0,ψ=φ


## 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 Figure6). 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°.Figure 6
Location of the high loess slope and terrain of Jingyang loess plateau.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 Figure7). 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.Figure 7
General view of landslide and cross-sections.Figure 8
Sketch of the cross-section of D-D′.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 indexIp 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.Table 2
Physical properties of soil for high loess slope.
Densityρ (kg/m3)Specific gravitydsWater contentwn (%)Void ratioePlasticity indexIp15402.7217.81.088.6


### 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 h 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 x-axis is recorded during the calculation process to analyze the slope stability. The calculation parameters are as follows: Young’s modulus E=5.0 MPa, Poisson’s ratio ν=0.33, and dilation angle ψ=10°. The calculations are performed under the groundwater level in 1976 (h=80 m) and 1992 (h=37 m) to compare the slope stability under different groundwater levels, respectively. The relationship between the displacement at node O in the x-direction and the unbalanced force rate with the different magnitude of b 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 1.0×10−5.Figure 9
Meshing model of slope before landslide event.Figure 10
Relationship betweenx-displacement of node O and unbalanced force rate with different b.Figure10 indicated that the slope changes towards instability with the decrease of the parameter b, and the slope failure occurs when the magnitude of b 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 b 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 b, and the values of FOS also decrease with the rising level of the groundwater. The magnitudes of b 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 h=37 m and b=0.5. As shown in Figure 12, the high loess slope becomes unstable when h=37 m and b=0.5.Figure 11
Relationship between FOS andb under the groundwater levels in 1976 and 1992.Figure 12
Pore-water pressure and displacement vector of slope inh=37 m and b=0.5.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 (h=35-50 m) of the groundwater level. Figure 13 shows the relationship between the displacement of node O in the x-direction and unbalanced force rate under different depths of the groundwater level when b=0.5. The results suggest that the displacement of node O in the x-direction increases significantly as the groundwater rises gradually. The slope failure occurred until the depth of the groundwater level is close to 35 m.Figure 13
Relationship between displacement of node O inx-direction and unbalanced force rate (b=0.5).


### 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 depthh=80 m, the calculation of seepage flow is performed with four scenarios: the rise of groundwater level Δh=30 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 thex-direction increases significantly when the unbalanced force rate is close to the target of convergence and the slope failure occurs suddenly when Δh 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 x-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.


## 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 h 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 x-axis is recorded during the calculation process to analyze the slope stability. The calculation parameters are as follows: Young’s modulus E=5.0 MPa, Poisson’s ratio ν=0.33, and dilation angle ψ=10°. The calculations are performed under the groundwater level in 1976 (h=80 m) and 1992 (h=37 m) to compare the slope stability under different groundwater levels, respectively. The relationship between the displacement at node O in the x-direction and the unbalanced force rate with the different magnitude of b 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 1.0×10−5.Figure 9
Meshing model of slope before landslide event.Figure 10
Relationship betweenx-displacement of node O and unbalanced force rate with different b.Figure10 indicated that the slope changes towards instability with the decrease of the parameter b, and the slope failure occurs when the magnitude of b 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 b 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 b, and the values of FOS also decrease with the rising level of the groundwater. The magnitudes of b 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 h=37 m and b=0.5. As shown in Figure 12, the high loess slope becomes unstable when h=37 m and b=0.5.Figure 11
Relationship between FOS andb under the groundwater levels in 1976 and 1992.Figure 12
Pore-water pressure and displacement vector of slope inh=37 m and b=0.5.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 (h=35-50 m) of the groundwater level. Figure 13 shows the relationship between the displacement of node O in the x-direction and unbalanced force rate under different depths of the groundwater level when b=0.5. The results suggest that the displacement of node O in the x-direction increases significantly as the groundwater rises gradually. The slope failure occurred until the depth of the groundwater level is close to 35 m.Figure 13
Relationship between displacement of node O inx-direction and unbalanced force rate (b=0.5).


## 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 depthh=80 m, the calculation of seepage flow is performed with four scenarios: the rise of groundwater level Δh=30 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 thex-direction increases significantly when the unbalanced force rate is close to the target of convergence and the slope failure occurs suddenly when Δh 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 x-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 (ψ=0). 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



---

*Source: 6637098-2021-01-11.xml*