1. INTRODUCTION

Stochastic simulations of fractured rock masses can

provide valuable information for the engineering

design of rock slopes, particularly when the natural

geologic discontinuities may form potential slope

failure modes. An essential component of such

engineering simulations is being able to compute

the probability of sliding for a given potential

failure mass once the geometry of that failure mass

has been identified through spatial rock-fracture

simulations. In cases where several thousand (or

more) simulations of possible failure geometries are

needed to provide a realistic representation of the

rock slope, computational efficiency is essential for

the repetitious protocol used to calculate the

probability of sliding.

An example of applying this type of geotechnical

approach to rock slope design was presented by

Miller and others [1] in a paper focused on the

design of catch benches for open pit mines and

quarries. This issue is important to NIOSH

(National Institute for Occupational Safety and

Health) as part of its research mission to improve

safety and health in the mining industry. Between

1995 and 2003 there were 42 reported fatalities due

to slope failure accidents at surface mines in the

United States, at least one and as many as eleven

each year. Additional accident statistics collected

by the Mine Safety and Health Administration

(MSHA) have shown that loose material from slope

and bench failures can pose significant safety

hazards to miners. To address these concerns the

NIOSH Spokane Research Laboratory has been

developing and testing rock slope stability software

over the past few years to provide advanced

technical tools for analyzing bench stability. A key

element of this software is a module used to

compute the probability of sliding for a given viable

failure geometry that has been simulated for the

bench under study.

Several different methods can be used to compute

the probability of sliding in a slope stability

analysis, including Monte Carlo simulation [2],

 

APPLICATIONS OF THE POINT ESTIMATION METHOD FOR

STOCHASTIC ROCK SLOPE ENGINEERING

Stanley M. Miller

Department of Civil Engineering, University of Idaho, Moscow, ID, USA

Jeffrey K. Whyatt & Edward L. McHugh

CDC-NIOSH Spokane Research Laboratory, Spokane, WA, USA

 

 

ABSTRACT: The point estimation method can be applied to the safety factor (SF) equation for any specified rock slope failure 

mode (such as plane shear, step path, or wedge) to obtain reliable estimates of the mean and standard deviation of the SF 

probability distribution. A gamma probability density function is recommended for modeling this probability distribution, because 

it allows only for positive values and is flexible enough to provide symmetrical shapes and right-skewed, exponential-type shapes 

for the SF distribution. The mean and standard deviation define this distribution, which then can be integrated numerically from 0 

to 1 to obtain the probability of sliding, PS (portion of the SF distribution where SF < 1.0). The overall probability of failure, PF, 

for the potential slope failure mass is the joint probability that the rock discontinuities are long enough to allow kinematic failure 

(PL) and that sliding occurs along the rock discontinuities (PS); that is, PF = PSPL. This method for estimating the probability of 

sliding is extremely efficient computationally, and thus, expedites slope stability simulation routines used by NIOSH software to 

stochastically describe rock slope behavior and assist the engineer in designing catch benches for large rock slopes. Enhanced 

bench design translates into increased operational efficiency and safer working conditions in open pit mines and quarries. 
Taylor series expansion [2], Fourier analysis [3],

and statistical point estimation [4]. Initial versions

of the NIOSH bench stability software relied on the

Fourier method [1] due to its numerical efficiency

and its capability to provide a discretized, general

output pdf (probability density function) for the

factor of safety, rather than relying on a specific

model for the pdf (e.g., a normal or lognormal pdf).

However, our recent experience with probabilistic

studies of rock slope stability has indicated that the

safety factor pdf tends to behave like a slightly

right-skewed gamma pdf or a left-truncated normal

pdf (truncated at zero, because the safety factor

realistically cannot take on negative values). The

right skew apparently is caused by a combination of

the positive-only gamma shape of the input pdf for

the shear strength (along the sliding plane) and the

exponential-type shape of the input pdf for the

fracture waviness. In this context, the waviness is

measured on a scale of about 1-2 meters and is

defined as the average dip of the fracture minus its

minimum dip, as presented by Call and others [5].

Thus, assuming that the output pdf for the safety

factor takes on the form of a gamma pdf, then the

point estimation method clearly has considerable

computational advantage even over the Fourier

method (which relies heavily on mathematical

manipulations of discretized pdf’s [3]).

2. POINT ESTIMATION METHOD

When a random variable of interest can be

expressed in an equation as the result of a

mathematical operation of other random variables,

then the point estimation method developed by

Rosenblueth [6] provides a direct computational

procedure to obtain moment estimates for that

random variable. In particular, these statistical

moments are the mean (i.e., the first moment about

the origin) and the variance (i.e., the second

moment about the mean). Geotechnical engineering

applications of this method have been around for

several decades, and recent publications [2, 4] have

clearly presented such work. The particular shape

of any pdf used for any input random variable is not

critical to the analysis, because the pdf is

represented by the mean and two hypothetical point

masses located at plus and minus one standard

