*Article* **Slope Failure Risk Assessment Considering Both the Randomness of Groundwater Level and Soil Shear Strength Parameters**

**Pu Peng <sup>1</sup> , Ze Li 1,\*, Xiaoyan Zhang <sup>1</sup> , Wenlian Liu <sup>2</sup> , Sugang Sui <sup>2</sup> and Hanhua Xu <sup>3</sup>**


**Abstract:** Conducting research on slope failure risk assessment is beneficial for the sustainable development of slopes. There will be various failure modes considering both the randomness of the groundwater level and soil shear strength parameters. Based on the integrated failure probability (*IFP*), the traditional failure risk analysis needs to count all failure modes, including the failure probability (*P<sup>f</sup>* ) and failure risk coefficient (*C*), one-by-one. A new slope failure risk assessment method that uses the sum of the element failure risk to calculate the overall failure risk is proposed in this paper and considers both the randomness of the groundwater level and soil shear strength parameters. The element failure probability is determined by their location information and failure situation; the element failure risk coefficient is determined by their area. It transforms the complex overall failure risk problem into a simple element failure risk problem, which simplifies the calculation process and improves the calculation efficiency greatly. The correctness is verified with the systematic analysis of a classical case. The results show that the slope failure probability and failure risk are greatly increased from 1.40% to 3.30% and 0.829 m<sup>2</sup> to 2.094 m<sup>2</sup> with rising groundwater level, respectively.

**Keywords:** failure risk; element failure probability; spatial variability; stochastic groundwater level; upper bound method

### **1. Introduction**

There are many factors that induce slope instability, among which the influences of the mechanical parameters of soil mass and groundwater are particularly important [1]. Due to long-term geological actions, such as sedimentation, post-sedimentation, chemical weathering, and physical denudation, the mechanical parameters of natural soil show certain spatial variability, which affects both the safety performance and failure mode of a slope. However, these effects are often ignored in practical engineering [2,3]. Groundwater mainly affects the safety performance from two aspects: one is making the seepage field change, affecting the stability performance; and the other is reducing shear strength parameters, affecting the safety performance [4]. Generally, the groundwater level inside the slope is not a certain value and is influenced by various uncertain factors. Therefore, it is essential to establish a failure risk analysis model of soil slopes that considers both the spatial variability and the stochasticity of the groundwater level so as to simulate the slope more in line with the real situation and achieve sustainable development of the slope.

In recent years, the slope reliability problem under consideration of spatial variability has received widespread attention by researchers. Using the limit equilibrium method

**Citation:** Peng, P.; Li, Z.; Zhang, X.; Liu, W.; Sui, S.; Xu, H. Slope Failure Risk Assessment Considering Both the Randomness of Groundwater Level and Soil Shear Strength Parameters. *Sustainability* **2023**, *15*, 7464. https://doi.org/10.3390/ su15097464

Academic Editor: Jianjun Ma

Received: 9 April 2023 Revised: 26 April 2023 Accepted: 28 April 2023 Published: 1 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(LEM), the probability density function (PDF) curve and cumulative probability density function (CDF) curve of the safety factor can be acquired directly by first assuming the sliding surface. However, the LEM only takes the one-dimensional spatial variability of the soil material into account and ignores the impact outside the sliding surface on the slope reliability. In addition, the LEM does not consider the constitutive relation of the soil, so the stress–strain relationship cannot be acquired [5–11]. Using the finite element method (FEM) without assuming the sliding surface, the constitutive relation and spatial variability of materials can be fully considered so as to obtain more reasonable calculation results compared with the LEM [12]. In addition, the FEM can analyze the stress–strain development and progressive slope failure process, but the calculation cost is high [12–16]. Using the limit analysis method (LAM), the range of true solutions and corresponding failure modes can be acquired efficiently and accurately with the upper bound method (UBM) and lower bound method (LBM) [17–24]. Chen Zhaohui et al. combined the stochastic FEM and LAM, taking the spatial variability into account for the slope reliability analysis, and gave the strict frequency distribution range of the safety factor [25]. Peng Pu et al. carried out a study on element failure probability and used affinity propagation (AP) cluster analysis to obtain the failure modes of soil slope [26]. The above research greatly promotes the application prospect of the LAM in slope reliability analysis; however, it is not always reasonable to use only the failure probability to evaluate the slope. For example, there are two failure modes of a slope, corresponding to shallow landslide (*Pf*1*, C*1) and deep landslide (*Pf*2*, C*2). *Pf*<sup>1</sup> is larger than *Pf*2, but *C*<sup>1</sup> is smaller than *C*2, so which failure mode has greater impact? To fully consider the impact of *P<sup>f</sup>* and *C*, researchers proposed the concept of slope failure risk (*R*) [27,28]. However, whether we use the LEM or FEM to calculate the slope failure risk, the failure modes of the slope need to be identified and counted first, and there is lower computational efficiency with large sample calculations. Zhang Xiaoyan et al., using the slope safety factor and velocity field information to calculate the element failure probability (EFP), provide a new idea for the calculation of slope failure risk [29].

Currently, the groundwater level is usually defaulted at a certain value, and the limit state function does not contain the groundwater level as a stochastic variable. However, due to the stochastic distribution of soil particles and pores inside the slope, the seepage field is uncertain. In addition, because of the uncertainty of the supply by precipitation, surface runoff, underground confluence, and the discharge by evaporation and pumping, the groundwater level is uncertain. It is well known that the groundwater level is directly related to the seepage field, so the stochastic groundwater level will lead to the uncertainty of the seepage field. In recent years, the research about the seepage field mainly focused on the change in the permeability coefficient and paid less attention to the stochastic groundwater level, which leads to the change in the seepage field [30,31]. In addition, autocorrelation functions, such as exponential, Gaussian, logarithmic, and triangular, are frequently applied to characterize the spatial variability of soil mass [32]. In view of this, a slope failure risk assessment method that considers both the randomness of the groundwater level and soil shear strength parameters is proposed in this paper to achieve sustainable development.

#### **2. Methodologies**

#### *2.1. Stable Seepage Field of Slope with Stochastic Groundwater Level*

In this paper, a slope failure risk assessment that considers both the randomness of the groundwater level and soil shear strength parameters is carried out. Figure 1 shows the stochastic groundwater level model for soil slopes, which includes the velocity of the triangular element, the model of the groundwater level distribution, the shear strength parameters' calculation, and the pore water pressure calculation. The soil mass is discretized by the triangular element, so the three nodes of the element *e* have horizontal velocity *u e i* vertical velocity *v e i* , cohesive force *c e i* , friction angle *ϕ e i* , and pore water pressure *p e i* , of which velocity *<sup>e</sup>*

pressure *<sup>e</sup>*

④ belong to element *b*.

the pore water pressure is *p h* =

*i* = (1, 2, 3) and *e* = (1, . . . , *Ne*). In addition, there is discontinuity between element *a* and element *b.* Nodes <sup>1</sup> and <sup>2</sup> belong to element *a;* Nodes <sup>2</sup> and <sup>4</sup> belong to element *b*. where *<sup>x</sup> k* and *<sup>y</sup> k* are the permeability coefficient in the *x* and *y* direction, respectively; and *Hw* is the stochastic groundwater head function at each point within the soil.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 3 of 26

γ

, of which

equation is used for the seepage analysis [31]. The specific equation is as follows:

*x y H H k k x y* ∂ ∂

*<sup>i</sup> u* , vertical velocity *<sup>e</sup>*

discretized by the triangular element, so the three nodes of the element *e* have horizontal

between element *a* and element *b.* Nodes ① and ② belong to element *a;* Nodes ② and

follows a truncated normal distribution, where min *H <sup>w</sup>* is the left-truncated tail of the stochastic groundwater level, max *H <sup>w</sup>* is the right-truncated tail of the stochastic groundwater level, and *mean H <sup>w</sup>* is the mean value of the stochastic groundwater level. In addition, this article makes the following assumptions: (1) the seepage field is a saturated and stable seepage field; (2) excess pore water pressure along with the effect on the soil shear strength parameters caused by a sudden rise or drop of groundwater level are not considered; (3) above the phreatic line, the pore water pressure is *p* = 0 ; and (4) below the phreatic line,

Similar to the literature [33], it is assumed that the stochastic groundwater level *Hw*

γ

