This paper quantifies gas-in-place, its uncertainty, and the displacements of major non-coal layers in the gas emission zones of longwall mines in a mining district in the Northern Appalachian Basin. These analyses will lead to improved design and application of methane capture and control systems.

In this study, the probabilistic gas emission zone prediction method presented in Karacan and Goodman (2011a) was combined with core data, including thickness and depth of major coal and non-coal formations, obtained from 276 exploration boreholes. Semivariogram analyses followed by sequential simulations were conducted to predict displacements in major non-coal layers and to predict gas-in-place in coals on a district-wide scale (in MMscf). Gas amounts were later specifically computed for a mined coal seam (Pittsburgh coal), in the caved zone, and in the fractured zone of the longwalls.

The results of simulations will be emphasized and analyzed in the Results and discussion section, although both kriging and simulations were performed and their results compared for variances and distributions. This preference is due to the inherent differences between these two techniques. Simulation differs from kriging in its superior ability to reproduce patterns of spatial continuity and produce realistic uncertainty modeling. Although there are no clear rules to exclude either kriging or simulation from a geostatistical analysis (Olea, 2009), kriging is recommended if the ultimate criterion is individual minimization of prediction error or the development of a smooth exploratory mapping of an unknown attribute. However, if the main objective is the correct assessment of confidence intervals or the correct modeling of spatial continuity, then simulation is the appropriate tool (Olea, 2009). Therefore, sequential Gaussian simulation was used as the preferred analysis technique in this study.

The geostatistical analyses used in this work gave not only the spatial correlations of thickness and gas content of coals (in scf/ton) that were used to calculate gas-in-place, but also enabled uncertainty analyses using the results of realizations at 5%, 50%, and 95% quantiles. Furthermore, grid blanking was used to isolate the volumes over the actual panels from the entire district to calculate gas-in-place and regions of higher rock displacement.

The studied longwall mining district is located in the southwestern Pennsylvania (in Greene County) section of the Northern Appalachian Basin (Fig. 1). The southwestern Pennsylvania section of the basin is a very important area for coal mining and for coalbed and coal mine methane capture and production. Markowski (1998) reported that 11 of the 24 intervals with coalbed methane in Pennsylvania are located in Greene County, Pennsylvania. Due to extensive coal mining in this area, these 11 gassy intervals (mostly coals) are closely associated with the mines operating in different coal seams and with their gas emission zones. Thus, these 11 coal intervals are potentially the major contributors to emissions in these mines.

The Monongahela Group of formations is within the gas emission zone when mining the Pittsburgh coal seam (Karacan, 2009b; Karacan and Goodman, 2011a). There are five main coals in the Monongahela Group of coal measures: the Pittsburgh coals, Redstone coal, Sewickley coal, Uniontown coal, and Waynesburg coals. The Pittsburgh coals consist of the main Pittsburgh seam and the Pittsburgh rider, which is usually located 1–3 ft above the Pittsburgh main coal. Overlying the Pittsburgh coal by 18–45 ft is the Redstone coal, which is not as extensive as the Pittsburgh. In general, the Redstone coal is 1–3 ft thick where present. The Sewickley coals include the Sewickley coal itself and any riders and splits of the main bench when it occurs. The deepest Sewickley coals are found in SW Pennsylvania at a depth of nearly 1170 ft. The Sewickley coal is laterally persistent and generally 1–6 ft thick. The Uniontown coal is not usually continuous, and when present, it is approximately 270–300 ft above the Pittsburgh coal with a thickness varying between 0.1 and 0.5 ft. The Waynesburg coals are made up of the Waynesburg coal, the Waynesburg A coal, and the Waynesburg B coal. The Waynesburg is laterally persistent and is usually multiple-bedded. These coal seams are shown in Fig. 1.

Pittsburgh coal is the main coal bed that is mined in the southwestern part of Pennsylvania. Due to their proximity to Pittsburgh coal, Pittsburgh riders and Redstone coals are usually within the caved volume and their gas is handled either by the ventilation system or by gob gas ventholes. Sewickley, Uniontown, and Waynesburg coals, on the other hand, are far enough from the Pittsburgh seam that they are part of the fractured section of the gas emission zone, as depicted in Fig. 1. These seams are believed to contribute to the gob gas in the fractured system (Karacan et al., 2007) and their gas emissions should be captured using gob gas ventholes.

