Six parameters for soil model.
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.
- slope stability
- finite element
- probabilistic methods
- dam failure
- risk assessment
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 . 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 . 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.
The Mohr-Coulomb failure criterion used in this program can be written as follows:
where and 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 cf′ and ϕf′ are therefore given by
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
|20°||10 kPa||0||10,000 kN/m2||0.3||20 kN/m3|
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.
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.
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
where derivatives are evaluated at (μX, μY).
To a first order of accuracy, the expected value of the function is given by
and the variance by
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,
In general, for a function of n uncorrelated random variables, the FOSM method tells us that
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:
According to Eqs. 4 and 7, the expect and variance of FS can be expressed as
Using a central difference estimate of the derivatives with perturbations of ±σ, then
Using program SLOPE64, FS calculated for each case is shown in Table 3.
|μϕ′, μ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|
So, the variance of FS can be calculated by
Assume that the FS probability density function is normal distribution (as shown in Figure 4).
Consider a “performance function” for this problem in which failure is defined when M < 0, the reliability index β in this case is given by .
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
so that failure occurs when M < 0 as before.
Once more assuming R and Q are uncorrelated, the FOSM method gives
3.1.2 Approach 2
In the classical “resistance” versus “load” problem, the performance function can be defined as
Assuming R and Q are uncorrelated, the FOSM method gives
Hence 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
so that failure occurs when M < 0 as before.
Once more assuming R and Q are uncorrelated, the FOSM method gives
Hence 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.
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 .
As shown before, the reliability index β is given as
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 (pf) is given by
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 . 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 .
In the general case, suppose that the safety margin M is a function of a sequence of random variables , that is, , and that the random variables X1, X2, . . . have covariance matrix C. Then, the Hasofer-Lind reliability index β is defined by .
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 . It is an iterative process to find β under this definition. On the curve M = 0, choose a value of x0 and compute β0, choose another point x1 on M = 0 and compute β1, 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/m3, α = 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.
FS for this slope can be expressed as follows :
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 . Further, according to FORM algorithm, .
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 .
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, X1 and X2. The response of the system to these inputs is then defined as a function g (X1, X2). Obviously, the function g (X1, X2) is also random because the input variables are random. Assume that system failure will occur when g(X1, X2) > gcrit, where gcrit represents the critical state. In the space of (X1, X2) values, there will be some region in which g (X1, X2) > gcrit, and the problem boils down to assessing the probability that the particular (X1, X2) which actually occurs will fall into the failure region. So the probability pf can be defined as
Figure 7 shows the 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 (νtanϕ′ = 0.3), which are the same with the previous example. It can be calculated that . Compared with the probability of failure calculated by FORM, pf 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 (pf) 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):
The lognormal parameters μlnc and σlnc given μc and σc are obtained via the transformations:
in which the coefficient of variation of c, νc, is defined as
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 . Consider a random process Z; Table 4 presents the local average procedure via the LAS method.
The algorithm proceeds as follows:
Generate a normally distributed random number with mean zero; the variance is obtained from the random field.
Subdivide into two equal parts, and ; the means and variances should be satisfied with three criteria:
Their variances meet the requirements of local averaging theory.
The relationship between and meets the requirements of local averaging theory.
The means of and are equal to the mean of , that is, .
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, and , should be properly correlated with and .
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
The mean of a lognormally random variable depends on both the mean and the variance of the underlying normal log variable:
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].
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 pf 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 pf 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.
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.
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 pf 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 pf 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, pf is relatively higher than the case that only c is the only random parameter.
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.
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.