Caving 2018 – Y Potvin and J Jakubec (eds)

© 2018 Australian Centre for Geomechanics, Perth, ISBN 978-0-9924810-9-4

https://papers.acg.uwa.edu.au/p/1815_09_Hebert/

Three-dimensional simulation of cave initiation, propagation

and surface subsidence using a coupled finite
difference–cellular automata solution

Y Hebert Itasca Australia Pty Ltd, Australia

G Sharrock Itasca Australia Pty Ltd, Australia

Abstract
This paper outlines a new methodology for modelling caveability and subsidence using bi-directional coupling
between the continuum code FLAC3D and the cellular automata code CAVESIM. FLAC3D, using the CaveHoek
constitutive model, simulates the progressive failure and disintegration of the rock mass from an
intact/jointed to a caved material. CAVESIM simulates gravity flow, in particular the collapse, bulking and
movement of caved rock. The coupled method captures many important aspects of caveability affecting cave
design such as hang-up formation, material recovery, timing of surface breakthrough or interaction with other
lifts, crater development, and surface subsidence. The key to improved modelling of many of these aspects is
the ability to accurately capture the impact of draw and gravity flow on cave propagation and subsidence.
Keywords: numerical modelling, caveability, subsidence, gravity flow

1 Introduction
Caveability affects many important aspects of cave design such as hang-up formation, material recovery,
timing of surface breakthrough or interaction with other lifts, crater development, and surface subsidence.
Prediction of caveability and subsidence is undertaken with empirical and numerical methods. Accurate
simulation of the caving process is complex, due to the wide range of mechanisms to be captured including
stress redistribution around the cave, rock mass failure in advance of the cave, reduction in strength from
peak to residual, dilation and bulking, material flow in the cave column, crater development, and slope
failure. The key to improved modelling of many of these mechanisms is capturing the impact of draw and
gravity flow on cave propagation.
As part of the Mass Mining Technology (MMT) project and over the past 20 years, an approach to estimating
rock mass material properties and simulating cave initiation and propagation has been developed. This paper
introduces new developments in the modelling method, whereby gravity flow in the cave column is
accounted for using a coupling technique between the finite difference program FLAC3D and the cellular
automata program CAVESIM. This allows for an improved representation of the impact of gravity flow and
production schedule on caveability and subsidence. In addition, the coupling technique combined with new
free rilling logic in CAVESIM has led to a new and improved way to estimate crater development.

2 Conceptual model of caving

In developing and interpreting numerical simulations of caving, it is useful to refer to a conceptual model of
the cave. The conceptual model of Duplancic and Brady (1999) has been employed to-date within the
International Caving Study (ICS) and MMT projects. In this model, a propagating cave consists of five main
behavioural regions as described below and illustrated in Figure 1(a).
The characteristics of each region are defined as follows:
 Pseudo-continuous domain or elastic region: The host rock mass around the caving region behaves
mainly elastically and has properties consistent with an ‘undisturbed’ rock mass.

Caving 2018, Vancouver, Canada 151

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

 Seismogenic zone: Microseismic (and sometimes seismic) activity is concentrated in this region
primarily due to slipping along pre-existing discontinuities and the initiation of new fractures.
 Zone of loosening or yielded zone: The rock mass in this region is fractured and has lost some or all
of its cohesive strength and provides minimal support to the overlying rock mass. Rock mass within
the yielded zone is subject to significant damage, i.e. open holes are cut off, time-domain
reflectometer (TDR) breakages are expected, and cracking is observable in infrastructure. Stress
components within this region are typically low in magnitude.
 Airgap.
 Caved zone or mobilised zone: This zone gives an estimate as to the portion of the orebody that has
moved in response to the production draw, and may be recoverable.
When considering surface subsidence, different terminology is often adopted. A review of large-scale surface
disturbances from block and panel caving mines was conducted by Lupo (1998), who found that the primary
surface features that develop as a result of block and panel caving include one or more of the following zones:
 A caved rock zone (crater).
 A large-scale surface cracking (fractured) zone.
 A small-scale surface displacement (continuous surface subsidence) zone.
 A stable zone.
Van As et al. (2003) proposed the terminology presented in Figure 1(b) to standardise the description of
subsidence features related to block and panel caving.

