Upper Bound Analysis of the Stability of 3D Slopes in the Saturated Soft Clay Subjected to Seismic Effect

In order to study the three-dimensional stability problem of the saturated soft clay slope under earthquake loads, based on the three-dimensional rotation failure model, the seismic force was introduced into the calculation by the quasi-static method. The work rate of external loads and the internal energy dissipation rate of the saturated soft clay slope were calculated using the upper bound method of limit analysis, and the analytical solution of stability coefficient of saturated soft clay slopes was derived based on the fictitious power principle. By virtue of the exhaust algorithm, the optimal solution of stability coefficient of saturated soft clay slopes was obtained. The influence of the slope angle and the horizontal and vertical seismic forces on the stability coefficient of saturated soft clay slope was analyzed. The results show that the slope angle has a great influence on the stability coefficient, and the relative difference is up to 35.7%. Therefore, the stability coefficient of saturated soft clay slopes can be effectively increased by a proper slope setting. The horizontal and vertical seismic forces also have a significant influence on the stability of saturated soft clay slopes. The relative differences of the stability coefficient under horizontal and vertical seismic forces are as high as 41 and 14.7%, respectively. If they are ignored, the stability coefficient of saturated soft clay slopes will be seriously overestimated. It is suggested that the effects of horizontal and vertical seismic forces must be considered simultaneously in the seismic design of saturated soft clay slopes.

1 Introduction

Safety of the slope has always been a concern in engineering. Soil slopes are usually made of clay, sand, silt, etc. These soil slopes tend to fail under the action of earthquake and other external loads, thus causing huge loss to related projects (Leong and Rahardjo 2012; Gofar and Rahardjo 2017; VandenBerge and McGuire 2019).

Alejano et al. (2011), Selcuk et al. (2015) concluded that considering the safety of slope under the action of earthquake, the limit analysis method and the limit equilibrium method are mainly adopted. Michalowski and Park (2020) believe that the principle of limit equilibrium method is simple, but it only considers the yield condition and stress equilibrium condition of soil, and ignores the constitutive relation between soils. The limit analysis method overcomes this defect and establishes the constitutive relation of soil through the flow law (Yang and Liu 2018). Compared with the limit equilibrium method, it is stricter in theory and more accurate in calculation, so it becomes an efficacious method to ensure the slope safety (Gischig et al., 2015; Rawat and Gupta 2016). The stability of the slope is considered by using the limit analysis method, and at present, the two-dimensional plane strain problem is gradually developed into a three-dimensional problem. In view of the complex soil slope, Jongmin et al. (2002) analyzed the slope stability and found that the method of limit analysis is better than the limit equilibrium method. Yao and Yang (2017) studied the stability of the unsaturated soil slope by introducing a uniform shear strength formula to consider the influence of intermediate principal stress. For the three-dimensional stability of the slope, Huang et al. (2002) expanded the 3D stability analysis method, making it possible to use “ two-directional force and moment equilibrium” in any shape of 3D failure mechanism. This method can accurately compute the sliding direction of the failure surface and reduce the tedious work of conventional methods. Gao et al. (2015) extended the plane strain analysis of the slope strengthened by pile arrangement to the three-dimensional situation and used the limit analysis upper limit method to determine the safety coefficient. Aiming at the safety of the unsaturated three-dimensional soil slope under the condition of steady seepage, Wang et al. (2019) studied the safety of the 3D slope under the suction induction effect and permeability. Han et al. (2014) analyzed the effect of heterogeneity and anisotropy on the safety of 3D slopes.

So far, there have been many results of using the upper limit method to discuss the safety of slopes under earthquake problems. In 2009, Michalowski and Drescher (2009), Michalowski (2010) put forward a 3D rotational failure mechanism of the slope and investigated the impact of seismic forces on the slope stability and safety coefficient by the upper limit method. On the basis of the upper limit theorem of limit analysis, Nian et al. (2016) introduced the quasi-static method to research the safety of the anchored slope under the action of seismic forces. Zhang et al. (2016) studied the impact of seismic forces on the safety of the 3D slope by using the limit analysis upper bound method. Sahoo et al. (2016) studied the seismic failure mechanism of the non-reinforced soil slope and reinforced soil slope in detail using the shaking table test. To sum up, there are many research studies on the seismic problem of the slope, but few studies on the simultaneous action of horizontal and vertical seismic forces and saturated soft clay. Therefore, in this work, the effect of horizontal and vertical seismic forces on the stability of the saturated soft clay slope was studied based on the 3D rotational failure mechanism of the slope and the limit analysis upper bound method combined with the quasi-static method, so as to provide reference for the seismic design of similar projects in the future.

2 Upper Bound Theorem of Limit Analysis

The aim of the upper bound theorem of limit analysis is to solve the limit failure load; for any hypothetical failure mode satisfying the allowable velocity field of motion, the upper limit of the ultimate load is the load obtained from the external power equal to the dissipated power of internal energy. The formula can be expressed as follows (Chen 2007; Yang and Wang 2018):

${\displaystyle {\int }_V{\sigma }_{ij}}{\dot{\epsilon }}_{ij}^{\text{*}}\text{d}V={\displaystyle {\int }_S{T}_i}{v}_i^{\text{*}}\text{d}S+{\displaystyle {\int }_V{F}_i}{v}_i^{\text{*}}\text{d}V,(1)$

where σ_ij and ɛ̇_ij are the stress state and the volume change rate at any point in the plastic failure zone, respectively; V and F_i are the volume and volume force of the unit in the plastic failure zone respectively; S and T_i are the surface area and surface force of the element in the plastic failure zone, respectively; and v_i^* refers to the velocity along the sliding surface.

3 3D Failure Mechanism

