# Analytical Computation of Sandy Slope Stability under the Groundwater Seepage and Earthquake Conditions

Authors: Yulin Lu; Xiaoran Chen; Li Wang
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Groundwater and earthquake are two of the factors behind the occurrence of slope instability. This paper focuses on the computation of the safety factor acting on a sandy slope when subjected to groundwater seepage and earthquake load. Based on the theory of limit equilibrium and a pseudodynamic method of analysis, a general solution for the dynamic factor of sandy slope containing two load conditions is proposed. In the solution, the effects of parameters such as water level, seismic acceleration coefficients, slope angle, and soil internal friction angle on the slope stability have been discussed. In order to evaluate the validity of the present formulation, the dynamic factors of safety computed by the present method are compared with those given by other cases. It is found that the results of the present method are in good agreement with those of the previous method. The results indicate that the dynamic of safety factor increases linearly with an increase in the soil internal friction angle but decreases with the slope angle and groundwater level and that the effects of direction and magnitude of the seismic inertia forces on the dynamic safety factor of sandy slope are different at different times.


---

## Body

## 1. Introduction

The effects of groundwater and earthquake on slope are a very important topic in geotechnical engineering. Therefore, the knowledge of slope stability is very essential under water seepage and seismic loading conditions [1–3]. The failure surface of a slope can be simplified as planar, arc, logarithmic spiral surface, and so on [4–6]. This paper presents an expression for the safety factor of slope with a planar failure surface subjected simultaneously to seismic forces and variable groundwater conditions, with the aim of emphasizing the fundamental importance of the complex conditions in slope stability analysis.Usually, groundwater is a factor behind the occurrence of slope instability, and the most unstable part of a sandy slope is the groundwater discharge zone, where the effect of seepage force occurs [7–10]. The critical role that the groundwater plays in the stability of slope was also recognized by Terzaghi [11] in the effective stress law. This instability is frequently associated with saturation of soil slopes, which increases the weight of soil and reduces the effective stress. A number of researchers have estimated the local stability of slopes under the effect of a seepage flow. Deere and Patton and Patton and Hendron discussed the implications of groundwater flow regimes for slope stability in a number of different slope environments [12–14]. The effects of groundwater on slope stability can also be seen by considering the slope under different degrees of submergence and drawdown. Lane and Griffiths introduced a finite-element method to the assessment of stability of slopes under drawdown conditions [15]. Nian et al. proposed a program based on the shear strength reduction finite-element method to compute the safety factor for the slopes under the joint action of surface water and groundwater [16]. Gao et al. investigated the influence of water drawdown on the stability of 3D slopes based on the kinematic approach of limit analysis [17]. Lu and Chen analyzed the slope stability subjected to the seepage of groundwater by using the simplified Bishop method [18]. Xue et al. also used the sequentially coupled method of seepage-softening-stability to analyze the loess slope stability with a rising groundwater level [19].Another important cause of slope instability is seismic shaking. Previously, analytical methods of slope stability include the pseudostatic method, time history analysis method, and Newmark’s sliding block method. Zhang (2020) used the time history analysis method and the Newmark method to evaluate the seismic stability of different slope shape under different amplitudes [20]. Ji et al. also proposed a simplified rotational permanent displacement calculation framework based on Newmark’s sliding block theory for seismic slope stability assessment [21]. However, the pseudostatic method used to be the most common method for evaluating the stability of landslides under seismic loads, and this pioneering work was first reported by Terzaghi et al. [22]. In the pseudostatic approach, the horizontal or vertical acceleration generated by earthquake shaking is considered to create seismic inertia force. Then, the pseudostatic safety factor is obtained by using the limit equilibrium method, which is related to the material behaviors of the sliding body, the shape of failure surface of slope, and the location and magnitude of earthquake. In the last decade, several researchers, such as Loukidis et al.; Choudhury et al.; Zhou et al., also implemented the pseudostatic method to analyze the stability of slope [23–25].Though the pseudostatic analysis is simple and straightforward, the representation of complex, dynamic effect of earthquake shaking by constant horizontal or vertical acceleration is actually pretty crude. The dynamic behavior of earthquake is considered only as time independent in this method, but it does not consider the effects of time, frequency, and body wave traveling the soil during the earthquake. To overcome this drawback, Steedman and Zeng introduced a pseudodynamic approach that considers the effect of time and phase difference due to the wave propagation behind the retaining wall [26]. One merit of the pseudodynamic method is that horizontal and vertical seismic displacement and finite shear and primary wave propagations through the soil with time variation can be taken into account. Another merit is its simplicity, with the ability to provide analytical solutions, which is convenient for analysis of the stability of slope. Recently, Zhou et al. and Hazari et al. proposed the pseudodynamic limit equilibrium method to analyze the stability of slope subjected to seismic loads [25, 27].Although there are many studies on the stability of slope under water seepage or earthquake, the evaluation of coupling effects on slope instability is rare at present. Loáiciga demonstrated the application of screening analysis for seismic slope stability through the results of numerous examples, providing a foundation to assist in designs of slope modification for groundwater control and to protect against earthquakes [28]. Song et al. analyzed the stability of sandy soil slope under the coupling of earthquake and groundwater based on the finite-element model of sandy soil slope [29]. Lai et al. used the strength reduction dynamic analysis method to analyze the dynamic stability of an artificial island under the coupled actions of gravity, seepage, and earthquake [30]. Obviously, most of the above research methods are numerical simulation and model test. However, the use of the pseudodynamic approach to evaluate the stability of slope under the groundwater seepage and seismic conditions in the framework of the limit equilibrium method has secured little attention in the past.In this paper, a pseudodynamic model is proposed to consider the effects of groundwater and earthquake on the stability of sandy slope, in which the horizontal and vertical seismic inertia forces and pore water pressure induced by earthquake acting on the potential slip body are taken into account. The limit equilibrium method, with a planar failure surface of sandy slope, is considered in the analysis. This study also provides a formulation for calculating the pseudodynamic factors of safety of slopes based on the groundwater seepage and earthquake shaking conditions and analyzes the effect of various parameters on the stability of slope.


## 2. Failure Mechanism of the Slope with Water Seepage and Earthquake

The stability analysis of slope based on the limit equilibrium theory is safe; as long as the resisting forces exceed the driving forces on the slip surface, the static equilibrium can be maintained. When the slope is subjected to groundwater seepage and earthquake shaking, the seismic forces and excess pore water pressure influence both the driving and resisting forces, leading to potential situations in which the driving forces at least temporarily exceed the resisting forces. Such situation will lead to relative motion between the slip body and the plane. In order to get a better understanding of the influence of the failure mechanism, the soil stress state with water seepage and earthquake is analyzed. It is generally considered that the stress state of the soil in Mohr–Coulomb theory is the total stress, as shown in Figure1(a). σ1 is the maximal principal stress, and σ3 is the minimum principal stress. When the water seepage is considered separately, the increase of pore water pressure will result in the decrease of effective stress of soil, which makes the Mohr circle of effective stress closer to the strength envelope of soil, and the soil is more likely to lose stability, as shown in Figure 1(b). The occurrence of earthquake will also aggravate the instability of slope. An actual slope is mostly affected by the double action of seismic force and water seepage force, and then the effective stress of the soil will decrease sharply under this complex stress condition, as shown in Figure 1(c).Figure 1
Stress state of soil. (a) Total stress state. (b) Single water seepage. (c) Double action of water seepage and earthquake.
(a)(b)(c)Where the pore water pressure isp, the horizontal and vertical stresses of the earthquake action are σh and σv, the effective stresses of soil are σ1′ and σ3′. It can be seen that the slope is prone to instability under the combined action of lateral seismic force or water pressure and seismic force, which is due to the significant loss of effective stress of soil.


## 3. Analytical Method of Sandy Slope with Water Seepage and Earthquake

### 3.1. Water Pressure Induced by Water Seepage

