Assessment of the application of point estimate methods in the probabilistic stability analysis of slopes

Computers and Geotechnics 69 (2015) 540–550

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Assessment of the application of point estimate methods in the probabilistic stability analysis of slopes Morteza Ahmadabadi ⇑, Rainer Poisel 1 Technical University of Vienna, Karlsplatz 13, 1040 Wien, Austria

a r t i c l e

i n f o

Article history: Received 26 August 2014 Received in revised form 18 June 2015 Accepted 19 June 2015

Keywords: Probabilistic analysis Point estimate method Monte Carlo Simulation Reliability index Probability of failure

a b s t r a c t This study examines the applicability of point estimate methods (PEMs) in the probabilistic stability analysis of slopes. For illustrative purposes, two slopes indicating two different mechanisms of failure are selected and probabilistic analyses are performed employing three PEMs that have been widely cited in different literature reports. To evaluate the benefits and shortcomings of each PEM, various cases are studied in which random parameters may be dependent or independent and can have arbitrary distributions. Furthermore, the efficiency and accuracy of the considered PEMs are assessed by comparing the results with those obtained by Monte Carlo Simulation (MCS). The results indicate that the considered PEMs are sufficiently precise and, in case an intricate geotechnical structure must be probabilistically analysed, they can readily be incorporated with numerical methods. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Numerical analyses are becoming increasingly popular in geomechanics for the investigation of the stability of geotechnical structures. However, soil and rock are known to be among the most variable of all engineering materials, leading to uncertainties in the outputs of various analyses; hence, probabilistic analyses are performed to take these uncertainties into account. The first step in performing such analyses is the selection of an appropriate probabilistic tool. Probabilistic tools, such as the first-order reliability method (FORM) [13], the second-order reliability method (SORM) [2], the Monte-Carlo Simulation (MCS) [37] and the point estimate methods (PEMs) [12,14,15,23,28,35,36,42,43] are some alternatives. In combination with numerical methods, the MCS and PEMs may be more practical because they do not require the definition of the performance function in a closed form, as is required when FORM and SORM are used. The MCS is based on several deterministic analyses, and better precision in the model output is achieved if the selected sample is of sufficient size; thus, when the model computation is time-consuming, the Monte Carlo Simulation would not be practical for uncertainty analysis. Furthermore, the MCS requires the specification of the marginal distributions of the random variables. Under such circumstances, ⇑ Corresponding author. Tel.: +43 6509440594. E-mail addresses: [email protected] (M. Ahmadabadi), rainer.poisel@ tuwien.ac.at (R. Poisel). 1 Tel.: +43 1 58801 20301. http://dx.doi.org/10.1016/j.compgeo.2015.06.016 0266-352X/Ó 2015 Elsevier Ltd. All rights reserved.

PEM often constitutes a more practical alternative because it requires less computation time along with requiring only statistical moments as inputs. In this study, the applicability of PEMs in probabilistically analysing the stability of slopes is investigated for three PEMs that have been widely discussed in different literature reports. In the following report, first, the basics of the considered point estimate methods are presented (i.e., Section 2). For some PEMs, shortcomings arise when random parameters are correlated or non-normally distributed. The proposed solution to these drawbacks is addressed in Section 3. Finally, two illustrative examples comprising two slopes with two different modes of failure are represented in Section 4. 2. Point estimate methods (PEMs) In general, PEMs consist of replacing a continuous density distribution function i.e., f x ðxÞ by specifically defining discrete probabilities that are supposed to model the same low-order moments of that distribution function. The determination of these moments is performed by adding up the weighted discrete realizations. In Fig. 1, this relationship is shown in a simplified way by two realizations at xþ and x , represented by the weights wþ and w , respectively. This method is straightforward, easy to use and has become a staple of geotechnical reliability analyses. This method was first proposed by Rosenblueth [35,36]. Several researchers in the field of reliability, e.g., Lind [23], Zhou and Nowak [43], Harr [12], Hong [15] and Zhao and Ono [42], have commented on the primary point estimate method (i.e., Rosenblueth’s PEM) and implemented

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

well. Rosenblueth’s algorithm considers two realization points for each random variable. When random variables are symmetric or normally distributed, the realization points are located at l  r, in which l and r are the mean value and standard deviation of random variable, respectively. By permuting the two points in each dimension, 2n (n: number of random parameters) combinations of realization points exist, i.e., the model should be evaluated 2n times. Describing the details of the concept and formulation of this method is not the objective of this paper; for further studies, refer to Refs. [35,36]

x+,w+

PDF

x-,w-

2.2. Harr’s method

μ Fig. 1. Simplified graphical representation of the PEM.

some improvements. Regardless of the concept, efficiency and accuracy of the above-mentioned PEMs, they mainly differ in the number and location of realization points and in the method of treating asymmetric or correlated variables. In this study, three PEMs (i.e., Rosenblueth’s, Harr’s and Hong’s methods) are chosen to model the uncertainties that exist in the stability analyses of the two slopes. Table 1 summarizes a general comparison of these three PEMs. This comparison provides an initial idea of the applicability, advantages and impediments of each method. For example, irrespective of the accuracy of each method investigated in this study, it is clear that when the number of random variables, i.e., n, is large, Rosenblueth’s method may be computationally inefficient because it must evaluate the model m ¼ 2n times. In addition, deficiencies in Harr’s and Hong’s approaches arise when they address asymmetric or correlated variables, respectively. In this study, the Nataf transformation [21,24,27] (see Section 3) is proposed to address the drawback of these methods with respect to the dependent non-normal variables. As is clear from the table above, contrary to Reosenblueth’s and Harr’s PEMs, the Hong’s PEM has the capability to consider more than two realization points for each random variable. In this study, Hong’s approach is categorized into four schemes depending on the number (i.e. k ¼ 2; 3; and 5) that is selected for each random parameter.