deviation (s) from the mean (μ).

Consequently, required inputs for a probabilistic

rock slope stability analysis are: 1) a defined

performance function (i.e., safety factor equation),

2) estimated value for each input attribute if it is

assumed to have negligible variability, and 3)

estimated mean μ and standard deviation s of each

input attribute treated as a random variable. Typical

attributes for a rock slope failure mass with a

defined geometry would include shear strength,

rock mass unit weight, and fracture waviness.

Calculation steps are presented below for the point

estimation method using two random variables X1

and X2 in a performance function to obtain the mean

and variance of F, the factor of safety.

1. Calculate the output value of F using the

performance function evaluated with the values of

mean-plus-one-s.d. for each of the two variables.

F++ = fn[(μ1 + s1), (μ2 + s2)] (1a)

Repeat for other combinations, as follows:

F- - = fn[(μ1 - s1), (μ2 - s2)] (1b)

F+ - = fn[(μ1 + s1), (μ2 - s2)] (1c)

F- + = fn[(μ1 - s1), (μ2 + s2)] (1d)

2. Calculate the point-mass “weights” [6].

P++ = P- - = (1/4)(1 + ρ12) (2a)

P+ - = P- + = (1/4)(1 - ρ12) (2b)

where ρ12 = correlation coefficient between

input variables X1 and X2.

3. Calculate the expectation (mean, μF) of F [6].

E(F) = P++ F++ + P- - F- - + P+ - F+ - + P- + F- + (3)

4. Calculate the variance (sF

2

) of F.

Var(F) = E(F

2

) – [E(F)]

2

(4)

where E(F

2

) is calculated using Eq. (3) with F

2

terms substituted for the F terms.

5. The standard deviation (sF) of F then is calculated by taking the square root of sF

2

.

For three input variables X1, X2, X3, there are eight

calculations in Step 1, and the point-mass weights

in Step 2 are given by [6]:

P+++ = P- - - = (1/8)(1 + ρ12 + ρ23 + ρ31) (5a)

P+ - - = P- ++ = (1/8)(1 - ρ12 + ρ23 - ρ31) (5b)

P++ - = P- - + = (1/8)(1 + ρ12 - ρ23 - ρ31) (5c)

P+ - + = P- + - = (1/8)(1 - ρ12 - ρ23 + ρ31) (5d)
Eq. (3) for E(F) is extended from a summation of

four terms to a summation of eight terms for this

case. The sF value subsequently can be calculated

using Eq. (4), after first using eight F

2

terms in the

extended Eq. (3) to calculate E(F

2

).

3. PLANE SHEAR ANALYSIS

A probabilistic procedure based on the point

estimation method can be applied to the twodimensional slope stability analysis used for the

plane shear failure mode (Fig.1). For a plane shear

mass with a defined geometry (i.e., slope angle δ,

dip of failure plane α, height of the failure mass H),

the SF equation is assumed to contain two random

variables, the shear strength and the waviness of the

geologic discontinuity. Other input terms, such as

the rock-mass density, which is used to calculate the

weight of the potential failure mass, are treated as

constants. The plane shear SF equation is:

sin( )

tan( )

α

τ σ

W

L L r

n F

+

=

( ) tan( )

sin( ) sin( )

r W

L

W

L n

⋅ ⋅

















= +

α

τ

α

σ

(6)

where: τ = shear strength, r = wavinessσn =

effective normal stress, α = average dip of sliding

surface, L = length of sliding surface, and W =

weight of slide mass. The point estimation method

can be applied to this expression if we have reliable

estimates of the mean and standard deviation of

each random variable.

L

W

H

α

geologic discontinuity

δ

Fig. 1. Example of plane shear failure mode.

3.1. Mean and S.D. of Shear Strength

The mean of τ is obtained from the linear or power

curve expression that relates the effective normal

stress σn to the shear strength [7], as shown by the

following options:

Linear: τ = c + σntanφ (7)

where: c = cohesion, φ = friction angle.

Power: τ = c + aσn

b

(8)

where: c = y-intercept of power curve, a,b = shape

parameters for power curve. This model reverts to a

linear model when b equals 1.

JRC Power: τ = σntan[JRC·log(JCS/σn) + φb] (9)

where: JRC= joint roughness coefficient, JCS= joint

wall compressive strength, φb = base friction angle.

For a simple 2d rock-slope bench geometry (shown

in Fig. 1), the effective normal stress is:

U L

W

n

−

⋅

= 1

cos(α)

σ (10)

where: W = total weight of the failure mass, U =

ground-water pore pressure on the sliding surface.

The standard deviation of τ is estimated using the

standard error (se) obtained from a least-squares

regression analysis of data obtained by laboratory

direct shear tests of natural fractures [7]. For a

linear regression model, we can estimate the

standard deviation of τ at a given value of effective

normal stress σn by the expression [8, p. 371]:

2

2

1 ( )

τ

τ

τ

σ

SSQ n m

m

n

