1 INTRODUCTION

Thorough engineering analyses of mine slopes cut in discontinuous rock masses should include investigations of 

bench stability and the likelihood of retaining serviceable 

catch benches during the mine life. When bench stability 

is controlled primarily by rock failures that slide along

natural fractures (such as plane shears and tetrahedral

wedges), a stochastic computer analysis can be used to 

evaluate the probability of retaining specified catchbench widths. If the original slope geometry plan and 

blasting layout are intended to produce catch benches of 

a certain width, it is unlikely that such width actually will 

be retained after blasting and excavation when kinematically viable rock failure modes are present in the

benches. Consequently, rockfall hazard assessment and 

related slope stability safety issues must consider the 

predicted, operational catch-bench geometry and not

the originally designed, ideal geometry.

 In fractured rock masses, bench failures most commonly occur in the upper portion of the bench, because 

the fracture lengths required for failure are shorter here. 

That is, small plane shears or wedges typically break out 

along the crest of a bench, due to the higher probability 

that natural fractures will be long enough here to form 

kinematically viable failure modes. This characteristic is 

observed in mine benches, and it should be reflected in 

the probabilistic outcome of a bench stability analysis. 

Thus, the probability of retaining a full width on the

bench is not as high as the probability of retaining, say, 

80% of the original bench width. Probabilities of retaining bench widths increase as the specified widths decrease. Only occasionally do longer fractures occur and 

allow for larger failures to affect much of the bench face 

and severely diminish the catch width.

Minor bench instabilities and rockfalls adversely impact mine safety in two key areas. First, as failures 

break back along the top of a bench, storage capacity 

for holding rockfall debris is significantly reduced, and 

falling rock from above may not be caught and retained 

on the bench. Second, as rockfall debris spills onto the 

bench below, it reduces the storage capacity of that 

bench and may even trigger multiple-bench failures. 

Thus, the capability to predict the size of bench failures 

and their impact on the catch width would provide key 

input to designing the bench geometry (i.e. bench height, 

face angle, and width).

The National Institute for Occupational Safety and 

Health (NIOSH) Spokane Research Laboratory focuses on safety and health issues in the mining industry. 

A project began several years ago aimed at mitigating 

rockfall hazards in open-pit mines and quarries. One 

aspect of this project has involved the development of 

computer software to analyze bench stability and backbreak characteristics in a hazard-based, stochastic

framework. One computer program analyzes planeshear failure modes in a two-dimensional framework by 

simulating plane-shear fractures in the bench and then 

calculating the probability of stability for each one, as 

well as identifying the corresponding backbreak distance 

on the bench. By repeating the simulation many times 

for a given bench, the probability of retaining various 

bench widths can be estimated. Another computer program analyzes three-dimensional wedges by simulating 

fractures from two fracture sets and conducting a similar 

back-break analysis. 

A stochastic (or probabilistic) approach is needed 

because rock fractures in a bench face are numerous 

and too difficult to analyze individually. For the computer model to generate realistic fracture patterns and 

subsequent slope failure modes, representative geotechnical input information must be available. This information often can be obtained readily from a thorough 

geotechnical site investigation.

Computer modeling of catch benches to mitigate rockfall hazards in open 

pit mines

Stanley M. Miller

Geological Engineering Program, University of Idaho, Moscow, Idaho

Jami M. Girard

NIOSH Spokane Research Laboratory, Spokane, Washington

Edward McHugh

NIOSH Spokane Research Laboratory, Spokane, Washington

ABSTRACT: A computer analysis of bench stability has been developed to account for multiple occurrences of potential slope-failure modes in discontinuous rock masses. Bench-scale plane shears and tetrahedral wedges are simulated and stochastically analyzed to estimate the probability of retaining specified catch-bench widths. This geotechnical information is useful in designing bench configurations to improve pit-slope stability and help alleviate rockfall 

hazards.
2 GEOTECHNICAL INPUT

Stability analysis of rock failure modes requires information on the slope geometry, the physical properties of 

rock discontinuities that define the modes, and local environmental conditions (such as ground water pore pressure). Slope geometry is specified by the engineer,

based on actual field conditions or on a proposed slope 

design plan. Other input data usually must be obtained 

by geotechnical site investigation procedures.

2.1 Mapping and analysis of rock discontinuities

Geotechnical data collection methods, such as scan-line 

