Stochastic modeling of gob gas venthole production performances in active

and completed longwall panels of coal mines

C. Özgen Karacan a,⁎,

a NIOSH, Office of Mine Safety and Health Research, Pittsburgh, PA, United States

⁎ Corresponding author. Tel.: +1 412 386 4008; fax: +1 412 386 6595.

E-mail address: cok6@cdc.gov (C.Ö. Karacan).

Kray Luxbacher b

b Virginia Tech, Dept. of Mining and Minerals Engineering, Blacksburg, VA, United States

Disclaimer: The findings and conclusions in this report are those of the authors and

do not necessarily represent the views of the National Institute for Occupational Safety

and Health. Commercial software mentioned in this paper is not endorsed by NIOSH.

abstract

Gob gas ventholes (GGVs) are an integral part of longwall coal mining operations, enhancing safety by

controlling methane in underground workings. As in many disciplines in earth sciences, uncertainties due to

the heterogeneity of geologic formations exist. These uncertainties, and the wide range of mining and

venthole operation parameters, lead to performance variability in GGVs. Random variations in parameters

affecting GGV performance and influencing parameters that cannot be quantified sufficiently due to lack of

information limit deterministic GGV models and even introduce error in severe cases. Therefore, evaluation of

GGV performance data and the uncertainty in input parameters is valuable for understanding the variability in

GGV production and for designing them accordingly.

This paper describes a practical approach for implementing stochastic determination of GGV production

performances and for generalizing the prediction capability of deterministic models. Deterministic sitespecific models were derived by using the GGV module in the recently developed MCP (Methane Control and

Prediction) software suite. These models were generated using multi-parameter regression techniques and

were then improved by inclusion of extra input parameters that eliminated the site dependency and

improved the predictions. Statistical distributions of input parameters in these models were quantified and

tested with the Kolmogorov–Smirnov goodness-of-fit technique. Next, Monte Carlo simulations were

performed using these distributions and generalized results for GGV performances were generated. The

results of this work indicate that this approach is a promising method of representing the variability in GGV

performances and to improve the limited and site-specific character of the deterministic models.

1. Introduction

Longwall gobs have high-permeability fractures that form open

pathways for gas transport towards the mining environment, making

them prime targets for gas control measures. Drilling vertical gob gas

ventholes (GGVs), which begin to vent gas as the longwall face

advances under their locations, is an effective technique for capturing

methane emissions within the overlying fractured strata before they

enter the work environment. The challenges involved in testing,

characterizing, modeling, and evaluating the reservoir properties of

the gob make it difficult to evaluate the performance and productivity

of GGVs and to design GGVs for the maintenance of underground

health and safety in the face of high methane emissions.

Characterization of the fractured zone and determination of its

reservoir properties are complicated. Because measured data are

scarce and formation deformation during subsidence is complex, it is

difficult to set up a numerical model of the gob for flow simulations.

Previously, Lunarzewski (1998) used boundary element and sequential bed separation methods for floor and roof strata relaxation and

immediate roof bending separation in addition to gas emission rate

calculations. Ren and Edwards (2002) used a computational fluid

dynamics (CFD) modeling approach to investigate gas capture from

surface GGVs. That paper introduced the use of CFD approach to

improve the design of surface gob wells for methane recovery while

minimizing the leakage of air into the gob. Tomita et al. (2003)

developed a three-dimensional (3-D) finite element model (FEM) to

predict the volume of methane gas emitted from surrounding coal and

rock layers based on stress distribution and permeability changes.

Geomechanical modeling is used to simulate the complex permeability, porosity, and elastic property changes that take place with the

fracturing of overburden strata, formation and recompaction of gob,

and desorption of methane from intact reservoirs. Esterhuizen and

Karacan (2007) coupled geomechanical modeling with reservoir

modeling, and Kelsey et al. (2003) coupled geomechanical modeling

with CFD to model drainage through the strata around longwall gob.

CFD is a robust method, but also computationally expensive,
especially for large scale, complex models, such as the longwall gob.

However, it can be used to understand the mechanisms governing

desorption and flow at micro- and macro-levels. Traditional reservoir

simulation, which can be accomplished with several numerical

methods, can also prove computationally expensive, and is primarily

used to understand the influence of reservoir and well parameters,

and forecast long term production. The two methods complement

each other by allowing for characterization of the mechanisms that

induce and influence flow, while also allowing for production

optimization.

In order to address some of the difficulties in modeling in a more

practical way, Karacan (2009a) developed an artificial neural network

(ANN)-based model using field data to predict GGV production rates

and methane concentrations. The total gas productions and methane

percentages of 10 ventholes located over 3 adjacent panels in the

southwest Pennsylvania section of the Northern Appalachian basin

were monitored at the wellhead using pressure, flow, temperature,

and methane sensors. The monitoring took place for more than

2 years, both during and after mining of these panels. The measurements, along with various spatial parameters related to venthole

location, borehole completion parameters, mining rate and panel

completion data, and exhauster operation data, were combined to

form an extensive database. The ANN model was successful in

predicting methane percentage and production rate with a correlation

coefficient, R, greater than 0.9. Sensitivity analyses about the mean of

the input variables were conducted using the ANN model to identify

which input variables had more influence on the performance of

GGVs.

Because methane production from a GGV can often be characterized as an event with multiple episodic phenomena controlled by

numerous variables associated with uncertainty and randomness,

deterministic methods of evaluation and their results may not always

be satisfactory. Stochastic methods, such as Monte Carlo, can support

deterministic methods. Monte Carlo (MC) methods produce an

approximation by calculating results from a large number of random

samples within the test domain.

General Monte Carlo methods have been proven a useful tool for

accurately estimating statistical uncertainties in standard errors and

confidence intervals in non-linear regression problems (Alper and

Gelb, 1990), and are applied frequently in studying particle physics

and in nuclear engineering (Haghighat and Wagner, 2003) to obtain a

probabilistic solution to a deterministic problem. MC methods are also

used in porous media flow and transport problems for ground water

contamination and remediation studies. Huang et al (2003) applied

MC methods to study groundwater flow and solute transport in

heterogeneous, dual-porosity media and compared the results with

analytical models. Morin and Ficarazzo (2006) used stochastic

techniques and MC simulations to predict fragmentation of rock

during blasting. They have shown that the results produced by the

simulator were comparable with the data obtained from a quarry, and

that the use of MC extended the understanding of the factors affecting

blast fragmentation. Sari (2009) demonstrated the use of MC

simulations to evaluate the strength and deformability of rock masses

by including the uncertainties of the intact rock strength and

discontinuity parameters. He concluded that the MC method provided

a viable means for assessing the variability of rock mass properties. Lu

and Zhang (2003) demonstrated the development of an important

sampling method to solve complicated problems with MC and applied

it to fluid transport problems in aquifers.

This paper describes a practical approach for implementing

stochastic determination of GGV production performances and for

generalizing the prediction capability of deterministic models. For this

approach, deterministic models were derived by using the GGV

module in the recently developed MCP (Methane Control and

Prediction) software suite (Karacan, 2010) and generalized to remove

site-specificity. Statistical distributions of input parameters in these

models were quantified and tested for goodness of fit. Monte Carlo

simulations were performed using these distributions and the

generalized deterministic models to assess the variability of GGV

production rates and methane percentages. While CFD and reservoir

simulation, coupled with geomechanical modeling can realistically

simulate the physics of single and multiphase flow they require

considerable knowledge of the reservoir, and significant computational time and expense. The approach presented here provides a

relatively fast estimation of performance for a single point in time and

space.

2. Approach to deterministic modeling of production rates and

methane concentrations

2.1. A brief introduction to Methane Control and Prediction (MCP)

software suite

The MCP software suite was developed to provide proxy solutions

to some of the problems related to controlling and capturing methane

from longwall mines. The software suite consists of dynamic link

library (DLL) extensions to MS-AccessTM, written in C++ and using

various mathematical approaches and artificial neural network (ANN)

methods.

The first version of MCP contains four main modules:

- Coal measure rock mechanical properties prediction. Predicts

dynamic elastic properties of coal-measure rocks for better roof

support and methane control.

- Mine ventilation emission prediction. Predicts ventilation air

methane (VAM) emissions from longwall mines.

- Degasification system selection. Recommends the best degasification choice for a given mine design, operational conditions, and

targeted VAM conditions.

- Gob gas venthole production performance prediction. Predicts

performance of GGVs in terms of total production rate and

methane concentration.

These models produce the desired predictions with known basic

log, mining, longwall panel, productivity, and coal bed characteristics

as input parameters. The application of these modules separately or in

combination for methane capture- and control-related problems can

improve the safety of coal mines. The theory and technical aspects of

developments of these software modules are given in detail in

Karacan (2008, 2009a,b,c), and the software suite is available for

download from http://www.cdc.gov/niosh/mining/products/product180.htm free of charge (Karacan, 2010).

Details of MCP and its possible applications are given in Dougherty

and Karacan (in press). This current paper deals with the deterministic and stochastic approaches to the prediction of GGV performance.

Readers are referred to Karacan (2009a) for the details of development of the GGV module and for a detailed discussion of input

parameters and their influences on production performances.

Fig. 1 shows two screen captures from the input parameter screen

of the GGV module in the MCP software suite. This screen requires

information on the status of the longwall panel, either active or

completed (Fig. 1-A), and on the status of longwall face, as advancing

or idle (Fig. 1-B). Other inputs include those related to borehole

drilling, completion, mine layout, mine operation and operation of an

exhauster at the surface. The module will generate prediction for two

panel and face situations: an advancing panel and active face (A-A) or a

completed panel and idle face (C-I). The A-A and C-I abbrevations will

be used throughout the rest of this paper to designate the two

circumstances that are common to the operation of gob gas ventholes:

when the mine is active (A-A), and after the panel is completed and/or

abandoned (C-I).
A B
 

