A new method to calculate permeability of gob for air leakage 

calculations and for improvements in methane control 

C.Ö. Karacan 

National Institute of Occupational Safety and Health (NIOSH), Office of Mine Safety and Health 

Research, Pittsburgh, Pennsylvania 

ABSTRACT: The porosity and permeability of the caved zone in a longwall operation impact many ventilation 

and methane control related issues, such as air leakage into the gob, the onset of spontaneous combustion, 

methane and air flow patterns in the gob, and the interaction of gob gas ventholes with the mining environment. 

Despite their importance, the gob is typically inaccessible for performing direct measurements of porosity and 

permeability. Thus, there has always been debate on the likely values of porosity and permeability of the caved 

zone and how these values can be predicted. This study demonstrates a fractal approach that allows the calculation 

of porosity and permeability from the size distribution of broken rock material in the gob, which can be 

determined from image analyses. The fractal approach presented here constructs flow equations and fractal 

crushing equations for granular materials to predict porosity for a completely fragmented porous medium. The 

virtual fragmented fractal porous medium so generated is exposed to various uniaxial stresses to simulate gob 

compaction and porosity and permeability changes during this process. 

1 Introduction 

Longwall mining is an underground mining method that 

can maximize coal production in coal beds that contain 

few major geological discontinuities, like faults. Longwall 

coal mining causes large scale disturbances of the 

surrounding rock mass due to fracturing and caving of the 

mine roof behind the shields as the mine face advances. 

The gob is highly fragmented since the overlying rock 

layers fall into the mine void and are broken into irregular 

shapes of various sizes. The gob can contain high void 

ratios due to fragmented rock pieces and may provide high 

permeability flow paths for any fluid flowing from 

surrounding formations into the mining environment. 

The porosity and permeability of the longwall gob can 

significantly affect methane and air flow patterns. Since 

methane is a major hazard in underground coal mining 

operations, extensive methane control techniques such as 

gob gas ventholes (GGV) are employed to supplement the 

existing mine ventilation system. Gob gas ventholes are 

drilled vertically into the overburden above longwall 

panels to capture the gas released by the fractured strata 

before it enters the mining environment. Most GGVs are 

drilled to within a short distance of the mined coalbed to 

stay out of the high-permeability gob so as not to draw 

excessive amounts of ventilation air along with the 

hazardous methane. In addition, the flow of mine air into 

the gob may increase the risk of spontaneous combustion. 

It has been shown that with an increase in gob 

permeability, airflow increases in the broken, coal-rich gob 

causing the onset temperature for spontaneous combustion 

to be reached more quickly (Yuan & Smith, 2008). 

As mining progresses, the gob gradually consolidates 

sufficiently to support large loads resulting from the 

overburden weight (Pappas & Mark, 1993). Consolidation 

results in a reduction in the void ratio (porosity) and the 

associated permeability. Although reduced to some degree 

due to compaction, prevailing high permeability pathways 

in the gob still affect the leakage of ventilation airflow 

from the face into the gob, the flow of methane from 

surrounding sources into the gob and into the mine, and the 

performance of methane extraction GGVs. Thus, an 

understanding of compaction phenomena and the resultant 

reservoir properties of gob material is very important for 

developing adequate methane control strategies. 

The properties of gobs and their influences on methane 

flow and control have been investigated through combined 

geomechanical and reservoir modeling (Esterhuizen & 

Karacan, 2005; Karacan et al., 2007) and through history 

matching of GGV production. In these studies, 

geomechanical models were used to characterize the stress 

changes due to mining, which were converted to porosity 

and permeability data. Esterhuizen & Karacan (2007) 

presented the development of a new methodology to 

determine directional variations in gob permeability based 

on a model of caving and block rotation which considered 

the effect of block dimensions and fall height on porosity. 

This work presents a fractal porous medium model to 

predict gob porosity and permeability for controlling 

methane flow and production from GGVs and to predict 

air leakage into the gob. For this approach, the particle size 

data for simulated gob material given by Pappas & Mark 

(1993) were used to calculate fragmentation fractal 

dimensions and particle size distributions before and after 

controlled loading tests. Fractal geometry concepts were 

used to construct porosity and flow equations for a 

completely fragmented porous medium through 

fragmentation, tortuosity and area fractal dimensions 

(Karacan & Halleck, 2003). The fragmented fractal porous 
 

medium model was exposed to various uniaxial 

compressive stresses to simulate gob compaction and to 

predict porosity-permeability changes in the gob. 

2 Laboratory Tests for Gob Material Simulation 

In this paper, rock fragment size distribution data was 

taken from Pappas & Mark (1993) for model validation 

purposes. In order to better understand the behaviour of the 

gob and to provide numerical modellers with more 

accurate data for simulating longwall mining conditions, 

they conducted a laboratory study to estimate the gradation 

of actual gob material. 

The first phase in the laboratory tests developed test 

materials with properties similar to those of actual gob 

material. The characteristics that were considered were the 

tensile and compressive rock strengths, rock density, 

surface roughness, rock shape, rock size, and size 

gradation. Most of these characteristics would be 

simulated by broken rock obtained from fresh roof falls. 

However, the rock size and gradation needed to be reduced 

to a laboratory scale by shifting the particle size 

distribution curve of the actual gob material determined 

from pictures taken from the headgate entries, where 

portions of the gob could be viewed. This shift was 

performed parallel to the horizontal axis of the particlesize curve until the desired maximum particle size of the 

laboratory sample was reached. The detailed procedure for 

this can be found in Pappas & Mark (1993). However, 

particle size affects the permeability of porous medium. 

Thus, the particle sizes used by Pappas & Mark (1993) 

were scaled up for use in the developed model to calculate 

permeability. 

The rock types that were obtained to simulate gob were 

shale, weak sandstone and strong sandstone. In this paper, 

the size distributions obtained with maximum particle size 