(detail-line) mapping and fracture-set mapping (Miller, 

1983), provide important information on fracture orientations, spacings, lengths, and roughness. Typical mapping sites in the project vicinity include natural rock outcrops (if the project is in initial development stages) or 

available rock slope cuts along roads or accessible mine 

benches.

The first step in analyzing such field data typically 

consists of plotting the poles to fractures on a lowerhemisphere stereonet in order to identify fracture sets, 

which appear as clusters of poles (Hoek and Bray

1978). The interaction of the proposed slope cut with 

the orientations of these fracture sets allow the engineer 

to identify potential slope failure modes (i.e. plane shears 

and wedges, for our particular study). It should be

noted that in any rock-slope stability evaluation, the general progression in the engineering analysis is:

1. Use fracture-set orientations to identify potential slope failure modes;

2. For the critically oriented sets, evaluate the likelihood of having sufficient fracture lengths to form 

kinematically viable failure blocks; and

3. For fracture sets with sufficient lengths, estimate 

the shear strength so that an engineering stability analysis can be conducted.

In the computer simulation of rock fractures in a set, 

one should strive to preserve the natural spatial dependence in fracture properties. Spatial covariance or semivariograms (Isaaks and Srivastava 1989) provide a statistical format for describing the spatial dependence in

fracture properties, which has been demonstrated by La 

Pointe (1980) and Miller (1979). Thus, instead of simulating fracture properties independently in space, the

measured spatial continuity can be incorporated using

methods described by Miller (1985). To conduct such 

a fracture-set simulation, each of the particular fracture 

properties need to be modeled by an appropriate semivariogram model using the “sill” (sample variance), the 

“nugget” value (i.e. the semi-variogram value at a separation distance of zero), and the spatial “range” of influence. A probability distribution model for each fracture 

property also is needed.

An essential input for the stability calculations is the 

mean length of fractures in a given set. An exponential 

pdf (probability density function) is assumed for fracture 

length, then the probability of a fracture being long

enough to form a viable failure path through the bench 

can be obtained directly from the exponential probability 

distribution. This pdf is a one-parameter distribution, 

being defined only by the mean value. See Section 3.1 

below. 

2.2 Shear strength

Shear strength along rock fractures typically is estimated 

in one of two ways: the JRC-JCS method proposed by 

Barton et al (1972), and by using laboratory directshear test data to describe either a linear MohrCoulomb failure envelope or a power-curve model

(Jaeger 1971). 

A general power-curve model has been adopted for 

use in the NIOSH bench analysis computer programs, 

given by the following expression:

τ = aσ

b

 + c (1)

where: τ = shear strength;

σ = effective normal stress; and 

a, b, c = model parameters.

This equation describes a power model with a yintercept. It reduces to a simple linear model when b 

equals 1.0, thus making “c” equal to cohesion, and “a” 

equal to the coefficient of friction (i.e. tanφ).

The variability of τ, given a predicted value of σ, also 

is needed in the bench stability analysis. Currently in the 

NIOSH codes, the shear strength is modeled with a 

gamma pdf with a standard deviation defined by a userspecified coefficient of variation (CV). This coefficient is 

given by:

CV = sτ / mτ or sτ = CV(mτ) (2)

where: sτ = standard deviation of τ; and

mτ = mean of τ given by Eq. (1).

Therefore, as the normal stress increases, so does the 

shear strength and so does the standard deviation of

shear strength. Typical values for the shear strength CV 

range from 0.15 to 0.35. Note that for small values of 

CV (i.e. less than 0.2), the gamma pdf begins to approximate a normal pdf. The key advantage in using a 

gamma pdf to describe shear strength is that this particular pdf is defined only for positive values, which means 

that τ in the computer analysis never can take on unrealistic negative values for low values of σ. Note that small 

normal stresses are common when analyzing small failure 

masses along bench crests.

In summary, the required geotechnical input needed 

for the NIOSH bench stability programs can be summarized as follows:
Bplane.exe (2-d analysis of plane shears)

Bench height (m) and width (m)

Number of back-break cells (typically set so cells 

are about 1-m wide)

Bench face angle (degrees)

Ground water height above bench toe (m)

Rock mass unit weight (tonne/cu.m): mean, sd

Fracture-set mean length (m)

Fracture-set dip (deg.): mean, sd, nugget value,

spatial range (no. of fractures)