(a) (b)
Figure 1 (a) Main behavioural regions of a propagating cave (after Duplancic & Brady 1999); and,
(b) Terminology used to describe subsidence features for block and panel cave mines (after
Van As et al. 2003)

3 Numerical model of caving

3.1 Critical factors influencing caveability

Pierce (2010) has defined six critical factors that impact the caveability of a rock mass. They are presented
schematically in Figure 2 and include:
 Cohesion and tension weakening: During caving, rock mass undergoes a reduction in strength from
its peak in situ value to a much lower residual value representative of a caved rock mass state. This

152 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

overall response is often termed a strain-softening process and is the result of strain-dependent
material properties.
 Post-peak brittleness: The rate at which the strength drops from peak to residual is referred to as
brittleness. Rocks that maintain their peak strength with continued loading are referred to as
perfectly plastic (ductile). Rock masses that instantaneously drop to residual strength properties
when they exceed their peak strength are referred to as perfectly brittle. In general, brittle rock
masses cave more readily than more ductile rock masses.
 Modulus softening: During caving, the rock mass increases in volume as intact rock blocks fracture,
separate and rotate during the yielding and mobilisation process. Along with this bulking, a
reduction in modulus is expected to occur. Representation of the decrease in stiffness is crucial for
assessing the evolving stress state around the cave, since, as rock mass dilates/bulks, its potential
to carry stress decreases.
 Dilatational behaviour: Dilation is the change in volume of a rock that occurs with shear distortion.
An accurate assessment and representation of the dilational behaviour of a jointed rock mass
during caving is essential in being able to accurately predict the correct volume increase during
plastic deformation. This will affect mass balance calculations and impact propagation through
confinement levels.
 Rock mass structure (faults, joints): The influence of geological structures on subsidence and caving
is recognised to be important. The presence of a steeply dipping fault can terminate the angle of
draw short of its normal value. If a gently dipping fault intersects the collapsing rock column, the
lateral extent of surface subsidence can increase outward to the place where the fault intersects
the ground surface. Using experimental data, Brown and Trollope (1970) demonstrate that even if
the rock mass structure and geomechanical properties are well conditioned for slip, the existence
of high horizontal in situ stresses can inhibit rock mass failure and, therefore, prevent block caving
initiation and development.
 Production draw schedule: The rate of draw and shape of the mining footprint has previously been
identified by Laubscher (1990, 1994) to impact the caveability of a rock mass.

Figure 2 Schematic diagram of six critical factors that impact caveability as defined by Pierce (2010)

The caving algorithm described in this paper captures the effect of all six factors listed above to simulate cave
initiation, propagation and surface subsidence.

Caving 2018, Vancouver, Canada 153

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

3.2 Modelling approach

The caving and stress-redistribution processes inherently involve large deformations, shear along pre-existing
joints and bedding surfaces, fracturing of intact rock blocks and fragmentation of the rock mass. Ideally, one
would model this process using a pure discontinuum approach with explicit breakage of rock blocks. However,
the computational size and time requirements to solve mine-scale problems make it currently impossible to
tackle the problem completely with the discontinuum approach (Sjöberg et al. 2017). Instead, an algorithm to
simulate caving has been developed within the concept of a continuum-based model.
Material flow within the cave as a response to draw is an important aspect of the caving process. The discrete
element method (DEM) is by far the most promising and widely used simulation method in granular
science. However, it is difficult to apply to mine-scale flow problems at the full resolution of the
fragmentation, largely because of computational limitations. Cellular automata (CA) methods are
well-suited to extremely large and complex problems in discrete systems. CA permits very large systems
to be studied while including a more limited physics. For this reason, CA has gained wide acceptance in
science and engineering fields concerned with discrete mechanics.
The caving process described in this paper is represented by a bi-directional coupling method between the
continuum code FLAC3D and the CA code CAVESIM:
 FLAC3D is a three-dimensional explicit finite difference program for geotechnical analysis that can
simulate the behaviour of continuous materials that undergo plastic flow when their yield limits are
reached. Materials are represented by polyhedral elements within a three-dimensional grid. Each
element behaves according to a prescribed linear or non-linear stress–strain law in response to
applied forces and boundary restraints. The explicit, Lagrangian calculation scheme and the
mixed-discretisation zoning technique ensure that plastic collapse and plastic flow are modelled
accurately. FLAC3D can also simulate discontinuities such as faults, joints, bedding planes, and
engineered boundaries along constructs.
 CAVESIM is a general purpose CA code for the simulation of gravity flow of caved rock in block and
