Comparative analysis of two models of progressive failure of slope

The physical significance of traditional slice method for slope stability analysis involves reducing the peak strength of slices to characterize the degree of strength degradation in the post-failure zone of geological materials. This study applies the partial strength reduction method, using the unbalanced thrust method as an example to extend its traditional assumptions. By employing both a perfect elasto-plastic model and a novel shear stress-strain complete process constitutive model method, an extended slice method is proposed to simulate the progressive failure process of landslides, the two methods for determining the critical stress state of slices are also presented. This approach tightly links traditional landslide stability analysis with deformation behavior, revealing the stepwise movement characteristics of critical states during landsliding. Landslide failure occurs when the final slice reaches the critical state. Utilizing the Daochi Village landslide in China as a case study, the unbalanced thrust slice method under two modeling frameworks reveals evolutionary characteristics of physical-mechanical parameters during progressive failure. Results indicate that both slice models can effectively describe the progressive failure process of the Daochi Village landslide. The stability coefficient derived from the traditional partial strength reduction perfect elasto-plastic model exceeds that of the complete process constitutive model.

1 Introduction

China is characterized by mountainous terrain and experiences substantial annual economic losses due to slope disasters. Geological disaster case studies show that factors such as earthquakes, rainfall, and excavation are the primary triggers for landslides and debris flows. Landslides typically undergo a process of crack initiation, local failure, and ultimately progressive failure leading to overall interconnection, with some evolving into debris flows. Therefore, it is crucial to study geological hazards and implement effective prediction and forecasting, with monitoring and early warning systems primarily based on deformation characteristics. Over the past century, more than ten types of limit equilibrium analysis methods have been developed, including Fellenius (1936), Bishop (1955), Spencer (1967), and Janbu (1973), the transfer coefficient method, Sarma and Tan (2006), and the finite element strength reduction method (FEM) (Kelesoglu, 2016; Nian et al., 2012; Zhu and Lee, 2002; Gharti et al., 2012; Tiwari et al., 2015; Wei et al., 2009). Limit equilibrium methods are widely used in engineering practice (Baker, 2004; Matthews et al., 2014; Faramarzi et al., 2017). However, for a given slip surface, different limit equilibrium methods make specific assumptions regarding the normal and tangential forces on the slice base, as well as the point and direction of thrust, leading to varying calculation results for the same slope. In two-dimensional (2D) landslide stability analysis, traditional slice methods typically divide the sliding mass into n slices, requiring the introduction of n−1 assumptions to make the problem statically determinate and solvable. These assumptions compromise theoretical rigor of slope stability analysis methods and rarely align with actual site conditions. Inappropriate method selection may produce erroneous results, potentially causing severe engineering consequences. These issues become even more complex when applying 2D methods to three-dimensional (3D) analysis, as the number of unknowns and required assumptions increases significantly, often rendering such methods impractical. Some studies have proposed a partial strength reduction method to simulate the progressive failure process of slopes (Lu et al., 2013; Lu et al., 2017; Lu, 2015; Lu et al., 2025).

The finite element method (FEM) is widely applied in numerical analysis in civil engineering. It was first proposed by Hrennikoff (1941) and has since been continuously refined and developed by numerous scholars (Oden et al., 1998; Wilson et al., 2013). Key code platforms in the field of slope engineering, such as ABAQUS, ANSYS, and FLAC, are primarily based on small deformation and continuum mechanics theories. In recent years, an improved finite element method has been proposed and applied in engineering practice (Mohammadi, 2008; Planas et al., 2003). Other computational methods, such as discrete element method (DEM) and discontinuous deformation analysis (DDA), are emerging and being actively employed in slope stability analysis. With the development of numerical analysis techniques, various computational methods have been introduced (Lu et al., 2013; Lu et al., 2017; Lu, 2015; Lu et al., 2025). More recently, a new numerical computation method has been proposed, which can account for the effects of various boundary conditions and possesses broad potential for extension and application (Zhang et al., 2023; Lu et al., 2023; Zhang et al., 2025).

With the improvement of understanding of natural materials, translational and retrogressive landslides can be classified along the slip surface into four zones: unstable, critical, metastable, and stable zone, i.e., the driving zone and the resisting zone. Within the driving zone, both the driving shear stress and the frictional resistance along the slip surface are discontinuous, and the corresponding displacements are also discontinuous. However, the normal stress remains continuous along the entire slip surface. This study characterizes critical state landslide behavior and establishes a slip surface boundary method for numerical analysis, resolving critical state identification and discontinuity challenges in stress strain solutions along failure surfaces. Based on the analysis of landslide deformation and failure mechanisms, five integrated approaches for describing the landslide failure process are defined: the comprehensive sliding-resisting method, the comprehensive displacement method, the main thrust method, the tensile failure method, and the surplus displacement method (Zhang et al., 2023; Lu et al., 2023; Zhang et al., 2025; Yingfa, 2016).

Traditional slice methods do not consider the progressive failure process of slopes, and the safety factor is only related to forces, not displacements. Therefore, it is necessary to extend these methods to characterize the progressive failure process and its close relation to slope deformation, enabling practical engineering applications. This study extends the traditional slice-based method for slope stability analysis by incorporating the perfect elasto-plastic model and the full-process constitutive model, thereby enabling the method to account for deformation effects. Comparative evaluation of both models provides a theoretical basis for the design of slope stabilization and treatment projects.

2 Overview of traditional unbalanced thrust methods

The traditional limit equilibrium slice method is widely applied in engineering practice, but it has some limitations, including simplifications of slices representation and force transmission compromising theoretical rigor, and neglect of deformation, which prevents the incorporation of field-measured displacements into the analysis.

Using the traditional unbalanced thrust method as an example, fundamental assumptions are:

2.1 Fundamental assumptions

1. Each slice is assumed to behave as a rigid body without deformation and is divided at fixed vertical intervals.

2. The force from the rear slice on the front slice is assumed to be parallel to the base of the rear slice, applied at the centroid of the front slice.

3. Rotational moments are neglected.

4. The base of each slice is in a critical state of shear stress.

Under these assumptions, slice division typically follows Figure 1. The fundamental equations of the unbalanced thrust method are given by Equations 1–9.

Figure 1

Figure 1. The scheme of slice block layout of slope stability analysis.

For the ith slice:

Normal force $N_i$:

$N_i=W_i\cos {\alpha }_i+P_{i-1}\sin \left({\alpha }_{i-1}-{\alpha }_i\right)+{\beta }_i{l}_i\cos {\alpha }_i\cos {\alpha }_i+{\Delta }_i{l}_i\cos {\alpha }_i\sin {\alpha }_i(1)$

Normal stress ${\sigma }_{i,n}$:

${\sigma }_{i,n}=N_i/l_i(2)$

Critical frictional stress ${\tau }_{i,peak}$:

${\tau }_{i,peak}=c_i+{\sigma }_{i,n}\mathit{tan}{\phi }_i(3)$

Frictional resistance $T_i^{crit}$:

$T_i^{crit}=c_i{l}_i+N_i\tan {\phi }_i(4)$

Reduced frictional stress for anti-sliding ${\tau }_{i,f}$:

${\tau }_{i,f}=\left(c_i+{\sigma }_{i,n}\mathit{tan}{\phi }_i\right)/f(5)$

Reduced frictional resistance after anti-sliding reduction $T_{i,F}^{crit}$:

$T_{i,F}^{crit}=T_i^{crit}/f(6)$

Driving force $P_i^S$:

$P_i^S=W_i\cos {\alpha }_i+P_{i-1}\cos \left({\alpha }_{i-1}-{\alpha }_i\right)+{\beta }_i{l}_i\cos {\alpha }_i\sin {\alpha }_i+{\Delta }_i{l}_i\cos {\alpha }_i\cos {\alpha }_i(7)$

Driving shear stress ${\tau }_i$:

${\tau }_i=P_i^S/l_i(8)$

Surplus driving force $P_i$:

$P_i=P_i^S-T_{i,F}^{crit}(9)$

Where:

$W_i$ is the weight of the ith slice; ${\beta }_i$ and ${\Delta }_i$ are vertical and horizontal uniform loads, respectively; $l_i$ is the base length of the slice; ${\alpha }_i$ is the angle between the base and the horizontal; $c_i$ is the cohesion; ${\phi }_i$ is the friction angle; $f$ is the strength reduction factor when the entire slip surface is at critical state; ${\sigma }_{i,n}$ is the normal stress on the base of the slice.

In engineering analysis, iterative computation of Equations 1–9 yields the reduction factor $f$ (i.e., stability coefficient). Notably, in this method, the frictional stress at the base of each slice is assumed to be at the ultimate limit state. However, a slip surface entirely under such critical stress conditions is rarely realized in natural conditions. The model becomes theoretically feasible only when the slip surface is in a residual stress state, allowing the physical and mechanical parameters of the base to be set to residual stress parameters. Nevertheless, after applying frictional resistance reduction, the calculated stress, except when $f=1$, logically cannot be compared with the actual stress state, even when the slip surface is in a residual stress state.

3 Partial strength reduction methods

For a general slope, the division of slices based on the unbalanced thrust method is illustrated in Figure 1. A common approach in this method is to consider all slices from 1 to n simultaneously in Equations 1–9 to derive the global stability coefficient of the slope. This coefficient essentially corresponds to the state where the driving force of the last slice equals its critical frictional resistance (i.e., the last slice attains critical state). Previous studies introduced the partial strength reduction method (Lu et al., 2013; Lu et al., 2017; Lu, 2015; Lu et al., 2025), which is explained as follows:

The implementation steps of the partial strength reduction method are as follows: Select slices 1 to m and compute via Equations 1–9, and when the calculated stability coefficient equals 1, the corresponding slice group is identified as being in the critical state (i.e., the mth slice is the critical slice). Then repeat the stability analysis progressively by expanding the range from slices 1 to m+1,. , up to 1 to n, thereby obtaining a series of partial strength reduction coefficients (stability coefficients, $f_i,i\in \left(m,n\right)$). The physical significance of this process is that the critical state point progressively moves forward along the slope from the mth slice to the (m+1)-th, . , until it reaches the nth slice. The results show that the stability coefficient is the lowest (equal to 1) when the mth slice is in the critical state, and it gradually increases as the critical state slice moves forward, reaching its maximum of $f_n$ when the nth slice is at the critical state. Evidently, the stability coefficient increases as the critical slice moves forward, which contradicts both field observations and rational analysis. By referencing concepts from the main thrust method and surplus displacement method proposed (Lu et al., 2013; Lu et al., 2017; Lu, 2015; Lu et al., 2025), the surplus stability coefficient ($f_{zs}^i$) for each critical slice slope can be computed. The surplus stability coefficient of the partial strength reduction method is defined as the difference between the global strength reduction stability coefficient of the slope ($f_n$) and the strength reduction stability coefficient of the ith critical slice ($f_i$). The surplus stability coefficient $f_{zs}^i$ of the slope corresponding to the ith critical state under the partial strength reduction method is defined as Equation 10:

$f_{zs}^i=f_n-f_i(10)$

Field application of this partial strength reduction surplus coefficient indicates that it is essential for researchers to determine the location of the critical slip surface in situ. This provides a theoretical basis for landslide predication, early warning, stabilization and treatment design.

4 Slice method based on the perfect elasto-plastic model

Employing the core concept of partial strength reduction, this section conducts progressive failure analysis using the perfect elasto-plastic model and proposes the following fundamental assumptions:

4.1 Fundamental assumptions

1. It is assumed that the slices (or the stress–strain relationship of geological materials) conform to the perfect elasto-plastic model, permit sufficient deformation, are divided vertically at fixed intervals, and that the peak shear strength is subject to strength reduction.

2. For slices within the failure zone, the shear strains of the ith and (i+1)-th slices satisfy the vector summation relationship, which involves shear strain components in the directions parallel and perpendicular to the base (Figure 2).

Traditional assumptions (2) and (3) from Section 1 remain operative. Based on these two fundamental assumptions, the slice division and force calculations consistent with the traditional method. In the failure zones, the shear strain of adjacent slices exhibits the following characteristics as shown in Equation 11.

${\gamma }_i^s={\gamma }_{i+1}^s+{\gamma }_{i+1}^{n}OR{\gamma }_{i+1}^s={\gamma }_i^s+{\gamma }_i^n(11)$

The physical implication of Assumption (6) is that, within the failure zone, there are two physical meanings:

First, the (i+1)-th slice advances, while the ith slice provides the deformation condition, as described in Equation 12:

${\gamma }_{i+1}^s={\gamma }_i^s\cos \left({\alpha }_i-{\alpha }_{i+1}\right)(12)$

Second, the ith slice advances, while the (i+1)-th slice provides the deformation condition, as described in Equation 13:

${\gamma }_i^s={\gamma }_{i+1}^s\cos \left({\alpha }_i-{\alpha }_{i+1}\right)(13)$

When the angles between the sliding direction and the horizontal are equal for both slices, their shear strains become consistent. Under these assumed conditions, the slice division is illustrated in Figure 1. The fundamental formulas for the partial strength reduction unbalanced thrust method are presented in Equations 1–9.

Figure 2

Figure 2. Shear strain and vector relationship between adjacent slices. (a) Shear strain between adjacent slices. (b) The ith block advances, and the conditions provided by the i+1-th block. (c) The i+1-th block advances, and the conditions provided by the ith block.

4.2 Perfect elasto-plastic model

In the progressive failure analysis based on the perfect elasto-plastic model, the partial strength reduction method is employed to reduce strength from a stress perspective. This reduction leads to an increase in surplus thrust within the post-failure zone, thereby forcing the critical stress state to gradually advance. The perfect elasto-plastic model is described as follows:

One of the traditional shear stress models is the perfect elasto-plastic model. In this model, shear stress is linearly related to shear strain before reaching the critical shear stress (which satisfies the Mohr-Coulomb criterion or other criteria), as shown in Figure 3.

${\tau }_i=G_i{\gamma }_i,when{\gamma }_i\le {\gamma }_{i,peak}(14)$

Figure 3

Figure 3. Perfect elasto-plastic model.

Where $G_i$ represents the shear modulus of the base of the ith slice (or element). From Equation 14, it is evident that when the shear stress reaches the critical stress state, the stress no longer uniquely determines the strain, indicating that the relationship between stress and strain becomes non-unique.

The traditional perfect elasto-plastic model presents several limitations in the context of slope stability analysis: (1) When the applied shear stress exceeds the critical shear stress, the shear stress can only be set to the critical shear stress and remains constant. As a result, the relationship between shear stress and shear strain becomes indeterminate, making it difficult to quantify the corresponding shear strain; (2) After the shear stress reaches the critical shear stress, it remains constant. This model cannot reflect the decreasing characteristic of shear stress (i.e., shear stress softening) when shear strain exceeds the critical shear strain. (3) Using the perfect elasto-plastic model, if a landslide calculation shows that the mth slice is in a critical state (Figure 1), the stability coefficient for slices 1 to m, calculated by the traditional frictional resistance reduction method, would be 1. However, as slip surface deformation increases, assuming the (m+1)-th slice’s deformation reaches a critical state, if the perfect elasto-plastic model is used without changing slice shape or mechanical parameters, the resisting force will increase (note: the resisting force of the (m+1)-th slice increases). Consequently, the calculated stability coefficient for slices 1 to (m+1) will be greater than that calculated for slices 1 to m, implying an increase in slope stability (this conclusion often contradicts actual engineering practice). These disadvantages limit the traditional model to reflect the progressive deformation and failure process of slopes, particularly the commonly recognized characteristics of decreasing resisting force and reduced stability in previously failed slices as displacement increases. To overcome these limitations, the partial strength reduction method is employed to reduce the frictional resistance of geological materials, thereby reflecting the softening characteristics of frictional resistance.

4.3 Determination of the critical state slice

For the perfect elasto-plastic model, the extended slice force-displacement method is used to calculate the frictional resistance, driving force, and surplus driving force along the slip surface. If, during the calculation, the surplus driving force for the (m-1)-th slice is $P_{m-1}>0$, but for the mth slice, the surplus driving force is $P_m<0$, then the critical state slice must be located within the mth slice. Since the base length of the mth slice can be relatively large, leading to $P_m<0$, it becomes necessary to subdivide the base length of the mth slice further. As slices in the slice force-displacement method are generally triangular, parallelogram, or trapezoidal, and a trapezoid can represent all three shapes, the following analysis focuses on the trapezoidal shape.

Assume that the four corner points of the mth trapezoidal slice are denoted as points $1\left(x_1^m,y_1^m\right)$, $2\left(x_2^m,y_2^m\right)$, 3 $\left(x_3^m,y_3^m\right)$, and 4 $\left(x_4^m,y_4^m\right)$, with an inclination angle of ${\alpha }_m$, and that the surplus driving force of (m−1)-th slice, $P_{m-1}>0$, and these parameters are all known. The mth slice is then decomposed into the mth and (m+1) slices, where the new mth slice is defined by points 1, 2, 5, and 6 (Figure 4). The objective is to determine the coordinates of points 5 and 6.

Figure 4

Figure 4. Scheme of slice decomposition.

The line formed by points 1 and 4 is given by Equation 15:

$y_{14}=k_{14}{x}_{14}+b_{14}(15)$

Where $k_{14}=\frac{y_4^m-y_1^m}{x_4^m-x_1^m}$ and $b_{14}=\frac{y_4^m-y_1^m}{x_4^m-x_1^m}{x}_1^m$.

Similarly, the equation of the line formed by points 2 and 3, The coordinates of points 5 and 6 are then $\left(x_5^m,k_{14}{x}_5^m+b_{14}\right)$ and $\left(x_5^m,k_{23}{x}_5^m+b_{23}\right)$, respectively. Thus, the only unknown is $x_5^m$. the weight $W_m^N$ and base length $l_m^N$ of trapezoid 1-2-5-6 can be calculated. The force expression for the redefined mth slice is then obtained by calculating the normal pressure $N_m^N$, Critical frictional resistance $T_m^{crit,N}$, and driving force $P_m^{s,N}$ using Equations 1, 4, 7, respectively.

According to the definition of the critical state slice, the driving force is set equal to the shear resistance, as expressed in Equation 16.

$T_m^{crit,N}=P_m^{s,N}(16)$

At this point, the resulting equation becomes a quadratic equation in a single variable, which can be solved for $x_5^m$. A solution satisfies the critical state condition if $x_5^m$ falls within the range [$x_1^m$, $x_3^m$]. When the slice is triangular, $y_3^m=y_4^m$; when the slice is a parallelogram, $\left|y_2^m-y_1^m\right|=\left|y_3^m-y_4^m\right|$.