Fracture-set spacing (m): mean, nugget value,

spatial range (no. of fractures)

Fracture-set waviness (deg.): mean, nugget value,

spatial range (no. of fractures)

(Note: waviness is the average dip minus the minimum dip of a fracture, and it represents a measure of

large-scale roughness)

Shear strength (tonne/sq.m) terms: a, b, c, CV 

Bwedge.exe (3-d analysis of wedges)

Bench height (m) and width (m)

Number of back-break cells (typically set so cells 

are about 1-m wide)

Bench face angle and dip direction (degrees)

Ground water height above bench toe (m)

Rock mass mean unit weight (tonne/cu.m)

The following input is needed for both the left fracture set and the right fracture set that form viable 

wedges:

Fracture-set mean length (m)

Fracture-set dip direction (deg.): mean, sd, nug- get 

value, spatial range (no. of fractures)

Fracture-set dip (deg.): mean, sd, nugget value,

spatial range (no. of fractures)

Fracture-set spacing (m): mean, nugget value,

spatial range (no. of fractures)

Shear strength (tonne/sq.m) terms: a, b, c, CV

(σ and τ expressed in tonne/sq.m)

3 STOCHASTIC MODELING CONCEPTS

The probability of retaining a specified bench width for 

given failure modes in a bench can be estimated by

simulating potential failure geometries and cataloging the 

back-break position of each one on the top of the

bench. Stability of a given failure geometry can occur 

two ways: 1) the failure length is not long enough to 

pass entirely through the bench, and 2) the failure length 

is long enough to pass through the bench, but sliding

does not occur (Miller 1983). The probability of stability for each geometry then is given by the sum of these 

two probability values:

Pstab = P(failure path not long enough) 

 + P(failure path long enough and no sliding)

Pstab = (1 - PL) + PL(1-PS) (3)

Thus, the probability of failure length and the probability 

of sliding must be computed for each potential failure 

mass generated in the bench simulation.

3.1 Probability of failure length

The probability that a given simulated fracture is long 

enough to pass entirely through the bench is computed 

as an exceedance probability using an exponential pdf 

model for the fracture-set lengths. The exponential cdf 

(cumulative distribution function) is a one-parameter cdf 

model given by (Devore 1995):

F(x) = 0, if x < 0

 = 1 - e-x/m, if x > 0 (4)

where: m = mean.

The length required for a through-going failure path 

for a plane-shear fracture is calculated by:

x = h/sin(D) (5)

where: h = vertical height of failure mass, as measured from the toe of the failure to the top of the bench; 

and

D = dip of failure plane (or wedge intersection line for wedge failures).

Thus, the probability that fracture length takes on a 

value greater than x is given by:

P(X > x) = 1 - P(X < x) = 1 - F(x) = 1 - (1 - e-x/m)

 = e-x/m = PL (6)

Example: for mean length = 1.6m and x = 3m,

P(X > 3) = e-3/1.6 = 0.153 = PL

In the case of three-dimensional wedges, which slide

along the line of intersection, the probability of length 

sufficient for failure is the joint probability that the left 

fracture is long enough and the right fracture is long

enough:

PL(wedge) = PL(left) x PL(right) (7)

After setting the length of the wedge intersection equal to 

x in Eq. (6), the corresponding PL(left) and PL(right) can 

be computed using the mean length for the left fracture 

set and the mean length for the right fracture set, respectively.
3.2 Probability of sliding

The probability of sliding for a given slope failure mode 

can be estimated by Monte Carlo methods applied to a 

limiting-equilibrium analysis, whereby a distribution (histogram) of safety factor values is generated by many repeated calculations using possible realizations of input 

values. The probability of sliding then is equal to the

fraction of safety factors that are less than 1.0. The

safety factor is defined as the ratio of resisting forces to 

driving forces, and a value of 1.0 indicates limiting equilibrium (i.e. the potential failure mass is just on the verge 

of sliding).

However, even after completing a Monte Carlo

simulation study using several thousand iterations, the resulting histogram of safety factors represents only one 

possible realization of the actual safety factor pdf. A 

slightly different distribution will result if the simulation is 

repeated using a different random seed starting value. 

Thus, questions always arise regarding the number of iterations to use and the repeatability of results.

Fourier analysis provides an alternative to Monte

Carlo simulation in estimating the probability distribution 