Fig. 1. GGV module input screens that show common input parameters and the selections for panel status (A) and longwall face status (B).

2.2. The use and application of the GGV production performance

prediction model for deterministic modeling

The GGV production performance prediction module in the current

version of the MCP calculates GGV total gas production rates and the

methane concentrations in the produced gas stream. The various input

parameters are complex, and since the original ANN model was

developed using the field data obtained from the Northern Appalachian

basin (Karacan, 2009a), in some cases the parameters are site-specific.

ANNs are adaptable deterministic systems that deliver one output

for a given set of inputs. They can determine relationships between

different sets of data to solve problems where conventional

deterministic models are inefficient or insufficiently accurate. Such

problems are usually complex and difficult to describe mathematically. Although the key advantages of neural networks are their

abilities to recognize patterns between input and output space and to

generalize solutions, the inherent disadvantage is that ANNs usually

do not offer a deterministic mathematical model that shows the

relations between inputs and outputs explicitly. In that respect, they

can be considered as “black-box” models (Maier and Dandy, 2000).

In this study, MCP's GGV performance prediction model was used

to develop proxy deterministic models for predicting methane

concentration and production rates. The proxy models quantify the

sensitivity of the production rate and methane concentration for

different mining and GGV parameters, and describe the performance

of GGVs with the range of input parameters that were used to develop

the original ANN model. Using these proxy models as the basis for

developing a stochastic approach allows for the statistical uncertainty

that is associated with each input parameter, for the addition of input

parameters, and for generalization of the site-specific nature

associated with some of the input parameters to an extended range

of situations. This approach eventually broadens the applicability of

the stochastic models to more general situations compared to the

deterministic models.

In order to obtain proxy deterministic models, the GGV performance prediction model was used to generate flow rate and methane

concentration output for A-A and C-I situations. Random number

generation between minimum and maximum limits of each input

parameter provided the values for the input parameters. These limits

were based on measured well data. For each of the A-A and C-I cases,

approximately 150 input–output patterns were generated. In this

study, one of the original input parameters in the A-A case, “% of panel

completed,” was replaced with a new parameter called “face past the

borehole location” (Table 1). This new parameter gives a true measure

of how far the face advanced compared to the borehole location and

also ensures that the GGV is operational at that face location. This

parameter is calculated as the difference of distance of the borehole

from panel start minus the product of the percent panel completed

and the panel length.

The input–output patterns were analyzed with multiple linear

regressioninorder to create a reasonably accurate andeasy tomanipulate

deterministic model for use in the stochastic approach. Since the input

values were known within a reasonable range, no elimination procedures

were necessary on the generated data. All of the inputs and their

contributions to the variance in the models were maintained.

Fig. 2-A and B gives the results for methane percentages obtained

from the MCP software plotted against those predicted by the proxy

models for A-A and C-I cases, respectively. As seen from these plots,

the data scatter around the coefficient of determination line, R2= 1,

within one standard deviation (shown by the gray lines). Fig. 2-C and

D shows that standardized residuals obtained by comparing MCP

results with the predictions of methane percentages in A-A and C-I

cases are mostly within +/− 2, which indicate that there is no

significant non-linearity in both cases (Devore and Peck, 2001).

Eqs. (1) and (2) represent the MCP calculations for production rate

and methane concentration in the production stream respectively.

Coefficients of determination (R2

) range from 0.69 to 0.88 and include

all the input parameters of the GGV module in MCP software, as

shown in Table 1. The variables “face past borehole location” and

“linear advance rate” are not included as input parameters in the C-I

case since panel and face are not active in this case.

GGV Production Rate ðscfmÞ = P1 × x1 + P2 × x2 + ::::: + Pn × xn ð1Þ

Methane ð%Þ = R1 × y1 + R2 × y2 + ::::: + Rn × yn ð2Þ

2.3. Improvement of deterministic GGV performance prediction models

(Table 1; Eqs. (1) and (2))

The model given in Karacan (2009a) and the equations that are

discussed in the previous section include most of the major
parameters that might influence GGV production rate and methane

concentration percentage. These parameters were discussed in

Karacan (2009a) in detail and will not be repeated here. It is

reasonable to assume that since the measured data used to develop

the prediction model came specifically from the Northern Appalachian basin, some parameters that affect other basins may not be

accounted for. Although most input parameters can be considered

common to other basins, there are a few that better reflect the

Northern Appalachian basin conditions. Furthermore, a few parameters were intentionally not included in the models, either because

they did not vary or because they were not measured. It is possible

that the intentionally and unintentionally excluded parameters could

have influenced the observed production rates and methane concentrations. These parameters were:

•	 Slotted casing length. The monitored GGVs were drilled and

completed with 200 ft of slotted casing. Therefore, this “variable”

was constant and was not included as an input parameter in the

models.

•	 Distance to the top of slotted casing. This information can be obtained

knowing the overburden depth, distance from bottom of slotted casing

to top of mined coal bed, and the slotted casing length. Since the casing

length was not a variable in the models, the distance to the top of

the slotted casing was not included either.

•	 Strata displacements and subsidence. The general stratigraphy in the

Northern Appalachian basin is almost uniform and has characteristic

formations that are present, although possibly with varying

thicknesses. The effect of varying strata and their properties are

potentially important for subsidence, or strata displacements. In the

absence of direct measurements, subsidence information could not

be included in the models as an input parameter.

•	 Gas content of the overlying formations. The coal beds above the

Pittsburgh seam, believed to be the source of methane produced

from the monitored GGVs in the Northern Appalachian basin, are

usually ranked as high volatile bituminous. These coal beds have

approximately the same amount of total gas content within the

areal scale of the monitored mining field. Since there were no gas

content measurements from the overlying strata for the studied

GGVs, this potential variable was not included in the original

models.

Table 1

Parameters and coefficients of multiple linear regression equations that calculate production rates and methane concentrations from GGVs for A-A and C-I cases.

Production rate (scfm) — Case: A-A Methane concentration (%) — Case: A-A

Parameter Coefficient Parameter Coefficient

Intercept 0.000E+00 Intercept 0.000E+00

Face past BH Face past BH