Other non-coal formations of the Monongahela Group that do not contain significant amounts of gas but may affect caving and fracturing of the gas emission zone are shales of various kinds, limestones, and sandstones. These formations are inter-bedded with the coal seams and exist in varying thicknesses depending on the location. Based on the identifications given in the exploration boreholes, the thicknesses and their relatively high rock strengths compared to other strata, the Fishpot limestone, Sewickley sandstone, and Uniontown sandstone were selected for further geostatistical analyses (Fig. 1).

The mining district that was modeled in this work hosted panels N-1 through K (Fig. 2). This mining district was selected for study because mining was already completed in this area, because of data availability from the exploration boreholes and due to the previous reservoir modeling study detailed in Karacan et al. (2007), which evaluated the productions of gob gas ventholes in panels G to K and determined those factors affecting their performances. The dimensions of the area shown in Fig. 2 are 20,500 ft in the y-direction and 18,900 ft in the x-direction.

In this district, overburden depths ranged between 500 and 900 ft. Longwall panels in the primary area were initially 830 ft wide and were increased to 1000 ft starting with F Panel. Thus, the panels are super-critical—i.e., the panel width is greater than the depth of the overburden, which results in a more complete caving of the overburden strata into the mined void. Several coal seams presented in the previous section are within the gas emission zone of the longwall panels shown in Fig. 2, and are believed to be the primary source of caved and fractured strata gases in the area.

Methane control in the study longwall district included the bleeder ventilation system, gob gas ventholes, and underground horizontal methane drainage boreholes drilled prior to mining. The bleeder ventilation system included the peripheral bleeder entries surrounding the panels, the former gateroads between the mined-out panels, and the associated bleeder fan shaft(s). The bleeder fans operated in this district were located at the top of 6-ft-diameter air shafts. Gob gas ventholes of the panels in this mining area were generally drilled to within 40–45 ft of the top of the Pittsburgh coal seam, cased with 7–8 in. of diameter pipe, and finished with 200 ft of slotted pipe on the bottom (Karacan et al., 2007).

Fig. 3 gives methane production rates from longwall bleeders of N–F panels (BF-2) and G–K panels (BF-3) during mining operation, with a missing period between the 300th and 700th days of mining. This figure also shows the methane production rates from gob gas ventholes drilled over these panels to remove gas from fractured strata. The data show that both bleeder rates and gob gas venthole rates differ between mining of these two panel districts. During mining of N–F panels, average methane production from GGVs is about 200 Mscfd, as opposed to 500–600 Mscfd from GGVs of G–K panels. This difference may be related to various factors discussed in Karacan (2009a) and differences in properties of the gob as a reservoir (Karacan, 2009c; Karacan and Goodman, 2011b). Similarly, methane rate from BF-2 averages around 2500–3000 Mscfd, where it is 1500–2000 Mscfd from BF-3. This may be related to the capacity of the bleeder fans. However, from both GGVs and bleeder fans, methane rates may be related to the availability of methane from various horizons and its spatial distribution as well, as explored in following sections.

The gas content of coal can be measured or estimated using various techniques. These techniques usually fall into two categories: (1) direct methods which actually measure the volume of gas released from a coal sample (preferably wire line core samples) sealed in desorption canisters, and (2) indirect methods based on empirical correlations or laboratory-derived gas storage capacity data from sorption isotherms. An extensive review of direct techniques for gas content measurement for coal has been published by Diamond and Schatzel (1998). Properly conducted direct-method testing of coal cores provides relatively accurate estimates of in-place gas contents of coals for most mine planning purposes and allows for resource evaluation at a reasonably low cost. A modified direct method (MDM) procedure (Ulery and Hyman, 1991) provides an increased level of accuracy, but requires a higher level of instrumentation sophistication, procedural complexity, and cost.