sublevel caves. The accurate representation and control of material bulking from the in situ to fully
bulked state allow the code to be used for calculation of accumulated airgap volume. In addition,
material flow against static faces and around or against disturbing features such as cave sidewalls,
overhangs, hang-ups and pendants can be modelled (Figure 3). In this method, the problem space
is firstly discretised into a fixed cubical lattice of dimension n.

Figure 3 Illustration of cave flow modelling using CAVESIM (after Sharrock et al. 2012)

154 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

3.3 Constitutive model

Within the continuum model, the constitutive relations provide a mathematical description of the rock mass
behaviour. As part of the ICS (Lorig 2000) and MMT (Chitombo & Pierce 2012) projects, several guidelines
have been developed to capture the key mechanisms associated with caving in a numerical model and are
summarised by Pierce (2013). The CaveHoek constitutive model has been developed by Itasca based on these
guidelines and aims to simulate the progressive failure and disintegration of the rock mass from an
intact/jointed to a caved material.
The peak-strength envelope is described by the Hoek–Brown failure criteria (Hoek et al. 2002). The numerical
implementation of the constitutive model uses a linear approximation, whereby the non-linear failure
surface is continuously approximated by the Mohr–Coulomb tangent, at the current zone stress level,
σ3 (Figure 4). The residual strength is typically that of a bulked rockfill (zero cohesion and a friction angle of
40–45°). As the material yields, the strength is degraded by linear interpolation between the linear
Mohr–Coulomb fit at peak and residual state based on the current stress state.

Figure 4 Rock mass failure envelope and strain-softening behaviour defined in the CaveHoek constitutive
model

Accumulated plastic shear strain (more specifically, the second invariant of the deviatoric plastic strain
tensor) is a common metric for irreversible shear strains in geomaterials and, in a more general sense, can
be considered as a measure of damage. The plastic shear strain required in going from peak to residual
strength (critical plastic strain, εpscrit) defines the brittleness of the rock mass failure, and impacts both the
caveability of a given unit, as well as the rate at which a cave will propagate for a given amount of draw.
An estimate of the relation between the critical strain and geological strength index (GSI) was determined by
a back-analysis of rock mass failure in caves and other openings (Lorig & Pierce 2000).
𝜀𝑝𝑠𝑐𝑟𝑖𝑡 = (12.5 − 0.125. 𝐺𝑆𝐼)/(100. 𝑑) (1)
where:
d = zone size.
During caving, the rock mass increases in volume (or bulks) due to dilation under shear, or due to volumetric
expansion under tension. In order to prevent zones from expanding to unrealistic levels, the dilation angle is
set to zero when the maximum porosity, i.e. maximum bulking potential, is reached. During bulking, a rock
mass is also expected to experience a significant drop in modulus. Representation of this drop is required in
order to account for stress shedding away from the cave and into the surrounding rock mass. In addition, it
is necessary to allow for modulus hardening that can occur when stresses are shed back onto exhausted or
undrawn parts of the cave. The rock mass modulus can be estimated from GSI, and the intact modulus by
using the relation developed by Hoek and Diederichs (2006). Based on laboratory testing results, Pappas and
Mark (1993) developed a set of equations to calculate the reduction in modulus as a function of the bulking
factor, uniaxial compressive strength of the intact rock fragments, and the aspect ratio of the fragments.

Caving 2018, Vancouver, Canada 155

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

3.4 Cellular automata flow simulation