s s

n

e

− ⋅

−

= + (11)

where: se = standard error of the regression; n =

number of data pairs used in the regression; mτ =

sample mean of the n τ values; and

= 

n

i

SSQ

2

τ (sum of the squared τ values)

For a power-curve regression model, Eq. (11) also

should give reasonable results if the se value from

that power-curve model is used. Example results

from laboratory direct-shear tests of natural rock

joints for actual rock-slope design projects are

summarized in Table 1. Stress units for the reported

shear strengths are given in tonnes per square meter

(tsm). As seen by the se values for the three rock

types, the regression errors are small for the latite
porphyry (a “tight” regression fit) and quite high for

the granodiorite (a “noisy” regression fit). These

examples in Table 1 were selected to show typical

se values for shear strength regression models.

Estimated standard deviations of τ reported at the

bottom of Table 1 were computed using Eq. (11).

The maximum normal stress acting on potential

sliding surfaces in a 25-m high bench is about 20

tsm; this justifies the 2, 10, and 18 tsm σn values in

the table. The representative sτ results shown in

Table 1 provide general guidance for reasonable

and acceptable values to be used in stochastic rock

slope studies, especially when direct-shear test

results are not available. If the engineer anticipates

little variability in discontinuity shear strength for

the critical failure mode in a mine bench, then a

reasonable value for sτ is 0.1 tsm (0.98 kPa). If a lot

of variability in shear strength is expected, then a

value up to 0.6 tsm (5.88 kPa) should be assumed.

Our experience has indicated that in most cases

reasonable “default” values are 0.3 to 0.4 tsm for sτ .

Table 1. Shear strength estimates for power-curve regressions.

Regression

Terms

Latite

Porphyry

Quartz

Monzonite

Granodiorite

a 1.3802 1.5412 0.8099

b 0.8377 0.7381 0.8468

c 0.0 0.0 0.0

se (tsm) 0.1823 0.3953 0.5718

Est. τ @ 2,

10, 18 tsm

2.47

9.50 15.54

2.57

8.43 13.01

1.46

5.69 9.36

SD of τ @ 2,

10, 18 tsm

0.118

0.105 0.202

0.260

0.213 0.382

0.383

0.294 0.482

3.2. Mean and S.D. of Waviness

The mean and standard deviation of tan(r) are

computed by applying a random variable transformation from waviness r to tan(r), then

numerically integrating the probability density

function of the new transformed variable.

Assuming the waviness can be modeled as an

exponentially distributed random variable, the pdf

of the new random variable V (i.e., v = tan(r)) is

given by [3]:







	










<

≥

+

=

−

0 0

0

(1 )

( )

2

[ arctan( )/ ]

if v

if v

v

e

f v

r

v

V

r

μ

μ

(12)

where: μr = mean waviness (in degrees).

Thus, for any given mean value of the fracture

waviness, the pdf of V is defined and can be

integrated numerically to obtain the mean and

variance of V (in other words, the mean and

variance of tan(r)). These procedures are outlined

below.

Mean of V:



∞

= = ⋅ ⋅

0

E[V] v f (v) dv μv V(13)

Variance of V:

Var[V] = E[V

2

] – μv

2

(14)

where:



∞

= ⋅ ⋅

0

2 2 E[V ] v f (v) dv V

Standard deviation of V:

s Var[V ] v = (15)

Numerical integrations needed to evaluate Eqs. (13)

and (14) have been completed for waviness values

from 1

o

to 12

o

, and the results are listed in Table 2.

Note that values for 0

o waviness are not needed,

because the expression in Eq. (6) reduces to only

one random variable for r = 0

o

.

Therefore, once a specified (simulated) mean

waviness value is known for a given plane shear

failure, the corresponding pdf of V (which has an

exponential-type shape) is defined by Eq. (12).

Waviness values should only be reported and used

as integer, whole numbers. Then, the appropriate

values for the mean and standard deviation of V can

be obtained from Table 2.

Table 2. Calculated means and standard deviations of random

variable V; note that v = tan(r).

Waviness (deg.) Mean of V S.D. of V

1 0.017464 0.014337

2 0.034992 0.035167

3 0.052654 0.053264

4 0.070523 0.072055

5 0.088691 0.091989
6 0.107277 0.114034

7 0.126444 0.140049

8 0.146394 0.172409

9 0.167338 0.212872

10 0.189444 0.261799

11 0.212803 0.318279

12 0.237413 0.380713

3.3. Example Calculations

To illustrate the point estimation procedures for

computing the probability of sliding for a 2-d plane

shear failure mode, a simple bench geometry is

considered with the following terms (see Fig. 1):

Bench face angle: δ = 65

o

Height of failure mass at face: H = 4.0 m

Dip of sliding surface: α = 32

o

Waviness of sliding surface: r = 3

o

Unit weight of rock mass: γ = 2.6 tcm (25.5 kN/m

3

)

Calculated length of sliding surface: L = 7.55 m