The typical instability of sandy slope is plane shear failure; thus, it is assumed that the slip surface is a planar surface, and the slip body is treated as a rigid body as shown in Figure2. A Cartesian coordinate system is established, which is centered at the toe of the slope A (0, 0). The height of slope is H, the angle of slope is θ, the angle of the failure plane with the horizontal ground surface is α, the unit weight of soil is γ, the saturation weight of soil is γsat, and the unit weight of water is W. The dotted line AM in the slope model is the phreatic surface, and the height of the initial groundwater level from the bottom is h. The outlet of the phreatic surface is located at the bottom of the slope, and hw is the intersection of phreatic surface and slip surface.Figure 2
Sandy slope calculated model.The slip body is divided inton slices, and the forces acting on any slice are plotted in Figure 3. The weight of soil above phreatic surface is dW1, the saturated weight of soil below the phreatic surface is dW2, the water pressure force of side AB is Pa, and the water pressure force of side CD is Pb. The water pressure force of side BC is dU, the normal force is dN, and the tangential force is dT. The angle of the failure plane with the horizontal ground surface is α, the angle of the phreatic surface with the horizontal ground surface is β, and the length of failure plane BC is dL.Figure 3
Forces acting on the slices subjected to seepage loads when the slice forces are neglected.The water pressure problem has been solved by Zheng et al. [10], and the water pressure force of BC is given by(1)dU=γwha+hbdL2cos2β.Equation (1) can be obtained by mathematical transformation as follows:(2)dU=γwha+hbdLcosα2cosαcos2β=γwdAwcos2βcosα=γwdAw1−sin2βcosα=γwdAw1−i2cosα.Therefore, the total water pressure force induced by the water seepage,U1, which is the result from the pore water as shown in equation (2), acting along the failure plane BC, is given by(3)U1=γwAw1−i2cosα,where the area of immersed soil is Aw, which is the area surrounded by the phreatic surface and failure plane. The hydraulic gradient is sinβ, and it can be simplified as average hydraulic gradient i.


### 3.2. Pore Water Pressure Induced by Earthquake

Under earthquake condition, an excess pore water pressure may be induced inside the slope. For a plane state, the excess pore water pressure △u can be obtained by a pseudostatic method of analysis and is given by Wang et al. [8] as follows:(4)Δu=γsatβkh+kv1+μ3+2αkh+kv2μ2−μ+1−3khkv,where μ is Poisson’s ratio of the soil, α and β are the coefficients of excess pore pressure, β = 1.0 for saturated soil, and kh and kv are the seismic acceleration coefficients in the horizontal and vertical directions.The total excess pore water pressure induced by the earthquake,U2, which is the result of the excess pore water, △u, shown in (4), acting alone the failure plane AF, is as follows:(5)U2=∫0hwsinαγsatβkh+kv1+μ3+2αkh+kv2μ2−μ+1−3khkvdz=γsathwβkh+kv1+μ3sinα+γsathw2αkh+kv2μ2−μ+1−3khkvsinα.


### 3.3. Force Analysis of Slip Body

The slip body is taken as the research object, as shown in Figure4. W is the total weight of the slip body, Eh is the horizontal seismic inertia force, Ev is the vertical seismic inertia force, U1 is the total water pressure caused by seepage, and U2 is the total excess pore water pressure caused by earthquake. N and S are the normal and tangential forces on the failure plane, respectively. The slip body is divided into three parts, containing AFA, AGFA, and GBCFG. The slip body is considered to be rigid in calculation, and the total force can be equivalent to the sum of three parts.Figure 4
Forces acting on the slip body.


#### 3.3.1. Calculation of Slip Body Weight

In the process of determining the phreatic surface, the physical approximations are the same as Dupuit assumed conditions, which are obtained by solving the Boussinesq equation. Based on the research results of seepage analysis, a parabola model can be used to calculate the phreatic surface induced by the groundwater seepage. This model has been applied to the dam, given by Zheng et al. [10], as shown in Figure 4. The phreatic surface equation is as follows:(6)z2=ax+b,where a and b are undetermined coefficients. The hydraulic radius of phreatic surface is R, and it is AE length as shown in Figure 2. It can be seen from Figure 2 that the coordinates of point A and point M are A (0, 0) and M (h, R), respectively, and both points are on the parabola line. So, the parameters a and b can be solved by taking point A and point B into (6) as follows:(7)b=0,a=h2R.At the same time, the water level at the intersection of the phreatic surface and slip surface is expressed as(8)hw=h2Rtanα.The hydraulic radiusR can be determined by the Gehart formula as follows:(9)R=10ΔHK,where △H is the water drawdown. K is the permeability coefficient. Then, the equation of phreatic surface is given by(10)z2=h2xR.It can be seen from Figure4 that the mathematical equation of failure plane AC can be expressed as(11)z=xtanα.Based on the above conclusions, the total weight of the slip body can be written as the sum of weight of geometryAFA, AGFA, and GBCFG.The weight of wedgeAFA is given by(12)WAFA=∫0h2/Rtanαγsatztanα−z2Rh2dz=γsath46R2tan3α.The weight of wedgeAGFA is given by(13)WAGFA=∫0h2/Rtanαγz2Rh2−ztanθdz=γh42tanθ−3tanα6R2tan3αtanθ.The weight of wedgeGBCFG is given by(14)WGBCFG=∫h2/RtanαHγcotα−cotθzdz=12γcotα−cotθH2−h4R2tan2α.The total weight of the slip body is given by(15)W=WAFA+WAGFA+WGBCFG.


#### 3.3.2. Calculation of Seismic Inertia Force of Slip Body

The seismic inertia force can be calculated by the pseudodynamic method, and the rationality of this method is proved by Bash and Babu [5]. It is assumed that the shear modulus is constant with depth through the soil and that only the phase, not the magnitude of acceleration, is varying. The major advantage of using sinusoidal acceleration in the pseudodynamic method is that it possesses most of the representative characters of the earthquake. Since the pseudodynamic method considers the effects of vibration characteristics, the pseudodynamic method will degenerate into the pseudostatic method when seismic wave velocity tends to infinity.The period of lateral shaking can be defined asT = 2π/w, where w is the angular frequency of the base shaking. In the present study, it is assumed that the slope is subjected to both horizontal and vertical sinusoidal acceleration with amplitude of acceleration khg and kvg, respectively, where g is the acceleration due to gravity. The variation of horizontal and vertical sinusoidal acceleration at any depth z below the top of the slope and time (t) with soil amplification factor (f) can be expressed as follows:(16)ahz,t=1+H−zHf−1khgsinwt−H−zvs,avz,t=1+H−zHf−1kvgsinwt−H−zvp,where vs is the shear wave velocity propagating through the soil and vp is the primary wave velocity propagating through the soil. ah is the horizontal seismic acceleration, and av is the vertical seismic acceleration. G is the shear modulus of soil, and g is the gravity acceleration. Based on the soil dynamics theory [3], the shear wave velocity and the primary wave velocity can be expressed as follows:(17)vs=Gρvp=G2−μρ1−2μ,where ρ is the density of soil mass and μ is Poisson’s ratio.The inertia force of different wedge can be calculated respectively. Because the value ofz in the pseudodynamic method is based on the free surface at the top of the slope, the value of z needs to be converted in Figure 4.In this study, the Cartesian coordinate system is established with the slope toe as the origin, as shown in Figure4. In order to simplify the calculation, the soil amplification factor f = 1 is assumed. Thus, the inertial forces of the slip body are calculated as follows:(1)
Horizontal and vertical inertial force acting on the wedgeAFA.The horizontal inertia force acting on the wedgeAFA can be written as follows:(18)Eh−AFA=∫0h2cotα/Rγsatgzcotα−z2Rh2khg⋅sinwt−zvsdz=m1γsatkhcoswt−h2cotαRvs+m2γsatkhsinwt−h2cotαRvs+m3γsatkhcoswt+m4γsatkhsinwt,where the coefficients can be expressed as follows:(19)m1=vs3w32Rh2,m2=−vs2w2cotα,m3=−vs3w32Rh2,m4=−vs2w2cotα.Similarly, the vertical inertia force acting on the wedgeAFA can be written as follows:(20)Ev−AFA=∫0h2cotα/Rγsatgzcotα−z2Rh2kvg⋅sinwt−zvpdz=m5γsatkvcoswt−h2cotαRvp+m6γsatkvsinwt−h2cotαRvp+m7γsatkvcoswt+m8γsatkvsinwt,where the coefficients can be expressed as follows:(21)m5=vp3w32Rh2,m6=−vp2w2cotα,m7=−vp3w32Rh2,m8=−vp2w2cotα.(2)
Horizontal and vertical inertial force acting on the wedgeAGFA.The horizontal inertia force acting on the wedgeAGFA can be written as follows:(22)Eh−AGFA=∫0h2cotα/Rγgz2Rh2−zcotθkhg⋅sinwt−zvsdz=m9γkhcoswt−h2cotαRvp+m10γkhsinwt−h2cotαRvp+m11γkhcoswt+m12γkhsinwt,where the coefficients can be expressed as follows:(23)m9=−vs3w32Rh2+vswh2cotαRcotα−cotθm10=vs2w22cotα−cotθm11=vs3w32Rh2m12=vs2w2cotθ.Similarly, the vertical inertia force acting on the wedge AGFA can be written as follows:(24)Ev−AGFA=∫0h2cotα/Rγgz2Rh2−zcotθkvg⋅sinwt−zvpdz=m13γkvcoswt−h2cotαRvp+m14γkvsinwt−h2cotαRvp+m15γkvcoswt+m16γkvsinwt,where the coefficients can be expressed as follows:(25)m13=−vp3w32Rh2+vpwh2cotαRcotα−cotθ,m14=vp2w22cotα−cotθ,m15=vp3w32Rh2,m16=vp2w2cotθ.(3)
Horizontal and vertical inertial force acting on the wedgeGBCFG.The horizontal inertia force acting on the wedge GBCFG can be written as follows:(26)Eh−GBCFG=∫h2cotα/RHγgcotα−cotθzkhgsinwt−zvsdz=m17γkhcoswt−h2cotαRvs+m18γkhsinwt−h2cotαRvs+m19γkhcoswt−Hvs+m20γkhsinwt−Hvs,where the coefficients can be expressed as follows:(27)m17=vswcotθ−cotαh2cotαR,m18=vs2w2cotθ−cotα,m19=vswcotθ−cotαH,m20=vs2w2cotθ−cotα.Similarly, the vertical inertia force acting on the wedgeGBCFG can be written as follows:(28)Ev−GBCFG=∫h2cotα/RHγgcotα−cotθzkvgsinwt−zvpdz=m21γkvcoswt−h2cotαRvp+m22γkvsinwt−h2cotαRvp+m23γkvcoswt−Hvp+m24γkvsinwt−Hvp,where the coefficients can be expressed as follows:(29)m21=vpwcotθ−cotαh2cotαR,m22=vp2w2cotθ−cotα,m23=vpwcotθ−cotαH,m24=vp2w2cotθ−cotα.Therefore, the total inertia force of the slip body is given by(30)Eh=Eh−AFA+Eh−AGFA+Eh−GBCFG,Ev=Ev−AFA+Ev−AGFA+Ev−GBCFG.