of 5.1 cm (2.0 inches) for shale were used for modeling. 

Table 1 gives the test summary of simulated shale gob 

material. 

Table 1 Summary of simulated materials (from Pappas 

and Mark, 1993) analyzed in this paper and their 

loading tests. 

Rock Max 

size 

cm 

Max 

stress 

MPa 

Initial 

por. 

Por. 

@ 5.4 

MPa 

Por. 

@ max. 

stress 

Shale 5.1 19.0 0.802 0.365 0.160 

Compressing a fragmented rock medium causes 

departures from the initial particle size distribution for 

both shale and sandstone material due to further crushing 

and fragmentation of the particles. This changes porosity 

and permeability. This will be discussed later in the rest of 

this paper. 

3	 Model Development 

3.1	 Fractal Fragmentation and Formulation of a 

Deforming Fractal Porous Medium 

Fragmentation may be considered as structural failure of a 

brittle material caused by multiple fractures of different 

lengths (Perfect, 1997). The disordered nature of pores and 

grains in fragmented porous medium suggests that the 

structure created by the fragmentation process shows 

scaling properties (Weiss, 2001). Fractal models not only 

can describe the scaling of mass and surface roughness of 

individual fragments, but also the fragment size 

distributions (Perfect & Kay, 1995). 

The power law scaling of the number-particle size 

relation given by Turcotte (1986): 

N(ω > Ω) = B Ω− DF 

F (1)

where: N is the cumulative number of particles of size 

greater than a characteristic size Ω, the exponent DF is 

referred to as the fragmentation fractal dimension and 

includes the information about the scale dependence of the 

number-size distribution of particles, and BF is a 

coefficient related to number of particles of unit diameter. 

In this study, the porous medium forming the gob was 

assumed to have pore sizes that can be considered as 

bundles of capillary tubes of different diameters as 

suggested by Yu & Cheng (2002). If the diameter of a 

capillary is assumed to be η and its tortuous length along 

the flow direction is Lt (η ), then the tortuous length will 

be longer than its representative length LO, which runs the 

shortest distance between two points. Following 

Wheatcraft & Tyler’s (1988) expression and scaling 

approach for flow in a heterogeneous medium, there is a 

fractal scaling between the diameter and length of the 

capillaries (Yu & Chang, 2002): 

L (η) η 1 DT

t = − L DTo (2)

In this equation, DT is the tortuosity fractal dimension 

which can take values between 1 and 2. DT represents the 

convoluteness of flow channels. When DT is 1, the flow 

channels are straight. The tortuosity increases with 

increasing values of DT. In the limiting case of DT =2, we 

have a highly tortuous line that fills a plane (Wheatcraft & 

Tyler, 1988). However, Wheatcraft & Tyler (1988) 

reported that DT values more than 1.5 do not have any 

physical significance. 

For the bundle of capillaries or flow channels, the 

number of channels with size η is also important. Since 

channels in a porous medium are analogous to the islands 

in a sea or spots on surfaces, the cumulative size 

distribution of pores can be written as: 

⎛η ⎞

DP 

N(L ≥η)
 =
 ⎜ max ⎜ ⎟

⎟ (3)

⎝ η ⎠
In this equation, N is the number of pores whose sizes 

(η ) are greater than a characteristic size, L and η max is the 

maximum pore size. Equation 3 can be differentiated to 

give: 

− dN = D ηDP η −( DP +1)

P max d (η) (4)

In Equation 4, DP is the pore area fractal dimension (1< 

DP <2). A DP value of 2 represents a regular pore area and 

the irregularity increases for DP <2. These equations 

describe the scaling relationship of cumulative pore 

populations. The total number of pores from the smallest 

diameter, η min, to the largest diameter, η max, thus can be 

obtained from the Equations 3 and 4 as: 

⎛η
 ⎞

DP

N t (L ≥η ) = max

min ⎜ ⎟ ⎜ ⎟ (5)

⎝ η
 min ⎠

Dividing Equation 4 by Equation 5 obtains an 

expression for the probability density function: 

dN − = D ηDP η −( DP +1)

P min d (η ) = f (η)d (η) (6) 

N t 

where 

f ( ) η = D D

Pη P 

min η
−( DP +1 ) (7)

The equations above can be used as a model of a 

capillary tube bundle of any shape. 

In order to characterize a pressure-driven, steady state 

flow of a fluid through long, straight, and rigid channels of 

any constant cross-sectional shape, with no-slip boundary 

conditions, the Hagen-Poiseuille (H-P) flow equation can 

be used. H-P equation for a circular area: 

πΔPη c 

Q i = (8) 16αμLt 

where η is the hydraulic diameter and Lt is the total 

length of the channel. In circular pipes, α=8 and the flowchannel diameter exponent “c” of the equivalent channel 

diameter is 4. For some known shapes such as rectangle, 

the shape correction factor (α) is given by Mortensen et al 

(2005): 

π 3γ 2 ⎛ − ∞ ⎛ n ⎞⎞ 1 γ
 2 πγ
 α
 =
 ⎜ 
 −
 ⎟ ⎜ ⎟ 8
 ∑ 4 5 5 tanh ⎜ ⎟ (9) 

⎝ n=1,3,5... πn π
 n ⎝ 2 ⎠⎠

where γ is the width-to-height (w/h) ratio of the 

rectangular channel. Tanh(x) can be approximated to 1 if x 

is very large which can be satisfied if γ is very large (slit 

shaped channels for example). In this case Equation 9 

becomes; 

12π 5γ 2 

α = 5 (10) 

π γ −186ζ (5) 

where ξ(5) results from truncating the summation after 5th 

iteration in Equation 9. 

In this paper, the average diameter exponent from Arya 

et al.’s work (1981) will be used in the H-P equation 

(Equation 8) to describe flow in single channels (3.531). 

The variable “α” (Equations 9 & 10) will be retained in the 