Calculated weight of slide mass: W = 23.588 tonne

Calculated effective normal stress: σn = 2.650 tsm

Shear strength parameters:

a = 0.6512 b = 0.988 c = 0.0

Calculated mean τ: μτ = 0.6512σn

0.988

= 1.706 tsm

Standard deviation of τ: sτ = 0.3 tsm

For 3

o waviness, μv = 0.052654 and sv = 0.053264

The first constant in Eq. (6) is computed as:

sin( )

1 W α

L

C = = 0.6039

The second constant in Eq. (6) is computed as:

sin( )

2

α

σ

W

L

C

n

= = 1.6003

Point estimate calculations from Eq. (1):

F++ = 0.6039(1.706 + 0.3) + 1.6003(0.052654 + 0.053264)

= 1.3807

F- - = 0.6039(1.706 - 0.3) + 1.6003(0.052654 - 0.053264)

= 0.8479

F+ - = 0.6039(1.706 + 0.3) + 1.6003(0.052654 - 0.053264)

= 1.2102

F- + = 0.6039(1.706 - 0.3) + 1.6003(0.052654 + 0.053264)

= 1.0184

Because the shear strength τ often is based on

small-scale laboratory test specimens (~15 cm) and

the waviness is measured over 1-2 m in the field, it

is reasonable to assume these two attributes are not

dependent on each other. Thus, the correlation

coefficient between the two is zero. The statistical

point-mass weights then are calculated from Eq. (2):

P++ = P- - = P+ - = P- + = (1/4)(1 + 0) = 0.25

Calculate the mean safety factor from Eq. (3):

μF = 0.25(1.3807 + 0.8479 + 1.2102 + 1.0184)

μF = 1.114

Then, calculate the variance of the safety factor

from Eq. (4):

sF

2

= 0.25(1.3807

2

+ 0.8479

2

+ 1.2102

2

+ 1.0184

2

) – 1.114

2

= 0.0401

Thus, the standard deviation of the safety factor is:

sF = 0.2002

A gamma pdf [8, p.161-163] for the safety factor

then can be defined and numerically integrated from

0 to 1 to estimate the probability of sliding for this

plane shear analysis. The u and z shape parameters

of the gamma pdf are computed directly from the

mean and variance.

u

1.114

2

0.0401

:= z

0.0401

1.114

:=

Ps

1

z

u ( ⋅ Γ (u)) 0

1

e x

x

z













−

x

u−1

⋅












:= d

Probability of sliding: Ps = 0.298

The safety factor gamma pdf is shown in Figure 2.

It is worthy to note that even though the mean

safety factor is 1.114 (indicating stability), this

plane shear actually has a 0.30 probability of

sliding, which results directly from the variability of

input values.
Fig. 2. Safety factor pdf for plane shear example based on

point estimation method; mean = 1.114, s.d.= 0.200, PS = 0.30.

4. WEDGE ANALYSIS

The probability of sliding also can be estimated for

a three-dimensional, tetrahedral wedge that may

form in the rock slope as a result of two intersecting

planar discontinuities (Fig. 3). The safety factor

equation for this type of wedge can be expressed as:

f

R R R R R

f

L L L L L

D

A A r

D

A A r

F

τ σ tan( ) τ σ tan( )

+

+

=

R

f

R

L

f

L

D

A

D

A

F τ τ

















+

















=

tan( ) tan( ) R R

f

R

L L

f

L

r

D

A

r

D

A

σ σ

















+

















+ (16)

where: AL = sliding area of left side of wedge, AR =

sliding area of right side of wedge, Df = driving

force of the wedge, τL = shear strength along left

side of wedge, τR = shear strength along right side

of wedge, σL = effective normal stress on left side

of wedge, σR = effective normal stress on right side

of wedge, rL = waviness of left side of wedge, and

rR = waviness of right side of wedge.

The shear strengths on either side of the wedge are

treated as random variables, and all other inputs are

considered as constants. Thus, the first two terms in

Eq. (16) are used in the point estimation method,

with the sum of the last two terms being added later

to the estimated mean of the safety factor. That is,

the last two terms in Eq. (16) effectively comprise a

waviness constant that increases the mean safety

factor but does not affect its variance.

δ

H

Fig. 3. Example of 3-d wedge failure mode.

4.1 Mean and S.D. of Shear Strength

Shear strengths for the left and right sides of the

wedge are handled in a similar fashion to that

presented in Section 3.1. In some cases (e.g., for a

symmetrical wedge), the mean shear strength will

be the same for both the left and right sides of the

wedge, but in the general case the two shear

strengths do not have to be the same. This is due to

the weight of the wedge being shared between the

two sliding surfaces. Thus, the effective normal

stress applied to each side of the wedge can vary

significantly, depending on the steepness of the

wedge plunge, the dihedral angle of the wedge, and

other geometric factors.

Also, application of the same type of shear strength

model to both sides of the wedge is not required for

this stability analysis. For example, the engineer