### 3.4. Stability Calculation of Slope

Limit equilibrium method is used to determine the safety factor against slope failure. The safety factor of the sandy slope can be expressed by the ratio of resisting force and sliding force on slip surface. The equilibrium equation of normal direction of slip surface (∑N = 0) is given by(31)N+U1+U2+Ehsinα−W−Evcosα=0.The equilibrium equation of tangential direction of slip surface (∑S = 0) is given by(32)Ehcosα+W−Evsinα−S=0.Theoretically, the cohesion of sandy soil is zero, but it is not zero in practice. In order to ensure the accuracy of the analytical calculation, the cohesion is assumed to exist in the resisting force calculation and then chosen according to its actual value. Therefore, the safety factor can be written as follows:(33)F=cl+W−Evcosα−Ehsinα−U1−U2tanφEhcosα+W−Evsinα,where c is the cohesion of soil mass and φ is the internal friction angle of soil mass.In (33), the relations of the terms b, l, α, and θ can be written as(34)b=Hcotα−cotθ,l=Hcscα.Based on the above work, when the seismic parameters, slope parameters, and hydraulic parameters are known, the safety factor is a function of the slip surface angleα and time t. The critical slip surface, whose safety factor is the minimum among all the available ones, can be determined by the objective function; therefore, the inclination angle extreme and the safety factor at each time can be obtained by the function extreme theory. The minimum safety factor of slope in the vibration time domain will be used to evaluate the stability of slope.It is generally believed that the cohesive force of sandy soil is 0 kPa. But in fact, the cohesive force of unsaturated sand and the cohesive force of silty sand are nonzero. The principle of selecting cohesive force is given by early researchers’ works [4]. If the cohesive force of sandy soil measured in the experiment is less than 6 kPa, it can be equivalent to pure sandy soil, and the calculated cohesive force is set to zero. When the cohesion is greater than 6 kPa, the cohesion force is generally reduced by half in calculation. The purpose of this processing in this study is to make better use of the classification experimental data of sandy soil for calculation.


## 3.1. Water Pressure Induced by Water Seepage

The typical instability of sandy slope is plane shear failure; thus, it is assumed that the slip surface is a planar surface, and the slip body is treated as a rigid body as shown in Figure2. A Cartesian coordinate system is established, which is centered at the toe of the slope A (0, 0). The height of slope is H, the angle of slope is θ, the angle of the failure plane with the horizontal ground surface is α, the unit weight of soil is γ, the saturation weight of soil is γsat, and the unit weight of water is W. The dotted line AM in the slope model is the phreatic surface, and the height of the initial groundwater level from the bottom is h. The outlet of the phreatic surface is located at the bottom of the slope, and hw is the intersection of phreatic surface and slip surface.Figure 2
Sandy slope calculated model.The slip body is divided inton slices, and the forces acting on any slice are plotted in Figure 3. The weight of soil above phreatic surface is dW1, the saturated weight of soil below the phreatic surface is dW2, the water pressure force of side AB is Pa, and the water pressure force of side CD is Pb. The water pressure force of side BC is dU, the normal force is dN, and the tangential force is dT. The angle of the failure plane with the horizontal ground surface is α, the angle of the phreatic surface with the horizontal ground surface is β, and the length of failure plane BC is dL.Figure 3
Forces acting on the slices subjected to seepage loads when the slice forces are neglected.The water pressure problem has been solved by Zheng et al. [10], and the water pressure force of BC is given by(1)dU=γwha+hbdL2cos2β.Equation (1) can be obtained by mathematical transformation as follows:(2)dU=γwha+hbdLcosα2cosαcos2β=γwdAwcos2βcosα=γwdAw1−sin2βcosα=γwdAw1−i2cosα.Therefore, the total water pressure force induced by the water seepage,U1, which is the result from the pore water as shown in equation (2), acting along the failure plane BC, is given by(3)U1=γwAw1−i2cosα,where the area of immersed soil is Aw, which is the area surrounded by the phreatic surface and failure plane. The hydraulic gradient is sinβ, and it can be simplified as average hydraulic gradient i.


## 3.2. Pore Water Pressure Induced by Earthquake

Under earthquake condition, an excess pore water pressure may be induced inside the slope. For a plane state, the excess pore water pressure △u can be obtained by a pseudostatic method of analysis and is given by Wang et al. [8] as follows:(4)Δu=γsatβkh+kv1+μ3+2αkh+kv2μ2−μ+1−3khkv,where μ is Poisson’s ratio of the soil, α and β are the coefficients of excess pore pressure, β = 1.0 for saturated soil, and kh and kv are the seismic acceleration coefficients in the horizontal and vertical directions.The total excess pore water pressure induced by the earthquake,U2, which is the result of the excess pore water, △u, shown in (4), acting alone the failure plane AF, is as follows:(5)U2=∫0hwsinαγsatβkh+kv1+μ3+2αkh+kv2μ2−μ+1−3khkvdz=γsathwβkh+kv1+μ3sinα+γsathw2αkh+kv2μ2−μ+1−3khkvsinα.


## 3.3. Force Analysis of Slip Body

The slip body is taken as the research object, as shown in Figure4. W is the total weight of the slip body, Eh is the horizontal seismic inertia force, Ev is the vertical seismic inertia force, U1 is the total water pressure caused by seepage, and U2 is the total excess pore water pressure caused by earthquake. N and S are the normal and tangential forces on the failure plane, respectively. The slip body is divided into three parts, containing AFA, AGFA, and GBCFG. The slip body is considered to be rigid in calculation, and the total force can be equivalent to the sum of three parts.Figure 4
Forces acting on the slip body.


### 3.3.1. Calculation of Slip Body Weight

In the process of determining the phreatic surface, the physical approximations are the same as Dupuit assumed conditions, which are obtained by solving the Boussinesq equation. Based on the research results of seepage analysis, a parabola model can be used to calculate the phreatic surface induced by the groundwater seepage. This model has been applied to the dam, given by Zheng et al. [10], as shown in Figure 4. The phreatic surface equation is as follows:(6)z2=ax+b,where a and b are undetermined coefficients. The hydraulic radius of phreatic surface is R, and it is AE length as shown in Figure 2. It can be seen from Figure 2 that the coordinates of point A and point M are A (0, 0) and M (h, R), respectively, and both points are on the parabola line. So, the parameters a and b can be solved by taking point A and point B into (6) as follows:(7)b=0,a=h2R.At the same time, the water level at the intersection of the phreatic surface and slip surface is expressed as(8)hw=h2Rtanα.The hydraulic radiusR can be determined by the Gehart formula as follows:(9)R=10ΔHK,where △H is the water drawdown. K is the permeability coefficient. Then, the equation of phreatic surface is given by(10)z2=h2xR.It can be seen from Figure4 that the mathematical equation of failure plane AC can be expressed as(11)z=xtanα.Based on the above conclusions, the total weight of the slip body can be written as the sum of weight of geometryAFA, AGFA, and GBCFG.The weight of wedgeAFA is given by(12)WAFA=∫0h2/Rtanαγsatztanα−z2Rh2dz=γsath46R2tan3α.The weight of wedgeAGFA is given by(13)WAGFA=∫0h2/Rtanαγz2Rh2−ztanθdz=γh42tanθ−3tanα6R2tan3αtanθ.The weight of wedgeGBCFG is given by(14)WGBCFG=∫h2/RtanαHγcotα−cotθzdz=12γcotα−cotθH2−h4R2tan2α.The total weight of the slip body is given by(15)W=WAFA+WAGFA+WGBCFG.


