Probabilistic Slope Stability Analysis for Embankment Dams

Abstract

Slope instability is one of the most common forms of dam failure. The commonly used slope stability analysis methods ignore the uncertainty and randomness of dam materials, which may overestimate the stability of dams. In this chapter, a deterministic slope stability analysis based on strength reduction finite-element method is introduced first. After that, the slope is investigated using simple probabilistic concepts and classical slope stability techniques, and the shear strength is treated as a single random variable. Further, the random finite-element method (RFEM) is shown, in which spatial correlation and local averaging are illustrated in detail. Finally, the RFEM is applied to slope stability risk assessment, and the results can lead to higher probabilities of failure.

1. Introduction

Slope instability is one of the most common forms of dam failure. Traditional slope stability analysis methods mainly depend on deterministic analysis, including limit equilibrium analysis and finite-element (FE) analysis. Equilibrium methods mainly include the ordinary method of slices, Bishop’s modified method, force equilibrium methods, Janbu’s generalized procedure of slices, Morgenstern and Price’s method, and Spencer’s method. All the equilibrium methods assume that the soil can be divided into slices, which is an artificial distinction. This assumption is the main characteristic that distinguishes different limit equilibrium methods. The main advantage of equilibrium methods is that they involve relatively simpler calculation, which leads to wide use [1, 2, 3, 4].

While the finite element method is another powerful approach for slope stability analysis, it can better reflect the stress–strain relationship of soils than the equilibrium methods. Slope failure in the finite-element model occurs naturally through the area in which the shear strength of the soil is insufficient to resist the shear stresses. There are several advantages of a FE approach to slope stability analysis over traditional limit equilibrium methods: (a) there is no assumption about the shape or location of the failure surface, (b) there are no slices and slice side forces, and (c) the FE method is able to monitor progressive failure up to and including overall shear failure [5, 6].

For a practical slope, not only the stress–strain relationship of soils but also the uncertainty of soil properties should be taken into consideration. However, traditional slope stability analysis methods always ignore the uncertainty and randomness of dam materials, which may overestimate the stability of dams. Attention was drawn to probabilistic slope stability analyses [7, 8]. Most probabilistic slope stability analyses continue to use classical slope stability analysis techniques which are mainly based on the equilibrium methods [9, 10, 11, 12]. An obvious deficiency of the traditional slope stability methods is that the shape of the failure surface is always fixed; therefore, the failure mechanism is not allowed to look for the most critical path through the soil. Besides, these traditional methods cannot take the importance of spatial correlation and local averaging of statistical geotechnical properties into consideration [13, 14, 15].

A more rigorous method, in which nonlinear finite-element methods are combined with random field generation techniques, called the random finite-element method (RFEM), was proposed by Griffiths and Fenton [16]. It can fully account for spatial correlation and averaging and is also a powerful slope stability analysis tool that does not require a priori assumptions relating to the shape or location of the failure mechanism.

In this chapter, a deterministic slope stability analysis based on strength reduction finite-element method is introduced first. After that, the slope is investigated using simple probabilistic methods, including first-order second-moment (FOSM) method, first-order reliability method (FORM), and Monte Carlo method. Further, RFEM is shown, in which spatial correlation and local averaging are illustrated in detail. Finally, the RFEM is applied to slope stability risk assessment, and the results can lead to higher probabilities of failure.

2. Deterministic slope stability analysis

Deterministic slope stability analysis in this chapter is based on FE analysis. The program used is called SLOPE64 [6]. This program is for two-dimensional plane strain analysis. The soil is assumed to follow a linear elastic-perfectly plastic behavior characterized by the Mohr-Coulomb shear failure criterion. In the gravity load generation, the stiffness matrix generation, and the stress redistribution procedure, the program uses eight-node quadrilateral elements with simplified integration (four Gauss points per element). Initially, the soil is assumed to be elastic, and the model generates normal and tangential stresses at all Gauss points in the grid. These stresses are then compared with the Mohr-Coulomb failure criterion. If the stress at a particular Gauss point is within the Mohr-Coulomb failure envelope, it is assumed that the position remains elastic. If the stress is on or outside the failure envelope, it is considered that the point is yielding. The yield stresses are redistributed in the whole grids by the viscoplastic algorithm. Overall shear failure occurs when a sufficient number of Gauss points have yielded to allow a mechanism to develop [5, 6].

2.1 Soil model

The soil model used in this program consists of six parameters, as shown in Table 1.

ϕ′	Friction angle
c′	Cohesion
ψ	Dilation angle
E′	Young’s modulus
ν′	Poisson’s ratio
γ	Unit weight

Table 1.

Six parameters for soil model.