can specify a power-curve regression model for the

left side and a linear model for the right side. Once

the effective normal stress is computed for each

sliding plane of the wedge, the corresponding mean

and standard deviation of τ for each side are

obtained independently.

4.2. Example Calculations

To illustrate the procedures for a 3-d wedge, the

following example analysis is provided for a simple

bench geometry:

Bench face angle: δ = 65

o

Bench face dip direction: Ddir = 160

o

Height of failure mass at face: H = 4.2 m

Dip direction and dip of left plane: 106

o

, 51

o

Dip direction and dip of right plane: 219

o

, 52

o

Effective waviness of left plane: 3.3

o

Effective waviness of right plane: 4.8

o

0 1 2 3 4 5

0

0.5

1

1.5

2

2.5

gamma pdf

f(x)

x
(Note: Effective waviness is the apparent angle of waviness

observed in the vertical plane in the direction of sliding rather

than in the direction of down-dip, so it is always less than the

true waviness angle.)

Shear strength parameters:

Left: a = 0.5543 b = 0.997 c = 0.0

Right: a = 0.5317 b = 0.988 c = 0.0

Unit weight of rock mass: γ = 2.6 tcm (25.5 kN/m

3

)

Calculated length of wedge intersection: L = 7.23 m

Calculated volume of slide mass: V = 14.37 m

3

Calculated weight of slide mass: W = 37.36 tonne

Calculated bearing and plunge of the wedge intersection: 162

o

, 36

o

Calculated driving force: Df = 21.73 tonne

Calculated area of left side: AL = 12.92 m

2

Calculated area of right side: AR = 12.17 m

2

Calculated effective normal stresses:

σL = 1.565 tsm σR = 1.652 tsm

Calculated mean τ:

Left: μτ L = 0.5543σ L

0.997

= 0.866 tsm

Right: μτ R = 0.5317σ R

0.988

= 0.910 tsm

Standard deviation of τ:

sτ L = 0.3 tsm sτ R = 0.5 tsm

The first two constants in Eq. (16) are computed:

AL /Df = 0.5945 AR /Df = 0.5600

The waviness constant is computed:

0.5945(1.565)tan(3.3

o

) + 0.5600(1.652)tan(4.8

o

)

= 0.132

Point estimate calculations:

F++ = 0.5945(0.866 + 0.3) + 0.5600(0.910 + 0.5)

= 1.4827

F- - = 0.5945(0.866 - 0.3) + 0.5600(0.910 - 0.5)

= 0.5660

F+ - = 0.5945(0.866 + 0.3) + 0.5600(0.910 - 0.5)

= 0.9227

F- + = 0.5945(0.866 - 0.3) + 0.5600(0.910 + 0.5)

= 1.1261

Calculate the mean safety factor from Eq. (3):

μF = 0.25(1.4827 + 0.5660 + 0.9227 + 1.1261)

= 1.024

Then, calculate the variance of the safety factor

from Eq. (4):

sF

2

= 0.25(1.4827

2

+ 0.5660

2

+ 0.9227

2

+ 1.1261

2

) – 1.024

2

= 0.111

Thus, the standard deviation of the safety factor is:

sF = 0.333

A gamma pdf for the safety factor distribution that

results from these stochastic procedures then can be

defined and numerically integrated from 0 to 1 to

estimate the probability of sliding for this 3-d

wedge analysis. The u and z shape parameters of

the gamma pdf are computed directly from the

mean and variance.

u

1.024

2

0.111

:= z

0.111

1.024

:=

Ps

1

z

u ( ⋅ Γ (u)) 0

1

e x

x

z













−

x

u−1

⋅












:= d

Probability of sliding: Ps = 0.514

The safety factor gamma pdf is shown in Figure 4.

It is worthy to note that even though the mean

safety factor is 1.024 (indicating stability), this 3-d

wedge actually has a 0.51 probability of sliding,

which results directly from the variability of input

values.

0 1 2 3 4 5

0

0.5

1

1.5

2

2.5

gamma pdf

f(x)

x

Fig. 4. Safety factor pdf for 3-d wedge example based on

point estimation method; mean = 1.024, s.d.= 0.333, PS = 0.51.
5. STEP-PATH ANALYSIS

Calculating the probability of sliding for step-path

failure modes (see Fig. 5) requires the application of

three variables in the safety factor equation. The

safety factor equation used here is based on the

work of Jaeger [9] and subsequent adaptations of

the step-path analysis to rock slope engineering by

others [10, 11]:

sin( )

tan( )

α

τ σ

W

L L r TL

n F

+ +

=

( ) ( ) T W

L

W

L

W

L r

n

⋅ ⋅ + ⋅

















= +

sin( ) sin( ) sin( )

tan( )

α α

τ

α

σ

(17)

where: L = effective sliding length of the potential

step-path failure; r = waviness of the master joint

set; α = average dip of the master joint set; W =

weight of the potential failure mass; and T is the

