DEEP COAL LONGWALL PANEL DESIGN FOR STRONG STRATA: 


THE INFLUENCE OF SOFTWARE CHOICE ON RESULTS 


By Mark K. Larson, Ph.D., P.E.,1

 and Jeffrey K. Whyatt, Ph.D., P.E.1 

 1Spokane Research Laboratory, National Institute for 

Occupational Safety and Health, Spokane, WA. 

ABSTRACT 

Software is often used to construct design models used 

in longwall panel design. These models necessarily reduce 

the abundant variability found in nature to a simplified 

representation that, ideally, captures the relevant characteristics of ground response to mining. The choice of stress 

analysis software is an important step in the modeling 

process. The importance of this step can be easily overlooked, yet assumptions inherent in modeling software can 

have a decisive influence. An appreciation of this step is 

important for both practitioners and users of design model 

results, particularly when decision-makers are integrating 

results into design decisions or evaluating the adequacy of 

design specifications. 

This paper examines four different stress analysis 

programs or tools for analyzing stress around a single 

longwall panel at 610 m (2,000 ft) depth beneath 

overburden containing some strong strata. The four tools 

are: an empirical model that underlies the Analysis 

of Longwall Pillar Stability (ALPS) program; the 

three-dimensional displacement-discontinuity program 

MULSIM/NL; the three-dimensional flexible overburden 

program LaModel, which has largely superseded 

MULSIM/NL; and the two-dimensional volume-element 

program FLAC. 

This study examines key aspects of panel simulation, 

including stress transfer through the gob, stress concentration in abutment ribs, and stress transfer distance 

into abutments. MULSIM/NL and FLAC results were the 

most similar. LaModel results varied greatly in terms of 

peak stress and stress transfer, while ALPS produced the 

least peak stress. This study also examines the transferability of calibrated input properties between MULSIM 

and LaModel. For instance, displacement-discontinuity 

codes have been replaced, at least in part, by LaModel, yet 

considerable experience exists in calibrated models with 

older tools. Whether and how this experience can be incorporated into analyses using other tools is a subject of some 

controversy and must be approached carefully. In all cases, 

selecting models and properties appropriate to site conditions is extremely important. 

INTRODUCTION 

Stress analysis tools such as ALPS, MULSIM/NL, 

LaModel, and FLAC have been used for some time to 

evaluate mine layout design. However, each of these tools 

has underlying assumptions that affect calculated results. 

This study examines the relative performance of these 

tools in building a generic design model for a single deep 

panel in geology typical of western U.S. longwall coal 

mines. The National Institute for Occupational Safety and 

Health has been evaluating these tools as part of a research 

project on control of bump hazards in deep western mines. 

The questions posed for this study are as follows: 

• 	 Are these tools appropriate for deep western coal 

mines where the overburden includes one or more 

strong members? 

• 	 What impacts do underlying assumptions have on 

results? 

• 	 Are there ways to translate experience gained with 

one tool into input for another? 

GENERIC MODEL 

Larson and Whyatt [2009] used a generic site model to 

compare the response of ALPS, LaModel, and FLAC to 

cases with a strong strata member in the overburden 

(Figure 1). The model’s stratigraphic column is typical of 

deep western coal mines. Elastic properties and strength 

properties used for each member were within the range of 

those found in the Wasatch Plateau and Book Cliffs region 

of Utah [Agapito et al. 1997; Haramy et al. 1988; Jones et 

al. 1990; Maleki 1995; Maleki 1988; Maleki et al. 1988; 

Maleki 2006; Pariseau 2007] and are presented in Table 1. 

These properties were used in the FLAC models. Elastic 

properties for the displacement-discontinuity codes were 

averaged from these properties according to the equivalent 

stiffness method and the weighted thickness method 

described by Larson and Whyatt [2009]. Only a single 

panel with 244-m (800-ft) width was considered. The 

study was expanded to include results of MULSIM/NL for 

this paper. 

In this study, Salamon’s [1966] nonlinear reconsolidation model with shale gob parameters as determined by 

Pappas and Mark [1993] was used as the constitutive law 

for gob in FLAC. In the case of MULSIM/NL and 

LaModel, the linearly hardening gob model was used, with 

parameters fit to the shale gob model used in FLAC. 
 Figure 1.—Stratigraphic column of generic model with 


indicated thickness of members used in models. 


Table 1.—Estimated properties of materials as used in the generic and FLAC models 

Property Soft shale Sandstone Shale Coal 

Young’s modulus, GPa 10.3 34.6 13.8 3.45 

Young’s modulus, million psi 

Poisson’s ratio 

1.50 

0.35 

5.00 

0.25 

2.00 

0.35 

0.50 

0.30 

Density, kg/m3

Density, lb/ft3

Cohesion, MPa 

2,310 

144 

20.5 

2,310 

144 

33.8 

2,310 

144 

20.5 

1,280 

80 

7.09 

Cohesion, psi 2,970 4,910 2,970 1,030 

Friction angle, ° 30 25 30 30 

Dilation angle, ° 

Tensile strength, MPa 

5 

2.07 

5 

5.03 

5 

6.89 

5 

2.07 

Tensile strength, psi 300 730 1,000 300 

Ubiquitous joint angle, ° 

Ubiquitous joint cohesion, MPa 

0 

1.4 

— 

— 

0 

1.4 

— 

— 

Ubiquitous joint cohesion, psi 200 — 200 — 

Ubiquitous joint friction angle, ° 

Ubiquitous joint dilation angle, ° 

25 

5 

— 

— 

25 

5 

— 

— 

Ubiquitous joint tensile strength, MPa 0.83 — 0.83 — 

Ubiquitous joint tensile strength, psi 120 — 120 — 

TOOL CAPABILITIES AND ASSUMPTIONS 

Numerical stress analysis tools are used in most cases 

to evaluate mining layout designs. For any given case 

where there is an excavation and adjacent pillars and 

abutments, three key points of information are determined 

by that model. Figure 2 illustrates these three points: 

(1) the fraction of overpanel weight that is transferred to 

