Probabilistic Risk Assessment of Soil Slope Stability Subjected to Water Drawdown by Finite Element Limit Analysis

Abstract

This study investigates the probabilistic stability of embankment slopes subjected to water level drawdown using the random field finite element method (RFEM) with strength reduction technology. The shear strength of soil properties was controlled by cohesion and internal friction angle for the slope shear failure. The cohesion and internal friction angle were modeled by a random field following the log-normal distribution. The factor of safety (FOS) for the embankment slope with random soil is calculated by strength reduction technology. During the numerical simulation, the limit analysis upper bound and lower bound method are applied to the finite element method, respectively, to obtain the upper bound and lower bound value of the FOS. Seepage action is also considered during the water drawdown by setting five different water levels (WLs). A total of 1000 Monte Carlo simulations are performed for each work condition, resulting in histograms of the FOSs. The results show that the FOSs obtained by the random field model are all lower than those by the deterministic method. Even if the FOSs obtained by the two methods are close, there still exists the possibility of slope failure. Compared to the deterministic results, the RFEM method is more reasonable for evaluating slope stability.

1. Introduction

In geotechnical engineering, embankment slope failure has always been an important issue within one of the most common failure modes. As one of the water-retaining structures, embankments play an important role in protecting the lives and properties further downstream during flood periods. The material of an embankment slope is always simplified as homogeneous, and the FOS is used to describe the failure mechanism in the conventional method. The limit equilibrium theory is used to evaluate slope stability in geotechnical engineering [1,2]. The shear strength parameters, e.g., cohesion and internal friction angle, are generally considered as deterministic parameters in the conventional method. Mahmood et al. [3] investigated the effects of rainfall intensity on the unsaturated slope. Ijaz et al. [4] studied the surficial stability of unsaturated expansive soil slope considering the hydro-mechanical effect. However, geotechnical material has obvious uncertainty and randomness as well as spatial correlation structure after long-term consolidation [5]. Phoon and Kulhawy [6] pointed out that there are three primary sources of geotechnical uncertainties: inherent soil, measurement errors and transformation uncertainties, and the vertical correlation length is less than the horizontal correlation length (l_x). Many previous studies have shown that the strength parameters, e.g., cohesion and internal friction angle, of geomaterials follow Gaussian distributions or log-normal distributions [7,8]. So far, probabilistic analysis has been widely adopted in geotechnical engineering to describe the uncertainty of geomaterials [9,10,11,12,13,14]. The reliability theory was used to investigate the structure failure in which the uncertainty of materials can be considered. The point estimation method was used to calculate the reliability index with the advantage of time saving and simplicity [15,16,17]. The first-order second-moment method was also used to study the reliability of geotechnical engineering [18,19]. Cho [20] studied the reliability index of two-layered slope stability considering the spatial variability of soil. The limit state function was used to calculate the failure probability and reliability index of engineering structures; however, the limit state function of the complex structures sometimes cannot be obtained easily. One of the obvious shortcomings of reliability theory is that the spatial correlation structure of the soil cannot be considered.

In recent decades, more advanced probabilistic analysis methods have been developed to study the failure of soil structures caused by the uncertainty of geomaterials. Ahmed and Soubra [21] analyzed the two neighboring strip footings by considering the spatial randomness of the soil using subset simulation coupled with the random field model. Griffiths and Fenton [22] first studied steady-state seepage by combining random field theory and the finite element method (RFEM), and since then, a lot of work using RFEM has been conducted to study geotechnical engineering, including slope stability [23,24,25,26]. This excellent work has attracted increasing attention and has been used to study various uncertainties in geotechnics. Based on different study cases, seepage analyses were also conducted using RFEM to describe the spatial variation [27,28,29,30]. Wang et al. [31] studied the multi-stage slopes using the bivariate RFEM and concluded that although the FOS obtained from deterministic results was close to the mean value based on random field results, there was still a high risk of instability. Li D et al. [32] investigated reservoir slope failure considering the spatial variability of the soil properties and concluded that the spatial variability has a strong effect on the stability.

In this study, the influence of spatial variability of the soil strength parameters, cohesion and internal friction angle, on the slope stability of embankments is studied by RFEM coupling with seepage. The influences of l_x and the coefficient of variation (COV) of shear strength on the probabilistic slope stability are highlighted. The FOS of the embankment subjected to slow drawdown is also studied. The numerical simulation method adopts finite element limit analysis.