effective tensile strength (i.e., resisting force) due to

intact rock bridges, given by:

T = Fir(ht)(To)/L (18)

where: Fir = fraction of intact rock contained in

tensile rock bridges along the step path; ht =

perpendicular height of the potential step-path

failure mass; and To = estimated tensile strength of

the intact rock. For a given step-path geometry, all

the terms in Eq. (18) are constants except for To, the

tensile strength, which typically is estimated by a

set of Brazilian disk tension tests. If this set

consists of at least four specimens, then a mean and

standard deviation can be computed for To, which

then are used to describe this input variable in the

safety factor computations. Thus, the third term in

Eq. (17) can be rewritten as:

( ) o

T W

I

⋅

sin(α)

where I = Fir(ht) .

The following intermediate calculations are useful

for helping to define the final constants to be used

in the point estimation method.

Effective length of sliding (side AC of triangle ABC

in Fig. 5; obtained by applying the law of sines):

sin( )sin( )

sin( )

β α α

α β

−

−

=

c

H c

L (19)

where: αc = average dip of cross joint set; β =

overall average angle of the step path.

β

α

master joint

cross joint

intact rock gap

h

H

t

δ

L

A

B

C

Fig. 5. Example of step-path failure mode.

Weight of the potential failure mass:

2sin( )sin( )

sin( )

2

δ β

γ δ − β

=

H

W (20)

where: γ = unit weight of the rock mass.

Perpendicular height of the potential failure mass:

sin( )

sin( )

β α

β

= −

H

ht(21)

5.1. Practical Considerations

Simulated step paths are generated by combining

realistic fractures from both the master joint set and

the cross joint set using the method presented in

[10]. A path is “continuous” if it continues uninterrupted from the daylight point on the slope face all

the way to the top of the bench. In this case, the Fir

value is zero, meaning that the third term in Eq. (17)

can be ignored. However, in many cases a

continuous step path does not form (e.g., a cross

joint length is not long enough to span the spacing

between adjacent master joints), and a new path is

initiated from the gap, leaving an intact rock

“bridge” that must fail in tension in order for the

entire failure mass to slide. The fraction of intact

rock is estimated as the summed length of all such

rock bridges divided by the value of ht

. Our

experience has shown that when this value exceeds

about 0.08 the probability of sliding is nil. Thus, in

computer routines the Ps is assigned a very small

value (say, 0.000001) when Fir > 0.08, and the point

estimation calculations are not needed.

If the master joint set is very planar with a waviness

angle of zero, then the second term in Eq. (17) can
be ignored, meaning that only two input variables

are analyzed. However, for the general case of

step-path analysis, all three variables are used in the

point estimation method as described earlier in

Section 2. These three variables are:

τ, shear strength along the master joint set;

tan(r), tangent of the waviness angle for the

master joint set;

To, tensile strength of the intact rock

Their respective constants (multipliers) are:

L /(Wsinα); (22a)

Lσn /(Wsinα); (22b)

Ι /(Wsinα). (22c)

It is worthy to note that engineering units of force

(actually, force per unit width for the 2-d analysis)

are preserved in the numerators for all three of the

terms in the safety factor equation (tonne/m). Such

units cancel the same units in the denominators (i.e.,

weight per unit width, or tonne/m) to provide a

unitless value for the safety factor. This was also

the case for the 2-d plane shear analysis presented

previously in Section 3.

5.2. Example Calculations

To demonstrate the point estimation procedures for

analyzing a 2-d step-path failure mode, a simple

bench geometry is considered with the following

terms (refer to Fig. 5):

Bench face angle: δ = 65

o

Height of failure mass at face: H = 6.1 m

Average dip of master joint set: α = 33

o

Waviness of master joint set: r = 4

o

Average dip of cross joint set: αc = 79

o

Unit weight of rock mass: γ = 2.6 tcm (25.5 kN/m

3

)

Mean tensile strength of intact rock: 270 tsm

S.D. of tensile strength of intact rock: 28 tsm

Step-path simulation provides a failure path with

one intact rock gap:

Overall step-path angle: β = 51

o

Total span of intact rock gap(s): I = 0.035 m

Calculated effective length of step path: L = 5.12 m

Calculated weight of slide mass: W = 16.62 tonne

Calculated effective normal stress: σn = 2.722 tsm

Shear strength parameters:

a = 0.5549 b = 0.988 c = 0.0

Calculated mean τ: μτ = 0.5549σn

0.988

= 1.492 tsm

Standard deviation of τ: sτ = 0.4 tsm

For 4

o waviness, μv = 0.070523 and sv = 0.072055

The first constant from Eq. (22a) is:

sin( )

1 W α

L

C = = 0.5656

The second constant from Eq. (22b) is:

sin( )

2

α

σ

W

L

C

n

= = 1.5396

The third constant from Eq. (22c) is:

sin( )

2 W α

I

C = = 0.0038

Point estimate calculations adapted from Eq. (1):