CAVESIM is a multithreaded, general purpose CA code for the simulation of gravity flow of caved rock in cave
mines. Coupled flow-deformation simulations often require the study of gravity flow in complex multi-cave
mining environments with existing voids, historic excavations such as open pits and open stopes, and
complex draw rules or flow behaviour. Within the FLAC3D–CAVESIM coupled methodology, the outcomes of
alternative extraction geometries and draw strategies on caveability and subsidence can be assessed. This is
achieved through accurate computation and tracking of extracted tonnes, grade, dilution (rock types),
bulking, and particle fragmentation. Furthermore, crater slope stability angles can be pre-defined within
geotechnical domains. This allows the emergence of craters with different critical state slope angles, which
can have a significant impact on subsidence. The accurate representation and control of material bulking
from the in situ to fully bulked state allows the code to be used for calculation of accumulated airgap volume
and changes in cave geometry, which can be fed into risk assessments and cave geometry estimation
processes. In addition, material flow against static faces and around or against disturbing features such as
cave sidewalls, overhangs, hang-ups and pendants can be modelled. CAVESIM can be used to construct
residual block models of caved material for use in multi-cave analyses or grade reconciliations.
Methods have been developed to interpolate block model grades into CAVESIM simulations, providing a
realistic spatial distribution of grade as an initial simulation condition. At the commencement of each
simulation, every particle is assigned a grade, which is then integrated at the instant of drawing to generate
an estimate of dilution and grade at each drawpoint. CAVESIM has two analytical engines (stochastic and
Newtonian), both based on the fundamental lattice cellular automaton concept first developed by John von
Neumann in 1947 (von Neumann 1966).
The stochastic engine treats gravity flow as a stochastic process; that is, “a process that can be described by
a random variable that depends on some stochastic parameter which may be discrete or continuous”
(Borowski & Borwein 1991). The stochastic parameter is “the probability of downward propagation of a
particle”, or, “upward propagation of a void”. Examples include Nilsson (1988), Bergmark (1975), and Heden
(1976). The physics of granular flow are almost completely ignored in stochastic methods (Nedderman 1992).
For example, fundamental physical laws such as Newton’s laws and cohesive-frictional behaviour are not
typically solved in stochastic simulations and particle ‘contacts’ do not explicitly exist. For example, lateral or
sideways movement of particles in a stochastic model are often constrained to a fixed direction and governed
by simple kinematic or probabilistic rules. Complex emergent flow phenomena such as arching, slow infilling
of voids, inrushes and piping cannot be explicitly simulated in stochastic models. However, emulation is
possible. The key to stochastic modelling is developing rules that provide valid emulation of physical
processes. Care must be taken to stay within the limits of these rules or unrealistic behaviour can result.
In the Newtonian engine, the problem space is firstly discretised into a fixed cubical lattice of dimension L.
Traditional approaches to modelling granular media with CA have assumed that particles are perfectly
centred on the nodes of a fixed, two-dimensional lattice. This lattice is commonly composed of equilateral
triangles, since the distance between each node on such a lattice is the same (Baxter & Behringer 1990).
Hence, in these traditional approaches, particles do not undergo smooth translations between lattice nodes,
but rather jumps, or step-wise translations. These step-wise translations create difficulties when solving the
equations of motion at the lattice nodes. One such difficulty is the possible occurrences of artificial shock
waves in the solutions. For this reason, Fox et al. (1994) introduced a method to reduce the problem of
stepwise translations, and hence create smoother solutions. In this approach, a particle still occupies a
specific node, but the actual location of the centre of mass of a particle is defined by a position offset from
that node which can change incrementally with time.

156 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

3.5 Simulation of cave initiation and propagation in a continuum model

An early attempt to simulate the caving process using a continuum model was performed by Lorig (2000) as
part of the International Caving Study (ICS 1), using a two-dimensional axisymmetric model in FLAC3D
(Figure 5(a)). A strain-softening Mohr–Coulomb constitutive model was used to represent the rock mass
behaviour, and the extraction process at the undercut level was simulated by reducing the vertical pressure
p(t) at this level (Figure 5(b)). For a specific undercut geometry and rock mass properties, the caveability
condition could be assessed by monitoring the cave height with continuous reduction of the undercut
support p(t). Using this approach, the hydraulic radius associated with cave initiation compared well to
Laubscher’s empirical cave stability chart over a range of mine rock mass rating (MRMR) values.
However, the cave shapes obtained with this method (flat cave back) did not reproduce actual shapes.

(a) (b)
Figure 5 Schematic representation of (a) The axisymmetric model; and, (b) The evolution of the undercut
pressure and cave height with numerical steps (Lorig 2000)