The Mohr-Coulomb failure criterion used in this program can be written as follows:

$$F=\frac{{\sigma }_1^{\prime }+{\sigma }_3^{\prime }}{2}\mathit{\text{sinϕ}}\prime -\frac{{\sigma }_1^{\prime }-{\sigma }_3^{\prime }}{2}-c\prime \mathit{\text{cosϕ}}\prime$$E1

where ${\sigma }_1^{\prime }$ and ${\sigma }_2^{\prime }$ are the major and minor principal effective stresses, respectively.

The failure function F can be described as follows:

F < 0 stresses inside the failure envelope (elastic).

F = 0 stresses on the failure envelope (yielding).

F > 0 stresses outside the failure envelope (yielding and must be redistributed).

2.2 Determination of the factor of safety (FS)

The FS of a soil slope is defined as the ratio between the strength of the soils and the actual load. It is exactly the same as that used in traditional limit equilibrium methods. The factored shear strength parameters c_f′ and ϕ_f′ are therefore given by

$$\begin{array}{l}{c}_f^{\prime }=c\text{'}/\mathit{FS}\\ {\varphi }_f^{\prime }=arctan\left(\frac{tan{\varphi }^{\prime }}{\mathit{FS}}\right)\end{array}$$E2

In this program, in order to find the actual FS, it is necessary to start a systematic search for FS values that will cause the slope to fail. This is achieved by the program that repeatedly solves problems using a sequence of user-specified FS values.

2.3 Slope stability analysis examples

Figure 1 shows a homogeneous slope with a foundation layer. The height of the slope (H) is 10 m, and the thickness of the foundation layer is H/2, 5 m. Soil parameters are shown in Table 2.

Figure 1.

Homogeneous slope with a foundation layer.

ϕ′	c′	ψ	E′	ν′	γ
20°	10 kPa	0	10,000 kN/m^2	0.3	20 kN/m^3

Table 2.

Soil parameters.

Figure 2 shows the undeformed mesh of the homogeneous slope. The slope is inclined at an angle of 26.578° (2:1). The left boundary is fixed horizontally but is free along the vertical direction, and the base boundary is fixed in both directions. Gravity loads were applied to the mesh, and the trial FS gradually increased until convergence could not be achieved within the iteration limit. The deformed mesh and the nodal displacement vectors are shown in Figure 3(a) and (b), respectively. The critical FS is calculated to be 1.34.

Figure 2.

Undeformed mesh of a homogeneous slope with a foundation layer.

Figure 3.

(a) Deformed mesh of a homogeneous slope with a foundation layer; (b) nodal displacement vectors.

3. Classical probabilistic slope stability analysis

In this section, a homogeneous slope and an infinite slope are investigated using simple and classical probabilistic slope stability methods, including FOSM, FORM, and Monte Carlo method. These methods are illustrated one by one in detail followed by a simple example, respectively.

3.1 FOSM

The FOSM method is a relatively simple method of including the effects of variability of input variables on a resulting dependent variable [17, 18]. It is basically a formalized methodology based on a first-order Taylor series expansion. This expansion is truncated after the linear term. The modified expansion is then used, along with the first two moments of the random variable(s), to determine the values of the first two moments of the dependent variable [19, 20, 21].

Consider a function f (X, Y) of two random variables X and Y. The Taylor series expansion of the function about the mean values (μ_X, μ_Y) gives

$$f\left(X,Y\right)\approx f\left({\mu }_X,{\mu }_Y\right)+\left(X-{\mu }_X\right)\frac{\partial f}{\partial x}+\left(Y-{\mu }_Y\right)\frac{\partial f}{\partial y}$$E3

where derivatives are evaluated at (μ_X, μ_Y).

To a first order of accuracy, the expected value of the function is given by

$$E\left[f\left(X,Y\right)\right]\approx f\left(E\left[X\right],E\left[Y\right]\right)$$E4

and the variance by

$$\mathit{Var}\left[f\left(X,Y\right)\right]\approx \mathit{Var}\left[\left(X-{\mu }_X\right)\frac{\partial f}{\partial x}+\left(Y-{\mu }_Y\right)\frac{\partial f}{\partial y}\right]$$E5

Hence,

$$\mathit{Var}\left[f\left(X,Y\right)\right]\approx {\left(\frac{\partial f}{\partial x}\right)}^2\mathit{Var}\left[X\right]+{\left(\frac{\partial f}{\partial y}\right)}^2\mathit{Var}\left[Y\right]+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\mathit{Cov}\left[X,Y\right]$$E6

where E[X] and E[Y] are the expected values of X and Y, respectively; Var[X] and Var[Y] are the variances of X and Y, respectively; Cov[X,Y] is the covariance of X and Y, and Cov[X,Y] = E[(X-E[X])(Y-E[Y])].

If X and Y are uncorrelated,

$$\mathit{Var}\left[f\left(X,Y\right)\right]\approx {\left(\frac{\partial f}{\partial x}\right)}^2\mathit{Var}\left[X\right]+{\left(\frac{\partial f}{\partial y}\right)}^2\mathit{Var}\left[Y\right]$$E7

In general, for a function of n uncorrelated random variables, the FOSM method tells us that

$$\mathit{Var}\left[f\left(X_1,X_2,\dots ,X_n\right)\right]\approx \sum _{i=1}^n{\left(\frac{\partial f}{\partial x_i}\right)}^2\mathit{Var}\left[X_i\right]$$E8

where the first derivatives are evaluated at the mean values (μ_X1, μ_X2,. .., μ_Xn).

Here is another example on the homogeneous slope presented in Section 2.3; a probabilistic analysis using FOSM is investigated on this slope. The shear strength parameters are as follows:

$$\begin{array}{l}{\mu }_{\varphi \text{'}}=20°,{\sigma }_{\varphi \text{'}}=3°\\ {\mu }_{c\text{'}}=10\mathit{kN}/m^2,{\sigma }_{c\text{'}}=3.0\mathit{kN}/m^2\end{array}$$

According to Eqs. 4 and 7, the expect and variance of FS can be expressed as

$$E\left[\mathit{FS}\right]\approx \mathit{FS}\left({\mu }_{\varphi \text{'}},{\mu }_{c\text{'}}\right)$$E9

$$\mathit{Var}\left[\mathit{FS}\right]\approx {\left(\frac{\partial \left(\mathit{FS}\right)}{\partial \varphi \text{'}}\right)}^2\mathit{Var}\left[\varphi \text{'}\right]+{\left(\frac{\partial \left(\mathit{FS}\right)}{\partial c\text{'}}\right)}^2\mathit{Var}\left[c\text{'}\right]$$E10

Using a central difference estimate of the derivatives with perturbations of ±σ, then

$$\mathit{Var}\left[\mathit{FS}\right]\approx {\left(\frac{\mathrm{\Delta }{\mathit{FS}}_{\varphi \text{'}}}{2}\right)}^2+{\left(\frac{\mathrm{\Delta }{\mathit{FS}}_{c\text{'}}}{2}\right)}^2$$E11

where

$$\begin{array}{l}\Delta {\mathit{FS}}_{\varphi \text{'}}=\mathit{FS}\left({\mu }_{\varphi \text{'}}+{\sigma }_{\varphi \text{'}},{\mu }_{c\text{'}}\right)-\mathit{FS}\left({\mu }_{\varphi \text{'}}-{\sigma }_{\varphi \text{'}},{\mu }_{c\text{'}}\right)\\ \Delta {\mathit{FS}}_{c\text{'}}=\mathit{FS}\left({\mu }_{\varphi \text{'}},{\mu }_{c\text{'}}+{\sigma }_{c\text{'}}\right)-\mathit{FS}\left({\mu }_{\varphi \text{'}},{\mu }_{c\text{'}}-{\sigma }_{c\text{'}}\right)\end{array}$$E12

Using program SLOPE64, FS calculated for each case is shown in Table 3.

	ϕ′	c′	FS
μϕ′, μc′	20.0	10.0	1.34	μFS = 1.34
μϕ′ + σϕ′, μc′	23.0	10.0	1.50	ΔFSϕ′ = 0.3
μϕ′ - σϕ′, μc′	17.0	10.0	1.20
μϕ′, μc′ + σc′	20.0	13.0	1.48	ΔFSc′ = 0.28
μϕ′, μc′ - σc′	20.0	7.0	1.20

Table 3.

Factor of safety for five cases.

So, the variance of FS can be calculated by

$$\mathit{Var}\left[\mathit{FS}\right]={\left(\frac{\Delta {\mathit{FS}}_{\varphi \text{'}}}{2}\right)}^2+{\left(\frac{\Delta {\mathit{FS}}_{c\text{'}}}{2}\right)}^2={\left(\frac{0.3}{2}\right)}^2+{\left(\frac{0.28}{2}\right)}^2=0.0421$$

Hence

$${\sigma }_{\mathit{FS}}=\sqrt{0.0421}=0.205$$

Assume that the FS probability density function is normal distribution (as shown in Figure 4).

$$p\left[\mathit{FS}<1\right]=\mathrm{\Phi }\left(\frac{1-1.34}{0.205}\right)=\mathrm{\Phi }\left(-1.66\right)=1-\mathrm{\Phi }\left(1.66\right)=1‐0.9515=0.0485\left(4.85\%\right)$$

Figure 4.

Normal distribution of FS.

Consider a “performance function” for this problem in which failure is defined when M < 0, the reliability index β in this case is given by $\beta =\frac{E\left[M\right]}{\sqrt{\mathit{Var}\left[M\right]}}$.

There are three different approaches calculating the reliability index β listed as follows.

3.1.1 Approach 1

For nonnegative loads and resistances (typical in geotechnical engineering), an alternative definition of the performance function could be

$$M=\frac{R}{Q}-1$$E13

so that failure occurs when M < 0 as before.

Once more assuming R and Q are uncorrelated, the FOSM method gives

$$E\left(M\right)=\frac{E\left[R\right]}{E\left[Q\right]}-1$$E14

$$\begin{array}{cc}\mathit{Var}\left[M\right]& \approx {\left(\frac{\partial M}{\partial R}\right)}^2\mathit{Var}\left[R\right]+{\left(\frac{\partial M}{\partial Q}\right)}^2\mathit{Var}\left[Q\right]\\ & =\frac{1}{E^2\left[Q\right]}\mathit{Var}\left[R\right]+\frac{E^2\left[R\right]}{E^4\left[Q\right]}\mathit{Var}\left[Q\right]\end{array}$$E15

Hence $\beta =\frac{{\mu }_Q\left({\mu }_R-{\mu }_Q\right)}{\sqrt{{\mu }_Q^2{\sigma }_R^2+{\mu }_R^2{\sigma }_Q^2}}$.

3.1.2 Approach 2

In the classical “resistance” versus “load” problem, the performance function can be defined as

$$M=R-Q$$E16

Assuming R and Q are uncorrelated, the FOSM method gives

$$E\left[M\right]=E\left[R\right]-E\left[Q\right]={\mu }_R-{\mu }_Q$$E17

And

$$\mathit{Var}\left[M\right]={\left(\frac{\partial M}{\partial R}\right)}^2\mathit{Var}\left[R\right]+{\left(\frac{\partial M}{\partial Q}\right)}^2\mathit{Var}\left[Q\right]={\sigma }_R^2+{\sigma }_Q^2$$E18

Hence $\beta =\frac{{\mu }_R-{\mu }_Q}{\sqrt{{\sigma }_R^2+{\sigma }_Q^2}}$ which is obviously different to Approach 1.

3.1.3 Approach 3

For nonnegative loads and resistances (typical in geotechnical engineering), an alternative definition of the performance function could be

$$M=ln\left(\frac{R}{Q}\right)$$E19

so that failure occurs when M < 0 as before.

Once more assuming R and Q are uncorrelated, the FOSM method gives

$$E\left[M\right]=ln\left(\frac{E\left[R\right]}{E\left[Q\right]}\right)=ln{\mu }_R-ln{\mu }_Q$$E20

$$\begin{array}{cc}\mathit{Var}\left[M\right]& \approx {\left(\frac{\partial M}{\partial R}\right)}^2\mathit{Var}\left[R\right]+{\left(\frac{\partial M}{\partial Q}\right)}^2\mathit{Var}\left[Q\right]\\ & =\frac{\mathit{Var}\left[R\right]}{R^2}+\frac{\mathit{Var}\left[Q\right]}{Q^2}\\ & =\frac{\mathit{Var}\left[R\right]}{E^2\left[Q\right]}+\frac{\mathit{Var}\left[Q\right]}{E^2\left[Q\right]}\\ & ={\nu }_R^2+{\nu }_Q^2\end{array}$$E21

Hence $\beta =\frac{\mathit{ln}\left({\mu }_R\right)-\mathit{ln}\left({\mu }_Q\right)}{\sqrt{{\nu }_R^2+{\nu }_Q^2}}$ which is clearly different to the results before.

Apparently, the reliability index β differs with the definition of the performance function, which is one of the major drawbacks of FOSM. Also, the method takes no account of the form of the probability density function describing the random variables, using only their mean and standard deviation, which ignores the effect of distribution of random variables to the results.

3.2 FORM

The major drawback to the FOSM method when used to compute probabilities relating to failure is that it can give different failure probabilities for the same problem [19, 22], which caused Hasofer and Lind to develop an improved approach, FORM [23].

As shown before, the reliability index β is given as

$$\beta =\frac{E\left[M\right]}{\sqrt{\mathit{Var}\left[M\right]}}$$E22

which measures how far the mean of the safety margin M is from zero (assumed to be the failure point) in units of number of standard deviations. The interesting point is on the probability that failure occurs, that is, M < 0. Therefore, a unique relationship between the reliability index (β) and the probability of failure (p_f) is given by

$$p_f=1-\mathrm{\Phi }\left(\beta \right)$$E23

where Φ is the standard normal cumulative distribution function. The point, line, or surface defined by M = 0 is called the failure surface.

The inconsistency of the FOSM method is due to that different definitions of margin M may have different mean estimates and different first derivatives. What the FOSM method does is calculating the distance from the average point to the failure surface in the gradient direction of the average point [18]. Hasofer and Lind solved the inconsistent problem by looking for the overall minimum distance between the average point and the failure surface, rather than just along the gradient direction [23].

In the general case, suppose that the safety margin M is a function of a sequence of random variables $X^T=\left\{X_1,X_2,\dots \right\}$, that is, $M=f\left(X_1,X_2,\dots \right)$, and that the random variables X₁, X₂, . . . have covariance matrix C. Then, the Hasofer-Lind reliability index β is defined by [23].

$$\beta =\underset{M=0}{min}\sqrt{{\left(x-E\left[X\right]\right)}^T{C}^{-1}\left(x-E\left[X\right]\right)}$$E24

which is the minimum distance between the failure surface (M = 0) and the mean point (E [X]) in units of number of standard deviations. For example, if M = f(X), then Eq. (24) simplifies to $\beta ={\mathit{min}}_x\left(x-{\mu }_X\right)/{\sigma }_X$. It is an iterative process to find β under this definition. On the curve M = 0, choose a value of x₀ and compute β₀, choose another point x₁ on M = 0 and compute β₁, and so on. The Hasofer-Lind reliability index is the minimum of all possible values of β_i. When the minimum reliability index β is determined, the probability of failure can be calculated by Eq. (23).

Figure 5 gives an example for an infinite slope. In this example, H = 5 m, γ = 20 kN/m³, α = 30°, c′ and tanϕ′ are lognormal random variables with μ_c′ = 10 kPa, σ_c′ = 3 kPa (ν_c′ = 0.3) and ϕ′ = 30°, μ_tanϕ′ = 0.5774, σ_tanϕ′ = 0.1732 (ν_tanϕ′ = 0.3); the logarithmic normal distributions of c′ and tanϕ′ are shown in Figure 6.

Figure 5.

Infinite slope.

Figure 6.

Logarithmic normal distributions of (a) c′ and (b) tanϕ′.

FS for this slope can be expressed as follows [24]:

$$\mathit{FS}=\frac{c\text{'}}{\mathit{\gamma H}sin\alpha cos\alpha }+\frac{tan\varphi \text{'}}{tan\alpha }$$E25

where H is the height of the slope, γ is the saturated unit weight, α is the slope angle to the horizontal direction, c′ is the effective cohesion, and ϕ′ is the effective friction angle.

Using Eq. (25), it can be calculated that $\overline{\mathit{FS}}=1.23$. Further, according to FORM algorithm, $\beta =0.835,p_f=20.20\%$.

In practical applications, there are many different complex optimization algorithms, usually involving the gradient of M, which can find the point where the failure plane is perpendicular to the origin. The distance between these two points is β [25, 26]. Now, many spreadsheet programs include algorithms that allow users to specify only the minimum equations and constraints on the solution. Unfortunately, nonlinear failure surfaces can sometimes have multiple local minima, with respect to the mean point, which further complicates the problem. In this case, techniques such as simulated annealing may be necessary, but which still do not guarantee finding the global minimum. Monte Carlo simulation is an alternative means of computing failure probabilities which is simple in concept. Furthermore, it is not limited to first order and can be extended easily to very difficult failure problems with only a penalty in computing time to achieve a high level of accuracy [16].

3.3 Monte Carlo method

The Monte Carlo method is a broad computational algorithm that relies on repeated random sampling to obtain numerical results. The basic concept is to use random numbers (sometimes pseudo-random numbers) to solve problems that might be deterministic in principle. This method was proposed in the 1940s and has been widely used in slope stability probability analysis [12, 27, 28, 29].

The idea of the Monte Carlo method is to randomly generate samples according to an input probability density function and evaluate the model response of each sample by a deterministic computational model. Consider the problem of determining the probability of failure of a system which has two random inputs, X₁ and X₂. The response of the system to these inputs is then defined as a function g (X₁, X₂). Obviously, the function g (X₁, X₂) is also random because the input variables are random. Assume that system failure will occur when g(X₁, X₂) > g_crit, where g_crit represents the critical state. In the space of (X₁, X₂) values, there will be some region in which g (X₁, X₂) > g_crit, and the problem boils down to assessing the probability that the particular (X₁, X₂) which actually occurs will fall into the failure region. So the probability p_f can be defined as

$$p_f=P\left[g\left(X_1,X_2\right)>g_{\mathit{\text{crit}}}\right]$$E26

Figure 7 shows the algorithm for Monte Carlo analysis of slope stability.

Figure 7.

Algorithm for Monte Carlo analysis of slope stability.

Consider the same infinite slope given in Figure 5, c′ and tanϕ′ are lognormal random variables with μ_c′ = 10 kPa, σ_c′ = 3 kPa (ν_c′ = 0.3) and ϕ′ = 30°, μ_tanϕ′ = 0.5774, σ_tanϕ′ = 0.1732 (νt_anϕ′ = 0.3), which are the same with the previous example. It can be calculated that $\overline{\mathit{FS}}=1.23,p_f=23.6\%$. Compared with the probability of failure calculated by FORM, p_f calculated using the Monte Carlo method is a little higher.

4. RFEM slope stability analysis

In this part, a new parameter spatial correlation and the local averaging method are illustrated first. After that, random finite-element method is presented. Finally, results from a full RFEM method are analyzed. Throughout this section, the probability of failure (p_f) is compared with the traditional FS that would be obtained from charts or classical limit equilibrium methods.

4.1 Spatial correlation

In probabilistic slope stability study, the shear strength c and ϕ are assumed to be characterized statistically by a normal distribution or lognormal distribution defined by means μ_c and μ_tanϕ and standard deviations σ_c and σ_tanϕ. The probability of the strength that is less than a given value can be found from standard normal distribution table. When the variables are characterized by lognormal distribution, the lognormal can be transformed to normal as follows (take c for example):

$$P\left[c<a\right]=P\left[lnc<lna\right]=P\left[Z<\frac{lna-{\mu }_{lnc}}{{\sigma }_{lnc}}\right]=\mathrm{\Phi }\left(\frac{lna-{\mu }_{lnc}}{{\sigma }_{lnc}}\right)$$E27

The lognormal parameters μ_lnc and σ_lnc given μ_c and σ_c are obtained via the transformations:

$$\begin{array}{l}{\sigma }_{lnc}^2=ln\left(1+{\nu }_c^2\right)\\ {\mu }_{lnc}=ln\left({\mu }_c\right)-\frac{1}{2}{\sigma }_{lnc}^2\end{array}$$E28

in which the coefficient of variation of c, ν_c, is defined as

$${\nu }_c=\frac{{\sigma }_c}{{\mu }_c}$$E29

Unlike the former simulation, another parameter, the spatial correlation length θ_c or θ_lnc, will be considered in the following study. The spatial correlation length describes the significant correlation distance between spatially random values in the Gaussian field. Thus, a small value of θ refers to a ragged field, while a large value refers to a smooth field. In practice, the spatial correlation length can be estimated from a set of shear strength data (c and ϕ) taken over some spatial region simply by performing the statistical analyses on the data.

It has been suggested that typical ν_c values for undrained shear strength lie in the range of 0.1–0.5. The spatial correlation length, however, is less well documented and may well exhibit anisotropy, especially when soils are typically horizontally layered. To simplify in this chapter, the spatial correlation will be assumed to be isotropic.

4.2 Local averaging

The local average subdivision (LAS) method is a fast and accurate method that produces realizations of a discrete local average random process [30]. Consider a random process Z; Table 4 presents the local average procedure via the LAS method.

Stage 0

$Z_1^0$

Stage 1

$Z_1^1$

$Z_2^1$

Stage 2

$Z_1^2$

$Z_2^2$

$Z_3^2$

$Z_4^2$

Stage 3

$Z_1^3$

$Z_2^3$

$Z_3^3$

$Z_4^3$

$Z_5^3$

$Z_6^3$

$Z_7^3$

$Z_8^3$

Stage 4

Table 4.

Procedure of the LAS method.

The algorithm proceeds as follows:

1. Generate a normally distributed random number $Z_1^0$ with mean zero; the variance is obtained from the random field.

2. Subdivide $Z_1^0$ into two equal parts, $Z_1^1$ and $Z_2^1$; the means and variances should be satisfied with three criteria:

   1. Their variances meet the requirements of local averaging theory.

   2. The relationship between $Z_1^1$ and $Z_2^1$ meets the requirements of local averaging theory.

   3. The means of $Z_1^1$ and $Z_2^1$ are equal to the mean of $Z_1^0$, that is, $Z_1^0=\frac{1}{2}\left(Z_1^1+Z_2^1\right)$.

3. Subdivide each cell in stage 1 into another two equal parts; the means and variances should be satisfied with the above three criteria, and another new requirement, $Z_1^2$ and $Z_2^2$, should be properly correlated with $Z_3^2$ and $Z_4^2$.

4. The above steps are repeated, and the cell is subdivided gradually until the size of the subunit reaches the expected requirement.

Using the RFEM approach to analyze a slope, each element is assigned a constant property, including the mean, standard deviation, and spatial correlation length of the shear strength, at each realization of the Monte Carlo process. The assigned property represents an average over the area of each finite element used to discretize the slope. If the point distribution is normal, local arithmetic averaging is used which results in a reduced variance but the mean is unaffected. In a lognormal distribution, however, local geometric averaging is used, and both the mean and the standard deviation are reduced by this form of averaging as is appropriate for situations in which low-strength regions dominate the effective strength. The reduction in both the mean and standard deviation is from

$${\mu }_X=E\left[X\right]=e^{{\mu }_{lnX}+\frac{1}{2}{\sigma }_{lnX}^2}$$E30

$${\sigma }_X^2=\mathit{Var}\left[X\right]={\mu }_X^2\left(e^{{\sigma }_{lnX}^2}-1\right)$$E31

The mean of a lognormally random variable depends on both the mean and the variance of the underlying normal log variable:

$${\sigma }_{lnX}^2=ln\left(1+\frac{{\sigma }_X^2}{{\mu }_X^2}\right)$$E32

$${\mu }_{lnX}=ln\left({\mu }_X\right)-\frac{1}{2}{\sigma }_{lnX}^2$$E33

Obviously, local averaging has a great influence on the form of a reduced mean and standard deviation. These adjustments are fully accounted for in the following RFEM analysis.

4.3 Random finite-element method

A powerful and general method of accounting for spatially random shear strength parameters and spatial correlation is the RFEM, which combines elastoplastic finite-element analysis with random field theory generated using the LAS method. Figure 8 shows a typical finite-element mesh for the test problem considered in this section. Most of the elements are square, and the elements adjacent to the slope are degenerated into triangles. Taking full account of element size in the local averaging process, the random field of shear strength values was generated and mapped onto the finite-element mesh. In a random field, the value assigned to each finite element is a random variable. The random variables can be correlated to one another by controlling the spatial correlation length θ_lnc as described previously. Figure 9a,b, and c shows the typical meshes corresponding to different spatial correlation lengths. Figure 9a shows a relatively low spatial correlation length of θ = 1, Figure 9b shows a medium spatial correlation length of θ = 5, and Figure 9c shows a relatively high spatial correlation length of θ = 10. In these figures, light regions represent weak- or low-strength soils, while dark regions represent strong- or high-strength soils. The shear strength distributions of these three cases come from the same lognormal distribution, and the only difference is the spatial correlation length. The slope stability analyses use the Tresca failure criterion which is an elastic-perfectly plastic stress–strain law. When the stresses exceed the yield stress, the program attempts to redistribute excess stresses to neighboring elements that still have reserves of strength. This process is iterative and will continue until the Tresca criterion and global equilibrium are satisfied at all points within the mesh under quite strict tolerances. Plastic stress redistribution is accomplished using a viscoplastic algorithm with eight-node quadrilateral elements and reduced integration in both the stiffness and stress redistribution parts of the algorithm [5, 6].

Figure 8.

Undeformed mesh of a homogeneous slope with a foundation layer.

Figure 9.

Deformed mesh at slope failure for three different spatial correlation lengths. (a) θ = 1; (b) θ = 5; (c) θ = 10.

4.4 Results of RFEM analyses

Figure 9 shows three typical random field realizations and the associated failure mechanisms for slopes with θ = 1, 5, and 10. It can be concluded that spatial correlation length has a great influence on the failure surface morphology. When θ is low, the shear strength between neighbored elements varies severely; when θ is high, similar properties can be found between neighbored elements. In the RFEM approach, the failure mechanism is free to seek out the weakest path through the soil. Thus, the failure surface will tend to pass through elements with weaker shear strength.

In the following part, the two parameters of shear strength, c and ϕ, are defined as the random variable, respectively, to investigate the influence of spatial correlation length and coefficient of variance on the probability of failure.

4.4.1 Define c as random

Defining friction angle as a deterministic parameter, ϕ = 20°, and then fixing the mean of cohesion with μ_c = 10 kPa, Figure 10 shows the probability of failure p_f as a function of the spatial correlation length θ for a range of coefficients of variation, with the mean cohesion fixed at μ_c = 10 kPa. Figure 11 shows the relationship between probability of failure p_f and the coefficient of variation ν_c with two different spatial correlation lengths. It can be seen from Figure 10 that the probability of failure can be divided into two branches, with the probability of failure tending to unity or zero for higher and lower values of ν_c, respectively. Figure 11 demonstrates that when θ becomes large, the probability of failure is overestimated (conservative) when the coefficient of variation is relatively small and underestimated (unconservative) when the coefficient of variation is relatively high. The RFEM results show that the inclusion of spatial correlation and local averaging in this case will always lead to a smaller probability of failure.

Figure 10.

Probability of failure versus spatial correlation length (the mean of cohesion is fixing at μ_c = 10 kPa).

Figure 11.

Probability of failure versus coefficient of variance (the mean of cohesion is fixing at μ_c = 10 kPa).

4.4.2 Define ϕ as random

Defining cohesion as a deterministic parameter, c = 10 kPa, and then fixing the mean of friction angle with μ_ϕ = 20°, Figures 12 and 13 show the effect of the spatial correlation length θ and the coefficient of variation ν_ϕ on the probability of failure for the test problem. It is obvious that Figures 12 and 13 show similar tendency with Figures 10 and 11. Comparing Figures 11 and 13, it can be concluded that the influence of spatial correlation length of ϕ on the probability of failure is less than that of c.

Figure 12.

Probability of failure versus spatial correlation length (the mean of friction angle is fixing at μ_ϕ = 20°).

Figure 13.

Probability of failure versus coefficient of variance (the mean of friction angle is fixing at μ_ϕ = 20°).

4.4.3 Define c and ϕ as random

Defining cohesion c and friction angle ϕ as random parameters, and then fixing the mean of cohesion with μ_c = 10 kPa and the mean of friction angle with μ_ϕ = 20°, Figure 14 shows the probability of failure versus spatial correlation length with different coefficients of variance of c and ϕ. Clearly, Figure 14 shows similar tendency with Figures 10 and 12. Figure 15 shows the probability of failure p_f as a function of coefficient of variance ν_c for two different θ = 2 and 10, with the mean cohesion fixed at μ_c = 10 kPa, the mean of friction angle fixing at μ_ϕ = 20° and ν_ϕ fixing at 1. Similarly, Figure 16 shows the probability of failure p_f as a function of coefficient of variance ν_ϕ for two different θ = 2 and 10, with the mean cohesion fixed at μ_c = 10 kPa, the mean of friction angle fixing at μ_ϕ = 20° and ν_c fixing at 1. Clearly, these two figures show a similar relationship with Figures 11 and 13. It is worth noting that defining ϕ as random has an apparent influence on the probability of failure versus coefficient of variance of cohesion. From Figure 15, p_f is relatively higher than the case that only c is the only random parameter.

Figure 14.

Probability of failure versus spatial correlation length (the mean of cohesion is fixing at μ_c = 10 kPa and the mean of friction angle is fixing at μ_ϕ = 20°).

Figure 15.

Probability of failure versus coefficient of variance of cohesion (the coefficient of variance of friction angle is fixing at ν_ϕ = 1).

Figure 16.

Probability of failure versus coefficient of variance of friction angle (the coefficient of variance of friction angle is fixing at ν_c = 1).

5. Conclusion

This chapter presents a deterministic slope stability analysis based on strength reduction finite-element method first. After that, three simple probabilistic methods, including FOSM, FORM, and Monte Carlo method, are introduced to perform a simple probabilistic slope stability analysis. Finally, the RFEM approach combining random field generation techniques and finite-element methods is shown and applied to slope stability risk assessment.

The elastoplastic finite-element slope stability method makes no a priori assumptions about the shape or location of the critical failure mechanism, offering significant benefits over traditional limit equilibrium methods on slope stability analysis.

FOSM, FORM, and Monte Carlo method are relatively basic and practical probabilistic analysis methods. Based on different algorithms, the uncertainty and randomness of the soil properties, especially the mean and standard deviation, can be taken into account from different views. However, there are some deficiencies, such as limit of accuracy and time-consuming on these methods.

The RFEM approach combines finite-element slope stability method and local averaging subdivision method, which can take full account of spatial correlation and local averaging. The influence of spatial correlation length and coefficient of variance on the probability of failure can be studied using a parametric approach. In the elastoplastic RFEM, the failure mechanism is free to seek out the weakest path through the soil, which leads to higher probabilities of failure than that conducted by finite-element local averaging alone.

In summary, simplified probabilistic analysis in which spatial variability is ignored can lead to unconservative estimates of the probability of failure, while the RFEM approach that considers spatial correlation and local averaging would be a practical method on slope stability risk assessment.

Acknowledgments

The author wishes to acknowledge the support from Professor D.V. Griffiths for his supervision during the period of the author’s visiting scholar at Colorado School of Mines. The author also acknowledges the support of the National Key Research and Development Program of China Grant No. 2018YFC1508602 and National Natural Science Foundation of China Grant No. 51539006.