An Improved Strength Reduction-Based Slope Stability Analysis

Abstract

Slope uncertainty predominantly originates from the imperfect analysis model and the inaccuracy and imprecision of the observations. The strength reduction method (SRM) is widely used to attain the safety factor (SF) of the slopes, which is similar to interpretation of the limit state (LS). In this paper, the spectral element method (SEM), using an elasto-plastic Mohr–Coulomb failure criterion, is employed to project the plausible LS of the soil slopes. An iterative SRM search method is proposed to evaluate the SF of the slopes regardless of the LS interpretation. The proposed SRM paradigm encompasses the design trigger to trace the uncertain parameters in decision-making. This method is applied to three numerical examples: (1) a homogeneous dry slope, (2) a dry slope with a weak layer, and (3) a partially-wet slope with a weak layer. It is shown that for the case study examples, the proposed SRM reasonably converges to the required precision. Results further are compared and contrasted with some of the conventional and standard techniques in slope stability. This hybrid procedure paves the road for fast and safe stability analysis of man-made and natural slopes.

1. Introduction

Slope stability analysis is one of the most important research fields in geoscience. Slope stability models are associated with many uncertainty sources, such as material properties [1], simulation methods [2], failure path, and applied loads. Numerous researchers have proposed a wide spectrum of analysis techniques (ranging from simple/approximate to detailed/complex numerical models) to evaluate the slope stability. From the perspective of the soil stability analysis, there are two widely-used analysis methods: (1) the limit equilibrium method (LEM) and (2) the strength reduction method (SRM).

Some of the fundamental research related to the LEM can be found in Spencer [3], Sarma [4]. The limitations of the LEM have been well discussed by Krahn [5]. The LEM is based on a predefined hypothetical failure surface with a rigid (i.e., unconformable) super-structure [6]. It is assumed that there are two forces acting on the rigid body, i.e., the driving force, which attempts to move the rigid body, and the resisting force which prevents this action. The slope is considered to be stable as long as the resisting forces are greater than or equal to the driving ones, and their ratio is called the safety factor (SF). Therefore, the LEM is deemed a deterministic approach for soil stability analysis [7].

On the other hand, the application of SRM has also gained popularity in slope stability assessment. It was employed by Dawson et al. [8], Griffiths and Lane [9], Matsui and San [10], and Ugai [11] for the calculation of the safety factor in soil slopes. In SRM, the shear strength of the soil is reduced to trigger the failure of the slope. The reduction factor imposed on the system at the failure threshold is designated as the safety factor. Although many strategies have been introduced for the definition of the trigger of the failure [12,13], the deterministic assumption and engineering judgment imply the model uncertainties.

Although many of the current slope stability applications are focused on the numerical discretization with the crude finite element method (FEM), several variations are also proposed, e.g., spectral element method (SEM), extended FEM (XFEM), boundary element method (BEM), smoothed FEM (SFEM), etc. Of particular interest is SEM to solve the boundary value problems. The SEM is, in fact, a high-order form of the FEM suitable for geotechnical engineering problems. Gharti et al. [14] introduced the SEM for the evaluation of 3D soil slope stability. Many other researchers have shown the application of SEM in different geomechanics and geotechnical engineering problems. Among them are: the effect of geometrics on 3D slope stability by Zhang et al. [15], the assessment of the vegetated and barren soil slopes by Tiwari et al. [16], the seismic slope stability of natural slopes by Tiwari et al. [17], the evaluation of long and steep slopes by Tiwari et al. [18], simulations of the elastic wave propagation in exploration by Zheng et al. [19], and 3D rock stability by Shen and Karakus [20].

In this study, the stability of the soil slope is evaluated based on the spectral element method. An approach rooted in the design criteria is proposed to find the safety factor under the strength reduction method. First, the fundamentals of the slope stability analysis, spectral element method and strength reduction method are discussed. Next, the proposed method is explained, and several case study problems are analyzed to show the advantages of this hybrid method.

2. Fundamental Theories

2.1. Spectral Element Method

SEM is constituted of three general parts: pre-, main, and post-processing. The pre-processing itself consists of geometry, material attributes, and boundary conditions. Nodes’ coordinates and elements’ connectivities are the most important parameters that need to be defined for a proper application. The unit weight, $\gamma$, Young’s modulus, E, internal friction angle, $\varphi$, dilation angle, $\psi$, and cohesion, c, are classified as material properties. Moreover, the boundary conditions must be applied carefully in order to develop a reliable simulation.