Location (ft) P1 −2.015E−02 x1 Location (ft R1 3.514E−04 y1

Linear adv. rate (ft/day) P2 7.644E−01 x2 Linear adv. Rate (ft/day) R2 1.261E−01 y2

Surface elevation (ft) P3 −3.062E−01 x3 Surface elevation (ft) R3 −1.505E−02 y3

Average OB depth (ft) P4 8.949E−01 x4 Average OB depth (ft) R4 6.039E−02 y4

Casing diameter (′′) P5 1.666E+02 x5 Casing diameter (′′) R5 −9.419E+00 y5

Slotted casing height from top of coal (ft) P6 −4.631E+00 x6 Slotted casing height from top of coal (ft) R6 2.122E−02 y6

Distance to TG (ft) P7 −1.091E+00 x7 Distance to TG (ft) R7 −2.025E−01 y7

Distance from start (ft) P8 1.837E−02 x8 Distance from start (ft) R8 −1.241E−03 y8

Panel length (ft) P9 −4.366E−02 x9 Panel length (ft) R9 3.491E−02 y9

Panel width (ft) P10 −5.221E−01 x10 Panel width (ft) R10 −1.262E−01 y10

Atm. press (in Hg) P11 1.268E+01 x11 Atm. press (in Hg) R11 −2.319E+00 y11

Ext vacuum (in water) P12 −5.938E−01 x12 Ext vacuum (in water) R12 −3.292E−02 y12

R2=0.81 R2=0.88

Production rate (scfm) — Case: C-I Methane concentration (%) — Case: C-I

Parameter Coefficient Parameter Coefficient

Intercept 0.000E+00 Intercept 0.000E+00

Surface elevation (ft) P1 1.256E−02 x1 Surface −levation (ft) R1 9.736E−03 y1

Average OB depth (ft) P2 1.560E−01 x2 Average OB depth (ft) R2 −8.627E−03 y2

Casing diameter (′′) P3 4.689E+01 x3 Casing diameter (′′) R3 1.192E−01 y3

Slotted casing height from top of coal (ft) P4 3.515E+00 x4 Slotted casing height from top of coal (ft) R4 −1.239E+00 y4

Distance to TG (ft) P5 −2.052E-01 x5 Distance to TG (ft) R5 −4.956E−02 y5

Distance from start (ft) P6 −2.051E-04 x6 Distance from start (ft) R6 9.779E−04 y6

Panel length (ft) P7 −9.452E−03 x7 Panel length (ft) R7 7.016E−03 y7

Panel width (ft) P8 −3.681E−01 x8 Panel width (ft) R8 −2.214E−02 y8

Atm. press (in Hg) P9 1.827E+00 x9 Atm. press (in Hg) R9 2.504E+00 y9

Ext vacuum (in water) P10 −3.384E+00 x10 Ext vacuum (in water) R10 −1.058E−01 y10

R2=0.69 R2=0.74

Because of the missing parameters and the uncertainty carried by

the included parameters, the developed models are also subject to

uncertainty. This makes a stochastic approach an attractive improvement method.

At this juncture, we proceed with the postulation that production

rates and methane concentrations observed from the GGVs are the

results of combinations of parameters, some of which could be

included in the models as inputs and others that could not. Therefore,

the missing inputs are embedded in the coefficients of the existing

ones in Eqs. (1) and (2) (Table 1). If the missing variables can be

included in the models as linear combinations of parameters, the

modified equation can be fitted to the MCP predictions of production

rate and methane concentration again, as in the previous section, to

redetermine the coefficients to all variables in the modified equation.

The following sections describe the potential effects of slotted

casing length (or distance to the top of slotted casing), strata

displacement (subsidence), and gas content data on rate and methane

concentration data, and how these were integrated into the

deterministic equations (Eqs. (1) and (2)). Other input parameters

will not be discussed here, but interested readers are referred to

Karacan (2009a).

2.3.1. Slotted casing length/distance to top of slotted casing and

flow fraction

In gob gas ventholes, the slotted casing is the production interval

along which all gas is expected to flow into the borehole. Thus, vertical

variation of gas-bearing strata associated with the various lengths of
slotted casing is expected to have a major effect on the methane

production. The inclusion of coal beds, gas-bearing sandstones, and

shales in the overlying strata, in the horizon of the slotted casing interval

may make a significant difference in the amount of gas production.

Therefore, the optimum length of slotted casing depends upon the

geological layers in the overlying strata and their gas emission potentials.

Although slotted casing length can significantly affect gas

production there are not sufficient field data to evaluate the effects

of various lengths of slotted casing. A numerical simulation-based

study investigated the changes in production performances of gob gas

ventholes when the slotted casing lengths were changed to 100 ft and

to 250 ft from their initial design length of 200 ft, while keeping other

completion parameters constant (Karacan et al., 2007). The modeling

results predicted that the cumulative methane production would

increase to 459.4 MMscf with a 250-ft casing compared to

391.8 MMscf with a 200-ft casing. This represented a 9.5% increase

in methane capture from the ventholes in the simulated area. When

the slotted casing length was shortened to 100 ft, the predicted

methane production decreased to 314.7 MMscf, an approximately

25% reduction from the 200-ft slotted casing case.

In the present study, gas flow entry locations in GGVs and the

fractions of the total flow as a function of depth were used to evaluate

the effect of casing length on GGV flow performance. A combined flow

profile was constructed by digitizing the data from measurements of

all ventholes given in Mazza and Mlinar (1977). The study was

conducted in the Appalachian Basin using a data logger and a gas gun

to release pulses of Kr-85 tagged nitrogen. The test system recorded

the transit time of a pulse of radioactive gas along a known interval,

while carried in the stream of produced gas. The velocities were then

determined with the known time and distance, from which flow rates

and flow percentages coming from each flow entry point in the

borehole could be calculated.

A C 

100 

MCP - Methane Concentration (%)

80 

60 

40 

20 

0 

Observations 

0 20 40 60 80 100 

Obs1 

Obs10 

Obs19 

Obs28 

Obs37 

Obs46 

Obs55 

Obs64 

Obs73 

Obs82 

Obs91 

Obs100 

Obs109 

Obs118 

Obs127 

-4 -2 0 2 4 

Equation - Methane Concentration (%) Standardized residuals 

B D 

100 

MCP - Methane Concentration (%)

80 

60 

40 

20 

0 

Observations 

0 20 40 60 80 100 -4 -2 0 2 4 

Obs1 

Obs9 

Obs17 

Obs25 

Obs33 

Obs41 

Obs49 

Obs57 

Obs65 

Obs73 

Obs81 

Obs89 

Obs97 

Obs105 

Obs113 

Equation - Methane Concentration (%) Standardized residuals 

Fig. 2. Regression plots of methane concentration expected from the GGVs for A-A (A) and C-I (B) cases, and the standardized residuals (C and D) from these predictions, respectively.

Fig. 3 shows the results of the combined flow profile. Only a minor

percentage of flow enters the ventholes within 50 ft of the mined coal

bed. A significant increase in the flow entry occurs at approximately

175–200 ft above the mined coal seam. Between 50 and 70% of the

total flow originates in this interval. As shown in Fig. 4 and discussed

in the next section, this interval corresponds to the location of major

displacements or strata separations in the overburden. It should be

noted that this interval may not necessarily correspond to the source

of the gas, but rather to the location of the major flow pathways.

Almost 100% of total flow enters the GGVs by about 300–350 ft from

the top of the mined coal bed.

The parameter “distance to the top of casing” can be used to

include the flow entry profile in modified deterministic equations

(Eqs. (1) and (2)). A new mathematical term “distance to the top of

slotted casing (D)” was defined. Distance to top of slotted casing (D) is

the average overburden depth (OB) minus the sum of slotted casing

length (SL) and the distance from the bottom of the slotted casing to

the top of the coal seam. This factor that modifies flow rates obtained

with 200-ft slotted casing for other casing lengths. An exponential

equation, y=41.991e2.0071×, was fitted to the data, as shown in Fig. 3.

Flow percents (FP) that can be obtained with casing lengths other

than 200 ft and the ones with 200 ft (for the same overburden and

distance of the casing bottom to top of coal bed) can be defined and

the flow-rate-change factor (FRF) can be calculated as:

=

0

ln FRF (

OB−D−SL)1,

41:991 FP j 200−ft ð3Þ

2:0071 @ A

This equation gives a flow rate scaling factor relative to 200-ft

casing length. The proposed approximation is applicable only if this
factor can be calculated for given values of SL, D, and OB, including

their uncertainities.

y = 41.991e2.0071x 

R2 = 0.8043 

400 

350 

300 

Distance from the top of the coal seam (ft) 250 

200 

150 

100 

50 

0 

0 0.5 1 1.5 

Fraction of flow 

Fig. 3. Combined flow profiles obtained from GGVs drilled over various longwall mining

operations and the exponential equation fitted to this data.

400

350

300

y = 350e-0.698x 

R2 = 0.402 

Distance from the top of the coal seam (ft) 250 

200 

150 

100

50 

0 

0 1 2 3 4 

Displacement (ft) 

Fig. 4. Combined strata separations obtained from GGVs drilled over various longwall

mining operations in the Northern Appalachian basin and the exponential equation

fitted to this data.

2.3.2. Strata displacements and subsidence

During longwall mining, vertical and horizontal fractures occur

during the movement of the overburden. Horizontal fractures in the

formation mostly occur along strong–weak rock layer interfaces and

influence the hydraulic conductivity of the overburden strata, which

creates methane emission pathways and determines methane

emissions into the mine (Karacan and Goodman, 2009). The ability

to estimate the locations and magnitudes of the fractures and their

hydraulic flow properties is important for placement of boreholes and,

consequently, for controlling methane more effectively.

Palchik (2003) has shown that the extent of the fractured zone

induced by mining in the Donetsk coal basin can be determined based

on the change in natural methane emissions from this zone. After

modifying the test system, he was able to locate the individual

fractures at bedding-plane separations and to determine their

apertures (Palchik, 2005, 2010). Further, the presence and absence

of estimated horizontal fractures was correlated with uniaxial

compressive strength and thickness of rock layers, distances from

the extracted coal seam to the rock layer interfaces, and the

thicknesses of extracted coal seams.

In addition, observations on the presence and absence of

horizontal fractures at different rock layer interfaces of the overburden showed that the probability of separation increased with

increasing compressive strength difference of neighboring rock layers,

and with decreasing distance of the layer from the mined coal bed

(Palchik, 2005; Karacan and Goodman, 2009). Therefore, beddingplane separations and resultant strata displacements are functions of

both distance from the coal seam and the type of the overlying

formations. The data suggest that subsidence observed at the surface

can be considered as a lumped parameter that results from the

interactions of these factors. In the absence of specific data regarding

the type and strength of overlying formations and their respective

depths, the average value of strata displacements or average

subsidence can be used as a lumped parameter with inherent

uncertainty.

In this study, combined displacement profiles were constructed,

and the magnitudes of strata separations in the profiles were used to

generate a function for the likelihoods of strata separations and

depths. The profiles were constructed from measurements of all

ventholes reported in Mazza and Mlinar (1977). The data were

measured using Cobalt-60 tagged shaped charges as formation

markers and gamma-ray logging. This system was used in the

ventholes to implant radioactive markers into the formations to

measure the subsidence and the strata separation during mining, and

the data are displayed in Fig. 4.

Fig. 4 shows that vertical displacements due to strata separations

can occur as much as 300–350 ft above the mined coal seam (in this

case, the Pittsburgh coal seam). The magnitudes of strata separations

at this distance to almost 200–250 ft above of the Pittsburgh seam are

generally less than 1 ft, and mostly around 0.5 ft. The displacement

trend in this interval is almost linear and not a function of depth. At

about 200 ft above the Pittsburgh coal bed, the separations start to

depart from low displacement values, and major strata separations on

the order of 2 ft can be measured. These high separations are mostly

100 ft above the coal seam, where the Sewickley sandstone occurs.

This observation and the displacement measurements are in line with
the flow entry measurements discussed in the previous section. The

displacement trend at and below 200 ft is nearly linear with an

approximate 45-degree relationship with depth.

In order to establish a correlation between strata displacement

(SD) and depth, an exponential function was fitted to field measurements (Fig. 4). This function was later used as a surrogate to

approximate surface subsidence, in accordance with treating average

subsidence as a lumped parameter. Representative average strata

displacement, used as a surrogate for subsidence, was calculated as

the sum of the displacements at depth intervals where the bottom and

top of the slotted casing were located with respect to the top of coal

bed. This relation is given as:

0 ( ) D ln 1

350 SD = − 0:698

1 0

ln(

OB−D2

1

B

B CC B 350 + B−

)

C

@ A @ C ð4Þ 0:698 A

where D1 is the distance from bottom of slotted casing to top of coal

bed, D2 distance from surface to the top of slotted casing, and OB is the

overburden depth. Vertical and horizontal anisotropies are not

considered independently, but are simply part of the lumped

parameter.

2.3.3. Gas content of overlying formations (coal beds)

Determining the source of longwall methane and of the gas

produced from gob gas ventholes is important since the gas source

may dictate the amount and composition of the gas. Diamond et al.

(1992) conducted a series of field tests in the Lower Kittanning coal

bed. Boreholes were drilled before and after mining a panel, in order

to obtain coal and rock samples from overlying strata, and to

determine their gas contents. The modified direct method (MDM)

was employed to determine the gas contents of recovered coal and

rock samples (Diamond and Schatzel, 1998). The results indicated

that approximately 90% of the gas removed from the overlying strata

came from coal beds. Material balance calculations were made to

compare the volume of gas produced from GGVs drilled over the panel

and gas removed by the mine's ventilation system with the volume of

gas removed from strata directly overlying the panel to a height of

275 ft.

Since most longwall gas comes from overlying coal beds, the gas

contents of the coals in the Pennsylvania (PA) section of the Northern

Appalachian basin were compiled from Diamond et al. (1986) and

plotted versus depth as shown in Fig. 5.

Total Gas Content (scf/ton)
 

Pennsylvania Coals
 

0 200 400 

Depth (ft) 

0 

200 

400 

600 

800 

1000 

1200 

1400 

1600 

y = 4.9631x 

R2 = 0.568 

Fig. 5. Total gas contents of high volatile bituminous (HV) coal beds of Pennsylvania and

their variations with overburden depth.

This plot established a relation

between gas content and depth, and was incorporated into the

deterministic equations as an input parameter affecting gas rate and

methane concentration from the GGVs. Here, the relation was

represented with a linear trend of gas content increase with

increasing depth.

Almost all the coals that are reported in Diamond et al. (1986) for

PA are high volatile bitumious A, with lesser amounts of B and C, but

this is not true in all coal basins. In order to include the gas contents of

coals of all ranks into the predictive equations, all HV (A, B and C),

medium- and low-rank coals, and sub-bituminous (A, B and C) coals

from all basins that are reported in Diamond et al. (1986) were

compiled and grouped into 3 main classes of rank and gas content

(HV, MV–LV, and sub-bituminous), as distributions. This compilation

is found in Section 3, Table 4.

2.3.4. Modified deterministic equations for gas production rate and

methane concentration prediction from gob gas ventholes

In order to modify the original deterministic equations and their

parameter coefficients, the parameters discussed in the previous

sections were added to original equations as linear combinations.

Since distance to the top of the slotted casing and casing length are

related and create multi-colinearity, only distance to the top of the

slotted casing was included in the new equations. Casing length can

easily be calculated with knowledge of overburden depth, distance to

top of slotted casing, and the distance from bottom of slotted casing to

top of coal bed.

The modified equations, shown as Eqs. (5) and (6), representing

well production rates and methane concentration, respectively, were

fitted to the same production flow rate and methane concentration

data obtained earlier from the MCP software, for both A-A and C-I

cases.

GGV Production Rate ðscfmÞ = P1 × x1 + P2 × x2 + ::::: + Pn × xn ð5Þ

Methane ð%Þ = R1 × y1 + R2 × y2 + ::::: + Rn × yn ð6Þ

These equations, when fitted to the gas production rate and

methane concentration data, resulted in regression coefficients (R2

) of

0.83 and 0.69 for prediction of production rate from A-A and C-I cases,

respectively. For prediction of methane concentrations, the regression

coefficients were 0.89 and 0.74 for A-A and C-I cases, respectively.

These values are very close to the ones obtained by using Eqs. (1) and

(2), shown in Table 1.

Table 2 gives the new coefficients of input parameters for Eqs. (5)

and (6), with the extra parameters from the modified equations in

bold. Analysis of parameter coefficients shows neither distance to top

of slotted casing (and thus slotted casing length) nor subsidence

produce statistically significant effects on methane concentration in

either A-A or C-I cases. Gas content of the coal bed had no statistically

significant effects on production flow rate of GGVs.
The results obtained using the original deterministic equations

(Table 1) and the modified equations (Table 2) were also compared

using values of descriptive statistics obtained from the predictions of

production rate and methane concentration. These data are given in

Table 3. Almost all descriptive values are very similar for predictions

of rates and methane percentages. In both A-A and C-I cases, absolute

value of the skewness decreased, indicating a more normal distribution of results with the modified models.

Table 2

Parameters and coefficients of the modified multilinear regression equations that were developed for calculating production rates and methane concentrations from GGVs for A-A

and C-I cases.

Production rate (scfm) — Case: A-A Methane concentration (%) — Case: A-A

Parameter Coefficient Parameter Coefficient

Intercept 0.000E+00 Intercept 0.000E+00

Distance to top of slotted casing (ft) P1 1.277E+00 x1 Distance to top of slotted casing (ft) R1 0.000E+00 y1

Face past BH Face past BH

Location (ft) P2 −2.051E−02 x2 Location (ft R2 3.334E−04 y2

Linear adv. rate (ft/day) P3 4.926E−01 x3 Linear adv. rate (ft/day) R3 1.254E−01 y3

Surface elevation (ft) P4 −2.301E−01 x4 Surface elevation (ft) R4 −1.468E−02 y4

Average OB depth (ft) P5 −1.084E+00 x5 Average OB depth (ft) R5 5.416E−02 y5

Casing diameter (′′) P6 1.549E+02 x6 Casing diameter (′′) R6 −9.379E+00 y6

Slotted Casing height from top of coal (ft) P7 −1.112E+01 x7 Slotted Casing height from top of coal (ft) R7 4.509E−02 y7

Distance to TG (ft) P8 −9.449E−01 x8 Distance to TG (ft) R8 −2.068E−01 y8

Distance from start (ft) P9 5.538E−03 x9 Distance from start (ft) R9 −1.226E−03 y9

Panel length (ft) P10 2.567E−02 x10 Panel Length (ft) R10 3.484E−02 y10

Panel width (ft) P11 3.579E−01 x11 Panel width (ft) R11 −1.276E−01 y11

Atm. press (in Hg) P12 1.574E+01 x12 Atm. press (in Hg) R12 −2.307E+00 y12

Ext vacuum (in water) P13 −1.611E+00 x13 Ext vacuum (in water) R13 −3.847E−02 y13

Coal gas content (scf/ton) P14 0.000E+00 x14 Coal gas content (scf/ton) R14 3.526E−02 y14

Subsidence (ft) P15 −2.493E+02 x15 Subsidence (ft) R15 0.000E+00 y15

R2=0.83 R2=0.89

Production rate (scfm) — Case: C-I Methane concentration (%) — Case: C-I

Parameter Coefficient Parameter Coefficient

Intercept 0.000E+00 Intercept 0.000E+00

Distance to top of slotted casing (ft) P1 6.463E+00 x1 Distance to top of slotted casing (ft) R1 0.000E+00 y1

Surface elevation (ft) P2 2.734E−02 x2 Surface elevation (ft) R2 9.734E−03 y2

Average OB depth (ft) P3 −6.222E+00 x3 Average OB depth (ft) R3 −1.194E−02 y3

Casing diameter (′′) P4 5.477E+01 x4 Casing diameter (′′) R4 1.642E−01 y4

Slotted casing height from top of Coal (ft) P5 1.594E+01 x5 Slotted casing height from top of Coal (ft) R5 −1.242E+00 y5

Distance to TG (ft) P6 −2.245E−01 x6 Distance to TG (ft) R6 −4.917E−02 y6

Distance from start (ft) P7 −3.187E−03 x7 Distance from start (ft) R7 9.722E−04 y7

Panel length (ft) P8 4.136E−03 x8 Panel length (ft) R8 6.985E−03 y8

Panel width (ft) P9 −2.769E−01 x9 Panel width (ft) R9 −2.234E−02 y9

Atm. press (in Hg) P10 7.302E+00 x10 Atm. press (in Hg) R10 2.513E+00 y10

Ext vacuum (in water) P11 −3.670E+00 x11 Ext vacuum (in water) R11 −1.066E−01 y11

Coal gas content (scf/ton) P12 0.000E+00 x12 Coal gas content (scf/ton) R12 1.430E−02 y12

Subsidence (ft) P13 1.380E+02 x13 Subsidence (ft) R13 0.000E+00 y13

R2=0.69 R2=0.74

Table 3

Descriptive statistics for the rate and methane concentration outputs obtained by two

sets of equations.

A-A case (active panel; Original multilinear eq. Modified multilinear eq.

advancing face) (Eqs. (1) and (2)) (Eqs. (5) and (6))

Rate (mean) 224.91 210.94

Rate (std. dev.) 189.98 157.42

Rate (median) 224.72 209.87

Rate (variance) 36094.22 24783.70

Rate (skewness) 0.140 -0.010

Methane % (mean) 57.04 56.18

Methane % (std. dev.) 27.29 27.29

Methane % (median) 57.29 56.52

Methane % (variance) 745.25 744.96

Methane % (skewness) -0.021 -0.021

C-I case (completed panel; Original multilinear eq. Modified multilinear eq.

idle face) (Eqs. (1) and (2)) (Eqs. (5) and (6))

Rate (mean) 183.19 187.61

Rate (std. dev.) 78.24 81.81

Rate (median) 186.02 189.38

Rate (variance) 6122.67 6694.48

Rate (skewness) -0.039 0.00007

Methane % (mean) 57.25 56.91

Methane % (std. dev.) 11.92 11.94

Methane % (median) 57.08 56.73

Methane % (variance) 142.23 142.66

Methane % (skewness) 0.017 0.016

3. Stochastic method for evaluating production rates and methane

concentrations using distributions of input parameters

3.1. Generating probability density functions (PDF) for distribution of

input values

To move from a deterministic model to a stochastic model,

probability distributions over the range of values for each parameter

must be determined. Monte Carlo simulation generates random values

for input parameters in a model whose parameters cannot be defined by

exact functions, but by probability density functions (PDFs). Describing

the population of any variable using probability density functions is

helpful when directintegration ofthe function ofthat variable in a model

is not practical or when the exact behavior of that variable is not known.

The PDFs that best describe these populations can then be used to

generate random values for each variable to be used in deterministic

models. Once the likely probability density function has been determined, a random number is fed into the inverse equation to determine

the value to be generated for the distribution of any input parameter. The
random number is generated from a uniform distribution between 0 and

1, so that it falls in an unbiased percentile range. The inversion is

employed in the sampling method stage (Ripley, 1987; Vose, 2000).

The output of a deterministic model that is solved this way is

reasonably accurate only if the number of simulations is very large.

The advantage of this method is that the complete probability

distribution of the model output can be obtained (Sari, 2009).

In this section, the distributions of input data used in the MCP

software for production rate and methane concentration for A-A and

C-I cases are used to make generalizations about the populations of

each of the input parameters. These generalizations are then used in

Monte Carlo simulations of the modified deterministic equations

(Eqs. (5) and (6)) for estimating and generalizing production rate and

methane percentages in A-A and C-I cases.

A 

0.0004 

0.00035 

0.0003 

0.00025 

B 

0.0007 

0.0006 

0.0005 

Density

0 2000 4000 6000 8000 10000 12000 0 2000 4000 6000 8000 10000 

Density

0.0004 

0.0003 

0.0002 

0.00015 

0.0002 0.0001 

0.00005 0.0001
 

0
 0 

GGV Distance from Panel Start (ft) GGV Distance from Panel Start (ft) 

Distance from Start (ft) Log-normal(7.415,0.806) Distance from Start (ft) Log-normal(7.061,0.714) 

C D
 

0.6 0.8 

0.7

0.5 

0.6 

0.4 0.5 

5 6 7 8 9 10 5 6 7 8 9 

Density

Density 0.4 

0.3 

0.3 

0.2 

0.2 

0.1 0.1
 

0
 0 

Casing Diameter (inch) Casing Diameter (inch) 

Casing Diameter (inch) Normal(7.416,0.705) Casing Diameter (inch) Normal(7.178,0.567) 

E F 

0.025 0.03 

0 

0.005 

0.01 

0.015 

0.02 

-200 -150 -100 -50 0 50 

0 

0.005 

0.01 

0.015 

0.02 

0.025 

-100 -50 0 

Density

Density 

Extraction Vacuum (in Water) Extraction Vacuum (in Water) 

Extraction Vacuum (in Water) Normal(-47.975,19.735) Extraction Vacuum (in Water) Normal(-38.519,15.747) 

Fig. 6. Example PDFs for the values used in running MCP for rate and methane percentage estimations from GGVs in A-A (A, C, and E) and C-I (B, D, and F) cases.

Fig. 6A–F shows example input variables “GGV distance from panel

start”, “casing diameter”, and “extraction vacuum” and their distributions for the A-A and C-I cases, used in generating rate and methane

percentage outputs from MCP. As mentioned in Section 2.3.3, gas
contents of coals of HV (high volatile), MV (medium volatile) and LV

(low volatile), and sub-bitiminous ranks were included in the

predictive equations of the A-A and C-I cases as 3 main classes

(Table 4).

Table 4

Generalized distributions, and their means and standard deviations, obtained for each of the input variables and the results of K–S statistics for hypothesis testing for goodness-of-fit.

A-A case: (active panel; advancing face)

Parameter Distribution (μ, σ) Kolmogorov–Smirnov test (D; p)

Distance to SC top (ft)

Face past BH location (ft)

Linear adv. rate (ft/day)

Surface elevation (ft)

Average OB depth (ft)

Casing diameter (′′)

Slotted casing height from top of coal (ft)

Distance to TG (ft)

Distance from start (ft)

Panel length (ft)

Panel width (ft)

Atm. press (in Hg)

Ext vacuum (in water)

Subsidence (ft)

N (502.982; 85.069)

N (3182.644; 1783.761)

N (30.123; 14.250)

N (1252.033; 111.941)

N (749.435; 67.908)

N (7.416; 0.705)

N (43.395; 6.34)

N (277.419; 27.191)

LN (7.415; 0.806)

N (10702.566; 548.847)

N (1282.733; 62.508)

N (26.071; 2.922)

N (−49.975; 19.735)

N (3.585; 0.395)

(0.018; 0.905)

(0.039; 0.103)

(0.022; 0.724)

(0.025; 0.596)

(0.019; 0.884)

(0.026; 0.531)

(0.013; 0.996)

(0.017; 0.942)

(0.032; 0.283)

(0.018; 0.921)

(0.024; 0.629)

(0.025; 0.610)

(0.024; 0.634)

(0.047; 0.331)

C-I case: (completed panel; idle face)

Parameter Distribution (μ, σ) Kolmogorov–Smirnov test (D; p)

Distance to SC top (ft)

Surface elevation (ft)

Average OB depth (ft)

Casing diameter (′′)

Slotted casing height from top of coal (ft)

Distance to TG (ft)

Distance from start (ft)

Panel length (ft)

Panel width (ft)

Atm. press (in Hg)

Ext vacuum (in water)

Subsidence (ft)

N (522.514; 68.786)

N (1176.299; 99.239)

N (741.300; 59.929)

N (7.178; 0.567)

N (44.522; 7.242)

N (266.941; 41.001)

LN (7.061; 0.714)

N (10350.710; 692.338)

N (1249.957; 103.286)

N (27.376; 1.776)

N (−38.519; 15.747)

N (3.715; 0.392)

(0.023; 0.690)

(0.026; 0.486)

(0.022; 0.706)

(0.020; 0.831)

(0.015; 0.979)

(0.024; 0.605)

(0.023; 0.687)

(0.023; 0.660)

(0.028; 0.423)

(0.023; 0.681)

(0.017; 0.944)

(0.037; 0.127)

Coal rank–gas content (scf/ton) Distribution (μ, σ) Kolmogorov–Smirnov test (D; p)

High vol. bit. (A, B, and C) coals

Low–medium volatile coals

Sub-bit. (A, B, and C) coals

High vol. bit. (A, B, and C) coals

Low–medium volatile coals

Sub-bit. (A, B, and C) coals

A-A

C-I

LN (4.270; 1.078)

N (356.957; 136.895)

LN (1.906; 0.918)

LN (4.270; 1.083)

N (364.496; 134.395)

LN (1.794; 0.853)

(0.043; 0.630)

(0.041; 0.660)

(0.038; 0.703)

(0.053; 0.314)

(0.024; 0.989)

(0.026; 0.975)

In order to generalize the populations of input values, they were

compared against different distribution functions, as shown in Fig. 6

for the example inputs, and as given in Table 7 for all input variables of

the modified equations. Densities of the values for each of the inputs

were statistically compared to the hypothetical distributions to

ensure that the “assumption” of a selected distribution would

represent the data distribution. This statistical procedure is called

the goodness-of-fit test.

For the goodness-of-fit test described above, the Kolmogorov–

Smirnov (K–S) test (Chakravarti et al., 1967) was used. The K–S test is

based on:

i−1 i D = max F Yi − ; −F Yi 7 1≤i≤N

(

ð Þ N N ð Þ)

ð Þ

In this equation, Yi are the data pairs ordered from smallest to

largest value (Y1, Y2... YN). N is the number of ordered data pairs. F is

the theoretical cumulative distribution of the sample distribution

being tested and it must be fully specified with mean and standard

deviation. The test starts with a hypothesis that the data follow a

specified distribution (H0) and a signi cance fi level (α). The hypothesis

regarding the proposed distribution is rejected if the test statistic, D, is

greater than the critical value, p, obtained for the corresponding

significance level, and if the critical value p is less than the significance

level. In that case, the use of an alternative distribution is suggested

(Ha). If a hypothesis is accepted, the value of p is the risk of rejecting

that hypothesis although it is true.

In this work, 0.01 was used as significance level in all the tests and

for all variables, from which the K–S statistics and critical values given

in Table 4 were obtained. Table 4 shows that distributions of most

values can be represented either by normal (N) or log-normal (LN)

distributions with specified means (μ) and standard deviations (σ).

The density functions for normal and log-normal distributions are

given in Eqs. (8) and (9), respectively.

ffiffiffiffiffiffi

ffiffiffiffiffiffi

1 2 x−μ

f xð Þ = p e

ð Þ

2 2σ ; σ N 0 ð8 σ 2π Þ

1 2 ðlnð Þx −μ

f x Þð Þ = e 2 2σ ; x; σ N 0 9 xσ

p

2π ð Þ

3.2. Monte Carlo (MC) simulations

RiskAMPTM (Structured Data, 2008) was used as the Monte Carlo

simulator in this work. This program provides utilization of Monte

Carlo functions and random distributions as MS-ExcelTM add-in

features and creates spreadsheet models that employ Monte Carlo

simulations.

In this study, 5000 simulations for each of the coal ranks (HV, MV–

LV and sub-bituminous) in each of the longwall cases (A-A and C-I) in

Table 4 were run using the Latin Hypercube sampling technique. Latin
Hypercube sampling is a sampling method usually utilized for close

representation of the parameter probability distributions and was first

described by McKay et al (1979) as an appropriate method for

selecting values of input parameters. These cases were run using the

generalized distributions given for each input in Table 4 in modified

linear regression equations, Eqs. (5) and (6), as predictors for the GGV

performance parameters.

It should be noted that Monte Carlo simulation cannot easily

include the covariance between input variables, resulting in shortcomings as far as joint uncertainty of input variables is concerned. This

problem is usually solved with Bayesian Markov Chain Monte Carlo

methods (Ades and Lu, 2003). However, this discussion is beyond the

scope of this paper.

3.2.1. Conditional arguments to MC simulations for refinement of

parameter ranges

Problems with unrealistic input values or unacceptable calculated

outputs can arise in random simulations. This usually occurs when

sampling values from either side of the distributions, and is due to the

departure of representative distribution functions from the actual

data. In order to correct this problem, a set of conditional arguments

can be defined.

Five thousand MC simulations were executed for each of the three

coal ranks, and each of the two longwall cases to produce 30,000

(6 × 5000) different possible combinations of production rate and

methane percentage outputs, each of which corresponded to a certain

combination of random inputs. Conditional arguments were used to

correct for unrealistic inputs and outputs.

The generated input data was run through a series of IF THEN

routines to control the validity of the input values. Input values that

violated the conditions would result in the entire row of inputs being

assigned null values. For example, any distance, depth, or gas content

value must be positive. So, if the value generated for the “distance to

the top of slotted casing” was negative, then those sets of inputs were

assigned null values. Similarly, the routines checked for impossible

combinations of values. For example, if the sum of the values

generated for “face past borehole location” and “distance of borehole

from panel start” was greater than the value generated for “panel

length,” then those inputs were assigned null values.

Another set of IF THEN statements was implemented for the

generated output values. Production rates and methane percentages

cannot be negative, and methane percentage cannot be lower than

10% (as an operational constraint) and cannot exceed 100%. Any

combinations of input values that gave outputs that violated these

conditions were also assigned null values and eventually removed

from input–output data.

Imposing conditions on inputs and outputs during MC simulations

results in defining the ranges of input values that gave reasonable

outputs. Once the populations and ranges of these input values were

defined, minimum, mean and maximum values for the ranges could

be determined. These values are given in Table 5.

4. Results: production rate and methane

concentration distributions

Conditional arguments imposed during MC simulations resulted in

ranges of input values that could be applicable to modeling GGV

productions (Table 5) and methane concentrations. Probability

density functions of rate and methane concentration distributions

computed using MC simulations within the range of values given in

Table 5 are given in Figs. 7 and 8 for A-A and C-I cases, respectively.

Figs. 7 and 8 also show that rate and methane concentration outputs

can be represented with normal distributions with varying mean and

standard deviation values depending on the rank of coals in overlying

formations.

Table 5

Value ranges that are generated as a result of conditional MC to use in modified

equations.

A-A case: (active panel; advancing face)

Parameter Min; mean; max

Distance to SC top (ft)

Face past BH location (ft)

Linear adv. rate (ft/day)

Surface elevation (ft)

Average OB depth (ft)

Casing diameter (′′)

Slotted casing height from top of coal (ft)

Distance to TG (ft)

Distance from start (ft)

Panel length (ft)

Panel width (ft)

Atm. press (in Hg)

Ext vacuum (in water)

Overlying coal gas content HV bit.

(scf/ton) LV–MV bit.

Sub-bit.

Subsidence (ft)

(201.9; 502.9; 790.1)

(9.4; 3182.6; 9128.3)

(0.1; 30.1; 107.4)

(918.3; 1232.0; 1573.8)

510.1; 743.4; 988.2)

(5.3; 7.4; 9.7)

(24.0; 43.4; 64.1)

(204.3; 277.4; 379.2)

(144.9; 2230.6; 10138.5)

(9010.6; 10702.6; 12215.5)

(1099.2; 1282.7; 1448.5)

(16.6; 26.1; 34.7)

(−110.5; −47.9; 17.8)

(1.8; 118.5; 699.6)

(29.4; 356.9; 715.7)

(0.7; 10.4; 135.9)

(2.6; 3.6; 5.1)

C-I case: (completed panel; idle face)

Parameter Min; mean; max

Distance to SC top (ft)

Surface elevation (ft)

Average OB depth (ft)

Casing diameter (′′)

Slotted casing height from top of coal (ft)

Distance to TG (ft)

Distance from start (ft)

Panel length (ft)

Panel width (ft)

Atm. press (in Hg)

Ext vacuum (in Water)

Overlying coal gas content HV bit.

(scf/ton) LV–MV bit.

Sub-bit.

Subsidence (ft)

(254.7; 522.5; 739.4)

(884.5; 1176.3; 1527.9)

(515.4; 741.3; 937.6)

(5.4; 7.2; 8.9)

(20.1; 44.5; 65.8)

(123.4; 266.9; 399.4)

(153.6; 1503.7; 9526.7)

(8301.3; 10350.7; 12830.5)

(926.5; 1249.9; 1599.6)

(21.1; 27.4; 33.5)

(-99.1; -38.5; 22.2)

(2.2; 115.2; 676.4)

(13.5; 364.5; 707.7)

(0.8; 8.7; 69.7)

(2.7; 3.7; 5.3)

Fig. 7 shows that, for a GGV operating in an A-A situation with the

overlying formation containing HV coals, the most probable GGV

production rate is 200–250 scfm (Fig. 7-A), although the rate can be as

high as 650 scfm. With different combinations of input values, the

most likely methane concentration is 55–60%. Concentrations can be

as high as 100% in the produced gas (Fig. 7-D).

Fig. 7-B, E, C and F shows the resultant distributions of rate and

methane concentration with medium–low volatile coals and subbituminous coals in the overburden, respectively. These distributions

show that medium-low volatile coals in the overburden most likely

produce 240 scfm (B) at a concentration of 40–80% (E). Subbituminous coals in the overburden, on the other hand, can result in

a most probable rate of 230 scfm (C) and a most probable methane

concentration around 50% (F).

Fig. 8 shows the possible values of production rate and methane

percentage that might be expected from aGGV operating in a completed

longwall panel. The figure shows that higher rates, in general, can be

achieved with different combinations of input values compared to the

A-A shown in Fig. 7. Since methane concentrations generally tend

towards lower values and the maximum methane concentrations are at

most 90%, this data may suggest that higher rates are due to drawing

more air, instead of methane, into the GGV for a given combination of

location, drilling, and operation parameters.

Table 6 shows the descriptive statistical values calculated using the

distributions of rates and methane concentrations obtained from MC

simulations. The values for the A-A case show that mean production

rates and methane concentrations are slightly higher in the case of

medium-low volatile coals present in overlying gas formations.
A D
 

0.0035 0.025 

0.003 

0.02 

0.0025 

Density 

0 200 400 600 0 20 40 60 8 100 120 

0.015 

Density

0.01

0.002 

0.0015 

0.001 

0.0005 

0.005 

0 0 

Production rate (scfm) Methane (%) 

Production rate (scfm) Normal(235.951,131.969) Methane (%) Normal(55.108,20.947) 

B E
 

0.0035 0.02 

0.003 

0.015 0.0025 

0.0005 

0 0 

0 20 40 60 80 100 120 

Production rate (scfm) Methane (%) Density 

0 200 400 600 

Density

0.002 

0.0015 

0.001 

0.01 

0.005 

Production rate (scfm) Normal(239.002,129.852) Methane (%) Normal(59.240,21.938) 

C F
 

0.004 0.025 

0.0035 

0.02 

0.003 

0.0025 

Density 

0 200 400 600 0 20 40 60 80 100 120 

0.015

Density

0.002 

0.0015 

0.001 

0.01 

0.005 

0.0005 

0 0 

Production rate (scfm) Methane (%) 

Production rate (scfm) Normal(231.082,123.965) Methane (%) Normal(53.245,20.519) 

Fig. 7. PDFs of rate (A, B, and C) and methane percentage (D, E, and F) outputs obtained from conditional MC simulations for A-A case, for HV coals (A and D), MV–LV coals (B and E),

and sub-bituminous coals (C and F), and the normal distributions fitted to those.

Predicted maximum rates are around 640 scfm with methane

concentrations as high as 100% in all cases. The C-I values show that

rates in that case have mean values of 330–340 scfm, with extreme

values of 750–800 sfcm. Maximum methane concentrations for the

C-I case do not exceed 90%.

A more efficient way of interpreting the results of an MC

simulation is to analyze the calculated percentiles. Fig. 9 shows the

percentiles corresponding to the different production rate (A) and

methane percentage outputs (B) from conditional MC simulations.

These results show that it is more probable for a GGV to produce at
higher rates in a C-I situation than in an A-A situation. However, the

probability of producing higher methane concentrations is greater in

an A-A situation. Further, in the C-I case, a large range of GGV

production percentiles (20th–70th) is windowed in a narrow

methane concentration (50–60%) range, indicating a relatively high

probability of producing gas within this methane percentage range.

This suggests that changes in input parameter values have a relatively

small effect on methane concentration for C-I cases.

A D 

0.0025 0.035 

0.03 

0.002 

0.025 

Density 

0 200 400 600 800 0 20 40 60 80 

0.0015

Density

0.02 

0.015 0.001 

0.01 

0.0005 

0.005
 

0
 0 

Production rate (scfm) Methane (%) 

Production rate (scfm) Normal(327.005,182.425) Methane (%) Normal(56.790,12.394) 

B E 

0.0025 0.04 

0.035 

0.002 

0.03 

0.025

Density 

0 200 400 600 800 20 40 60 80 100 

0.0015

Density

0.02 

0.015 

0.01 

0.0005 

0.005
 

0
 0 

Production rate (scfm) Methane (%) 

0.001 

Production rate (scfm) Normal(337.767,189.055) Methane (%) Normal(58.401,11.541) 

C F
 

0.0025 0.035

0.03 

0.002 

0.025 

Density 

0 200 400 600 800 0 20 40 60 80 100 

0.0015 

0.001

Density

0.02 

0.015 

0.01 

0.0005 

0.005
 

0
 0 

Production rate (scfm) Methane (%) 

Production rate (scfm) Normal(339.106,187.622) Methane (%) Normal(54.148,12.235) 

Fig. 8. PDFs of rate (A, B, and C) and methane percentage (D, E, and F) outputs obtained from conditional MC simulations for C-I case, for HV coals (A and D), MV–LV coals (B and E),

and sub-bituminous coals (C and F). Also shown are the fitted normal distributions.

The results shown in Fig. 9 are presented in Table 7 in tabulated

form. In this table, the rates and methane concentrations
corresponding to 20th, 50th and 80th percentiles are shown in red to

highlight values which exhibit low-, mean- and high-possibilities for

the total number of iterations. For the A-A case, 80% of total iterations

resulted in production rates of 350 scfm or lower and methane

concentrations of 73% or lower when there were HV coals in the

fractured overburden. The numbers for the same percentile for the C-I

case and HV coals in the overlying strata as the main methane source

are 504 scfm or less and 68% or less. Similarly, analyses can be

conducted for other coal ranks and at different percentiles.

Table 6

Summary of descriptive statistics performed on the MC results.

A-A case: (active panel; advancing face)

HV coals MV-LV coals Sub-bit. coals

Parameter Rate Methane

(scfm) (%)

Rate

(scfm)

Methane

(%)

Rate

(scfm)

Methane

(%)

Mean

Standard deviation

Skewness

Minimum

Maximum

235.95 55.11

131.97 20.95

0.57 0.07

12.33 10.33

638.92 99.69

239.00

129.85

0.27

10.28

632.01

59.24

21.93

-0.18

10.68

100.00

231.08

123.96

0.59

10.34

644.47

53.24

20.52

0.09

11.51

99.51

C-I case: (completed panel; idle face)

HV coals MV–LV coals Sub-bit. coals

Parameter Rate Methane

(scfm) (%)

Rate

(scfm)

Methane

(%)

Rate

(scfm)

Methane

(%)

Mean

Standard deviation

Skewness

Minimum

Maximum

327.00

182.42

0.25

10.79

752.27

56.79

12.39

-0.11

15.84

85.59

337.76

189.05

0.20

10.29

799.73

58.40

11.54

0.01

21.63

90.27

339.11

187.62

0.27

11.26

753.13

54.15

12.23

-0.09

18.88

90.18

Table 7

Tabulated percentiles of the data presented in Fig. 9.

A-A case: (active panel; advancing face)

Coal rank HV bit. coals MV–LV coals Sub-bit. coals

Percentile Prod. rate Methane

(scfm) (%)

Prod. rate Methane

(scfm) (%)

Prod. rate Methane

(scfm) (%)

100% 638.93 99.7 632.02 100.0

90% 412.83 83.0 423.64 89.2

80% 349.24 72.8 356.71 80.7

70% 304.92 67.4 306.35 73.6

60% 258.08 61.5 264.17 66.5

50% 222.12 54.5 230.37 59.5

40% 187.89 48.9 200.85 53.3

30% 146.89 43.8 163.69 46.2

20% 112.62 36.2 118.27 39.3

10% 84.61 27.9 63.77 27.7

0% 12.33 10.3 10.29 10.7

C-I case: (completed panel; idle Face)

644.47 99.5

396.45 80.6

337.02 70.5

287.97 64.4

247.27 60.2

220.25 53.7

194.37 47.0

159.21 40.8

115.13 35.5

76.10 24.9

10.35 11.5

Coal rank HV bit. coals MV–LV coals Sub-bit. coals

Percentile Prod. rate Methane Prod. rate Methane

(scfm) (%) (scfm) (%)

Prod. rate Methane

(scfm) (%)

100%

90%

80%

70%

60%

50%

40%

30%

20%

10%

0%

752.27

588.18

504.49

419.28

368.88

313.42

261.46

219.31

149.38

94.22

10.79

85.6

72.7

67.9

62.9

60.2

57.4

54.3

50.0

46.3

41.3

15.8

799.73

603.49

526.59

456.71

385.81

328.16

264.30

215.96

155.91

85.18

10.30

90.3

72.7

67.7

64.1

61.5

58.5

55.7

52.1

48.7

44.2

21.6

753.13

617.09

521.88

437.03

387.35

327.82

259.03

214.21

165.39

97.34

11.26

90.2

68.7

65.7

60.5

57.2

54.1

51.5

48.1

43.9

38.0

18.9

The numbers in bold correspond to 20%, 50% and 80% percentiles, which were selected

as the low, mean and high values within the range.

Table 8

Parameters and their value ranges that are changed to evaluate responses of production

rate and methane concentration from a GGV in active and completed panels.

A-A case: (active panel; advancing face)

Parameter

Distance to SC top (ft)

Min; mean; max

(878.5; 1155.4; 1419.2)

Surface elevation (ft)

Average OB depth (ft)

(2095.0; 2502.7; 2855.4)

(1296.0; 1499.7; 1743.9)

C-I case: (completed panel; idle face)

Parameter Min; mean; max

Distance to SC top (ft) (966.3; 1190.5; 1423.5)

Surface elevation (ft)

Average OB depth (ft)

(2197.3; 2499.9; 2803.2)

(1338.6; 1503.7; 1700.7)

A 

0 

100 

200 

300 

400 

500 

600 

700 

800 

900 

100%

20% 70%80%90%

Production Rate (scfm) 

A-A:HV A-A:MV-LV 

A-A:Sub-Bit. C-I:HV 

C-I:MV-LV C-I:Sub-Bit. 0% 10% 30% 40%

50%60% 

Percentile 

B 

110 

A-A:HV A-A:MV-LV 

A-A:Sub-Bit. C-I:HV 90

%) 