### 3.3.2. Calculation of Seismic Inertia Force of Slip Body

The seismic inertia force can be calculated by the pseudodynamic method, and the rationality of this method is proved by Bash and Babu [5]. It is assumed that the shear modulus is constant with depth through the soil and that only the phase, not the magnitude of acceleration, is varying. The major advantage of using sinusoidal acceleration in the pseudodynamic method is that it possesses most of the representative characters of the earthquake. Since the pseudodynamic method considers the effects of vibration characteristics, the pseudodynamic method will degenerate into the pseudostatic method when seismic wave velocity tends to infinity.The period of lateral shaking can be defined asT = 2π/w, where w is the angular frequency of the base shaking. In the present study, it is assumed that the slope is subjected to both horizontal and vertical sinusoidal acceleration with amplitude of acceleration khg and kvg, respectively, where g is the acceleration due to gravity. The variation of horizontal and vertical sinusoidal acceleration at any depth z below the top of the slope and time (t) with soil amplification factor (f) can be expressed as follows:(16)ahz,t=1+H−zHf−1khgsinwt−H−zvs,avz,t=1+H−zHf−1kvgsinwt−H−zvp,where vs is the shear wave velocity propagating through the soil and vp is the primary wave velocity propagating through the soil. ah is the horizontal seismic acceleration, and av is the vertical seismic acceleration. G is the shear modulus of soil, and g is the gravity acceleration. Based on the soil dynamics theory [3], the shear wave velocity and the primary wave velocity can be expressed as follows:(17)vs=Gρvp=G2−μρ1−2μ,where ρ is the density of soil mass and μ is Poisson’s ratio.The inertia force of different wedge can be calculated respectively. Because the value ofz in the pseudodynamic method is based on the free surface at the top of the slope, the value of z needs to be converted in Figure 4.In this study, the Cartesian coordinate system is established with the slope toe as the origin, as shown in Figure4. In order to simplify the calculation, the soil amplification factor f = 1 is assumed. Thus, the inertial forces of the slip body are calculated as follows:(1)
Horizontal and vertical inertial force acting on the wedgeAFA.The horizontal inertia force acting on the wedgeAFA can be written as follows:(18)Eh−AFA=∫0h2cotα/Rγsatgzcotα−z2Rh2khg⋅sinwt−zvsdz=m1γsatkhcoswt−h2cotαRvs+m2γsatkhsinwt−h2cotαRvs+m3γsatkhcoswt+m4γsatkhsinwt,where the coefficients can be expressed as follows:(19)m1=vs3w32Rh2,m2=−vs2w2cotα,m3=−vs3w32Rh2,m4=−vs2w2cotα.Similarly, the vertical inertia force acting on the wedgeAFA can be written as follows:(20)Ev−AFA=∫0h2cotα/Rγsatgzcotα−z2Rh2kvg⋅sinwt−zvpdz=m5γsatkvcoswt−h2cotαRvp+m6γsatkvsinwt−h2cotαRvp+m7γsatkvcoswt+m8γsatkvsinwt,where the coefficients can be expressed as follows:(21)m5=vp3w32Rh2,m6=−vp2w2cotα,m7=−vp3w32Rh2,m8=−vp2w2cotα.(2)
Horizontal and vertical inertial force acting on the wedgeAGFA.The horizontal inertia force acting on the wedgeAGFA can be written as follows:(22)Eh−AGFA=∫0h2cotα/Rγgz2Rh2−zcotθkhg⋅sinwt−zvsdz=m9γkhcoswt−h2cotαRvp+m10γkhsinwt−h2cotαRvp+m11γkhcoswt+m12γkhsinwt,where the coefficients can be expressed as follows:(23)m9=−vs3w32Rh2+vswh2cotαRcotα−cotθm10=vs2w22cotα−cotθm11=vs3w32Rh2m12=vs2w2cotθ.Similarly, the vertical inertia force acting on the wedge AGFA can be written as follows:(24)Ev−AGFA=∫0h2cotα/Rγgz2Rh2−zcotθkvg⋅sinwt−zvpdz=m13γkvcoswt−h2cotαRvp+m14γkvsinwt−h2cotαRvp+m15γkvcoswt+m16γkvsinwt,where the coefficients can be expressed as follows:(25)m13=−vp3w32Rh2+vpwh2cotαRcotα−cotθ,m14=vp2w22cotα−cotθ,m15=vp3w32Rh2,m16=vp2w2cotθ.(3)
Horizontal and vertical inertial force acting on the wedgeGBCFG.The horizontal inertia force acting on the wedge GBCFG can be written as follows:(26)Eh−GBCFG=∫h2cotα/RHγgcotα−cotθzkhgsinwt−zvsdz=m17γkhcoswt−h2cotαRvs+m18γkhsinwt−h2cotαRvs+m19γkhcoswt−Hvs+m20γkhsinwt−Hvs,where the coefficients can be expressed as follows:(27)m17=vswcotθ−cotαh2cotαR,m18=vs2w2cotθ−cotα,m19=vswcotθ−cotαH,m20=vs2w2cotθ−cotα.Similarly, the vertical inertia force acting on the wedgeGBCFG can be written as follows:(28)Ev−GBCFG=∫h2cotα/RHγgcotα−cotθzkvgsinwt−zvpdz=m21γkvcoswt−h2cotαRvp+m22γkvsinwt−h2cotαRvp+m23γkvcoswt−Hvp+m24γkvsinwt−Hvp,where the coefficients can be expressed as follows:(29)m21=vpwcotθ−cotαh2cotαR,m22=vp2w2cotθ−cotα,m23=vpwcotθ−cotαH,m24=vp2w2cotθ−cotα.Therefore, the total inertia force of the slip body is given by(30)Eh=Eh−AFA+Eh−AGFA+Eh−GBCFG,Ev=Ev−AFA+Ev−AGFA+Ev−GBCFG.


## 3.3.1. Calculation of Slip Body Weight

In the process of determining the phreatic surface, the physical approximations are the same as Dupuit assumed conditions, which are obtained by solving the Boussinesq equation. Based on the research results of seepage analysis, a parabola model can be used to calculate the phreatic surface induced by the groundwater seepage. This model has been applied to the dam, given by Zheng et al. [10], as shown in Figure 4. The phreatic surface equation is as follows:(6)z2=ax+b,where a and b are undetermined coefficients. The hydraulic radius of phreatic surface is R, and it is AE length as shown in Figure 2. It can be seen from Figure 2 that the coordinates of point A and point M are A (0, 0) and M (h, R), respectively, and both points are on the parabola line. So, the parameters a and b can be solved by taking point A and point B into (6) as follows:(7)b=0,a=h2R.At the same time, the water level at the intersection of the phreatic surface and slip surface is expressed as(8)hw=h2Rtanα.The hydraulic radiusR can be determined by the Gehart formula as follows:(9)R=10ΔHK,where △H is the water drawdown. K is the permeability coefficient. Then, the equation of phreatic surface is given by(10)z2=h2xR.It can be seen from Figure4 that the mathematical equation of failure plane AC can be expressed as(11)z=xtanα.Based on the above conclusions, the total weight of the slip body can be written as the sum of weight of geometryAFA, AGFA, and GBCFG.The weight of wedgeAFA is given by(12)WAFA=∫0h2/Rtanαγsatztanα−z2Rh2dz=γsath46R2tan3α.The weight of wedgeAGFA is given by(13)WAGFA=∫0h2/Rtanαγz2Rh2−ztanθdz=γh42tanθ−3tanα6R2tan3αtanθ.The weight of wedgeGBCFG is given by(14)WGBCFG=∫h2/RtanαHγcotα−cotθzdz=12γcotα−cotθH2−h4R2tan2α.The total weight of the slip body is given by(15)W=WAFA+WAGFA+WGBCFG.


## 3.3.2. Calculation of Seismic Inertia Force of Slip Body

