Surficial stability analysis of soil slope under seepage based on a novel failure mode

Abstract

Normally, the edge effects of surficial landslides are not considered in the infinite slope method for surficial stability analysis of soil slopes. In this study, the limit stress state and discrimination equation of an infinite slope under saturated seepage flow were analyzed based on the Mohr-Coulomb strength criterion. Therefore, a novel failure mode involving three sliding zones (upper tension zone, middle shear sliding zone, and lower compression zone) was proposed. Accordingly, based on the limit equilibrium analysis, a semi-analytical framework considering the edge effect for the surficial stability of a soil slope under downslope seepage was established. Subsequently, the new failure mode was verified via a numerical finite element analysis based on the reduced strength theory with ABAQUS and some simplified methods using SLIDE software. The results obtained by the new failure mode agree well with those obtained by the numerical analysis and traditional simplified methods, and can be efficiently used to assess the surficial stability of soil slopes under rainwater seepage. Finally, an evaluation of the infinite slope method was performed using the semi-analytical method proposed in this study. The results show that the infinite slope tends to be conservative because the edge effect is neglected, particularly when the ratio of surficial slope length to depth is relatively small.

Graphical abstract

Keywords

soil slope / seepage / surficial failure mode / stress state / edge effects

1 Introduction

Surficial failures of soil slopes in tropical zones are commonly triggered by heavy rainfall because rainwater infiltration into unsaturated soil will reduce or eliminate matric suction, which can decrease the soil shear strength [1]. Furthermore, heavy rainfall may cause positive pore water pressure in the upper layers of a soil slope when an impervious layer exists at some depth of the slope, further reducing the soil shear strength and making the slope more susceptible to failure [2].

Surficial failure is most frequently addressed by the well-known infinite slope method, which is widely used in practice [3,4]. This method satisfies the condition of complete static force equilibrium under infinite length assumptions, and involves straightforward calculation procedures [5]. Although surficial landslides assumed to be infinitely long are usually valid as long as the length of the sliding block is sufficiently large compared with its buried depth, this is rarely justified, particularly when applying an analytical method to evaluate the ratio of length to depth (L/z_w), which is suitable for calculating the safety factor. This is because the resistance effects at the edges play a leading role when L/z_w is small. Much attention has been focused on analyzing pore water pressure distribution due to different rainfall intensities, and proposing repair measurements for surficial slope failures, which are established using the infinite slope method [6–8]. However, the boundary conditions of surficial landslides have rarely been considered. This can be mainly because the assumption of an infinitely long length employed in the infinite slope method is considered valid in most cases of surficial landslides by actual observations and justification [9–11]. This problem has been addressed using the finite element method (FEM) to calibrate the infinite slope method, and finally proved that when L/z_w is small, the infinite slope method produces large errors owing to the irrationality of the infinite length assumption [12]. Specifically, when L/z_w exceeds 16, the analysis result of the infinite slope method is acceptable. The validity of an infinite length assumption was reviewed using the FEM to the possible range of natural soil slopes using random parameter exploration [13]. It was found that when L/z_w exceeds 25, the relative difference in the safety factor between the infinite slope method and FEM always converges to within 5%. However, interpreting the influence of individual variables on the ratio of length and critical depth is difficult, as the FEM only provides general results [13]. In addition, rainfall leading to water accumulation in tension cracks at the top of a slope is unfavorable to the stability of the slope. A series of achievements have been made to reveal the mechanism of hydraulic fracturing and crack propagation using the phase-field method in porous media [14–19].

This paper proposes a semi-analytical stability method based on the classical infinite slope method for the surficial failure of soil slopes subjected to saturated seepage parallel to the slope surface. The semi-analytical method based on limit equilibrium theory considers the surficial landslide edge resistance effect, which is neglected in the infinite slope method.

2 Infinite slope method

When heavy rain infiltrates and saturates the upper layers of a soil slope and produces seepage parallel to the surface by an impervious layer at some depth z_w, slope failure begins to occur. The surficial stability model is illustrated in Fig. 1.

The classical infinite slope method established under the condition that the steady saturated seepage flow is parallel to the slope surface has been described in a considerable number of studies [20–22]. However, in practical engineering, when rainwater infiltrates the surficial layer of a soil slope, the angle between the streamline and the horizontal line usually varies from 0° to 90°; nevertheless, the infinite slope method assuming seepage parallel to the slope surface is still widely used [23–25]. To deepen the understanding of semi-analytical analysis, a brief review is presented here. One of the critical prerequisites of the infinite slope equation is that the slope is assumed to be infinitely long relative to the vertical depth z_w of the sliding surface.

Figure 1(a) shows the geometrical parameters and forces acting on the unit length of a soil slice, which is located on an infinite slope with an inclined angle α. The depth of water infiltration is denoted by z_w. The side forces P_l and P_r are earth pressures at the two ends of the slice, and both are directed parallel to the slope surface and equal in magnitude but in opposite directions. No consideration needs to be taken regarding P_l and P_r because they do not show any effect under an infinite slope condition. Figure 1(b) illustrates the flow net, as well as its expression for calculating the piezometric heads at any depth and the hydraulic gradient between equipotential lines. The water pressure u at vertical depth z_w is:

(1)u={\gamma }_{\mathrm{w}}{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha ,

where γ_w is the unit weight of water.

Along the flow direction, the seepage force J acting on the slice can be expressed as follows:

(2)J={\gamma }_{\mathrm{w}}i{z}_{\mathrm{w}},

where i is the hydraulic gradient and is equal to sinα on the same flow line between any equipotential lines because of the same piezometric heads at the same depth.

Therefore, if we define the safety factor f_s of the infinite slope as the ratio of the total resistance force to the driving force of the slice, f_s can be given by:

(3)f_{\mathrm{s}}={\displaystyle \frac{R}{T}}={\displaystyle \frac{c^{\prime }+{\gamma }^{\prime }{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha \tan {\phi }^{\prime }}{W^{\prime }\cos \alpha \sin \alpha \mathrm{+}J\cos \alpha }},

where cʹ is the effective cohesion, φʹ is the effective inner friction angle, γ′ is the buyout unit weight of soil, and W′ is the effective weight of the slice.

The classical infinite slope analytical expression can be obtained by substituting Eqs. (1) and (2) into Eq. (3):

(4)f_{\mathrm{s}}={\displaystyle \frac{c^{\prime }+{\gamma }^{\prime }{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha \tan {\phi }^{\prime }}{{\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\cos \alpha \sin \alpha }},

where γ_sat is the saturated unit weight of soil.

Simultaneously, the pore water pressure coefficient r_u can be defined as follows:

(5)r_{\mathrm{u}}={\displaystyle \frac{u}{{\sigma }_{\mathrm{v}}}},

where σ_v is equal to z_wγ_sat.

By substituting Eq. (5) into Eq. (4), the expression for the safety factor can be obtained as follows:

(6)f_{\mathrm{s}}={\displaystyle \frac{c^{\prime }}{{\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\cos \alpha \sin \alpha }}\mathrm{+(}\mathrm{1}-{\displaystyle \frac{r_{\mathrm{u}}}{{\cos }^{\mathrm{2}}\alpha }}){\displaystyle \frac{\tan {\phi }^{\prime }}{\tan \alpha }}.

For the total stress analysis, the strength parameters φ_u = 0 and c_u≠ 0 for saturated soil. Equation (6) can be simplified as follows:

(7)f_{\mathrm{s}}={\displaystyle \frac{c_{\mathrm{u}}}{{\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\cos \alpha \sin \alpha }}.

The failure model of the unit length slice can represent the entire infinite slope because vertical slices of soil of any length within the infinite slope are the same as the other vertical slices in all respects. Meanwhile, Eq. (4) is established based on whether the shear stress at the bottom of the unit length slice reaches the shear strength, which neglects the resistance effects from the upper and lower edges existing in actual surficial failures of natural landslides.

3 Modified method for surficial stability analysis

3.1 Upper and lower edge failure surfaces

The inclination angle of the infinite slope is denoted as α, as shown in Fig. 2(a). At depth z_w, the stress state of point B can be expressed easily under no seepage flow conditions [21,25]. The pressures p_v and p_l are the vertical stress and the lateral stress parallel to the slope surface, respectively, and they are conjugate stresses [21]. Similarly, p_v may also be substituted by the normal stress σ_n and shear stress τ_n [26–28]. Accordingly, the stress state under downslope seepage flow is influenced by the pore water pressure. The effective normal stress σ_n′ and shear stress τ_n, which acts on plane BB′, can be calculated using:

(8){\sigma }_n^{\prime }={\gamma }^{\prime }{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha ,

(9){\tau }_n={\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\cos \alpha \sin \alpha .

The lateral stress p_l′, which is conjugated to p_v′, acts on the plane parallel to the direction of p_v′, and the direction and magnitude of the resultant pressure p_v′ on plane BB′ has already changed when compared with the vertical pressure p_v, as shown in Fig. 2(c). According to Eqs. (8) and (9), the resultant pressure p_v′ can be expressed as follows:

(10){p^{\prime }}_{\mathrm{v}}=\sqrt{{({\sigma }_{\mathrm{n}}{}^{\prime })}^{\mathrm{2}}+{({\tau }_{\mathrm{n}})}^{\mathrm{2}}}={\gamma }^{\prime }{z}_{\mathrm{w}}\cos \alpha \sqrt{\mathrm{1}+{\sin }^{\mathrm{2}}\alpha \left[{\left({\displaystyle \frac{{\gamma }_{\mathrm{sat}}}{{\gamma }^{\prime }}}\right)}^{\mathrm{2}}-\mathrm{1}\right]}.

According to Eq. (10), p_v′ has a linear relationship with z_w and is represented by the distance OA in Fig. 3. The distance OB₁ represents the vertical stress p_v with its obliquity α, and the distance OB represents the pressure p_v′ with obliquity α′. The functional relationship between α′ and α is given by:

(11)\tan {\alpha }^{\prime }\mathrm{=}{\displaystyle \frac{{\gamma }_{\mathrm{sat}}}{{\gamma }^{\prime }}}\tan \alpha .

Considering Eq. (11), α′ is always greater than α. However, the hydraulic fracturing phenomenon induced by the seepage force will occur as the soil mass cannot sustain tension conditions, and α′ exceeds α₃ (α₃ = 45° + φ′/2). Hence, the geometry and mechanical condition of the entire surficial failure of the soil slope are α′<45° + φ′/2.

The corresponding depths of points B and B₁ are given by the critical depth z_cr of the infinite slope under the seepage condition parallel to the slope surface and no seepage flow, respectively. At point B, the shear stress τ_z is equal to the shear strength τ_f and the critical depth expression z_cr, shown in Eq. (12), can be derived, which is the same as Eq. (4) with a factor of safety f_s = 1. The z_cr value of point B is closer to the slope surface than that of point B₁.

(12)z_{\mathrm{cr}}={\displaystyle \frac{c^{\prime }}{{\gamma }_{\mathrm{sat}}\cos \alpha \sin \alpha -{\gamma }^{\prime }{\cos }^{\mathrm{2}}\alpha \tan {\phi }^{\prime }}}.

For a given finitely long slope, the failure surface must include the upper and lower edges; otherwise, surficial landslides cannot slip out. The rainwater infiltration to the depth z_w = z_cr is just a necessary condition, but not sufficient for surficial soil slope failure. This is because the unit length slice is always stable at any depth above z_w. Nonetheless, z_w is less than z_cr, and surficial failure occurs only on an infinitely long slope condition when depth z_w reaches z_cr. The vertical soil slice produces surplus slide forces to promote the upper and lower edges, resulting in expansion and compression failures when z_w is greater than z_cr.

For points A or A′ (Figs. 2(a) and 4) at the lower or upper edge, the soil will reach a limited state if Eqs. (13a) and (13b) are satisfied, which are the Mohr stress circle equation and limit state equation based on the Mohr-Coulomb failure criterion, respectively.

(13a){[{{\sigma }^{\prime }}_{\mathrm{z}}-({\displaystyle \frac{{{\sigma }^{\prime }}_{\mathrm{1}}+{{\sigma }^{\prime }}_{\mathrm{3}}}{\mathrm{2}}}\left)\right]}^{\mathrm{2}}+{\tau }_z^{\mathrm{2}}={\left({\displaystyle \frac{{{\sigma }^{\prime }}_{\mathrm{1}}-{{\sigma }^{\prime }}_{\mathrm{3}}}{\mathrm{2}}}\right)}^{\mathrm{2}},

(13b){{\sigma }^{\prime }}_{\mathrm{3}}={{\sigma }^{\prime }}_{\mathrm{1}}{K}_{\mathrm{a}}-\mathrm{2}{c}^{\prime }\sqrt{K_{\mathrm{a}}}.

Substituting Eqs. (8) and (9) into Eqs. (13a) and (13b), the two failure criterion equations of the soil mass at point A or A′ corresponding to the limit state are given by:

(14a){({p^{\prime }}_{\mathrm{v}})}^{\mathrm{2}}={\gamma }^{\prime }{\mathrm{z}}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha [{{\sigma }^{\prime }}_{\mathrm{3}}(\mathrm{1}+K_{\mathrm{p}})+\mathrm{2}{c}^{\prime }\sqrt{K_{\mathrm{p}}}]+(K_{\mathrm{p}}{{\sigma }^{\prime }}_{\mathrm{3}}^{\mathrm{2}}+\mathrm{2}{c}^{\prime }{{\sigma }^{\prime }}_{\mathrm{3}}\sqrt{K_{\mathrm{p}}}),

(14b){({p^{\prime }}_{\mathrm{v}})}^{\mathrm{2}}={\gamma }^{\prime }{\mathrm{z}}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha [{{\sigma }^{\prime }}_{\mathrm{1}}(\mathrm{1}+K_{\mathrm{a}})-\mathrm{2}{c}^{\prime }\sqrt{K_{\mathrm{a}}}]-(K_{\mathrm{a}}{{\sigma }^{\prime }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{2}{c}^{\prime }{{\sigma }^{\prime }}_{\mathrm{1}}\sqrt{K_{\mathrm{a}}}),

where σ₁′ and σ₃′ are the effective maximum and minimum principal stresses, respectively; and K_a = tan²(45° – φ′/2) and K_p = tan²(45° + φ′/2) are the active and passive earth pressure coefficients, respectively.

According to Eqs. (14a) and (14b), two limit Mohr circles can be obtained, such as C₁ and C₃ shown in Fig. 4. The circle with its center at point C₂ is an intermediate case that requires its center on the stress σ′ axis and cannot cross the strength envelope. Point A represents the stress p_v′ acting on the plane BB′ (Fig. 2(a)), the AB line has an angle α′ with the σ′ axis, and points G, E, and H represent the stress p_l′ (Fig. 4). Owing to the expansion by the slice of the middle shear sliding zone, the lateral pressure p_l′ reduces to the minimum (p_l′)_min, while point A at the upper edge reaches an actively tensile limited state. Because it is compressed by the slice of the middle shear sliding zone, p_l′ increases to the maximum (p_l′)_max. Meanwhile, point A at the lower edge reaches a passively compressive limited state.

In the triangles AC₁M and AC₃N, the expressions η and β with z_w can be derived based on the cosine law as follows:

(15a)\cos \beta =\mathrm{1}-\left({\displaystyle \frac{{[{{\sigma }^{\prime }}_z-{\displaystyle \frac{(\sqrt{\mathrm{\Delta }}-b)-\mathrm{2}{c}^{\prime }\tan {\phi }^{\prime }}{\mathrm{2}(\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }^{\prime })}}]}^{\mathrm{2}}+{[{\tau }_z-{\displaystyle \frac{(\sqrt{\mathrm{\Delta }}-b)\tan {\phi }^{\prime }+\mathrm{2}{c}^{\prime }}{\mathrm{2}(\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }^{\prime })}}]}^{\mathrm{2}}}{\mathrm{2}({[{{\sigma }^{\prime }}_z-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\sqrt{\mathrm{\Delta }}-b)]}^{\mathrm{2}}+{\tau }_z^{\mathrm{2}})}}\right),

(15b)\cos \eta =\mathrm{1}-\left({\displaystyle \frac{{[{{\sigma }^{\prime }}_z+{\displaystyle \frac{(\sqrt{\mathrm{\Delta }}+b)+\mathrm{2}{c}^{\prime }\tan {\phi }^{\prime }}{\mathrm{2}(\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }^{\prime })}}]}^{\mathrm{2}}+{[{\tau }_z+{\displaystyle \frac{(\sqrt{\mathrm{\Delta }}+b)\tan {\phi }^{\prime }-\mathrm{2}{c}^{\prime }}{\mathrm{2}(\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }^{\prime })}}]}^{\mathrm{2}}}{\mathrm{2}({[{{\sigma }^{\prime }}_z+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\sqrt{\mathrm{\Delta }}+b)]}^{\mathrm{2}}+{\tau }_z^{\mathrm{2}})}}\right),

where

(16)–b=\mathrm{2}c\prime \tan \phi \prime +\mathrm{2}{\sigma }_{\mathrm{z}}\prime (\mathrm{1}+{\tan }^{\mathrm{2}}\phi \prime ),

(17)\sqrt{\mathrm{\Delta }}=\sqrt{\mathrm{4}{[c^{\prime }\tan {\phi }^{\prime }+{{\sigma }^{\prime }}_z(\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }^{\prime })]}^{\mathrm{2}}-\mathrm{4}[(\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }^{\prime })({{\sigma }^{\prime }}_z^{\mathrm{2}}+{\tau }_z^{\mathrm{2}})-{c^{\prime }}^{\mathrm{2}}]},

and β ranges from 0° to (90° – φ′). If the influence of cracks on the soil slope is not considered due to tensile pressures at the upper edge, then η is located between 0° and (90° + φ′). Otherwise, the stress σ₃′ must exceed zero and η is between 0° and (90° + φ′ – 2α′).

Substituting σ₃′ = 0 into Eq. (14a) provides the expression of the crack depth z₀:

(18)z_{\mathrm{0}}={\displaystyle \frac{\mathrm{2}{c}^{\prime }}{{\gamma }^{\prime }\sqrt{K_{\mathrm{a}}}}}\cdot {\displaystyle \frac{\mathrm{1}}{{\cos }^{\mathrm{2}}\alpha +{[\sin \alpha ({\gamma }_{\mathrm{sat}}/{\gamma }^{\prime })]}^{\mathrm{2}}}}.

For the total stress analysis, Eq. (17) recovers to z_{\mathrm{0}}=\mathrm{2}c/\gamma \sqrt{K_{\mathrm{a}}} under no seepage flow conditions.

The surficial sliding mode under downslope seepage flow is shown in Fig. 5(a). The surficial landslide comprises three sliding zones: the upper tension zone, middle shear sliding zone, and lower compression zone. Because the surficial sliding mode mainly involves a planar translational slide, it can be defined as a “downslope curve” failure mode along a composite sliding surface.

3.2 Finite difference method

The finite difference method was used to express the shape of the upper and lower edge sliding surfaces, as shown in Fig. 5(b).

Taking the upper edge for example, the partial derivative of the function y = f(z_w, x) at point i can be represented as follows:

(19){\displaystyle \frac{\partial f(z_{\mathrm{w}},x_i)}{\partial x_i}}=-\cot (\mathrm{0}\mathrm{.5}{\eta }_i).

Neglecting the higher derivatives, Eq. (19) can be written in the first derivative difference form:

(20){\displaystyle \frac{\partial f(z_{\mathrm{w}},x_i)}{\partial x_i}}\approx {\displaystyle \frac{f(z_{\mathrm{w}},x_{i\mathrm{+1}})-f(z_{\mathrm{w}},x_{i-\mathrm{1}})}{\mathrm{2}\mathrm{\Delta }x}}.

The initial conditions for the difference in Eq. (20) are expressed as Eqs. (21) and (22).

(21){f(z_{\mathrm{w}},x)|}_{x=x_n}=\mathrm{0},

(22){{\displaystyle \frac{\partial f(z_{\mathrm{w}},x)}{\partial x}}|}_{x=x_i}=-\cot (\mathrm{0}\mathrm{.5}{\eta }_i).

Substituting Eqs. (21) and (22) into (20) leads to the slide length of the upper tension zone along the slope surface f(z_w, x₀) as follows:

(23)f(z_{\mathrm{w}},x_{\mathrm{0}})=\mathrm{2}\mathrm{\Delta }x[\cot (\mathrm{0.5}{\eta }_{\mathrm{1}})+\cot (\mathrm{0.5}{\eta }_{\mathrm{3}})+L+\cot (\mathrm{0.5}{\eta }_{\mathrm{n-3}})+\cot (\mathrm{0.5}{\eta }_{\mathrm{n-1}})].

Taking the lower edge as an example, the slide length of the lower compression zone along the slope surface f(z_w, y₀) can also be obtained as follows:

(24)f(z_{\mathrm{w}},y_{\mathrm{0}})=\mathrm{2}\mathrm{\Delta }y[\cot (\mathrm{0.5}{\beta }_{\mathrm{1}})+\cot (\mathrm{0.5}{\beta }_{\mathrm{3}})+L+\cot (\mathrm{0.5}{\beta }_{\mathrm{n-3}})+\cot (\mathrm{0.5}{\beta }_{\mathrm{n-1}})].

Thus, a set of slip lines is located in the upper tension and lower compression zones, as shown in Fig. 6(a).

3.3 Log spiral expression of sliding surface at the upper and lower edges

Although the sliding surface dimensions of the upper and lower edges can be calculated using the finite difference method, the stability analysis cannot be performed owing to a lack of functional expressions. The log spiral trace is well-known as “theoretically the best trace for slope stability analysis,” and its adequateness has been demonstrated for homogeneous slopes [29]. According to the sliding theory, the angle between two intersecting log spiral lines assumes a log spiral trace corresponding to the potential failure surface, and the formulation of the failure and the notation mechanism are shown in Fig. 6(b). The log-spiral equation can be expressed as follows:

(25)r=r_{\mathrm{0}}{\mathrm{e}}^{\theta \tan {\phi }^{\prime }},

where r and r₀ are the length variables of a log spiral, and θ is the angular variable of a log spiral.

Based on the log spiral geometry in Fig. 6(b), the length and angular variables of the log spiral functions r₁, r₂, r₃, θ₁, θ₂, and θ₃ can be solved. The total length of the surficial landslide mass is L = L_d + L_c + L_u. The lengths L_d and L_u represent the ranges of the upper and lower compression zones, respectively; and L_uc is the length of the upper tension zone, considering the influence of the crack (Fig. 6(a)).

3.4 Parametric analysis and examples

Using the finite difference method to verify the acceptability of the log spiral trace, Fig. 7 shows the sliding surfaces at the upper and lower edges changing with the effective cohesion c′ and the effective friction angle φ′, with the constant slope ratio of 1:1.5 and γ_sat = 20 kN/m³. In addition, the value of c′ is 0–20 kPa, and that of φ′ is 0–40°. First, the log spiral trace is reasonable, regardless of how the shear strength variables change. Second, when φ′ is constant, the potential failure surfaces are parallel to expand, r₃/r₁ is equal in magnitude, and the curve morphology does not change (Fig. 7(a)). If c′ is constant, r₃/r₁ increases with an increase in φ′, and the curve tends to develop linearly (Fig. 7(b)). This can be reasonably accounted for by the phenomenon that the slope failure surface of cohesionless soil develops along an approximate plane, but for purely cohesive soil, the slope failure surface develops along a circular surface.

For an impervious layer depth z_w greater than z_cr, the point stress of the location z_w has already yielded under the assumption of an infinite slope (Fig. 3), and the analytical equations in (15a) and (15b) about the relationship between the failure surface and ground surface would not be established. According to unsaturated soil mechanics [30], rainwater infiltration only reduces the matric suction strength, which is considered a part of cohesion when the deformation of the soil skeleton is neglected. Figure 7(a) illustrates the geomorphology of the potential failure surface and why r₃/r₁ does not vary with c′. Therefore, it is assumed that the log spiral variables r₃/r₁ at depth z_w are the same as those at depth z_cr.

4 Limit equilibrium analysis with three slide components

As described herein, the slide blocks in the upper and lower zones were studied on a differential slice within the length dx, as shown in Fig. 8(a). Similar to the Fellenius method [31] assumption, the differential slice side forces are parallel to the failure surface with the same magnitude but in the opposite direction. Moreover, the force was located at an elevation of h_i/3, as shown in Fig. 8(b). For the middle shear sliding zone, the lateral forces of the slice with length L_c are represented as P_n and P_{n+ 1}, as shown in Fig. 8(c).

The following equations can calculate the resistance force R_i and the normal force N_i along the base of slice i:

(26)R_i=(c\cdot ds+N_i\tan \phi )/F_{\mathrm{s}},

(27)N_i=W_i^{\prime }\cos {\delta }_i+J_i\sin (\alpha -{\delta }_i).

Around the pole of the log spiral, the two-moment equilibrium equations for the upper and lower compression zones can be expressed as follows, respectively:

(28){\displaystyle \sum M_{{\mathrm{W}}_i{}^{\prime }}^{\mathrm{u}}}+{\displaystyle \sum M_{{\mathrm{J}}_i}^{\mathrm{u}}}-{\displaystyle \sum M_{{\mathrm{N}}_i}^{\mathrm{u}}}-{\displaystyle \sum M_{{\mathrm{P}}_{\mathrm{n}}}^{\mathrm{u}}}-{\displaystyle \sum M_{{\mathrm{R}}_i}^{\mathrm{u}}}=\mathrm{0},

(29){\displaystyle \sum M_{{\mathrm{W}}_i{}^{\prime }}^{\mathrm{d}}}+{\displaystyle \sum M_{{\mathrm{J}}_i}^{\mathrm{d}}}-{\displaystyle \sum M_{{\mathrm{N}}_i}^{\mathrm{d}}}+{\displaystyle \sum M_{{\mathrm{P}}_{\mathrm{n+1}}}^{\mathrm{d}}}-{\displaystyle \sum M_{{\mathrm{R}}_{{}_i}}^{\mathrm{d}}}=\mathrm{0},

where the superscripts u and d represent the upper tension and lower compression zones, respectively, and i is the ID number of the slices. When i ranges from 1 to n, it represents the slice located at the upper tension zone; when i is located between n + 1 and m, it is for the slice located at the lower compression zone. M_{{\mathrm{W}}_i^{\prime }} is the moment due to W_i^{\prime }, M_{{\mathrm{N}}_i} is the moment due to N_i, M_{{\mathrm{J}}_i} is the moment due to J_i, M_{{\mathrm{R}}_i} is the moment due to R_i, M_{{\mathrm{P}}_{\mathrm{1}}} is the moment due to P_{\mathrm{n}}, and M_{{\mathrm{P}}_{\mathrm{3}}} is the moment due to P_{\mathrm{n+1}}.

Using Eqs. (28) and (29), the corresponding lateral forces P_{\mathrm{n}} and P_{\mathrm{n+1}} are easily obtained and can be represented by the following equations:

(30){\displaystyle \frac{{\displaystyle \sum J_{i}r\sin (\theta +\alpha )}+{\displaystyle \sum W_i{}^{\prime }r\cos {\phi }^{\prime }\sin {\delta }_i}-{\displaystyle \sum R_{i}r\cos {\phi }^{\prime }}}{d}}=P_{\mathrm{n}},

(31){\displaystyle \frac{{\displaystyle \sum R_{i}r\cos {\phi }^{\prime }}-{\displaystyle \sum J_{i}r\sin (\theta +\alpha )}-{\displaystyle \sum W_i{}^{\prime }r\cos {\phi }^{\prime }\sin {\delta }_i}}{d}}=P_{\mathrm{n}+\mathrm{1}},

where d is the perpendicular distance from P_{\mathrm{n}} or P_{\mathrm{n+1}} to the pole of the log spiral, which is equal to r_{\mathrm{2}}\cos {\phi }^{\prime }-\mathrm{1}/\mathrm{3}{z}_{\mathrm{w}}\cos \alpha.

The middle shear sliding zone force equilibrium equation can be expressed as (Fig. 8(c)):

(32){\displaystyle \frac{L_{\mathrm{c}}}{\cos \alpha }}(c^{\prime }+{\gamma }^{\prime }{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha \tan {\phi }^{\prime })/F_{\mathrm{s}}\mathrm{+}{P}_{\mathrm{n}+\mathrm{1}}={\displaystyle \frac{L_{\mathrm{c}}}{\cos \alpha }}{\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\cos \alpha \sin \alpha +P_{\mathrm{n}}.

Combining Eqs. (30)–(32), a generalized expression for the safety factor with the composite sliding mode included in the three zones can be expressed as follows:

(33)F_S={\displaystyle \frac{(c^{\prime }+{\gamma }^{\prime }{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha \tan {\phi }^{\prime })L_{\mathrm{c}}d/\cos \alpha +{\displaystyle \sum (c^{\prime }ds+N_i\tan \phi )r\cos {\phi }^{\prime }}}{L_{\mathrm{c}}d{\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\sin \alpha +{\displaystyle \sum W_i{}^{\prime }r\cos {\phi }^{\prime }\sin {\delta }_i}\mathrm{+}{\displaystyle \sum J_{i}r\sin (\theta +\alpha )}}}.

As the length L_c gradually increases to infinity, the above equation will recover to the infinite slope equation, as follows:

(34)\underset{L_{\mathrm{c}}\to \infty }{\lim }{F}_{\mathrm{s}}={\displaystyle \frac{c^{\prime }+{\gamma }^{\prime }{z}_{\mathrm{w}}{\cos }^{\mathrm{2}}\alpha \tan {\phi }^{\prime }}{{\gamma }_{sat}{z}_{\mathrm{w}}\cos \alpha \sin \alpha }}=f_{{}_{\mathrm{S}}}.

Thus, the infinite slope method is a particular case of the semi-analytical slope stability method for a composite sliding surface.

As L_c gradually decreases to zero, the “downslope curve” failure mode along the composite sliding surface degenerates to the entire log spiral failure mode. The equation for the safety factor now becomes:

(35)F_s={\displaystyle \frac{{\displaystyle \sum (c^{\prime }ds+N_i\tan {\phi }^{\prime })r\cos {\phi }^{\prime }}}{{\displaystyle \sum W_i{}^{\prime }r\cos {\phi }^{\prime }\sin {\delta }_i}\mathrm{+}{\displaystyle \sum J_{i}r\cos {\phi }^{\prime }\cos (\alpha -{\delta }_i)}}},

where

(36a){\displaystyle \sum (c^{\prime }ds+N_i\tan {\phi }^{\prime })r\cos {\phi }^{\prime }}={\displaystyle {\int }_{{\theta }_{\mathrm{1}}}^{{\theta }_{\mathrm{3}}}(c^{\prime }+h_i\sin (\theta -{\phi }^{\prime })[{\gamma }^{\prime }\sin (\theta -{\phi }^{\prime })\mathrm{+}{\gamma }_{\mathrm{w}}\sin \alpha \sin (\alpha -{\delta }_i)]\tan {\phi }^{\prime })r^{\mathrm{2}}d\theta },

(36b){\displaystyle \sum W_i{}^{\prime }r\cos {\phi }^{\prime }\sin {\delta }_i}={\displaystyle {\int }_{{\theta }_{\mathrm{1}}}^{{\theta }_{\mathrm{3}}}{h}_i{\gamma }^{\prime }\sin (\theta -\phi \text{'})\cos (\theta -{\phi }^{\prime })r^{\mathrm{2}}d\theta },

(36c){\displaystyle \sum J_{i}r\cos {\phi }^{\prime }\cos (\alpha -{\delta }_i)}={\displaystyle {\int }_{{\theta }_{\mathrm{1}}}^{{\theta }_{\mathrm{3}}}{h}_i{\gamma }_{\mathrm{w}}\sin \alpha \sin (\theta -\phi \text{'})\cos (\alpha -\delta )r^{\mathrm{2}}d\theta }.

For the total stress analysis: φ_u = 0; c_u soil, Eq. (33) is simplified to:

(37)F_S={\displaystyle \frac{c_{\mathrm{u}}d{L}_{\mathrm{c}}/\cos \alpha +{\displaystyle \sum c_{\mathrm{u}}rds}}{L_{\mathrm{c}}d{\gamma }_{\mathrm{sat}}{z}_{\mathrm{w}}\sin \alpha +{\displaystyle \sum W_{i}r\sin {\delta }_i}}}.

The downslope-curve composite sliding surface method is semi-analytical. This method starts with the stress state of an infinite-slope soil mass. It is based on both the Mohr-Coulomb strength criterion and the finite difference method to determine the dangerous sliding surface and apply the logarithmic spiral function to express the sliding surface. With this method, the expression of the safety factor is established such that the sliding bodies of the upper tension zone and lower compression zone meet the torque balance, and the sliding body of the middle shear sliding zone meets the force equilibrium. From analyzing the collapse of the surficial part of the soil slope, the semi-analytical method can obtain relatively accurate results because of the consistency of the upper and lower sliding body stress states with the theoretical analysis.

For soil slopes that generally have boundaries in the upper and lower parts, as shown in Fig. 8(a), the composite sliding surface position is difficult to determine using the actual states of stress. However, the fact that the sliding surface passes through the toe of the slope has already been accepted, and it is assumed that the lower edge compression zone shear outlet is at the toe of the slope. For the downslope-curve surficial slide, the entire slope surface slide is more unfavorable than the local surface slide by observing the model test [19]. It was demonstrated that the tension zone along the crack always reached the top of the slope. Thus, the tension slope can be precisely determined with the x- or y-coordinates of known B′. Moreover, starting point B′ of the tension zone can be estimated via a simple trial. Therefore, the position of the upper edge tension zone was determined using the single-variable search method. This approximate treatment can quickly evaluate the location of the downslope-curve failure surface.

5 Verification of the new failure mode

To verify the rationality of the semi-analytical method and the novel failure mode proposed in this paper, comparative analyses of the FEM using ABAQUS and a single circle search technique derived from the commercial software SLIDE were performed.

5.1 Comparison with numerical FEM analysis based on reduced-strength method

To facilitate the comparison, a case history reported by Griffiths et al. [12], as shown in Fig. 9, was utilized in this study, and the slope geometry and soil physical and mechanical parameters are set to the same as those in [12]. Specifically, the slope geometry parameters are z_w = 2.5 m, α = 30°, and L/z_w = 4; and the physical and mechanical parameters for the saturated zone are γ_sat = 20 kN/m³, c_u = 25 kPa, and φ_u = 0. Meanwhile, the model boundaries of the semi-analytical method are consistent with those of the model in [12]; namely, D = 2z_w takes bilateral boundaries, and the distance between the top and the bottom boundary is 2H, which can meet the requirements of calculation accuracy [32]. The slope-geometry boundary conditions are shown in Fig. 9.

Based on the semi-analytical method proposed herein, the safety factor of the slope model in Fig. 9 can be obtained as 1.903 using Eq. (37). In addition, a dangerous sliding surface can be deduced, as shown in Fig. 9. The upper tension zone is L_u = 2.400 m, the middle shear sliding zone is L_c = 3.422 m, and the lower compression zone L_d = 4.178 m.

A numerical study implemented using ABAQUS based on the finite-element reduced-strength method was conducted to demonstrate the comparison. Therefore, the slope model and setup for the numerical FEM are the same as those in the semi-analytical method. Other parameters used in the numerical analysis are listed in Table 1.

The base of the numerical model is fully fixed, and the extreme vertical boundaries to the left and right allow only vertical movement. A typical finite element mesh with 4-node bilinear plane-strain quadrilateral elements is shown in Fig. 10. In the numerical analysis, the elastic-plastic constitutive model, which obeys the Mohr-Coulomb yield criterion, is attributed to the slope soil. In addition, the unassociated flow rule (dilation angle ψ = 0.1°) was used to establish a finite element numerical solution of the slope stability.

Figure 11 shows the relationship between the displacement of the slope top point δ and the reduction factor. It can be observed that δ develops gradually initially and increases abruptly at some point. Based on the displacement mutation criterion, the safety factor F_s can be acquired as 1.938 based on the numerical analysis, which is very close to 1.927, as described by Griffiths et al. [12]. Therefore, the numerical analysis results are reliable and can be used for the comparison analysis.

For varying L/z_w, Table 2 lists the calculated results of the safety factor obtained using different methods, including FEM numerical analysis, semi-analytical method, and Griffiths et al.'s method, as well as the infinite slope theory. The results of the FEM numerical analysis, semi-analytical method, and Griffiths et al.'s method decrease gradually with the increment of L/z_w. In contrast, the result of the infinite slope method remains constant because it assumes that the slope length is infinite, which indicates that the results derived using the infinite slope theory tend to be conservative when L/z_w is relatively small.

Figures 12(a) and 12(b) illustrate the contour results of the equivalent plastic strain and displacement, respectively. It can be seen that the numerical failure mode appears to be a three-zone sliding mode, which agrees well with the result deduced by the semi-analytical method, as shown in Fig. 9. Therefore, the sliding surface derived from the semi-analytical method is reasonable.

5.2 Comparison of results derived from SLIDE

A single-center search strategy applies only one center to search for composite sliding surfaces in the SLIDE software, as shown in Fig. 13. When one sliding surface is produced, the intersections B and B′ of the circle line and the interfacial line of soil layers will be determined, and the arc line BB′ will be replaced by straight-line BB′. The traditional limit equilibrium method (Fellenius method [32] and simplified Bishop method [33]) is used to determine the safety factor of the composite surface with the combination lines of CB-BB′-C′B′ (i.e., the highlighted blue line in Fig. 10).

The soil properties were set as follows: γ_sat = 18.8 kN/m³, cʹ = 6.3 kPa, φʹ = 28.4° for silty sand soil located at University Technology Malaysia campus [2], z_w = 0.5~1.5 m for a typical surficial failure depth of soil slope caused by rainfall and H = 6 m with four slope angles: α = 30, 34 (slope ratio= 1.0:1.5), 40, and 45°.

Figure 14 only compares the dangerous sliding surfaces between the semi-analytical method and SLIDE search method under different infiltration depths when the slope angle α is 34°. Similar to Fig. 12(a), there is only a slight difference at the upper edge of the tensile zone. The calculated results for the safety factors are shown in Fig. 15. The results indicate that the safety factor decreases with the increase in infiltration depth z_w, and the difference between the results derived from the semi-analytical method and the SLIDE search methods, which involve the simplified Bishop method, Fellenius method, and simplified Janbu method [34], is small.

6 Evaluation of the infinite slope method

The key assumption of the classical infinite slope method is that the slope length L is infinite, and the safety factor is calculated using the change in slope height H, as shown in Fig. 16(a). This indicates that the safety factor obtained from the infinite slope method is independent of H (Table 1 and Fig. 13(a)), while that calculated using the semi-analytical method shows a gradual decrease with the increase in H, and is always higher than that of the former. When H tends to infinity, the semi-analytical method is reduced to an infinite slope method (as shown in Eq. (34)); however, when H is small, the results of the two methods are quite different. For instance, when H = 4 m and z_w = 1 m, the safety factor obtained using the latter method is 1.42 times that of the former. Therefore, to further examine the effectiveness of the infinite slope assumption, the relative difference ε between the predictions, which can represent the error of the infinite slope model due to neglecting the edge effects, is defined as Eq. (38).

(38)\epsilon ={\displaystyle \frac{F_{\mathrm{s}}-f_{\mathrm{s}}}{F_{\mathrm{s}}}}\times \mathrm{100}\%.

The relative difference between the safety factors from the infinite slope method is illustrated in Fig. 16(b) (α = 34°). The relationship between ε and λ is almost identical for different depths z_w. According to the project's needs, 5% or 10% were set as the allowable error, and the corresponding λ is denoted by λ_ε. As the relative difference ε equals 10%, the corresponding λ_ε ranges from 14 to 16. Then, ε increases steeply with a decrease in λ. When ε decreased to 5%, the corresponding λ_ε = 28–33. The semi-analytical method may then be recognized to recover to the infinite slope model, and the infinite length assumption is valid.

Other results are shown in Table 3, where ε = 5% and 10%, the range of the validity of the infinite slope assumptions λ_ε increases with a decrease in the slope angle α. At α = 45°, the ratio λ_ε is between 21~24 and 10~12. When the ratio λ_ε is less than 10, it will be too conservative to adopt the infinite slope method. Meanwhile, λ_ε is not constant but is also affected by the soil slope and geometric parameters.

7 Conclusions and discussion

A composite sliding mode, including the upper tension zone, middle shear sliding zone, and lower compression zone, is presented. According to the finite differential method, the upper and lower edge failure surfaces were determined and verified based on an established log spiral trace. Analysis of the failure surface morphology with a range of soil shear strengths illustrated that the failure surface of a soil slope is approximately straight for cohesionless soil and a circle for purely cohesive soil.

For the surficial stability of a soil slope with the consideration of edge effects, a semi-analytical framework has been proposed using composite sliding surface-based limit equilibrium analysis. The proposed semi-analytical method requires that the upper and lower sliding masses satisfy the moment equilibrium, and the middle sliding mass satisfies the force equilibrium. With the increase in the middle shear sliding zone length to infinitely long, the method recovers to the infinite slope method. However, with the decrease in the length of the middle shear sliding zone to zero, it degenerates to a log-spiral model. The semi-analytical method was well-verified via numerical FEM analysis and some simplified methods within the SLIDE software from safety factors and search for a dangerous sliding surface. Furthermore, the slip surface search pattern in the semi-analytical method only requires one variable to determine the critical slip surface, which is brief and precise and can be efficiently used to evaluate the surficial stability of soil slopes subjected to rainfall infiltration.

To assess the validity of the infinite slope method, a typical slope with silty sand soil was tested by assuming the relative difference between the proposed semi-analytical method and the infinite slope method mentioned later. The results indicate that the infinite slope method is conservative, resulting from neglection of the edge shear resistance effects, particularly when the length-to-depth ratio (L/z_w) is less than 10. When L/z_w increases to 21–38, the relative difference between predictions falls within 5%, and the edge effects relative to large ranges of the middle shear sliding zone can be neglected gradually.

References

[1] Wang X, Ma C, Wang Y, Wang Y, Li M, Dai Z, Li M. Effect of root architecture on rainfall threshold for slope stability: Variabilities in saturated hydraulic conductivity and strength of root-soil composite. Landslides, 2020, 17(8): 1965–1977

[2] Lee L M, Gofar N, Rahardjo H. A simple model for preliminary evaluation of rainfall-induced slope instability. Engineering Geology, 2009, 108(3–4): 272–285

[3] Si A, Zhang J, Zhang Y, Kazuva E, Dong Z, Bao Y, Rong G. Debris flow susceptibility assessment using the integrated random forest based steady-state infinite slope method: A case study in Changbai mountain, China. Water (Basel), 2020, 12(7): 2057

[4] Urciuoli G, Pirone M, Picarelli L. Considerations on the mechanics of failure of the infinite slope. Computers and Geotechnics, 2019, 107(3): 68–79

[5] Travis Q B, Houston S L, Marinho F A M, Schmeeckle M. Unsaturated infinite slope stability considering surface flux conditions. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(7): 963–974

[6] Simoni S, Zanotti F, Bertoldi G, Rigon G. Modelling the probability of occurrence of shallow landslides and channelized debris flows using GEOtop–FS. Hydrological Processes, 2008, 22(4): 532–545

[7] Feng S, Liu H W, Ng C W. Dimensional analysis of pore-water pressure response in a vegetated infinite slope. Canadian Geotechnical Journal, 2019, 56(8): 1119–1133

[8] Day R W. Design and repair for surficial slope failures. Practice Periodical on Structural Design and Construction, 1996, 1(3): 83–87

[9] Montgomery D R, Dietrich W E. A physically based model for the topographic control on shallow landsliding. Water Resources Research, 1994, 30(4): 1153–1171

[10] Haneberg W C. A rational probabilistic method for spatially distributed landslide hazard assessment. Environmental & Engineering Geoscience, 2004, 10(1): 27–43

[11] Warburton J, Milledge D G, Johnson R. Assessment of shallow landslide activity following the January 2005 storm. In: Cumberland Geological Society proceedings. Northern Cumbria: Cumberland Geological Society, 2008, 7: 263–283

[12] Griffiths D V, Huang J, Dewolfe G F. Numerical and analytical observations on long and infinite slopes. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(5): 569–585

[13] Milledge D G, Griffiths D V, Lane S N, Warburton J. Limits on the validity of infinite length assumptions for modeling shallow landslides. Earth Surface Processes and Landforms, 2012, 37(11): 1158–1166

[14] Zhuang X, Zhou S, Sheng M, Li G. On the hydraulic fracturing in naturally-layered porous media using the phase field method. Engineering Geology, 2020, 266: 105306

[15] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 169–198

[16] Zhou S, Zhuang X, Rabczuk T. Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation. Computer Methods in Applied Mechanics and Engineering, 2019, 355: 729–752

[17] Zhou S, Rabczuk T, Zhuang X. Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies. Advances in Engineering Software, 2018, 122: 31–49

[18] Zhou S, Zhuang X, Rabczuk T. A phase-field modeling approach of fracture propagation in poroelastic media. Engineering Geology, 2018, 240: 189–203

[19] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192

[20] Wiberg N E, Koponen M, Runesson K. Finite element analysis of progressive failure in long slopes. International Journal for Numerical and Analytical Methods in Geomechanics, 1990, 14(9): 599–612

[21] Taylor D W. Fundamentals of Soil Mechanics. New York: Wiley, 1948

[22] Day R W, Axten G W. Surficial stability of compacted clay slopes. Journal of Geotechnical Engineering, 1989, 115(4): 577–580

[23] Pradel D, Raad G. Effect of permeability on surficial stability of homogeneous slopes. Journal of Geotechnical Engineering, 1993, 119(2): 315–332

[24] Rosso R, Rulli M C, Vannucchi G. A physically based model for the hydrologic control on shallow landsliding. Water Resources Research, 2006, 42(6): 770–775

[25] Lade P V. The mechanics of surficial failure in soil slopes. Engineering Geology, 2010, 114(1–2): 57–64

[26] Wang L Z, Guo D J. Secondary stress analysis of slope tunnel. Mechanical Engineering, 2000, 22(4): 25–28 (in Chinese)

[27] Li W H, Peng L M, Lei M F, An Y L. Secondary and primary stress distribution of terrain bias tunnel. Journal of Railway Science and Engineering, 2012, 9(4): 63–69 (in Chinese)

[28] Gao X, Wang H N. Analytical solution of ground response induced by excavation of shallow-embedded bias tunnel under slope. Water Resources and Hydropower Engineering, 2019, 50(1): 1–9 (in Chinese)

[29] Zhao L, Cheng X, Li L, Chen J, Zhang Y. Seismic displacement along a log-spiral failure surface with crack using rock Hoek-Brown failure criterion. Soil Dynamics and Earthquake Engineering, 2017, 99: 74–85

[30] Fredlund D G, Rahardjo H. Soil Mechanics for Unsaturated Soils. New York: John Wiley & Sons Inc., 1993

[31] Zheng Y, Zhao S, Li A, Tang X. FEM Limit Analysis and Its Application in Slope Engineering. Bejing: China Communication Press, 2011 (in Chinese)

[32] Fellenius W. Calculation of the stability of earth dams. In: Transactions of the 2nd Congress on Large Dams. Washington D.C.: International Commission on Large Dams (ICOLD) Paris, 1936, 445–463

[33] Ji J, Zhang W, Zhang F, Gao Y, Lu Q. Reliability analysis on permanent displacement of earth slopes using the simplified Bishop method. Computers and Geotechnics, 2020, 117(1): 103286

[34] Janbu N. Stability Analysis of Slopes with Dimensionless Parameters. Cambridge: Harvard University, 1959

RIGHTS & PERMISSIONS

Higher Education Press