Figure 1 shows the 3D rotational failure mechanism for saturated soft clay slopes. Rotating the failure surface around O point with angular velocity ω, it passes through A point at the top of the slope and C point at the toe of the slope, forming two logarithmic spiral curves of AC and A´C´, where the angle between OA and OA´ and the horizontal direction is θ_A, the lengths are r₀ and r₀´, respectively, the radius of C is r_h, the angle of OB is θ_B, the angle between the direction diameter and the horizontal direction of any point on arcs AC and A´C´ is θ, and the lengths are r and r´, respectively. As shown in Figure 1, the plane perpendicular to the paper surface at O point intersects with the 3D rotational failure mechanism to form a circle; a rectangular coordinate system is established on the circle, with the positive direction of the Y-axis along the radial direction. Here, R is the radius of the circle and r_m is the radius of the center of the circle.

FIGURE 1

FIGURE 1. 3D rotational failure mechanism.

According to geometry,

$r=r_{\text{0}}\exp \left[\left(\theta -{\theta }_0\right)\tan \phi \right](2)$

$r^{\prime }=r_0^{\prime }\exp \left[-\left(\theta -{\theta }_0\right)\tan \phi \right](3)$

$r_{\text{m}}=\frac{r+r^{\prime }}{2}=r_{\text{0}}{f}_1(4)$

$R=\frac{r-r^{\prime }}{2}=r_{\text{0}}{f}_2(5)$

$\frac{H}{r_0}={\mathrm{e}}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }\cdot \sin \left({\theta }_{\text{h}}\right)-\sin \left({\theta }_{\text{0}}\right),(6)$

where f₁ and f₂ are, respectively:

$f_1=\frac{1}{2}\left[{\mathrm{e}}^{\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }+\frac{r_0^{\prime }}{r_0}{\mathrm{e}}^{-\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }\right](7)$

$f_2=\frac{1}{2}\left[{\mathrm{e}}^{\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }-\frac{r_0^{\prime }}{r_0}{\mathrm{e}}^{-\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }\right].(8)$

In order to reflect the state of saturated soft clay slope failure more truly, a block with width b is inserted into the middle of the 3D rotational failure mechanism, as shown in Figure 2; as b→+∞, the 3D failure mechanism degenerates to a 2D plane strain condition. The combined failure mechanism can be used to analyze the 3D stability of saturated soft clay slopes of any width; in the analysis of this study, let the maximum width of the 3D rotational failure mechanism with insert block be b₂, and the maximum width of the failure mechanism without insert block be b₁, then b + b₁ = b₂.

FIGURE 2

FIGURE 2. 3D failure mechanism. (A) 3D failure mechanism without plane insert, (B) 3D failure mechanism with plane insert.

4 Calculation of the Work Rates

In this study, the saturated soft clay slope with a single soil layer is calculated, and the following assumptions are made: 1) the soil mass satisfies the ideal elastic–plastic model and complies with the Mohr–Coulomb failure criterion and the associated flow law, and 2) for saturated soft clay, the calculation assumes that its internal friction angle value is close to 0.

4.1 The External Work Rates

In the calculation, the 3D rotational failure mechanism of the saturated soft clay slope is divided into two parts: the 3D rotational failure mechanism without the plane insert and with the plane insert, respectively; finally, the two are superimposed. The external work rates are composed of the work rates of soil weight and seismic forces.

1) The work rates of soil weight

The work rates of soil weight W_γ are divided into two parts for the calculation: the first part is the 3D rotational failure mechanism without plane insert W_γ-3D.