The methodology adopted to simulate the caving process in the continuum model is based on the principle
described above, adapted to a full 3D model and modified to include a better representation of the cave
propagation phase in order to produce more realistic cave shapes. By monotonically reducing the support
pressure at the undercut level, cave initiation occurs when the critical hydraulic radius is reached for the rock
mass located immediately above the undercut. During cave propagation, and to simulate the reduction in
support provided by the muck pile on the cave back as a result of material flow inside the cave, the vertical
pressure at the boundary of the entire cave (Figure 6) is reduced during the simulation.

Figure 6 Example of cave boundary identified in the FLAC3D model during cave propagation

Caving 2018, Vancouver, Canada 157

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

Typically, areas of the cave with high flow provide less support to the cave back than static areas with little
to no flow. To simulate this effect, a local reduction factor (RFlocal) is applied to the vertical stress at the cave
boundary, meaning that the reduction in support pressure is non-uniform along the boundary of the cave
(Figure 7). RFlocal is calculated as a function of the volumetric flow rate provided by the CAVESIM model.
Controlling the pressure applied by the caved material on the cave back directly at the cave boundary thus
allows for the local effects of flow on cave propagation to be taken into account.
As part of the flow calculations, airgap development is also simulated within the CAVESIM model. In the
continuum model, the airgap is simulated by reducing the support pressure provided by the air on the cave
back in the same way as described above.

Figure 7 Illustration of support pressure reduction during the simulation in different areas of the cave
boundary, as a function of the volumetric flow rate calculated in CAVESIM

3.6 Bi-directional coupling between continuum and flow models

Production advance is simulated in CAVESIM, where material is drawn from the drawpoints based on the
input draw schedule. Material flow is simulated within CAVESIM, and the resulting airgap (volume and
location) and volumetric flow rate in the cave are communicated to the FLAC3D model. Within FLAC3D, the
stress state at the cave boundary is reduced to account for the airgap development and the reduction in
support provided by the muck pile during flow. As a result, stress redistribution occurs leading to propagation
of the yielded zone in the model. The new caved zone (material at residual strength with a vertical
displacement greater than 1 m) is identified and sent to CAVESIM. This cycle is repeated until the end of the
mining sequence (Figure 8). Communication between the two codes is performed using socket connections
(as used for TCP/IP transmission over the Internet), which allows fast transmission of data in binary with no
loss of precision.

158 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

CAVESIM

Draw and
Flow Calculation

Air Gap Volume

New Cave Shape
Flow Rate in the cave

Cave Propagation

FLAC3D

Figure 8 Conceptual illustration of the bi-directional coupling method to simulate the caving process

Ideally, the continuum and flow models should be run in parallel and at the same resolution, small enough
to capture in detail the different mechanisms associated with both cave propagation and flow. However, the
computing intensity required would be too intense to solve real problems in a reasonable amount of time.
Typically, a lattice cell size of 1 to 2 m gives a good representation of the flow behaviour using a CA method.
However, this resolution is difficult to achieve in continuum models at the mine scale due to the computing
intensity required to solve the strain–stress relations. Two simplifications are made in the coupling algorithm
to reduce the simulation time.
The mining sequence is divided into periods of time, and communication between the flow and continuum
phase occurs at the end of each period. Depending on the accuracy required, time periods ranging from one
day to one year can be defined. Back-analysis of actual cave propagation undertaken by Sharrock et al. (2012)
using this method indicate that weekly or two-weekly periods provide a good compromise between accuracy
and speed. It is noted that time is not simulated explicitly in the models and is only tracked during the
simulation as a function of production in CAVESIM. In the continuum model, mass balance calculations are
carried out during the simulation to ensure that the cave volume never exceeds the maximum volume
allowed for a given extracted mass and bulking factor. During a given period, cave propagation in FLAC3D
stops when the maximum cave volume allowed is reached, or when the cave stalls.
To account for the difference in model resolution between the continuum and flow models, a cell mapping
technique has been developed whereby each lattice cell within the CAVESIM model is associated to one
FLAC3D zone based on its location in space (Figure 9). This methodology allows for improved representation
of blocks and faults in the flow phase. The volumetric flow rate in each lattice cell is transmitted to the
continuum model and averaged out in each zone.

Figure 9 Conceptual illustration of the cell mapping technique used to associate each lattice cell in the
CAVESIM model (right) with each FLAC3D zone (left)