equation to run parametric analyses related to channel 

shape. Thus, H-P equation from a single channel becomes: 

Q i () π ΔP η3.531η = (11) 16 Lt (η)α μ 

If the total flow rate is calculated for each channel by 

integrating between the minimum and maximum sizes over 

the entire range of pore sizes, the total flow rate is obtained 

as;

η max π ΔP L − DT D

Q o P

T = − ∫ Q i (η) dN (η) = 

η min 

16 μ α 2.531 + D T − DP 

⎡ 2.531+ D T −DP ⎤

× η2.531+ DT ⎛η max ⎢1 − ⎜ min ⎞

⎟ ⎥ ⎜ ⎟ ⎢ ⎝η max ⎠ ⎥ ⎣ ⎦

 (12) 

For the special case of circular channels, this equation is: 

π ΔPL− DT D QT = o P ×

128 μ 3 + DT − DP

⎡ ⎛η ⎞

3+DT −DP

 (13) ⎤

η 3+DT

max ⎢1 − ⎜ min ⎟ ⎥ ⎜ ⎟ ⎢ ⎝η ⎣ max ⎠ ⎥

⎦

Equations 12-13 honour the properties of H-P equation 

where the streamlines are formed by individual particle 

sizes that create tortuous pores of a particular dimension 

and shape and those they do not really mix. Combining 

kA ΔP

Equations 12 and 13 with Darcy’s Law ( Q T 

T =
 ),

μLo 

where k is the permeability and AT is the total area to flow, 

yields the definition of permeability for a porous medium 

composed of tortuous channels of varying diameters and 

shapes: 

π L

1−DT

o D k = P η2.531+DT

16
 A α
 max

T 2.531
 +
 DT −
DP

2.531+D T −D (14) ⎡ ⎛ ⎤ η
 ⎞ p 

×
 ⎢1 −
⎜ min ⎟ ⎥ ⎜ ⎟ ⎢ ⎝η
 ⎠ ⎥ ⎣ max ⎦
For circular channels; 

+ D −Dp

π L 3 1−D T T

o D ⎡ ⎤ P 3 +DT ⎛η k = η ⎢1− ⎜ min ⎞ ⎥ (15) ⎜ 128 AT 3 + D max ⎟

⎟ T − DP ⎢ η max ⎠ ⎥ ⎣ ⎝ ⎦

In these equations, total area of flow (AT) is the 

combination of areas for both pores and the rock surface 

exposed to the flow. Since the gob material in the caved 

zone is highly fragmented and the total channel area is the 

same as the effective porosity where the flow occurs, total 

channel area (AP) is used for AT. In order to determine the 

value of AP, the areas of individual flow channels can be 

summed mathematically to give (Karacan, 2009): 

s πΩ2 ⎞ A =
 A =
 ∑⎛

⎜ i φ N1−D

⎜ T

T p i ⎟ (1

i=1 6 i ⎟ 6)

⎝ ⎠

as a combination of virtual pore-size fractions (Rieu & 

Sposito, 1991) for total porosity and formed by broad 

range of pore sizes. Since the pores are assumed to be 

formed in a packing arrangement of the grains and the 

relation of ηi (mean pore diameter in the ith virtual pore 

size fraction) with Ωi ,then the number of pore-size classes 

can assumed to be equal to the number of particle size 

classes (“s” in Equation 16). Total porosity in terms of 

particle sizes is defined as: 

⎛ Ω ⎞

3−DF 

φ T = 1−
 ⎜ min ⎟ ⎜ ⎟ (17)

⎝ Ωmax ⎠

where DF is the fragmentation fractal dimension that can 

be determined from the slope of the best linear fit to the 

number-size relation given in Equation 1. 

In the caved zone of longwall gobs, the degree of 

compaction is a function of initial porosity and strength of 

the rock fragments. The initial porosity varies with the 

shape of rock fragments, their sizes and size distribution 

(Yavuz, 2004). He notes that the initial compaction of 

caved material will be larger during the initial loading 

(compression) stage. As compression increases, the broken 

material will stiffen due to the densely compacted material 

and a reduction in porosity in the caved zone will result by 

further crushing or fragmentation of the initial rock pieces. 

Crushing is accompanied by changes in particle size 

distribution. The crushing strength of a particle is a 

function of its size and its co-ordination number (number 

of contacts with neighbours). In the case of a large 

population of granular materials subjected to the above 

processes, the change in compaction is important for 

evaluating porosity. 

Using the principles of grain crushing and a particlesize distribution modeled as a self-similar pattern, the 

plastic reduction in porosity with a one-dimensional 

applied stress increment is given by (McDowell et al., 

1996): 

(1−DF ) 

2 ε dφ T = −βf (2
 − D )
 (1−κ )σΩ
 F 

max (18) m(1−DF ) m( DF −1) −1

× m σ 2 σ
 2 0 dσ

In this equation DF is the fragmentation fractal 

dimension defined by Equation 1, ε is the surface energy in 

linear-elastic fracture mechanics (Griffith, 1920), m is the 

Weibull modulus based on the observations that particlesurvival probability follows a Weibull distribution, σ 0 is

the tensile strength of the grains, σ is the applied 

macroscopic stress and f is fracturing probability term 

(McDowell et al., 1996). This equation can be written in a 

more compact form by separating it into the plastic 

compressibility index (Λ) and stress increment ( dσ ) 

terms: 

m(DF −1)−2

dφ T = −Λσ 2 dσ (19)

where Λ is defined as: 

(1− DF ) 

2 ε m (1− DF )

Λ = βf (2 − D )mσ 2 (20) 

(1 − κ )σΩ F 0 

max

Equation 19 can be integrated, assuming Λ is constant 

at each stress level, between two successive stress 

increments (i to i+1) to obtain (Karacan, 2009): 

⎡ m( D (D −1) −1 ⎤ σ
 F −1) m