F+++= 0.5656(1.492 + 0.4) + 1.5396(0.070523 + 0.072055)

+ 0.0038(270 + 28) = 2.4220

F- - - = 0.5656(1.492 - 0.4) + 1.5396(0.070523 - 0.072055)

+ 0.0038(270 - 28) = 1.5349

F++ -= 0.5656(1.492 + 0.4) + 1.5396(0.070523 + 0.072055)

+ 0.0038(270 - 28) = 2.2092

F+ - -= 0.5656(1.492 + 0.4) + 1.5396(0.070523 - 0.072055)

+ 0.0038(270 - 28) = 1.9873

F- ++= 0.5656(1.492 - 0.4) + 1.5396(0.070523 + 0.072055)

+ 0.0038(270 + 28) = 1.9694

F- - += 0.5656(1.492 - 0.4) + 1.5396(0.070523 - 0.072055)

+ 0.0038(270 + 28) = 1.7477

F+ - += 0.5656(1.492 + 0.4) + 1.5396(0.070523 - 0.072055)

+ 0.0038(270 + 28) = 2.2002

F - + -= 0.5656(1.492 - 0.4) + 1.5396(0.070523 + 0.072055)

+ 0.0038(270 - 28) = 1.7567

Assuming independence between each pair of the

three variables, the correlation coefficients are zero.

The statistical point-mass weights then are given by

Eq. (5): (1/8)(1 + 0) = 0.125 .

Calculate the mean safety factor from an expanded

form of Eq. (3):

μF = 0.125(2.422 + 1.5349 + 2.2092 + 1.9873

+ 1.9694 + 1.7477 + 2.2002 + 1.7567)

= 1.978
Then, calculate the variance of the safety factor by

applying Eq. (4):

sF

2

= 0.125(2.422

2

+ 1.5349

2

+ 2.2092

2

+ 1.9873

2

+ 1.9694

2

+ 1.7477

2

+ 2.2002

2

+ 1.7567

2

) – 1.978

2

= 0.0765

Thus, the standard deviation of the safety factor is:

sF = 0.2766

A gamma pdf for the safety factor distribution that

results from these stochastic procedures then can be

defined and numerically integrated from 0 to 1 to

estimate the probability of sliding for this step-path

analysis. The u and z shape parameters of the

gamma pdf are computed directly from the mean

and variance.

u

1.978

2

0.0765

:= z

0.0765

1.978

:=

Ps

1

z

u ( ⋅ Γ (u)) 0

1

e x

x

z













−

x

u−1

⋅












:= d

Probability of sliding: Ps = 0.000074

The safety factor gamma pdf is shown in Figure 6.

The calculated probability of sliding is small here,

primarily due to the relatively high mean value of

the safety factor.

0 1 2 3 4 5

0

0.5

1

1.5

2

2.5

gamma pdf

f(x)

x

Fig. 6. Safety factor pdf for step-path example based on point

estimation method; mean = 1.978, s.d.= 0.277, PS = 0.00007.

For this same step-path geometry with no intact

rock gaps (i.e., I = 0.0), we can neglect the

contribution from tensile strength and use only the

first two terms in Eq. (17). This situation leads to a

calculated mean safety factor of 0.952 and s.d. of

0.2537, which yield a value of 0.60 for the

probability of sliding. This simple example clearly

shows the significant increase in stability provided

by small intact rock bridges along a step-path

failure surface.

The driving force for this 2-d step path example is

Wsin(α) = 9.05 tonne. The expected contribution to

stability from the small intact rock bridge is found

by: (0.035 m)(1 m)(270 tsm) = 9.45 tonne, which is

the tensile force that can be carried by the intact

rock bridge without rupturing. Thus, we see why

the mean safety factor can increase from 0.952 to

1.978 when the effect of the intact rock bridge is

included in the stability analysis. If observed field

conditions do not indicate solid, intact zones of rock

substance (e.g., rock masses that are highly

fractured or that contain weathered or altered zones)

that can form such rock bridges, then the engineer

should consider reducing the input value for the

mean tensile strength. This is especially important

when the tensile strength has been estimated from

small, intact specimens used in laboratory testing.

6. DISCUSSION

The point estimation method provides a computationally efficient method to calculate the mean and

standard deviation of the safety factor for common

rock slope failure modes. The shape of the safety

factor pdf must be specified by the engineer in order

to estimate the probability of sliding by numerically

integrating the area under this pdf from zero to one.

This pdf shape will closely resemble that of the pdf

chosen for the shear strength. It is prudent to use a

pdf model for shear strength that only takes on

positive values, such as the gamma pdf. We

recommend this choice, which rationally leads to an

assumption of a gamma pdf for the safety factor as

well. For this model, the shape of the safety factor

distribution appears fairly bell-shaped for small

values of the standard deviation, but becomes more

right-skewed as the variance increases (Fig. 7).

For the plane-shear and wedge failure modes, the

probability that the failure surface is long enough to