At this stage, the frictional stress along the base of slices 1 to m is the stress after strength reduction (calculated using Equation 5). The normal stress is calculated using Equation 2, and the driving shear stress is calculated using Equation 8. Slices from m+1 to n are in a pre-peak shear stress state. Here, the frictional stress is numerically equal to the driving force divided by the base length, the normal stress is calculated using Equation 2, and the driving shear stress is calculated using Equation 8. Thus, the current frictional shear stress, driving shear stress, and normal stress for all slices from 1 to n are fully determined, and denoted as ${\tau }_i^X$, ${\tau }_i^u$, ${\sigma }_{i,n}^X$, respectively, $i\in \left(1,n\right)$.

4.4 Determination of slip surface displacement for the perfect elasto-plastic model

For the perfect elasto-plastic model, once the critical stress state slices have been determined as described above, the stress distribution along the slip surface is also established. The determination of the current shear strain for slices 1 to n using the slice method is as follows:

For slices from m+1 to n, which are in a linear shear stress state, the shear strain is calculated as shown in Equation 17.

${\gamma }_i=P_i^s/\left(l_i{G}_i\right),iϵ\left(m+1,n\right)(17)$

For the mth slice, the shear strain is given by Equation 18:

${\gamma }_{m,peak}=T_m^{peak}/\left(l_m{G}_m\right)(18)$

Regarding the treatment of deformation in the failure zone, Assumption (6) is proposed. Based on this assumption, the specific strain solutions are as follows:

For the failure zone (including slices or elements at the critical state), the relationship between the base shear strains of the ith and (i+1)-th slices can be categorized into two cases:

Case 1. The advancing strain of ith slice is the contribution strain (${\gamma }_i^s$) from the failed (i+1)-th slice to the ith slice (see Equations 12, 13).

${\gamma }_i={\gamma }_i^s$

However, if the contribution strain (${\gamma }_i^s$) is less than the critical strain (${\gamma }_{i,peak}$), then the critical strain (${\gamma }_{i,peak}$) is adopted, as specified in Equation 19.

Case 2. The shear strain of the ith slice is the sum of the critical shear strain before reduction (${\gamma }_{i,peak}$), the strain induced by strength reduction (${\gamma }_{i,peak}^{*}$), and the contribution strain from the failed (i+1)-th slice to the ith slice (${\gamma }_i^s$), as formulated in Equation 20:

${\gamma }_i={\gamma }_{i,peak}+{\gamma }_i^s-{\gamma }_{i,peak}^{*}$

The frictional stress (${\tau }_i^z$) after reducing the resisting stress is given in Equations 21, 22:

${\tau }_{i,f}={\tau }_{i,peak}/f_i$

${\gamma }_{i,peak}^{*}={\tau }_{i,f}/G_i$

For critical state slices (or elements) as shown in Equation 23.

${\gamma }_i={\gamma }_{i.,peak},{\gamma }_{i,peak}={\tau }_{i,peak}/G_i$

Under these two assumed conditions, the non-unique relationship between stress and strain in the peak stress and post-failure zones is fixed, which facilitates subsequent stability evaluation.

Based on Sections 3.1–3.4, solutions for the strain problem can be obtained under the assumptions. Utilizing these solutions, the stability evaluation of landslide progressive failure can be performed. The calculated stresses and displacements are the current shear stresses and displacements. Furthermore, these results can be directly compared with in situ monitoring data of stress and displacement, allowing for calibration and refinement of corresponding model parameters.

4.5 Progressive failure analysis using partial strength reduction

The objective of failure analysis is to determine the final failure state of the landslide. With a known failure path, the partial strength reduction method, based on the perfect elasto-plastic model, is applied to progressively reduce the strength of slices. This ensures that each slice undergoes a critical stress state, until the last slice (or element) reaches its critical condition, obtaining the stress-strain distribution characteristics along the slip surface. Consequently, at the point of failure, the frictional shear stress, driving shear stress, normal stress, and shear strain along the slip surface are fully determined, denoted as: ${\tau }_i^{p,b}$, ${\tau }_i^{u,p,b}$, ${\sigma }_{i,n}^{p,b}$, and ${\gamma }_i^{p,b}$, respectively, $i\in \left(1,n\right)$.

5 Slice method based on the shear stress-strain complete process constitutive model

The slice method based on the shear stress-strain complete process constitutive model aims to overcome several limitations of the traditional slice method. For instance: 1) It does not reduce the resisting force, which leads to a more realistic determination of forces. 2) The novel complete process constitutive theory can reflect the characteristic strength degradation of materials as displacement increases, and it can also incorporate field-measured deformations into the stability analysis. Taking the surplus thrust method as an example, the slice method based on the stress-strain complete process constitutive model is described as follows:

5.1 Fundamental assumptions

1. It is assumed that the slices possess sufficient deformation capacity, and that they are vertically divided at fixed intervals.

2. Resisting stress at the base of each slice satisfies the complete process constitutive equation.

Assumptions (2), (3), and (6) remain valid. Based on these assumptions, the calculation procedure using the shear stress-strain complete process constitutive model proceeds as follows:

The normal force $N_i$, normal stress ${\sigma }_{i,n}$, driving force $P_i^S$, driving shear stress ${\tau }_i$, and surplus driving force $P_i$ for the ith slice are calculated using Equations 1, 2, 7–9. However, the shear stress ${\tau }_i$ and shear resistance $T_i$ along the base are calculated using Equations 24, 25 as follows:

${\tau }_i=G_i{\gamma }_i{\left[1+{\gamma }_i^{q_i}/P_i\right]}^{{\zeta }_i}(24)$

$T_i=G_i{\gamma }_i{\left[1+{\gamma }_i^{q_i}/P_i\right]}^{{\zeta }_i}{l}_i(25)$

5.2 Determination of critical state for the shear stress-strain complete process constitutive model (D-CPCM)

The shear stress-strain complete process constitutive model can describe the softening mechanical characteristics of the post-failure zone. Under the condition that the critical frictional stresses all satisfy a strength criterion (e.g., Mohr-Coulomb criterion (MC) or other criteria), the critical state slice corresponding to the shear stress-strain complete process constitutive model will occur after the mth slice. The calculation steps are as follows:

Based on the calculations from the perfect elasto-plastic model (with a reduction coefficient equal to 1), the frictional shear stress, driving shear stress, and normal stress for slices 1 to m are obtained and denoted as ${\tau }_i^X$, ${\tau }_i$, ${\sigma }_{i,n}^X$, and ${\gamma }_i$, respectively, $i\in \left(1,m\right)$. For slices 1 to m, calculations based on the shear stress-strain complete process constitutive model are performed as follows: For any ith slice, the normal force $N_i$, normal stress ${\sigma }_{i,n}$, critical shear stress ${\tau }_{i,peak}$, and critical shear resistance $T_i^{crit}$ are calculated using Equations 1–4, respectively. The corresponding critical strain is calculated by Equation 26:

${\gamma }_{i,peak}=a_1^0+a_2^0{\sigma }_{i,n}(26)$

Where $a_0^1$ is dimensionless and $a_0^2$ has dimension $1/kPa$.

For the dimensionless softening parameter $\xi$, it can be expressed as shown in Equation 27:

${\xi }_i={\xi }_{i,0}/\left(1+{\left({\xi }_{i,0}/{\xi }_{i,c}-1\right)\left({\sigma }_{i,n}/{\sigma }_{i,n}^c\right)}^{{\zeta }_i}\right)(27)$

Where ${\xi }_{i,0}$ is the value of ${\xi }_i$ when the normal stress ${\sigma }_{i,n}=0$ , and ${\xi }_{i,0}=-0.9999$. ${\xi }_{i,c}$ is the value of ${\xi }_i$ when ${\sigma }_{i,n}$ equals to ${\sigma }_{i,n}^c$, and ${\zeta }_i$ is a constant coefficient. This relationship can be obtained from shear test curves under different normal pressures.

The model parameter $q_i$ is given in Equation 28:

$q_i=\frac{{\left[{\tau }_{i,peak}/\left({\gamma }_{i,peak}{G}_i\right)\right]}^{1/{\xi }_i}}{{\xi }_i\left[1-{\left({\tau }_{i,peak}/\left({\gamma }_{i,peak}{G}_i\right)\right)}^{1/{\xi }_i}\right]}(28)$

The model parameter $p_i$ can be calculated as shown in Equation 29:

$p_i={-\left({\gamma }_{i,peak}\right)}^{q_i}\left(1+q_i{\xi }_i\right)(29)$

Furthermore, the frictional stress ${\tau }_i$, frictional resistance $T_i$, driving force $P_i^S$, driving shear stress ${\tau }_i^u$, and surplus driving force $P_i$ are calculated using Equations 7–9, respectively. According to assumption (6), the critical shear strain for each slice from 1 to m is calculated, completing the first step of strain evaluation for slices 1 to m. The calculation is then repeated for slices 1 to m using the equations described above. When the calculation reaches the mth slice, its surplus driving force is $P_m>0$. This implies that, under the new model conditions, the mth slice is not yet in a critical state. Therefore, calculations must continue for slice m+1, with the following steps: 1) Calculate the critical strain using Equations 26, 2) Re-calculate the shear strain for slices from m+1 back to slice 1, based on assumption (6).