2 i +1 σ
 F −1

φ ⎢ 2 i ⎥ Ti +1 −
 φ
Ti = −Λ
 m( 1) −
 m(D −1) (21) ⎢ DF − −1 F −1 ⎥ ⎣ 2 2 ⎦

In order to solve Equations 19-21, the plastic 

compressibility index (Λ) should be defined by either 

substituting typical values of physical parameters in 

Equation 20 or by using a typical value of 0.1 for sands 

and clays (McDowell et al., 1996). 

3.2 How to Determine Model Parameters 

In effect, Equations 19-21 evaluate the change in porosity 

as a function of increasing vertical stress and 

fragmentation fractal dimension. In each iteration, or 

incremental vertical stress, the new porosity replaces the 

porosity term (Ø) in Equation 16 to modify total flow area 

(AT). Finally, the updated flow area is inserted into 

Equation 14, for rectangular channel openings or 

Equation 15 for cylindrical channels, to estimate the 

permeability at a particular stress and fragmentation level 

(Karacan, 2009). 

The general solution procedure stated above requires 

determination of these parameters: representative length 

“Lo”, maximum and minimum particle dimensions “Ωmin 

and Ωmax” in the range(s) that shows fractal behaviour, 

Weibull modulus “m”, pore, tortuosity and fragmentation 

fractal dimensions “DP, DT and DF” respectively, virtual 
 

 

 

 

pore size classes “s”, minimum and maximum pore 

(channel) diameters “ηmin and ηmax”. 

Weibull modulus (m) is a measure of the variability in 

strength of the material. For chalk, brick, stone and 

cement, m is between 1 and 5 (McDowell et al., 1996). 

Engineering ceramics, for instance, have values for m of 

about 10 due to much less variation in strength compared 

to natural rocks or compared to composite materials with 

relatively high inter- and intra-particle porosities and 

flaws. Thus, a low m-value around 2 should be expected 

for natural materials. 

Representative length (Lo) is defined as the shortest 

distance between two points where there is a pressure 

differential that drives the flow. These two points may be 

on either side of the measurement interval within the gob 

or from a point in the gob to the gateroads. Therefore, the 

representative lengths change based on the location of the 

evaluation points in the gob. 

Fragmentation fractal dimension (DF) defined by

Turcotte (1986) showed that particle sizes of geologic 

materials exhibit fractal behaviour as shown in Equation 1. 

The cumulative number versus particle size relation

(Figure 1) that is used to determine DF can also be used to 

find the maximum and minimum particle size. Maximum 

and minimum particle diameters (Ωmin and Ωmax) are 

approximately the cut-off limits of deviation from linearity 

in the number-size relation. 

1.E+00 

1.E+01 

1.E+02 

1.E+03 

1.E+04 

1.E+05 

0.1 1

Particle Size (cm) 

N (>Omega) 

After Compaction 

Before Compaction 

Figure 1	 Log-log plots for number-particle size data 

(Equation 1) for shale test before and after 

loading. DF’s for this test - before compaction 

and after compaction are 2.19 and 2.89, 

respectively. 

The fragmentation fractal 

dimensions and maximum and minimum particle sizes 

before and after loading tests are listed in Table 2. 

10 

Table 2 Results of DF, DT, DP and maximum (Ωmax) and 

minimum (Ωmin) particle sizes (up-scaled) 

determined before and after loading tests using 

number-size relationship 

Before Loading 

s DF Ωmin 

(cm) 

Ωmax 

(cm) 

9 2.19 15 51 

Before Loading 

DT 

(straight) 

DT 

(tortuous) 

DP 

(uniform) 

DP 

(irregular) 

1 1.10 2 1.19 

After Loading 

s DF Ωmin 

(cm) 

Ωmax 

(cm) 

10 2.89 5 51 

After Loading 

DT 

(straight) 

DT 

(tortuous) 

DP 

(uniform) 

DP 

(irregular) 

1 1.26 2 1.89 

The number of virtual pore size classes “s” is equal to 

the number of different particle size fractions that are 

within the cut-off limits of the data set where fractal 

scaling is observed. Since a linear trend is observed 

between all particle size classes, the number of particle 

size classes used will be the same as the virtual pore size 

fractions (Table 2). 

Fractal dimension of the one-dimensional trace of pore 

channel (DT), or tortuosity fractal dimension, describes 

how the pores are convoluted in relation to Lo. However, it 

is not always easy to determine the tortuosity of a channel 

to find its fractal dimension. In this study, the “slit-island” 

theorem of Mandelbrot et al. (1994) was used to determine 

DT, since there is really no practical means of measuring 

this parameter. This approach has been used in the 

literature (Sammis & Steacy, 1995) to estimate the fractal 

dimension of lower order geometries. 

In order to apply the slit-island theorem to estimate the 

tortuosity fractal dimension (DT) in the crushed material, 

the fractal increment values reported by various 

researchers on different type of soil and porous rocks were 

evaluated. Jacquin & Adler (1988) observed fractal 

increments of 0.23 in pore structures of dolomitic 

limestone. Tyler & Wheatcraft (1989) reported fractal 

increments between 0.14 and 0.43, with an average of 0.23 

when obtained from soil particle size data. Karacan & 

Halleck (2003) determined that an average DT increment 

of 0.29 for incompletely fragmented porous media using 

shock waves. These increments can be used to calculate DT 

values of 1.23 and 1.29 from these references for highly 

compacted and incompletely fragmented porous media, 

respectively. The average of these two values (1.26) will 

be taken as the representative value for the compacted gob. 

For the initial gob material, 1.10 will be used as the 

average of DT. This value is based on the averages 

calculated by Wheatcraft & Tyler (1988) as 1.08 using 

Monte-Carlo simulations, values reported by Yu & Cheng 

(2002), and Yu & Liu (2004) for pores in 52% porosity 

medium between 1.10-1.12. These values have physical 