2. Analysis Method of Seepage and Stability

2.1. Darcy Seepage

The steady-state flow through porous media can be described by the general Darcy law in two-dimensional space.

$$q_x=k_x\frac{\partial h}{\partial x}\mathrm{d}y$$

(1a)

$$q_y=k_y\frac{\partial h}{\partial y}\mathrm{d}x$$

(1b)

where $q_x$, $q_y$ is the fluid velocity along the x and y directions, respectively, ${h=h}^{*}-p_s/{\mathit{\gamma }}_w$ is the total head, $h^{*}$ is the position head, $p_s$ is the static water pressure, ${\mathit{\gamma }}_w$ is the water gravity, and $k_x$ and $k_y$ are saturated hydraulic conductivity along the horizontal direction and vertical direction. Considering the mass conservation equation and assuming $k_x=k_y$ lead to the following classic steady-state seepage Laplace equation.

$$\frac{{\partial }^{2}h}{\partial x^2}+\frac{{\partial }^{2}h}{\partial y^2}=0$$

(2)

2.2. Equilibrium Equation and Strain-Displacement Relationships

The static equilibrium equations for an objective and the boundary conditions and strain–displacement relationship under the small deformation assumption are given by

$$\{\begin{array}{l}{\nabla }^{\mathrm{T}}\sigma +b=0 in V\\ P^{\mathrm{T}}\sigma =\alpha t on S_{\mathsf{\sigma }}\\ u=u_{\mathrm{S}} on{S}_u\\ \epsilon =\nabla u\end{array}$$

(3)

where ${\nabla }^{\mathrm{T}}$ is the usual equilibrium operator, $\sigma$ is the total stress tensor, $b$ is the body force vector, which is usually the gravity vector, $V$ is the configuration of interest, $P$ is the projection matrix and $\mathsf{\alpha }$ is a factor of load. $t$ are the tractions applied on the surface, $S$ is the boundary, $\epsilon$ is the strain tensor, $u$ is the displacement vector, $u_{\mathrm{S}}$ is specified boundary displacement vector and $\nabla$ is the strain–displacement operator which is dual to the equilibrium operator.

In recent years, finite element limit analysis has become a popular method in geotechnical engineering [33,34,35,36,37,38], including lower bound limit analysis and upper bound limit analysis, which are different from the limit equilibrium method. The core question of limit analysis is to determine the maximum magnitude of the external load that can be sustained without the structure suffering collapse [36,39].

The scaled loads are applied to the slope during the analysis, and the above governing equations can be rewritten as an optimization problem of variation principle given by [37,39]:

$$\begin{array}{c}\underset{u}{\min }\underset{\sigma ,\mathsf{\alpha }}{\max }\mathsf{\alpha }+{\displaystyle \underset{V}{\int }{\sigma }^{\mathrm{T}}\nabla \dot{u}}\mathrm{d}V-{\displaystyle \underset{V}{\int }{b}^{\mathrm{T}}\dot{u}}\mathrm{d}V-{\displaystyle \underset{S}{\int }\mathsf{\alpha }{t}^{\mathrm{T}}\dot{u}\mathrm{d}S}\\ s. t. F\left(\sigma \right)\le 0\end{array}$$

(4)

2.3. Constitutive Relationship in Effective Stress Space

2.3.1. Effective Stress Principle

One of the most significant features for a porous soil-water structure is that the effective stresses cannot be ignored. Terzaghi’s effective stress principle is adopted in this paper to describe the stress state of porous soil material:

$${\sigma }^{\prime }=\sigma -\omega p_{\mathrm{s}}$$

(5)

where ${\sigma }^{\prime }$ are effective stresses, $\sigma$ are total stresses, $\omega ={\left(1,1,0\right)}^{\mathrm{T}}$ responding to the plane–strain problem, and $p_s$ is the pore pressure due to the seepage. That is, only the static water pressure is considered without considering the excess pore pressure.

2.3.2. Rigid-Plastic Constitutive Relationship

The material behavior is assumed to be rigid-plastic, and the deformation at incipient damage is of a magnitude based on the small strain theory. This means that there is no deformation before the yield occurs. The constitutive relation using effective stress is given as follows:

$$\{\begin{array}{l}F\left({\sigma }^{\prime }\right)\le 0\\ {\dot{\epsilon }}^{\mathrm{p}}= \dot{\mathsf{\lambda }}\partial G/\partial \sigma \\ \dot{\mathsf{\lambda }}F\left({\sigma }^{\prime }\right)=0 \dot{\mathsf{\lambda }}\ge 0\end{array}$$

(6)

where $F$ represents the yield functions, ${\dot{\epsilon }}^{\mathrm{p}}$ is the plasticity strain rate tensor, $G$ is the plastic potential function and $\dot{\lambda }$ is the plastic multiplier.

2.4. Finite Element Limit Analysis with Strength Reduction Technology

2.4.1. Lower Bound Limit Analysis

Figure 1 shows a typical triangular element adopted in the lower bound limit analysis [33,38]. The stresses obey equilibrium and satisfy both the boundary condition and the yield function. The objective function corresponding to the governing equations considering the effective stresses in terms of the lower bound limit analysis is given as below:

$$\begin{array}{lll}\mathrm{Maximize}& & \mathsf{\alpha }\\ \mathrm{Subject}\mathrm{to}& & \{\begin{array}{l}{\nabla }^{\mathrm{T}} \left({\sigma }^{\prime }+\omega p_s\right)+b=0 in V\\ P^T \left({\sigma }^{\prime }+\omega p_s\right)=\mathsf{\alpha }t+t_s on S\\ F^T {\sigma }^{\prime }-k+\kappa =0 s\ge 0\end{array}\end{array}$$

(7)

where $t_s$ represents the seepage pressures. The yield function $F\left({\sigma }^{\prime }\right)\le 0$ is replaced by a set of linear constrains $f_i^T{\sigma }^{\prime }-k_i\le 0, \left(i=1\sim n\right)$ and then can be further rewritten as $F^T {\sigma }^{\prime }-k+\kappa =0,\kappa \ge 0$, where $F$ and $k$ collects the contributions, $f_i$ and $k_i$, respectively, and $\kappa$ represents slack variables.

2.4.2. Upper Bound Limit Analysis

Figure 2 shows a typical triangular element adopted in the upper bound limit analysis [34,38]. Every node has two velocity components, and the element is associated with a vector of three unknown stresses. Velocity fields follow an associated flow rule and satisfy the boundary conditions. The objective function corresponding to the governing equations considering the effective stress in terms of the upper bound limit analysis is given as below:

$$\begin{array}{lll}\mathrm{Minimize}& & {\displaystyle \underset{V}{\int }{k}^{\mathrm{T}}\dot{\mathsf{\lambda }}\mathrm{d}V}-{\displaystyle \underset{V}{\int }{\left[{\nabla }^{\mathrm{T}} \left(\omega p_s\right)+b\right]}^{\mathrm{T}}\dot{u}\mathrm{d}V}+{\displaystyle \underset{S}{\int }{t}_s^T\dot{u}\mathrm{d}S}\\ \mathrm{Subject\; to}& & \{\begin{array}{l}{\nabla }^{\mathrm{T}}u=F\dot{\mathsf{\lambda }} \dot{\mathsf{\lambda }}\ge 0\\ {\displaystyle \underset{V}{\int }{t}^T\dot{u}\mathrm{d}S}=1\end{array}\end{array}$$

(8)

In this paper, the yield criterion and plastic potential function G are given in terms of the Mohr–Coulomb relations based on the plane strain assumption [40,41],

$$\sqrt{{\left({\mathit{\sigma }}_{xx}-{\mathit{\sigma }}_{yy}\right)}^2+4{\mathit{\tau }}_{xy}^2}+\left({\mathit{\sigma }}_{xx}+{\mathit{\sigma }}_{yy}\right)\sin \mathit{\phi }\le 2c\cos \mathit{\phi }$$

(9)

$$G=\sqrt{{\left({\mathit{\sigma }}_{xx}-{\mathit{\sigma }}_{yy}\right)}^2+4{\mathit{\tau }}_{xy}^2}+\left({\mathit{\sigma }}_{xx}+{\mathit{\sigma }}_{yy}\right)\sin \mathit{\psi }$$

(10)

where $c$ is the cohesion, $\mathit{\phi }$ is the internal friction angle, and $\mathit{\psi }$ is the dilation angle with $\mathit{\psi }=0$ implying incompressibility.

The optimization calculation process is solved employing the second-order cone programming method that has been embedded in the Optum G2 program [42], which is commercially available and facilitates the straightforward use of this finite element limit analysis and strength reduction technology.

2.4.3. Strength Reduction Technology

The two most common methods for slope stability analysis are the limit equilibrium method and the finite element strength reduction technology. Due to the clear concept of the limit equilibrium method, it has been widely used by engineers and researchers to assess slope stability. However, there are some shortcomings when using the limit equilibrium method, which gives different FOSs based on different assumptions. More comparisons between the limit equilibrium method and finite element strength reduction method can be found in Liu et al. [43].

The FOS of this embankment slope is calculated by the strength reduction method. Comparing with the limit equilibrium method, the strength reduction technology has the advantage that it does not need to assume the position of the shear failure surface in advance. The shear failure surface can be searched automatically during the process of strength reduction. In this method, to obtain the minimum strength reduction factor $f_s$ of the embankment slope, a series of $f_s$ values were set to determine the actual FOS until the slope failed. For simplicity, the Mohr–Coulomb yield model is used to describe the failure behavior. The strength reduction process in the Mohr–Coulomb model is shown in Figure 3 and follows the Formulas (11) and (12).

$$c^{\prime }=c/f_s$$

(11)

$$\mathit{\phi } \prime =\arctan (\tan \mathit{\phi }/f_s)$$

(12)

where $c\prime$ and $\mathit{\phi }\prime$ are the strength parameters after reduction and $f_s$ is the strength reduction factor. The slope is stable if $f_s>1$, in a critical state if $f_s=1$ and unstable if $f_s<1$.

3. Random Field Generation

3.1. Spatial Correlation Model

The most common spatial correlation models of one-dimensional space are the exponential model, Gauss model, spherical model and second-order autoregressive model. The exponential model (also called Markovian decaying function) is used in this paper to describe the spatial correlation structure of the randomness of the soil. The formulas in one-dimensional and two-dimensional expressions are as follows:

$$\mathit{\rho }\left(h\right)=\exp \left(-2\frac{d}{l}\right)$$

(13)

$$\rho (d_x,d_y)=\exp \left[-\sqrt{{\left(\frac{2d_x}{l_x}\right)}^2+{\left(\frac{2d_y}{l_y}\right)}^2}\right]$$

(14)

where $\rho$ is the correlation coefficient of different positions, $d$ is the separated distance in the random field, $l$ is the correlation length, $d_x$, $d_y$ are the separated distances along the horizontal direction and vertical direction, respectively, and $l_x$, $l_y$ are the horizontal correlation length and vertical correlation length in two-dimensional space.

A plot of this exponential model in two-dimensional space is given in Figure 4, indicating the properties at two points separated by l_x and l_y.

Figure 4 shows that the soil properties will be almost the same if the distance between two points is very close. A large correlation length implies that the soil property is correlated over a large spatial extent, resulting in a smooth variation within the soil profile. On the other hand, a small correlation length indicates that the fluctuation of the soil property is serious.

3.2. Generation of Random Fields Using the Karhunen–Loeve Series Expansion Method

The spatial variability of random soils is represented by a random field model. Several common models have been developed and used to generate random fields such as the covariance matrix decomposition method [44,45,46], the local average subdivision method [9,47,48] and the Karhunen–Loeve (K-L) series expansion method [27,49,50,51]. This paper applies the following K-L series expansion method to generate the Gauss random field because of its efficiency and precision.

$${\widehat{H}}_i\left(x\right)={\mu }_{X_i}+{\mathrm{s}}_{X_i}{\widehat{H}}_i{}^D\left(x\right)={\mu }_{X_i}+{\mathrm{s}}_{X_i}{\displaystyle \sum _{j=1}^M\sqrt{{\mathsf{\lambda }}_j}{f}_j\left(x\right){\chi }_{ji}}$$

(15)

where ${\widehat{H}}_i\left(x\right)$ is the desired random field considering the mean value ${\mu }_{X_i}$ and standard deviation ${\mathrm{s}}_{X_i}$, i is the i-th random field, $x$ is the spatial position, ${\widehat{H}}_i{}^D\left(x\right)$ is a standard Gaussian random field after discretizing the space, $M$ is the number of series expansion terms of a random field, ${\lambda }_j$ and $f_j\left(x\right)$ are the $j^{\mathrm{th}}$ eigenvalue and $j^{\mathrm{th}}$ eigenfunction of spatial correlation model $\rho$, respectively, $\chi$ is a stochastic matrix of size $M\times N$ considering the autocorrelation of the random field, $\chi =\xi L^{\mathrm{T}}$, where $\xi$ is a stochastic vector of size $M\times N$ in standard Gaussian space, and $L$ is the lower triangular matrix obtained by the Cholesky decomposition of matrix $\rho$ in standard Gaussian space, and its size is $M\times N$.

3.3. Probability Distribution of Strength Parameters

The cohesion and internal friction angle are considered as random variables in this numerical model of slope stability for this test embankment.

Based on previous studies, the shear strength parameters, e.g., cohesion and internal friction angle in geotechnical engineering can be modeled as Gaussian distributions or log-normal distributions [5,31]. Since the shear strength parameters in geotechnical engineering are always positive, in this study, $c$ and $\mathit{\phi }$ are assumed to be characterized statistically by a log-normal distribution defined by two mean values, ${\mathit{\mu }}_c$ and ${\mathit{\mu }}_{\mathit{\phi }}$, and two standard deviations $s_c$ and $s_{\mathit{\phi }}$, respectively. Once the mean and standard deviation are expressed in terms of the dimensionless COV, defined as $\mathrm{COV}=s/\mathit{\mu }$, the mean value and standard deviation of the underlying Gauss distribution of $\ln c$ and $\ln \mathit{\phi }$ are given by [28,31,49],

$$s_{\ln c}=\sqrt{\ln \left(1+{\mathrm{COV}}_c^2\right)}$$

(16a)

$${\mathit{\mu }}_{\ln c}=\ln {\mathit{\mu }}_c-0.5s_c^2$$

(16b)

and

$$s_{\ln \mathit{\phi }}=\sqrt{\ln \left(1+{\mathrm{COV}}_{\mathit{\phi }}^2\right)}$$

(17a)

$${\mathit{\mu }}_{\ln \mathit{\phi }}=\ln {\mathit{\mu }}_{\mathit{\phi }}-0.5s_{\mathit{\phi }}^2$$

(17b)

4. Results and Discussion

4.1. Monte Carlo Simulation

There are many stochastic methods for probabilistic slope analysis in previous studies, but the Monte Carlo simulation is still the most common technique. Even though a large amount of computing time is required to use this technology, the mathematical formulation of the Monte Carlo method is relatively simple.

During the process of probabilistic analysis, the finite element method with strength reduction technology is repeated enough times to reach a stable value. Figure 5 shows the convergence of the estimated mean value and standard deviation of FOS from 2000 Monte Carlo simulations for the test embankment slope. The result converges to a stable value when the number of Monte Carlo simulations reaches about 1000 from Figure 5. The 1000 random field will be generated during the 1000 Monte Carlo simulations; then, 1000 FOSs could be obtained based on the mean value, standard deviation and the spatial correlated structural of the strength parameters. The average of the 1000 FOSs will be used as the result to evaluate the stability. Indeed, there is still an error compared with the true solution; however, it is acceptable for geotechnical engineering. Therefore, all the results of FOS are obtained by performing Monte Carlo simulations 1000 times, for the balance of computing time and numerical error. These results are used to describe the probabilistic stability of the slope.

4.2. The Characters of a Test Embankment Slope

The instability of an embankment slope is a very complex physical process, which involves seepage, boundary conditions, pore water pressure, stress history of the soil and so on. Herein, the FOS of the embankment slope subjected to slow drawdown (implying that the excess pore water pressure is neglected) is studied considering the seepage and spatial variability of the soil.

The numerical modeling is conducted on a test embankment slope as shown in Figure 6. This embankment constructed on a 5.5 m thick soil foundation that is 5.5 m high, with upstream and downstream slopes of 2h:1v, and the water level (WL) in the river is 5 m above the foundation. To evaluate the stability of this embankment slope, the FOS is calculated by the finite element limit analysis with the strength reduction method mentioned in Section 2.4.

The soil properties are displayed in

Table 1. Statistical properties of soil parameters.

Parameter	Mean	Coefficient of Variation	Horizontal Correlation Length
c (kPa)	15	0.1, 0.2, 0.3, 0.4, 0.5	1 m, 30 m, 40 m, 50 m, 500 m
$\phi (°)$	20	0.1, 0.2, 0.3, 0.4, 0.5	1 m, 30 m, 40 m, 50 m, 500 m
E (MPa)	40	-	-
$\mathit{\nu }$	0.3	-	-
${\mathit{\gamma }}_{\mathrm{d}}$$(\mathrm{kN}/{\mathrm{m}}^3$)	18	-	-
${\mathit{\gamma }}_{\mathrm{sat}}$$(\mathrm{kN}/{\mathrm{m}}^3$)	20	-	-
ks (m/s)	2 × 10−6	-	-

. The Young’s modulus $E$, Poisson ratio $\mathsf{\nu }$ and hydraulic conductivity $k_s$ are assumed to be deterministic, as they have turned out to have little influence on the results of the slope stability analysis. Several previous studies [6,7] pointed out that l_y is much smaller than l_x in geotechnics, so that l_y is fixed to 1 meter. The l_x of the cohesion and internal friction angle is set at five different levels. In particular, the dilation angle is set to zero, corresponding to no volume change during the process of yielding.

4.3. Deterministic Analysis Results

A series of deterministic analyses are carried out in OptumG2 before the probability analysis during the process of the slow drawdown of the WL. The finite element limit analysis with the mesh adaptivity technique is used for this numerical model. Figure 7 shows that the number of the adaptive finite elements in 2000 is suitable for obtaining stable FOSs.

Several different WLs are set to simulate the slow water drawdown, because the pore water pressure can completely dissipate during the process of water falling slowly. The FOSs of the embankment slope at six different WLs of 5 m, 4 m, 3 m, 2 m, 1 m, and 0.5 m are calculated to describe the slow water drawdown. As shown in Figure 8, the FOSs increase with the slow water drawdown. The reason is that it is assumed to be a steady-state seepage during the slow drawdown. The effective stress increases with the dissipation of pore water pressure, and the matric suction also increases during the slow drawdown. These two aspects will strengthen the stability of the embankment slope.

As mentioned in Section 2.4.3 of this paper, the upper and lower bound values of FOSs are calculated by strength reduction technology. The true solution then falls somewhere between the upper bound and lower bound values, and thus, the mean value of the upper and lower values is regarded as the true solution, as shown in Figure 8. The probabilistic results described in the following sections are also calculated using this method.

4.4. Analysis of Random Field Results

In this section, there are in total 150 random field analysis results obtained from five different l_x and five different COVs under six different WLs. The detailed computational setup can be found in the Appendix A, and corresponding FOSs are displayed in Figure 9.

Furthermore, to compare with the deterministic results in Figure 8, during the slow water drawdown, one group of random field results with the parameter set, COV = 0.1, l_x = 50 m, corresponding to six water levels are plotted as shown in Figure 10. It shows the variation of the statistical FOSs during the slow drawdown of the WL, during which the FOSs present a gradual upward trend with the decrease in the WL. Such a result comes from two reasons. The first one is that the effective stress of soil increases with the decrease in WL because the pore water pressure disappears. The other one is that the increase in mechanical bite force leads to the soil slope being more stable.

With the same WL, all FOSs obtained from the random field analysis are smaller than the results obtained from the deterministic analysis. In other words, the deterministic analysis method may sometimes overestimate the stability of the embankment slope, having potential risk. It implies that it is more reasonable to adopt the random field analysis, especially when the objective is the natural geomaterials. It is interesting to note that all the other results from the different parameter sets have a similar trend to Figure 10.

4.4.1. Effect of Different COV

Some literature studied that the shear strength parameters, cohesion and internal friction angle, in geotechnical materials varied from 0.1 to 0.5 [6]. The FOS is calculated, and the parametric study of the spatial variation of shear strength parameters is investigated during the slow water drawdown. All the results of the random field have similar properties with the histograms of the FOSs depicted in Figure 11, which correspond to the parameter set (l_x = 50 m, WL = 5 m) versus five different COVs. For comparison, the same deterministic results were 1.684, as marked with an arrow in Figure 11. Figure 11 also shows that all the mean values in the five subfigures are lower than the deterministic value of 1.684. The differences between the mean value and deterministic value become larger with the increase in COV of shear strength. The statistical result from Figure 11 shows that although the FOS obtained from deterministic FOS is close to the mean value of that based on random field analysis, still nearly half of the FOSs are less than the deterministic value. The result from Figure 11e shows that there is a large difference between the mean value and the deterministic value, and the failure probability was about 23.3%. This result presents a high risk of landslides for the overall slope system, despite the mean value from the random field analysis and the deterministic value both being located at a relatively safe level. Figure 11 also shows that the range of FOSs increases with the increase in the COV, having similar properties to the previous study. The result also indicates that the random field method is a more reasonable method for assessing the stability of the geotechnical slope.

Figure 12 shows the contour of the shear dissipation energy at the last stage of the failure slope for this test embankment, in which the subfigures correspond to the subfigures in Figure 11. The shear strength of the soil set in this test embankment is strong because the FOS is relatively high, resulting in the failure of circular slip surface through the toe shown in Figure 12. The failure circular slip surface moved toward the bottom of the foundation with the COV of the shear strength becoming larger. In other words, the magnitude of the failure soil increased with the COV of the shear strength becoming larger.

Figure 13 shows that the mean FOS decreases as the COV increases. The reason is that when the COV becomes greater, the fluctuation of the uncertainty of the soil strength becomes more serious. It is more possible to take much weaker strength parameters, cohesion and internal friction angle, so that the FOS becomes lower. This is consistent with the real situation that there are strong and weak areas in an actual slope system. The failure of the embankment slope occurs easily in weak soil areas.

4.4.2. Effect of Different l_x

Because of the long-term depositional process and load history, the spatial correlation structure of geomaterials should be considered as a basic property that cannot be ignored. In previous studies, the magnitude of the vertical correlation length l_y which is from 1 to 2 m, is always lower than the horizontal correlation length l_x, which is usually from 30 to 60 m [6]. In this paper, the l_y is assumed to be 1 m as a fixed value, and the influence on the FOS of five different l_x is mainly studied. Figure 14 shows five different random fields of cohesion generated by the K-L series expansion method. The obvious spatial correlation structure can be observed by the contours.

Figure 15 shows the variation of mean FOSs with l_x for different COVs at the 5 m high water level. The mean FOSs slightly decrease with the increase in the l_x for the range 1–50 m, and this is in good agreement with Zhu et al. [28]. The mean FOSs reach a minimum when the l_x ranges from 30 to 50 m, because the width of the foundation of this test embankment in this paper is 36 m. The stability of this test model will almost not be influenced when the l_x exceeds 50 m. Hence, the end part of the curve in Figure 11 tends to be horizontal. This is also consistent with real engineering, and the spatial correlation length is constant for a specific case, and the maximum length is no more than the outline. However, the standard deviation of FOS increases with the increase in l_x. The increasing trend is obvious when the COV is high.

5. Conclusions

In this study, the finite element limit analysis was used to study the stability of embankment slopes subject to the slow water drawdown, considering the spatial variability of geomaterials. We assumed that the pore water pressure in the soil body has been dissipated during the slow drawdown. Hence, five different WLs were set to represent that the water level continues to drop.

The spatial variability of geomaterials is described by a random field model, which is determined by the mean, standard deviation and spatial correlation length of the shear strength cohesion and internal friction angle. Both the cohesion and internal friction angle are assumed to follow a log-normal distribution. The range of the COV is set from 0.1 to 0.5. The vertical correlation length is fixed at 1 m, and five different horizontal correlation lengths are set, namely 1 m, 30 m, 40 m, 50 m, and 500 m. All the mean values of FOS obtained from 1000 Monte Carlo simulations for different working conditions are lower than the corresponding deterministic analysis. The FOS increases with the decrease in the WL and decreases with the increase in COV. This is consistent with the actual situation that the shear failure surface is more likely to occur at the position where the soil strength is weaker. Although the mean value of FOS was greater than 1 for this test embankment, the probability of failure is still high, e.g., the probability of failure is 23.3% in Figure 11e.

The study shows that the probability of potential risk is high because the FOS obtained from the deterministic analysis is higher than the random field analysis result. Hence, it is more reasonable to evaluate the stability of an embankment slope using the probability framework rather than the deterministic analysis. Further work will focus on the influence of rapid water drawdown on slope stability as well as the random field action.