SEM is a flexible technique, which can adapt a variety of the element types and the arbitrary order of integration points via Gaussian quadrature from FEM [14]. Figure 1 illustrates the fundamental difference between SEM and FEM. Both of them employ the same number of Lagrangian interpolation points in identical locations within the element, while the integration points are different. The readers should note that the provided mass matrices are used only in the dynamic analysis, and they are shown for the completeness of the comparison between SEM and FEM.

SEM utilizes Gauss–Legendre–Lobatto (GLL) points as a nodal quadrature for which interpolation nodes coincide with the integration points. In particular, additional integration points in FEM, which are called Gauss–Legendre (GL) points, are generated at different positions compared to the interpolation points. This innovative change in the configuration of the SEM improves its capabilities in two ways. First, removing the interpolation to calculate the nodal quantities from quadrature points makes the problem computationally simple. Second, interpolating functions become orthogonal on quadrature points, resulting in a diagonal mass matrix (used only in the dynamic analysis). Both changes increase the accuracy and decrease the computational effort. From the constitutive relation point of view, the elasto-plastic Mohr–Coulomb failure criterion is used to simulate the material behavior.

The unstructured hexahedron mesh technique is used to prepare the desired geometry for subsequent simulation. Then, the GLL nodes are introduced in each element. The nodes, interpolation points, and GLL points of an element are shown in Figure 2. The first step in mesh generation based on a hexahedral element is to divide the full geometrical model into several elements with common surface and corner nodes; see Figure 2a. Next, the interpolation and GLL points are mapped by SEM for each element; see Figure 2b,c. Finally, the boundary conditions along with external loads are applied.

2.2. Strength Reduction Method

Slope stability requires a series of nonlinear elasto-plastic simulations based on Mohr–Coulomb criteria. In the strength reduction method, the actual strength parameters are reduced by a certain reduction factor as [21]:

$$\begin{array}{ccc}\hfill c^{\prime }& =& \frac{c}{SRF}\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\varphi }^{\prime }& =& arctan\left(\frac{tan\varphi }{SRF}\right)\hfill \end{array}$$

(2)

where SRF is a strength reduction factor.

To obtain an accurate safety factor, it is important to trace the SRF gradually until the reduced strength parameters, $c^{\prime }$ and ${\varphi }^{\prime }$, yield the predefined limit state (LS). This traditional approach makes the problem iterative and subsequently increases the analysis time and effort. However, in the proposed method, an intelligent method is used to update the new SRF based on the old values. The SRF corresponding to an LS is designated as SF for the slope (i.e., SF = SRF). Therefore, a legitimate interpretation of the LS is important. There are various LS criteria to identify the slope failure such as experimental tests [22], limiting the shear stress on the potential failure surface [23], divergence of the computation problem [8], calculation of the total equivalent plastic strain zone [24], running-through of the yielding zone, or a prescribed generalized shear strain developed from the toe through the crest of the slope [25], a rapid increase in displacement [9], and the formulation of an initial value problem to trace the critical slide line [26].

The impacts of the material properties on the overall slope stability are also studied. Zheng et al. [27] and Duan et al. [26] investigated the effects of Young’s modulus and Poisson’s ratio on the soil stability analysis. Dawson et al. [8] claimed that the value of elasticity parameters (i.e., E and $\nu$) had little impact on the SF. Therefore, the effects of elasticity parameters will not be discussed in this paper.

3. Methodology

The strategy to implement the SEM on the soil slopes is illustrated in Figure 3. This procedure can be categorized into four steps. First, the geometrical model and the material properties are defined as input parameters, which can be obtained by surveying techniques and the in situ tests/measurements, respectively. Next, the computational model is built to obtain a realistic response in which the SEM and the Mohr–Coulomb relation are employed. Different outputs can be extracted; however, the maximum displacement is the most palpable metric. Finally, the SRM is iteratively repeated to find the target SF.

The proposed improved strength reduction method is presented in Figure 4. This flowchart shows that the main matrices in the SEM processing are displacement, $\left[U\right]$, stiffness, $\left[K\right]$, and the force vector, $\left[F\right]$, which are developed in the pre-processing stage. Obviously, the processing stage is the most computationally demanding part, which is associated with the number of plastic nodes. The outputs are nodal displacement, deformed shapes, and the number of iterations.