If the surplus driving force of the (L−1)-th slice is greater than zero while that of the Lth slice is less than zero (as shown in Figure 4), the Lth slice must be subdivided according to the method described in Section 3.3, and the above steps are repeated to determine the critical slice for the new shear stress constitutive model (e.g., the Lth slice).

According to progressive failure analysis, assuming the Lth slice develops into a critical stress state, the surplus thrust from the (L-1)-th slice to Lth slice can be calculated. The length of the Lth slice is then calculated as follows:

According to Section 4.2, the surplus driving force of the (L−1)-th slice $P_{L-1}$ can be calculated. Setting the surplus thrust of the Lth slice $P_L=0$, the $P_{L-1}$ can be calculated as shown in Equations 30, 31:

$P_{L-1}=\frac{X_x}{X_y}(30)$

Where:

$\left.\begin{array}{l}{X}_x=W_L\left(\mathit{sin}{\alpha }_L-\cos {\alpha }_L\tan {\phi }_L\right)+\\ {\beta }_L{l}_L\cos {\alpha }_L\left(\sin {\alpha }_L-\cos {\alpha }_L\tan {\phi }_L\right)+\\ {\Delta }_L{l}_L\cos {\alpha }_L\left(\cos {\alpha }_L-\cos {\alpha }_L\tan {\phi }_L\right)-c_L{l}_L\\ X_y=\mathit{sin}\left({\alpha }_{L-1}-{\alpha }_L\right)\tan {\phi }_L-\cos \left({\alpha }_{L-1}-{\alpha }_L\right)\end{array}\right\}(31)$

Thus, the base length $l_L$ of the critical slice L can be determined based on the surplus thrust received from slice (L−1).

At this point, the frictional stress along the base of slices 1 to L is calculated using Equation 24, the normal stress using Equation 2, and the driving shear stress using Equation 8. Slices from L+1 to n are in a pre-peak shear stress state. Here, the frictional stress is numerically equal to the driving force divided by the base length, the normal stress and the driving shear force are calculated using Equations 2, 8. Thus, the current frictional shear stress, driving shear stress, and normal stress for all slices from 1 to n are fully determined, and denoted as ${\tau }_i^X$, ${\tau }_i$, and ${\sigma }_{i,n}^X$, respectively, $i\in \left(1,n\right)$.

6 Stability analysis

For a specific slope, the global traditional stability coefficient of a landslide can be calculated using Equations 1–9. To evaluate the landslide characteristics at different stages and assess its stability, the current frictional shear stress, driving shear stress, normal stress, and shear strain are obtained using both the perfect elasto-plastic model and the shear stress-strain complete process constitutive model, denoted as ${\tau }_i^X$, ${\tau }_i$, ${\sigma }_{i,n}^X$, and ${\gamma }_i$, respectively, $i\in \left(1,n\right)$. Under failure conditions, the corresponding parameters are ${\tau }_i^{p,b}$, ${\tau }_i^{u,p,b}$, ${\sigma }_{i,n}^{p,b}$, and ${\gamma }_i^{p,b}$, respectively, $i\in \left(1,n\right)$. Under two-dimensional conditions, the overall stability evaluation of the sliding mass is expressed as follows:

6.1 Stability coefficient based on the comprehensive sliding–resisting method ($F_{CSRM}$)

The horizontal and vertical components, as well as the comprehensive vector sum, of the current driving force $P_x^S$ for each slice are given as follows, as shown in Equations 32–34:

$P_x^S={\displaystyle \sum _{i=1}^n}{\tau }_i{l}_i\cos {\alpha }_i(32)$

$P_y^S={\displaystyle \sum _{i=1}^n}{\tau }_i{l}_i\sin {\alpha }_i(33)$

$P^S=\sqrt{{\left(P_x^S\right)}^2+{\left(P_y^S\right)}^2}(34)$

The vector sum is denoted as $P^S$, and the minimum angle ${\alpha }_s$ it forms with the horizontal axis is given by Equation 25:

${\alpha }_s=\arctan \left(P_y^S/P_x^S\right)(35)$

Resistant shear force: For each slice along the slip surface in a failure state, the components of the frictional resistance and normal force in the horizontal and vertical directions, as well as their comprehensive vector sum, are given as follows, as shown in Equations 36–38:

$T^x={\displaystyle \sum _{i=1}^n}{\tau }_i^{p,b}{l}_i\cos {\alpha }_i+{\displaystyle \sum _{i=1}^n}{\sigma }_{i,n}^{p,b}{l}_i\cos {\alpha }_i(36)$

$T^y={\displaystyle \sum _{i=1}^n}{\tau }_i^{p,b}{l}_i\sin {\alpha }_i+{\displaystyle \sum _{i=1}^n}{\sigma }_{i,n}^{p,b}{l}_i\sin {\alpha }_i(37)$

$T=\sqrt{{\left(T^x\right)}^2+{\left(T^y\right)}^2}(38)$

The comprehensive vector sum $T$ forms a minimum angle α with the horizontal axis, given by Equation 39:

${\alpha }_f=\arctan \left(T^y/T^x\right)(39)$

The stability coefficient in the horizontal direction is calculated by Equation 40:

$F_{CSRM}^x=\left|T^x/P_x^S\right|(40)$

The stability coefficient in the vertical direction is described in Equation 41:

$F_{CSRM}^y=\left|T^y/P_y^S\right|(41)$

The stability coefficient along the direction of the driving shear force is given in Equation 42:

$F_{CSRM}^s=\left|T\cos \left({\alpha }_f-{\alpha }_s\right)/P^S\right|(42)$

6.2 Stability coefficient based on the main thrust method ($F_{MTM}$)

For slices ranging from the first slice to the critical state slice (i.e., the mth slice), the vector sum of the surplus driving force $P_i$ in the horizontal and vertical directions is calculated. The horizontal and vertical vector sums obtained are $P_{xm}$ and $P_{ym}$, respectively. Their comprehensive vector sum is $P^m$, and the minimum angle it forms with the horizontal axis is ${\alpha }^m$. They can be calculated as shown in Equations 43–46:

$P_{xm}={\displaystyle \sum _{i=1}^m}{P}_i\cos {\alpha }_i(43)$

$P_{ym}={\displaystyle \sum _{i=1}^m}{P}_i\sin {\alpha }_i(44)$

$P^m=\sqrt{{\left(P_{xm}\right)}^2+{\left(P_{ym}\right)}^2}(45)$

${\alpha }^m=\arctan \left(P_{ym}/P_{xm}\right)(46)$

For slices from m+1 to n ($i\in \left(m+1,n\right)$), the difference between the frictional stress at failure (${\tau }_i^{p,b}$) and the current frictional stress (${\tau }_i^X$) for each slice along the slip surface, in terms of their horizontal and vertical components and their resultant vector sum, are respectively shown in Equations 47–49:

$T^{xm}={\displaystyle \sum _{i=m+1}^n}\left({\tau }_i^{p,b}-{\tau }_i^X\right)l_i\cos {\alpha }_i(47)$

$T^{ym}={\displaystyle \sum _{i=m+1}^n}\left({\tau }_i^{p,b}-{\tau }_i^X\right)l_i\sin {\alpha }_i(48)$

$T^m=\sqrt{{\left(T^{xm}\right)}^2+{\left(T^{ym}\right)}^2}(49)$

The resultant vector $T^m$ forms a minimum angle ${\alpha }_f^m$ with the horizontal axis, given by Equation 50:

${\alpha }_f^m=\arctan \left(T^{ym}/T^{xm}\right)(50)$

The main thrust stability coefficient in the horizontal direction is given in Equation 51:

$F_{MTM}^x=\left|T^{xm}/P_{xm}\right|(51)$

The main thrust stability coefficient in the vertical direction is given in Equation 52:

$F_{MTM}^y=\left|T^{ym}/P_{ym}\right|(52)$

The stability coefficient along the direction of the main thrust is shown in Equation 53:

$F_{MTM}^s=\left|T^m\cos \left({\alpha }_f^m-{\alpha }^m\right)/P^m\right|(53)$

6.3 Stability coefficient based on the comprehensive displacement method ($F_{CDM}$) (Note: displacement caused by normal stress is not considered in this method)

The current shear strain ${\gamma }_i$ ($i\in \left(1,n\right)$) of each slice is resolved into horizontal and vertical components, and the comprehensive resultant vector are respectively shown in Equations 54–56:

$S_x^s={\displaystyle \sum _{i=1}^n}{\gamma }_i{l}_i\cos {\alpha }_i(54)$

$S_y^s={\displaystyle \sum _{i=1}^n}{\gamma }_i{l}_i\sin {\alpha }_i(55)$

$S^s=\sqrt{{\left(S_x^s\right)}^2+{\left(S_y^s\right)}^2}(56)$

The resultant vector $S^s$ forms a minimum angle ${\alpha }_s^s$ with the horizontal axis, given by Equation 57:

${\alpha }_s^s=\arctan \left(S_y^s/S_x^s\right)(57)$

Displacement at Failure: For each slice along the slip surface, its displacement components in the horizontal and vertical directions, and their resultant vector sum (Note: normal displacement is ignored here) are respectively shown in Equations 58–60:

$S_x={\displaystyle \sum _{i=1}^n}{\gamma }_i^{p,b}{l}_i\cos {\alpha }_i(58)$

$S_y={\displaystyle \sum _{i=1}^n}{\gamma }_i^{p,b}{l}_i\sin {\alpha }_i(59)$