Alternatively, regional gas content data on individual coal samples can be used along with auxiliary data on coal rank and/or depth to construct curves, as shown in Fig. 4, for estimating in-place gas contents for a specific area. Such curves are generated for a particular coal seam or closely associated group of coal beds and can be used to estimate gas content values only if the rank or depth is known for the coal bed of interest. The graph given in Fig. 4 presents such curves for bituminous coals of the Black Warrior Basin, Alabama (McFall et al., 1986). However, in considering Fig. 4, it should be kept in mind that lithotype characteristics of coal also affect the methane content. Gas content generally varies positively with the amounts of vitrinite and liptinite, which usually offer high methane storage capacity. No obvious relationship is observed with inertinite content. Therefore, the smooth graphs in Fig. 4 should be interpreted and used with caution.

The geostatistical estimates of gas-in-place presented in this paper relate to coals of the Lower Monongahela Group. Site-specific information used to determine gas content of the coals and the gas content estimation technique were based on bivariate normal distributions presented in Karacan and Goodman (2011a).

In Karacan and Goodman (2011a), the gas contents of all seams given in Diamond et al. (1986) in Pennsylvania were compiled in a database, classifying all the coal beds above and including the Pittsburgh seam. Means (μ) and standard deviations (σ) of gas content and overburden depth for each coal studied in this work and the correlation coefficients (ρ) between variable pairs that were used in calculating joint probability distributions for gas content and depth are given in Table 1. All coals within the interval of interest in this study area are high volatile bituminous-A type coals with the gas content-depth relationship generally showing an increasing trend with depth. However, it should also be noted that gas contents show significant variations even at the same depth (Fig. 5) to establish a clear regression.

The bivariate normal probability calculations and representations of probabilities for obtaining gas contents less than or equal to 350 scf/ton with four-parameter logistic functions were used to estimate gas contents. These analyses showed that gas content of coals, most of which were from southwestern Pennsylvania, could be represented by the values of a dependent variable at 50% value of the independent variables (50% quantile—Q50) of these sigmoid functions, up to an overburden of 1400 ft. The polynomial representation of these coefficients gives a Pearson’s correlation coefficient of 0.9969 and is also given in Table 1. This relation was later used in this study to calculate the gas content of coals at Q50 as a function of depth.

Monte Carlo (MC) simulations were performed in order to determine similar relationships for 5% and 95% quantiles (Q5 and Q95) for each coal bed. These simulations were conducted on the models that aimed systematically to generate the distribution of original gas contents shown in Fig. 5, as well as the distribution of overburden values in a bivariate normal distribution. In this process, the correlation coefficient between the data pairs was matched, and the values of Q5 and Q95 were determined from the resultant distribution of population. After this step, the coefficients of the quadratic equations, which were initially arbitrary for Q5 and Q95, were iteratively determined by changing overburden depths and by targeting the Q5 and Q95 values in a goal-seek manner to minimize the standard errors of predictions. For these runs, 20,000 simulations were performed. The results of these analyses and the equations for Q5 and Q95 are also given in Table 2.

One of the key steps in forecasting gas emissions during and after mining is to calculate the volume of GIP that will potentially migrate to the underground mining environment. The simplest method for calculating the GIP for coal seams is based on commonly available geologic mapping data for the mine site and using site-specific gas content data (Diamond, 1982), as follows: (1)GIP=(ρ×h×A)GC where GIP is gas-in-place in coal; ρ is coal density; h is coal thickness; A is the area that coal occupies; GC is the gas content of coal (volume-to-mass ratio), in ft³ gas/ton of coal.

A more elaborate approach is a reservoir analysis-based volumetric approach (Karacan and Diamond, 2006) given in Eq. (2). Volumetric-based calculations may also be used to determine the gas-in-place for coal beds in a specific geographic area. This approach both relates the volume of gas in the reservoir at reservoir conditions to the volume at surface conditions, treating separately the volume of free gas in fracture porosity and adsorbed gas in bulk volume.

In this work, a pseudo-volumetric method was used to calculate gas-in-place for coal seams. The final form of the equation, without including moisture content (f_m), which is used after geostatistical modeling is: (3)GIP=AhGCρ(1-fa).