$\begin{array}{l}{W}_{\gamma -3\text{D}}=2\omega \gamma \cdot [{\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{x_1^{\ast }}{\displaystyle {\int }_{d_1}^{y^{\ast }}{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta +{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{x_2^{\ast }}{\displaystyle {\int }_{d_2}^{y^{\ast }}{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta ]\\ =\gamma \omega r_{\text{0}}^4{g}_1,\end{array}(9)$

where γ is the soil weight, and g₁ is the intermediate variable.

The other part is the work rate of soil weight with plane insert W_γ-insert.

$\begin{array}{l}{W}_{\gamma -\text{insert}}=2\omega \gamma \cdot [{\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{\frac{b}{2}}{\displaystyle {\int }_{d_1}^R{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta +{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{\frac{b}{2}}{\displaystyle {\int }_{d_2}^R{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta ]\\ =\gamma \omega r_{\text{0}}^4{g}_2\end{array}(10)$

where g₂ is the intermediate variable.

In summary, the work rate of the total soil weight W_γ in the 3D rotational failure mechanism with a plane insert is as follows:

$W_{\gamma }=W_{\gamma -3\text{D}}+W_{\gamma -\text{insert}},(11)$

where

$x_i^{\ast }=\sqrt{R^2-d_i^2}\left(i=\mathrm{1,2}\right)(12)$

$y^{\ast }=\sqrt{R^2-x^2}(13)$

$d_1=\frac{r_0\sin {\theta }_0}{\sin \theta }-r_{\text{m}}=r_0{f}_3(14)$

$\begin{array}{l}{d}_2=\frac{r_0\sin \left({\theta }_{\text{h}}+\beta \right)}{\sin \left(\theta +\beta \right)}{e}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }-r_{\text{m}}\\ =r_0{f}_4\end{array}(15)$

${\theta }_{\text{B}}=\arctan \frac{\sin {\theta }_{\text{0}}}{\cos {\theta }_0-l}(16)$

$l=\frac{\sin \left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)}{\sin {\theta }_{\text{h}}}-\frac{\sin \left({\theta }_{\text{h}}+\beta \right)}{\sin {\theta }_{\text{h}}\sin \beta }\left(\sin {\theta }_{\text{h}}{e}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }-\sin {\theta }_{\text{0}}\right),(17)$

where f₃ and f₄ are dimensional functions of 1.

2) The work rates of seismic force

In this article, the pseudo-static method is used to analyze the stability of the three-dimensional slope. According to the pseudo-static method, the horizontal seismic force F_kh = k_hG and the vertical seismic force F_kv = k_vG, where G is soil gravity and k_h is the horizontal seismic coefficient. According to the Code for Seismic Design of Buildings GB50011-2010 (China), its value generally ranges from 0 to 0.3; k_v is the vertical seismic coefficient, k_v = ζk_h, and ζ is the vertical earthquake proportional coefficient. Because the effect of vertical earthquake is less than that of horizontal earthquake, its value generally ranges from −1 to 1. When ζ > 0, the vertical earthquake force is downward.

The work rate by the horizontal seismic forces in the 3D rotational failure mechanism without the plane insert $W_{k_{\text{h}}-3\text{D}}$ is as follows:

$\begin{array}{l}{W}_{k_{\text{h}}-3\text{D}}=2k_{\text{h}}\omega \gamma \cdot [{\displaystyle {\int }_{{\theta }_0}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{x_1^{\ast }}{\displaystyle {\int }_{d_1}^{y^{\ast }}{\left(r_{\text{m}}+y\right)}^2\sin \theta \text{d}x\text{d}y}}}\text{d}\theta +{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{x_2^{\ast }}{\displaystyle {\int }_{d_2}^{y^{\ast }}{\left(r_{\text{m}}+y\right)}^2\sin \theta \text{d}x\text{d}y}}}\text{d}\theta ]\\ =k_{\text{h}}\gamma \omega r_{\text{0}}^4{g}_3,\end{array}(18)$

where g₃ is the intermediate variable.

Similarly, the work rate generated by the horizontal seismic forces acting on plane insert $W_{k_{\text{h}}-\text{insert}}$ is as follows:

$\begin{array}{l}{W}_{k_{\text{h}}-\text{insert}}=2k_{\text{h}}\omega \gamma \cdot [{\displaystyle {\int }_{{\theta }_0}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{\frac{b}{2}}{\displaystyle {\int }_{d_1}^R{\left(r_{\text{m}}+y\right)}^2\sin \theta \text{d}x\text{d}y}}}\text{d}\theta +{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{\frac{b}{2}}{\displaystyle {\int }_{d_2}^R{\left(r_{\text{m}}+y\right)}^2\sin \theta \text{d}x\text{d}y}}}\text{d}\theta ]\\ =k_{\text{h}}\gamma \omega r_{\text{0}}^4{g}_4,\end{array}(19)$

where g₄ is the intermediate variable.

In summary, the work rates generated by the total horizontal seismic forces in the 3D rotational failure mechanism without the plane insert $W_{k_{\text{h}}}$ are as follows:

$W_{k_{\text{h}}}=W_{k_{\text{h}}-3\text{D}}+W_{k_{\text{h}}-\text{insert}}.(20)$

Similar to the work done by the horizontal seismic forces, the work rates of vertical seismic forces are also composed of two parts: the 3D rotational failure mechanism without plane insert $W_{k_{\text{v}}-3\text{D}}$ is given as follows:

$\begin{array}{l}{W}_{k_{\text{v}}-3\text{D}}=2\zeta k_{\text{h}}\omega \gamma \cdot [{\displaystyle {\int }_{{\theta }_0}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{x_1^{\ast }}{\displaystyle {\int }_{d_1}^{y^{\ast }}{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta +{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{x_2^{\ast }}{\displaystyle {\int }_{d_2}^{y^{\ast }}{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta ]\\ =\zeta k_{\text{h}}\gamma \omega r_{\text{0}}^4{g}_1,\end{array}(21)$

the work rate generated by the horizontal seismic forces acting on plane insert $W_{k_{\text{v}}-\text{insert}}$ is given as follows:

$\begin{array}{l}{W}_{k_{\text{v}}-\text{insert}}=2\zeta k_{\text{h}}\omega \gamma \cdot [{\displaystyle {\int }_{{\theta }_0}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{\frac{b}{2}}{\displaystyle {\int }_{d_1}^R{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta +{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{\frac{b}{2}}{\displaystyle {\int }_{d_2}^R{\left(r_{\text{m}}+y\right)}^2\cos \theta \text{d}x\text{d}y}}}\text{d}\theta ]\\ =\zeta k_{\text{h}}\gamma \omega r_{\text{0}}^4{g}_2,\end{array}(22)$

Therefore, the total work rates of vertical seismic forces, which act on the 3D rotational failure mechanism without plane insert $W_{k_{\text{v}}}$, are as follows:

$W_{k_{\text{v}}}=W_{k_{\text{v}}-3\text{D}}+W_{k_{\text{v}}-\text{insert}},(23)$

In summary, the work rates of soil weight and seismic forces in the 3D rotational failure mechanism with plane insert, namely, the external work rates can be expressed as follows:

$W=W_{\gamma }+W_{k_{\text{h}}}+W_{k_{\text{v}}},(24)$

4.2 Internal Energy Dissipation

The internal energy dissipation rates of the 3D rotational failure mechanism without the plane insert are calculated, and the saturated soft clay slope is divided into AB and BC parts along the extension line of OB, the internal energy dissipation rates are D_AB-3D and D _BC-3D, respectively:

$\begin{array}{l}{D}_{\text{AB}-3\text{D}}=2\omega cR\cdot {\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{B}}}{\displaystyle {\int }_0^{{\alpha }_1^{\ast }}{\left(r_{\text{m}}+R\cos \alpha \right)}^2\text{d}\alpha \text{d}\theta }}\\ =\omega c{r}_0^3{g}_5\end{array}(25)$

$\begin{array}{l}{D}_{\text{BC}-3\text{D}}=2\omega cR\cdot {\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}{\displaystyle {\int }_0^{{\alpha }_2^{\text{*}}}{\left(r_{\text{m}}+R\cos \alpha \right)}^2\text{d}\alpha \text{d}\theta }}\\ =\omega c{r}_0^3{g}_6,\end{array}(26)$

where c is the cohesion of soil, and g₅ and g₆ are intermediate variables, respectively.

The internal energy dissipated rates acting on the plane insert D _insert are as follows:

$\begin{array}{l}{D}_{\text{insert}}=b{\displaystyle {\int }_{{\theta }_0}^{{\theta }_{\text{h}}}cv\cos \phi \frac{r}{\cos \phi }\text{d}\theta }\\ =\omega c{r}_0^3{g}_7,\end{array}(27)$

where v is the velocity and g₇ is the intermediate variable.

In summary, the total internal energy dissipated rates can be expressed as follows:

$D=D_{\text{AB}-3\text{D}}+D_{\text{BC}-3\text{D}}+D_{\text{insert}},(28)$

where

${\alpha }_1^{\ast }=\arccos \left(d_1/R\right),(29)$

${\alpha }_2^{\ast }=\arccos \left(d_2/R\right),(30)$

$v=\omega \rho =\left(r_{\text{m}}+y\right)\omega ,(31)$

where ρ is the direction diameter of any point and ω is the angular velocity at this point.

4.3 Upper Bound Solution

As the internal energy dissipation rate and external work rate are equal, the critical height of the saturated soft clay slope can be obtained as follows:

$H_{\text{c}}=\frac{c\cot \phi }{\gamma }\cdot \left[{\mathrm{e}}^{({\theta }_{\text{h}}-{\theta }_{\text{0}})\tan \phi }\sin {\theta }_{\text{h}}-\sin {\theta }_{\text{0}}\right]\cdot \frac{g_5+g_6+g_7}{\left(1+\zeta k_{\text{h}}\right)\left(g_1+g_2\right)+k_{\text{h}}\left(g_3+g_4\right)},(32)$

In order to analyze the stability of the 3D saturated soft clay slope, the stability factor N_L is introduced in this article. Under the Mohr–Coulomb failure criterion, the stability factor of the saturated soft clay slope can be defined as follows:

$N_{\text{L}}=\gamma H_{\text{c}}/c,(33)$

The stability of the slope is different with the sliding surface. Therefore, we need to take different θ₀ and θ_h to calculate different N_L according to Eqs. 32, 33. The minimum N_L is the stability factor of the slope, and the corresponding θ₀ and θ_h determine the slip surface of the slope. The formula to solve the stability factor N_L is complicated, and it needs to be calculated by MATLAB software.

In addition, in the calculation, f₃∼ f₁₀ and g₁ ∼ g₇ are functions of dimension 1, and their expressions are, respectively, given as follows:

$f_3=\frac{\sin {\theta }_0}{\sin \theta }-\frac{1}{2}\left[{\mathrm{e}}^{\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }+\frac{r_0^{\prime }}{r_0}{\mathrm{e}}^{-\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }\right](34)$

$f_4=\frac{\sin \left({\theta }_{\text{h}}+\beta \right)}{\sin \left(\theta +\beta \right)}{e}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }-\frac{1}{2}\left[{\mathrm{e}}^{\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }+\frac{r_0^{\prime }}{r_0}{\mathrm{e}}^{-\left(\theta -{\theta }_{\text{0}}\right)\tan \phi }\right](35)$

$f_5=\frac{1}{3\left(1+9{\tan }^2\phi \right)}\cdot [\left(3\tan \phi \cos {\theta }_{\text{h}}+\sin {\theta }_{\text{h}}\right)e^{3\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }-\left(3\tan \phi \cos {\theta }_{\text{0}}+\sin {\theta }_{\text{0}}\right)](36)$

$f_6=\frac{1}{6}\left(2\cos {\theta }_{\text{0}}-\frac{L}{r_0}\right)\frac{L\sin {\theta }_{\text{0}}}{r_0}(37)$

$f_7=\frac{H}{3r_0}\cdot [\left({\cos }^2{\theta }_{\text{h}}+\sin {\theta }_{\text{h}}\cos {\theta }_{\text{h}}\cot \beta \right){\mathrm{e}}^{2\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }+\frac{H}{2r_0}\left(\cos {\theta }_{\text{h}}\cot \beta +\sin {\theta }_{\text{h}}{\cot }^2\beta \right){\mathrm{e}}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }](38)$

$f_8=\frac{1}{3\left(1+9{\tan }^2\phi \right)}\cdot [\left(3\tan \phi \sin {\theta }_{\text{h}}-\cos {\theta }_{\text{h}}\right)e^{3\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }-3\tan \phi \sin {\theta }_{\text{0}}+\cos {\theta }_{\text{0}}](39)$

$f_9=\frac{1}{3}\frac{L}{r_0}{\sin }^2{\theta }_{\text{0}}(40)$

$f_{10}=\frac{H}{3r_0}\cdot [\left(\cos {\theta }_{\text{h}}\sin {\theta }_{\text{h}}+{\sin }^2{\theta }_{\text{h}}\cot \beta \right){\mathrm{e}}^{2\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }-\frac{H}{2r_0}\left(\cos {\theta }_{\text{h}}+\sin {\theta }_{\text{h}}\cot \beta \right){\mathrm{e}}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }](41)$

$\begin{array}{l}{g}_1=2{\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{B}}}[\left(\frac{f_2^2{f}_3}{8}-\frac{f_3^3}{4}-\frac{2f_1{f}_3^2}{3}-\frac{f_3{f}_1^2}{2}+\frac{2f_1{f}_2^2}{3}\right)\cdot \sqrt{f_2^2-f_3^2}+\left(\frac{f_2^4}{8}+\frac{f_2^2{f}_1^2}{2}\right)\cdot \arcsin \left(\frac{\sqrt{f_2^2-f_3^2}}{f_2}\right)]\text{cos}\theta }\text{d}\theta \\ +2{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}[\left(\frac{f_2^2{f}_4}{8}-\frac{f_4^3}{4}-\frac{2f_1{f}_4^2}{3}-\frac{f_4{f}_1^2}{2}+\frac{2f_1{f}_2^2}{3}\right)\cdot \sqrt{f_2^2-f_4^2}+\left(\frac{f_2^4}{8}+\frac{f_2^2{f}_1^2}{2}\right)\cdot \arcsin \left(\frac{\sqrt{f_2^2-f_4^2}}{f_2}\right)]\text{cos}\theta }\text{d}\theta \end{array}(42)$

$g_2=\frac{b}{H}\left(f_5-f_6-f_7\right)\cdot \left[{\mathrm{e}}^{\left({\theta }_{\text{h}}-{\theta }_{\text{0}}\right)\tan \phi }\sin {\theta }_{\text{h}}-\sin {\theta }_{\text{0}}\right](43)$

$\begin{array}{l}{g}_3=2{\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{B}}}[\left(\frac{f_2^2{f}_3}{8}-\frac{f_3^3}{4}-\frac{2f_1{f}_3^2}{3}-\frac{f_3{f}_1^2}{2}+\frac{2f_1{f}_2^2}{3}\right)\cdot \sqrt{f_2^2-f_3^2}+\left(\frac{f_2^4}{8}+\frac{f_2^2{f}_1^2}{2}\right)\arcsin \left(\frac{\sqrt{f_2^2-f_3^2}}{f_2}\right)]\text{sin}\theta }\text{d}\theta \\ +2{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}[\left(\frac{f_2^2{f}_4}{8}-\frac{f_4^3}{4}-\frac{2f_1{f}_4^2}{3}-\frac{f_4{f}_1^2}{2}+\frac{2f_1{f}_2^2}{3}\right)\cdot \sqrt{f_2^2-f_4^2}+\left(\frac{f_2^4}{8}+\frac{f_2^2{f}_1^2}{2}\right)\arcsin \left(\frac{\sqrt{f_2^2-f_4^2}}{f_2}\right)]\text{sin}\theta }\text{d}\theta \end{array}(44)$

