Computers and Geotechnics 37 (2010) 1015–1022
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Efficient Monte Carlo Simulation of parameter sensitivity in probabilistic slope stability analysis Yu Wang ⇑, Zijun Cao 1, Siu-Kui Au 2 Department of Building and Construction, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong
a r t i c l e
i n f o
Article history: Received 21 April 2010 Received in revised form 12 July 2010 Accepted 23 August 2010 Available online 18 September 2010 Keywords: Probabilistic failure analysis Slope stability Monte Carlo Simulation Subset Simulation Hypothesis tests Bayesian analysis
a b s t r a c t Monte Carlo Simulation (MCS) method has been widely used in probabilistic analysis of slope stability, and it provides a robust and simple way to assess failure probability. However, MCS method does not offer insight into the relative contributions of various uncertainties (e.g., inherent spatial variability of soil properties and subsurface stratigraphy) to the failure probability and suffers from a lack of resolution and efficiency at small probability levels. This paper develop a probabilistic failure analysis approach that makes use of the failure samples generated in the MCS and analyzes these failure samples to assess the effects of various uncertainties on slope failure probability. The approach contains two major components: hypothesis tests for prioritizing effects of various uncertainties and Bayesian analysis for further quantifying their effects. Equations are derived for the hypothesis tests and Bayesian analysis. The probabilistic failure analysis requires a large number of failure samples in MCS, and an advanced Monte Carlo Simulation called Subset Simulation is employed to improve efficiency of generating failure samples in MCS. As an illustration, the proposed probabilistic failure analysis approach is applied to study a design scenario of James Bay Dyke. The hypothesis tests show that the uncertainty of undrained shear strength of lacustrine clay has the most significant effect on the slope failure probability, while the uncertainty of the clay crust thickness contributes the least. The effect of the former is then further quantified by a Bayesian analysis. Both hypothesis test results and Bayesian analysis results are validated against independent sensitivity studies. It is shown that probabilistic failure analysis provides results that are equivalent to those from additional sensitivity studies, but it has the advantage of avoiding additional computational times and efforts for repeated runs of MCS in sensitivity studies. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Various uncertainties exist in slope engineering, such as inherent spatial variability of soil properties, subsurface stratigraphy, simplifications and approximations adopted in geotechnical models. Effects of these uncertainties on probability of slope failure are often significant, and insight on these effects is of great value for understanding failure mechanisms and designing slope remedial measures. Several probabilistic methodologies have been developed to incorporate these uncertainties in slope stability analysis, such as the First Order Second Moment (FOSM) method [1–4], First Order Reliability Method (FORM) [5–9], Monte Carlo Simulation (MCS) method [10–20], and other methods based on artificial neural network, support vector machine or random set
⇑ Corresponding author. Tel.: +852 2788 7605; fax: +852 2788 7612. E-mail addresses:
[email protected] (Y. Wang),
[email protected]. edu.hk (Z. Cao),
[email protected] (S.-K. Au). 1 Tel.: +852 3442 6492; fax: +852 2788 7612. 2 Tel.: +852 2194 2769; fax: +852 2788 7612. 0266-352X/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2010.08.010
[21–25]. Among these methods, MCS method that integrates with different deterministic analysis methods (e.g., limit equilibrium methods [10,11,13], finite element methods [12,14,15]) are gaining popularity due to their robustness and conceptual simplicity. However, as pointed out by Baecher and Christian [26], although MCS method provides a robust and simple way to assess failure probability, it does not offer insight into the relative contributions of various uncertainties to the failure probability and suffers from a lack of resolution and efficiency at small probability levels. This paper develops a probabilistic failure analysis approach that makes uses of the failure samples generated in the MCS and analyzes these failure samples to assess the effects of various uncertainties on slope failure probability. An advanced Monte Carlo Simulation called Subset Simulation (Subsim) [27,28] is employed to improve efficiency of generating failure samples in MCS and resolution at small failure probability levels. The paper starts with mathematical formulation of the approach, including hypothesis tests for prioritizing effects of various uncertainties and Bayesian analysis for further quantifying their effects. Subset Simulation and its implementation in a commonly-available spreadsheet environment are briefly described. As an illustration,
1016
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
the proposed approach is applied to study a design scenario of James Bay Dyke [2,10,14]. 2. Probabilistic failure analysis approach Probabilistic failure analysis is similar to back-analysis [29–31], which is a common analysis procedure in geotechnical engineering that intends to find a set of model parameters that would result in the observed performance of geo-structures. Similarly, probabilistic failure analysis aims to identify a group of uncertain parameters that would significantly affect the slope performance (i.e., the probability of slope failure). The back-analysis, however, relies on the observed performance and it is inapplicable when observed performance of the geo-structures of interest is unavailable (e.g., during design analysis of a new geo-structure). On the other hand, probabilistic failure analysis makes uses of the failure samples generated in the MCS, and it is readily applicable in design analysis for evaluating the effects of various uncertainties. The proposed probabilistic failure analysis approach contains two major components: hypothesis tests for prioritizing effects of various uncertainties and Bayesian analysis for further quantifying their effects, which are described in the following two subsections, respectively. 2.1. Hypothesis tests The effects of various uncertainties on the probability of slope failure are prioritized by comparing, statistically, failure samples with their respective nominal (unconditional) samples. When the uncertainty of an uncertain system parameter has a significant effect on the probability of slope failure, the mean l of failure samples of the parameter differs significantly from the mean l0 of its unconditional samples. The statistical difference between l and l0 is evaluated by hypothesis tests. A null hypothesis H0 and alternative hypothesis HA are defined as [32]
H0 : l ¼ l0 HA : l–l0
ð1Þ
Then, a hypothesis test statistic ZH of the parameter is formulated as
ZH ¼
l l0 pffiffiffi r= n
ð2Þ
where r is standard deviation of the uncertain parameter and n is the number of failure samples. Based on Central Limit Theorem, ZH follows the standard Normal distribution when n is large (e.g., n P 30) [32]. When the failure sample mean l deviates statistically from the unconditional mean l0 of the parameter, the absolute value of ZH is relatively large. As the absolute value of ZH increases, the statistical difference between l and l0 becomes growingly significant, and the effect of the uncertain parameter on failure probability also becomes growingly significant. The absolute value of ZH is therefore formulated in this paper as an index to measure the effects of the uncertain parameters on failure probability and to prioritize their relative effects on failure probability. Using the absolute value of ZH, the uncertain parameters that have significant effects on failure probability are selected, and their effects are further quantified using a Bayesian analysis approach described in the next subsection. 2.2. Bayesian analysis The failure samples generated in the MCS are further analyzed by a Bayesian analysis to quantify effects of various uncertainties. Let h denote an uncertain parameter selected based on hypothesis tests. In the context of the Bayesian Theorem [33,34]
PðFjhÞ ¼ PðFÞPðhjFÞ=PðhÞ
ð3Þ
where P(F|h) is the conditional probability density function (PDF) of slope failure for a given h value; P(F) is the probability of slope failure; P(h|F) is the conditional PDF of h given that the slope has failed; and P(h) is the unconditional PDF of h. As both P(F) and P(h|F) are estimated from failure samples of MCS and P(h) is given before MCS, Eq. (3) can be used to estimate P(F|h) using P(h) and P(h|F) obtained from analysis of failure samples. Note that P(F|h) is a variation of failure probability as a function of h, and it can be considered as results of a sensitivity study of h on slope failure probability. In other words, the probabilistic failure analysis approach presented in this paper, which makes use of the failure samples generated in a single run of MCS for assessment of failure probability, provides results that are equivalent to those from a sensitivity study, which frequently includes many repeated runs of MCS with different given values of h in each run. Additional computational times and efforts for repeated runs of MCS in the sensitivity study can be avoided using the probabilistic failure analysis approach described herein. In addition, Eq. (3) implies that comparison between the conditional probability P(h|F) and its unconditional one P(h) provides an indication of the effect of the uncertain parameter h on failure probability. In general, P(F|h) changes as the values of the uncertain parameter h changes. However, when P(h|F) is similar to P(h), P(F|h) remains more or less constant regardless of the values of h. This implies that the effect of h on the slope failure probability is minimal. Such implication can be used to validate the prioritization obtained from hypothesis tests, as shown in the example of James Bay Dyke later. The resolution of P(F) and P(h|F) is pivotal to obtain P(F|h), and it depends on the number of failure samples generated in MCS. As the number of failure samples increases, the resolution improves. For a given slope stability problem, the value of P(F) is constant, although unknown before MCS. In this case, increasing the number of failure samples necessitates an increase in the total number of samples in MCS. One possible way to improve the resolution is, therefore, to increase the total number of samples in MCS at the expense of computational time. Alternatively, advanced MCS methods can be employed to improve efficiency and resolution at small failure probability levels. An advanced MCS called Subset Simulation (Subsim) [27,28,35,36] is used in this paper to calculate the failure probability and generate failure samples efficiently for probabilistic failure analysis. 3. Subset Simulation 3.1. Algorithm Subset Simulation is an adaptive stochastic simulation procedure for efficiently generating failure samples and computing small tail probability [27,28]. It expresses a small probability event as a sequence of intermediate events {F1, F2, . . ., Fm} with larger conditional probability and employs specially designed Markov chains to generate conditional samples of these intermediate events until the final target failure region is achieved. For the slope stability problem, the factor of safety (FS) is the key parameter, and let Y = 1/FS be the critical response [20,35]. The probability of Y = 1/FS larger than a given value y (i.e., P(Y = 1/FS > y)) is of interest and let 0 < y1 < y2 < < ym1 < ym = y be an increasing sequence of intermediate threshold values. The sequence of intermediate events {F1, F2, . . ., Fm} are chosen as Fi = {Y > yi, i = 1, 2, . . ., m} for these intermediate threshold values. By sequentially conditioning on the event {Fi, i = 1, 2, . . ., m}, the failure probability can be written as
PðFÞ ¼ PðF m Þ ¼ PðF 1 Þ
m Y
PðF i jF i1 Þ
ð4Þ
i¼2
where P(F1) is equal to P(Y > y1) and P(Fi|Fi1) is equal to {P(Y > yi|Y > yi1): i = 2, . . ., m}. In implementations, y1, y2, . . ., ym
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
are generated adaptively using information from simulated samples so that the sample estimate of P(F1) and {P(Fi|Fi1): i = 2, . . ., m} always corresponds to a common specified value of conditional probability p0 (p0 = 0.1 is found to be a good choice) [35,36]. The efficient generation of conditional samples is pivotal in the success of Subset Simulation, and it is made possible through the machinery of Markov Chain Monte Carlo (MCMC). In MCMC, Metropolis algorithm [37] is used and successive samples are generated from a specially designed Markov chain whose limiting stationary distribution tends to the target PDF as the length of Markov chain increases. Details of the computational process of Subset Simulation and its application to slope stability problems are referred to Au and Beck [27,28] and Wang et al. [20], respectively. The Subset Simulation provides much more failure samples than Direct MCS under the same total number of samples, especially when the failure probability is relatively small. When compared with failure samples generated in direct MCS where each sample carries equal weight in the calculation of P(F) and P(h|F), the samples generated by Subset Simulation are conditional samples and carry different weights for different intermediate events Fm. Thus, when using these conditional failure samples collected from Subset Simulation to construct the conditional PDF P(h|F) required in Eq. (3), a weighted summation by Total Probability Theorem [33,34] is necessary, which is described in the following subsection. 3.2. Estimation of P(h|F) based on conditional failure samples Consider a Subset Simulation that performs m + 1 levels of simulations. The first level of Subset Simulation is direct MCS, and samples of the next level are generated conditional on the samples collected from the previous level. The intermediate threshold values {yi, i = 1, 2, . . ., m} divide the sample space X of an uncertain parameter h into m individual sets {Xi, i = 0, 1, 2, . . ., m}. According to the Total Probability Theorem [34], the failure probability can be written as
PðFÞ ¼
m X
PðFjXi ÞPðXi Þ
ð5Þ
i¼0
where X0 = {Y 6 y1}; Xi, i = 1, . . ., m 1 is equal to Fi Fi+1 (i.e., Xi = {yi 6 Y 6 yi+1}); Xm is equal to Fm (i.e., Xm = {Y P ym}); P(F|Xi) is the conditional failure probability given sampling in Xi; P(Xi) is the probability of the event Xi. P(F|Xi) is estimated as the fraction of the failure samples in Xi. The failure samples are collected from samples generated by Subset Simulation and are based on the performance failure criteria (i.e., Y = 1/FS P 1 for a slope stability problem). P(Xi) is calculated as
PðX0 Þ ¼ 1 p0 PðXi Þ ¼ pi0 p0iþ1
ði ¼ 1; . . . ; m 1Þ
ð6Þ
PðXm Þ ¼ pm 0 P Note that P(Xi \ Xj) = 0 for i – j and m i¼0 PðXi Þ ¼ 1. When P(F), P(F|Xi) and P(Xi) are obtained, the conditional probability P(Xi|F) is calculated using the Bayesian Theorem
PðXi jFÞ ¼ PðFjXi ÞPðXi Þ=PðFÞ
ð7Þ
Then the conditional PDF P(h|F) of an uncertain parameter h is given by the Total Probability Theorem as
PðhjFÞ ¼
m X
PðhjXi \ FÞPðXi jFÞ
ð8Þ
i¼0
where P(h|Xi \ F) is the conditional probability of h estimated from failure samples that lie in Xi. In this paper, P(h|Xi \ F) is estimated
1017
from the failure sample histogram in Xi. The number of bins k in the failure sample histogram is estimated as [38] n
k ¼ 1 þ log2i
ð9Þ
where ni is the number of the failure samples in Xi. Using Eqs. (5) and (8), P(F) and P(h|F) can be calculated from Subset Simulation, and subsequently, P(F|h) is calculated using Eq. (3). 4. Implementation of Subset Simulation in spreadsheet environment The Subset Simulation described above has been implemented in a commonly-available spreadsheet environment by a package of worksheets and functions/Add-In in EXCEL with the aids of Visual Basic for Application (VBA) [35,36]. The package is divided into three parts: deterministic model worksheet, uncertainty model worksheet and Subset Simulation Add-In, which are described briefly in the following three subsections, respectively. 4.1. Deterministic model worksheet For a slope stability problem, deterministic model analysis is the process of calculating factor of safety for a given nominal set of values of system parameters. The system parameters include the geometry information of the slope and the slip surface, soil properties and profile of soil layers. In this paper, limit equilibrium methods (e.g., Swedish Circle method, Simplified Bishop method and Spencer method) [39] are employed to calculate the factor of safety for the critical slip surface. The calculation process of deterministic analysis is implemented in a series of worksheets assisted by some VBA functions/Add-In [35,36]. From an input–output perspective, the deterministic analysis worksheets take a given set of values as input, calculate the factor of safety and return the factor of safety as an output. 4.2. Uncertainty model worksheet An uncertainty model worksheet is developed to generate random samples of uncertain system parameters that are treated as random variables in the analysis. The uncertain worksheet includes detailed information of random variables, such as statistics, distribution type and correlation information. The generation of random samples starts with an EXCEL built-in function ‘‘RAND()” for generating uniform random samples, which are then transformed to random samples of the target distribution type (e.g., normal distribution or lognormal distribution). If the random variables are considered correlated, Cholesky factorization of the correlation matrix is performed to obtain a lower triangular matrix, which is used in the transformation to generate correlated random samples. Details of the random sample generation process are referred to Au et al. [36]. From the input–output perspective, the uncertainty model worksheet takes no input but returns a set of random samples of the uncertain system parameters as its output. When deterministic model worksheet and uncertainty model worksheet are developed, they are linked together through their input/output cells to execute probabilistic analysis of slope stability. The connection is carried out by simply setting the cell references for nominal values of uncertain parameters in deterministic model worksheet to be the cell references for the random samples in the uncertainty model worksheet in EXCEL. After this task, the values of uncertain system parameters shown in the deterministic model worksheet are equal to that generated in the uncertainty model worksheet, and so the values of the safety factor calculated in the deterministic modeling worksheet are random.
1018
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
4.3. Subset Simulation Add-In
5. The James Bay Dyke As an illustration, the probabilistic failure analysis approach is applied to analyze a design scenario of the James Bay Dyke. The James Bay Dyke is a 50 km-long earth dyke of the James Bay hydroelectric project in Canada. Soil properties and various design scenarios of the dyke were studied by Ladd et al. [40], Soulié et al. [41], Christian et al. [2], El-Ramly [10], El-Ramly et al. [11] and Xu and Low [14]. As shown in Fig. 2, the embankment is 12 m high with a 56 m-wide berm at mid-height. The slope angle of the embankment is about 18.4° (3H:1V). The embankment is overlying on a clay crust with a thickness Tcr. The clay crust is underlain by a layer of 8.0-m thick sensitive marine clay and a layer of lacustrine
Elevation (m)
When the deterministic analysis and uncertainty model worksheets are completed and linked together, Subset Simulation procedure is invoked for uncertainty propagation. In this paper, Subset Simulation is implemented as an Add-In in EXCEL [35,36]. The userform of the Add-In is shown in Fig. 1. The upper four input fields of the userform (i.e., number of Subset Simulation runs, number of samples per level N, conditional probability from one level to next p0, the highest Subset Simulation level m) control the number of samples generated by Subset Simulation. The total number of samples per Subset Simulation run is equal to N + mN(1 p0). The lower four input fields of the userform record the cell references of the random variables, their PDF values, and the cell references of the system response (e.g., Y = 1/FS in this paper) and other variables V (e.g., the samples generated in this paper) of interest, respectively. After each simulation run, the Add-In provides the complementary cumulative density function (CDF) of the driving variable versus the threshold level, i.e., estimate for P(Y > y) versus y = 1/FS in this paper, into a new spreadsheet and produces a plot of it [36]. Then, based on the output information obtained, the CDF, histograms or conditional counterparts (e.g., P(h|F)) of uncertain parameters of interest are calculated using the procedures and equations described in the previous section.
56.0m 45
(x, y)
3 1 Emb m ankment γFill; φFill Claay Cr Cl Cru ust 15 Marine Clay y SuM S y uL 5 Lacustrine Clay Tiill T -5 0 80 20 60 40 35 25
3
100
120
1
12.0m Tcr 8.0m DTill TL Critical slip p circle
140
160
Distance (m) Fig. 2. The James Bay Dyke (modified after [10]).
clay with a thickness TL. The undrained shear strength (i.e., SuM and SuL) of the marine clay and the lacustrine clay were measured by field vane tests [2,10,40,41]. The lacustrine clay is overlying on a stiff till layer, the depth to the top of which is DTill. Six uncertain system parameters have been considered in literature (e.g., El-Ramly [10], El-Ramly et al. [11] and Xu and Low [14]), including the friction angle /Fill and unit weight cFill of embankment material, the thickness Tcr of clay crust, the undrained shear strength SuM of the marine clay, the undrained shear strength SuL of the lacustrine clay, and the depth of the till layer DTill. During the probabilistic failure analysis of the dyke, the six uncertain parameters are represented by six independent Gaussian random variables [10], respectively. Statistics [i.e., mean, standard deviation and coefficient of variation (COV)] of these six random variables are summarized in Table 1. These statistics are used to generate random samples for each random variable in uncertain model worksheet. Note that thickness of the lacustrine clay TL is an uncertain variable that depends on Tcr and DTill and has a mean of about 6.5 m (see Fig. 2). In addition to these uncertain parameters, other system parameters are taken as deterministic, including an undrained shear strength of 41 kPa for the clay crust and unit weights of 19 kN/m3, 19 kN/m3 and 20.5 kN/m3 for the clay crust, marine clay and lacustrine clay, respectively [10]. Using the soil properties described above, El-Ramly [10] and Xu and Low [14] employed direct MCS methods integrating with limit equilibrium methods and response surface method, respectively, to evaluate failure probability of the dyke. To enable a consistent comparison with the analyses by El-Ramly [10], the critical slip surfaces recommended by El-Ramly [10] are adopted in this paper, which is circular and always tangential to top of the till layer and pass through the point (x = 4.9 m, y = 36.0 m). The x coordinate of the center is fixed at 85.9 m. For each set of random samples, the critical slip surface is specified uniquely by the value of DTill. In this paper, the safety factor of the critical slip surface is calculated by Simplified Bishop method, and two Subset Simulation runs are performed. One has the highest simulation level m = 3 and sample number N = 1000 per each level, as opposed to m = 4 and N = 10,000 per each level in the other run.
Table 1 Soil properties of the James Bay Dyke (modified after [10]). Soil layers
Uncertain system parameters*
Mean
Standard deviation
Coefficient of variation (%)
Embankment
/Fill (°) cFill (kN/m3) Tcr (m) SuM (kN/m2) SuL (kN/m2) DTill (m)
30.0 20.0 4.0 34.5 31.2 18.5
1.79 1.10 0.48 3.95 6.31 1.00
6.0 5.5 12.0 11.5 20.2 5.4
Clay crust Marine clay Lacustrine clay Till *
Fig. 1. The userform of Subset Simulation Add-In.
All parameters are modeled as Gaussian random variables. Thickness of the Lacustrine clay layer TL is an uncertain variable that depends on Tcr and DTill and has a mean of about 6.5 m.
1019
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
5.1. Simulation results Fig. 3 shows a typical complementary CDF for Y = 1/FS (i.e., P(Y > y) versus y), from two Subset Simulation runs with a total sample number NT = 1000 + 3 1000 (1–0.1) = 3700 (i.e., Subsim Run 1) and NT = 10,000 + 4 10,000 (1–0.1) = 46,000 (i.e., Subsim Run 2), respectively. For comparison, the result from a direct MCS with 20,000 samples is also plotted. Three consistent failure probabilities P(Y P 1) = 0.22%, 0.23%, and 0.25% are estimated, respectively, from direct MCS and two runs of Subset Simulations. In addition, the two runs of Subset Simulations provide results that are consistent even at low probability levels [e.g., P(Y P y) = 0.01%] where the CDF curve from direct MCS becomes erratic. Table 2 compares the simulation results with those reported by El-Ramly [10] and Xu and Low [14]. El-Ramly [10] performed direct MCS with the Simplified Bishop method (i.e., the same limit equilibrium method used in this paper) and obtained a failure probability Pf = 0.24%. This Pf value is almost identical to the average of the three Pf values obtained in this paper (i.e., 0.22%, 0.23%, and 0.25% in Columns 4–6 of Table 2). In addition, Xu and Low [14] combined MCS with response surface method to estimate the Pf of the James Bay Dyke and obtained a Pf value of 0.33%. Although different deterministic slope stability analysis methods were used, the Pf values obtained compare favourably with each other. This implies that the probabilistic analysis models for the James Bay Dyke presented in this paper work properly. However, note that the Pf values summarized in Table 2 are obtained under the assumption of circular slip surfaces. In contrast, non-circular failure mechanism has been identified in recent finite element analyses incorporating spatially varying soil properties [14,19]. When the searching of non-circular slip surfaces is included in the MCS, the corresponding Pf values will increase. Table 2 also compares the number of failure samples in direct MCS (e.g., Columns 2 or 4) with that in Subset Simulations (e.g.,
100
P(Y>=y) (%)
10
1
6. Probabilistic failure analysis results With the large number of failure samples generated from Subset Simulations, probabilistic failure analysis are performed for the James Bay Dyke, including hypothesis tests for identifying key uncertainties that have significant effects on slope failure probability and Bayesian analysis for further quantifying their effects. 6.1. Hypothesis test results Based on the failure samples generated from Subsim Run 1, the hypothesis test statistics ZH defined by Eq. (2) are calculated and shown in Fig. 4 for all uncertain parameters. The absolute values of ZH varies from less than 2 for the thickness of clay crust Tcr to about 95 for undrained shear strength of the lacustrine clay SuL. The decreasing order of the ZH absolute values is: SuL, DTill, cFill, SuM, /Fill, and Tcr. This implies that the uncertainty of SuL has the most significant effects on the slope failure probability, while the uncertainty of Tcr contributes the least to the failure probability. This result is consistent with that reported by El-Ramly [10] who employed an Excel spreadsheet – based MCS program @RISK [42] and compared the spearman rank correlation coefficients for various uncertain parameters to show that SuL and Tcr are the most and least influential parameter, respectively. The results can be validated by an independent sensitivity studies on SuL and Tcr, which are described in next subsection. 6.2. Validation of hypothesis test results
0.1 Subsim Run 1
3700 Samples
Subsim Run 2 46000 Samples
0.01
Direct MCS
0.001 0.5
Columns 5 or 6). For a total sample number NT = 20,000, direct MCS leads to only 48 or 44 failure samples. In contrast, Subset Simulations with m = 3 and NT = 3700 (i.e., Run 1 in Column 5) or m = 4 and NT = 46,000 (i.e., Run 2 in Column 6) result in a failure sample number of 1128 and 20,482, respectively. This comparison clearly shows that Subset Simulations significantly improve the efficiency of generating failure samples, which enables generation of a large number of failure samples with relative ease and makes probabilistic failure analysis feasible. Table 2 also shows that, as the value of m increases (e.g., from 3 in Column 5 to 4 in Column 6), the efficiency increases as well (e.g., the percentage of failure sample increases from 30.5% to 44.5%).
0.6
0.7
20000 Samples
0.8
0.9
1.0
1.1
1.2
1.3
y=1/FS Fig. 3. Complementary CDF plot from simulations.
Sensitivity studies are performed to explore the effect of SuL and Tcr uncertainties on slope failure probability. The coefficient of variation (COV) of SuL and Tcr is varied in the sensitivity studies from 0.15 to 0.50 and from 0.05 to 0.25, respectively. The SuL COV range adopted in this study follows the typical range of COV of undrained shear strength of clay measured by vane shear tests [43], and the COV of Tcr varies from half to about twice of the value reported by El-Ramly [10]. Other parameters, including the means of SuL and Tcr, remain unchanged in the sensitivity studies. About 40 additional Subset Simulation runs are performed to validate the hypothesis test results, and their results are shown in Fig. 5a and
Table 2 Comparison of simulation results.
a b
Simulation results
El-Ramly [10]
Xu and Low [14]
Direct MCS with EXCEL
Subsim with EXCEL Run 1
Subsim with EXCEL Run 2
Failure probability Pf (%) Number of failure samples, NF Number of Total Samples, NT Percentage of failure samples, NF/NT (%)
0.24a 48 20,000 0.24
0.33b N/A N/A N/A
0.22a 44 20,000 0.22
0.23a 1128 3700 30.5
0.25a 20,482 46,000 44.5
Pf is calculated by MCS methods integrating with the simplified Bishop method. Pf is calculated by MCS methods integrating with the response surface method.
1020
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
γ Fill φ Fill T cr S uM S uL D Till 0
10
20
30
40
50
60
70
80
90
100
Absolute Value of Z H Fig. 4. Summary of absolute values of ZH.
b for Tcr and SuL, respectively. In addition, sensitivity studies on COVs of Tcr and SuL are also carried out using the Excel spreadsheet – based MCS program @RISK [42] and a FORM calculation spreadsheet developed by Low and his coworkers [5–8]. Fig. 5 includes the results from @RISK and FORM for comparison. Fig. 5a shows that, when the COV of Tcr varies from 0.05 to 0.25, the slope failure probability fluctuates between 0.19% and 0.43%. For comparison, the baseline failure probability (i.e., about 0.24% corresponding to the values summarized in Table 1) is also included in the figure. The failure probabilities from sensitivity study using Subset Simulations fall around the horizontal line of 0.24%,
and the failure probability is insensitive to the uncertainty on Tcr. The results from Subset Simulations are in good agreement with those from @RISK and FORM. This validates the results from hypothesis tests that the Tcr uncertainty has the least effect on the failure probability. It is interesting to note that the probability of slope failure should, theoretically, increase as the COV of Tcr increases from 0.05 to 0.25, as shown by the FORM results [i.e., a rather slight increase of the Pf values shown by the open circles in Fig. 5a when the COV of Tcr increases from 0.05 to 0.25]. The effect is, however, so minimal that it is dominated by the MCS ‘‘noise” (i.e., the random fluctuations of the Pf values obtained from Subset Simulations and @RISK). On the other hand, Fig. 5b shows that the slope failure probability increases as the COV of SuL increases. The results from Subset Simulations are in good agreement with those from @RISK and FORM. When the COV of SuL is 0.15, the probability is less than 0.1%. When the COV of SuL increases to 0.50, the slope failure probability increases by two orders of magnitude (i.e., increases to about 10%). The failure probability varies significantly with the change of the SuL COV. This agrees well with the results from hypothesis tests that the uncertainty of SuL has significant effect on the slope failure probability. Such agreement further validates that the hypothesis test procedure proposed in this paper properly prioritizes effects of various uncertainties on failure probability.
6.3. Bayesian analysis results Based on the failure samples generated from Subsim Run 2, a Bayesian analysis is performed using Eqs. (3), (7), and (8) accord-
100
Subsim FORM @ Risk isk
Failure Probability (%)
Failure Probability, P(F) (% )
100
10
1 0.24%
0.1
1
0.1
0.01 Parametric P(F|SuL) estimated from sensitivity study P(F|SuL) estimated from Bayesian analysis Series4
COV=12.0% used by El-Ramly [10]
0.01 0.05
0.10
0.15
0.20
0.001
0.25
COV of T cr
12
Failure Probability (%)
10
1 0.24%
0.1
0.30
0.35
20
22
24
10
1 0.1
0.01
COV=20.2% used by El-Ramly [10]
0.25
18
100
Subsim FORM Risk @ Risk
0.20
16
SuL(kN/m ) (a) θ = undrained shear strength of the lacustrine clay, SuL
100
0.01 0.15
14
2
(a) Thickness of clay crust, Tcr
Failure Probability, P(F) (% )
10
0.40
0.45
0.50
COV of S uL (b) Undrained shear strength of the lacustrine clay, SuL Fig. 5. Effects of coefficient of variation (COV) for different system parameters.
0.001 17
Parametric P(F|DTill) estimated from sensitivity study Series4 P(F|DTill) estimated from Bayesian analysis
18
19
20
21
22
23
DTill (m) (b) θ = depth of the till layer, DTill Fig. 6. Bayesian analysis results for different system parameters h.
24
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
ingly. Fig. 6 shows the Bayesian analysis results by open circles for SuL and DTill, which have been identified from the Hypothesis tests (see Section 6.1) as the two most influential uncertain parameters. Note that the conditional probability [i.e., P(F|SuL) in Fig. 6a and P(F|DTill) in Fig. 6b] obtained from Bayesian analysis is a variation of failure probability as a function of SuL or DTill. Fig. 6a shows that, as SuL increases from 12 kPa to 24 kPa, the slope failure probability decreases from more than 10% to less than 0.1%. Similarly, Fig. 6b shows that, as DTill increases from about 18 m to 23 m, the slope failure probability increases from about 0.1% to about 10%. It is obvious that the values of SuL and DTill have significant effects on slope failure probability, and such effects can be quantified explicitly from the Bayesian analysis of failure samples.
Probability (%)
100
10
1
0.1 P(Tcr|F) Series4 Series1 P(Tcr)
0.01
2
3
4
5
6
Tcr(m) (a) θ = thickness of clay crust, Tcr
Probability (%)
100
6.4. Validation of Bayesian analysis results Note that variations of failure probability as a function of SuL or DTill shown in Fig. 6 can also be obtained from sensitivity studies on SuL or DTill, which frequently include many repeated runs of MCS with different given values of SuL or DTill in each run. Therefore, to validate the Bayesian analysis results, about 40 additional Subset Simulation runs are performed with different given values of SuL or DTill in each run. Fig. 6 also includes results from these additional sensitivity runs by open triangles. The open triangles follow trends similar to the open circles (i.e., the Bayesian analysis results) and plot closely to the open circles as well. This validates that a Bayesian analysis of failure samples generated in MCS or Subset Simulations provides results that are equivalent to those from additional sensitivity studies. In addition, the Bayesian analysis has the advantage of avoiding additional computational times and efforts for repeated runs of MCS or Subset Simulations in the sensitivity studies. As mentioned before, Eq. (3) implies that comparison between the conditional probability P(h|F) and its unconditional one P(h) provides an indication of the effect of the uncertain parameter h on failure probability. This offers a means to check the Bayesian analysis results with those from hypothesis tests. Fig. 7 compares the conditional probabilities [i.e., P(Tcr|F), P(SuL|F) and P(DTill|F)] for Tcr, SuL and DTill and their unconditional ones [i.e., P(Tcr), P(SuL) and P(DTill)]. Fig. 7a shows that P(Tcr|F) and P(Tcr) are almost identical for different Tcr values, which is consistent with the hypothesis test results that Tcr has the least effect on failure probability. In contrast, the hypothesis tests show that SuL and DTill are the two most influential uncertain parameters. Fig. 7b and c show that P(SuL|F) and P(DTill|F) differ significantly from their respective unconditional one [i.e., P(SuL) and P(DTill)]. The Bayesian analysis results agree well with the hypothesis test results.
10
7. Conclusions 1
0.1 Series4 P(S uL|F) P(SuL) Series1
0.01 12
14
16
18
20
22
24
2
SuL(kN/m ) (b) θ = undrained shear strength of the lacustrine clay, SuL 100
10
Probability (%)
1021
1
0.1
0.01 Series6 P(DTill|F) Series4 P(DTill)
0.001 17
18
19
20
21
22
23
24
D Till (m)
(c) θ = depth of the till layer, DTill Fig. 7. Conditional probability density function (PDF) (P(h|F)) for different system parameters h.
This paper developed a probabilistic failure analysis approach that makes use of the failure samples generated in the MCS and analyzes these failure samples to assess the effects of various uncertainties on slope failure probability. The approach contains two major components: hypothesis tests for prioritizing effects of various uncertainties and Bayesian analysis for further quantifying their effects. A hypothesis test statistic ZH was formulated to evaluate the statistical difference of failure samples with their respective nominal (unconditional) samples. As the absolute value of ZH increases, the statistical difference between failure samples and their respective nominal samples becomes growingly significant, and the effect of the uncertain parameter on failure probability also becomes growingly significant. Therefore, the absolute value of ZH is used as an index to measure the effects of the uncertain parameters on failure probability and to prioritize their relative effects on failure probability. A Bayesian analysis approach was developed to further quantify effects of the uncertain parameters that have been identified from the hypothesis tests as influential parameters. Equations were derived for estimating conditional PDF [i.e., P(F|h)] of slope failure for a given value of uncertain parameter h. As P(F|h) is a variation of failure probability as a function of h, it can be considered as results of a sensitivity study of h on slope failure probability. In other words, a Bayesian analysis of the failure samples provides results that are equivalent to those from additional sensitivity studies. In addition, it has the advantage of avoiding additional computational times and efforts for repeated runs of MCS or Subset Simulations in the sensitivity studies. Furthermore, it was shown that comparison
1022
Y. Wang et al. / Computers and Geotechnics 37 (2010) 1015–1022
between the conditional probability P(h|F) and its unconditional one P(h) provides an indication of the effect of the uncertain parameter h on failure probability. This offers a means to check the Bayesian analysis results with those from hypothesis tests. The resolution of P(F) and P(h|F) is pivotal to obtain P(F|h), and it depends on the number of failure samples generated in MCS. An advanced Monte Carlo Simulation called Subset Simulation was employed to improve efficiency of generating failure samples in MCS and resolution at small failure probability levels. Subset Simulation algorithm and its implementation in a commonly-available spreadsheet environment were briefly described. As an illustration, the proposed probabilistic failure analysis approach was applied to study a design scenario of James Bay Dyke. The hypothesis tests show that the uncertainty of SuL has the most significant effect on the slope failure probability, while the uncertainty of Tcr contributes the least to the failure probability. The hypothesis test results are very consistent with results from independent sensitivity studies. Such agreement validates that the hypothesis test procedure proposed in this paper properly prioritizes effects of various uncertainties on failure probability. A Bayesian analysis was performed to quantify explicitly the effects of SuL and DTill, which have been identified from the Hypothesis tests as the two most influential uncertain parameters. It is shown that the slope failure probability change significantly as the values of SuL and DTill change. The Bayesian analysis results have also been validated against those from independent sensitivity studies. In addition, a cross-check between the hypothesis test results and the Bayesian analysis results shows that they agree well with each other. It is worthwhile to note that, although the proposed approach was developed together with a slope stability analysis problem, the approach is general and applicable to other types of geotechnical analyses and engineering problems. Acknowledgements The work described in this paper was supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China [Project No. 9041484 (CityU 110109)] and a strategic research grant from City University of Hong Kong (Project No. 7002455). The financial support is gratefully acknowledged. Reference [1] Tang WH, Yucemen MS, Ang HS. Probability based short-term design of slope. Can Geotech J 1976;13:201–15. [2] Christian JT, Ladd CC, Baecher GB. Reliability applied to slope stability analysis. J Geotech Eng 1994;120(12):2180–207. [3] Hassan AM, Wolff TF. Search algorithm for minimum reliability index of earth slopes. J Geotech Geoenviron Eng 1999;125(4):301–8. [4] Suchomel R, Masin D. Comparison of different probabilistic methods for predicting stability of a slope in spatially variable. Comput Geotech 2010;37(1–2):132–40. [5] Low BK, Tang WH. Reliability analysis of reinforced embankments on soft ground. Can Geotech J 1997;34(5):672–85. [6] Low BK, Gilbert RB, Wright SG. Slope reliability analysis using generalized method of slices. J Geotech Geoenviron Eng 1998;124(4):350–62. [7] Low BK. Practical probabilistic slope stability analysis. In: Proceeding of 12th Panamerican conference on soil mechanics and geotechnical engineering and 39th US rock mechanics symposium, M.I.T., Cambridge (MA): Verlag Gluckauf Gmbh Essen; 2003. p. 2777–84. [8] Low BK, Tang WH. Efficient spreadsheet algorithm for first-order reliability method. J Eng Mech 2007;133(2):1378–87. [9] Hong HP, Roh G. Reliability evaluation of earth slopes. J Geotech Geoenviron Eng 2008;134(12):1700–5. [10] El-Ramly H. Probabilistic analysis of landslide hazards and risks bridging theory and practice. PHD Thesis of University of Alberta; 2001.
[11] El-Ramly H, Morgenstern NR, Cruden DM. Probabilistic slope stability analysis for practice. Can Geotech J 2002;39(3):665–83. [12] Griffiths DV, Fenton GA. Probabilistic slope stability analysis by finite elements. J Geotech Geoenviron Eng 2004;130(5):507–18. [13] El-Ramly H, Morgenstern NR, Cruden DM. Probabilistic assessment of stability of a cut slope in residual soil. Geotechnique 2005;55(1):77–84. [14] Xu B, Low BK. Probabilistic stability analyses of embankments based on finite element method. J Geotech Geoenviron Eng 2006;132(11):1444–54. [15] Griffiths DV, Jinsong Huang, Fenton GA. Influence of spatial variability on slope reliability using 2-D random fields. J Geotech Geoenviron Eng 2009;135(10):1367–78. [16] El-Ramly H, Morgenstern NR, Cruden DM. Lodalen slide: a probabilistic assessment. Can Geotech J 2006;43(9):956–68. [17] Cho SE. Effects of spatial variability of soil properties on slope stability. Eng Geol 2007;92(3–4):97–109. [18] Grilli ST, Taylor ODS, Baxter CDP, Maretzki SA. Probabilistic approach for determining submarine landslide tsunami hazard along the upper east coast of the United States. Mar Geol 2009;264(1–2):74–97. [19] Cho SE. Probabilistic assessment of slope stability that considers the spatial variability of soil properties. J Geotech Geoenviron Eng 2010;136(7):975–84. [20] Wang Y, Cao Z, Au SK. Practical reliability analysis of slope stability by advanced Monte Carlo simulations in spreadsheet. Can Geotech J, accepted for publication. [21] Rubio E, Hall JW, Anderson MG. Uncertainty analysis in a slope hydrology and stability model using probabilistic and imprecise information. Comput Geotech 2004;31(7):529–36. [22] Schweiger HF, Peschl GM. Reliability analysis in geotechnics with the random set finite element method. Comput Geotech 2005;32(6):422–35. [23] Zhao HB. Slope reliability analysis using a support vector machine. Comput Geotech 2008;35(3):459–67. [24] Cho SE. Probabilistic stability analyses of slopes using the ANN-based response surface. Comput Geotech 2009;36(5):787–97. [25] Ching JY, Phoon KK, Hu YG. Efficient evaluation of reliability for slopes with circular slip surfaces using importance sampling. J Geotech Geoenviron Eng 2009;135(6):768–77. [26] Baecher GB, Christian JT. Reliability and statistics in geotechnical engineering. New York: John Wiley; 2003. [27] Au SK, Beck JL. Estimation of small failure probabilities in high dimensions by Subset Simulation. Probab Eng Mech 2001;16(4):263–77. [28] Au SK, Beck JL. Subset Simulation and its application to probabilistic seismic performance assessment. J Eng Mech 2003;129(8):1–17. [29] Luckman PG, Der Kiureghian A, Sitar N. Use of stochastic stability analysis for Bayesian back calculation of pore pressures acting in a cut at failure. In: Proceeding of 5th international conference on application of statistics and probability in soil and structure engineering. University of British Columbia, Vancouver, Canada; 1987. [30] Gilbert RB, Wright SG, Liedtke E. Uncertainty in back-analysis of slopes: Kettleman Hills case history. J Geotech Geoenviron Eng 1998;124(12): 1167–76. [31] Zhang J, Tang WH, Zhang LM. Efficient probabilistic back-analysis of slope stability model parameters. J Geotech Geoenviron Eng 2010;136(1):99–109. [32] Walpole RE, Myers RH, Myers SL. Probability and Statistics for Engineers and Scientists. 6th ed. Upper Saddle River (New Jersey): Prentice Hall; 1998. [33] Au SK. Reliability-based design sensitivity by efficient simulation. Comput Struct 2005;83(14):1048–61. [34] Ang HS, Tang WH. Probability concepts in engineering: emphasis on applications to civil and environmental engineering. 2nd ed. New York: John Wiley and Sons; 2007. [35] Au SK, Wang Y, Cao Z. Reliability analysis of slope stability by advanced simulation with spreadsheet. In: Proceedings of The Second International Symposium on geotechnical safety and risk. Gifu (Japan): CRC Press; 2009. p. 275–80. [36] Au SK, Cao Z, Wang Y. Implementing advanced Monte Carlo simulation under spreadsheet environment. Struct Saf 2010;32:281–92. [37] Metropolis N, Rosenbluth A, Rosenbluth M, Teller A. Equations of state calculations by fast computing machines. J Chem Phys 1953;21(6):1087–92. [38] Sturges H. The choice of a class-interval. J Am Stat Assoc 1926;21:65–6. [39] Duncan JM, Wright SG. Soil strength and slope stability. 1st ed. New York: John Wiley and Sons; 2005. [40] Ladd CC, Dascal O, Law KT, Lefebrve G, Lessard G, Mesri G, et al. In: Report of the subcommittee on embankment stability – annex II. Committee of specialists on sensitive clays on the NBR complex. Société d’Energie de la BaieJames, Montréal, Que; 1983. [41] Soulié M, Montes P, Silvestri V. Modelling spatial variability of soil parameters. Can Geotech J 1990;27:617–30. [42] Palisade Corporation. @Risk and the Decision Tools Suite Version 5.5.1. http:// www.palisade.com.2010. [43] Phoon KK, Kulhawy FH. Evaluation of geotechnical property variability. Can Geotech J 1999;36(4):625–39.