of the safety factor, provided that the safety factor equation can be written as the sum of independent pdf’s. A 

computationally efficient way to estimate the actual pdf 

of the safety factor relies on discrete Fourier methods, 

which take advantage of the computing speed of the fast 

Fourier transform (Miller 1982). As presented by Feller 

(1966), the sum of independent pdf’s in the “space”

domain is analogous to the product of their Fourier

transforms in the “frequency” domain.

For our case, if fracture shear strength is assumed to 

have a gamma pdf and fracture waviness is assumed to 

have an exponential pdf (which is a special form of a 

gamma pdf), then the output safety factor pdf can be 

described as a gamma pdf. The probability of sliding 

(PS) is computed by numerically integrating the area under the discretized pdf of the safety factor to the left of 

safety factor = 1.0. That is,

PS = P(SF < 1.0) (8)

Additional information on this analytical method

based on Fourier convolution of pdf’s was reported by 

Miller (1982).

4 BENCH STABILITY ANALYSIS

The concept of bench back-break cells is illustrated in 

Figures 1 and 2. For the plane-shear analysis (Fig. 1), a 

random starting point is selected near the bench toe and 

then fracture locations up the bench are simulated by

generating spatially dependent fracture spacings. Fracture dips and waviness values also are generated using 

spatial dependence and assigned to individual fractures 

previously located on the slope face.

Figure 1. Simulated plane shears (Miller 1983).

Figure 2. Simulated 3-d wedges (Miller 1983).

By simulating many realizations of a given bench, 

each of which contains multiple occurrences of the particular failure mode, the probability of stability for any 

specified back-failure cell can be estimated as follows 

(Miller 1983):

 N Ji