$S=\sqrt{{\left(S_x\right)}^2+{\left(S_y\right)}^2}(60)$

The resultant vector $S$ forms a minimum angle ${\alpha }_s^{crit}$ with the horizontal axis, given by Equation 61:

${\alpha }_s^{crit}=\arctan \left(S_y/S_x\right)(61)$

The stability coefficient in the horizontal displacement direction is described as Equation 62:

$F_{CDM}^x=\left|S_x/S_y\right|(62)$

The stability coefficient in the vertical displacement direction is given by Equation 63:

$F_{CDM}^y=\left|S_x/S_y^s\right|(63)$

The stability coefficient along the sliding displacement direction is given by Equation 64:

$F_{CDM}^s=\left|Scos\left({\alpha }_s^{crit}-{\alpha }_s^s\right)/S^s\right|(64)$

6.4 Stability coefficient based on the surplus displacement method ($F_{SDM}$)

Similar to the definition of the surplus force stability coefficient, a surplus displacement stability coefficient can be obtained. For a selected sliding mass, within the range of slices 1 to n, the displacement $S_i$ of each element along the slip surface at failure is calculated. Similarly, within the range of slices from m+1 to n, when the sliding mass undergoes global failure along the potential slip surface, the difference between the failure displacement $S_{crit}^i$ and the current displacement $S^i$ for each slice under the possible failure mode along the slip surface is calculated. The vector sums of these displacement differences are then determined in both the horizontal and vertical directions, yielding the horizontal and vertical displacement difference vector sums as $S_{crit}^x$ and $S_{crit}^y$, respectively, as shown in Equations 65–68. Their resultant vector sum is $S_{crit}$, and the minimum angle it forms with the horizontal axis is ${\alpha }_s^{crit}$.

$S_{crit}^x={\displaystyle \sum _{i=m+1}^n}\left({\gamma }_i^{p,b}-{\gamma }_i^u\right)l_i\cos {\alpha }_i(65)$

$S_{crit}^y={\displaystyle \sum _{i=m+1}^n}\left({\gamma }_i^{p,b}-{\gamma }_i^u\right)l_i\sin {\alpha }_i(66)$

$S_{crit}=\sqrt{{\left(S_{crit}^x\right)}^2+{\left(S_{crit}^y\right)}^2}(67)$

${\alpha }_s^{crit}=\arctan \left(S_{crit}^y/S_{crit}^x\right)(68)$

Surplus coefficient in the horizontal displacement direction is given by Equation 69:

$F_{SDM}^x=\left|S_{crit}^x/S_m^x\right|(69)$

Surplus coefficient in the vertical displacement direction is given by Equation 70:

$F_{SDM}^y=\left|S_{crit}^y/S_m^y\right|(70)$

Surplus coefficient along the main sliding displacement direction is calculated as Equation 71.

$F_{SDM}^s=\left|S_{crit}\cos \left({\alpha }_s^{crit}-{\alpha }_{sm}^x\right)/S_m\right|(71)$

Note: In the above equations, m refers to the critical slice (mth slice). For the multi-parameter stability description of progressive failure based on the shear stress-strain complete process constitutive model, the relationships are consistent with those presented above, and the critical stress state described above varies with displacement.

From the above theory, the object of this study is geomaterials, hose failure mode is compressive–shear, making it suitable for describing slope failure processes without significant rotational movement. To perform the evaluation of slope stability in progressive failure process, an algorithm was shown in Figure 5.

Figure 5

Figure 5. Algorithm of comparative analysis of two models of progressive failure of slope.

7 Case study

On 4 July 2016, the Huangyantang Landslide in Daochi Village, Hongtu Township, Enshi City was triggered by continuous heavy rainfall, resulting in the collapse of rock and soil mass. The collapsed slide mass had an average width of 240 m, a length of 200 m, and a volume of approximately 720,000 m³. Numerous surface cracks appeared on the slope, with their orientations generally parallel to the landslide direction, mostly being tensile and shear-offset cracks. The tensile cracks typically had openings of 50–100 mm, with a maximum opening of 300 mm and detectable depths of 200–400 mm; the shear-offset cracks generally had vertical displacements of 30–200 mm, causing cracks in the buildings above. The lower section of the old retaining wall fractured, while the unfractured sections developed cracks that penetrated the entire retaining wall. The rear edge of the landslide was closely connected to Maoping Highway, a two-lane road that serves as the main route connecting Enshi City to Jianshi and Hefeng Counties, posing a threat to 3 households near the rear edge, 10 households at the front edge, and Maoping Highway.

7.1 Geological overview

The Huangyantang Landslide is in the tectonic denudation hilly area of southwestern Hubei Province. The bedrock strikes northeast, forming a monoclinal landslide with a northwest dip direction. The cross-section exhibits a curvilinear morphology with slope gradients ranging 30°–38°. The slope surface is densely vegetated, and at the toe of the slope lies the Huangyantang pond. This rock-soil slope has surface elevations between 780.30 and 816.11 m and total vertical relief of 32.35–35.81 m. A platform area is present on the slope: it is designated for construction use and measures approximately 400 m in length and 30 m in width, with several existing buildings, showing approximately rectangular planar shape. The plan-view shape of the platform is roughly rectangular. The terrain is gentle with slope gradient less than 5° and ground elevation of 780.30–782.00 m, having relative relief of 1.70 m (Figure 6).

Figure 6

Figure 6. Pane view of Huangyantang Landslide.

According to field geological surveys and drilling data, the exposed strata within the investigation area primarily include miscellaneous fill soil (Q4m1) and the Lower Triassic Daye Formation (T1d). These strata, described from the youngest to the oldest, are as follows:

1. Miscellaneous fill soil (Q4m1): The lithology appears yellowish-brown, loose, and slightly moisture. The material composition consists mainly of silty clay mixed with fragments of yellowish-brown silty mudstone, mainly composed of silty clay with yellowish-brown argillaceous siltstone fragments. It is primarily distributed on the slope and lower platform of the Huangyantang landslide, with a general thickness of 0.40–12.20 m.

2. Strongly weathered limestone of the Lower Triassic Daye Formation (T1d): The rock stratum attitude is 290°–310°∠25°–35°. The lithology is yellowish-brown with argillaceous texture (high clay content) and soft quality. This layer is highly sensitive to water, showing characteristics of swelling, shrinkage, and large deformation, and it has poor resistance to weathering. Upon water exposure, it easily softens; when dehydrated, it readily cracks. After modification by shallow surface processes, the shallow rock mass becomes severely weathered and fragmented, with highly developed weathering fractures, poor rock mass integrity, low strength, and water-saturated soil properties, classifying it as a highly permeable layer. Dominated by strongly weathered limestone, it appears angular to subangular and flaky, in strongly to completely weathered states, forming weathered surplus blocks. This layer is mainly distributed around the Huangyantang Landslide area, relatively uniform in distribution and medium-thick layered, with Drilling revealing a general thickness of 6.90–8.80 m.

3. Moderately weathered limestone of the Lower Triassic Daye Formation (T1d): he moderately weathered limestone is gray to light gray. Exploration indicates this layer has a thickness of 7.60–31.40 m with relatively uniform distribution. The rock mass is intact, and according to existing geological data, its standard value of laboratory saturated uniaxial compressive strength is 73.92 MPa.

7.2 Hydrogeological conditions and deformation characteristics of the landslide

The landslide area exhibits underdeveloped groundwater systems, however, during the flood season, the region experiences abundant rainfall. The loosely structured miscellaneous fill distributed across the slope is prone to water absorption and saturation, forming pore water. Influenced by topography and soil-rock structure, this pore water generally migrates along the soil-rock interface. Long-term wet-dry cycling at these interfaces leads to progressive strength reduction, thereby increasing the possibility of slope failure and collapse of the overlying soil mass.

The underlying bedrock consists of the Lower Triassic Daye Formation (T1d), primarily composed of strongly weathered limestone. The lithology is soft and poorly resistant to weathering. This stratum primarily receives pore water recharge through shallow weathering fractures, forming bedrock fissure water that migrates along joint and fracture planes. Under long-term groundwater erosion, the slope is further destabilized, making landslide reactivation more likely.

7.3 Numerical calculation and analysis

Based on the I-I′ cross-section shown in the plane view (Figure 6), the slice-based calculation diagram is constructed as shown in Figure 7. The unit weight of the landslide mass is taken as 18 kN/m³. The slice division, as well as the base angles and lengths of each slice, are also indicated in Figure 7.

Figure 7

Figure 7. Slice division diagram of Huangyantang Landslide.

Based on laboratory testing and field experience, the model parameters are assigned as follows:

Cohesion $C$ = 24 kPa, friction angle $\phi$ = 24°, shear modulus $G$ = 3,000 kPa, ${\rho }_{i,c}$ = −2.96, ${\sigma }_i^{n,c}$ = 600 kPa, ${\zeta }_i$ = 0.3, $b_1$ = 50, and $b_2$ = 0 kPa^-1.

To ensure equivalence between the perfect elasto-plastic model and full-range shear stress-strain model at peak stress states, the critical shear strain is defined by linear formulas: ${\alpha }_1^0=c_i/G_i$, ${\alpha }_2^0=\tan {\phi }_i/G_i$.