In this equation, GC is the gas content of coal (determined as explained in Section 3.1 for each coal seam). “A” is the area, “h” is thickness of a particular coal bed of interest at a specific point, f_a is the ash content of coal and ρ is its bulk density.

Ash in coal is a combination of all non-organic matter including clays and minerals. Ash content is strongly dependent on the primary paleo-depositional environment, the changes in the climate and various precipitations during the peat-forming process, and the present hydrogeology of the area at a particular depth. Although there is no a priori reason that ash content of coals should correlate with depth, it has been observed in West Virginia that ash content does indeed decrease on average for coals in the Monongahela Group from the Waynesburg coal to the Pittsburgh coal (West Virginia Geological and Economic Survey, 2005), similar to Fig. 6.

Ash is considered as a constituent of coal that does not have a significant methane sorption capacity. Therefore, its amount should be determined for GIP calculations (Eqs. (2) and (3)). In addition, ash has a higher density than organic matter (112.4–137.4 lb/ft³; 1.8–2.2 g/cm³) and thus ash content is one of the factors that most influences coal density.

Ash contents of coal seams within a thick stratigraphic column usually follow a scattered trend, making it difficult to establish a simple correlation (Fig. 6). In addition, there are variations within each seam as a function of overburden depth, as observed from scatter of the data given in Fig. 6 and from the correlation coefficients presented in Table 2. Therefore, a probabilistic approach of uncertainty based on ash-content distributions of the individual seams over the depth they exist seems to be the best method to determine ash contents of individual coals. This method requires bivariate distribution analyses of two random functions for determining joint probabilities of occurrences of different values. In this study, these functions were assumed to be normal for distributions of depth and ash contents of individual seams.

Data for ash contents for the coals of this study were taken from Diamond et al. (1986). The data were partitioned for each seam and was analyzed using bivariate normal distributions to determine ash content at Q5, Q50, and Q95. For these analyses, a similar technique described in Section 3.1 for gas contents was performed using a combination of bivariate distribution analyses and Monte Carlo simulations. These analyses effectively treated for each seam such that at each depth ash content had a distribution function, from which probable values could be determined at Q5, Q50, and Q95 to generate functional relationships. The values of ash at Q5, Q50, and Q95, as well as their standard errors from MC simulations and the equations that represent ash values as a function of depth, are given in Tables 3 and 4.

Ash content values calculated for individual seams at different quantiles were later converted to their bulk densities using the correlation that was established for these same coals from the data given in Fig. 7. The correlation between ash content and bulk density is given in Eq. (4). Ash and density correlations were built and used in the absence of spatial data of these attributes from analyses of cores of the exploration boreholes.

As discussed in Section 2, major non-coal strata within the possible gas emission zone in the area of study are the Fishpot limestone (FPLS), Sewickley sandstone (SWSS) and Uniontown sandstone (UNSS). Properties of these rock units are important for movement of gas within the gas emission zone. If these strata subside and fracture vertically and horizontally (termed “displacement” collectively in this paper), then all sources of gas-in-place will be in communication with each other, and potentially will be in communication with the mine as well. If these strata do not displace significantly, then gas sources may be isolated from each other. Correlating maximum displacement locations of major non-coal strata with high gas-in-place locations in the gas emission zone helps in predicting those locations of high gas inflows into the mines and can best locate gas capture boreholes.

Displacements were globally described in accordance with Karacan and Goodman (2011a) using the field data presented in Mazza and Mlinar (1977) for southwestern Pennsylvania. These data, which show that strata displacements are a function of distance from the top of the mined coal seam (the Pittsburgh seam in this case) are given in Fig. 8. This figure shows that displacements are at a maximum of ~2–3 ft close to the mining zone and decrease with increasing distance from the mined seam. Eventually, fracturing and separations cease around 350 ft above the mining horizon.

Strata displacements and their probabilities were described in Karacan and Goodman (2011a). Means (μ) and standard deviations (σ) of displacements and distances from the Pittsburgh seam and the correlation coefficients (ρ) between these variable pairs that were used in calculating joint probability distributions for these variables are given in Table 5.