$g_4=\frac{b}{H}\left(f_8-f_9-f_{10}\right)\cdot \left[{\mathrm{e}}^{({\theta }_{\text{h}}-{\theta }_{\text{0}})\tan \phi }\sin {\theta }_{\text{h}}-\sin {\theta }_{\text{0}}\right](45)$

$g_5=2{\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{B}}}[f_2\left(f_1^2+\frac{f_2^2}{2}\right)\arccos \left(\frac{f_3}{f_2}\right)}+f_2\left(2f_1+\frac{f_3}{2}\right)\sqrt{f_2^2-f_3^2}]\text{d}\theta (46)$

$g_6=2{\displaystyle {\int }_{{\theta }_{\text{B}}}^{{\theta }_{\text{h}}}[f_2\left(f_1^2+\frac{f_2^2}{2}\right)\arccos \left(\frac{f_4}{f_2}\right)}+f_2\left(2f_1+\frac{f_4}{2}\right)\sqrt{f_2^2-f_4^2}]\text{d}\theta (47)$

$g_7=\frac{b}{H}\left[{\mathrm{e}}^{({\theta }_{\text{h}}-{\theta }_{\text{0}})\tan \phi }\sin {\theta }_{\text{h}}-\sin {\theta }_{\text{0}}\right]\cdot {\displaystyle {\int }_{{\theta }_{\text{0}}}^{{\theta }_{\text{h}}}{\mathrm{e}}^{2(\theta -{\theta }_{\text{0}})\tan \phi }}\text{d}\theta ,(48)$