For the Huangyantang landslide, the traditional critical method (TCM) was applied. Under saturated conditions, the calculated stability coefficient is 1.50025, indicating that the slope is in a stable state. At this stability coefficient, the relationships between the driving force ($P_i^{STCM}$), frictional resistance ($T_i^{TCM}$), and surplus driving force ($P_i^{TCM}$) for each slice and the slice number (SN) are illustrated in Figure 8. When the stability coefficient is 1.50025, the maximum driving force and maximum frictional resistance both occur at slice 11, after which they gradually decrease. By the final slice, the surplus frictional resistance reduces to zero.

Figure 8

Figure 8. Forces under the traditional critical method.

According to the partial strength reduction method (PSRM), when the stability coefficient equals 1.0, the critical slice is identified as slice 10. As the critical state progressively advances, the corresponding partial strength reduction stability coefficient gradually increases. When the last slice reaches the critical state, the partial strength reduction stability coefficient equals to the traditional strength reduction stability coefficient.

The relationship between the partial strength reduction stability coefficient and the critical slice number (CSN) is illustrated in Figure 9. The surplus stability coefficients ($f_{zs}^l$) corresponding to different critical states under the partial strength reduction method are shown in Figure 10. Figures 9, 10 demonstrate that as the critical state progressively advances forward, the overall stability coefficient of the landslide increases, while the surplus stability coefficient decreases. This indicates that the progressive forward movement of the critical stress state corresponds to a gradual decline in the overall stability of the slope.

Figure 9

Figure 9. The partial strength reduction stability coefficient.

Figure 10

Figure 10. The partial strength reduction surplus stability coefficient.

Based on the partial strength reduction method using the perfect elasto-plastic model (PEPM), a progressive failure analysis was conducted. During this process, five selected critical state slices were analyzed, and their corresponding driving forces, shear resistances, and surplus driving forces are shown in Figures 11–15, respectively. As the critical state slices progressively advances, the variations of the four slope stability coefficients ($F_{CRSM}$, $F_{MTM}$, $F_{CDM}$, and $F_{SDM}$) are illustrated in Figures 16–19.

Figure 11

Figure 11. Force distribution characteristics along slip surface for all slices when Slice 12 reaches critical state in perfect elasto-plastic model.

Figure 12

Figure 12. Force distribution characteristics along slip surface for all slices when Slice 15 reaches critical state in perfect elasto-plastic model.

Figure 13

Figure 13. Force distribution characteristics along slip surface for all slices when Slice 18 reaches critical state in perfect elasto-plastic model.

Figure 14

Figure 14. Force distribution characteristics along slip surface for all slices when Slice 21 reaches critical state in perfect elasto-plastic model.

Figure 15

Figure 15. Force distribution characteristics along slip surface for all slices when Slice 24 reaches critical state in perfect elasto-plastic model.

Figure 16

Figure 16. Stability coefficient based on CRSM using the perfect elasto-plastic model.

Figure 17

Figure 17. Stability coefficient based on MTM using the perfect elasto-plastic model.

Figure 18

Figure 18. Stability coefficient based on CDM using the perfect elasto-plastic model.

Figure 19

Figure 19. Stability coefficient based on SDM using the perfect elasto-plastic model.

As shown in Figures 11–15, with the progressive advancement of the critical state, the locations of maximum driving force, frictional resistance, and surplus thrust evolve respectively. Specifically, the slice with the maximum driving force shifts from slice 5 to slices 6, 8, 11, and 11, the slice with the maximum frictional resistance remains at slice 8 consistently, and the slice with the maximum surplus thrust changes from slice 5 to slices 6, 10, 11, and 12. Overall, the maximum driving force increases from 2,934 kN to 5,435 kN and the surplus thrust grows from 1,459 kN to 4,475 kN, whereas the maximum slice frictional resistance decreases from 1,563 kN to 1,116 kN. This decrease is attributed to the reduction in frictional resistance caused by an increasing reduction coefficient. Regarding the global stability coefficient evaluation methods (CRSM, MTM, CDM, and SDM), during the transition of the critical state from slice 10 to the final slice 24: the CRSM values decrease from 1.81, 4.88, 2.43 to 0.55, 2.14, 0.73 in X, Y and resultant directions respectively; the MTM values decline from 1.29, 0.32, 1.03 to 0.0, 0.0, 0.0; the CDM values reduce from 5.43, 5.35, 5.41 to 1.0, 1.0, 1.0; and the SDM values drop from 0.58, 0.12, 0.42 to 0.0, 0.0, 0.0. The complete evolutionary processes are illustrated in Figures 15–18 respectively.

Under the condition that the peak shear stress and the corresponding shear strain are equal to those defined in theory, the newly developed complete process constitutive model (CPCM) is employed to simulate the progressive failure process. The initial critical-state slice is identified as slice 11. With increasing deformation at the landslide’s rear edge, the critical state progressively advanced along the slip surface. When the final slice reached critical state, the landslide entered global failure. During this progressive failure process, five critical-state slices were selected, and their corresponding driving force, frictional resistance, and surplus thrust were plotted, as shown in Figures 20–24. As displacement continues to increase and the critical state slice changes, the evolution patterns of four slope stability coefficients ($F_{CRSM}$, $F_{MTM}$, $F_{CDM}$, and $F_{SDM}$) are illustrated in Figures 25–28, respectively.

Figure 20

Figure 20. Force distribution of all slices when slice 12 is in the critical state under the complete process constitutive model.

Figure 21

Figure 21. Force distribution of all slices when slice 15 is in the critical state under the complete process constitutive model.

Figure 22

Figure 22. Force distribution of all slices when slice 18 is in the critical state under the complete process constitutive model.

Figure 23

Figure 23. Force distribution of all slices when slice 21 is in the critical state under the complete process constitutive model.

Figure 24

Figure 24. Force distribution of all slices when slice 24 is in the critical state under the complete process constitutive model.

Figure 25

Figure 25. CRSM stability coefficient under the complete process constitutive model.

Figure 26

Figure 26. MTM stability coefficient under the complete process constitutive model.

Figure 27

Figure 27. CDM stability coefficient under the complete process constitutive model.

Figure 28

Figure 28. SDM stability coefficient under the complete process constitutive model.

As shown in Figures 20–24, with the progressive advancement of critical states, the maximum driving force shifts from slice 6 to slices 6, 7, and 8; the maximum frictional resistance remains consistently at slice 8; while the surplus thrust progresses from slice 4 to slices 5, 6, 7, and 11. Overall, the maximum driving force increases from 3,344 kN to 4,790 kN and the surplus thrust grows from 1,944 kN to 3,595 kN, whereas the maximum frictional resistance decreases from 1,641 kN to 1,330 kN due to frictional resistance reduction calculated by the new constitutive model with increasing displacement.

Using the new evaluates model global stability coefficient of slope (CRSM, MTM, CDM, SDM), during the progression from critical Slice 11 to final Slice 24: the CRSM values decrease from 1.55, 4.15, 2.08 to 0.51, 2.13, 0.71 in X, Y and resultant directions respectively; the MTM values decline from 0.81, 0.18, 0.63 to 0.0, 0.0, 0.0; the CDM values reduce from 5.25, 2.97, 3.72 to 1.0, 1.0, 1.0; and the SDM values drop from 0.14, 0.02, 0.08 to 0.0, 0.0, 0.0. The complete evolutionary processes are illustrated in Figures 25–28 respectively.

From the analysis based on the traditional slope stability coefficient method, the partial strength reduction method (PSRM), the perfect elasto-plastic model (PEPM), and the complete process constitutive model (CPCM), the following insights can be drawn: the stability factor obtained by the traditional slice method represents the overall stability factor of the sliding mass; the surplus stability factor obtained by the partial strength reduction method reflects the evolution pattern of the surplus stability degree of the sliding mass as the critical state shifts; the main thrust method and surplus displacement method represent the surplus degree of the overall force and displacement of the sliding mass, respectively; the comprehensive sliding force–resisting force method and comprehensive displacement method describe the evolution patterns of the overall force and displacement of the landslide as displacement changes.

Regarding the mechanisms described by these evolution patterns, traditional methods primarily use force to characterize the stability of the sliding mass. A comparative analysis of the new methods (partial strength reduction method, ideal elasto-plastic model method, andcomplete process constitutive model method) yields the following results under current conditions: the stability coefficient obtained by the traditional strength reduction method is 1.50, whereas the perfect elasto-plastic partial reduction model generates respective stability coefficients of (1.81, 4.88, 2.43) for CRSM, (1.29, 0.32, 1.03) for MTM, (5.43, 5.35, 5.41) for CDM, and (0.58, 0.12, 0.42) for SDM. The complete process constitutive model produces corresponding values of (1.55, 4.15, 2.08), (0.81, 0.18, 0.63), (5.25, 2.97, 3.72), and (0.14, 0.02, 0.08). Results indicate higher stability coefficients from the perfect elasto-plastic partial reduction method relative to the complete process constitutive model. The conventional stability coefficient of 1.5 suggests that the slope is stable, yet field observations reveal cracks at both the front and rear edges of the slope. From the multi-parameter stability evaluation analysis of the slope, the stability coefficients in the X-direction for the main thrust method and the surplus displacement method are 0.32 and 0.18, respectively, indicating that the slope is unstable.

8 Conclusion

To investigate the process of landslide evolution from initial cracking to overall failure, this study conducted a detailed analysis of the slope stability evolution using both the partial strength reduction method with an perfect elasto-plastic model and the complete process constitutive model. By extending the traditional slice method, the stability coefficient is now closely linked to deformation. Under the condition that the peak stress and the corresponding strain are equal for both models, the following conclusions are drawn:

1. Based on the traditional strength reduction method, this study further clarifies the computational steps and physical interpretation of the partial strength reduction method. It demonstrates that the partial strength reduction method can approximately describe the evolution of forces during the progressive failure of a slope.

2. The traditional assumptions of the slice method are expanded by incorporating both ideal elasto-plastic and complete process constitutive models into the slice analysis framework. The partial strength reduction perfect elasto-plastic model and complete process constitutive model provide a comprehensive description of slope progressive failure. The complete process shear stress-strain constitutive model effectively captures the force and displacement characteristics during progressive failure. Notably, all four stability coefficients derived from the partial strength reduction perfect elasto-plastic model (comprehensive sliding-resisting force method, main thrust method, comprehensive displacement method, and surplus displacement method) exceed those obtained from the complete process shear stress-strain constitutive model.

3. Two critical-state slice determination methods are proposed: the slice division method and the surplus force calculation method. The deformation relationship between adjacent slices in the failure zone is examined, leading to the formulation of deformation compatibility conditions for interconnected slices. The generalized slice analysis method for progressive slope failure developed in this study can be extended to other slice methods, including Janbu’s method, Sarma’s method, and the simplified Bishop method.

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

SH: Data curation, Formal Analysis, Methodology, Writing – original draft. LL: Investigation, Methodology, Project administration, Supervision, Writing – review and editing. LY: Investigation, Software, Writing – original draft. YL: Investigation, Validation, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research received support from the National Natural Science Foundation of China (Grants No. 42071264 and 41641027), for which the authors are grateful.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

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

Baker, R. (2004). Nonlinear mohr envelopes based on triaxial data. J. Geotechnical Geoenvironmental Engineering 130 (5), 498–506. doi:10.1061/(ASCE)1090-0241(2004)130:5(498)

CrossRef Full Text | Google Scholar

Bishop, A. W. (1955). The use of the slip circle in the stability analysis of slopes. Geotechnique 5 (1), 7–17. doi:10.1680/geot.1955.5.1.7

CrossRef Full Text | Google Scholar

Faramarzi, L., Zare, M., Azhari, A., and Tabaei, M. (2017). Assessment of rock slope stability at cham-shir dam power plant pit using the limit equilibrium method and numerical modeling. Bull. Eng. Geol. Environ. 76 (2), 783–794. doi:10.1007/s10064-016-0870-x

CrossRef Full Text | Google Scholar

Fellenius, W. (1936). Calculation of the stability of earth dams. Proc. Second Congress Large Dams 4, 445–463.

Google Scholar

Gharti, H. N., Komatitsch, D., Oye, V., Martin, R., and Tromp, J. (2012). Application of an elastoplastic spectral-element method to 3d slope stability analysis. Int. J. Numer. Methods Eng. 91 (1), 1–26. doi:10.1002/nme.3374

CrossRef Full Text | Google Scholar

Hrennikoff, A. (1941). Solution of problems of elasticity by the framework method. J. Appl. Mech. 8, A169–A175. doi:10.1115/1.4009129

CrossRef Full Text | Google Scholar

Janbu, N. (1973). Slope stability computations. Wiley (John) and Sons.

Google Scholar

Kelesoglu, M. (2016). The evaluation of three-dimensional effects on slope stability by the strength reduction method. KSCE J. Civ. Eng. 20 (1), 229–242. doi:10.1007/s12205-015-0686-4

CrossRef Full Text | Google Scholar

Lu, Y. (2015). Deformation and failure mechanism of slope in three dimensions. J. Rock Mech. Geotechnical Eng. 7 (2), 109–119. doi:10.1016/j.jrmge.2015.02.008

CrossRef Full Text | Google Scholar

Lu, Y., Yang, L., and Liu, D. (2013). A new joint constitutive model and several new methods of stability coefficient calculation of landslides. Chin. J. Rock Mech. Eng. 32 (12) 2431–2438.

Google Scholar

Lu, Y., Huang, X., and Liu, D. (2017). Distribution characteristics of force and stability analysis of slope. Chin. J. Geotechnical Eng. 30 1322–1330.

Google Scholar

Lu, Y., Sun, W., Yang, H., Jiang, J., and Lu, L. (2023). A new calculation method of force and displacement of retaining wall and slope. Appl. Sci. 13 (9), 5806. doi:10.3390/app13095806

CrossRef Full Text | Google Scholar

Lu, L., Yang, H., Lu, Y., and Zhang, S. (2025). Multi-criterion stability evaluation and surface boundary method for slopes. Adv. Civ. Eng. 2025 (1), 6588377. doi:10.1155/adce/6588377

CrossRef Full Text | Google Scholar

Matthews, C., Farook, Z., and Helm, P. (2014). Slope stability analysis–limit equilibrium or the finite element method. Ground Eng. 48 (5), 22–28. Available online at: https://www.geplus.co.uk/technical-paper/technical-paper-slope-stability-analysis-limit-equilibrium-or-the-finite-element-method-01-05-2014/.

Google Scholar

Mohammadi, S. (2008). Extended finite element method: for fracture analysis of structures. doi:10.1002/9780470697795

CrossRef Full Text | Google Scholar

Nian, T.-K., Huang, R.-Q., Wan, S.-S., and Chen, G.-Q. (2012). Three-dimensional strength-reduction finite element analysis of slopes: geometric effects. Can. Geotechnical J. 49 (5), 574–588. doi:10.1139/t2012-014

CrossRef Full Text | Google Scholar

Oden, J. T., Duarte, C., and Zienkiewicz, O. C. (1998). A new cloud-based hp finite element method. Comput. Methods Applied Mechanics Engineering 153 (1-2), 117–126. doi:10.1016/S0045-7825(97)00039-X

CrossRef Full Text | Google Scholar

Planas, J., Elices, M., Guinea, G., Gómez, F., Cendón, D., and Arbilla, I. (2003). Generalizations and specializations of cohesive crack models. Eng. Fracture Mechanics 70 (14), 1759–1776. doi:10.1016/S0013-7944(03)00123-1

CrossRef Full Text | Google Scholar

Sarma, S., and Tan, D. (2006). Determination of critical slip surface in slope analysis. Geotechnique 56 (8), 539–550. doi:10.1680/geot.2006.56.8.539

CrossRef Full Text | Google Scholar

Spencer, E. (1967). A method of analysis of the stability of embankments assuming parallel inter-slice forces. Geotechnique 17 (1), 11–26. doi:10.1680/geot.1967.17.1.11

CrossRef Full Text | Google Scholar

Tiwari, R. C., Bhandary, N. P., and Yatabe, R. (2015). 3d sem approach to evaluate the stability of large-scale landslides in Nepal himalaya. Geotechnical Geol. Eng. 33 (4), 773–793. doi:10.1007/s10706-015-9858-8

CrossRef Full Text | Google Scholar

Wei, W., Cheng, Y. M., and Li, L. (2009). Three-dimensional slope failure analysis by the strength reduction and limit equilibrium methods. Comput. Geotechnics 36 (1-2), 70–80. doi:10.1016/j.compgeo.2008.03.003

CrossRef Full Text | Google Scholar

Wilson, Z. A., Borden, M. J., and Landis, C. M. (2013). A phase-field model for fracture in piezoelectric ceramics. Int. J. Fract. 183 (2), 135–153. doi:10.1007/s10704-013-9881-9

CrossRef Full Text | Google Scholar

Yingfa, L. U. (2016). A new constitutive model and its parameter calibration. Rock Soil Mech. 37 (o.8), 2138–2144.

Google Scholar

Zhang, Y., Lu, Y., Zhong, Y., Li, J., and Liu, D. (2023). Stability analysis of landfills contained by retaining walls using continuous stress method. CMES-Computer Model. Eng. and Sci. 134 (1), 357–381. doi:10.32604/cmes.2022.020874

CrossRef Full Text | Google Scholar

Zhang, S., Lu, Y., and Lu, L. (2025). A new method for evaluating the stability of retaining walls. Buildings 15 (10), 1732. doi:10.3390/buildings15101732

CrossRef Full Text | Google Scholar

Zhu, D., and Lee, C. (2002). Explicit limit equilibrium solution for slope stability. Int. J. Numer. Anal. Methods Geomechanics 26 (15), 1573–1590. doi:10.1002/nag.260

CrossRef Full Text | Google Scholar

Nomenclature

Latin letters

CDM comprehensive displacement method

CRSM comprehensive sliding force - resisting force method

CSN critical state slice block

CPCM shear stress-strain complete process constitutive model

D-CPCM determination of critical state for the shear stress-strain complete process constitutive model

${\mathit{F}}_{\mathit{C}\mathit{S}\mathit{R}\mathit{M}}^{\mathit{i}}$ safety factor components of CRSM in the X, Y and their resultant’s directions, respectively

${\mathit{F}}_{\mathit{C}\mathit{D}\mathit{M}}^{\mathit{i}}$ safety factor components of CDM in the X, Y and their resultant’s directions, respectively

${\mathit{F}}_{\mathit{M}\mathit{T}\mathit{M}}^{\mathit{i}}$ safety factor components of MTM in the X, Y and their resultant’s directions, respectively

${\mathit{F}}_{\mathit{S}\mathit{D}\mathit{M}}^{\mathit{i}}$ safety factor components of SDM in the X, Y and their resultant’s directions, respectively