It is noteworthy that the open-source software SPECFEM3D-GEOTECH was used in the present study, which can be found via Computational Infrastructure for Geodynamics (www.geodynamics.org). The software is parallelized on the basis of domain decomposition [14]. It can handle the material heterogeneity and complex models with significant topography. Both simple and variable free surface profiles may be used on the basis of a simple hydrostatic relation. To be used in the context of this paper, the above-mentioned software is coupled with [28] to perform the iterative simulations and automatically post-process the results.

As previously stated, the maximum displacement is used to evaluate SF based on the SRF. The proposed method is an automatic search technique, which is started with three SRFs, i.e., $SR{F}_1=1.0$, $SR{F}_2=1.5$, $SR{F}_3=2.0$. Then, the maximum displacements, $U_i$, corresponding to $SR{F}_i$, are calculated. Figure 5a shows the variation of the maximum displacement to initialize the search parameters. The required precision (RP) is defined in order to automatize the process of SF determination. In other words, the iterative process will be repeated until the pre-defined RP is satisfied (it is set to 0.01 in this study). The SF corresponds to a rapid increase in the displacement response. The gradient of the displacement vs. SRF is considered as a design trigger (DT) to identify the LS. The DT is used to alleviate the aleatory selection of SF. Theoretically, the DT ranges from zero to infinity. However, in this study, three values are adopted, i.e., 30%, 50%, and 70%. The left and right gradients (denoted as MLand MR) are computed by Equation (3) to satisfy the DT constraint; see Figure 5b.

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill ML& =\hfill & \frac{U_2-U_1}{SR{F}_2-SR{F}_1}\\ \hfill MR& =\hfill & \frac{U_3-U_2}{SR{F}_3-SR{F}_2}\end{array}\end{array}$$

(3)

4. Case Study Models

Several examples have been simulated to evaluate and validate the slope stability analysis problems [14,29,30]. In this study, three slope stability problems are discussed: (1) homogeneous dry slope, (2) dry slope with a weak layer, and (3) partially-wet slope with a weak layer. All three case studies are symmetric in the XZ-plane [14]. Boundary conditions (BCs) are applied for half of the model in SEM; see Figure 6a. All the side faces are constrained in the perpendicular direction. The symmetry plane (back face) is free to move in the X and Z directions. The front face is further constrained in the vertical direction (due to a hypothetical valley). Since the bottom face is connected to the Earth, it is fully constrained.

5. Results and Discussion

5.1. Case 1

The dimensions and geometry of the first case study are illustrated in Figure 7a, which consist of 1670 elements and 2200 nodes. This model represents the dry region with material properties of $\gamma$ = 18.8 kN/m${}^3$, $\varphi$ = 20${}^{°}$, $c=$ 29 kN/m${}^2$, $\psi$ = 0, E = ${10}^5$ kN/m${}^2$, and $\nu =$ 0.3.

5.2. Cases 2 and 3

The dimensions and geometry of Case 2 and Case 3 are shown in Figure 7b, which consist of two different material types with 1890 elements and 2453 nodes. The highlighted weak layer is made of material properties with $\gamma$ = 18.8 kN/m${}^3$, $\varphi$ = 0${}^{°}$, $c=$ 10 kN/m${}^2$. The slope material properties are identical to Case 1.

5.3. Safety Factor Determination

The proposed algorithm in Section 3 was applied to the three case studies. In each case, three different DT margins were adopted, which were subsequently used to update the ML and MR values. The iteration process was continued until the RP reached 0.01.

In general, for the smaller SRF values, a simulation requires less iterations to fulfill the process (the system still behaves in the linear range). Increasing the SRF increases the number of plastic nodes. For example, in Case 1, with DT = 30%, the SEM required only 120 iterations for the system with SRF of 2.16 (elastic material), while it required more than 350 iterations for the system with SRF of 2.18 (which underwent plastic deformations). In Case 2 (which had a weak layer), the differences between the elastic and plastic iterations were even greater. For Case 2 with DT = 30%, the SEM required 69 and 1311 iterations for an SRF of 1.3 and 1.7, respectively.

Figure 7 compares the slope stability analysis of the three case studies. In each case, the results of the numerical simulations are presented in terms of the maximum displacement, ${\Delta }_{max}$, and the percentage of the gradient. The linear part of the system exhibited the elastic behavior of the system, while the nonlinear range presented the plasticity. Moreover, the identified SF in the ${\Delta }_{max}$-SRF coordinate system is shown with a dashed line for each DT value. Several observations can be drawn:

• The displacement-based curves were smoother than the gradient-based ones.
• In displacement-based curves, the three plots resulting from three DT values were very close to each other. However, they had meaningful differences in the gradient-based plots.
• The range of SRF in which the system moved from the linear to plastic phase was longer for Cases 2 and 3 than Case 1.
• In all case studies, the safety factor for the system with DT = 30% was the minimum, and for the system with DT = 70%, it was the maximum.
• The safety factor of the homogeneous model (i.e., Case 1) was more than the two other cases (with a weak layer). This is consistent with the physics of the problem.
• Comparing Cases 2 and 3 revealed that the safety factor of the dry model (Case 2) was higher than the partially-wet slope (i.e., Case 3).
• The number of required simulations depended on the LS and pre-defined accuracy.

5.4. Validation

The three case studies in this paper were already examined by other researchers to find the SF.

Table 1. Comparison of the safety factor for the three case studies using different techniques.

Methods	Case 1	Case 2	Case 3
Xing [32]	2.122	1.553	1.441
Bishop’s modified method (Lam and Fredlund [31])	-	1.607	1.511
Janbu’s simplified method (Lam and Fredlund [31])	-	1.558	1.481
CLARA (Lam and Fredlund [31])	-	1.62	1.54
CLE (Lam and Fredlund [31])	-	1.603	1.508
Chen et al. [30]	2.262	1.717	-
Chen et al. [30]	2.187	1.603	-
FEM (Griffiths and Marquez [33])	2.17	1.58	-
SEM (Gharti et al. [14])	2.18	1.57	1.49
SEM (DT = 30%)—present study	2.17	1.59	1.46
SEM (DT = 50%)—present study	2.18	1.61	1.48
SEM (DT = 70%)—present study	2.19	1.62	1.5

compares the SFs obtained by different methods: the limit equilibrium method of Lam and Fredlund [31], Xing [32], the limit analysis method of Chen et al. [29,30], the finite element method of Griffiths and Marquez [33], the spectral element method of Gharti et al. [14], and the present study.

The results of the LEM varied owing to the conceptual assumptions, the complexity of the failure surface, and the force decomposition assumption on the column direction, particularly in three-dimensional slope stability analysis. This is due to the deterministic assumption of the LEM, which was discussed in the “Introduction” Section. The following conclusions can be deduced from the results.

• The SFs obtained by Gharti et al. [14] were approximately equal to the proposed hybrid SEM-SRM with DT of 50%, 10%, and 60% for Case Studies 1, 2, and 3, respectively. This indicates that the inconsistency in SF resulted from the conventional SRM.
• The differences between SFs were equal to 0.02, 0.03, and 0.04 within DT of 30–70% for Case Studies 1, 2, and 3, respectively. This reveals that the influence of DT on SF was numerically significant for multilayer soil slopes. One should note that from the practical point of view, these differences are negligible.
• The number of required analyses to evaluate a particular example that considered both design and precision parameters was about 10 using the proposed method. However, it was a considerable number in other methods (depending on the technique).

6. Summary

In this paper, a strength reduction-based technique is presented for slope stability analysis. The key parameters affecting the stability analysis of soil slopes in conventional methods such as LEM and SRM are discussed. The spectral element method as a high-order form of the finite element method was used to attain the slope deformation responses. The soil stress-strain relation was simulated based on the elasto-plastic Mohr–Coulomb failure criterion. An iterative search method was proposed to evaluate the safety factor irrespective of the engineer’s judgment about the limit state criteria.

The proposed search method took into account both the design and precision parameters. It was used to evaluate three case studies, i.e., a homogeneous dry slope, a dry slope with a weak layer, and a partially-wet slope with a weak layer. Soil parameters were considered constant in all of the examples. The advantages of the proposed hybrid method were well-revealed in the simulation of the real-world problems, where it reduced the computational effort compared to the traditional methods.

It was shown that the proposed method converged to the pre-defined precision parameter (i.e., 0.01) with only a few iterations. The computed safety factors were compared with other conventional methods, and an inconsistency in the conventional SRM was found in the estimation of the safety factor. Moreover, the significance of the design trigger for the safety factor of multilayer soil slopes was discussed. Finally, it can be concluded that the proposed method was successful at defining unique limit state criteria for slope stability analysis.

Future work can be directed at the application of the random filed theory in conjunction with the spectral element method in order to compute a more precise safety factor for heterogeneous soil slops. Furthermore, this technique can be used in slope reliability analysis and determination of the failure probability. Last, but not least, this technique can be extended for complex case studies with many degrees of freedom and random variables.