2.1. Rosenblueth’s method The primary point estimate method proposed by Rosenblueth [35,36] is widely used in geotechnical practice for reliability calculations [6]. An advantage of this method is its simplicity in comparison with the other PEMs, although it has some disadvantages as

Table 1 Comparison of the three employed PEMs (n: number of variables). Author’s method

Rosenblueth Harr Hong 2n scheme 2n + 1 scheme 3n scheme 4n + 1 scheme

No. of realization points for each variable (k)

No. of model evaluations (m)

2 2

2n 2n

2

2n

3

2n þ 1

3

3n

5

4n þ 1

541

Capability for variables with

Asymmetry

Correlation

U  U

U U 

Using Harr’s algorithm (1989) [12] to solve problems involving n normally distributed and correlated random variables; a hyperpffiffiffi sphere with a radius n centered at the origin of standardized parameters space is defined as representing the correlation matrix. The realization points at which the model should be evaluated are located at intersections of the hypersphere and eigenvectors of the correlation matrix. Because each of the n eigenvectors has two intersections with the hypersphere, 2n realization points are produced, i.e., the model should be evaluated 2n times. Although Harr’s PEM is computationally more efficient compared with Rosenblueth’s PEM, Harr’s PEM has the disadvantage that it cannot model non-normally distributed variables. In addition, when the number of variables n is large, the realization points are located farther from the mean, thereby causing the inaccuracy or inapplicability of Harr’s PEM [6]. Further explanations of Harr’s PEM can be found in Refs. [1,4,6,12,14,28]. 2.3. Hong’s method When using Hong’s method [15], more than two realization points can be selected. In addition, Hong’s method takes into account higher-order moments of random variables; as a result, it is possible to calculate higher-order moments of the model output. In the case that, for each variable, k estimation points are considered, the model should be evaluated a total of n  k times. In the present study, this method is categorized into four schemes based on the number of realization points considered. For example, the 4n þ 1 scheme corresponds to the case that for each variable, five realization points are selected, such that one point is located at mean value and two points are located at each side of the mean. Note that Hong’s PEM also has some shortcomings. An important drawback of Hong’s PEM is that it is unable to model correlated random parameters. Furthermore, sometimes some schemes of this method are unusable because realization points are located beyond the marginal distribution or they are imaginary numbers. This alternative is not discussed in detail, and further details can be found in Refs. [3,6,15]. 3. Nataf transformation As mentioned in Section 2, Harr’s PEM is incapable to consider asymmetric variables; conversely, Hong’s schemes are only applicable for uncorrelated variables. To circumvent such limitations, the Nataf transformation [27] is implemented when dealing with correlated non-normal variables. In the following, the Nataf transformation is briefly described. Let the correlated random variables vector X be

X ¼ ðx1 ; x2 ; . . . ; xn Þ

ð1Þ

Assuming the marginal CDF F xi ðxi Þ; i ¼ 1; . . . ; n; of each random variable xi ; i ¼ 1; . . . n; of each random variable X i and the correlation matrix q ¼ ðqij Þnn are known, and applying the

542

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

Fig. 2. Geometry of slope I.

Table 2 Material properties for slope I. Material

f x ðxÞ ¼

Parameters

MH Weathered A M Weathered A Sed. Marine. Fresh Andesite Mat 7 Mat 8

c ðkN=m3 Þ

rt (kPa)

c (kPa)

/ (°)

E (MPa)

m (–)

21 21 19 21 19 21

8 7 4 100 0 100

8 7 4 100 0 100

34 24 17 30 37 30

50 50 50 50 50 50

0.4 0.4 0.4 0.4 0.4 0.4

Generally, this distribution model is referred to as the Nataf distribution. To solve q0 ¼ ðq0ij Þnn in Eq. (4), any two random variables

qij ¼

Uðyi Þ ¼ F xi ðxi Þ

ð2Þ

yi ¼ U1 ½F xi ðxi Þ

in which UðÞ and U1 ðÞ are the standard normal CDF and the inverse standard normal, respectively, with joint probability density function Un ðy; q0 Þ having zero means, unit standard deviation and a correlation matrix q0 ¼ ðq0ij Þnn that is given by

1

  1 T 1 y q0 y 2

Un ðy; q0 Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ð2pÞn jq0 j

ð4Þ

ðxi ; xj Þ are considered. The linear correlation between these random variables is calculated by

isoprobabilistic transformation, the random vector X can be transformed to the standardized normal random vector Y ¼ ðy1 ; y2 ; . . . ; yn Þin y space as follows:

(

f X1 ðx1 Þf X2 ðx2 Þ; . . . ; f Xn ðxn Þ /n ðy; q0 Þ /ðy1 Þ/ðy2 Þ; . . . ; /ðyn Þ

ð3Þ

Then, according to Nataf transformation theory and the differentiation rules of implicit functions, the approximate joint probability density function f x ðÞ in x space is expressed as

Z

þ1

Z

1

þ1

1

F 1 i ½Uðyi Þ  lxi

!

rxi

F 1 j ½Uðyj Þ  lxj

rxj

! /2 ðyi ; yj ; q0ij Þdyi dyj ð5Þ

where lxi and lxj are the means of correlated random variables xi and xj , respectively; rxi and rxj are the standard deviation of correlated random variables xi and xj , respectively. Given the marginal distributions and the correlation for a pair of random variables, the equivalent correlation coefficient q0ij can be determined by solving the nonlinear equation, as shown in Eq. (5). Because a practical solution of this equation is not easily obtained, Liu and Der Kiureghian [24] developed a set of semi-empirical formulae. However, such empirical formulae may not always be valid for random variables following arbitrary distributions. In this study, the solution of this equation, i.e., Eq. (5), is performed using a trial and error procedure by means of a MATLAB code that is presented in Appendix A. By obtaining the equivalent correlation matrix q0 , which can be decomposed as

q0 ¼ C0 CT0

(where C0 is the lower triangular matrix obtained from Cholesky decomposition), the dependent standard normal vector Y can be obtained from the independent standard normal random vector U as follows:

Table 3 Statistical information of random parameters of slope I. Material

Random parameters c (kPa)

MH Weathered A M Weathered A Sed. Marine.

ð6Þ

/ (°)

Y ¼ C0 U

l

r

l

r

8 7 4

2 2 1

34 24 17

3 3 2

ð7Þ

The original non-normal random variables are given by

xi ¼ F 1 i ½Uðyi Þ;

i ¼ 1; 2; . . . ; n

ð8Þ

543

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

Fig. 3. A typical result of the model evaluation of slope I (i.e., safety factor) using the shear strength reduction method.

Table 4 Results of the probabilistic stability analysis of slope I using PEMs. Employed numerical code

FLAC

PEM

No. model evaluations m

Rosenblueth Harr & Hong’s 2n scheme Hong 2n + 1 & 3n schemes 4n + 1 scheme

Factor of safety

Reliability index

Probability of failure

l

r

ðbÞ

P f ð%Þ

64 12

1.036 1.038

0.0893 0.0889

0.4009 0.4252

34.42 33.53

13 & 18 25

1.029 1.021

0.0873 0.0842

0.3276 0.2458

37.16 40.29

Table 5 Results of probabilistic stability analysis of slope I performed by Rocscience [31]. Employed numerical code

Probabilistic tool

Factor of safety

Probability of failure

l

r

P f ð%Þ

64

1.05

0.092

30.40

10,000

1.07

0.091

22.0

No. model evaluations m

different mechanisms of failure incorporating above-mentioned PEMs with different numerical methods.

the

4.1. Example 1  slope I with circular failure mode

4. Illustrative examples

This example is taken from Rocscience [31] and aims to reveal how uncertainties that exist in the shear strength of slope materials can be taken into consideration by means of the aforementioned PEMs. In addition, the examined slope indicates the case that a cut is made in weak materials, such as highly weathered or closely fractured rock and rock fills. In such materials, a strongly defined structural pattern (wedge or toppling) no longer exists, and the slide surface is free to find the line of least resistance through the slope. Observations of slope failures in these materials suggest that this slide surface generally takes the form of a circle [41]. In this case, the failure of rock slope is not controlled by geological features, such as bedding planes and joints that divide the rock into discontinues mass; thus, the slope is supposed to be as continuum and it is analysed by FLAC [16] using the strength reduction technique [39].

In this section, the applicability and usefulness of the aforementioned PEMs in probabilistic stability analysis of slopes is considered. The examples are used to investigate probabilistically two

4.1.1. Slope geometry The examined slope depicted in Fig. 2 has a height of H ¼ 60 m and a general face inclination of wf  58 . It is assumed that water

Finite element (phase 2) Analytical (Slide)

Rosenblueth’s PEM MCS

Table 6 Number of simulations required for calculating pf ¼ 35% with 10% error using MCS in example 1. Standard deviation (d)

Confidence (%)

No. simulations

1.28 1.64 1.96

80 85 95

305 500 714

544

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

table is very low, meaning that its effects could be disregarded in the stability calculations. Note that for the sake of consistency and to compare the results of analyses with those presented by [31], the slope geometry and the properties of slope materials are represented without any manipulation. 4.1.2. Material properties In Fig. 2, slope consists of six materials. The properties of the slope materials are represented in Table 2, in which c; rt ; c; /; E and m denote unit weight, tensile strength, cohesion, internal friction angle, elasticity module and Poisson’s ratio, respectively. As mentioned above, the shear strength of the slope materials is assumed to be uncertain. In addition, the behavior of the materials is assumed to be modeled by the Mohr–Coulomb failure criterion, meaning that cohesion ðcÞ and friction angle ð/Þ are treated as random variables, and 6  2 ¼ 12 random variables should be considered in probabilistic stability calculations; to make the probabilistic analyses more efficient, the friction angle and cohesion of three materials termed Mat 8, Mat 7, and Fresh Andesite are assumed to be deterministic because they have higher values in comparison to other materials, and they do not determine the location of critical failure surface (i.e., with minimum safety factor). Table 3 contains the statistical information of the random parameters that are assumed to be uncorrelated and normally distributed. 4.1.3. Probabilistic stability analyses As mentioned earlier, probabilistic stability analyses of slope are performed by incorporating the above-mentioned point estimate methods with FLAC. Fig. 3 displays a sample of the model evaluations in terms of the safety factor; for clarity, the velocity vectors and maximum shear strains are also displayed. Note that in the numerical model, the boundary conditions correspond to fixed movement in the horizontal direction at the right and the left end points of the model (i.e., x ¼ 0 m & x ¼ 200 m), and the fixed movements in the horizontal and vertical directions at the lowest part of the model (i.e. y ¼ 40 m). 4.1.4. Results Table 4 summarizes the results of probabilistic stability analyses of slope using the above-mentioned PEMs. These results include the mean ðlÞ and standard deviation ðrÞ of the safety factor, reliability index ðbÞ and probability of failure ðPf Þ. Note that for calculating Pf and b, it is necessary to know (or estimate) the

distribution of the model output (i.e., safety factor). In this example, it is supposed that safety factor has a normal distribution; hence, Pf and b are computed using Eq. (9).

b¼

n¼

ð10Þ

100 Number of simulations Ns=500

90

Number of simulations Ns=1500

80 70

50

Frequency

Frequency

 2 d 1p p a

Table 6 represents the number of simulations required to calculate the probability of failure Pf ¼ 35% with 10% error with different confidences. According to the table below, to calculate the probability of failure with 10% error and 85% confidence (i.e. standard deviation d = 1.64) with MCS, sample size should be Ns = 500 (see row 2 of the Table 6). It is apparent that MCS is much more inefficient in comparison with PEMs when it is coupled with numerical modeling codes (e.g., FLAC). To better understand the effect of sample size in the results of MCS, two sets of probabilistic stability analyses are conducted for samples with sizes N s1 ¼ 500 and N s2 ¼ 1500. The results of these two probabilistic analyses are compared in Fig. 4 in terms of the histogram of safety factor. In Fig. 5, the histograms of input random parameters for a samples with size of Ns2 = 1500 are presented.

60

40 30

60 50 40 30

20

20

10 0

ð9Þ

where UðÞ denotes the cumulative density function (CDF) of the standard normal distribution. Note that when random variables are uncorrelated and normally distributed, Hong’s 2n scheme and Harr’s algorithm are identical. Similarly, Hong’s 2n þ 1 scheme and Hong’s 3n scheme are also the same. A comparison of the results of different point estimate methods indicates that there is good agreement among these methods; however the results of probabilistic stability analyses performed by Rocscience [31] using finite element [32] and analytical methods [34] are provided in Table 5 for further comparison. To highlight the accuracy of the PEMs, the stability of slope is also probabilistically analysed using Monte Carlo Simulation (MCS). As mentioned before, the accuracy of MCS is directly correlated with the number of simulations; nevertheless, the number of simulations was chosen to be as small as possible because model evaluation with a PC that is equipped with a quad core processor takes approximately 25 min, i.e., if a large sample is selected, then probabilistic analyses by means of MCS are infeasible. To estimate the number of simulations required, Eq. (10) is used [7]. This formula indicates that, for calculating the probability of failure p with a% error, n simulations are required.

80 70

l1 ; Pf ¼ 1  UðbÞ; r

10 0.85

0.9

0.95

1

1.05

1.1

Safety factor

1.15

1.2

1.25

0

0.9

0.95

1

1.05

1.1

Safety factor

Fig. 4. Histogram of the safety factors for slope I with MCS, left graph N s1 ¼ 500, right graph N s2 ¼ 1500.

1.15

545

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

100 90

120 Number of samples=1500

Number of samples=1500

100

80

80

60

Frequency

Frequency

70

50 40

60

40

30 20

20

10 0

0

4000 5000 6000 7000 8000 9000 10000110001200013000

30

32

Cohesion1 (Pa)

Number of samples=1500

90

Frequency

Frequency

40

6000

7000

Number of samples=1500

70

60 50 40 30

60 50 40 30

20

20

10

10 13

14

15

16

17

18

19

20

0

21

1000

2000

3000

4000

5000

Cohesion2 (Pa)

Friction angle 2 90

120

Number of samples=1500

Number of samples=1500

80

100

70 60

Frequency

80

Frequency

38

80

70

0

36

100

90 80

34

Friction angle1

60

50 40 30

40

20 20 10 0

2000 3000 4000 5000 6000 7000 8000 9000 1000011000

Cohesion 3 (Pa)

0

19

20

21

22

23

24

25

26

27

28

29

Friction angle 3

Fig. 5. Histograms of the random parameters of slope I, with normal distribution fits.

In addition, a quantitative comparison is also demonstrated in Table 7. An increase in the number of simulations is found to lead to a considerable change in P f , although the mean value of the safety factor changes slightly. From the data in Table 7, and by substituting n ¼ 1500; a ¼ 10% and p ¼ 21% in Eq. (10), we obtain

d ¼ 1:99, which means that with a sample size of n ¼ 1500, the probability of failure, pf , would be 21%, with more than 95% confidence (see Table 6). Note that these sample sizes are large for this type of numerical calculation, and this makes the probabilistic stability analyses expensive.

546

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

Table 7 Results of the probabilistic stability analysis of slope I with MCS. Employed numerical code

No. simulations

Factor of safety

l

r

Probability of failure P f ð%Þ

FLAC

N s1 ¼ 500 N s2 ¼ 1500

1.035 1.040

0.0815 0.0524

32 21.3

and wp is the dip of the sliding plane. b denotes the location of tension crack behind the slope crest, and B is the location of tension crack behind the slope toe. The following assumptions are made in the plane failure analysis:  Both sliding surface and tension cracks strike parallel to the slope.  The tension crack is vertical and is filled with water to a depth zw .  Water enters the sliding surface along the base of the tension crack and seeps along the sliding surface, escaping at atmospheric pressure where the sliding surface daylights in the slope face.  The shear strength s of the sliding surface is defined by cohesion c and friction angle /, which are related by the Mohr– Coulomb failure criterion. 4.2.2. Material properties Again, for the sake of consistency and to compare the results of the probabilistic analyses with those presented by Wyllie & Mah [41], the input parameters and variables are considered without any manipulation. The deterministic and statistic parameters that are utilized in the probabilistic stability analysis of slope II are listed in Tables 8 and 9. The statistical parameters of slope 2 are listed in Table 9. In this table, wp , /, c and zw symbolise the dip of the sliding plane, the friction angle, the cohesion of the joint and the depth of water in the tension crack, respectively.

Fig. 6. Geometry of slope II with plane failure.

Table 8 Deterministic parameters of slope II. Parameter

Unit

Value

H

m

cw crock

kN m3 kN m3

30.50 10

wf

deg. m m m

b B z

27 58 12.50 31.56 19.01 m

4.2. Example 2  rock slope II with plane failure This example is taken from chapter 6 of Wyllie & Mah [41] and aims to probabilistically investigate the stability of a rock slope in the context of plane failure. Although plane failure is a comparatively rare sight in rock slopes, it would not be appropriate to ignore the two-dimensional case because there are many valuable lessons to be learned from a consideration of the mechanics of this simple failure mode. Plane failure is particularly useful for demonstrating the sensitivity of the slope to changes in shear strength and ground water conditions – changes that are less obvious when dealing with the more complex mechanics of three-dimensional slope failure [41]. 4.2.1. Slope geometry The slope geometry and ground water condition considered in this analysis are defined in Fig. 6. In this figure, z is the depth of the tension crack, H is the slope height, wf is the slope face angle

4.2.3. Probabilistic stability analyses Before evaluating the model, the realization points should be computed based on the desired PEMs. In this example, the random variables are assumed to be uncorrelated, but two random variables, i.e., c and zzw have non-normal distributions; consequently, as mentioned in Sections 2 and 3, Harr’s PEM should be manipulated by means of the Nataf transformation, i.e., the non-normal variables are first transformed to the normal space and then Harr’s PEM is used to obtain the realization points. In this research, this procedure that is a combination of Harr’s PEM and a normal transformation is called the modified Harr’s PEM. Table 10 presents a comparison among the realization points that are calculated using the modified Harr’s and the primary Harr’s PEMs. Patently, in the case that the non-normal random parameter has a substantial effect on the model output, the modification of Harr’s method is crucial. 4.2.4. Results Although, in a 2D model of rock slope with plane failure, there is an analytical solution for the safety factor (see chapter 6 of [41]), probabilistic stability analyses of rock slopes are accomplished by coupling the aforesaid point estimate methods with UDEC [17]. The distinct element model for the rock slope is shown in Fig. 7. In the figure, the results of the model evaluation in terms of safety factor along with the vectors of displacement of the sliding block are also observed.

Table 9 Statistical parameters of slope II. Variable

Distribution

Mean ðlÞ

Standard deviation ðrÞ

Lowe bound

Upper bound

wp (deg.) / (deg.) c (kPa) zw z (%)

Normal Normal Triangular Triangular

20 20 80 50

2.4 2.7 – –

– – 40 0

– – 130 100

547

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550 Table 10 Comparison of the realization points that are computed by Harr’s and modified Harr’s PEMs in the probabilistic stability analyses of slope II. Realization points

1 2 3 4 5 6 7 8

Harr’s PEM

Modified Harr’s PEM

wp (deg.)

/ (deg.)

c (kPa)

zw z

24.8 15.2 20 20 20 20 20 20

20 20 25.4 14.6 20 20 20 20

83.33 83.33 83.33 83.33 120.15 46.51 83.33 83.33

50 50 50 50 50 50 90.8 9.20

(%)

wp (deg.)

/ (deg.)

c (kPa)

zw z

24.8 15.2 20 20 20 20 20 20

20 20 25.4 14.6 20 20 20 20

82.56 82.56 82.56 82.56 119.88 49.05 82.56 82.56

50 50 50 50 50 50 89.33 10.66

(%)

Fig. 7. Distinct element model of slope II with UDEC.

Table 11 Results of the probabilistic analyses of slope II using PEMs. Probabilistic tool

No. model evaluations

Safety factor

Model output distribution Normal

Rosenblueth Harr Hong 2n scheme Hong 2n + 1 scheme Hong 3n scheme Hong 4n + 1 scheme

16 8 8 9 12 17

Log-normal

Mean (l)

Std. ðrÞ

Safety index ðbÞ

P f (%)

P f (%)

1.3581 1.3604 1.3583 1.3590 1.3535 1.3364

0.2561 0.2463 0.2513 0.2524 0.2462 0.2531

1.3983 1.4633 1.4258 1.4223 1.4358 1.3291

8.10 7.17 7.70 7.75 7.56 9.20

6.13 5.22 5.73 5.78 5.62 7.34

548

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

Table 12 Results of the probabilistic stability analyses of slope II with MCS. Numerical code

No. simulations

Safety factor

RocPlane [18,38]

10,000

1.3559

0.2558

6.26

UDEC

1000 6000

1.3797 1.3435

0.2563 0.2260

4.80 5.60

Mean (l)

P f (%) Std. ðrÞ

As in the first example, the mean and standard deviation of the safety factor and the safety index ðbÞ and the probability of failure ðP f Þ are calculated using three PEMs, the results of which are included in Table 11. It is remarkable that although the two variables are non-normal, all of the methods produce different realization points, and they should be considered separately in the

probabilistic stability calculations. Again, note that to interpret the results of probabilistic analyses by means of PEMs in terms of probability of failure, an appropriate distribution of the model output (safety factor) should be considered. In this example, because two random variables have non-normal distributions, it might be thought that they affect the distribution of the safety factor; as a result, a lognormal distribution of the safety factor is also included in the analyses. Note that probability of failure for lognormal distribution is calculated using Eq. (11).

l2Fos

mu ¼ ln qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; l2Fos þ r2Fos

700 600

No. of simulations Ns=1000

120 100

Frequency

Frequency

lFos

140 No. of simulations Ns=6000

500 400 300

80 60

200

40

100

20

0

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 ffi lFos þ r2Fos sigma ¼ ln 2

0.6

0.8

1

1.2

1.4

1.6

1.8

0

2

0.8

1

1.2

1.4

1.6

1.8

2

2.2

Safety factor

Safety factor

Fig. 8. Histograms of the safety factor for slope II with MCS: left graph N s1 ¼ 6000, right graph N s2 ¼ 1000.

Table 13 qc;/ Variation of F ¼ q 0 for different quantities of qc;/ in the probabilistic stability analysis of slope II. c;/

qc;/

qc;/

qc;/

qc;/0

F¼ q 0 c;/

qc;/

qc;/0

F¼ q 0 c;/

0.10 0.20 0.30 0.40 0.50

0.1005 0.2009 0.3013 0.4018 0.5022

1.0050 1.0045 1.0043 1.0045 1.0044

0.60 0.70 0.80 0.90 1.0

0.6026 0.7030 0.8035 0.9039 1.0043

1.0043 1.0043 1.0044 1.0043 1.0043

Table 14 Comparison of the realization points that are produced by the original and the modified Hong’s PEM (3n scheme) in the probabilistic stability analysis of slope II. Realization points

1 2 3 4 5 6 7 8 9 10 11 12

Original Hong’s PEM ðqc:/ ¼ 0Þ

Modified Hong’s PEM (qc:/ ¼ 0:50)

wp (deg.)

/ (deg.)

c (kPa)

zw (%)

wp (deg.)

/ (deg.)

c (kPa)

zw (%)

20 24.16 15.84 20 20 20 20 20 20 20 20 20

20 20 20 20 24.67 15.32 20 20 20 20 20 20

83.33 83.33 83.33 83.33 83.33 83.33 113.46 56.81 79.04 83.33 83.33 83.33

50 50 50 50 50 50 50 50 50 50 81.60 18.39

20 24.16 15.84 20 20 20 20 20 20 20 20 20

20 20 20 20 24.67 15.32 20 20 20 20 20 20

83.23 83.23 83.23 83.23 67.03 101.24 59.90 79.04 83.23 83.23 83.23 83.23

50 50 50 50 50 50 50 50 50 50 81.60 18.39

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

1.6

/ are negatively correlated. A series of correlation coefficients, qc;/ ¼ 0 1:0; is used to model the common shear test results, in which the cohesion generally decreases as the friction angle increases, and vice versa [25,26]. The Nataf transformation is used to consider the correlation where the reliability analyses are performed by means of Hong’s method. Note that, although cohesion, c, has a triangular distribution, because its distribution is approximately symmetric about the mean value, the coefficient factor

1.5

F¼

2

Rosenblueth's PEM Harr's PEM Hong 2n scheme Hong 2n+1 scheme Hong 3n scheme Hong 4n+1 scheme Monte Carlo Simulation.

Safety Index (β)

1.9 1.8 1.7

1.4 1.3 0

-0.1

-0.2

-0.3

-0.4

-0.5

-0.6

-0.7

-0.8

-0.9

-1

Correlation coefficient (ρc,φ) Fig. 9. Variation of the reliability index versus the correlation coefficient in the probabilistic stability analyses of slope II.

8.2

Probability of failure (Pf ) %

549

Rosenblueth's PEM Harr's PEM Hong 2n scheme Hong 2n+1 scheme Hong 3n scheme Hong 4n+1 scheme Monte Carlo Simulation

7.2 6.2 5.2 4.2 3.2

2.2

-0.05

-0.15

-0.25

-0.35

-0.45

-0.55

-0.65

-0.75

-0.85

-0.95

q0ij qij ,

as it is described by [21,24], would be equal to 1.0. This

matter is readily proved by solving Eq. (5), the results of which are given in Table 13. According to this table, when the random variables are transformed to the normal space, the variation in correlation coefficient ðq0ij  qij Þ is less than 0.50%. The table below, Table 14, indicates a comparison between the original Hong’s PEM and the modified Hong’s PEM in terms of the realization points that are used for the model evaluations. The importance of this manipulation by means of the NATAF transformation is better appreciated by examining Figs. 9 and 10, in which the variation of the safety index, b, and the probability of failure, Pf , versus the correlation coefficient are shown, respectively. It is evident from the figures above that the results obtained using the different methods are not much different, demonstrating the accuracy of proposed solution by means of the Nataf transformation. In these figures, the Monte Carlo Simulations are performed using [33] a sample size of n ¼ 100; 000. Evidently, by considering the correlation coefficient between cohesion c and the internal friction angle of rock joints /, the safety index, b, and the probability of failure, P f , change considerably. In addition, in this example, disregarding the correlation coefficient between cohesion c and internal friction angle of rock joints / clearly leads to a conservative geotechnical design because the probability of failure decreases as the correlation coefficient increases.

Correlation coefficient (ρc,φ)

5. Summary and conclusions Fig. 10. Variation of the probability of failure versus the correlation coefficient in the probabilistic stability analyses of slope II.

Pf

Lognormal

¼ logncdf ð1; mu; sigmaÞ

ð11Þ

where logncdf ðÞ denotes the cumulative distribution function corresponding to a lognormal distribution, which can be calculated using MATLAB code. To estimate the accuracy of the results of the PEMs, the results of three probabilistic analyses based on MCS are given in Table 12. These analyses are conducted with different sample sizes and with different numerical methods. The first results are those of Wyllie & Mah [41], in which the stability of the discussed slope is analysed using an analytical method [33]. The other two Monte Carlo Simulations are performed with UDEC for N s1 ¼ 1000 and N s2 ¼ 6000 simulations. Fig. 8 exhibits a comparison between the histograms of the safety factor based on the results of Monte Carlo Simulations with UDEC. A comparison of the results of Monte Carlo Simulations and the results calculated by the PEMs proves the accuracy of the PEMs; furthermore, it illustrates that assuming a log-normal distribution for the safety factor is more accurate. 4.2.5. Correlation of the random variables As mentioned in Section 2, Hong’s approach was originally established for uncorrelated random variables. This section represents how to treat correlated random parameters when using Hong’s method; moreover, the effect of correlation coefficient variation on reliability index, b, and probability of failure, pf , is studied. Among the four uncertain parameters (see Table 9), it is assumed that cohesion c and internal friction angle of rock joints

Application of point estimate methods in the reliability analysis of slopes was examined. Three PEMs were investigated using two examples, which include two slopes with two modes of failure. The following conclusions can be made from the results of this study. 1. Although Rosenblueth’s PEM is applicable for both skewed and asymmetric variables, it may be infeasible due to the computational effort required when the number of random parameters is large. 2. In problems dealing with normally distributed random variables and when the number of random variables is not large, Harr’s algorithm is not only very efficient but also provides sufficiently accurate results. For problems involving random parameters with a non-normal distribution, the modified Harr’s PEM is proposed as an alternative. In the second example, the asymmetry of the variables changes the probability of failure from 7.0% to 7.17%. However, undoubtedly, when non-normal variables play an important role in the performance function, considering the asymmetry of non-normal variables may drastically change the results. 3. In contrast to Rosenblueth’s and Harr’s methods, Hong’s PEM is capable of theoretically considering more than two realization points for each random parameter. In the case that all variables are uncorrelated and normally distributed, Harr’s PEM and Hong’s 2n scheme, as well as Hong’s 3n and 2n þ 1 schemes, are the same. 4. When Hong’s PEM is utilized for problems involving correlated variables, the Nataf transformation is proposed to resolve the drawbacks of Hong’s PEM with respect to the correlated

550

M. Ahmadabadi, R. Poisel / Computers and Geotechnics 69 (2015) 540–550

variables. In the second example, the probability of failure is observed to vary from Pf ¼ 8:0% to P f ¼ 2:8% when the correlation coefficient qc;£ varies from 0 to 1. 5. The more accurate is assumed the distribution of performance function, the more precise are the results obtained by the PEMs; as a result, if an analytical or approximate solution of the examined problem is possible, running the MCS is recommended to achieve a better estimation of the distribution of the performance function.

Acknowledgements The authors thank the anonymous reviewers for their constructive comments, which improved this paper. The authors acknowledge the funding of this work by the Austrian Society for Geomechanics. The first author is grateful to Professors A. Rohatsch and A. Preh for their ongoing supports, advice and suggestions. Appendix A. ‘‘MATLAB code for solving the Eq. (5)’’

clc; % Input information of random variables mu1 = 20; std1 = 2.70; a2 = 40; mu2 = 80; b2 = 130; delta = 0.0001; R = 0.90; pd1 = makedist(‘normal’,’mu’,mu1,’sigma’,std1); std1 = std(pd1); mu1 = mean(pd1); pd2 = makedist(‘triangular’, ’a’, a2, ’b’, mu2, ’c’, b2); std2 = std(pd2); mu2 = mean(pd2); % This array includes the outputs count = 0; a1 = 7; a2 = 7;% Integral Limits q = 1; while q > 1e5 count = count + 1; R0 = R + delta ⁄ count fun = @(z1, z2) (((icdf(pd1,normcdf(z1,0,1)))m u1)./std1).⁄ ... (((icdf(pd2, normcdf(z2, 0,1)))mu2)./std2).⁄ ... 1./(2 ⁄ pi. ⁄ sqrt(1R0^2)).⁄ ... exp(1 ⁄ (z1.^2 + z2.^2-2. ⁄ R0. ⁄ z1. ⁄ z2)./... (2. ⁄ (1R0.^2))); q = integral2(fun, a1, a1, a2, a2)R end F = R0/R

References [1] Baecher GB, Christian JT. Reliability and statistics in geotechnical engineering. Wiley; 2003. [2] Breitung K, Hohenbichler M. Asymptotic approximations for multivariate integrals with an application to multinormal probabilities. J Multivariate Anal 1989;30(1):80–97. [3] Caramia P, Carpinelli G, Varilone P. Point estimate schemes for probabilistic three-phase load flow. Electr Power Syst Res 2010;80(2):168–75. [4] Che-Hao C, Yeou-Koung T, Jinn-Chuang Y. Evaluation of probability point estimate methods. Appl Math Model 1995;19(2):95–105. [6] Christian JT, Baecher GB. The point-estimate method with large numbers of variables. Int J Numer Anal Meth Geomech 2002;26(15):1515–29.

[7] Gibson W. Probabilistic methods for slope analysis and design. Austral Geomech 2011;46(3):29–40. [12] Harr ME. Probabilistic estimates for multivariate analyses. Appl Math Model 1989;13(5):313–8. [13] Hasofer AM, Lind NC. An exact and invariant first-order reliability format. J Eng Mech, ASCE 1974;100(1):111–21. [14] He J, Sällfors G. An optimal point estimate method for uncertainty studies. Appl Math Model 1994;18(9):494–9. [15] Hong HP. An efficient point estimate method for probabilistic analysis. Reliab Eng Syst Saf 1998;59(3):261–7. [16] ITASCA consulting group Inc. Fast lagrangian analysis of continua (FLAC) ver.7.0: two-dimensional explicit finite difference program for engineering mechanics computation; 2015. [17] ITASCA consulting group Inc. Universal distinct element code (UDEC) ver.5.0: advanced, two-dimensional distinct element modeling for geotechnical analysis of rock, soil and structural support; 2015. [18] Jiang S, Li D. Reliability analysis of unsaturated slope with spatially correlated soil properties. In: Vulnerability, uncertainty, and risk. American Society of Civil Engineers; 2014. p. 2429–38. [21] Lemaire M. Structural reliability. John Wiley & Sons; 2009. [23] Lind NC. Modelling of uncertainty in discrete dynamical systems. Appl Math Model 1983;7(3):146–52. [24] Liu P-L, Der Kiureghian A. Multivariate distribution models with prescribed marginals and covariances. Probab Eng Mech 1986;1(2):105–12. [25] Low BK. Reliability analysis of rock slopes involving correlated nonnormals. Int J Rock Mech Min Sci 2007;44(6):922–35. [26] Low BK. Efficient probabilistic algorithm illustrated for a rock slope. Rock Mech Rock Eng 2008;41(5):715–34. [27] Nataf A. Determination des distributions de probabilites dont les marges sont donnes. Comptes Renud de l’Academie des Sci 1962;225(5):42–3. [28] Panchalingam G, Harr ME. Modelling of many correlated and skewed random variables. Appl Math Model 1994;18(11):635–40. [31] Rocscience Inc. Probabilistic slope stability analysis. Available: . [32] Rocscience Inc. Phase 2 ver.8.0: finite element analysis for excavations and slopes; 2015. [33] Rocscience Inc. RocPlane ver.3.0: planar sliding stability analysis for rock slopes; 2015. [34] Rocscience Inc. Slide ver.6.0: 2D limit equilibrium slope stability analysis; 2015. [35] Rosenblueth E. Point estimates for probability moments. Proc Natl Acad Sci USA 1975;72(10):3812–4. [36] Rosenblueth E. Two-point estimates in probabilities. Appl Math Model 1981;5(5):329–35. [37] Rubinstein RY, Kroese DP. Simulation and the Monte Carlo method. 2nd ed. Wiley; 2008. [38] Shamekhi E, Tannant DD. Probabilistic assessment of rock slope stability using response surfaces determined from finite element models of geometric realizations. Comput Geotech 2015;69:70–81. [39] Matsui T, San K-C. Finite element slope stability analysis by shear strength reduction technique. Soils Found 1992;32(1):59–70. [41] Wyllie DC, Mah C. Rock slope engineering. 4th ed. Taylor & Francis; 2004. [42] Zhao Y-G, Ono T. New point estimates for probability moments. J Eng Mech 2000;126(4):433–6. [43] Zhou J, Nowak AS. Integration formulas to evaluate functions of random variables. Struct Saf 1988;5(4):267–84.

Further reading [5] Christian JT, Baeche GB. Point estimate method as numerical quadrature. Geotech Geoenviron Eng 1999;125(9):779–86. [8] Griffiths DV, Huang J, Fenton GA. Probabilistic infinite slope analysis. Comput Geotech 2011;38(4):577–84. [9] Griffiths VD, Fenton GA. Probabilistic methods in geotechnical engineering. Springer; 2007. [10] Haldar A, Mahadevan S. Probability, reliability, and statistical methods in engineering design. John Wiley; 2000. [11] Haldar A, Mahadevan S. Reliability assessment using stochastic finite element analysis. Wiley; 2000. [19] Kim JM, Sitar N. Reliability approach to slope stability analysis with spatially correlated soilproperties. Soils Found 2013;53(1):1–10. [20] Le TMH. Reliability of heterogeneous slopes with cross-correlated shear strength parameters. In: Geotechnical safety and risk IV. CRC Press; 2013. p. 161–7. [22] Li D, Chen Y, Lu W, Zhou C. Stochastic response surface method for reliability analysis of rock slopes involving correlated non-normal variables. Comput Geotech 2011;38(1):58–68. [29] Phoon KK. Reliability-based design in geotechnical engineering: computations and applications. Taylor & Francis; 2008. [30] Reale C, Xue J, Pan Z, Gavin K. Deterministic and probabilistic multi-modal analysis of slope stability. Comput Geotech 2015;66:172–9. [40] Wang Y, Tung Y-K. Improved probabilistic point estimation schemes for uncertainty analysis. Appl Math Model 2009;33(2):1042–57.