The seismic inertia force can be calculated by the pseudodynamic method, and the rationality of this method is proved by Bash and Babu [5]. It is assumed that the shear modulus is constant with depth through the soil and that only the phase, not the magnitude of acceleration, is varying. The major advantage of using sinusoidal acceleration in the pseudodynamic method is that it possesses most of the representative characters of the earthquake. Since the pseudodynamic method considers the effects of vibration characteristics, the pseudodynamic method will degenerate into the pseudostatic method when seismic wave velocity tends to infinity.The period of lateral shaking can be defined asT = 2π/w, where w is the angular frequency of the base shaking. In the present study, it is assumed that the slope is subjected to both horizontal and vertical sinusoidal acceleration with amplitude of acceleration khg and kvg, respectively, where g is the acceleration due to gravity. The variation of horizontal and vertical sinusoidal acceleration at any depth z below the top of the slope and time (t) with soil amplification factor (f) can be expressed as follows:(16)ahz,t=1+H−zHf−1khgsinwt−H−zvs,avz,t=1+H−zHf−1kvgsinwt−H−zvp,where vs is the shear wave velocity propagating through the soil and vp is the primary wave velocity propagating through the soil. ah is the horizontal seismic acceleration, and av is the vertical seismic acceleration. G is the shear modulus of soil, and g is the gravity acceleration. Based on the soil dynamics theory [3], the shear wave velocity and the primary wave velocity can be expressed as follows:(17)vs=Gρvp=G2−μρ1−2μ,where ρ is the density of soil mass and μ is Poisson’s ratio.The inertia force of different wedge can be calculated respectively. Because the value ofz in the pseudodynamic method is based on the free surface at the top of the slope, the value of z needs to be converted in Figure 4.In this study, the Cartesian coordinate system is established with the slope toe as the origin, as shown in Figure4. In order to simplify the calculation, the soil amplification factor f = 1 is assumed. Thus, the inertial forces of the slip body are calculated as follows:(1)
Horizontal and vertical inertial force acting on the wedgeAFA.The horizontal inertia force acting on the wedgeAFA can be written as follows:(18)Eh−AFA=∫0h2cotα/Rγsatgzcotα−z2Rh2khg⋅sinwt−zvsdz=m1γsatkhcoswt−h2cotαRvs+m2γsatkhsinwt−h2cotαRvs+m3γsatkhcoswt+m4γsatkhsinwt,where the coefficients can be expressed as follows:(19)m1=vs3w32Rh2,m2=−vs2w2cotα,m3=−vs3w32Rh2,m4=−vs2w2cotα.Similarly, the vertical inertia force acting on the wedgeAFA can be written as follows:(20)Ev−AFA=∫0h2cotα/Rγsatgzcotα−z2Rh2kvg⋅sinwt−zvpdz=m5γsatkvcoswt−h2cotαRvp+m6γsatkvsinwt−h2cotαRvp+m7γsatkvcoswt+m8γsatkvsinwt,where the coefficients can be expressed as follows:(21)m5=vp3w32Rh2,m6=−vp2w2cotα,m7=−vp3w32Rh2,m8=−vp2w2cotα.(2)
Horizontal and vertical inertial force acting on the wedgeAGFA.The horizontal inertia force acting on the wedgeAGFA can be written as follows:(22)Eh−AGFA=∫0h2cotα/Rγgz2Rh2−zcotθkhg⋅sinwt−zvsdz=m9γkhcoswt−h2cotαRvp+m10γkhsinwt−h2cotαRvp+m11γkhcoswt+m12γkhsinwt,where the coefficients can be expressed as follows:(23)m9=−vs3w32Rh2+vswh2cotαRcotα−cotθm10=vs2w22cotα−cotθm11=vs3w32Rh2m12=vs2w2cotθ.Similarly, the vertical inertia force acting on the wedge AGFA can be written as follows:(24)Ev−AGFA=∫0h2cotα/Rγgz2Rh2−zcotθkvg⋅sinwt−zvpdz=m13γkvcoswt−h2cotαRvp+m14γkvsinwt−h2cotαRvp+m15γkvcoswt+m16γkvsinwt,where the coefficients can be expressed as follows:(25)m13=−vp3w32Rh2+vpwh2cotαRcotα−cotθ,m14=vp2w22cotα−cotθ,m15=vp3w32Rh2,m16=vp2w2cotθ.(3)
Horizontal and vertical inertial force acting on the wedgeGBCFG.The horizontal inertia force acting on the wedge GBCFG can be written as follows:(26)Eh−GBCFG=∫h2cotα/RHγgcotα−cotθzkhgsinwt−zvsdz=m17γkhcoswt−h2cotαRvs+m18γkhsinwt−h2cotαRvs+m19γkhcoswt−Hvs+m20γkhsinwt−Hvs,where the coefficients can be expressed as follows:(27)m17=vswcotθ−cotαh2cotαR,m18=vs2w2cotθ−cotα,m19=vswcotθ−cotαH,m20=vs2w2cotθ−cotα.Similarly, the vertical inertia force acting on the wedgeGBCFG can be written as follows:(28)Ev−GBCFG=∫h2cotα/RHγgcotα−cotθzkvgsinwt−zvpdz=m21γkvcoswt−h2cotαRvp+m22γkvsinwt−h2cotαRvp+m23γkvcoswt−Hvp+m24γkvsinwt−Hvp,where the coefficients can be expressed as follows:(29)m21=vpwcotθ−cotαh2cotαR,m22=vp2w2cotθ−cotα,m23=vpwcotθ−cotαH,m24=vp2w2cotθ−cotα.Therefore, the total inertia force of the slip body is given by(30)Eh=Eh−AFA+Eh−AGFA+Eh−GBCFG,Ev=Ev−AFA+Ev−AGFA+Ev−GBCFG.


## 3.4. Stability Calculation of Slope

Limit equilibrium method is used to determine the safety factor against slope failure. The safety factor of the sandy slope can be expressed by the ratio of resisting force and sliding force on slip surface. The equilibrium equation of normal direction of slip surface (∑N = 0) is given by(31)N+U1+U2+Ehsinα−W−Evcosα=0.The equilibrium equation of tangential direction of slip surface (∑S = 0) is given by(32)Ehcosα+W−Evsinα−S=0.Theoretically, the cohesion of sandy soil is zero, but it is not zero in practice. In order to ensure the accuracy of the analytical calculation, the cohesion is assumed to exist in the resisting force calculation and then chosen according to its actual value. Therefore, the safety factor can be written as follows:(33)F=cl+W−Evcosα−Ehsinα−U1−U2tanφEhcosα+W−Evsinα,where c is the cohesion of soil mass and φ is the internal friction angle of soil mass.In (33), the relations of the terms b, l, α, and θ can be written as(34)b=Hcotα−cotθ,l=Hcscα.Based on the above work, when the seismic parameters, slope parameters, and hydraulic parameters are known, the safety factor is a function of the slip surface angleα and time t. The critical slip surface, whose safety factor is the minimum among all the available ones, can be determined by the objective function; therefore, the inclination angle extreme and the safety factor at each time can be obtained by the function extreme theory. The minimum safety factor of slope in the vibration time domain will be used to evaluate the stability of slope.It is generally believed that the cohesive force of sandy soil is 0 kPa. But in fact, the cohesive force of unsaturated sand and the cohesive force of silty sand are nonzero. The principle of selecting cohesive force is given by early researchers’ works [4]. If the cohesive force of sandy soil measured in the experiment is less than 6 kPa, it can be equivalent to pure sandy soil, and the calculated cohesive force is set to zero. When the cohesion is greater than 6 kPa, the cohesion force is generally reduced by half in calculation. The purpose of this processing in this study is to make better use of the classification experimental data of sandy soil for calculation.


## 4. Result and Discussion

### 4.1. Verification of the Present Formulation

To verify the accuracy of the present solution, the safety factors of the two cases obtained in the present study were compared with the values reported in studies by Song et al. [29] and Guan [31] based on different methods.Case 1: The slope has weightH = 10 m and gradient tanθ = 1 : 1.5. The characteristics of the soil material are as follows: c = 3 kPa, φ = 40°, γ = 20.5 kN/m3, and γsat = 22 kN/m3. The shear modulus of soil is G = 19.5 MPa, and Poisson’s ratio of soil is μ = 0.3.Under earthquake condition, the horizontal seismic acceleration is 0.1 g, and no vertical acceleration is assumed to exist in this example. The height of groundwater level from the bottom of slope is 2 m, 6 m, and 10 m, respectively. The shear wave velocity and primary wave velocity can be expressed as equation (3), so in this present analysis, vs = 97.53 m/s and vp = 201.06 m/s. The angular frequency w can be calculated by the relation provided by Kramer, as T = 2π/w = 4H/vs; thus, w = 15.32 rad/s and T = 0.41 s. The horizontal seismic acceleration coefficient is kh = 0.1, and the pore water pressure coefficient is α = 0.75 and β = 1.0, respectively. In the computation, the cohesion of sandy soil is set to zero. Table 1 is the result of the safety factor computed by different methods.Case 2: The slope has weightH = 24 m and slope angle θ = 30°. The characteristics of the soil material are as follows: c = 7 kPa, φ = 36°, γ = 20.5 kN/m3, and kh = 0.21. The shear modulus of soil is G = 23.8 MPa, Poisson’s ratio of soil is μ = 0.3, and the period of sinusoidal acceleration is T = 0.6 s. The height of groundwater level is 3 m, 6 m, and 9 m, respectively. Table 2 is the result of the safety factor computed by different methods.Results presented in Tables1 and 2 indicate that the present solution is in reasonably good agreement with the case values.Table 1
Safety factor by different methods.
MethodFh = 2 mh = 6 mh = 10 mDynamic time history method [29]1.340.870.25Standard method1.361.010.82Analytical method1.340.970.74Table 2
Safety factor by different methods.
MethodFh = 3 mh = 6 mh = 9 mDynamic time history method1.221.091.02Standard method [31]1.271.221.17Analytical method1.201.040.96