5 Comparison

In order to simplify the study, the research object of this article is the same as Michalowski and Drescher (2009), and the friction problem of complex clays is not considered. Without considering the seismic load (k_h = 0, k_v = 0), the stability factor N_L of the saturated soft clay slope (φ = 0) is calculated. Table 1 shows the comparisons between the stability factors calculated in the work and the results of Michalowski and Drescher (2009). As can be seen from Table 1, the results of the two methods are quite similar, and the maximum error is only 10.2%. It shows that calculation in this article is effective.

TABLE 1

TABLE 1. Comparing the results between this study and Michalowski and Drescher (2009).

6 Results Analysis

6.1 Influence of Slope Width-to-Height Ratio and Slope Angle

Considering the seismic load (k_h = 0.2, ζ = 0.5), the influence of B/H and β on stability factors is shown in Figure 3. As shown in Figure 3A, the slope stability factor N_L decreases with the increase in B/H, and the decreasing trend becomes more gentle and finally tends to level. As can be seen from Table 2, when B/H ≥ 8 and β = 45°, 60°, 75°, and 90°, the relative error between stability factors is only 1%. Assume that the relative error allowed is 1%, and the conclusions are as follows:

1) When B/H = 1 or 2, the three-dimensional effect of the slope is very prominent, and the relative error between stability factors can reach 21.4%.

2) As B/H increases from 1 to 7, the stability coefficient of the slope decreases gradually, and the three-dimensional effect of the slope becomes less and less prominent.