extend all the way to the top of the bench must be

considered in the final calculation to obtain the
probability of failure (note: this probability of sufficient length is not needed for the step-path mode,

which presumes a continuous failure path including

any small intact rock bridges present). The length

required for a continuous plane shear is shown by L

in Fig. 1. The length required for a wedge failure is

measured along the intersection line of the wedge,

and both planes must have trace lengths at least as

long as this intersection for the wedge to be viable.

A reasonable and convenient model for the length

distribution in a fracture set is the exponential pdf

[1, 5], which is a one-parameter pdf defined by its

mean. The cumulative distribution function (cdf;

denoted by F(x)) of an exponential variable can be

used directly to obtain the probability that length in

a fracture set will take on a value at least as long as

that required for the failure mode to be viable:

PL = 1 – F(x) = 1 - (1 - e

-x/m

) = e

-x/m

(23)

where: m = mean length of the fracture set. Thus,

for the plane shear mode, the probability of failure

is: PF = PSPL , and for the wedge mode it is given

by: PF = PS(PL1PL2), which reflects the joint probability of sufficient length for the two planar

discontinuities that form the wedge. There are

accepted methods for estimating the mean trace

length of a fracture set based on field data obtained

by scan-line or window mapping [12]. Such

procedures provide the mean length to be used in

Eq. (23) to obtain the needed PL values.

These methods for estimating the probability of

failure have been programmed into several software

0 1 2 3 4 5

0

0.8

1.6

2.4

3.2

4

A. mean=1.1, sd=0.12, Ps=0.21

f(x)

x

0 1 2 3 4 5

0

0.4

0.8

1.2

1.6

2

B. mean=1.1, sd=0.3, Ps=0.40

f(x)

x

0 1 2 3 4 5

0

0.4

0.8

1.2

1.6

2

C. mean=1.1, sd=0.6, Ps=0.51

f(x)

x

Fig. 7. Comparison of safety factor gamma pdf’s for different

standard deviation values; the mean SF in all figures is 1.1.

components used in NIOSH bench simulation

software [1], which can be used to help engineers

stochastically describe rock slope behavior. The

tools also can be applied to other aspects of rock

slope engineering whenever estimates of the

probability of failure are desired.

REFERENCES

1. Miller, S.M., J.M. Girard, and E. McHugh. 2000.

Computer modeling of catch benches to mitigate

rockfall hazards in open pit mines. In Pacific Rocks

2000 (Proc. of Fourth N. Amer. Rock Mech. Symp.),

eds. J. Girard, M. Liebman, C. Breeds, T. Doe, 539-

545. Rotterdam: Balkema.

2. Abramson, L.W., T.S. Lee, S. Sharma, and G.M.

Boyce. 2002. Slope Stability and Stabilization Methods,

2nd Ed., 427-454. New York: John Wiley & Sons,

Inc.

3. Miller, S.M. 1982. Fourier analysis for estimating

probability of sliding for the plane shear failure mode.

In Issues in Rock Mechanics (Proc. of 23rd U.S. Symp

on Rock Mech.), eds. R.E.Goodman, F.E. Hueze, 124-

131. New York: SME-AIME.

4. Wolff, T.F. 1996. Probabilistic slope stability in

theory and practice. In Uncertainty in the Geologic

Environment: From Theory to Practice, 419-433.

ASCE Geotech. Speciality Publ. No. 58.

5. Call, R.D., J.P. Savely, and D.E. Nicholas. 1976.

Estimation of joint set characteristics from surface

mapping data. In Monograph on Rock Mechanics
Applications in Mining, eds. W.S. Brown, S.S. Green,

W.A. Hustrulid, 65-73. New York: SME-AIME.

6. Rosenblueth, E. 1975. Point estimates for probability

moments. In Proc. National Academy of Sci., 72:10

(Mathematics), 3812-3814.

7. Miller, S.M. 1988. Modeling shear strength at low

normal stresses for enhanced rock slope engineering.

In Proc. of 39th Highway Geology Symp., 346-356.

Salt Lake City: Utah DOT.

8. Hines, W.W. and D.C. Montgomery. 1980. Probability

and Statistics in Engineering and Management Science,

2nd Ed., 634 p. New York: John Wiley & Sons, Inc.

9. Jaeger, J.C. 1971. Friction of rocks and stability of

rock slopes. Geotechnique, 21: 97-134.

10. Call, R.D. and D.E. Nicholas. 1978. Prediction of step

path failure geometry for slope stability analysis. In

Proc. of 19th U.S. Symp. on Rock Mechanics, 8 p.

Lake Tahoe, NV.

11. Miller, S.M. 1984. Probabilistic rock slope engineering. Publ. No. GL-84-8, 75 p. USAE-WES, Vicksburg, MS.

12. Pahl, P.J. 1981. Estimating the mean length of discontinuity traces. Intl. J. Rock Mech. Min. Sci. &

Geomech. Abstr., 18: 221-228.