sense since compaction or incomplete fragmentation 

compared to complete fragmentation increases the 

tortuosity of pore channels. 

If the pore area represents a smooth object (circle, 

square, triangle etc.) covering the whole area, then DP can 

be taken as 2. If the pores are irregular in shape and 

distribution, then their area-fractal dimension will be less 
than 2 (Ahmed & Dryzmala, 2005). In case of irregular 

pore shapes and distributions, slit-island theorem can be 

used by employing the DF values (Table 2) obtained in this 

study. Since particle size distribution is a 3-dimensional 

entity (Tyler & Wheatcraft, 1989), its Euclidean dimension 

is 3. Then, the calculated increments (DF-3) are added to 

the Euclidean dimension for area to obtain DP (Table 2). 

Tortuosity (DT) and pore-area (DP) fractal dimensions 

for straight and regular pores as well as the values 

determined using the above approaches for the convoluted 

and irregular pores are given in Table 2. These values will 

be used in model calculations later in this paper. 

In order to calculate porosity and permeability using 

actual gob-material sizes and their distributions, the 

gradation curves of laboratory-simulated gob material 

were up-scaled by 10 times to their original sizes. Upscaling, in effect, changes particle sizes and eventually the 

channel sizes that will be calculated based on particle size. 

However, as long as the gradation curve is held constant, 

up-scaling does not change the fragmentation fractal 

dimensions since this will not change the slope. Thus, the 

fractal dimensions reported in Table 2 will still be valid. 

The size of the grains in a fragmented porous media 

will directly affect the pore flow area and pore diameters 

forming between the grains. In order to calculate 

permeability under increasing loading, up-scaled grain 

sizes will be used in area- and pore-size calculations. In 

order to incorporate the effect of loading on permeability, 

a fractal crushing model for porosity will be integrated into 

calculations. 

Equation 16 is used to calculate the total flow area in 

the fragmented gob material. This equation can 

alternatively be represented with an integral: 

π Ω max 

∫Ω φ N1−DT dΩ (22) 3 Ω min 

In this integral, maximum and minimum particle sizes 

are the up-scaled values and the number (N) is a function 

of grain size. In order to find this function, the numbers of 

up-scaled particles for each up-scaled size were plotted 

(Figure 2). 

 

0 

10 

20 

30 

40 

50 

60 

70 

80 

90 

0 20 40 60 80 100 

Grain diamter (cm) 

Number of Particles 

Figure 2 Number-size distributions for the up-scaled 

diameters for the initial gob material. 

By up-scaling the material mass (assuming 

density stays constant) by 1000 times due to a 10 fold 

increase in particle diameter, the functional relationship 

given in Equation 23. Figure 2 shows that as diameters 

increase, the number of particles decreases, as expected. 

More importantly, the figure shows that the data can be 

successfully represented by power functions. 

In Equation 16, Ø-values are the partial porosities for 

each infinitesimal size interval. Integration between Ωmin 

and Ωmax will result in the total fractal porosity (ØT), which 

will not be dependent on Ω. With these treatments, the 

fractal flow area (Equation 16) becomes: 

π
 Ω max

φ
T ∫ (10326Ω 1− −1 ) .74 DTΩ dΩ (23) 3
 Ωmin 

After rearrangement, the integral shown below is obtained: 

10326 1−D Ω T π φ T ∫

max

Ω −0.74+1.74DT dΩ (24) 3 Ωmin 

Solving the integral gives the fractal pore area for 

initial gob material as: 

⎡ Ω 0.26+1.74DT ⎤ max − 10326 1−D ⎢ ⎥ T π
 ⎢1.74DT + 0.26 ⎥ AT = AP = φ
T (25) 

3 ⎢ 0.26+1.74DT ⎥ Ω⎢ min ⎥

⎣⎢1.74DT + 0.26 ⎦⎥

This equation can directly be used in Equations 14 or 

15 for “area“once the fractal bounds of particle sizes and 

the size of the flow channels are determined. In Equation 

25, ØT is the porosity which changes as uniaxial stress 

progressively increases according to Equations 19-21 

(Karacan, 2009). 

In order to determine maximum and minimum channel 

sizes, Equation 26 given below is used. 

⎛ 1/ 2 2

η
 = 1 DT

i i ⎜ φi N − ⎞ Ω i ⎟ (26) ⎝ 3 ⎠

This equation gives the sizes of flow channels in each 

partial volume of a fractal porous medium. In this 

equation, the number term (N) is replaced by its functional 

form given in the integrand in parenthesis of Equation 23. 

The partial porosity term in the equation is, on the other 

hand, the amount of porosity in that partial volume (Rieu 

& Sposito, 1991). In order to obtain partial porosity, the 

total porosity at a particular loading stage can be divided 

by the “virtual pore size fractions” or the number of 

particle size classes “s” given in Table 2. This value was 

around 9-10 in this study as tabulated in Table 2. With 

these additions, the maximum and minimum channel 

diameters (ηmax and ηmin) as a function of Ω can be 

calculated. These relations are given in Table 3. 
 

 

Table 3 Minimum and maximum sizes (ηmin and ηmax) of 

tortuous flow channels as a function of particle 

size. Particle diameters to be used are the upscaled values. The units of both particle sizes and 

the resulting channel sizes are in cm. 

ηmax 

0.82 T 

s 

φ ( ) D T DT 0.26 1.74 10236 max

1 − + Ω 

ηmin 

0.82 T 

s 

φ ( ) D T DT 0.26 1.74 10236 min

1 − + Ω 

4 Model Application to Predict Gob Porosity and 

Permeability under Loading 

4.1 Porosity of Initial and Compacted Gob Material 

In order to compare the initial porosities given by the 

experiment at 0.0 MPa load (Table 1) with the initial 

porosities calculated using particle size distribution and 

fractal model equations, Equation 17 was used first. This 