3) When B/H ≥ 8, the stability factor of the slope is approximate to a certain value, and the three-dimensional slope can be simplified as a two-dimensional problem.

FIGURE 3

FIGURE 3. Effects of the ratio B/H and slope angle β on stability factor N_L: (A) B/H; (B) β.

TABLE 2

TABLE 2. Effects of the ratio B/H on stability factor N_L (k_h = 0.2, ζ = 0.5).

As shown in Figure 3B, when β increases, the slope stability factor N_L decreases, and when B/H is larger (B/H = 10), N_L decreases more. Compared with the original slope, when the slope angle increases, the sliding surface does not change, but the soil mass of the slope increases. The internal energy dissipation rate remains unchanged, and the external work rate increases, so the slope stability factor N_L decreases. The larger the B/H is (such as B/H = 10), the more soil mass increases, and the greater the slope stability factor decreases.

As can be seen from Table 3, when the slope angle β decreases by 5° each time from 90° to 45°, the stability factor N_L approximately increases by 1/3. The slope angle has a significant effect on slope stability, and reasonable design of a smaller slope angle is particularly important, which greatly increases the safety of the slope and reduces the possibility of slope landslide.

TABLE 3

TABLE 3. Effects of slope angle β on stability factor N_L (k_h = 0.2, ζ = 0.5).

In conclusion, B/H and β significantly affect the stability factor N_L. The ratio B/H represents the three-dimensional effect degree of the slope. When B/H is large enough (B/H ≥ 8), the slope stability can be simplified into a two-dimensional problem.

${\Delta }_2=\left|\frac{N_{\text{L}}\left(B/H=n+1\right)-N_{\text{L}}\left(B/H=n\right)}{N_{\text{L}}\left(B/H=n\right)}\right|\times 100\%$

${\Delta }_3=\left|\frac{N_{\text{L}}\left(\beta =n\right)-N_{\text{L}}\left(\beta =90°\right)}{N_{\text{L}}\left(\beta =90°\right)}\right|\times 100\%.$

6.2 Influence of Horizontal Seismic Force

It can be seen from Figure 4, the stability factor N_L of the saturated soft clay slope decreases with the increase in k_h. As the critical height H_c decreases, both the rate of external work and the rate of internal energy dissipation decrease, but the work rate of external force decreases more. When the horizontal seismic force increases, the work rate of the three-dimensional slope increases. According to the conservation of energy, the critical height H_c of the slope decreases, so the stability factor N_L decreases.

FIGURE 4

FIGURE 4. Effects of the horizontal seismic effect coefficient k_h on stability factor NL: (A) B/H = 1, β = 60°; (B) B/H = 2, β = 60°; (C) B/H = 5, β = 60°; (D) B/H = 10, β = 60°; (E) B/H = 5, β = 45°; (F) B/H = 5, β = 75°.

As can be seen from Table 4, taking B/H = 1, β = 45° as an example, without considering the effect of horizontal seismic force (k_h = 0), the safety factor is 9.06. When k_h increased from 0.1 to 0.3, the safety factor N_L decreased from 7.98 to 6.33. Compared with k_h = 0, the minimum relative error was 11.9% and the maximum was 30.2%.

TABLE 4

TABLE 4. Effects of the horizontal seismic effect coefficient k_h on stability factor N_L (ζ = 0).

Taking B/H = 10, β = 45° as an example, without considering the effect of horizontal seismic force (k_h = 0), the safety factor is 6.15. When k_h increased from 0.1 to 0.3, the safety factor N_L decreased from 5.20 to 3.63. Compared with k_h = 0, the minimum relative error was 15.5% and the maximum was 41.0%.

It shows that the horizontal seismic force obviously affects the safety factors of saturated soft clay slope. Especially when B/H is large, the slope stability problem can be regarded as a two-dimensional problem. If the influence of the horizontal seismic force is ignored, the stability factor of the slope will be significantly overestimated.

$\begin{array}{l}{\Delta }_4=\left|\frac{N_{\text{L}}\left(k_{\text{h}}=0.1\right)-N_{\text{L}}\left(k_{\text{h}}=0\right)}{N_{\text{L}}\left(k_{\text{h}}=0\right)}\right|\times 100\%,{\Delta }_5=\left|\frac{N_{\text{L}}\left(k_{\text{h}}=0.2\right)-N_{\text{L}}\left(k_{\text{h}}=0\right)}{N_{\text{L}}\left(k_{\text{h}}=0\right)}\right|\times 100\%,\\ {\Delta }_6=\left|\frac{N_{\text{L}}\left(k_{\text{h}}=0.3\right)-N_{\text{L}}\left(k_{\text{h}}=0\right)}{N_{\text{L}}\left(k_{\text{h}}=0\right)}\right|\times 100\%\end{array}$

6.3 Influence of Vertical Seismic Force

It can be seen from Figure 5, when the vertical seismic effect proportion coefficient ζ increases, the stability factor N_L of the saturated soft clay slope decreases. When ζ > 0, the vertical seismic force is downward, which is not conducive to the stability of the slope. When the vertical seismic force increases, the work rate of the three-dimensional slope increases. According to the conservation of energy, the critical height H_c of the slope decreases, so the stability factor N_L decreases.

FIGURE 5

FIGURE 5. Effects of the vertical seismic effect proportion coefficient ζ on stability factor N_L. (A) B/H = 1, β = 60°; (B) B/H = 2, β = 60°; (C) B/H = 5, β = 60°; (D) B/H = 10, β = 60°; (E) B/H = 5, β = 45°; (F) B/H = 5, β = 75°.