### 4.2. Different Methods versus Analytical Method

Again, to aid the understanding of the accuracy of the present solution for groundwater and seismic cases, the safety factorF is compared with the values reported from the different methods; for h = 2 m, 6 m, and 10 m, kh = 0.1 and kv = 0 and for the h = 3 m, 6 m, and 9 m, kh = 0.21 and kv = 0. The values of safety factor predicted by the analytical method are lower than those of the standard methods, and the difference increases for the higher values of h. The reasons resulting in the differences come from the difference in computing the excess pore water pressure induced by the earthquake shaking. In the present solutions, the excess pore water pressure is assumed to be induced by both the horizontal and vertical earthquake forces, but in the standard method, it is assumed to be only related to the dynamic pore water pressure model. The results of case 1 computed by the dynamic time history method are lower than the values computed by the analytical method, which is restricted by strict geometric boundaries due to the consideration of the real seismic wave as input condition. It can be concluded that the dynamic characteristics of real seismic wave are more significant than those of the pseudodynamic method, such as the effect of time and phase difference due to shear wave and primary wave velocities.


### 4.3. The Dynamic Safety Factor under Water Seepage and Earthquake

In Figure5, the curve represents the case 2 results of the analytical method, and the solid line stands for the results of the standard method. The dotted line stands for the critical value of safety factor (F = 1). It can be seen that the values of safety factor from the analytical method change with time. It shows obvious periodicity, and the period of the safety factor F is similar to that of the seismic wave from the pseudodynamic method. In this study, the minimum safety factor is considered as the dynamic safety factor. When the dynamic safety factor F is larger than the critical value 1, the slope is stable. On the contrary, when the dynamic safety factor F is less than the critical value 1, the slope starts to slide. It is found from Figure 5 that the dynamic safety factor is less than the standard method.Figure 5
Safety factor comparison with different methods. (a)h = 3 m, kh = 0.21, (b) h = 6 m, kh = 0.21, and (c) h = 9 m, kh = 0.21.
(a)(b)(c)According to the different water level slope, the seismic loads acting on the sliding body vary with time, which leads to a change of location of the slip surface with time. Therefore, in order to find out the relationship between time and the safety factor, it is necessary to determine the slip surface at different times. The slip surface location can be expressed by the angle of the failure plane with the horizontal ground surfaceα, in accordance with (33), and the slip surface can be obtained. Figure 6 shows the relationship between safety factor and angle α. The safety factor curve is larger than the critical value 1, and the slope is safety. The slip surface corresponds to the minimum safety factor; it can be concluded that inclination angle α is also different under different water levels. Figure 6 also shows that the higher the water level is, the more serious the loss of the safety factor is and the easier the slope is to lose stability.Figure 6
Relationship between safety factorF and inclination angle α. (a) Results of case 1. (b) Results of case 2.
(a)(b)


### 4.4. The Effects of Parameters on the Stability of Slope under Water Seepage and Earthquake

#### 4.4.1. The Effects of Internal Friction Angleφ on the Safety Factor F

The focus of this study is on sandy slope because the cohesion of sandy soil is very small; only the influence of internal friction angle on safety factor is considered. The range and values of parameters involved in case 2 arec = 0 kPa, kh = 0.21, and φ = 5°∼35°. As the critical slip surface differs from parameter to parameter, it is necessary to find out the critical slip surface for each of different sets of parameters in the computations. It can be observed from Figure 7 that the internal friction angle φ plays an important role in the stability of slope. The safety factor linearly increases as the internal friction angle φ increases. Additionally, the safety factor obviously increases with decreases in water level h. And the higher the water level h is, the faster the safety factor F decreases under the same internal friction angle φ. For example, with a constant value of φ = 35°, the safety factors for water level conditions with h = 6 m and h = 9 m are about 15.03% and 35.07%, respectively, smaller than the water level h = 3 m.Figure 7
Effect of water levelh and internal friction angle φ on safety factor (kh = 0.21).


#### 4.4.2. The Effects of Water Levelh and Slope Angle θ on the Safety Factor F

The range and values of parameters involved in case 1 are the water levelh = 0.5∼2.0 m and the slope angle θ = 15°∼35° and in case 2 are the water level h = 1.0∼4.0 m and the slope angle θ = 15°～30°. It is clearly found from Figure 8 that the safety factor decreases nonlinearly with increasing water level and slope angle. It indicates that the higher the water level is, the greater the slope angle is and the more likely the slope is to be unstable. For instance, in case 1, for h = 2 m, when slope angle increases from 15° to 35°, the value of the safety factor decreases by 60% (from 4.06 to 1.61).Figure 8
Effect of water levelh and slope angle θ on safety factor. (a) Result of case 1 (kh = 0.1). (b) Result of case 2 (kh = 0.21).
(a)(b)


#### 4.4.3. The Effects of Seismic Acceleration Coefficient on the Safety FactorF

During an earthquake, the inertia forces’ direction acting on the slip body varies with time. Here, the range and values of parameters involved in case 1 and case 2 arekh = −0.4∼0.4 and kh = −0.4∼0.4. It is seen from Figure 9 that the effects of seismic acceleration coefficients kh and kv on the slope stability are similar, and the safety factor changes nonlinearly with an increase of seismic acceleration coefficients. When the horizontal seismic acceleration coefficient kh increases, it corresponds to the horizontal inertia force whose direction points to the outside of the slope, and then the safety factor F decreases. And similarly, when the vertical seismic acceleration coefficient kv decreases, it corresponds to the vertical inertia force whose direction points to the downside of the slope, and then the safety factor F decreases. Therefore, the seismic inertia forces pointing to the outside and downside of the slope will attenuate the slope stability. The influence of the horizontal seismic acceleration coefficient on the slope stability is significantly greater than that of the vertical seismic acceleration coefficient, which is also the reason why the horizontal seismic action is mainly considered in stability analysis.Figure 9
Effect of seismic acceleration coefficient on safety factor. (a) Horizontal seismic acceleration coefficient with safety factor. (b) Horizontal seismic acceleration coefficient with safety factor.
(a)(b)


#### 4.4.4. Effects of Seismic Acceleration Coefficient on Safety FactorF and Inclination Angle α

The effects of the horizontal seismic acceleration coefficient on the safety factor and inclination angle are depicted in Figure10. The range and values of parameters involved in case 1 and case 2 are kh = −0.4∼0.4 and kh = −0.3∼0.3. It is observed from Figure 10 that the values of the safety factor and inclination angle decrease with an increase of the horizontal seismic acceleration coefficient, and the safety factor is highly sensitive to the change in kh values. Hence, the horizontal seismic inertia force plays a critical role in the instability of the slope subjected to seismic loading. As for an illustration in case 1, when kh increases from −0.4 to 0.4, the value of the safety factor decreases by 85% (from 5.81 to 0.85). Thus, it can be concluded that the slip body volume shows a significant increase with an increase in the value of the horizontal seismic acceleration coefficient whose direction points to the outside of slope.Figure 10
Effect of horizontal seismic acceleration coefficient on safety factorF and inclination angle α. (a) Result of case 1 (h = 2 m). (b) Result of case 2 (h = 3 m).
(a)(b)Figure11 illustrates the effects of the vertical seismic acceleration coefficient on the safety factor and inclination angle. The range and values of parameters involved in case 1 and case 2 are kv = −0.4∼0.4 and kv = −0.3∼0.3. It is found from Figure 11 that the values of the safety factor and inclination angle increase with an increase of the vertical seismic acceleration coefficient. It can also be concluded that the slip body volume shows a significant increase with an increase in value of the vertical seismic acceleration coefficient whose direction points to the downside of slope. Compared with Figure 10, it can be seen that the effect of the vertical seismic acceleration coefficient on the safety factor is lower than that of the horizontal seismic acceleration coefficient. As for an illustration in case 1, when kv increases from −0.4 to 0.4, the value of safety factor increases only by 20% (from 1.67 to 2.01).Figure 11
Effect of vertical seismic acceleration coefficient on safety factorF and inclination angle α. (a) Result of case 1 (h = 2 m). (b) Result of case 2 (h = 3 m).
(a)(b)The size and position of slip body under the different horizontal seismic acceleration coefficientskh can be seen in Figure 12. The slip volume increases with an increase of the horizontal seismic acceleration coefficient, and it is easily observed that the horizontal seismic force obviously affects the state of slope.Figure 12
Slip body state under earthquake action. (a) Result of case 1. (b) Result of case 2.
(a)(b)