equation requires the upper and lower bounds (in terms of 

particle size) of distribution across which a fractal scaling 

is observed and the calculated fragmentation fractal 

dimension. Since the size distributions given in Figure 1 

reveal reasonably good linear relations for the initial gob 

material, the upper- and lower-size limits were selected to 

capture the whole particle-size range in each test (Table 2). 

In addition, the DF’s given in Table 2 were used for the 

same cases. 

The porosity values obtained using Equation 17 and the 

values determined in the laboratory before loading are 0.63 

and 0.8, respectively. Before calculating the porosities 

using particle size distributions and the theory described in 

previous sections, porosities that should be experienced in 

the gob were calculated using the experimental stressstrain curve (Pappas & Mark, 1993) and the initial porosity 

determined from Equation 17 (0.63). The results of this 

calculation are shown in Figure 3. 

0.00 

0.10 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 

0.80 

0.90 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

Porosity 

Figure 3	 Porosity predicted for compacting gob materials 

using Equations 17-21. 

The porosity evaluation during uniaxial compression of 

simulated gob material shows that at earlier loadings, the 

material could be compressed more when compared to 

later stages. At around a stress level of 19 MPa (2750 psi), 

the predicted final porosity of the material was about 20­

22%. A difference of 4-5% between the values given in 

Figure 3 and the values in Table 1 can be considered close 

for this analysis. As an inner check point, the porosity 

values reported in Table for 5.4 MPa (800 psi) were 

compared with the corresponding data shown in Figure 3. 

These values are also within 2-3% of each other. 

As the longwall face advances during mining and the 

overburden starts settling on the caved material, the 

initially caved gob material is subjected to progressively 

increasing stresses. These stresses decrease the initially 

observed high porosities by crushing the gob material into 

smaller pieces as shown in the previous sections 

(Figure 3). 

The porosity of simulated gob material after laboratory 

loading test was calculated by using the particle size 

distributions data after the test (Figure 1) and the fractal 

model given in Equation 17. In summary, a similar 

procedure that was followed to calculate initial porosities 

was applied. The porosities were calculated using the 

fragmentation fractal dimensions given in Table 2 within 

the linear range of particle size data (Figure 1) where 

fractal behaviour is observed. For these calculations, the 

lower bound of the particle diameter after loading was 

determined as 5 cm (for up-scaled diameters). In this 

calculation, a porosity value of 0.23 was predicted, 

compared to 0.16 that was obtained during laboratory tests. 

Again, this difference between predictions and 

reported/calculated values can be regarded as acceptable 

predictions for practical purposes. 

The porosity calculated using fractal methods and the 

number-size relations discussed so far and their 

comparison with laboratory measured or calculated data 

shows that the gob material forms a fractal porous 

medium. The porosity of this medium, before or after 

being subjected to a vertical stress, can be calculated 

within an acceptable error using number-size distribution 

data available for the fragmented gob material. The 

evolution of the porosity during the entire loading period 

can also be estimated by using Equations 17-21, as 

presented in Figure 3. This suggests that the fractal 

crushing equation can be used to predict evolving 

porosities in a gob as a result of compression. One can use 

the initial number-size distribution of a granular media and 

its fragmentation fractal dimension to calculate changing 

porosities as stress is progressively increased. 

4.2	 Calculation of Permeability of Initial and 

Compacted Gob Material 

In order to calculate gob permeability, Equations 14 or 15 

can be used depending on the expected shape of the flow 

channels. Certainly, this is a hard judgment to make for a 
complex fragmented medium. However, in any case, in 

this study the channels were taken as rectangular shaped 

pores, channels, and even as slit-like rectangular 

geometries. Mathematically, this transition in the shape 

can be created by changing the aspect ratio (Equations 9­

10) between the two-axes of the rectangular geometry for 

use in the H-P type flow and permeability prediction. 

In order to calculate permeability of a gob composed of 

broken rocks and rectangular flow channels, Equation 14 

was used. In this equation, total pore area for flow was 

calculated using Equation 25 and the flow channel sizes 

were taken as in Table 3. As it can be noticed, the flow 

area is a function of particle size and tortuosity and area 

fractal dimensions (DT and DP). Flow channel sizes are 

also a function of DT. In the permeability calculations, the 

DT=1.26 was used as tortuosity of the rectangular 

channels, whose areas were irregular (DP=1.19). The upscaled maximum and minimum rock sizes (Ωmax and Ωmin) 

for initial gob material were based on Table 2. 

Table 4 shows the aspect ratios used for the rectangular 

channels and the calculated shape factors (Equation 9).

Table 4 Aspect ratios used and the calculated shape 

factors for rectangular flow channels. 

Aspect 

Ratio 

Shape Factor 

10 149.32 

20 266.13 

100 1224.04 

250 3023.75 

750 9023.63 

5000 60023.57 

 As 

this table indicates as the aspect ratio increases, i.e. as the 

rectangular opening becomes more of a slit-like structure, 

the shape factor increases. Also, since shape factor (α) is in 

the denominator of Equation 14, the value of permeability 

is expected to decrease. This is a situation that is expected 

in reality too. 

Figures 4 and 5 show the predicted permeability of 

rectangular channels of various aspect ratios in gob as a 

function of uniaxial compressive load.

0.00E+00 

1.00E-05 

2.00E-05 

3.00E-05 

4.00E-05 

5.00E-05 

6.00E-05 

7.00E-05 

8.00E-05 

9.00E-05 

1.00E-04 

1.10E-04 

1.20E-04 

0 5 10 15 20

Uniaxial Loading (MPa) 

Permeability (cm 2 ) 

10:1 

20:1 

100:1 

Figure 4 Calculated permeability change for a shale-rich 

gob, where porous space for gas migration is 

formed with irregular rectangular channels of 

geometrical aspect ratios of 10:1, 20:1, 100:1. 

0.00E+00 

5.00E-07 