the abutment versus the gob, (2) the distance into the 

abutment that the overpanel weight is transferred, and 

(3) the peak stress in the remaining pillars or abutment and 

its location. 

 Figure 2.—Vertical cross-section across the width of a 

panel showing important concepts of stress redistribution 

resulting from excavation of a single panel. 

Each tool simplifies the problem with underlying 

assumptions. For any specific case, the user must evaluate 

whether the underlying assumptions of a tool are appropriate. If a software tool is used without considering the 

underlying assumptions, the tool may give erroneous 

results. A brief summary of each tool with its capabilities 

and assumptions follows. 
 

ALPS 

Mark [1987] used case studies to empirically calibrate 

a simple estimate of stress distribution around a retreating 

coal panel. This method, called ALPS, considers a longwall panel across its width in vertical cross-section. 

Analysis of Retreat Mining Pillar Stability (ARMPS) is its 

equivalent for room-and-pillar retreat mining [Mark and 

Chase 1997]. ALPS considers no geology, only width of 

the excavation and overburden height to determine a stress 

distribution on the coal seam. Load on pillars and their 

stability factors are then calculated. Mark [1990] defines 

the stability factor as the load-bearing capacity of the pillar 

system divided by the design loading. A database of stability factors for various cases and their classification of 

“satisfactory” or “not satisfactory” allows the user to compare the case at hand with many others. Thus, deficiencies 

in estimation are “corrected” by using experience to define 

the critical stability factor. 

The original ALPS did not take into account the condition of the roof. Molinda and Mark [1994] developed the 

Coal Mine Roof Rating (CMRR). ALPS stability information was then modified to incorporate roof conditions. 

Mark et al. [1994] noted that the line 

 ALPS SF = 1.76 - 0.014 CMRR (1) 

separated “successful” cases from “not successful” tailgate 


cases in 82% of the case histories in the database. 


Advantages of this tool are: 

•	 It is quick and easy to calculate the stability 

factors. 

• 	 It has a large database for comparison with other 

cases. 

The assumptions of this tool are: 

•	 Caving and load transfer fit a simple model. Figure 3 depicts supercritical and subcritical vertical 

sections showing a wedge volume with unit 

thickness of the overpanel weight that is transferred to the abutment. That wedge is defined by 

the angle, f. No differences in geology are directly 

considered. The overpanel strata between the triangles are assumed to cave, and their weight is 

fully supported by the gob. 

• 	 Mark [1990] found for six cases in the Eastern 

United States, f ranged from 10.7° to 25.2° and 

recommended that f be assumed as 21°. That 

assumption is constant in the ALPS database. 

 Figure 3.—Vertical cross-section across the width of mined 

panels showing geometry of supercritical (left) and subcritical 

(right) panels (after Heasley [2008a]). The subcritical geometry 

is used by Mark [1987] in ALPS. 

Maximum load transfer distance is represented by the 

following equation [Pariseau 2007]: 

 D = 9.3 H , (2)

where D = maximum load transfer distance (ft); 

and H = overburden height (ft). 

 The vertical stress profile on the seam is represented 

by 

 

 3L


σ
 

a =
 

 

 

 

 

 

( )2 D − s x
 , 
 (3)
 3 D
 