The bivariate normal probability calculations and representations of probabilities of observing displacements ≤4 ft with four-parameter logistic functions were used to quantify displacements. Results showed that displacements in southwestern Pennsylvania, as a result of mining of the Pittsburgh seam, could be represented by the values of the dependent variable at 50% value of the independent variable up to distances of 350 ft. The polynomial representation of these coefficients gives a Pearson’s correlation coefficient of 0.9952 (Table 5). This relation was used to calculate displacements of FPLS, SWSS, and UNSS as a function of their distances from the Pittsburgh seam at exploration borehole locations, which later were used for geostatistical modeling of strata displacements in these three major rock layers.

The data used in geostatistical modeling were obtained from 276 vertical exploration boreholes drilled over the mining area shown in Fig. 2. These boreholes and their spatial locations are shown in Fig. 9. For modeling that will be discussed in detail in the forthcoming sections, these data were assigned to a 100 × 100 Cartesian grid in which each grid was 189 ft in the x-direction and 205 ft in the y-direction. In preparing the data, the thicknesses of each identified layers of each of the coal and non-coal units were combined, so that each unit would be represented as a single layer. Accordingly, the average of the depths of these layers was taken as the representative depth of each unit, as were their depth-dependent attributes.

All empirical data from the boreholes had identified coal and non-coal formations drilled along the entire length of the boreholes. Since the top of the gas emission zone for mines operating in the Pittsburgh seam was around 350 ft above the top of the Pittsburgh coal (Karacan and Goodman, 2011a; Mazza and Mlinar, 1977), the data beyond this interval were excluded from further analyses. Thus, Fig. 9 shows only the coal and non-coal layers within the expected gas emission zone. In this figure, each data point represents an identified formation, with different colors representing their elevation using the top of the Pittsburgh seam as a reference datum.

To calculate gas-in-place in coals and in different sections of the gas emission zone, the coal beds within the 350-ft interval above the Pittsburgh seam, and including the Pittsburgh seam, were identified and isolated from the rest of the data as separate datasets. These coal beds were the Pittsburgh seam (main bench—PCMB), Pittsburgh rider seams (PCR), Redstone coal (RDC), Sewickley coal (SWC), Uniontown coal (UNC), and Waynesburg coals (both upper and lower—WBC). These data were further analyzed carefully to extract erroneous data points and to see whether they appeared as a single layer or as splits. In case of splits or multiple layers of the same seam, data were combined as noted before. Therefore, the maximum number of data points for each unit was equal to the number of boreholes when a particular seam was present in all boreholes, and the number of data points was less than the number of boreholes when the seam was missing at a particular location (Table 6).

For coal seams of interest, three attributes were determined at each spatial position for geostatistical modeling and for calculation of gas-in-place using Eq. (3). These attributes were overburden depth, thickness of coal, and gas content. Results of univariate statistical analyses of these data for each coal seam are given in Table 6.

A similar data extraction and preparation procedure was also performed for non-coal layers as well. For these layers, two attributes were determined at each spatial position for geostatistical modeling. These attributes were overburden depth and displacements, calculated using distance from the Pittsburgh seam at each borehole location. Results of univariate statistical analyses of these data for each of the non-coal layers are given in Table 7.

The frequency distributions of the data (histograms) of each of these attributes for both coal seams and non-coal layers were also checked for their Gaussian behavior. Olea (2009) and Remy et al. (2009) state that for some of the geostatistical modeling and simulation techniques to be applicable without the need of transformation of the data, the data should follow a Gaussian (normal) distribution. Plotting histograms of depth, thickness, and gas content data for coals and of depth and displacements for rock layers showed that the data were not exactly Gaussian in all cases. Therefore, automatically, all distributions in this work were transformed to normal scores by targeting a Gaussian distribution with mean 0 and variance 1. Semivariograms were then modeled on the normal-score data. However, the values of the attributes were transformed back to the original space during simulations by targeting the original distribution.

Examples of these histograms for thickness data of Waynesburg (WBC) and depth of Sewickley (SWC) coals, and normal-score transformation of these attributes are shown in Fig. 10.