1.00E-06 

1.50E-06 

2.00E-06 

2.50E-06 

3.00E-06 

3.50E-06 

4.00E-06 

4.50E-06 

5.00E-06 

5.50E-06 

6.00E-06 

0 5 10 15 20 

Unia xi al Loa di ng (MPa ) 

Perm eability (cm 2 ) 

250: 1 

750: 1 

5000: 1 

Figure 5 Calculated permeability change for a shale-rich 

gob, where porous space for gas migration is 

formed with irregular rectangular channels of 

geometrical aspect ratios of 250:1, 750:1, 

5000:1. 

 Both figures show 

that the permeability is higher initially and decreases as a 

result of crushing and compaction during loading. The 

predicted permeabilities are changing from 0.00011 cm2 

(11×103

 Darcy) for a gob initially composed of bigger 

rocks forming irregular, meandering tortuous channels of 

rectangular shape and 10:1 aspect ratio to 0.00004 cm2 

(4×103

 Darcy) for the same gob under compressional 

loading. When the initial aspect ratio is 20:1, the initial 

permeability is around 0.00006 cm2

 (6×103 Darcy), which 

decreases to 0.000025 cm2

 (2.5×103 Darcy) when the gob 

is compressed (Figure 4). As it can be seen from Figures 4 

& 5, as the aspect ratio increases the initial and final 

permeabilities decrease for each of the six cases modeled 

for this paper. 

The obvious question, of course, is how these 

permeabilities compare with the measurements.

Unfortunately, there are not any field measurements of 

permeability of gobs due to their inaccessibility for 

performing direct measurements. However, there are some 

reported values calculated based on volumetric strains and 

geomechanical calculations. For instance, Brunner (1985) 

used gob permeability values in ventilation network 

models as 1×10-7 m2

 to 1×10-5 m2 (1×105 to 1×107 Darcy). 

Ren et al. (1997) estimated the permeability in the gob to 

be on the order of 1×10-10 m2

 (1×102 Darcy) in the 

compacted region. Wendt & Balusu (2002) used maximum 

values around 1×10-9 m2

 (1×103 Darcy). Whittles et al. 

(2006) reported calculated values in the range of 5×10-7 m2 

(5×105

 Darcy) to 1×10-8 m2 (1×104 Darcy) for a gob. 

Esterhuizen & Karacan (2007) reported calculated values 

of 1×10-9 m2

 (1×103 Darcy). The differences in reported 

permeabilities may be related to different coal-seam 

geologies, different panel layouts and to differences in 

caving characteristics. Nevertheless, it should be noted that 

all these reported values are within similar orders of 

magnitude and close to the permeability values calculated 

by the fractal flow models presented in this study. Thus, 

the fractal-fragmentation and fractal-crushing driven 

methodologies and the flow models presented in these 
analyses can predict permeabilities close to the values 

reported in other studies. 

5 Summary and Conclusions 

This study developed a predictive approach using the 

principles of scaling and fractal porous medium combined 

with fluid flow. The fractal approach constructed flow 

equations and fractal crushing equations for granular 

materials to predict porosity and permeability for a 

completely fragmented porous medium. The virtual 

fragmented fractal porous medium so generated was 

exposed to various uniaxial stresses to simulate gob 

compaction and porosity and permeability changes during 

this process for various flow channel complexities and 

irregularities. 

Particle size distributions of gob material determined 

from digital pictures or from simulated tests can be used to 

find the fragmentation fractal dimension of gob material 

based on a number-size relationship, which then can be 

used to find porosity. Fluid flow in a fractal porous media, 

on the other hand, can be described using fractal scaling 

between pore and grain sizes. These information can be 

used in the Hagen-Poiseuille (H-P) equation to describe 

flow rate, which then can be integrated with the Darcy 

equation to describe permeability of a fractal porous media 

in terms of textural and fractal properties. Moreover, 

fractal crushing of granular materials can be used in the 

definition of flow area to simulate compaction of a gob. In 

this study, the porosities predicted using a fractal crushing 

model were found to be in good agreement with the 

porosity values calculated using stress-strain data and 

particle size data. 

Geometrically, circular and smooth channels should 

give the highest permeability. This study was conducted on 

rectangular-shaped geometries as a more realistic approach 

to channel shapes in the gob. Based on this approach, it 

was shown that as the irregularity and aspect ratio of the 

channels increase, permeabilities decrease as expected. 

This paper demonstrated a new approach to calculate 

porosity and permeability of gob. Using the approach 

developed in this paper, it is possible to make predictions 

of gob porosity and permeability. This will lead to 

predictions regarding the flow amount and flow patterns in 

the gob. The ability of improving predictions in the flow 

amounts and flow patterns in the gob using this simple 

technique will lead to performing leakage calculations and 

methane control projections. However, one should exercise 

caution in treating the data from either digital images or 

from particle size experiments as the technique is closely 

related to geometrical and scaling parameters. Therefore, 

results can be sensitive to errors in data analyses, which 

can lead to spurious porosity and permeability values. 

Nomenclature 

AG area of grains 

AP area of pores 

AT total area 

DF fragmentation fractal dimension 

DP pore area fractal dimension 

DT tortousity fractal dimension 

f fracturing probability term 

h thickness of the analysis location 

k permeability 

Lt tortuous length 

Lo representative length 

s number of virtual pore size fractions 

Ni number of particles in the ith size class 

ΔP differential pressure 

m Weibull modulus 

Q total flow rate 

V volume of porous medium 

Symbols 

α shape factor in H-P equation 

η pore diameter 

ηi pore diameter of the ith size class 

φ porosity of fragmented porous medium 

φi porosity in partial volume in fractal 

porous media 

φT total porosity 

μ viscosity 

Ωi particle diameter of the ith size class 

Ω average particle diameter 

Π hydraulic resistance 

Ψ arbitrary cross sectional area 