C-I:MV-LV C-I:Sub-Bit. 

n (

70 

tratio

50 

Concen

30 

ane

10 

Meth

-10

20%70%80%90%

0%

10%30%40%50%60%

% 100

Percentile 

Fig. 9. Percentiles corresponding to different production rate (A) and methane

percentage outputs (B) from conditional MC simulations.

The ranges of GGV and mine operation parameters given in Table 5

are suitable to cover most longwall operations. The parameters most

likely to fall outside these ranges may be highly location-dependent

parameters such as the depths and surface elevations, and thus the

length of the slotted casing. In order to evaluate the effects on predicted

GGV production rates and methane concentrations in situations where

distance to slotted casing top, surface elevation, and overburden depths

are larger than given in Table 5, the allowable ranges of these

parameters were raised, following normal distributions. This gave the

minimum, mean and maximum values tabulated in Table 8 for A-A and

C-I conditions. Changing only these parameters, while keeping the

others constant, increases the length of the slotted casing.

MC simulations were performed with the data given in Table 8 for

A-A and C-I cases, keeping the values of the other parameters as in

Table 5. Simulation results are shown as percentiles in Fig. 10-A and B
for production rates and methane concentrations, respectively. The

simulation percentile, production rate, and methane concentration

data are given in Table 9.

Table 9

Tabulated percentiles of the data presented in Fig. 10.

A-A case: (active panel; advancing face)

Coal rank	 HV bit. coals MV–LV coals Sub-bit. coals

Percentile	 Prod. rate Methane

(scfm) (%)

Prod. rate Methane

(scfm) (%)

Prod. rate Methane

(scfm) (%)

100% 1137.38 99.7 986.19 100.0

90% 728.18 87.7 703.87 89.2

80% 629.79 81.0 626.74 83.5

70% 572.05 73.7 576.87 77.1

60% 525.25 68.1 533.82 72.6

50% 481.46 61.7 483.81 66.7

40% 445.42 56.3 446.70 60.9

30% 396.43 48.5 410.91 55.0

20% 343.01 41.0 346.38 48.0

10% 260.08 30.1 265.40 37.7

0% 13.90 10.7 46.45 10.7

C-I case: (completed panel; idle face)

1051.50 99.2

744.68 84.9

654.30 75.6

589.47 68.8

542.33 63.8

491.89 58.0

450.15 51.4

402.14 45.0

346.72 37.7

250.50 29.0

21.66 10.9

Coal rank HV bit. coals MV–LV coals Sub-bit. coals

Percentile Prod. rate Methane Prod. rate Methane

(scfm) (%) (scfm) (%)

Prod. rate Methane

(scfm) (%)

100%

90%

80%

70%

60%

50%

40%

30%

20%

10%

0%

817.73

451.11

399.32

342.22

269.67

213.56

173.12

139.68

96.79

44.46

0.11

93.9

73.6

69.6

65.8

62.3

58.3

55.2

52.3

49.0

42.2

28.0

672.94

436.42

371.87

309.27

254.34

213.33

167.12

121.84

83.09

41.00

1.46

97.4

79.4

73.1

69.4

65.4

62.2

59.4

56.4

52.7

48.0

30.7

718.53

488.99

403.91

326.96

279.85

230.81

177.39

140.95

89.37

48.97

0.21

85.7

71.7

67.3

63.5

60.1

57.8

55.0

52.3

47.9

43.1

27.9

The numbers in bold correspond to 20%, 50% and 80% percentiles, which were selected

as the low, mean and high values within the range.

Production rates (Fig. 10-A) showed a drastic

increase, especially for the A-A case, with the increased values of

distance to slotted casing, surface elevation, and average overburden

compared to Fig. 9-A.

A	 

0 

200 

400 

600 

800 

1000 

1200 

Production Rate (scfm) 

A-A:HV A-A:MV-LV 

A-A:Sub-Bit. C-I:HV 

C-I:MV-LV C-I:Sub-Bit. 0% 10% 20% 30% 40%

50%60% 

70% 70% 80% 80% 90% 90% 

100% 100% 

Percentile 

B	 

110 

A-A:HV A-A:MV-LV 

Methane Concentration (%)

A-A:Sub-Bit. C-I:HV 90 

C-I:MV-LV C-I:Sub-Bit. 

70 

50 

30 

10 

0%

10%20%30%40%50%60%

-10 

Percentile 

Fig. 10. Percentiles corresponding to different production rates (A) and methane

percentage outputs (B) from conditional MC simulations using increased depths

(Table 8) and the other parameters in Table 5.

Comparing the percentile values in Table 9 to

those in Table 7 shows that the rates were almost doubled. By

contrast, Fig. 10-B shows only slight increases in methane concentrations for both A-A and C-I cases compared to the results with the

lower range (Fig. 9-B).

5. Summary and conclusions

GGVs are integral to the safety of many underground longwall coal

mines, but their performances are difficult to forecast and, as a result,

it is a challenge to optimize their designs. The complex interactions

between GGV operating parameters, mining operating parameters

and a host of dynamic geologic conditions result in GGV performance

data that are difficult to understand and model, making stochastic

models a good choice for representing the data.

In this work, a set of deterministic models was derived using the

previously developed Methane Control and Prediction (MCP) software suite. The input parameters were assessed and adjusted so that

the models could be generalized to sites other than the northern

Appalachian basin, where the field data were collected. Methane

concentrations and production rates were generated using MCP with

random input values from a realistic range. Multilinear regressions

were performed on the data. Reasonable coefficients of determination

(R2

), from 0.69 to 0.88, were obtained. Next, it was determined that

slotted casing length, strata displacement, and gas content of the

overlying strata were important parameters for generalizing the

models away from site-specific data. Subsidence is difficult to

generalize across many sites, but it is a parameter that can

significantly influence gas production by increasing conductivity to

other gas bearing seams. Under this more generalized model, strata

displacement was considered instead, allowing for these mechanisms

to be taken into account. Strata displacement was only considered

over the length of the slotted casing, which is crucial because the

slotted casing promotes collection of gas from these fractured strata.

Their inclusion in the deterministic models resulted in coefficients of

determination ranging from 0.69 to 0.89. Comparison of statistical

data between the two deterministic models indicated good agreement, except for skewness in production rates. In both A-A and C-I

cases, absolute value of the skewness decreased, representing a more

normal distribution with the modified models.

Probability distributions for the input data were created for use in

Monte Carlo simulations of the modified deterministic models

presented above. The Kolmogorov–Smirnov (K–S) goodness-of-fit

test determined that the distributions could be represented with

normal and log-normal distributions. 5000 simulations for each of the

coal ranks (HV, MV–LV and sub-bituminous) in each of the longwall

cases (A-A and C-I) were run for a total of 30,000 simulations.

Conditional arguments were then applied to the resulting data in order

to check that input parameters and outputs were valid. Simulations

with invalid parameters or outputs were removed from the dataset.

Analysis of the output from the Monte Carlo simulations indicated

that overall production rates were generally higher for GGVs over

completed, inactive (C-I) longwall panels, but that methane concentrations were likely to be slightly lower in the produced gas stream for C-I

panels than for advancing, active (A-A) panels. Further, the results

indicated that in the 20th to 70th percentile range for the C-I case,
methane concentration ranged from 50 to 60%, a relatively small

variation. This indicates that changes in well operating parameters are

not as important during the completed and idle phase ofGGVproduction.

This research quantifies the likelihood of operating a GGV within a

given range of production rates and methane concentrations, using a

stochastic modeling method. This allows operators to quickly

determine if GGVs are operating within normal ranges and to explore

why they may operate outside of these ranges. These data can also be

used for preliminary design of GGVs. The advantage of this method is

that input parameters that cannot be accounted for in deterministic

modeling are incorporated into the stochastic model, so that a more

comprehensive model is generated. The detailed discussion of

development of the stochastic model allows for the procedure to be

transferred to specific sites or basins.

This work indicates that the stochastic approach is a promising

method to represent the randomness in GGV performances and to

improve the limited and site-specific character of deterministic

models. Future work may include analyzing data from other basins

to prove that this stochastic method can be generalized to other sites.

The initial deterministic model may be improved to generate a more

comprehensive data distribution for Monte Carlo simulation. The

stochastic model presented here for forecasting GGV performance can

be used both in design of longwall GGVs and in optimizing operating

parameters of existing GGVs to more efficiently extract methane from

working areas of underground mines.

List of symbols

Pi Production parameters (listed in Table 1)

xi Production parameter coefficients (listed in Table 1)

Ri Coefficients of determination (listed in Table 1)

yi Methane concentration parameters (listed in Table 1)

FRF Flow rate change factors

FP Flow percent

OB Overburden depth

D Distance to the top of the slotted casing

SL Slotted casing length

SD Strata displacement

μ Mean for normal and log-normal distributions

σ Standard deviation for normal and log-normal distributions

Conversion Table (English to SI units)

1 ft= 0.3048 m

1 MMscf=28316 m3

1 scfm= 0.0004719 m3/s

References

Ades, A.E., Lu, G., 2003. Correlations between parameters in risk models: estimation

and propagation of uncertainty by Markov Chain Monte Carlo. Risk Analysis 23,

1165–1172.

Alper, J.S., Gelb, R.I., 1990. Standard errors and confidence intervals in nonlinear

regression: comparison of Monte Carlo and parametric statistics. Journal of Physical

Chemistry 94, 4747–4751.

Chakravarti, I.M., Laha, R.G., Roy, J., 1967. Handbook of Methods of Applied Statistics,

vol. I. John Wiley and Sons, pp. 392–394.

Devore, J., Peck, R., 2001. Statistics: the Exploration and Analysis of Data. Duxbury

Thomson Learning Inc., Pacific Grove, California.

Diamond, W.P., Schatzel, S.J., 1998. Measuring the gas content of coal: a review.

International Journal of Coal Geology 35, 311–331.

Diamond, W.P., La Scola, J.C., Hyman, D.M. 1986. Results of direct-method determination of the gas content of the US coalbeds. Information Circular No. 9067, US Dept.

of Interior, US Bureau of Mines, Pittsburgh, PA.

Diamond, W.P., Ulery, J.P., Kravitz, S.J. 1992. Determining the source of longwall gob

gas: Lower Kittanning Coalbed, Cambria County, PA. Bureau of Mines, Information

Circular No: 9430.

Dougherty, H., Karacan, C.Ö., in press. A new methane control and prediction (MCP)

software suite for longwall mines. Computers and Geosciences.

Esterhuizen, G., Karacan, C., 2007. A methodology for determining gob permeability

distributions and its application to reservoir modeling of coal mine longwalls.

2007 SME Annual Meeting and Exhibit, Februrary 25–28, Denver, Colorado,

preprint 07-78.

Haghighat, A., Wagner, J.C., 2003. Monte Carlo variance reduction with deterministic

importance functions. Progress in Nuclear Energy 42, 25–53.

Huang, H., Hassan, A.E., Hu, B.X., 2003. Monte Carlo study of conservative transport in

heterogeneous dual-porosity media. Journal of Hydrology 275, 229–241.

Karacan, C.Ö., 2008. Modeling and prediction of ventilation methane emissions of U.S.

longwall mines using supervised artificial neural networks. International Journal of

Coal Geology 73, 371–387.

Karacan, C.Ö., 2009a. Forecasting gob gas venthole production performances using

intelligent computing methods for optimum methane control in longwall coal

mines. International Journal of Coal Geology 79, 131–144.

Karacan, C.Ö., 2009b. Elastic and shear moduli of coal-measure rocks derived from basic

well logs using fractal statistics and radial basis function networks. International

Journal of Rock Mechanics and Mining Sciences 46, 1281–1295.

Karacan, C.Ö., 2009c. Degasification system selection for U.S. longwall mines using an

expert classification system. Computers and Geosciences 35, 515–526.

Karacan, C.Ö., 2010. Methane control and prediction (MCP) software. http://www.cdc.

gov/niosh/mining/products/product180.htm.

Karacan, C.Ö., Goodman, G.V.R., 2009. Hydraulic conductivity changes and influencing

factors in longwall overburden determined by slug tests in gob gas ventholes.

International Journal of Rock Mechanics and Mining Sciences 46, 1162–1174.

Karacan, C.Ö., Esterhuizen, G.S., Schatzel, S.J., Diamond, W.P., 2007. Reservoir

simulation-based modeling for characterizing longwall methane emissions and

gob gas venthole production. International Journal of Coal Geology. 71 (2–3),

225–245.

Kelsey, A., Lea, C.J., Lowndes, I.S., Whittles, D., Ren, T.X., 2003. CFD modelling of methane

movement in mines. The 20th International Conference of Safety in Mines Research

Institues, South African Institute of Mining and Metallurgy.

Lu, Z., Zhang, D., 2003. On the importance sampling Monte Carlo approach to

uncertainty analysis for flow and transport. Advances in Water Resources 26,

1177–1188.

Lunarzewski, L., 1998. Gas emission prediction and recovery in underground coal

mines. International Journal of Coal Geology 35, 117–145.

Maier, H.R., Dandy, G.C., 2000. Neural networks for the prediction and forecasting of

water resources variables: a review of modeling issues and applications.

Environmental Modeling and Software 15, 101–124.

Mazza, R.L., Mlinar, M.P., 1977. Reducing methane in coal mine gob areas with vertical

boreholes. U.S. Department of Interior, Bureau of Mines Research Report No. 1607­

1-77.

McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for

selecting values of input variables in the analyses of output from a computer code.

Technometrics 21, 239–245.

Morin, A.M., Ficarazzo, F., 2006. Monte Carlo simulation as a tool to predict blasting

fragmentation based on the Kuz–Ram model. Computers and Geosciences 32,

352–359.

Palchik, V., 2003. Formation of fractured zones in overburden due to longwall mining.

Environmental Geology 44, 28–38.

Palchik, V., 2005. Localization of mining-induced horizontal fractures along rock layer

interfaces in overburden: field measurements and prediction. Environmental

Geology 48, 68–80.

Palchik, V., 2010. Experimental investigation of apertures of mining induced horizontal

fractures. International Journal of Rock Mechanics and Mining Sciences 47,

502–508.

Ren, T.X., Edwards, J.S., 2002. Goaf gas modeling techniques to maximize methane capture

from surface gob wells. In: De Souza, Euler (Ed.), Mine Ventilation, pp. 279–286.

Ripley, B.D., 1987. Stochastic Simulation. Wiley, New York.

Sari, M., 2009. The stochastic assessment of strength and deformability characteristics

for a pyroclastic rock mass. International Journal of Rock Mechanics and Mining

Sciences 46, 613–626.

Structured data, 2008. RiskAMP: risk analysis and modeling program. http://www.

riskamp.com.

Tomita, S., Deguchi, G., Matsuyama, S., Li, H., Kawahara, H., 2003. Development of a

simulation program to predict gas emission based on 3D stress analysis. 30th

International Conference of Safety in Mines Research Institutes, South African

Institute of Mining and Metallurgy, pp. 69–76.

Vose, D., 2000. Risk Analysis: a Quantitative Guide, 2nd Ed. Wiley, Chichester.