PCS = [(NT- N)/NT] + (1/NT) Σ { Π [(1 - PLj) | si

 i=1 j=1

 + PLj(1 - PSj) | si } (9)
where: PCS = probability of cell stability;

NT = total number of bench simulations;

 N = number of bench simulations that

have at least one failure path in the specified cell;

si = the i-th bench simulation;

Ji = number of failure paths in the specified cell for the i-th bench simulation;

PLj = probability of sufficient length for

the j-th failure path;

PSj = probability of sliding for the j-th

failure path;

For simulating three-dimensional wedges in a bench, 

a standard length along the bench face must be specified 

to define an area for probability accumulations. This 

length typically is set equal to the bench height to provide for “square” units that can be analyzed along the 

bench face (Fig. 2). The number and size of simulation 

windows depend on fracture set spacings, lengths, and 

on engineering judgment (Miller 1983).

5 IMPLICATIONS FOR OVERALL SLOPES

Results from a bench stability simulation study can be 

used to help select interramp slope angles and overall 

slope angles. Bench geometry has a direct influence on 

the overall slope angle as expressed in the following

form:

tan(A) = 1 / [(W/H) + (1/tanB)] (10)

where: A = overall (average) slope angle;

B = bench-face angle;

H = vertical height of bench; and

W = horizontal width of bench.

Example: for H = 15m, W = 8m, B = 64 deg.,

 A = arctan{1 / [(8/15) + (1/tan64)] = 44 deg.

If an overall steeper angle is desired, then the W:H ratio 

of benches must be decreased or the bench faces cut at 

a steeper angle. However, for typical applications, results from the stochastic bench simulations will guide the 

engineer in selecting the overall slope angle to minimize 

extensive loss of catch bench width and thus minimize 

subsequent rockfall hazards. Relationships between

bench geometry, catch-bench width, and overall slope 

angle can be displayed in graphs, which can be used to 

predict overall slope angles when the probability of retaining a given catch-width has been specified. Examples are shown in Figure 3.

A. Plane shears; bench angle of 0.5:1 (64 deg.)

B. Plane shears; bench angle of 0.2:1 (79 deg.).

 Figure 3. Examples of catch-bench stability graphs.

Failure of a total bench width (i.e. related to the 

probability of retaining a bench width of zero) will have 

significant effects on interramp or overall slope stability, 

because multiple-bench failures often result when sections of benches are lost due to excessive failures of underlying benches. Thus, the design of overall slope angles should provide for very high probabilities (greater 

than 0.95) of retaining bench widths less than 1 meter.

One way to interpret the probability of retaining a 

given bench width is to directly relate this probability to 

the percentage of bench run in a given segment of the pit 

wall that will display that given bench width some time 

after bench excavation. For example, a 0.80 probability 

of retaining 4m wide catch benches can be interpreted 
as expecting about 80% of the bench run (as viewed 

horizontally) to have retained widths of at least 4m. 

Recent experience in using these bench stability computer programs has indicated that the number of bench 

simulations for plane shears should be greater than 100 

in most cases to provide consistent probabilistic results 

when repeating the analysis for different random seed 

values. The user-interface input screen for Bplane.exe 

currently lists a default value of 100 simulations, with a 

maximum of 200 allowed. Due to the extensive computational effort in the wedge simulations, Bwedge.exe currently lists a default value of 50, with a maximum of 100 

allowed.

Estimated volumes of rockfall debris can be related 

directly to the probability of retaining given catch-bench 

widths. If the probability associated with a specified

back-break cell is calculated as Pi, then the average failure volume associated with Pi

 can be estimated as follows:

First, calculate h, the vertical height of an average

failure:

h = Ci(sinDsinB) / sin(B-D) (11)

where: Ci = back-break distance to center of cell 

with probability Pi;

 B = bench face angle; and

 D = average dip of plane-shears or average plunge of wedges in the simulation.

Then, calculate the unit-width area (i.e. the area associated with a 1-m increment along a run of bench):

Ai = 0.5h[Ci + (h/tanB)] (12)

The associated intact volume of rock prior to the failure 

is then: Vi = Ai(1m) cubic meters of rock per 1m run 

of bench. A bulking factor (usually, 0.15 to 0.25) then 

is multiplied by this volume to estimate the volume of

loose rock debris lost from the bench crest and which 

must be contained on the catch bench below. If this 

volume exceeds the expected storage volume on the

lower bench, then the debris can be expected to cascade farther down the overall slope.

6 SUMMARY AND CONCLUSIONS

Computer software for a PC platform has been developed to stochastically analyze rock slope stability, particularly aimed at benches in an open-pit mine. The 

analysis also could be used for large rock slopes constructed for civil projects. The computer programs

simulated potential plane-shear or wedge failure modes 

and calculates to probability of retaining specified widths 

on the affected catch benches. Probabilistic estimates of 

potential failure volumes also can be obtained from this 

analysis.

Such information is useful in the design and selection 

of bench geometries and overall slope angles so as to 

minimize rockfall hazards while maintaining adequate

catch-bench widths. Field studies are underway to 

evaluate and verify the results provided by this type of 

stochastic analysis. In situations where blast damage or 

highly-fractured rock will control the back-break of

benches, the analysis of structurally controlled failures 

(plane shears and wedges) will not be adequate to describe bench stability. Also, no allowance has been

made in these computer codes for tension cracks that 

may truncate the failure paths, so the stochastic results 

may be approximate.

REFERENCES

Devore, J.L. 1995. Probability and statistics for engineering and the sciences, 4 ed. Belmont, CA: 

Wadsworth, 793 p.

Feller, W. 1966. An introduction to probability theory and its applications, v. II. New York: Wiley, 

626 p.

Hoek, E. & J.W. Bray. 1977. Rock slope engineering, 2nd ed. London: Inst. of Mining & Metallurgy, 

402 p. 

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

rock slopes. Geotechnique. 21: 97-134.

La Pointe, P.R. 1980. Analysis of the spatial variation in 

rock mass properties through geostatistics. In Proc. 

of 21st U.S. Symp. on rock mechanics, Rolla, 

MO: 570-580.

Miller, S.M. 1979. Geostatistical analysis for evaluating 

spatial dependence in fracture set characteristics. In 

Proc. of 16th Intl. Symp. APCOM, New York:

SME-AIME: 537-545.

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

probability of sliding for the plane shear failure mode. 

In R. Goodman & F. Hueze (eds), Issues in rock 

mechanics. New York: SME-AIME: 124-131.

Miller, S.M. 1983. Probabilistic analysis of bench

stability for use in designing open pit mine slopes. In 

Proc. of 24th U.S. Symp. on Rock Mechanics, College Station, TX: 621-629.

Miller, S.M. 1984. Probabilistic rock slope engineering, 

Publ. No. GL-84-8. Vicksburg, MS: USAE-WES: 

75 p.

Miller, S.M. & L.E. Borgman 1985. Spectral-type 

simulation of spatially correlated fracture set properties. 

Math Geology 17: 41-52.