As can be seen from Table 5, taking B/H = 1, β = 45° as an example, without considering the effect of vertical seismic force (ζ = 0), the safety factor is 7.08. When ζ increased from 0.5 to 1, the safety factor N_L decreased from 6.58 to 6.14. Compared with ζ = 0, the minimum relative error was 7.1% and the maximum was 13.2%.

TABLE 5

TABLE 5. Effects of the vertical seismic effect proportion coefficient ζ on stability factor N_L (k_h = 0.2).

Taking B/H = 10, β = 45° as an example, without considering the effect of vertical seismic force (ζ = 0), the safety factor is 4.36. When ζ increased from 0.5 to 1, the safety factor N_L decreased from 4.09 to 3.86. Compared with ζ = 0, the minimum relative error was 6.1% and the maximum was 11.5%.

It shows that the vertical seismic force also has great influence on the stability of the saturated soft clay slope. If the vertical seismic force effect is ignored, the accuracy of calculation results will be directly affected.

Horizontal and vertical seismic loads have significant effects on slope stability. Therefore, in practical engineering, the influence of earthquake on slope stability needs to be considered, and the slope structure must be designed according to relevant seismic design code in order to ensure the stability of the slope.

$\begin{array}{l}{\Delta }_7=\left|\frac{N_{\text{L}}\left(\zeta =-1\right)-N_{\text{L}}\left(\zeta =0\right)}{N_{\text{L}}\left(\zeta =0\right)}\right|\times 100\%,{\Delta }_8=\left|\frac{N_{\text{L}}\left(\zeta =-0.5\right)-N_{\text{L}}\left(\zeta =0\right)}{N_{\text{L}}\left(\zeta =0\right)}\right|\times 100\%\\ {\Delta }_9=\left|\frac{N_{\text{L}}\left(\zeta =0.5\right)-N_{\text{L}}\left(\zeta =0\right)}{N_{\text{L}}\left(\zeta =0\right)}\right|\times 100\%,{\Delta }_{10}=\left|\frac{N_{\text{L}}\left(\zeta =1\right)-N_{\text{L}}\left(\zeta =0\right)}{N_{\text{L}}\left(\zeta =0\right)}\right|\times 100\%\end{array},$

7 Conclusion

1) The safety factors calculated in this article are very close to the results of Michalowski and Drescher (2009). The maximum error is only 7.7%, which shows that the calculation method in this article is effective.

2) The ratio B/H represents the three-dimensional effect degree of the saturated soft clay slope. Considering the effect of seismic load (k_h = 0.2, ζ = 0.5), the conclusions are as follows:

(a) When B/H = 1 or 2, the three-dimensional effect of the slope is very obvious.

(b) When B/H increases from 1 to 7, the three-dimensional effect of the slope gradually weakens.

(c) When the B/H ≥ 8, the three-dimensional slope stability problem can be simplified to a two-dimensional problem.

3) The stability factor N_L increases with the decrease in slope angle β. When the slope angle decreases from 90° to 45°, the stability coefficient N_L increases approximately 1/3. It is suggested that slope angle should be considered in the design of saturated soft clay slope.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author Contributions

The overarching research goals were developed by BZ. BZ established the models and analyzed the formal. YJ analyzed the results and wrote the initial draft of the manuscript. BZ and ZL revised and polished the manuscript. HC offered proposals and polished the manuscript. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This study was funded by the National Natural Science Foundation of China (52004088 and 52074116) and the Science Foundation of Hunan University of Science and Technology (E52076).

Conflict of Interest

Author HC is employed by Foshan Transportation Science and Technology Limited Company.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alejano, L. R., Ferrero, A. M., Ramírez-Oyanguren, P., and Álvarez Fernández, M. I. (2011). Comparison of Limit-Equilibrium, Numerical and Physical Models of wall Slope Stability. Int. J. Rock Mech. Mining Sci. 48, 16–26. doi:10.1016/j.ijrmms.2010.06.013

CrossRef Full Text | Google Scholar

Alemdag, S., KayaMustafa, A. K., Karadag, M., Gurocak, Z., and Bulut, F. (2015). Utilization of the Limit Equilibrium and Finite Element Methods for the Stability Analysis of the Slope Debris: An Example of the Kalebasi District (NE Turkey). J. Afr. Earth Sci. 106, 134–146. doi:10.1016/j.jafrearsci.2015.03.010

CrossRef Full Text | Google Scholar

Ausilio, E., Conte, E., and Dente, G. (2000). Seismic Stability Analysis of Reinforced Slopes. Soil Dyn. Earthquake Eng. 19, 159–172. doi:10.1016/S0267-7261(00)00005-1

CrossRef Full Text | Google Scholar

Chen, W. F. (2007). Limit Analysis and Soil Plasticity. Florida: J. Ross Publishing.

Google Scholar

Gao, Y.-f., Ye, M., and Zhang, F. (2015). Three-dimensional Analysis of Slopes Reinforced with Piles. J. Cent. South. Univ. 22, 2322–2327. doi:10.1007/s11771-015-2757-6

CrossRef Full Text | Google Scholar

Gischig, V. S., Eberhardt, E., Moore, J. R., and Hungr, O. (2015). On the Seismic Response of Deep-Seated Rock Slope Instabilities - Insights from Numerical Modeling. Eng. Geology. 193, 1–18. doi:10.1016/j.enggeo.2015.04.003

CrossRef Full Text | Google Scholar

Gofar, N., and Rahardjo, H. (2017). Saturated and Unsaturated Stability Analysis of Slope Subjected to Rainfall Infiltration. MATEC Web Conf. 101, 05004. doi:10.1051/matecconf/201710105004

CrossRef Full Text | Google Scholar