the vertical distance from the point on the phreatic line that is on the same equipotential line as the calculation point to the calculation point. The two-dimensional stable seepage

2 2

2 2 0 *w w*

+ =

*<sup>i</sup> p* , of which (1, 2, 3) *i* = and (1,..., ) *<sup>e</sup> e N* = . In addition, there is discontinuity

*<sup>i</sup> c* , friction angle *<sup>e</sup>*

ϕ

is the volume weight of water, and *h* is

∂ ∂ (1)

*<sup>i</sup>* , and pore water

*<sup>i</sup> v* , cohesive force *<sup>e</sup>*

**Figure 1.** Stochastic groundwater level model for soil slope. **Figure 1.** Stochastic groundwater level model for soil slope.

*2.2. Stochastic Field of Parameter Spatial Variability*  Considering that slope reliability is sensitive to autocorrelation length and insensitive to the form of the autocorrelation function [34,35], this paper selects the exponential autocorrelation function with a simpler mathematical expression for research, which can be expressed with the following equation [35]: [ ] ( , );( , ) exp 2 2 *aa bb ab ab xy xy h v xx yy L L* ρ − − =− − (2) Similar to the literature [33], it is assumed that the stochastic groundwater level *H<sup>w</sup>* follows a truncated normal distribution, where *Hmin <sup>w</sup>* is the left-truncated tail of the stochastic groundwater level, *Hmax <sup>w</sup>* is the right-truncated tail of the stochastic groundwater level, and *Hmean <sup>w</sup>* is the mean value of the stochastic groundwater level. In addition, this article makes the following assumptions: (1) the seepage field is a saturated and stable seepage field; (2) excess pore water pressure along with the effect on the soil shear strength parameters caused by a sudden rise or drop of groundwater level are not considered; (3) above the phreatic line, the pore water pressure is *p* = 0; and (4) below the phreatic line, the pore water pressure is *p* = *γh*, of which *γ* is the volume weight of water, and *h* is the vertical distance from the point on the phreatic line that is on the same equipotential line as the calculation point to the calculation point. The two-dimensional stable seepage equation is used for the seepage analysis [31]. The specific equation is as follows:

$$k\_x \frac{\partial^2 H\_w}{\partial x^2} + k\_y \frac{\partial^2 H\_w}{\partial y^2} = 0 \tag{1}$$

where *k<sup>x</sup>* and *k<sup>y</sup>* are the permeability coefficient in the *x* and *y* direction, respectively; and *H<sup>w</sup>* is the stochastic groundwater head function at each point within the soil.

### *2.2. Stochastic Field of Parameter Spatial Variability*

Considering that slope reliability is sensitive to autocorrelation length and insensitive to the form of the autocorrelation function [34,35], this paper selects the exponential autocorrelation function with a simpler mathematical expression for research, which can be expressed with the following equation [35]:

$$\rho\_{\left[ (\mathbf{x}\_{a,y\_a}); (\mathbf{x}\_b, y\_b) \right]} = \exp \left[ -2 \left( \frac{\mathbf{x}\_a - \mathbf{x}\_b}{L\_h} \right) - 2 \left( \frac{y\_a - y\_b}{L\_v} \right) \right] \tag{2}$$

where *xa*, *ya*, *x<sup>b</sup>* , and *y<sup>b</sup>* are the coordinate components of the spatial coordinates; and *L<sup>h</sup>* and *L<sup>v</sup>* are the fluctuation range in the *x* and *y* direction, respectively.

Due to the existence of certain mutual correlations between soil parameters and certain autocorrelations of the soil parameters themselves thus involving the discrete process of the associated non-Gaussian fields, this paper adopts a method similar to that of the literature [35] using the midpoint method of Cholesky decomposition for the stochastic field generation, and for the exponential autocorrelation functions, the associated logarithmic stochastic field of the soil materials is as follows:

$$c\_i^{\varepsilon} = \exp(u\_{\ln \varepsilon} + \sigma\_{\ln \varepsilon} \cdot c\_i^{\varepsilon D}) \tag{3}$$

$$\varphi\_{\rm i}^{\varepsilon} = \exp(u\_{\ln \varphi} + \sigma\_{\ln \varphi} \cdot \varphi\_{\rm i}^{\varepsilon D}) \tag{4}$$