Λ plastic compressibility index 

v(x, y) fluid velocity perpendicular to x-y plane 

ε surface energy in linear-elastic fracture 

mechanics 

σ 0 tensile strength of the grains 

σ applied macroscopic stress 

Permeability conversion factor 

1 Darcy 10-12 m2

 = 10-8 cm2 

References 

Ahmed, A.M.H. & Drzymala, J. (2005). Two-dimensional fractal 

linearization of distributed curves. Physicochemical 

Problems of Mineral Processing 39, 129-139. 

Arya, L.M., Paris, J.F., (1981). A physicoemprical model to 

predict the soil moisture characteristics from particle-size 

distribution and bulk density data. Soil. Sci. Soc. Am. 

Journal, 45, 1023-1030. 

Brunner, D.J. (1985). Ventilation models for longwall leakage 

simulation. Proc. 2nd US Mine Ventilation Symposium, 23­

25 September, Reno, NV. 

Esterhuizen, G.S. & Karacan, C.Ö. (2005). Development of 

numerical models to investigate permeability changes and 

gas emission around longwall mining panels. The 

AlaskaRocks 2005, 40th US Symposium on Rock 

Mechanics, Anchorage, Alaska, 25-26 June. 

Esterhuizen G.S. & Karacan C.Ö. (2007). A methodology for 

determining gob permeability distributions and its 

application to reservoir modeling of coal mine longwalls. 

2007 SME Annual Meeting, Denver, CO. 
Griffith, A.A. (1920). The phenomena of rupture and flow in 

solids. Philos. Trans. R. Society of London A221, 163-198. 

Jacquin, C. G. & Adler, P. M. (1988). Fractal porous media II: 

geometry of porous geologic structures. Transport in Porous 

Media, 2, 571-596. 

Karacan, C.Ö. (2009). Prediction of porosity and permeability of 

caved zone in longwall gobs. Transport in Porous Media (in 

press).

Karacan, C.Ö. & Halleck, P.M. (2003). A fractal model for 

predicting permeability around perforation tunnels using 

size distribution of fragmented grains. J. Petroleum Sci. and 

Eng. 40, 159-176. 

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. 

Mandelbrot, B.B., Passaja, D.E. & Paulley, A.J. (1994). Fractal 

character of fracture surfaces of metals. Nature, 21, 308. 

McDowell, G.R., Bolton, M.D. & Robertson, D. (1996). The 

fractal crushing of granular materials. Journal of the 

Mechanics and Physics of Solids 44, 2079-2101. 

Mortensen, N.A., Okkels, F. & Bruss, H. (2005). Reexamination 

of Hagen-Poiseuille flow: Shape dependence of the 

hydraulic resistance in microchannels. Physical Review E 

71, 057301, 1-4. 

Pappas, D.M., & Mark, C. (1993). Behavior of simulated 

longwall gob material. Report of Investigations RI No. 

9458, US Dept. of Interior, US Bureau of Mines. 

Perfect, E. & Kay, B.D. (1995). Applications of fractals in soil 

and tillage research: a review. Soil Tillage Res. 36, 1-20.
 

Perfect, E. (1997). Fractal models for the fragmentation of rocks 

and soils: a review. Engineering Geology, 48, 185-198. 

Ren, T.X., Edwards, J.S. & Josefowicz, R.R. (1997). CFD 

modeling of methane flow around longwall coal faces. Proc 

6th International Mine Ventilation Congress, 17-22 May, 

Pittsburgh, PA. 

Rieu M. & Sposito, G. (1991). Fractal fragmentation, soil 

porosity and soil water properties: theory, Soil. Sci. Soc. 

Am. J., 55, 1231-1238. 

Sammis, C.G. & Steacy, S.J. (1995). Fractal fragmentation in 

crustal shear zones. In: Fractals in the Earth Sciences. 

Barton C., La Pointe, P.R. (Eds.), Pleanum Publishing Co., 

New York, NY. 265 pp. 

Turcotte, D.L. (1986). Fractals and fragmentation. J. Geophysical 

Res. 91, 1921-1926. 

Weiss, J. (2001). Fracture and fragmentation of ice: a fractal 

analysis of scale invariance. Engineering Fracture 

Mechanics, 68, 1975-2012. 

Tyler, S.W. & Wheatcraft, S.W. (1989). Application of fractal 

mathematics to soil water retention estimation. Soil. Sci. 

Soc. Am. J., 53, 987-996. 

Wendt, M. & Balusu, R. (2002). CFD modeling of longwall goaf 

gas flow dynamics. Coal and Safety 20, 17-34. 

Wheatcraft, S.W. & Tyler, S.W. (1988). An explanation to scaledependent dispersivity in heterogeneous aquifers using 

concepts of fractal geometry. Water Res. Research, 24, 566­

578. 

Whittles, D.N., Lowndes, I.S., Kingman, S.W., Yates, C., & 

Jobling, S. (2006). Influence of geotechnical factors on gas 

flow experienced in a UK longwall coal mine panel. 

International Journal of Rock Mechanics and Mining 

Sciences 43, 369-387. 

Yavuz, H. (2004). An estimation of cover pressure reestablishment distance and pressure distribution in the goaf 

longwall coal mines. International J. of Rock Mech. and 

Mining Sciences 41, 193-205. 

Yu, B.M. & Cheng, P. (2002). A fractal model for bi-dispersed 

porous media. Int. J. of Heat and Mass Transfer, 45 (11), 

2983-2993. 

Yu, B.M. & Liu, W. (2004). Fractal analysis of permeabilities for 

porous media. AIChE Journal 50(1), 46-57. 

Yuan, L. & Smith, A.C., (2008). Effects of ventilation and gob 

characteristics on spontaneous heating in longwall gob 

areas. Proc. 12th U.S./North American Mine Ventilation 

Symposium, pp. 141-147, Reno, NV, June 9-11.