Han, C.-y., Chen, J.-j., Xia, X.-h., and Wang, J.-h. (2014). Three-dimensional Stability Analysis of Anisotropic and Non-homogeneous Slopes Using Limit Analysis. J. Cent. South. Univ. 21, 1142–1147. doi:10.1007/s11771-014-2047-8

CrossRef Full Text | Google Scholar

Huang, C.-C., Tsai, C.-C., and Chen, Y.-H. (2002). Generalized Method for Three-Dimensional Slope Stability Analysis. J. Geotech. Geoenviron. Eng. 128128, 83610–83848. doi:10.1061/(ASCE)1090-024110.1061/(asce)1090-0241(2002)128:10(836)

CrossRef Full Text | Google Scholar

Jongmin, K., Rodrigo, S., and Lee, J. (2002). Stability Analysis of Complex Soil Slopes Using Limit Analysis. J. Geotechnical Geoenvironmental Eng. 128128, 5467–5557. doi:10.1061/(ASCE)1090-0241

CrossRef Full Text | Google Scholar

Leong, E. C., and Rahardjo, H. (2012). Two and Three-Dimensional Slope Stability Reanalyses of Bukit Batok Slope. Comput. Geotechnics 42, 81–88. doi:10.1016/j.compgeo.2012.01.001

CrossRef Full Text | Google Scholar

Li, A. J., Lyamin, A. V., and Merifield, R. S. (2009). Seismic Rock Slope Stability Charts Based on Limit Analysis Methods. Comput. Geotechnics 36, 135–148. doi:10.1016/j.compgeo.2008.01.004

CrossRef Full Text | Google Scholar

Michalowski, R. L., and Drescher, A. (2009). Three-dimensional Stability of Slopes and Excavations. Géotechnique 59 (10), 839–850. doi:10.1680/geot.8.p.136

CrossRef Full Text | Google Scholar

Michalowski, R. L., and Martel, T. (2011). Stability Charts for 3D Failures of Steep Slopes Subjected to Seismic Excitation. J. Geotech. Geoenviron. Eng. 137, 183–189. doi:10.1061/(asce)gt.1943-5606.0000412

CrossRef Full Text | Google Scholar

Michalowski, R. L., and Park, D. (2020). Stability Assessment of Slopes in Rock Governed by the Hoek-Brown Strength Criterion. Int. J. Rock Mech. Mining Sci. 127, 1–12. doi:10.1016/j.ijrmms.2020.104217

CrossRef Full Text | Google Scholar

National Standard of the People’s Republic of China (2010). Code for Seismic Design of Building (GB50011-2010). Beijing, China: China Architecture and Building Press.

Google Scholar

Nian, T.-K., Jiang, J.-C., Wang, F.-W., Yang, Q., and Luan, M.-T. (2016). Seismic Stability Analysis of Slope Reinforced with a Row of Piles. Soil Dyn. Earthquake Eng. 84, 83–93. doi:10.1016/j.soildyn.2016.01.023

CrossRef Full Text | Google Scholar

Rawat, S., and Gupta, A. K. (2016). Analysis of a Nailed Soil Slope Using Limit Equilibrium and Finite Element Methods. Int. J. Geosynth. Ground Eng. 2, 34. doi:10.1007/s40891-016-0076-0

CrossRef Full Text | Google Scholar

Selcuk, A., Ayberk, K., and Mustafa, K (2015). Utilization of the limit equilibrium and finite element methods for the stability analysis of the slope debris: An example of the Kalebasi District (NE Turkey). Afr. Earth Sci.. doi:10.1016/j.jafrearsci.2015.03.010

CrossRef Full Text | Google Scholar

Sahoo, S., Manna, B., and Sharma, K. G. (2016). Seismic Stability Analysis of Un-reinforced and Reinforced Soil Slopes. Geo-china. Int. Conf.. doi:10.1061/9780784480007.009

CrossRef Full Text | Google Scholar

VandenBerge, D. R., and McGuire, M. P. (2019). Practical Use of Modified Hoek-Brown Criterion for Soil Slope Stability Analysis. Geotech Geol. Eng. 37 (6), 5441–5455. doi:10.1007/s10706-019-00991-1

CrossRef Full Text | Google Scholar

Wang, L., Hu, W., Sun, D. A., and Li, L. (2019). 3D Stability of Unsaturated Soil Slopes with Tension Cracks under Steady Infiltrations. Int. J. Numer. Anal. Methods Geomech 43 (6), 1184–1206. doi:10.1002/nag.2889

CrossRef Full Text | Google Scholar

Yang, X. L., and Liu, Z. A. (2018). Reliability Analysis of Three-Dimensional Rock Slope. Geomechanics Eng. 15 (6), 1183–1191. doi:10.12989/gae.2018.15.6.1183

CrossRef Full Text | Google Scholar

Yang, X. L., and Wang, H. Y. (2018). Catastrophe Analysis of Active-Passive Mechanisms for Shallow Tunnels with Settlement. Geomechanics Eng. 15 (1), 621–630. doi:10.12989/gae.2018.15.1.621

CrossRef Full Text | Google Scholar

Yao, C., and Yang, X. (2017). Limit Analysis of Unsaturated Soil Slope Stability Considering Intermediate Principal Stress and Strength Nonlinearity. Geotech Geol. Eng. 35, 2053–2063. doi:10.1007/s10706-017-0226-8

CrossRef Full Text | Google Scholar

Zhang, F., Gao, Y., Wu, Y., Zhang, N., and Qiu, Y. (2016). Effects of Vertical Seismic Acceleration on 3D Slope Stability. Earthq. Eng. Eng. Vib. 15, 487–494. doi:10.1007/s11803-016-0338-9

CrossRef Full Text | Google Scholar