Caving 2018, Vancouver, Canada 159

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

3.7 Surface subsidence

Crater development is simulated using a bi-directional coupling method between CAVESIM and FLAC3D
(Figure 10). The free-surface rilling logic in CAVESIM allows a rill slope to form at a pre-defined angle of repose
through the avalanching mechanism. This allows for crater development to be simulated within CAVESIM as
material is being drawn. The crater shape obtained after each production advance is then initialised within
the FLAC3D model in order to assess the stability of the slopes. If a volume of material becomes unstable in
the FLAC3D model (due to rock mass failure and/or slipping along a fault), it is released in the CAVESIM model
and allowed to flow into the crater.

Figure 10 (a) Conceptual illustration of the CAVESIM–FLAC3D coupling methodology to simulate crater
development; (b) Isometric view; and, (c) North–south section of an actual crater geometry
obtained using this method

4 Example application
A series of conceptual models have been developed to illustrate the capability of the algorithm to simulate
typical behaviours that can be expected during cave propagation. All the models consider a 200 × 200 m
footprint located 600 m below the ground surface. The pre-mining state of stress is assumed lithostatic, with
a constant ratio of horizontal to vertical stress (K0 = 2). Three sets of material properties, provided in
Figure 11, have been used to represent different material strengths (good, fair and poor rock mass strength).
Two production schedules have been developed for the purpose of this exercise. In both cases, production
is simulated in weekly increments, and the undercut is opened at a rate of 250 m2/day from the southwest
corner of the footprint. During each period, the same mass of material is drawn in both schedules, the only
difference being the repartition of draw:
 A uniform draw schedule, whereby each drawpoint is drawn uniformly at a rate of 50 t/day after
complete establishment of the undercut.
 A non-uniform draw schedule that emphasises draw in the southwest corner of the footprint.

160 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

A summary of the different model permutations is provided in Table 1.

Figure 11 Rock mass strength and properties adopted to simulate good, fair and poor rock mass strengths

Table 1 Summary of model permutations

Caving 2018, Vancouver, Canada 161

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

Figure 12 illustrates the evolution of the cave for each model permutation, from cave initiation to surface
breakthrough. Figure 13 provides a contour of major principal stress around the cave (output from FLAC3D)
and cumulative vertical displacement of the caved material (output from CAVESIM) after 18 months of
production on an east–west cross-section through the centre of the cave for all model permutations. In all
cases, cave initiation occurs after four months of undercutting, when 72 drawpoints are established
(hydraulic radius of 27 m). Following initiation, the caving rate and cave shape varies in each model. Figure 14
illustrates the evolution of the mobilised zone, yielded zone, and seismogenic zone on a vertical cross-section
in the early stages of cave propagation for model 5.
Models 1, 3, and 4; moderate draw ramp-up in fair and poor rock: Model 1 represents the aim of most caving
operations, namely successful caving of the ore column with a broad cave with a fully developed stress arch.
The cave shape is initially slightly asymmetric in the fair rock mass, in response to the draw profile and
ramp-up starting location. However, the relatively low and uniform draw strategy results in the successful
formation of a wide stress arch in the cave crown and vertical caving of the sidewalls with no overhangs.
Finally, breakthrough to surface is rapid through poor near-surface materials with surface breakthrough after
130 weeks of production. By contrast, model 3 represents the same draw schedule but with strong rock at
the base of the cave (east side) resulting in an ore overhang forming. High stresses develop on the east side
of the cave, but are still too low to allow cave propagation in the strong material for the span. The overhang
induces faster vertical propagation, cave narrowing and potential stall (or slow caving). The key driver is that
the cave sidewalls propagate at a slower rate than forced by draw, resulting instead in void diffusion to the
cave crown. In such a scenario, both the draw schedule and footprint geometry could be altered to either
eliminate the overhang, and/or reduce cave narrowing during propagation. The impact of a major fault can
be similar to an ore overhang, as demonstrated in model 4. Here, a steeply dipping fault results in local
de-stressing of the cave sidewalls and more rapid vertical cave propagation. In such scenarios, substantial
rilling of material is indicated by the flow phase, and draw should be reduced under the overhang to prevent
void diffusion and dilution ingress.
Models 2 and 5; aggressive draw ramp-up in fair and poor rock: In these cases, an aggressive start-up draw
strategy leads to rapid cave propagation in the southwest corner of the footprint, where more material is
being drawn. In model 2, the cave narrows as it propagates due to the cave sidewalls propagating at a slower
rate than forced by draw, resulting instead in void diffusion to the cave crown. As a result of the cave
narrowing, cave stall or slow caving could occur if a strong rock mass or a fault is encountered in the cave
crown, leading to insufficient span and low crown and sidewall stress states. This example shows the risks of
aggressive early draw resulting in cave narrowing and potential stall, or slow caving, and demonstrate the
importance of capturing gravity flow in caveability models. Another outcome of the aggressive draw strategy
is breakthrough to surface occurs 20 weeks earlier than in the base case scenario. Model 5 shows the impact
of an aggressive schedule, combined with an overhang and fault, resulting in the formation of a conical cave
geometry with a narrow crown, with high potential for cave stall or slow caving if a strong rock mass or fault
is encountered in the cave crown, leading to insufficient span. In this case, surface breakthrough is 40 weeks
or 30% earlier than a moderate draw option with no overhangs or faults.
The above examples demonstrate the significant impact of draw and rock mass strength/stress/structure on
cave shape and caveability risks such as stall/slow caving, stranded reserves in overhangs, and airgap/airblast
risks. In addition, if draw is not adjusted to account for non-vertical cave geometry, the potential exists for
ore loss from high levels of internal mixing, promotion of fines or mud ingress, or poor connections with
other caves. The coupled method allows assessment of both caveability and flow risks within the same
modelling framework.