The Stanford University Geostatistical Modeling Software (SGeMS) was used for spatial correlation analyses and for stochastic (sequential Gaussian) simulations. SGeMS implements several geostatistics algorithms for modeling of earth systems and phenomena that exhibit space-time distributions (Remy et al., 2009). The software also has a graphical user interface and a script compiler that gives the user ease and flexibility to implement complex problems. It includes several of the algorithms that are in the GSLIB—Geostatistical Software Library (Deutsch and Journel, 1998) in its platform.

All modeling studies and their results require some level of verification. In this study, the results of sequential Gaussian simulations of modeled attributes were compared with the original data before proceeding with calculations of GIP and the associated uncertainties. These comparisons required development of histograms and Q–Q plots of hard data, along with SGSIM realizations and calculation of elementary descriptive analyses. The analysis procedure is schematically shown in Fig. 13 for thickness of the Pittsburgh seam and was applied to each of the simulated parameters for evaluations of results.

Basic statistical analyses of the data obtained from E-maps of simulated attributes of coals and non-coal units, with the hard data given in Tables 6 and 7, showed that almost all of the basic statistical parameters are similar. This indicates that average results from simulations have similar properties in relation to the ranges and limits of the hard data. However, it should be recognized that the variances obtained from E-data are less because the point-wise averages of simulated data in E-maps have less dispersion compared to both hard data and to individual realizations. For instance, Fig. 14 shows example histograms of Waynesburg coal thickness and Sewickley coal depth plotted from the 50th realization and from the E-type maps (point-wise average of 100 realizations). Histograms of the hard data for the same properties and formations are given in Fig. 10. Comparison of distributions from the presented realizations and E-type maps with the hard data shows that the histograms of individual realizations are very similar to the hard data, since they are generated to closely match the cumulative probability distribution of the original data. On the other hand, E-type data have less dispersion due to averaging of results from a large number of realizations compared to the original data and to individual realizations, as expected. This was the case for all modeled attributes. Therefore, individual realizations and properties derived from those attributes were used to calculate gas-in-place and uncertainties.

In addition to comparisons of basic statistical parameters outlined above, Q–Q plots were prepared for each modeled attribute for their hard data versus E-type data quantiles and for hard data versus individual realizations. Q–Q plots are used to compare probability distributions. A straight line is an indication of equality between the distributions being compared and that the data in both axes have similar quantile values (Krishnamoorthy, 2006).

The Q–Q plot analyses conducted in this work gave acceptable linear trends between hard data and realizations, as well as between E-type maps. Fig. 15 shows, as an example, this comparison for thickness of the Waynesburg coal. This figure shows that comparisons of both types of SGSIM data give nearly linear relations with the hard data from which they are generated, and thus the probability distributions of all data are almost the same. Note the smoothing in the E-type map, revealed by higher modeled values when the empirical values are low, and lower modeled values when the empirical values are high. An E-type map shares a similar smoothing effect as a kriging map.

Since spatial features and uncertainty play an important role in predictions of GIP and methane emissions that result from the gas emission zone, individual realizations, rather than E-type maps, of modeled attributes from SGSIM were used to create realizations of GIP for further analyses and computations.

GIP and displacement calculations for the gas emission zone and the uncertainty evaluation strategy were as described below.

In addition to assessing these three distinct zones of deformation for the entire model area, GIP in the following longwall panel series was also calculated for N–F panels, and G–K panels. The reason why these two panel series were evaluated separately is that the methane from the ventilation system was controlled with two separate bleeder fans (Figs. 2 and 3).

The values of GIP for individual coals and longwall gob zones, as well as displacements from realizations, were calculated at 5%, 50%, and 95% quantiles (Q5, Q50, Q95). These quantiles represent the values where each estimated value is ranked at 5%, 50%, and 95% in the distribution of results, and the probability of the estimated values is lower than the actual unknown values. Therefore, the data of any attribute at 5% quantile are the estimated values at 5% probability, and can be considered as the lower limit. Similarly, Q95 is the data with 95% chance of being higher than the actual value and is considered as the upper limit. Between these, Q50 has a special meaning, representing the median of the possible population distribution. Q50 is called M-type where each estimated value has a 50% chance to be higher than the actual value.