## 4.1. Verification of the Present Formulation

To verify the accuracy of the present solution, the safety factors of the two cases obtained in the present study were compared with the values reported in studies by Song et al. [29] and Guan [31] based on different methods.Case 1: The slope has weightH = 10 m and gradient tanθ = 1 : 1.5. The characteristics of the soil material are as follows: c = 3 kPa, φ = 40°, γ = 20.5 kN/m3, and γsat = 22 kN/m3. The shear modulus of soil is G = 19.5 MPa, and Poisson’s ratio of soil is μ = 0.3.Under earthquake condition, the horizontal seismic acceleration is 0.1 g, and no vertical acceleration is assumed to exist in this example. The height of groundwater level from the bottom of slope is 2 m, 6 m, and 10 m, respectively. The shear wave velocity and primary wave velocity can be expressed as equation (3), so in this present analysis, vs = 97.53 m/s and vp = 201.06 m/s. The angular frequency w can be calculated by the relation provided by Kramer, as T = 2π/w = 4H/vs; thus, w = 15.32 rad/s and T = 0.41 s. The horizontal seismic acceleration coefficient is kh = 0.1, and the pore water pressure coefficient is α = 0.75 and β = 1.0, respectively. In the computation, the cohesion of sandy soil is set to zero. Table 1 is the result of the safety factor computed by different methods.Case 2: The slope has weightH = 24 m and slope angle θ = 30°. The characteristics of the soil material are as follows: c = 7 kPa, φ = 36°, γ = 20.5 kN/m3, and kh = 0.21. The shear modulus of soil is G = 23.8 MPa, Poisson’s ratio of soil is μ = 0.3, and the period of sinusoidal acceleration is T = 0.6 s. The height of groundwater level is 3 m, 6 m, and 9 m, respectively. Table 2 is the result of the safety factor computed by different methods.Results presented in Tables1 and 2 indicate that the present solution is in reasonably good agreement with the case values.Table 1
Safety factor by different methods.
MethodFh = 2 mh = 6 mh = 10 mDynamic time history method [29]1.340.870.25Standard method1.361.010.82Analytical method1.340.970.74Table 2
Safety factor by different methods.
MethodFh = 3 mh = 6 mh = 9 mDynamic time history method1.221.091.02Standard method [31]1.271.221.17Analytical method1.201.040.96


## 4.2. Different Methods versus Analytical Method

Again, to aid the understanding of the accuracy of the present solution for groundwater and seismic cases, the safety factorF is compared with the values reported from the different methods; for h = 2 m, 6 m, and 10 m, kh = 0.1 and kv = 0 and for the h = 3 m, 6 m, and 9 m, kh = 0.21 and kv = 0. The values of safety factor predicted by the analytical method are lower than those of the standard methods, and the difference increases for the higher values of h. The reasons resulting in the differences come from the difference in computing the excess pore water pressure induced by the earthquake shaking. In the present solutions, the excess pore water pressure is assumed to be induced by both the horizontal and vertical earthquake forces, but in the standard method, it is assumed to be only related to the dynamic pore water pressure model. The results of case 1 computed by the dynamic time history method are lower than the values computed by the analytical method, which is restricted by strict geometric boundaries due to the consideration of the real seismic wave as input condition. It can be concluded that the dynamic characteristics of real seismic wave are more significant than those of the pseudodynamic method, such as the effect of time and phase difference due to shear wave and primary wave velocities.


## 4.3. The Dynamic Safety Factor under Water Seepage and Earthquake

In Figure5, the curve represents the case 2 results of the analytical method, and the solid line stands for the results of the standard method. The dotted line stands for the critical value of safety factor (F = 1). It can be seen that the values of safety factor from the analytical method change with time. It shows obvious periodicity, and the period of the safety factor F is similar to that of the seismic wave from the pseudodynamic method. In this study, the minimum safety factor is considered as the dynamic safety factor. When the dynamic safety factor F is larger than the critical value 1, the slope is stable. On the contrary, when the dynamic safety factor F is less than the critical value 1, the slope starts to slide. It is found from Figure 5 that the dynamic safety factor is less than the standard method.Figure 5
Safety factor comparison with different methods. (a)h = 3 m, kh = 0.21, (b) h = 6 m, kh = 0.21, and (c) h = 9 m, kh = 0.21.
(a)(b)(c)According to the different water level slope, the seismic loads acting on the sliding body vary with time, which leads to a change of location of the slip surface with time. Therefore, in order to find out the relationship between time and the safety factor, it is necessary to determine the slip surface at different times. The slip surface location can be expressed by the angle of the failure plane with the horizontal ground surfaceα, in accordance with (33), and the slip surface can be obtained. Figure 6 shows the relationship between safety factor and angle α. The safety factor curve is larger than the critical value 1, and the slope is safety. The slip surface corresponds to the minimum safety factor; it can be concluded that inclination angle α is also different under different water levels. Figure 6 also shows that the higher the water level is, the more serious the loss of the safety factor is and the easier the slope is to lose stability.Figure 6
Relationship between safety factorF and inclination angle α. (a) Results of case 1. (b) Results of case 2.
(a)(b)


## 4.4. The Effects of Parameters on the Stability of Slope under Water Seepage and Earthquake

### 4.4.1. The Effects of Internal Friction Angleφ on the Safety Factor F

The focus of this study is on sandy slope because the cohesion of sandy soil is very small; only the influence of internal friction angle on safety factor is considered. The range and values of parameters involved in case 2 arec = 0 kPa, kh = 0.21, and φ = 5°∼35°. As the critical slip surface differs from parameter to parameter, it is necessary to find out the critical slip surface for each of different sets of parameters in the computations. It can be observed from Figure 7 that the internal friction angle φ plays an important role in the stability of slope. The safety factor linearly increases as the internal friction angle φ increases. Additionally, the safety factor obviously increases with decreases in water level h. And the higher the water level h is, the faster the safety factor F decreases under the same internal friction angle φ. For example, with a constant value of φ = 35°, the safety factors for water level conditions with h = 6 m and h = 9 m are about 15.03% and 35.07%, respectively, smaller than the water level h = 3 m.Figure 7
Effect of water levelh and internal friction angle φ on safety factor (kh = 0.21).


### 4.4.2. The Effects of Water Levelh and Slope Angle θ on the Safety Factor F

The range and values of parameters involved in case 1 are the water levelh = 0.5∼2.0 m and the slope angle θ = 15°∼35° and in case 2 are the water level h = 1.0∼4.0 m and the slope angle θ = 15°～30°. It is clearly found from Figure 8 that the safety factor decreases nonlinearly with increasing water level and slope angle. It indicates that the higher the water level is, the greater the slope angle is and the more likely the slope is to be unstable. For instance, in case 1, for h = 2 m, when slope angle increases from 15° to 35°, the value of the safety factor decreases by 60% (from 4.06 to 1.61).Figure 8
Effect of water levelh and slope angle θ on safety factor. (a) Result of case 1 (kh = 0.1). (b) Result of case 2 (kh = 0.21).
(a)(b)


### 4.4.3. The Effects of Seismic Acceleration Coefficient on the Safety FactorF

During an earthquake, the inertia forces’ direction acting on the slip body varies with time. Here, the range and values of parameters involved in case 1 and case 2 arekh = −0.4∼0.4 and kh = −0.4∼0.4. It is seen from Figure 9 that the effects of seismic acceleration coefficients kh and kv on the slope stability are similar, and the safety factor changes nonlinearly with an increase of seismic acceleration coefficients. When the horizontal seismic acceleration coefficient kh increases, it corresponds to the horizontal inertia force whose direction points to the outside of the slope, and then the safety factor F decreases. And similarly, when the vertical seismic acceleration coefficient kv decreases, it corresponds to the vertical inertia force whose direction points to the downside of the slope, and then the safety factor F decreases. Therefore, the seismic inertia forces pointing to the outside and downside of the slope will attenuate the slope stability. The influence of the horizontal seismic acceleration coefficient on the slope stability is significantly greater than that of the vertical seismic acceleration coefficient, which is also the reason why the horizontal seismic action is mainly considered in stability analysis.Figure 9
Effect of seismic acceleration coefficient on safety factor. (a) Horizontal seismic acceleration coefficient with safety factor. (b) Horizontal seismic acceleration coefficient with safety factor.
(a)(b)


### 4.4.4. Effects of Seismic Acceleration Coefficient on Safety FactorF and Inclination Angle α