162 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

Figure 12 Cave propagation for all the model permutations

Caving 2018, Vancouver, Canada 163

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

Figure 13 Contour of major principal stress around the cave and cumulative vertical displacement of the
caved material after one year of production on an east–west cross-section through the centre of
the cave for all model permutations

Figure 14 Evolution of the mobilised zone (caved rock zone), yielded zone and seismogenic zone in the
early stages of cave propagation on a vertical cross-section (model 5)

164 Caving 2018, Vancouver, Canada

 Cave mechanics and numerical modelling

5 Conclusion
A methodology to simulate cave initiation, propagation, surface subsidence and gravity flow has been
developed using a coupled finite difference–cellular automata solution. The algorithm captures important
aspects of caveability such as overhang formation, influence of the draw strategy on cave propagation and
material recovery, potential for cave stall, airgap development and associated airblast risks, timing of surface
breakthrough or interaction with other lifts, crater development, and surface subsidence. Capturing the impact
of draw and gravity flow on cave propagation and subsidence is the key to improved modelling of these aspects.
Mechanisms captured in the algorithm include progressive failure and disintegration of the rock mass from an
intact/jointed to a caved material through the CaveHoek material model, complex gravity flow inside the cave
with CAVESIM, influence of geological structures, and crater development with FLAC3D. A strength of the
method is the assessment of caveability, subsidence and flow within the same modelling framework.

Acknowledgement
The work presented in this paper is the result of many years of development within Itasca and is largely
inspired by the work of Loren Lorig and Matthew Pierce from Itasca Consulting Group. The authors thank Ian
Brunton from Itasca Australia for his help in developing the method and his invaluable knowledge on
caveability and gravity flow.