${\mathit{G}}_{\mathit{i}}$ shear modulus

GM rigid model

L critical state slice block

MC Mohr-Coulomb criterion

MTM main thrust method

${\mathit{N}}_{\mathit{i}}$ normal pressure of the ith unit slice block

${\mathit{N}}_{\mathit{m}}^{\mathit{N}}$ normal pressure of the refined mth unit slice block

PEPM perfect elasto-plastic model

${\mathit{P}}_{\mathit{i}}$ unbalanced thrust of the ith unit slice block

${\mathit{P}}_{\mathit{x}\mathit{m}},{\mathit{P}}_{\mathit{y}\mathit{m}}$ surplus sliding force components along the sliding surface in the X and Y directions under the current stress state, respectively

${\mathit{P}}^{\mathit{m}}$ resultant of $P_{xm},P_{ym}$

PSRM partial strength reduction method

${{\mathit{P}}_{\mathit{i}}}^{\mathbf{S}}$ sliding force of the ith unit slice block

${{\mathit{P}}_{\mathit{m}}}^{\mathbf{S},\mathbf{N}}$ sliding force of the redefined mth unit slice block

${\mathit{P}}_{\mathit{i}}^{\mathit{T}\mathit{C}\mathit{M}}$ unbalanced thrust of the ith slice block for TCM

${\mathit{P}}_{\mathit{i}}^{\mathit{S}\mathit{T}\mathit{C}\mathit{M}}$ sliding force of the ith slice block for TCM

${\mathit{P}}_{\mathit{x}}^{\mathit{S}},{\mathit{P}}_{\mathit{y}}^{\mathit{S}}$ sliding force components along the sliding surface in the X and Y directions under the current stress state, respectively

${\mathit{P}}^{\mathit{S}}$ resultant of $P_x^S,P_y^S$

SDM surplus displacement method

${\mathit{S}}_{\mathit{x}}^{\mathit{s}},{\mathit{S}}_{\mathit{y}}^{\mathit{s}},{\mathit{S}}^{\mathit{s}}$ displacement components along the sliding surface in the X, Y and their resultant directions under the current state, respectively

${\mathit{S}}_{\mathit{x}},{\mathit{S}}_{\mathit{y}},\mathit{S}$ displacement components along the sliding surface in the X, Y and resultant’s directions under the failure state, respectively

SN Number of slice block

${\mathit{S}}_{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}^{\mathit{x}},{\mathit{S}}_{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}^{\mathit{y}},{\mathit{S}}_{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}$ displacement difference components along the sliding surface in the X, Y and resultant’s directions, respectively

SRM Strength reduction method

TCM traditional critical state slice method

${\mathit{T}}_{\mathit{i}}^{\mathit{T}\mathit{C}\mathit{M}}$ frictional resistance of the ith slice block for TCM

${\mathit{T}}_{\mathit{i}}$ anti-sliding force of the ith unit slice block

${{\mathit{T}}_{\mathit{i}}}^{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}$ peak anti-sliding force of the ith unit slice block

${{\mathit{T}}_{\mathit{m}}}^{\mathit{c}\mathit{r}\mathit{i}\mathit{t},\mathit{N}}$ peak anti-sliding force of the redefined mth unit slice block

${\mathit{T}}_{\mathit{i},\mathit{F}}^{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}$ reduced anti-sliding force of the ith unit slice block

${\mathit{T}}^{\mathit{x}\mathit{m}},{\mathit{T}}^{\mathit{y}\mathit{m}},{\mathit{T}}^{\mathit{m}}$ surplus anti-sliding force components along the sliding surface in the X and Y directions and their resultant under the failure stress state, respectively

${\mathit{T}}^{\mathit{x}},{\mathit{T}}^{\mathit{y}}$ anti-sliding components along the sliding surface in the X and Y directions under the failure stress state

$\mathit{T}$ resultant of $T^x,T^y$

UTM unbalanced thrust slice method

${\mathit{W}}_{\mathit{i}}$ weight of the ith unit slice block

${\mathit{W}}_{\mathit{m}}^{\mathit{N}}$ weight of the redefined mth unit slice block

${\mathit{a}}_{\mathbf{1}}^{\mathbf{0}},{\mathit{a}}_{\mathbf{2}}^{\mathbf{0}}$ constant coefficients

${\mathit{c}}_{\mathit{i}}$ cohesive force of the ith unit slice block(or element)

$\mathit{f}$ strength reduction coefficient along the entire sliding surface

${\mathit{f}}_{\mathit{i}}$ partial strength reduction coefficient

${\mathit{f}}_{\mathit{z}\mathit{s}}^{\mathit{i}}$ partial strength reduction surplus stability coefficient

$\mathit{i}\left({\mathit{x}}_{\mathit{i}}^{\mathit{m}},{\mathit{y}}_{\mathit{i}}^{\mathit{m}}\right),\mathit{i}\in \left(\mathbf{1},\mathbf{6}\right)$ coordinate of the mth unit slice block at i point

${\mathit{l}}_{\mathit{m}}^{\mathit{N}}$ length of the refined mth unit slice block

${\mathit{l}}_{\mathit{i}}$ length of the bottom edge of the ith unit slice block

n tatol number of element (or slice block)

${\mathit{p}}_{\mathit{i}},{\mathit{q}}_{\mathit{i}},{\mathit{\varsigma }}_{\mathit{i}}$ parameters of constitutive equation of the ith unit slice block

${\mathit{\xi }}_{\mathit{i}}$ softening parameters of constitutive equation of the ith unit slice block

${\mathit{\xi }}_{\mathit{i},\mathit{c}},{\mathit{\xi }}_{\mathit{i},\mathit{o}}$ constant coefficients

${\mathit{\beta }}_{\mathit{i}}$ uniformly distributed vertical load of the ith unit slice block on the surface of the earth

${\mathbf{\Delta }}_{\mathit{i}}$ uniformly distributed horizontal load of the ith unit slice block on the surface of the earth

x coordinate in X direction

y coordinate in Y direction

Greek letters

${\mathit{\alpha }}_{\mathit{i}}$ angle between the bottom edge and the horizontal of the ith unit slice block

${\mathit{\alpha }}_{\mathit{f}}$ angle of T with respect to X-axis

${\mathit{\alpha }}_{\mathit{f}}^{\mathit{m}}$ angle of T^m with respect to X-axis

${\mathit{\alpha }}_{\mathit{S}}$ angle of P^S with respect to X-axis

${\mathit{\alpha }}_{\mathit{s}}^{\mathit{s}}$ angle of $S_x^s$ with respect to $S_y^s$

${\mathit{\alpha }}_{\mathit{s}}^{\mathit{c}\mathit{r}\mathit{i}\mathit{t}}$ angle of S with respect to X-axis

${\mathit{\alpha }}^{\mathit{m}}$ angle of $P^m$ with respect to X-axis

${\mathit{\gamma }}_{\mathit{i}}$ shear strain along the bottom edge of the ith unit slice block

${\mathit{\gamma }}_{\mathit{i}}^{\mathit{S}}$,${\mathit{\gamma }}_{\mathit{i}}^{\mathit{n}}$ shear strain along the bottom edge and the horizontal of the ith unit slice block

${\mathit{\gamma }}_{\mathit{i},\mathit{p}\mathit{e}\mathit{a}\mathit{k}}^{*}$ critical shear strain induced by strength reduction of the ith unit slice block

${\mathit{\gamma }}_{\mathit{i},\mathit{p}\mathit{e}\mathit{a}\mathit{k}}$ critical shear strain at the peak stress of the ith unit slice block

${\mathit{\gamma }}_{\mathit{i}}^{\mathit{p},\mathit{b}}$ failure shear strain of the ith unit slice block

${\mathit{\sigma }}_{\mathit{i},\mathit{n}}$ normal stress of the ith unit slice block

${\mathit{\sigma }}_{\mathit{i},\mathit{n}}^{\mathit{c}}$ constant normal stress of the ith unit slice block

${\mathit{\sigma }}_{\mathit{i},\mathit{n}}^{\mathit{p},\mathit{b}}$ normal stress of the ith unit slice block under failure stress state

${\mathit{\sigma }}_{\mathit{i},\mathit{n}}^{\mathit{X}}$ current normal stress of the ith unit slice block

${\mathit{\sigma }}_{\mathit{i},\mathit{n}}^{\mathit{p},\mathit{b}}$ failure normal stress of the ith unit slice block

${\mathit{\tau }}_{\mathit{i}}$ shear stress of the bottom edge of the ith unit slice block

${\mathit{\tau }}_{\mathit{i},\mathit{f}}$ reduced stress of the ith unit slice block

${\mathit{\tau }}_{\mathit{i},\mathit{p}\mathit{e}\mathit{a}\mathit{k}}$ peak anti-sliding stress of the ith unit slice block

${\mathit{\tau }}_{\mathit{i}}^{\mathit{X}}$ current frictional shear stress of the ith unit slice block

${\mathit{\tau }}_{\mathit{i}}^{\mathit{p},\mathit{b}}$ failure frictional shear stress of the ith unit slice block

${\mathit{\tau }}_{\mathit{i}}^{\mathit{u},\mathit{p},\mathit{b}}$ failure driving shear stress of the ith unit slice block

${\mathit{\phi }}_{\mathit{i}}$ internal friction angle of the ith unit slice block