where *c e i* and *ϕ e i* are the soil cohesion and internal friction angle of nodes *i* upon element *e* in the non-Gaussian stochastic field, respectively; *c eD i* and *ϕ eD i* are the soil cohesion and internal friction angle of nodes *i* upon element *e* in the Gaussian stochastic field, respectively; *σ*ln *<sup>c</sup>* = p ln(1 + (*σ<sup>c</sup>* <sup>2</sup>/*u<sup>c</sup>* <sup>2</sup>)); *σ*ln *<sup>ϕ</sup>* = q ln(1 + (*σ<sup>ϕ</sup>* <sup>2</sup>/*u<sup>ϕ</sup>* <sup>2</sup>)); *<sup>u</sup>*ln *<sup>c</sup>* = ln *<sup>u</sup><sup>c</sup>* − *<sup>σ</sup>*ln *<sup>c</sup>* <sup>2</sup>/2; *u*ln *<sup>ϕ</sup>* = ln *u<sup>ϕ</sup>* − *σ*ln *<sup>ϕ</sup>* <sup>2</sup>/2; *u<sup>c</sup>* and *σ<sup>c</sup>* are the mean and standard deviation of the log-normal distribution of *c*, respectively; *u<sup>ϕ</sup>* and *σ<sup>ϕ</sup>* are the mean and standard deviation of the lognormal distribution of *ϕ*, respectively; *u*ln *<sup>c</sup>* and *σ*ln *<sup>c</sup>* are the mean and standard deviation of the corresponding normal variable of ln *c*, respectively; and *u*ln *<sup>ϕ</sup>* and *σ*ln *<sup>ϕ</sup>* are the mean and standard deviation of the corresponding normal variable of ln *ϕ*, respectively.

### *2.3. Stochastic Programming Model*

Sloan et al. [20] discretized the soil mass with triangular elements (as shown in Figure 1) and constructed the kinematically admissible velocity fields (KAVF). It is easy to discern from the UBM that the KAVF is one-to-one, corresponding to the external load, where the minimum is infinitely close to the limit load; therefore, the UBM can be understood for solving the minimization problem. On the basis of previous studies [20,26,29,36], the stochastic programming model that considers both the randomness of the groundwater level and soil shear strength parameters established in this paper is as follows:

$$\begin{cases} Z = k\_m - 1\\ Mminimize: \quad k\_\gamma = \mathbf{W}\_\varepsilon + \mathbf{W}\_d - \mathbf{W}\_\varepsilon^p - \mathbf{W}\_d^p\\ Ssubject\ to: \quad \mathbf{a}\_e^1 \mathbf{u}\_e - \mathbf{a}\_e^2 \chi\_\varepsilon = 0\\ \quad \mathbf{a}\_d^1 \mathbf{u}\_d - \mathbf{a}\_d^2 \chi\_d = 0\\ \quad \mathbf{a}\_b \mathbf{u}\_b = 0\\ \mathbf{W}\_\mathbf{G} = 1; \chi\_\varepsilon \ge 0; \chi\_d \ge 0 \end{cases} \tag{5}$$

where *Z* is the limit state function; *k<sup>m</sup>* and *k<sup>γ</sup>* are the safety factor and volume weight overload factor, respectively; *χ<sup>e</sup>* and *χ<sup>d</sup>* are the plastic multiplier of finite element *e* and node *d*, respectively; *a* 1 *e* , *a* 2 *<sup>e</sup>* and *a* 1 *d* , *a* 2 *d* are the plastic flow constraint matrix of the finite element *e* and velocity discontinuity *d*, respectively; **u***<sup>e</sup>* = - *u e* 1 *v e* 1 *u e* 2 *v e* 2 *u e* 3 *v e* 3 *T* , of which *u e i* and *v e i* are the horizontal and vertical velocity components of nodes *i* upon element *e*, respectively; **u***<sup>d</sup>* = - *u d* 1 *v d* 1 *u d* 2 *v d* 2 *u d* 3 *v d* 3 *u d* 4 *v d* 4 *T* , of which *u d i* and *v d i* are the horizontal and vertical velocity components of nodes *i* upon velocity discontinuity *d*, respectively; *a<sup>b</sup>* = cos *θ* sin *θ* − sin *θ* cos *θ* , of which *θ* is the rotation angle; **u***<sup>b</sup>* = h *u b* 1 *v b* <sup>1</sup> K K *u b j v b j* i , of which *u b j* and *v b j* are the horizontal and vertical velocity components on the boundary, respectively; **W***<sup>e</sup>* and **W***<sup>d</sup>* are the internal power of finite element *e* and velocity discontinuity *d*, respectively; **W** *p <sup>e</sup>* and **W** *p d* are the external power of finite element *e* and velocity discontinuity *d*, respectively; and **W**<sup>G</sup> is the external power of self-weight [20].

#### **3. Solution Strategy**

The *a* 2 *e* , *a* 2 *d* , **W***<sup>e</sup>* , **W***<sup>d</sup>* , **W** *p e* , and **W** *p <sup>d</sup>* matrices in Equation (5) are all associated with the groundwater level as well as the shear strength parameters. Currently, there is no better method for solving this kind of problem; in view of this, on the basis of the Monte Carlo simulation, an iterative method is proposed for solving this issue. The detailed implementation process is as follows:

(1) Assuming that the stochastic groundwater level follows a truncated normal distribution, which is generated with the Monte Carlo simulation method:

$$\begin{cases} \begin{array}{l} H\_{\overline{w}}(T\_{\overline{w}}) = \text{Random}(\text{Normal}, \mu\_{\overline{w}}, \sigma\_{\overline{w}}, 1, \text{N}\_{\overline{w}})\\ H\_{\overline{w}}^{\text{min}} \le H\_{\overline{w}}(T\_{\overline{w}}) \le H\_{\overline{w}}^{\text{max}} \end{array} \end{cases} \tag{6}$$

where *Hw*(*Tw*) is the *Tw*th groundwater level; *T<sup>w</sup>* = (1, . . . , *Nw*), *N<sup>w</sup>* is the number of groundwater level; and *µ<sup>w</sup>* and *σ<sup>w</sup>* are the mean and standard deviation of the groundwater level, respectively.

(2) Assuming that the autocorrelation functions of the soil materials are exponential type, using Equations (3) and (4), the midpoint method of Cholesky decomposition for the stochastic field generation *c e* (*Tm*) and *ϕ e* (*Tm*), of which *T<sup>m</sup>* = (1, . . . , *Nm*); *N<sup>m</sup>* is the number of random fields for the shear strength parameters.

In this paper, the slope is allowed to reach the limit state by means of capacitive overload, and the capacitive overload factor of the slope when taking the parameters' spatial variability into account under the effect of the stochastic groundwater level is as follows:

$$k\_{\gamma}(T\_{\mathfrak{m}}, T\_{\mathfrak{w}}) = \frac{\gamma\_{\mathfrak{c}}(\mathfrak{c}^{\mathfrak{c}}(T\_{\mathfrak{m}}), \mathfrak{q}^{\mathfrak{c}}(T\_{\mathfrak{m}}), H\_{\mathfrak{w}}(T\_{\mathfrak{w}}))}{\gamma\_{\mathfrak{a}}} \tag{7}$$

where *kγ*(*Tm*, *Tw*) is the capacitive overload factor; *γc*(*c e* (*Tm*), *ϕ e* (*Tm*), *Hw*(*Tw*)) is the capacity of the soil in the ultimate state, which is related to *c e* (*Tm*), *ϕ e* (*Tm*), and *Hw*(*Tw*); and *γ<sup>a</sup>* is the actual capacity of the soil.

The slope safety factor when taking the parameters' spatial variability into account under the effects of the stochastic groundwater level is as follows:

$$k\_m(T\_m, T\_w) = \frac{c^\varepsilon(T\_m)}{c^{\prime \varepsilon}(T\_m)} = \frac{q^\varepsilon(T\_m)}{q^{\prime \varepsilon}(T\_m)}\tag{8}$$

where *km*(*Tm*, *Tw*) is the safety factor; and *c* 0*e* (*Tm*) and *ϕ* 0*e* (*Tm*) are the cohesion and internal friction angle of finite element *e* in the *Tm*th non-Gaussian stochastic field after the strength reduction, respectively.


*T N w w* = ; ( ) *<sup>e</sup>*

*<sup>m</sup> c T* , ( ) *<sup>e</sup>* ϕ

*N N w m* × numbers of capacity overload factors [ ( , )] *m w kTT*

(5) Calculation of the slope safety factor and plot the related curve.

**Figure 2.** Specific numerical solution flow for slope. **Figure 2.** Specific numerical solution flow for slope.

#### **4. Reliability Index for Slope**

**4. Reliability Index for Slope**  The traditional *IFP* uses the threshold value of the safety factor to determine whether the slope will fail or not and then conducts the slope reliability analysis accordingly [28]. The traditional *IFP* uses the threshold value of the safety factor to determine whether the slope will fail or not and then conducts the slope reliability analysis accordingly [28]. Based on the *IFP*, the slope failure function is as follows:

$$I(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) = \begin{cases} \ \ 0, \text{ if } k\_{\mathfrak{m}\prime}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) \ge 1\\ \ \ 1, \text{ if } k\_{\mathfrak{m}\prime}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) < 1 \end{cases} \tag{9}$$

*Tm* from 1 *Tm* = to *T N m m* = cycles, the number of *Nm* sto-

γ

while, at the same time,

chastic fields are brought into Equation (5) and use the dual simplex method to obtain

use the bisection method to obtain *N N w m* × numbers of the slope safety factor

0, if ( , ) 1 (,) 1, if ( , ) 1 *m w mmw IT T kTT* <sup>≥</sup> <sup>=</sup> <sup>&</sup>lt; (9) where *I*(*Tm*, *Tw*) is the failure function of the slope corresponding to the *Tm*th stochastic field under *Tw*th groundwater level acts.

*mmw*

where (,) *m w IT T* is the failure function of the slope corresponding to the *Tm* th stochastic field under *Tw* th groundwater level acts. According to the Equation (8) slope failure function, the IFP under *Tw*th groundwater level acts can be acquired as follows:

$$P\_f^{IFP}(T\_w) = \frac{\sum\_{w=1}^{N\_w} I(T\_{W\_1}, T\_w)}{N\_m} \times 100\% \tag{10}$$

(,) *m N m w IT T* (10) where *P IFP f* (*Tw*) is the *IFP* of the slope under *Tw*th groundwater level acts.

1 ( ) 100% *m IFP T P T* <sup>=</sup> = × Further, the *IFP* of the slope is as follows:

*f w*

$$P\_f^{IFP} = \sum\_{T\_W=1}^{N\_W} P\_f^{IFP}(T\_W) \tag{11}$$

Further, the *IFP* of the slope is as follows: where *P IFP f* is the *IFP* of the slope under all potential groundwater level acts.

The integrated failure risk (*IFR*) of the slope under *Tw*th groundwater level acts is as follows:

$$R^{IFP}(T\_w) = P\_f^{IFP}(T\_w) \mathbb{C} \tag{12}$$

where *R IFP*(*Tw*) is the *IFR* of the slope under *Tw*th groundwater level acts, and *C* is the area of the slope failure mode.

Based on the *IFP*, the *IFR* of the slope is calculated as follows [27,28]:

$$R^{IFP} = \sum\_{T\_W=1}^{N\_w} R^{IFP}(T\_W) \tag{13}$$

where *R IFP* is the *IFR* of the slope.

The Equation (13) default is that only one failure mode exists when conducting the slope failure risk analysis, which is contrary to the fact that slopes have multiple failure modes, considering both the randomness of the groundwater level and soil shear strength parameters. For that reason, Huang et al. [37] gave a calculation equation, which considered all failure modes of the slope simultaneously as follows:

$$R\_k^{IFP}(T\_w) = P\_{fk}^{IFP}(T\_w) \mathbb{C}\_k \tag{14}$$

where *R IFP k* (*Tw*) is the *IFR* of the *k*th failure mode of the slope under *Tw*th groundwater level acts; *P IFP f k* (*tw*) is the *IFP* of the *k*th failure mode of the slope under *Tw*th groundwater level acts; and *C<sup>k</sup>* is the area of the *k*th slope failure mode.

The IFR of the slope is as follows:

$$\mathcal{R}^{IFP} = \sum\_{T\_w=1}^{N\_w} \sum\_{k=1}^{N} \mathcal{R}\_k^{IFP}(T\_w) \tag{15}$$

Compared with Equations (13) and (15), it can perform both single and multiple failure mode slope risk assessments; however, both the LEM and FEM require prior work on complex failure mode classification. Zhang Xiaoyan et al. [29] proposed a new concept of *EFP* on the basis of the UBM. They used the slope safety factor and the KVAF to judge the failure elements. The failure function, considering the stochasticity of the groundwater level and spatial variability of soil material, is as follows:

$$I\_{\varepsilon}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) = \begin{cases} \begin{array}{l} 0 \text{ if } k\_{\mathrm{m}}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) \ge 1 \\ 0 \text{ if } \mu^{\varepsilon}\_{\mathrm{c}}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) = 0 \text{ and } k\_{\mathrm{m}}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) < 1 \\ 1 \text{ if } \mu^{\varepsilon}\_{\mathrm{c}}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) > 0 \text{ and } k\_{\mathrm{m}}(T\_{\mathfrak{m}\prime}T\_{\mathfrak{w}}) < 1 \end{array} \end{cases} \tag{16}$$

where *Ie*(*Tm*, *Tw*) is the failure function of the element *e* corresponding to the *Tm*th stochastic field under *Tw*th groundwater level acts, and *u e* c (*Tm*, *Tw*) is the velocity of the element *e* corresponding to the *Tm*th stochastic field under *Tw*th groundwater level acts.

$$u\_{\mathbf{c}}^{\varepsilon}(T\_{\mathbf{m}\prime}T\_{\mathbf{w}}) = \sqrt{\left(\sum\_{i=1}^{3} \frac{u\_i^{\varepsilon}(T\_{\mathbf{m}\prime}T\_{\mathbf{w}})}{\mathfrak{3}}\right)^2 + \left(\sum\_{i=1}^{3} \frac{v\_i^{\varepsilon}(T\_{\mathbf{m}\prime}T\_{\mathbf{w}})}{\mathfrak{3}}\right)^2} \tag{17}$$

where *u e i* (*Tm*, *Tw*) is the horizontal velocity component of nodes *i* upon element *e* corresponding to the *Tm*th stochastic field under *Tw*th groundwater level acts, and *v e i* (*Tm*, *Tw*) is the vertical velocity component of nodes *i* upon element *e* corresponding to the *Tm*th stochastic field under *Tw*th groundwater level acts.

According to the slope element failure function in Equation (16), the slope element failure probability (*EFP*) under *Tw*th groundwater level acts is as follows:

$$P\_{fc}^{EFP}(T\_w) = \frac{\sum\_{m=1}^{N\_m} I^c(T\_{m\prime}, T\_w)}{N\_m} \times 100\% \tag{18}$$

where *P EFP f e* (*Tw*) is the *EFP* under *Tw*th groundwater level acts.

The slope *EFP* under all potential groundwater level acts is as follows:

$$P\_{fc}^{EFP} = \frac{\sum\_{w=1}^{N\_w} \sum\_{T\_m=1}^{N\_m} I^e(T\_{m\prime}, T\_w)}{N\_w \times N\_m} \times 100\% \tag{19}$$

The element failure risk (*EFR*) of the slope under *Tw*th groundwater level acts is as follows:

$$R\_{\varepsilon}^{EFP}(T\_w) = P\_{f\varepsilon}^{EFP}(T\_w) \mathcal{C}\_{\varepsilon} \tag{20}$$

where *R EFP e* (*Tw*) is the failure risk of element *e* under *Tw*th groundwater level acts, and *C<sup>e</sup>* is the area of element *e*.

The element failure risk (*EFR*) of the slope under all potential groundwater level acts is as follows:

$$R\_{\varepsilon}^{EFP} = \sum\_{T\_W=1}^{N\_W} R\_{\varepsilon}^{EFP}(Tw) \tag{21}$$

where *R EFP e* is the failure risk of element *e* under all potential groundwater level acts.

The *IFR* based on the *EFP* of the slope under *Tw*th groundwater level acts is calculated in the following equation [20]:

$$\mathcal{R}^{EFP}(T\_w) = \sum\_{\varepsilon=1}^{N\_\varepsilon} \mathcal{R}^{EFP}\_\varepsilon(T\_w) \tag{22}$$

where *R EFP*(*Tw*) is the *IFR* under *Tw*th groundwater level acts.

The *IFR* based on the *EFP* of the slope under all potential groundwater level acts is calculated in the following equation:

$$R^{EFP} = \sum\_{T\_{\mathcal{W}}=1}^{N\_{\mathcal{W}}} R^{EFP}(T\_{\mathcal{W}}) \tag{23}$$

where *R EFP* is the *IFR* under all potential groundwater level acts.

A new slope failure risk assessment under consideration of the stochastic groundwater level involves *EFP* and *C<sup>e</sup>* . The *EFP* can be easily obtained by solving the stochastic programming model, and *C<sup>e</sup>* is the area of element *e*, which is constant compared to *C<sup>k</sup>* .

According to the solution strategy of the stochastic programming model, the *N<sup>w</sup>* × *N<sup>m</sup>* slope safety factors *km*(*Tm*, *Tw*) can be easily acquired. Using statistical knowledge, the mean and standard deviation of the slope safety factor can be calculated as follows:

$$\mu\_k(T\_w) = \frac{\sum\_{m=1}^{N\_m} k\_m(T\_{m\nu}T\_w)}{N\_m} \tag{24}$$

$$\mu\_k = \frac{\sum\_{w=1}^{N\_w} \sum\_{m=1}^{N\_m} k\_m (T\_{m\nu} T\_{w})}{(N\_w \times N\_m)} \tag{25}$$

$$\sigma\_k(T\_w) = \sqrt{\frac{\sum\_{T\_w=1}^{N\_m} \left(k\_m(T\_{\rm inv}, T\_w) - \mu\_k(T\_w)\right)}{N\_m - 1}}\tag{26}$$

μ

−

$$\sigma\_k = \sqrt{\frac{\sum\_{\text{w}}^{N\_{\text{w}}} \sum\_{\text{n}=1}^{N\_{\text{w}}} \left(k\_{\text{w}}(T\_{\text{w}}, T\_{\text{w}}) - \mu\_k\right)}{N\_{\text{w}} \times N\_{\text{m}} - 1}} \tag{27}$$

where *T<sup>w</sup>* = (1, . . . , *Nw*); *µ<sup>k</sup>* (*Tw*) and *σ<sup>k</sup>* (*Tw*) are the mean and standard deviation of the slope safety factor under *Tw*th groundwater level acts, respectively; and *µ<sup>k</sup>* and *σ<sup>k</sup>* are the mean and standard deviation of the slope safety factor under all potential groundwater level acts, respectively. **5. Calibration and Application**  The UBM program is compiled, and a classic slope calculation example is calculated and analyzed. Comparing the result with the calculation result of the LEM, we verified

1

*m*

1 1

*w m*

*T T*

= = <sup>=</sup> <sup>×</sup>

1

1 1

<sup>=</sup> × −

*w m*

*T T*

= =

*w m*

*N N*

*m*

*T*

*m*

*N*

( ) <sup>1</sup>

<sup>=</sup> <sup>−</sup>

*w m*

*N N*

*T*

*m*

*N*

( )

<sup>=</sup> =

*T*

*k w*

*k*

μ

*k w*

σ

*T*

*k*

σ

<sup>=</sup>

μ

(,)

(,)

( )

*kTT T*

*mmw kw*

( )

*mmw k*

1

(,)

*m*

*N*

*kTT*

*w m*

*N N*

(,) ()

−

μ

*mmw*

*kTT*

(24)

(25)

(26)

(27)

*<sup>k</sup>* are

*mmw*

*kTT*

*m*

*N*

( )

*N N*

*w m*

#### **5. Calibration and Application** the correctness of the calculation method.

μ

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 9 of 26

The UBM program is compiled, and a classic slope calculation example is calculated and analyzed. Comparing the result with the calculation result of the LEM, we verified the correctness of the calculation method. *5.1. Numerical Simulations*  Figure 3 shows the homogeneous slope calculation model. The height is 5 m, the width

#### *5.1. Numerical Simulations* of the top is 10 m, and the ratio is 1:2. Li Dian qing et al. have calculated the failure prob-

where (1,..., ) *T N w w* = ; ( )

ter level acts, respectively.

Figure 3 shows the homogeneous slope calculation model. The height is 5 m, the width of the top is 10 m, and the ratio is 1:2. Li Dian qing et al. have calculated the failure probability of the slope [35] but have not conducted systematic research on the slope failure risk. In view of this, on the basis of the method proposed in this paper, the slope failure risk is studied. Using the triangular element to discretize the slope, 989 finite elements, 2967 nodes, and 1536 discontinuities are acquired. In addition, there are three pore water pressure monitoring points: P1 (5,5), P2 (10,5), and P3 (15,5). The set soil volume weight *γ* is 20.0 kN/m<sup>3</sup> , and the permeability coefficient is *<sup>K</sup>* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−7m/s, both of which are determined values. See Table 1 for the other calculation parameters [26]. ability of the slope [35] but have not conducted systematic research on the slope failure risk. In view of this, on the basis of the method proposed in this paper, the slope failure risk is studied. Using the triangular element to discretize the slope, 989 finite elements, 2967 nodes, and 1536 discontinuities are acquired. In addition, there are three pore water pressure monitoring points: P1 (5,5), P2 (10,5), and P3 (15,5). The set soil volume weight *γ* is 20.0 kN/m3, and the permeability coefficient is <sup>7</sup> *K*=5 10 m/s <sup>−</sup> × , both of which are determined values. See Table 1 for the other calculation parameters [26].

**Figure 3.** Homogeneous slope calculation model. **Figure 3.** Homogeneous slope calculation model.

**Table 1.** Statistics of the soil parameters.


The mean and standard deviation of the stochastic groundwater level is 7.5 m and 2.25, respectively, on the basis of the measured groundwater in many projects. The upper boundary of the groundwater level is 8.5 m. The lower boundary of the groundwater level is 6.5 m. The quantity of the stochastic groundwater levels is 50. The groundwater level at

the right slope toe is 5 m. According to Equation (6), 50 stochastic numbers of groundwater levels are generated, and their distribution is shown in Figure 4. The stochastic number distribution of the groundwater is relatively close to the mean value of the groundwater level, and relatively small near the relatively high and low groundwater levels. the right slope toe is 5 m. According to Equation (6), 50 stochastic numbers of groundwater levels are generated, and their distribution is shown in Figure 4. The stochastic number distribution of the groundwater is relatively close to the mean value of the groundwater level, and relatively small near the relatively high and low groundwater levels. is 6.5 m. The quantity of the stochastic groundwater levels is 50. The groundwater level at the right slope toe is 5 m. According to Equation (6), 50 stochastic numbers of groundwater levels are generated, and their distribution is shown in Figure 4. The stochastic number distribution of the groundwater is relatively close to the mean value of the groundwater

The mean and standard deviation of the stochastic groundwater level is 7.5 m and

2.25, respectively, on the basis of the measured groundwater in many projects. The upper

is 6.5 m. The quantity of the stochastic groundwater levels is 50. The groundwater level at

*Lv* <sup>=</sup> 4 m *ρc,<sup>φ</sup>* <sup>=</sup> <sup>−</sup>0.5

*Lv* <sup>=</sup> 4 m *ρc,<sup>φ</sup>* <sup>=</sup> <sup>−</sup>0.5

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 10 of 26

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 10 of 26

*φ*(°) 30 0.2 Lognormal

*φ*(°) 30 0.2 Lognormal

*c*(kPa) 10 0.3 Lognormal *Lh* = 40 m

level, and relatively small near the relatively high and low groundwater levels.

*c*(kPa) 10 0.3 Lognormal *Lh* = 40 m

2.25, respectively, on the basis of the measured groundwater in many projects. The upper boundary of the groundwater level is 8.5 m. The lower boundary of the groundwater level

**Figure 4.** Stochastic groundwater levels distribution histogram. **Figure 4.** Stochastic groundwater levels distribution histogram. Figure 5a–c shows the slope stable seepage fields when =1 *Tw* (*Hw* = 6.8006 m),

Figure 5a–c shows the slope stable seepage fields when =1 *Tw* (*Hw* = 6.8006 m), =25 *Tw* (*Hw* = 7.4848 m), and =50 *Tw* (*Hw* = 8.1951 m), respectively. It can be observed that the contours of the pore water pressure become steeper, and the saturated area inside the slope increases as the groundwater level rises. Figure 6 shows the key points pore water pressure. Under the same groundwater level act, the pore water pressure decreases gradually as the coordinates of the key points move to the right. At P1, the mean and standard deviation of the pore water pressure are −19.77 kPa and 2.31 kPa, respectively. At P2, the Figure 5a–c shows the slope stable seepage fields when *T<sup>w</sup>* = 1 (*H<sup>w</sup>* = 6.8006 m), *T<sup>w</sup>* = 25 (*H<sup>w</sup>* = 7.4848 m), and *T<sup>w</sup>* = 50 (*H<sup>w</sup>* = 8.1951 m), respectively. It can be observed that the contours of the pore water pressure become steeper, and the saturated area inside the slope increases as the groundwater level rises. Figure 6 shows the key points pore water pressure. Under the same groundwater level act, the pore water pressure decreases gradually as the coordinates of the key points move to the right. At P1, the mean and standard deviation of the pore water pressure are −19.77 kPa and 2.31 kPa, respectively. At P2, the mean and standard deviation of the pore water pressure are −14.59 kPa and 1.78 kPa, respectively. At P3, the mean and standard deviation of the pore water pressure are −8.94 kPa and 1.13 kPa, respectively. =25 *Tw* (*Hw* = 7.4848 m), and =50 *Tw* (*Hw* = 8.1951 m), respectively. It can be observed that the contours of the pore water pressure become steeper, and the saturated area inside the slope increases as the groundwater level rises. Figure 6 shows the key points pore water pressure. Under the same groundwater level act, the pore water pressure decreases gradually as the coordinates of the key points move to the right. At P1, the mean and standard deviation of the pore water pressure are −19.77 kPa and 2.31 kPa, respectively. At P2, the mean and standard deviation of the pore water pressure are −14.59 kPa and 1.78 kPa, respectively. At P3, the mean and standard deviation of the pore water pressure are −8.94 kPa and 1.13 kPa, respectively.

**Figure 5.** *Cont*.

**Figure 5.** Steady seepage **Figure 5.** fields of slope: ( Steady seepage fields of slope: ( **a**) *Tw* = 1; (**b**) *Tw***a**= 25; (**c**) *Tw* = 50. ) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50.

**Figure 6.** The key points pore water pressure.

**Figure 6.** The key points pore water pressure.

relation with internal friction angle in space.

Assuming that the autocorrelation functions of the soil materials are exponential type, we used the midpoint method of Cholesky decomposition for the stochastic field generation. The quantity of the stochastic fields of the soil parameters is 2000; according to Equations (3) and (4), 2000 stochastic fields of the slope shear strength parameters are generated. Figure 7a,b shows the stochastic fields of *ce*(500) and *ϕe*(500), respectively. The

generated. Figure 7a,b shows the stochastic fields of (500) *<sup>e</sup> c* and (500)

Assuming that the autocorrelation functions of the soil materials are exponential type, we used the midpoint method of Cholesky decomposition for the stochastic field generation. The quantity of the stochastic fields of the soil parameters is 2000; according to Equations (3) and (4), 2000 stochastic fields of the slope shear strength parameters are

The maximum value of soil cohesion and internal friction angle are 18.9084 kPa and 23.3219°, respectively; The minimum value of soil cohesion and internal friction angle are 3.1696 kPa and9.5851°, respectively. In addition, soil cohesion has a certain negative cor-

ϕ

*<sup>e</sup>* , respectively.

maximum value of soil cohesion and internal friction angle are 18.9084 kPa and 23.3219◦ , respectively; The minimum value of soil cohesion and internal friction angle are 3.1696 kPa and 9.5851◦ , respectively. In addition, soil cohesion has a certain negative correlation with internal friction angle in space. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 13 of 26

**Figure 7.** Stochastic fields of the slope shear strength parameters: (**a**) *cr* (500); (**b**) *φ<sup>r</sup>* (500). **Figure 7.** Stochastic fields of the slope shear strength parameters: (**a**) *c r* (500); (**b**) *ϕ r* (500).

LEM 1.2622 0.1461 2.70

LEM 1.2310 0.1432 4.65

#### *5.2. Main Results of the Research 5.2. Main Results of the Research*

The calculation results show the following:

This paper selects low, medium, and high groundwater levels ( =1 *Tw* , =25 *Tw* , and =50 *Tw* ) for comparative analysis to verify the effectiveness of the calculation method. Table 2 shows the calculation results of the reliability index comparison between the UBM and LEM under three groundwater level acts: =1 *Tw* , =25 *Tw* , and =50 *Tw* . Figure 8 shows the distribution characteristics of the slope safety factor under three groundwater level acts. Figure 9 shows the distribution histograms of the slope safety factor under three This paper selects low, medium, and high groundwater levels (*T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50) for comparative analysis to verify the effectiveness of the calculation method. Table 2 shows the calculation results of the reliability index comparison between the UBM and LEM under three groundwater level acts: *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50. Figure 8 shows the distribution characteristics of the slope safety factor under three groundwater level acts. Figure 9 shows the distribution histograms of the slope safety factor under three groundwater level acts.


groundwater level acts. **Table 2.** Calculation results of the reliability index.

is larger than that of the LEM, but the error is small, which conforms to the features of the upper bound solution. In addition, the slope safety factors acquired with the UBM and LEM methods decrease as the groundwater level rises. The slope failure probability increases acquired with the UBM and LEM decrease as the groundwater level rises. On the basis of the upper bound theorem, the slope safety factor acquired

*Tw* <sup>=</sup> <sup>50</sup>UBM 1.2357 0.1365 3.30

with the UBM must be greater than the real solution. Therefore, the UBM will slightly

(2) Figure 8a,b shows the PDF and CDF curves of the slope safety factors acquired with the UBM and LEM under three groundwater level acts, =1 *Tw* , =25 *Tw* , and =50 *Tw* , respectively. It is not difficult to see that the PDF and CDF curves of the slope safety factors acquired with the UBM and LEM are very close with small errors. In addition, the PDF and CDF curves gradually move to the left as the groundwater level rises. (3) Figure 9a–c are the distribution histograms of the slope safety factors acquired with the UBM under the three groundwater level acts, =1 *Tw* , =25 *Tw* , and =50 *Tw* , respectively. It is not difficult to see that the distribution of the slope safety factors is similar

underestimate the failure slope probability.

to the stochastic groundwater levels.

**Figure 8.** Distribution characteristics of the slope safety factors: (**a**) PDF curve of the safety factor; (**b**) CDF curve of the safety factor. **Figure 8.** Distribution characteristics of the slope safety factors: (**a**) PDF curve of the safety factor; (**b**) CDF curve of the safety factor.

The calculation results show the following:


**Figure 9.** Distribution histograms of the slope safety factors: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50. **Figure 9.** Distribution histograms of the slope safety factors: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50.

On the basis of the stochastic programming model and solution strategy, the quantity of 100,000 slope safety factors was acquired. Figures 10 and 11 are the number of 50 PDF and CDF curves of the slope safety factors, respectively. Figure 12 reflects the mean and standard deviation of the slope safety factors versus the groundwater level. Figure 13 On the basis of the stochastic programming model and solution strategy, the quantity of 100,000 slope safety factors was acquired. Figures 10 and 11 are the number of 50 PDF and CDF curves of the slope safety factors, respectively. Figure 12 reflects the mean and standard deviation of the slope safety factors versus the groundwater level. Figure 13

shows the PDF and CDF curves of the slope safety factors under all potential groundwater level acts acquired. The following rules can be acquired through analysis: *Sustainability* **2023**, *15*, x FOR PEER REVIEW 16 of 26


$$
\mu\_k(T\_w) = -0.0439 H\_w + 1.5963\tag{28}
$$

$$
\sigma\_k(T\_w) = -0.0039 H\_w + 0.1685\tag{29}
$$

(3) The quantity of 100,000 slope safety factors was acquired from 2000 stochastic fields under 50 groundwater level acts to perform the statistical analysis of all the acquired data; under 50 groundwater level acts, the mean and standard deviation of the slope safety factors are 1.2664 and 0.11, respectively. *kw w T H* =− + (29) (3) The quantity of 100,000 slope safety factors was acquired from 2000 stochastic fields under 50 groundwater level acts to perform the statistical analysis of all the acquired data; under 50 groundwater level acts, the mean and standard deviation of the slope

( ) 0.0039 0.1685

σ

safety factors are 1.2664 and 0.11, respectively.

**Figure 10.** PDF curve of the slope safety factors. **Figure 10.** PDF curve of the slope safety factors.

Based on 100,000 Monte Carlo simulations of the slope, the LEM default is that only one failure mode exists (as shown in the LEM sliding surface in Figure 14). The biggest advantage of the UBM compared with the LEM in this paper is that all slope failure modes can be captured according to the failure information of the elements (as shown in the pink failure area in Figure 14) and then the reliability index for the slope can be calculated. Similar to the method in reference 37, the failure areas were used to classify the failure modes; there are six failure modes of the slope under 50 groundwater level acts (as shown in Figure 14). Table 3 lists the failure risk of the above six slope failure modes. Table 4 lists the failure risk corresponding to the failure modes when *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50. The results show that failure modes 1 and 2 have a small failure area that belongs to a shallow landslide, and the failure area is between 28.68 and 58.91 m<sup>2</sup> (failure times is 997). Failure modes 5 and 6 have a large failure area that belongs to a deep landslide, and the failure area is between 89.17 and 119.38 m<sup>2</sup> (failure times is 160). Failure modes 3 and 4 are between a shallow landslide and a deep landslide, and the failure area is between 58.91 and 89.04 m<sup>2</sup> (failure times is 989). When the groundwater level is *T<sup>w</sup>* = 1, the slope has only

five failure modes; when the groundwater level is *T<sup>w</sup>* = 25 and *T<sup>w</sup>* = 50, the fifth failure modes occur. The phreatic line moves up, correspondingly, and the saturated area inside the slope increases as the groundwater level rises, thus increasing the probability of the failure mode. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 17 of 26 *Sustainability* **2023**, *15*, x FOR PEER REVIEW 17 of 26

**Figure 11.** CDF curve of the slope safety factors. **Figure 11.** CDF curve of the slope safety factors. **Figure 11.** CDF curve of the slope safety factors.

**Figure 12. Figure 12.** Slope safety factors versus groundwater level: ( Slope safety factors versus groundwater level: ( **<sup>a</sup>**) *uk* (*T***a***w*); (**b**) *σ<sup>k</sup>* (*Tw*). ) *u<sup>k</sup>* (*Tw*); (**b**) *σ<sup>k</sup>* (*Tw*).

acts.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 18 of 26

**Figure 13.** PDF and CDF curves of the slope safety factors under all potential groundwater level **Figure 13.** PDF and CDF curves of the slope safety factors under all potential groundwater level acts. failure mode.

**Figure 13.** PDF and CDF curves of the slope safety factors under all potential groundwater level

Based on 100,000 Monte Carlo simulations of the slope, the LEM default is that only

one failure mode exists (as shown in the LEM sliding surface in Figure 14). The biggest

landslide, and the failure area is between 28.68 and 58.91 m2 (failure times is 997). Failure modes 5 and 6 have a large failure area that belongs to a deep landslide, and the failure area is between 89.17 and 119.38 m2 (failure times is 160). Failure modes 3 and 4 are be-**Figure 14.** Schematic diagram results of the slope failure modes: (**a**) Failure mode 1; (**b**) Failure mode 2; (**c**) Failure mode 3; (**d**) Failure mode 4; (**e**) Failure mode 5; (**f**) Failure mode 6. **Figure 14.** Schematic diagram results of the slope failure modes: (**a**) Failure mode 1; (**b**) Failure mode 2; (**c**) Failure mode 3; (**d**) Failure mode 4; (**e**) Failure mode 5; (**f**) Failure mode 6.


tween a shallow landslide and a deep landslide, and the failure area is between 58.91 and **Table 3.** Failure risk of the slope.

> The failure mode acquired with the LEM is only consistent with failure mode 3 in this paper. The main reason is that, when the LEM is adopted to calculate the slope safety factor, the initial slip surface is assumed in advance, then a constantly repeated search based on the mean of the shear parameters to acquire the critical slip crack surface is performed, and

**Figure 14.** Schematic diagram results of the slope failure modes: (**a**) Failure mode 1; (**b**) Failure mode

2; (**c**) Failure mode 3; (**d**) Failure mode 4; (**e**) Failure mode 5; (**f**) Failure mode 6.

then the slope safety factor is calculated based on the unique critical slip surface. Thus, a difference exists between the critical slip surface obtained with the LEM and the actual slip surface. The UBM constructs a stochastic programming model through finite elements, and then uses the method of stochastic mathematical programming to search the instability area of the slope, so the result is more in line with the real situation. The calculation indicates that the traditional LEM will ignore some failure modes and may miscalculate the slope failure risk, considering both the randomness of the groundwater level and soil shear strength parameters. The UBM can ignore the constitutive relation of the materials and acquire the slope safety factor and failure modes. Its applicability and calculation efficiency are both high.


**Table 4.** Failure risk of the slope when *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50.

Table 5 is the statistical table of the slope failure risk acquired with the three methods. The LEM default is that only one failure mode exists when *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50; the slope failure risk according to the LEM with Equation (12) is 1.121 m<sup>2</sup> , 1.593 m<sup>2</sup> , and 2.325 m<sup>2</sup> , respectively. All failure modes can be acquired with the UBM. When *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50, the slope failure risk according to the UBM with Equation (14) and the UBM Equation (22) are 0.829 m<sup>2</sup> , 1.302 m<sup>2</sup> , and 2.094 m<sup>2</sup> , respectively. It should be noted that Equation (14) has a difference calculation principle from Equation (22). All slope failure modes are required to be counted when Equation (14) is used to calculate the slope failure risk; the EFP for all elements is easy to acquire by solving Equation (5), and the element area is fixed when using Equation (22) to calculate the slope failure risk. From the calculation principle, the proposed method will simplify the calculation process and make the calculation more efficient.

**Table 5.** Failure risk statistical table of the slope (Unit: m<sup>2</sup> ).


 **Method** 

calculation more efficient.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 20 of 26

**Table 5.** Failure risk statistical table of the slope (Unit: m2).

**Groundwater LEM with Equation (12) UBM with Equation (14) UBM with Equation (22)** 

*Tw* = 1 1.121 0.829 0.829 *Tw* = 25 1.593 1.302 1.302

Equation (14) has a difference calculation principle from Equation (22). All slope failure modes are required to be counted when Equation (14) is used to calculate the slope failure risk; the EFP for all elements is easy to acquire by solving Equation (5), and the element area is fixed when using Equation (22) to calculate the slope failure risk. From the calculation principle, the proposed method will simplify the calculation process and make the

The EFP of the slope under *Tw*th groundwater level acts when *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50 obtained with Equation (18) is shown in Figure 15. The EFP of the slope under all potential groundwater level acts obtained with Equation (19) is shown in Figure 16. It can be observed that the groundwater level influences the EFP of the slope. The white area in the figure is the non-failure element, and the blue area is the failure element. In addition, according to the theory of EFP, the darker the color, the greater the EFP. The EFP of the slope under *Tw* th groundwater level acts when =1 *Tw* , =25 *Tw* , and =50 *Tw* obtained with Equation (18) is shown in Figure 15. The EFP of the slope under all potential groundwater level acts obtained with Equation (19) is shown in Figure 16. It can be observed that the groundwater level influences the EFP of the slope. The white area in the figure is the non-failure element, and the blue area is the failure element. In addition, according to the theory of EFP, the darker the color, the greater the EFP.

**Figure 15.** EFP of the slope under *Tw*th groundwater level acts: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50. **Figure 15.** EFP of the slope under *Tw*th groundwater level acts: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50.

The EFR of the slope under *Tw*th groundwater level acts when *T<sup>w</sup>* = 1, *T<sup>w</sup>* = 25, and *T<sup>w</sup>* = 50 obtained with Equation (20) is shown in Figure 17. The EFR of the slope under all potential groundwater level acts obtained with Equation (21) is shown in Figure 18. It

The EFR of the slope under *Tw* th groundwater level acts when =1 *Tw* , =25 *Tw* , and

=50 *Tw* obtained with Equation (20) is shown in Figure 17. The EFR of the slope under all potential groundwater level acts obtained with Equation (21) is shown in Figure 18. It can be observed that the groundwater level influences the EFR of the slope. When =1 *Tw* , the frequency of EFR between 0 and 0.00062 m2 is 675, the frequency of EFR between 0.00062 and 0.00124 m2 is 36, the frequency of EFR between 0.00124 and 0.00186 m2 is 58, the frequency of EFR between 0.00186 and 0.00248 m2 is 69, the frequency of EFR between 0.00248 and 0.00310 m2 is 140, the frequency of EFR between 0.00310 and 0.00372 m2 is 11, and the slope failure risk is 0.829 m2. When =25 *Tw* , the frequency of EFR between 0 and 0.00082 m2 is 670, the frequency of EFR between 0.00082 and 0.00164 m2 is 32, the frequency of EFR between 0.00164 and 0.00246 m2 is 54, the frequency of EFR between 0.00246 and 0.00328 m2 is 56, the frequency of EFR between 0.00328 and 0.00410 m2 is 165, the frequency of EFR between 0.00410 and 0.00492 m2 is 12, and the slope failure risk is 1.302 m2. When =50 *Tw* , the frequency of EFR between 0 and 0.00150 m2 is 666, the frequency of EFR between 0.00150 and 0.00300 m2 is 41, the frequency of EFR between 0.00300 and 0.00450 m2 is 48, the frequency of EFR between 0.00450 and 0.00600 m2 is 62, the frequency of EFR between 0.00600 and 0.00750 m2 is 158, and the frequency of EFR between 0.00750 and 0.00900 m2 is 14. Under all potential groundwater level acts, the frequency of EFR between 0 and 0.04820 m2 is 674, the frequency of EFR between 0.04820 and 0.09640 m2 is 30, and the frequency of EFR between 0.09640 and 0.14460 m2 is 57. The frequency of EFR between 0.014460 and 0.19280 m2 is 59, the frequency of EFR between 0.19280 and 0.24100 m2 is 157, the frequency of EFR between 0.24100 and 0.28892 m2 is 12, and the slope failure risk is 2.094 m2. The element failure risk comprehensively reflects the contribution of the EFP and the element failure risk coefficient, which can make a quantitative judgment on the slope failure risk of each part. The slope failure risk is 1.332 m2, obtained by the sum of all element failure risks. It is observed that the slope failure risk

**Figure 16.** EFP of the slope under all potential groundwater level acts.

can be observed that the groundwater level influences the EFR of the slope. When *T<sup>w</sup>* = 1, the frequency of EFR between 0 and 0.00062 m<sup>2</sup> is 675, the frequency of EFR between 0.00062 and 0.00124 m<sup>2</sup> is 36, the frequency of EFR between 0.00124 and 0.00186 m<sup>2</sup> is 58, the frequency of EFR between 0.00186 and 0.00248 m<sup>2</sup> is 69, the frequency of EFR between 0.00248 and 0.00310 m<sup>2</sup> is 140, the frequency of EFR between 0.00310 and 0.00372 m<sup>2</sup> is 11, and the slope failure risk is 0.829 m<sup>2</sup> . When *T<sup>w</sup>* = 25, the frequency of EFR between 0 and 0.00082 m<sup>2</sup> is 670, the frequency of EFR between 0.00082 and 0.00164 m<sup>2</sup> is 32, the frequency of EFR between 0.00164 and 0.00246 m<sup>2</sup> is 54, the frequency of EFR between 0.00246 and 0.00328 m<sup>2</sup> is 56, the frequency of EFR between 0.00328 and 0.00410 m<sup>2</sup> is 165, the frequency of EFR between 0.00410 and 0.00492 m<sup>2</sup> is 12, and the slope failure risk is 1.302 m<sup>2</sup> . When *T<sup>w</sup>* = 50, the frequency of EFR between 0 and 0.00150 m<sup>2</sup> is 666, the frequency of EFR between 0.00150 and 0.00300 m<sup>2</sup> is 41, the frequency of EFR between 0.00300 and 0.00450 m<sup>2</sup> is 48, the frequency of EFR between 0.00450 and 0.00600 m<sup>2</sup> is 62, the frequency of EFR between 0.00600 and 0.00750 m<sup>2</sup> is 158, and the frequency of EFR between 0.00750 and 0.00900 m<sup>2</sup> is 14. Under all potential groundwater level acts, the frequency of EFR between 0 and 0.04820 m<sup>2</sup> is 674, the frequency of EFR between 0.04820 and 0.09640 m<sup>2</sup> is 30, and the frequency of EFR between 0.09640 and 0.14460 m<sup>2</sup> is 57. The frequency of EFR between 0.014460 and 0.19280 m<sup>2</sup> is 59, the frequency of EFR between 0.19280 and 0.24100 m<sup>2</sup> is 157, the frequency of EFR between 0.24100 and 0.28892 m<sup>2</sup> is 12, and the slope failure risk is 2.094 m<sup>2</sup> . The element failure risk comprehensively reflects the contribution of the EFP and the element failure risk coefficient, which can make a quantitative judgment on the slope failure risk of each part. The slope failure risk is 1.332 m<sup>2</sup> , obtained by the sum of all element failure risks. It is observed that the slope failure risk assessment method proposed in this paper can avoid the screening and statistical work of failure modes compared with the failure risk calculation results in Table 3. *Sustainability* **2023**, *15*, x FOR PEER REVIEW 21 of 26 **Figure 15.** EFP of the slope under *Tw*th groundwater level acts: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50.

**Figure 16.** EFP of the slope under all potential groundwater level acts. **Figure 16.** EFP of the slope under all potential groundwater level acts.

The EFR of the slope under *Tw* th groundwater level acts when =1 *Tw* , =25 *Tw* , and =50 *Tw* obtained with Equation (20) is shown in Figure 17. The EFR of the slope under all potential groundwater level acts obtained with Equation (21) is shown in Figure 18. It can The safety performance decreases, and the IFP of the slope increases from 1.40% to 3.30% with the gradual increase in the groundwater level. Through cubic polynomial fitting, the IFP of the slope versus the groundwater level acquired (as shown in Figure 19) is as follows:

$$P\_f^{IFP}(T\_w) = 0.5271H\_w^3 - 11.2957H\_w^2 + 81.5794H\_w - 196.7446\tag{30}$$

frequency of EFR between 0 and 0.00062 m2 is 675, the frequency of EFR between 0.00062 and 0.00124 m2 is 36, the frequency of EFR between 0.00124 and 0.00186 m2 is 58, the frequency of EFR between 0.00186 and 0.00248 m2 is 69, the frequency of EFR between 0.00248 and 0.00310 m2 is 140, the frequency of EFR between 0.00310 and 0.00372 m2 is 11, The EFR of the slope increases from 0.829 m<sup>2</sup> to 2.094 m<sup>2</sup> with the gradual increase in the groundwater level. Through cubic polynomial fitting, the EFR of the slope versus the groundwater level acquired (as shown in Figure 20) is as follows:

$$R^{\rm EFP}(T\_w) = 0.1251 H\_w^3 - 2.4622 H\_w^2 + 16.6665 H\_w - 37.9602\tag{31}$$

quency of EFR between 0.00150 and 0.00300 m2 is 41, the frequency of EFR between 0.00300 and 0.00450 m2 is 48, the frequency of EFR between 0.00450 and 0.00600 m2 is 62, the frequency of EFR between 0.00600 and 0.00750 m2 is 158, and the frequency of EFR between 0.00750 and 0.00900 m2 is 14. Under all potential groundwater level acts, the frequency of EFR between 0 and 0.04820 m2 is 674, the frequency of EFR between 0.04820 and 0.09640 m2 is 30, and the frequency of EFR between 0.09640 and 0.14460 m2 is 57. The frequency of EFR between 0.014460 and 0.19280 m2 is 59, the frequency of EFR between 0.19280 and 0.24100 m2 is 157, the frequency of EFR between 0.24100 and 0.28892 m2 is 12, and the slope failure risk is 2.094 m2. The element failure risk comprehensively reflects the contribution of the EFP and the element failure risk coefficient, which can make a quantitative judgment on the slope failure risk of each part. The slope failure risk is 1.332 m2, obtained by the sum of all element failure risks. It is observed that the slope failure risk

0.00082 m2 is 670, the frequency of EFR between 0.00082 and 0.00164 m2 is 32, the frequency of EFR between 0.00164 and 0.00246 m2 is 54, the frequency of EFR between 0.00246 and 0.00328 m2 is 56, the frequency of EFR between 0.00328 and 0.00410 m2 is 165, the frequency of EFR between 0.00410 and 0.00492 m2 is 12, and the slope failure risk is

assessment method proposed in this paper can avoid the screening and statistical work of

failure modes compared with the failure risk calculation results in Table 3.

**Figure 17.** EFR of the slope under *Tw*th groundwater level acts: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50. **Figure 17.** EFR of the slope under *Tw*th groundwater level acts: (**a**) *Tw* = 1; (**b**) *Tw* = 25; (**c**) *Tw* = 50.

**Figure 18.** EFR of the slope under all potential groundwater level acts. **Figure 18.** EFR of the slope under all potential groundwater level acts.

*Sustainability* **2023**, *15*, x FOR PEER REVIEW 23 of 26

**Figure 19.** IFP of the slope versus the groundwater level (Unit: %). **Figure 19.** IFP of the slope versus the groundwater level (Unit: %).

groundwater level and soil shear strength parameters is proposed in this paper. The corresponding calculation program is compiled. The major conclusions are the following: (1) When the randomness of the groundwater level and soil shear strength parameters are considered comprehensively, the traditional LEM will ignore multiple failure modes and may miscalculate the slope failure risk. However, all failure modes can be acquired with the UBM for seeking the minimum value of the KAVF. Thus, the result is more consistent with the real situation. In addition, the traditional LEM only judges the slope stability by the safety factor, which only reflects the degree of the IPF. The EFP is used to calculate the EFR of the slope, which cannot only reflect the degree of the IFP but, also, the slope failure risk can be accurately acquired. It should

be noted that this calculation method can greatly reduce the calculation cost. (2) The IFP and EFR of the slope are increasing from 1.40% to 3.30% and 0.829 m2 to 2.094 m2 with the rise of the groundwater level, respectively. Based on the EFP, the proposed method can accurately obtain the EFR of the slope under each groundwater level act by using the element's location information and failure situation. This will provide engineers with realistic reference values for the slope reinforcement design

(3) Groundwater level and earthquakes are two important causes of slope instability and failure. However, this study does not consider the impact of earthquakes on slope reliability. Therefore, relevant studies on seismic slope stability will be carried out in the future. In addition, according to the upper bound theory, the upper bound solution is inevitably greater than the true solution. Therefore, the failure probability will be underestimated when using the UBM for slope reliability analysis. To solve this problem, there is a necessity to study the slope reliability calculation method on the basis of the lower bound theory in future research work. The solution of slope failure probability with the UBM and LBM can be obtained at the same time, so the interval range of the real failure probability can be accurately judged, and the reliability index

**Author Contributions:** Conceptualization, P.P., Z.L., X.Z., W.L., S.S., and H.X.; methodology, P.P.; software, Z.L.; validation, W.L.; formal analysis, Z.L.; investigation, Z.L.; resources, W.L.; data curation, W.L.; writing—original draft preparation, Z.L.; writing—review and editing, P.P.;

**Figure 20.** EFR of the slope versus the groundwater level (Unit: m<sup>2</sup> ).

**Figure 20.** EFR of the slope versus the groundwater level (Unit: m2).

to achieve sustainable development.

of the slope can be quantified more accurately.

**6. Conclusions** 

#### **6. Conclusions**

A slope failure risk assessment method that considers both the randomness of the groundwater level and soil shear strength parameters is proposed in this paper. The corresponding calculation program is compiled. The major conclusions are the following:


**Author Contributions:** Conceptualization, P.P., Z.L., X.Z., W.L., S.S. and H.X.; methodology, P.P.; software, Z.L.; validation, W.L.; formal analysis, Z.L.; investigation, Z.L.; resources, W.L.; data curation, W.L.; writing—original draft preparation, Z.L.; writing—review and editing, P.P.; visualization, H.X.; supervision, W.L. and Z.L.; project administration, S.S. and Z.L.; funding acquisition, X.Z. and Z.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China grant number 12162018 and 12262016; the Research Foundation of Science and technology plan project of Yunnan Provincial Department of Science and Technology grant number 202201AT070283 and 202305AD160064.

**Institutional Review Board Statement:** Not applicable for studies not involving humans or animals.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