References
Baxter, GW & Behringer, RP 1990, ‘Cellular automata models of granular flow’, Physical Review A, vol. 42, no. 2, pp. 1017–1020.
Bergmark, J 1975, The Calculation of Drift Spacing and Ring Burden for Sublevel Caving, report no. RU 75-16, Luossavaara-Kiirunavaara
Aktiebolag, Kiruna, in Swedish.
Borowsk, EJ & Borwein, JM 1991, The Harper Collins Dictionary of Mathematics, Harper Collins, New York.
Brown, ET & Trollope, DH 1970, ‘Strength of a model of jointed rock’, Journal of Soil Mechanics and Foundations Division, vol. 96,
pp. 685–704.
Chitombo, G & Pierce, ME 2012, A Practical Guide on the Use of the MMT Tools and Methodologies - Block Panel and Sublevel Caving,
report to Mass Mining Technology project.
Duplancic, P & Brady, BH 1999, ‘Characterization of caving mechanisms by analysis of seismicity and rock stress’, Proceedings of the
9th International Congress on Rock Mechanics, vol. 2, A.A. Balkema, Rotterdam, pp. 1049–1053.
Fox, GC, Williams, RD & Messima, PC 1994, Parallel Computing Works, Morgan Kaufmann, San Francisco.
Heden, H 1976, ‘Physical models for describing gravity flow’, in H Heden (ed.), Proceedings of the Seminar on Sublevel Caving in
Kiruna, Luossavaara-Kiirunavaara Aktiebolag, Kiruna, pp. 1–9, in Swedish.
Hoek, E, Carranza-Torres, C & Corkum, B 2002, ‘Hoek-Brown failure criterion - 2002 edition’, in R Hammah (ed.), Proceedings of the
5th North American Rock Mechanics Symposium and the 17th Tunnelling Association of Canada Conference, University of
Toronto, Toronto, pp. 267–273.
Hoek, E & Diederichs, MS 2006, ‘Empirical estimation of rock mass modulus’, International Journal of Rock Mechanics and Mining
Sciences, vol. 43, pp. 203–215.
Laubscher, DH 1990, ‘A geomechanical classification system for the rating of rock mass in mine design’, Journal of the South African
Institute of Mining and Metallurgy, vol. 90, no. 10, pp. 257–273.
Laubscher, DH 1994, ‘Cave mining – the state of the art’, Journal of the South African Institute of Mining and Metallurgy, vol. 94,
no. 10, pp. 279–293.
Lorig, L 2000, The Role of Numerical Modelling in Assessing Caveabilty, Itasca Consulting Group Inc., Minneapolis, report to the
International Caving Study, ICG00-099-3-16.
Lorig, L & Pierce, M 2000, Methodology and Guidelines for Numerical Modelling of Undercut and Extraction Level Behaviour in Caving
Mines, Itasca Consulting Group Inc., Minneapolis, report to the International Caving Study.
Lupo, JF 1998, ‘Large-scale surface disturbances resulting from underground mass mining’, International Journal of Rock Mechanics
and Mining Sciences, vol. 35, no. 4–5, paper no. 25.
Nedderman, RM 1992, Statics and Kinematics of Granular Materials, Cambridge University Press, Cambridge.
Nilsson, L 1988, A Theoretical Model of Gravity Flow in Sublevel Caving, report no. 88-740, Luossavaara-Kiirunavaara Aktiebolag,
Kiruna, in Swedish.
Pappas, D & Mark, C 1993, Behaviour of Simulated Longwall Gob Material, Report of Investigations 9458, United States Department
of the Interior, Bureau of Mines, Washington DC.
Pierce, M 2010, CaveHoek Constitutive Model: Theory and Implementation, Itasca technical memorandum.
Pierce, M 2013, ‘Numerical modeling of rock mass weakening, bulking and softening associated with cave mining’, ARMA
e-Newsletter, Spring 2013, no. 9.

Caving 2018, Vancouver, Canada 165

Three-dimensional simulation of cave initiation, propagation and surface subsidence Y Hebert and G Sharrock
using a coupled finite difference–cellular automata solution

Sharrock, GB, Beck, DA, Capes, G & Brunton, I 2012, ‘Applying coupled Newtonian cellular automata-discontinuum finite element
models to simulate propagation of Ridgeway Deeps Block Cave’, Proceedings of MassMin 2012, Canadian Institute of Mining,
Metallurgy and Petroleum, Westmount.
Sjöberg, J, Perman, F, Lope Álvarez, D, Stöckel, B-M, Mäkitaavola, K, Storvall, E & Lavoie, T 2017, ‘Deep sublevel cave mining and
surface influence’, in J Wesseloo (ed.), Proceedings of the Eighth International Conference on Deep and High Stress Mining,
Australian Centre for Geomechanics, Perth, pp. 357–372.
Van As, A, Davison, J & Moss, A 2003, Subsidence Definitions for Block Caving Mines, Rio Tinto technical report, Rio Tinto Technical
Services.
von Neumann, J 1966, Theory of Self-reproducing Automata, University of Illinois Press, Urbana.

166 Caving 2018, Vancouver, Canada