where (a = abutment stress distribution function; 

 x = distance from the edge of the panel; 

and Ls = total side abutment load. 

MULSIM/NL 

Boundary-element techniques, pioneered in the late 

1960s and early 1970s, have a computational advantage 

over volume-element techniques for problems in infinite or 

semi-infinite domains in that the system of equations to 

solve is much smaller for the same problem [Crouch and 

Starfield 1983]. However, the equations are not sparse, 

as with volume-element tools, meaning that there are not 

many zero coefficients in the system of equations that must 

be solved. The displacement-discontinuity method is a 

subset of the boundary-element method, which solves the 

problem of a discontinuity in displacement between opposite surfaces of a crack over a finite length in an infinite 

elastic medium. Such a solution can be applied to a tabular 

deposit, such as a coal seam, where the behavior of the 

deposit is simulated with the crack. The seam is represented by a grid, and each square or block in that grid is 

assigned its own set of properties, strengths, and constitutive law. Constitutive laws for seam elements include 

linear elastic, strain softening, elastic-plastic, bilinear 
hardening, strain hardening, and linear elastic gob. The 

off-seam material is isotropic elastic only. 

Crouch and Fairhurst [1973] developed software 

implementing the technique for mine structural analysis. 

St. John [1978] used the technique in his code, 

EXPAREA. Sinha [1979] used the technique in the 

development of three computer programs, one of which 

(MULSIM) analyzes cases of multiple, parallel seams. He 

included the ability to subdivide coarse blocks into finer 

mesh so that the scheme was computationally more efficient. Beckett and Madrid [1986, 1988] developed additional features for MULSIM (their version was called 

MULSIM/BM), such as additional seam materials like 

gob, pack walls, and cribs; graphical development of grids; 

and an increase in the number of coarse blocks and the 

number of blocks that could be subdivided into finer mesh. 

Donato [1992] converted MULSIM/BM to a PC environment. Itasca Consulting Group added the ability to consider multiple mining steps [Zipf 1992b]. Zipf [1992a,b] 

added nonlinear seam materials such as strain-softening, 

elastic-plastic, bilinear hardening, and strain-hardening 

(MULSIM/NL). 

MULSIM/NL can handle up to four parallel seams dipping at some angle to horizontal as specified by the user. 

Initial three-dimensional stress conditions and stress gradients with depth are specified. Once the system of equations 

is solved iteratively, the full stress tensor and displacement 

vector components are output for each element and at userspecified locations in the surrounding ground. 

The DOS version (or Zipf version) of the tool uses 

coarse mesh to streamline calculations. For problems in this 

study (68 × 68 coarse blocks with fine mesh in the middle 

50 × 50 of the coarse blocks), a stress threshold of 1 psi 

served as the equilibrium convergence criterion so that the 

number of iterations typically ranged from 250 to 350. 

The Windows version (or Heasley version) of the tool 

lacks the coarse mesh, but allows a 400 × 400 fine mesh. 

Meshes of this size were initially used in this study, but all 

iterations were stopped at 90. Use of this version in this 

study was eventually abandoned because of long run time. 

Advantages of this tool are: 

• 	 Calculation time with the Zipf version is relatively 

short. 

• 	 The full stress tensor is output. 

• 	 Nonlinear in-seam behavior is available, as mentioned earlier. 

Assumptions of this tool are: 

•	 Overburden and underburden behavior can be 

adequately simulated with a one-material, elastic 

medium. 

• 	 Interaction between an elastic off-seam material 

and a nonlinear seam will adequately and realistically simulate gob-roof displacement and caving 

 

behavior. In short, the roof does not fail and cave. 

Instead, elastic sag of the overburden and appropriate in-seam gob constants can provide realistic 

load to the gob. 

LAMODEL 

Heasley [1998] developed LaModel, a displacementdiscontinuity modeling tool that uses the thin-plate or 

lamination formulation [Salamon 1991]. Each layer is 

separated by parallel, frictionless joints in the overburden 

and underburden at even intervals specified by the user. 

Heasley also included the same nonlinear models used in 

MULSIM/NL. The frictionless joints make the overburden 

less stiff, increasing closure of excavated areas of the 

seam. Heasley found this formulation tracked stress distributions in a coal seam and provided a better match to 

subsidence observations than an elastic model [Heasley 

1998]. However, he cautioned that the user must calibrate 

for displacement or stress, and calibrating for one entity 

may not provide realistic results for the other [Heasley 

2008b]. 

Solutions are completed about the same or faster than 

the Zipf version of MULSIM/NL. Off-seam calculations 

take a much longer time. Horizontal stresses are not 

computed. In addition, the seam can only be horizontal. 

Overburden and underburden are assumed to consist of 

layers of elastic material interspersed with horizontal, 

frictionless, cohesionless joints. The result is an increase in 

mechanical flexibility and the amount of predicted surface 

subsidence compared to the MULSIM/NL model, which 

assumes elastic behavior of the overburden. 

Heasley [2008b] suggested a method for calibrating a 

LaModel simulation. He recommended adjusting the overburden Young’s modulus so that 90% of the load transferred to the abutment lies within 

D = 5 H , (4) .9 

where H = seam depth (ft). Then Heasley recommended 

using a layer thickness according to 

 

 

 

 

 

 

 

2E 12(1−ν 2) 5 H − d 

2 

, (5) 

E h ln(0.1) 

t = s 

where E = elastic modulus of the overburden; 

v	 = Poisson’s ratio of the overburden; 

Es	 = elastic modulus of the seam; 

h	 = seam thickness (ft); 

d	 = extent of the coal yielding at the gob 

 edge (ft); 

and H = seam depth (ft), as in Equation 4. 
If Es , v, h, H, and d are maintained constant, then 

Equation 5 can be reformed to find the product, tE, that is: 

 tE	 = constant (6) 

In fact, this product is the actual overburden “property,” i.e., calculated results will be the same as long as 

this product is constant. 

Advantages of the tool are: 

• 	 Calculation time is relatively short. 

• 	 Nonlinear in-seam behavior is available. 

Assumptions of this tool are: 

•	 Overburden and underburden behavior can be 

adequately simulated with a one-material, elastic 

medium with embedded frictionless, cohesionless 

joints. 

• 	 Interaction between an elastic off-seam material 

with embedded frictionless, cohesionless joints and 

a nonlinear seam will adequately and realistically 

simulate gob-roof displacement and caving 

behavior. In short, the roof does not fail and cave. 

Instead, elastic sag of the overburden and appropriate in-seam gob constants can provide realistic 

load to the gob. 

EQUIVALENCY OF MULSIM/NL AND LAMODEL 

The transition from an elastic overburden mass in 

MULSIM/NL to overburden with frictionless, cohesionless 

joints in LaModel begs the question whether MULSIM/NL 

and LaModel can produce equivalent or similar results for 

appropriate inputs. Gates et al. [2008, Appendix V], based 

on Heasley’s [1998] formulation, proposed defining 

equivalent properties by matching midpanel closure. Solutions for the two methods, assuming cracks rather than 

coal seams, are: 

q s

2 2 2 

h ( ) x = 4(1 − υ ) (L − x ) an

E , d (7)

q sl ( ) 12 1 − υ2

x =	 (L 2 − x 

2 ),

 t E 

( ) (8)

where sh = seam convergence of the homogeneous 

 (MULSIM/NL) case; 

sl = seam convergence of the laminated case; 

x	 = distance from the panel centerline; 

v = rock mass Poisson’s ratio; 

t = layer or lamination thickness; 

q = overburden stress; 

E = rock Young’s modulus; 

and L = half-width of longwall panel. 

Combining these, Gates et al. [2008, Appendix V] found 

3 E

t = homogeneous L 

4 E 1−υ2 , laminated 

(9)

or 

tElaminated = kEhomogeneous (10) 

Gates et al. [2008, p. V-1] present Equation 9 as a 

method for finding the “required thickness” for translating 

a calibrated elastic overburden modulus to LaModel overburden properties. More specifically, they state that Equation 9 “could be used to estimate properties that would 

equate the laminated strata behavior with the homogeneous 

rock mass used in other boundary element programs” 

[Gates et al. 2008, p. 115]. This suggests that Equation 9 

may provide equivalence beyond closure at centerline of 

the panel. Such an equivalence would be a valuable link 

between tools, even though limited to a particular panel 

width (2L). 

A simple test was devised to explore the extent of this 

equivalence. An elastic model was constructed for both 

MULSIM/NL and LaModel with layer thickness calculated according to Equation 9. Young’s modulus of coal 

was set high (207 GPa (30,000,000 psi)) to minimize seam 

contribution to closure. Figure 4 shows the stress and closure profiles for a case with panel half-width set at 122 m 

(400 ft), Poisson’s ratio set at 0.35, and Ehomogeneous = 

Elaminated = 10.3 GPa (1.5 million psi); and thus t = k = 

113 m (370 ft). At midpanel, LaModel calculated closure 

to be 9.8% higher than MULSIM/NL. This difference 

likely is a result of effects from element size and edges 

[Heasley 1998, 2009; Zipf 1992b]. However, the stress 

profile within 46 m (150 ft) of the panel differed significantly. Thus, this “equivalence” is extremely limited and 

should not be used if abutment stresses are a concern, 

which is usually the case. 

FLAC OR VOLUME ELEMENT 

The concept of volume-element discretization for 

stress analysis has been around for a long time, but has 

become increasingly important since the advent of computers. For example, the finite-element method was 

developed in the aerospace industry in the 1950s 

[Segerlind 1976]. Turner et al. [1956] are credited with 

being the first to use the method in solid mechanics in 

1956. Finite-difference concepts have also been long 

known to the mathematical world [Dahlquist and Björck 

1974]. Dr. Peter Cundall first used the technique in his 

FLAC computer code to solve solid mechanics problems 

specifically for geomaterials in 1986. 
 Figure 4.—A, Vertical stress profile; B, closure profile for 

a purely elastic case with large Young’s modulus for coal. 

The closures at centerline of panel were expected to be 

equivalent. 

FLAC, as with most volume-element codes, calculates 

stress for each discretized volume element and displacement at each grid point that defines the volume elements. 

Each element has its own constitutive law and properties, 

including the possibility of material yielding and plastic 

flow. With FLAC version 5.0, used in this study, care must 

be taken in choosing mesh size and strength properties. If 

an element exhibits too much localized plastic flow, 

numerical instability may result. To detect such an 

impending situation, FLAC stops calculations if any triangular subzone of an element has an area below a 

threshold fraction of the entire element volume. However, 

FLAC version 6.0 includes the capability of dynamic 

meshing, thus eliminating this instability under large plastic deformations. 

Volume-element models discretize all modeled space, 

so boundaries are needed to limit model size. As a result, 

model boundaries must be far away from the area of 

interest to avoid influencing results, but close enough to 

limit problem size. Generally, this method needs a large 

number of elements and, thus, the size of the model can be 

very large. 

Advantages of FLAC are: 

• 	 Each member of the stratigraphic column can be 

represented according to an appropriate constitutive law, specific elastic properties, and specific 

strength properties. 

• 	 Failure of elements is determined by the code, not 

by the user. 

• 	 The code has an embedded simple computer language, FISH, that permits the user much versatility 

in model construction, model running, and inclusion of user-defined constitutive laws. 

• 	 Complex boundary conditions and initial conditions can be input to the model. 

Assumptions of this tool are: 

• 	 In the case of a two-dimensional model, the vertical cross-section of the model is far enough from 

panel ends that a plain-strain condition exists. 

• 	 The first-order volume-element response (in the 

elastic formulation, terms with exponents greater 

than 2 are neglected) adequately represents material behavior. If beamlike behavior is important in 

an analysis, then the grid must be fine to calculate 

accurate deformations and stresses. 

GENERIC MODEL RESULTS 

The generic model study with each of the four tools 

considered estimated three aspects of stress redistribution 

resulting from panel mining, as shown in Figure 2. These 

are (1) the fraction of overpanel weight that is transferred 

to the abutment versus the gob, (2) the distance into the 

abutment that the overpanel weight is transferred, and 

(3) the peak stress in the remaining pillars or abutment and 

its location. 

Overpanel Weight Transfer to Abutment 

or to Gob 

Figure 5 shows the range of overpanel weight fraction 

transferred to the abutment for roughly “equivalent” 

properties. The range for MULSIM/NL was very small— 

approximately 0.94–0.96 for the whole range of overburden properties used. MULSIM/NL and FLAC results 

are similar. The larger range of FLAC results was likely 

caused by various degrees of failure in the set of models. 

The full range of LaModel results with the same 

overburden elastic properties (E, v) was less than that of 

MULSIM/NL. Stress transfer with LaModel using a 

lamination thickness on the order of the overburden 

thickness was clearly not the same as that calculated by 

MULSIM/NL. The fraction of overpanel weight transferred to the abutment by ALPS for f = 21° was 

approximately 0.73—clearly below that of FLAC and 

MULSIM/NL. The upper value of 38° was chosen to test 

the impact of increased bridging of strata on results. 

 Figure 5.—Chart showing the range of proportion of overpanel weight shifted to abutments for three numerical modeling tools and ALPS for a panel width of 244 m (800 ft) and an 

overburden thickness of 610 m (2,000 ft). The locations of 

numbers in the LaModel column represent average results 

for that layer interval, where the interval is in feet. 

Peak Stress and Location 

Figure 6 shows the vertical stress profile on the coal in 

the first 46 m (150 ft) from the edge of the abutment for 

two generic model cases.

 Figure 6.—Vertical stress profile on abutment for the 

generic model case of 61 m (200 ft) of roof sandstone. Gob is 

modeled. A, No immediate roof shale; B, immediate roof shale 

thickness is 15 m (50 ft). 

 Results from LaModel were calculated with 1.5-m (5-ft) elements with eight yield rings, 

while the results from MULSIM/NL were calculated with 

3.0-m (10-ft) elements and four yield rings. Results from 

FLAC were calculated using 0.76-m (2.5-ft) elements. 

While these element sizes are not the same, they represent 

the best results one can get for a single panel model within 

the constraints of each numerical tool. 

The peak stresses and their locations are significantly 

affected by element size. To adequately compare results 

using the Mark-Bieniawski formula (MULSIM/NL and 

LaModel) and an equivalent-strength, elastic-perfectlyplastic constitutive law (FLAC), the elements would need 

to be the same size and probably smaller (0.76-m (2.5-ft) 

or smaller). 

For cases where only the seam is extracted (no immediate roof shale caves (e.g., Figure 6A)), FLAC tends to 

calculate higher peak stress than MULSIM/NL and 

LaModel. However, FLAC, MULSIM/NL, and LaModel 

with layer thickness set at overburden height calculate 

similar peak stresses and stress profiles. 

FLAC, MULSIM/NL, and LaModel handle material 

failure near the rib differently. In FLAC, zones stressed 

beyond the elastic limit undergo plastic flow, thus reducing 

the stiffness near the rib and shifting overburden weight to 

stiffer coal farther from the rib. MULSIM/NL and 

LaModel impose an elastic limit according to the MarkBieniawski formula and allow more closure of the seam in 

an element if stresses are above the limit. These two 

methods are not equivalent and thus can affect the amount 

of peak stress calculated. It is, therefore, not surprising 

when FLAC peak stresses are significantly different from 

either MULSIM/NL or LaModel. ALPS, of course, has no 

ability to simulate failure of seam material near a rib. 

Where it is assumed that immediate roof shale caves and 

forms gob (FLAC results in Figure 6B), the location of the 

top of the seam with respect to the geometry of the 

opening affects the amount of peak stress. Such geometry 

is not possible in the boundary-element codes, where fullheight extraction of the seam only is assumed. 
Stress Transfer To Abutment 

Figures 7 and 8 show the range of locations in the 

abutment where total vertical stress returns to 150% and 

200%, respectively, of premining vertical stress for 

reasonable ranges of input for each tool.

 Figure 7.—Distance from gob where vertical stress on coal 

is 150% of premining vertical stress. Results include immediate shale thickness of 0, 3.0, and 15 m (0, 10, and 50 ft) and all 

LaModel layer intervals modeled. 

 Figure 8.—Distance from gob where vertical stress on coal 

is 200% of premining vertical stress. Results include immediate shale thickness of 0, 3.0, and 15 m (0, 10, and 50 ft) and all 

LaModel layer intervals modeled. 

 It is, in essence, 

a map of accessible solutions. LaModel lamination thickness significantly affects stress transfer distance. Results 

with lamination thickness of 7.6-m (25-ft) plot at the 

bottom of the LaModel range, while results with lamination thickness of 610-m (2,000-ft) plot at the top of the 

LaModel range. MULSIM/NL stress transfer distance 

seems to be at approximately midrange of the LaModel 

results. MULSIM/NL stress does not seem to transfer as 

far into the abutment as ALPS, FLAC, or LaModel with 

layer thickness set at the overburden thickness. However, 

FLAC results with smaller roof sandstone thickness plots 

near corresponding MULSIM/NL results. 

Equivalency of MULSIM/NL and LaModel: 


Generic Model Results 


Equivalence of MULSIM/NL and LaModel tools was 

also tested for the generic model to see if some combination of input parameters might produce essentially equivalent stress distributions. This was addressed by applying 

generic model E and varying t in LaModel. 

Figure 9 shows stress profiles for the generic case 

where immediate roof shale thickness is 0 m and roof 

sandstone thickness is 61 m (200 ft).

 Figure 9.—Vertical stress profiles across half-width of 

panel and abutment for MULSIM/NL and LaModel for the 

same overburden properties and various layer thicknesses 

for LaModel. Properties were determined by the weighted 

thickness method for shale = 0 m and sandstone = 61 m 

(200 ft). 

 Stress profile equivalency does not seem possible for this practical example 

with inelastic rib properties. When layer thickness was set 

to overburden thickness, peak stress calculated by the two 

tools was nearly the same. However, LaModel results 

showed more load closer to the opening than those of 

MULSIM/NL. It seems that no value of the constant tE 

will produce a LaModel stress profile equivalent to that 

from MULSIM/NL. 
Discussion 

These results show that the underlying assumptions of 

each of these stress analysis tools (ALPS, MULSIM/NL, 

LaModel, and FLAC) significantly influenced results. Of 

these, FLAC and MULSIM/NL with averaged overburden 

properties provide the most similar results. The main 

difference is that FLAC transfers overpanel weight slightly 

farther into the abutment. The recommended f = 21° for 

ALPS, which controls the portion of overpanel weight 

transferred to the abutment, underestimates abutment load 

transfer relative to FLAC and MULSIM/NL. Larger values 

of f provide closer results, but nonstandard input to ALPS 

makes comparison with the ALPS database invalid. 

LaModel produces the widest range of possible 

solutions and thus is most sensitive to input parameters. 

Heasley [2008b] addresses this by recommending a 

specific input development process. However, in deep 

western conditions where the stratigraphic column includes a strong, stiff member, overpanel weight may be 

transferred farther [Barron 1990; DeMarco et al. 1995; 

Gilbride and Hardy 2004; Goodrich et al. 1999; Kelly 

1999; Maleki 2006] than Equation 2 would indicate. In 

such a case, observations or measurements of load transfer 

would be necessary to determine the maximum load 

transfer distance. 

In this study, the gob model was not varied. However, 

accurate simulation of gob behavior is perhaps the input 

that most significantly affects stress transfer to the abutment. All four tools may be able to adequately simulate 

site-specific gob behavior where detailed observations or 

measurements have been made of strata behavior, subsidence, extent of cave, gob loading, etc. However, not all 

tools can simulate important mechanisms. For example, 

a strong, stiff overburden member may result in arching of 

stress over the excavated panel and subsequent sudden 

collapse. Displacement-discontinuity models cannot simulate that mechanism or its effects. Only models that 

simulate the behavior of individual stratigraphic members 

can simulate and/or predict this mechanism. 

When moving between MULSIM/NL and LaModel, 

it may be possible to get equivalent closures at midpanel 

by limiting error resulting from edge and element size 

effects, but equivalent stress profiles cannot be achieved in 

this case. Equivalence seems to be approached, but not 

attained, as lamination thickness is increased to 

overburden depth. The user’s choice between these two 

tools should be motivated by measured behavior or 

characteristics of the stratigraphic column, such as the 

strength of beds, bedding plane properties, thickness of 

beds, etc. 

CALIBRATION TO STRESS MEASUREMENTS 

Comparison to actual conditions is the only way to 

evaluate the validity of a model. It can also be used to 

“calibrate” model input, i.e., results of the previous section 

show only relative performance, not which is “best.” 

Figures 7 and 8 show ranges of stress transfer distance that 

each tool can achieve for a given geology. Therefore, it is 

best, where possible, to assess whether a tool is appropriate 

for a specific site. 

As an example, Larson and Whyatt [2009] compared 

stress in the abutment calculated with ALPS, LaModel, 

and FLAC with borehole pressure cell measurements of 

stress induced by mining at two sites. Figures 10 and 11 

show the same results with MULSIM/NL calculations 

added. Gate road entries were not included in the models 

for convenience. However, aside from local perturbations 

around entries, each tool can be optimized to provide the 

best possible approximation of the measured stress 

distribution. 

While the borehole pressure cells were not located 

close enough to the rib to capture the actual peak stress, 

the ALPS stress distributions do not simulate the rapid 

change in stress near the gob. Calibrated MULSIM/NL 

stress change profiles reasonably matched the changes in 

stress measurements. LaModel stress profiles varied 

greatly, with either the peak stress too high or the stress 

not decreasing as quickly as the measurements with 

distance from the gob. FLAC models were not built to 

closely simulate stratigraphy because stratigraphic members or their properties was not sufficiently described 

[Barron 1990; Koehler et al. 1996]. Instead, generic model 

geometries with actual panel widths and overburden 

heights were used to approximate actual site conditions. In 

these cases, calculated vertical stress profiles with roof 

sandstone thicknesses of 3 m (10 ft) and 15 m (50 ft) are 

reasonably close to the changes in stress measurements, 

suggesting that calculated results likely would match stress 

change measurements closely if the actual stratigraphic 

columns were modeled. 
 Figure 10.—Mining-induced stress around the 6th Right 

gate roads at a Utah Mine [Koehler et al. 1996]. A, ALPSassumed stress distribution function and measurements; 

B, MULSIM/NL and LaModel results; C, FLAC results with 

measurements. (BPC = borehole pressure cell; Sh = shale; 

SS = sandstone; YM = Young’s modulus.) 

 Figure 11.—Mining-induced stress around the 9th East 

gate roads at a Utah Mine [Barron 1990]. A, ALPS-assumed 

stress distribution function and measurements; 

B, MULSIM/NL and LaModel results; C, FLAC results 

with measurements. (BPC = borehole pressure cell; 

Sh = shale; SS = sandstone; YM = Young’s modulus.) 

The importance of calibrating tool input properties to 

site conditions is emphasized by these results. This point 

was also underscored in a recent Program Information 

Bulletin from the Mine Safety and Health Administration 

[Skiles and Stricklin 2009]. In the bulletin, an eight-step 

general process is outlined for successful calibration to 

specific sites. The bulletin states, in summary:

“Successful numerical simulation requires a substantial 

effort including the observation of in-mine conditions 

in many areas and the often repetitive process of 

calibrating model parameters...It cannot be overemphasized, however, that in order to be of value, 

a numerical model must be validated and provide a 

realistic representation of the underground environment 

for which it is applied.” 
CONCLUSIONS 

Empirical, boundary-element, and volume-element 

tools are often used to evaluate mine plans. These tools are 

not interchangeable. Site-specific information must be 

used to calibrate and assess the appropriateness of each 

tool. For instance, LaModel does not consider tectonic 

stress. Thus, it might be a poor choice if high horizontal 

stresses are present and can influence outcomes in the 

mine. For empirical tools, application site conditions 

should be compared with conditions at underlying cases to 

determine whether the tool might be justified. The published database includes a wide variety of geologic and 

stress conditions, but it is possible that cases exist outside 

of database conditions. If possible, additional cases based 

on local experience should be compared to the published 

database to establish validity. 

Generally, a tool’s input properties should be calibrated to site-specific conditions and observations. This 

point is even more important when a strong, stiff member 

is present in the stratigraphic column. Calibration procedures should include achieving the correct abutment stress 

profile and reasonable stress transfer distance. Empirical 

tools are often an exception, as the reference database typically includes a variety of site conditions. However, if 

local behavior or cases depart significantly from the 

underlying database, the user may need to make adjustments. These can include revising coal strength and cave 

angle, for example. Alternatively, a revised critical 

stability factor might be proposed. In either case, 

adjustments delink results from the underlying empirical 

database and the established success criterion, creating, 

in essence, a new empirical method that must be justified 

on its merits. 

Results of the generic model study show that assumptions and features of the tools studied differ too markedly 

for a model constructed with one tool to be “converted” 

into another through application of the same input 

parameters. Thus, care must be taken when taking input 

parameters from past analyses using different modeling 

programs. New models should be calibrated to field observations and, ideally, measurements of critical behavior. 

Results from this study show that if FLAC input 

accounts for individual stratigraphic members and their 

properties, then the seam stress profile that it calculates is 

most similar to results calculated by MULSIM/NL, where 

overburden properties were determined by averaging individual member properties by some reasonable method. 

However, the user should evaluate the stratigraphic column and the properties of its individual members with 

respect to the caving mechanism of the site. An analysis 

tool should be able to simulate the most important effects 

of that mechanism. 

LaModel results are highly sensitive to input parameters and, therefore, careful model calibration, adjusting tE 

to fit observed or measured stress conditions, is required. 

LaModel results did not fit the example cases very well. 

The quick decrease of stress near the rib requires tE to 

be relatively small, but this increases peak stress 

unrealistically. 

ALPS stress distribution significantly underestimates 

peak stress. However, using the default value of f may, 

in many cases, be a reasonable estimate of total load transferred to the abutment. 

In this study, we assumed a model for gob (shale from 

Pappas and Mark [1993]) and kept that constitutive model 

constant for all models. If a stiffer gob were used, less 

stress would be transferred to the abutment. Determining 

the correct gob stiffness relationship is one of the most 

important parts of model calibration because that relationship is the most significant factor in determining amount of 

load transfer to the abutment. 

The user must be very careful when selecting the 

proper analysis tool. If caving behavior can be affected 

significantly by variability in the strength and stiffness of 

stratigraphic column members, then a tool should be used 

that can take that influence into account. This means that 

tools that assume a single off-seam material that is elastic 

may miss details of off-seam strata behavior that may be 

significant. 

ACKNOWLEDGMENTS 

The authors would like to thank Keith A. Heasley, 

Ph.D., West Virginia University, Department of Mining 

Engineering, for discussions that strengthened this paper. 

REFERENCES 

Agapito JFT, Goodrich RR, Moon M [1997]. Dealing 

with coal bursts at Deer Creek. Min Eng 49(7):31–37. 

Barron LR [1990]. Longwall stability analysis of a 

deep, bump-prone western coal mine: case study. In: Peng 

SS, ed. Proceedings of the Ninth International Conference 

on Ground Control in Mining. Morgantown, WV: West 

Virginia University, pp. 142–149. 

Beckett LA, Madrid RS [1986]. Practical application 

of MULSIM/BM for improved mine design. In: Proceedings of the Third Conference on the Use of Computers in 

the Coal Industry (Morgantown, WV, July 28–30, 1986), 

pp. 209–219. 

Beckett LA, Madrid RS [1988]. MULSIM/BM: a structural analysis computer program for mine design. Denver, 

CO: U.S. Department of the Interior, Bureau of Mines, 

IC 9168. NTIS No. PB 88-237565. 

Crouch SL, Fairhurst C [1973]. The mechanics of coal 

mine bumps and the interaction between coal pillars, mine 

roof, and floor. Minneapolis, MN: University of Minnesota, Department of Civil and Mineral Engineering. U.S. 

Bureau of Mines contract No. H0101778. NTIS No. PB 

222 898. 
Crouch SL, Starfield AM [1983]. Boundary element 

methods in solid mechanics. London: George Allen & 

Unwin. 

Dahlquist G, Björck Å [1974]. Numerical methods. 

Englewood Cliffs, NJ: Prentice-Hall. 

DeMarco MJ, Koehler JR, Maleki H [1995]. Gate road 

design considerations for mitigation of coal bumps in 

western U.S. longwall operations. In: Maleki H, Wopat PF, 

Repsher RC, Tuchman RJ, eds. Proceedings: Mechanics 

and Mitigation of Violent Failure in Coal and Hard-Rock 

Mines. Spokane, WA: U.S. Department of the Interior, 

Bureau of Mines, SP 01-95, pp. 141–165. NTIS No. PB95­

211967. 

Donato DA [1992]. MULSIM/PC: a personal 

computer-based structural analysis program for mine 

design in deep tabular deposits. Denver, CO: U.S. Department of the Interior, Bureau of Mines, IC 9325. 

Gates RA, Gauna M, Morley TA, O’Donnell JR Jr., 

Smith GE, Watkins TR, Weaver CA, Zelanko JC [2008]. 

Report of investigation: underground coal mine, fatal 

underground coal burst accidents, August 6 and 16, 2007, 

Crandall Canyon mine, Genwal Resources, Inc., Huntington, Emery County, Utah, ID No. 42-01715. Arlington, 

VA: U.S. Department of Labor, Mine Safety and Health 

Administration. 

Gilbride LJ, Hardy MP [2004]. Interpanel barriers for 

deep western U.S. longwall mining. In: Peng SS, Mark C, 

Finfinger GL, Tadolini SC, Heasley KA, Khair AW, eds. 

Proceedings of the 23rd International Conference on 

Ground Control in Mining. Morgantown, WV: West Virginia University, pp. 35–41. 

Goodrich RR, Agapito JFT, Pollastro C, LaFrentz L, 

Fleck K [1999]. Long load transfer distances at the Deer 

Creek mine. In: Amadei B, Kranz RL, Scott GA, Smeallie 

PH, eds. Rock Mechanics for Industry. Proceedings of the 

37th U.S. Rock Mechanics Symposium (Vail, CO, June 6– 

9, 1999). Rotterdam, Netherlands: A. A. Balkema, pp. 517– 

523. 

Haramy K, Magers JA, McDonnell JP [1988]. Mining 

under strong roof. In: Peng SS, ed. Proceedings of the 

Seventh International Conference on Ground Control in 

Mining. Morgantown, WV: West Virginia University, pp. 

179–194. 

Heasley KA [1998]. Numerical modeling of coal 

mines with a laminated displacement-discontinuity code 

[Dissertation]. Golden, CO: Colorado School of Mines, 

Department of Mining and Earth Systems Engineering. 

Heasley KA [2008a]. Back analysis of the Crandall 

Canyon mine using the LaModel program. Appendix S in: 

Gates RA, Gauna M, Morley TA, O’Donnell JR Jr., Smith 

GE, Watkins TR, Weaver CA, Zelanko JC. Report of 

investigation: underground coal mine, fatal underground 

coal burst accidents, August 6 and 16, 2007, Crandall 

Canyon mine, Genwal Resources, Inc., Huntington, Emery 

County, Utah, ID No. 42-01715. Arlington, VA: U.S. 

Department of Labor, Mine Safety and Health Administration. 

Heasley KA [2008b]. Some thoughts on calibrating 

LaModel. In: Peng SS, Tadolini SC, Mark C, Finfinger 

GL, Heasley KA, Khair AW, Luo Y, eds. Proceedings of 

the 27th International Conference on Ground Control in 

Mining. Morgantown, WV: West Virginia University, pp. 

7–13. 

Heasley KA [2009]. Personal communication, March 

17. 

Jones RE, Pariseau WG, Payne V, Takenaka G [1990]. 

Sandstone escarpment stability in vicinity of longwall 

mining operations. In: Hustrulid WA, Johnson GA, eds. 

Rock Mechanics Contributions and Challenges: Proceedings of the 31st U.S. Symposium (Golden, CO, June 18– 

20, 1990). Rotterdam, Netherlands: A. A. Balkema, pp. 

555–562. 

Kelly M [1999]. 3D aspects of longwall geomechanics. 

In: Ground behaviour and longwall faces and its effect on 

mining. Exploration and Mining Report 560F, ACARP 

project C5017. 

Koehler JR, DeMarco MJ, Marshall RJ, Fielder J 

[1996]. Performance evaluation of a cable bolted yieldabutment gate road system at the Crandall Canyon No. 1 

mine, Genwal Resources, Inc., Huntington, Utah. In: 

Ozdemir L, Hanna K, Haramy KY, Peng S, eds. Proceedings of the 15th International Conference on Ground Control in Mining. Golden, CO: Colorado School of Mines, 

pp. 477–495. 

Larson MK, Whyatt JK [2009]. Critical review of 

numerical stress analysis tools for deep coal longwall 

panels under strong strata. SME preprint 09-011. Littleton, 

CO: Society for Mining, Metallurgy, and Exploration, Inc. 

Maleki H [1995]. An analysis of violent failure in U.S. 

coal mines: case studies. In: Maleki H, Wopat PF, Repsher 

RC, Tuchman RJ, eds. Proceedings: Mechanics and Mitigation of Violent Failure in Coal and Hard-Rock Mines. 

Spokane, WA: U.S. Department of the Interior, Bureau of 

Mines, SP 01-95, pp. 5–25. NTIS No. PB95-211967. 

Maleki H [1988]. Ground response to longwall mining: 

a case study of two-entry yield pillar evolution in weak 

rock. Colo Sch Mines Q 83(30):1–60. 

Maleki H [2006]. Caving, load transfer, and mine 

design in western U.S. mines. In: Yale DP, Holtz SC, 

Breeds C, Ozbay U, eds. Proceedings of the 41st U.S. 

Rock Mechanics Symposium (Golden, CO, June 17–21, 

2006). Alexandria, VA: American Rock Mechanics 

Association. 

Maleki H, Agapito JFT, Moon M [1988]. In-situ pillar 

strength determination for two-entry longwall gates. In: 

Peng SS, ed. Proceedings of the Seventh International Conference on Ground Control in Mining. Morgantown, WV: 

West Virginia University, pp. 10–19. 

Mark C [1987]. Analysis of longwall pillar stability 

[Dissertation]. University Park, PA: The Pennsylvania 

State University. 
Mark C [1990]. Pillar design methods for longwall 

mining. Pittsburgh, PA: U.S. Department of the Interior, 

Bureau of Mines, IC 9247. NTIS No. PB 90-222449. 

Mark C, Chase FE [1997]. Analysis of retreat mining 

pillar stability (ARMPS). In: Mark C, Tuchman RJ, eds. 

Proceedings: New technology for ground control in retreat 

mining. Pittsburgh, PA: U.S. Department of Health and 

Human Services, Centers for Disease Control and Prevention, National Institute for Occupational Safety and Health, 

IC 9446, pp. 17–34. 

Mark C, Chase FE, Molinda GM [1994]. Design of 

longwall gate entry systems using roof classification. In: 

Mark C, Tuchman RJ, Repsher RC, Simon CL, eds. New 

Technology for Longwall Ground Control. Proceedings: 

U.S. Bureau of Mines Technology Transfer Seminar. 

Pittsburgh, PA: U.S. Department of the Interior, Bureau 

of Mines, SP 01–94, pp. 5–17. NTIS No. PB95-188421. 

Molinda GM, Mark C [1994]. Coal mine roof rating 

(CMRR): a practical rock mass classification for coal 

mines. Pittsburgh, PA: U.S. Department of the Interior, 

Bureau of Mines, IC 9387. NTIS No. PB94-160041. 

Pappas DM, Mark C [1993]. Behavior of simulated 

longwall gob material. Pittsburgh, PA: U.S. Department of 

the Interior, Bureau of Mines, RI 9458. NTIS No. PB93­

198034. 

Pariseau WG [2007]. Finite element analysis of interpanel barrier pillar width at the Aberdeen (Tower) mine. 

Department of Mining Engineering, University of Utah, 

and Bureau of Land Management, Salt Lake City, UT. 

Salamon MDG [1966]. Reconsolidation of caved 

areas. Transvaal and Orange Free State Chamber of Mines 

Research Organization, project No. 801/66, research report 

No. 58/66. 

Salamon MDG [1991]. Deformation of stratified rock 

masses: a laminated model. J S Afr Inst Min Metall 91(1): 

9–26. 

Segerlind LJ [1976]. Applied finite element analysis. 

New York: John Wiley & Sons, Inc. 

Sinha KP [1979]. Displacement discontinuity technique for analyzing stress and displacements due to mining 

in seam deposits [Dissertation]. University of Minnesota. 

Skiles ME, Stricklin KG [2009]. General guidelines 

for the use of numerical modeling to evaluate ground control aspects of proposed coal mining plans. Arlington, VA: 

Mine Safety and Health Administration, Program Information Bulletin No. P09-03, March 16, 2009. Available at: 

http://www.msha.gov/regs/complian/PIB/2009/pib09-03.pdf 

St. John CM [1978]. EXPAREA: a computer code for 

analyses of test scale underground excavations for disposal 

of radioactive waste in bedded salt deposits. St. Paul, MN: 

University of Minnesota, Department of Civil and Mineral 

Engineering. DOE contract No. W-7405-ENG-26, report 

No. Y/OWI/SUB-7118/2. 

Turner MJ, Clough RW, Martin HC, Topp LJ [1956]. 

Stiffness and deflection analysis of complex structures. 

J Aeronaut Sci 23(9):805–823, 854. 

Zipf RK Jr. [1992a]. MULSIM/NL: application and 

practitioner’s manual. Pittsburgh, PA: U.S. Department of 

the Interior, Bureau of Mines, IC 9322. NTIS No. PB93­

131993. 

Zipf RK Jr. [1992b]. MULSIM/NL: theoretical and 

programmer’s manual. Pittsburgh, PA: U.S. Department of 

the Interior, Bureau of Mines, IC 9321. NTIS No. PB2004­

105476.