The effects of the horizontal seismic acceleration coefficient on the safety factor and inclination angle are depicted in Figure10. The range and values of parameters involved in case 1 and case 2 are kh = −0.4∼0.4 and kh = −0.3∼0.3. It is observed from Figure 10 that the values of the safety factor and inclination angle decrease with an increase of the horizontal seismic acceleration coefficient, and the safety factor is highly sensitive to the change in kh values. Hence, the horizontal seismic inertia force plays a critical role in the instability of the slope subjected to seismic loading. As for an illustration in case 1, when kh increases from −0.4 to 0.4, the value of the safety factor decreases by 85% (from 5.81 to 0.85). Thus, it can be concluded that the slip body volume shows a significant increase with an increase in the value of the horizontal seismic acceleration coefficient whose direction points to the outside of slope.Figure 10
Effect of horizontal seismic acceleration coefficient on safety factorF and inclination angle α. (a) Result of case 1 (h = 2 m). (b) Result of case 2 (h = 3 m).
(a)(b)Figure11 illustrates the effects of the vertical seismic acceleration coefficient on the safety factor and inclination angle. The range and values of parameters involved in case 1 and case 2 are kv = −0.4∼0.4 and kv = −0.3∼0.3. It is found from Figure 11 that the values of the safety factor and inclination angle increase with an increase of the vertical seismic acceleration coefficient. It can also be concluded that the slip body volume shows a significant increase with an increase in value of the vertical seismic acceleration coefficient whose direction points to the downside of slope. Compared with Figure 10, it can be seen that the effect of the vertical seismic acceleration coefficient on the safety factor is lower than that of the horizontal seismic acceleration coefficient. As for an illustration in case 1, when kv increases from −0.4 to 0.4, the value of safety factor increases only by 20% (from 1.67 to 2.01).Figure 11
Effect of vertical seismic acceleration coefficient on safety factorF and inclination angle α. (a) Result of case 1 (h = 2 m). (b) Result of case 2 (h = 3 m).
(a)(b)The size and position of slip body under the different horizontal seismic acceleration coefficientskh can be seen in Figure 12. The slip volume increases with an increase of the horizontal seismic acceleration coefficient, and it is easily observed that the horizontal seismic force obviously affects the state of slope.Figure 12
Slip body state under earthquake action. (a) Result of case 1. (b) Result of case 2.
(a)(b)


## 4.4.1. The Effects of Internal Friction Angleφ on the Safety Factor F

The focus of this study is on sandy slope because the cohesion of sandy soil is very small; only the influence of internal friction angle on safety factor is considered. The range and values of parameters involved in case 2 arec = 0 kPa, kh = 0.21, and φ = 5°∼35°. As the critical slip surface differs from parameter to parameter, it is necessary to find out the critical slip surface for each of different sets of parameters in the computations. It can be observed from Figure 7 that the internal friction angle φ plays an important role in the stability of slope. The safety factor linearly increases as the internal friction angle φ increases. Additionally, the safety factor obviously increases with decreases in water level h. And the higher the water level h is, the faster the safety factor F decreases under the same internal friction angle φ. For example, with a constant value of φ = 35°, the safety factors for water level conditions with h = 6 m and h = 9 m are about 15.03% and 35.07%, respectively, smaller than the water level h = 3 m.Figure 7
Effect of water levelh and internal friction angle φ on safety factor (kh = 0.21).


## 4.4.2. The Effects of Water Levelh and Slope Angle θ on the Safety Factor F

The range and values of parameters involved in case 1 are the water levelh = 0.5∼2.0 m and the slope angle θ = 15°∼35° and in case 2 are the water level h = 1.0∼4.0 m and the slope angle θ = 15°～30°. It is clearly found from Figure 8 that the safety factor decreases nonlinearly with increasing water level and slope angle. It indicates that the higher the water level is, the greater the slope angle is and the more likely the slope is to be unstable. For instance, in case 1, for h = 2 m, when slope angle increases from 15° to 35°, the value of the safety factor decreases by 60% (from 4.06 to 1.61).Figure 8
Effect of water levelh and slope angle θ on safety factor. (a) Result of case 1 (kh = 0.1). (b) Result of case 2 (kh = 0.21).
(a)(b)


## 4.4.3. The Effects of Seismic Acceleration Coefficient on the Safety FactorF

During an earthquake, the inertia forces’ direction acting on the slip body varies with time. Here, the range and values of parameters involved in case 1 and case 2 arekh = −0.4∼0.4 and kh = −0.4∼0.4. It is seen from Figure 9 that the effects of seismic acceleration coefficients kh and kv on the slope stability are similar, and the safety factor changes nonlinearly with an increase of seismic acceleration coefficients. When the horizontal seismic acceleration coefficient kh increases, it corresponds to the horizontal inertia force whose direction points to the outside of the slope, and then the safety factor F decreases. And similarly, when the vertical seismic acceleration coefficient kv decreases, it corresponds to the vertical inertia force whose direction points to the downside of the slope, and then the safety factor F decreases. Therefore, the seismic inertia forces pointing to the outside and downside of the slope will attenuate the slope stability. The influence of the horizontal seismic acceleration coefficient on the slope stability is significantly greater than that of the vertical seismic acceleration coefficient, which is also the reason why the horizontal seismic action is mainly considered in stability analysis.Figure 9
Effect of seismic acceleration coefficient on safety factor. (a) Horizontal seismic acceleration coefficient with safety factor. (b) Horizontal seismic acceleration coefficient with safety factor.
(a)(b)


## 4.4.4. Effects of Seismic Acceleration Coefficient on Safety FactorF and Inclination Angle α

The effects of the horizontal seismic acceleration coefficient on the safety factor and inclination angle are depicted in Figure10. The range and values of parameters involved in case 1 and case 2 are kh = −0.4∼0.4 and kh = −0.3∼0.3. It is observed from Figure 10 that the values of the safety factor and inclination angle decrease with an increase of the horizontal seismic acceleration coefficient, and the safety factor is highly sensitive to the change in kh values. Hence, the horizontal seismic inertia force plays a critical role in the instability of the slope subjected to seismic loading. As for an illustration in case 1, when kh increases from −0.4 to 0.4, the value of the safety factor decreases by 85% (from 5.81 to 0.85). Thus, it can be concluded that the slip body volume shows a significant increase with an increase in the value of the horizontal seismic acceleration coefficient whose direction points to the outside of slope.Figure 10
Effect of horizontal seismic acceleration coefficient on safety factorF and inclination angle α. (a) Result of case 1 (h = 2 m). (b) Result of case 2 (h = 3 m).
(a)(b)Figure11 illustrates the effects of the vertical seismic acceleration coefficient on the safety factor and inclination angle. The range and values of parameters involved in case 1 and case 2 are kv = −0.4∼0.4 and kv = −0.3∼0.3. It is found from Figure 11 that the values of the safety factor and inclination angle increase with an increase of the vertical seismic acceleration coefficient. It can also be concluded that the slip body volume shows a significant increase with an increase in value of the vertical seismic acceleration coefficient whose direction points to the downside of slope. Compared with Figure 10, it can be seen that the effect of the vertical seismic acceleration coefficient on the safety factor is lower than that of the horizontal seismic acceleration coefficient. As for an illustration in case 1, when kv increases from −0.4 to 0.4, the value of safety factor increases only by 20% (from 1.67 to 2.01).Figure 11
Effect of vertical seismic acceleration coefficient on safety factorF and inclination angle α. (a) Result of case 1 (h = 2 m). (b) Result of case 2 (h = 3 m).
(a)(b)The size and position of slip body under the different horizontal seismic acceleration coefficientskh can be seen in Figure 12. The slip volume increases with an increase of the horizontal seismic acceleration coefficient, and it is easily observed that the horizontal seismic force obviously affects the state of slope.Figure 12
Slip body state under earthquake action. (a) Result of case 1. (b) Result of case 2.
(a)(b)


## 5. Conclusions

In this paper, a pseudodynamic model is proposed to determine the dynamic safety factor of sandy slope subjected to groundwater seepage and seismic loads. The effects of variation of different parameters on the stability of slope, namely, water level, seismic acceleration coefficients, slope angle, and soil internal friction angle, are investigated here. After comparisons between the results of the present and previous methods, it is found that the results of the present method are in good agreement with those from the previous ones. The major findings of the present investigation are summarized as follows:(1)
The safety factor shows periodicity. The period of the safety factor is similar to that of seismic wave, and the minimum value of the safety factor is treated as the dynamic safety factor.(2)
The groundwater level plays a significant role in the stability of slope subjected to the groundwater seepage. The value of the safety factor decreases with an increase in water level.(3)
The dynamic safety factor increases linearly with an increase in the soil internal friction angle and decreases nonlinearly with the slope angle.(4)
The dynamic safety factor is sensitive to seismic acceleration coefficientskh and kv, and the effects of those coefficients on the safety factor are different at different times, but the influence of the horizontal seismic acceleration coefficient on the slope stability is significantly greater than that of vertical seismic acceleration coefficient.The values obtained from the present analysis are reasonably comparable; thus, they can be used for the analysis of the sandy slope under complex loading conditions.



---

*Source: 5286878-2022-05-25.xml*