Prediction of Porosity and Permeability of Caved Zone

in Longwall Gobs

C. Özgen Karacan

C. Ö. Karacan

NIOSH, Pittsburgh Research Laboratory, Pittsburgh, PA 15236, USA

e-mail: cok6@cdc.gov

Abstract The porosity and permeability of the caved zone (gob) 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 its importance, the gob is

typically inaccessible for performing directmeasurements 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 predictive approach

that combines fractal scaling in porous medium with principles of fluid flow. The approach

allows the calculation of porosity and permeability from the size distribution of broken rock

material in the gob, which can be determined from image analyzes of gob material using the

theories on 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. The results suggest that the gob

porosity and permeability values can be predicted by this approach and the presented models

are capable to produce values close to values documented by other researchers.

List of Symbols

AG Area of grains, L2
 

A Area of pores, L2
 P

AT Total area, L2
 

DF Fragmentation fractal dimension (fractal dimension)
 

DP Pore-area fractal dimension (fractal dimension)
 

DT Tortuosity fractal dimension (fractal dimension)
 
f Fracturing probability term (dimensionless)

h Thickness of the analysis location, L

k Permeability, L2

Lt Tortuous length, L

Lo Representative length, L

s Number of virtual pore-size fractions, number, N

Ni Number of particles in the ith size class, number, N

IP Differential pressure, P

r Similarity coefficient (dimensionless)

m Weibull modulus (dimensionless)

Q Total flow rate, L3/t

V Volume of porous medium, L3

α Shape factor in H–P equation

η Pore diameter, L

ηi Pore diameter of the ith size class, L

Porosity of fragmented porous medium, L3 φ L−3

p b

φi Porosity in partial volume in fractal porous media, L3L−3

p b

L−3 φT Total porosity, L3

p b , dimensionless

r Pore coefficient, dimensionless

μ Viscosity, ML−1t

−1

Qi Particle diameter of the ith size class, L

Q ̄ Average particle diameter, L

I Hydraulic resistance, L−2

' Arbitrary cross-sectional area, L2

A Plastic compressibility index (dimensionless)

ν(x, y) Fluid velocity perpendicular to x–y-plane, L/t

ε Surface energy in linear-elastic fracture mechanics, E

σ0 Tensile strength of the grains, P

σ ̄ Applied macroscopic stress, P

Permeability conversion factor 1 Darcy=10−12m2 = 10−8cm2

1 Introduction

Longwall mining is an underground mining method that can maximize coal production in

coal beds, which contain few geological discontinuities. In these operations, a mechanical

shearer progressively mines a large block of coal, called a panel, which is outlined with development entries or gate roads. Asthe coal is extracted, the supports automatically advance and

the roof strata allowed to cave behind the supports (Karacan et al. 2007) creating large-scale

disturbances of the surrounding rock mass due to fracturing and caving of mine roof behind

the shields. The caved zone (gob) created by longwall mining could reach 4–11 times the

thickness of the mining height where overburden rocks are weak and porous. The gob material is highly fragmented and is 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 within the zone or flowing from surrounding formations into the

mining environment. Figure 1 shows example pictures from the caved zones of two Eastern

Kentucky coal mines.
Fig. 1 Gob materials photographed in two Eastern Kentucky coal mines (from Pappas and Mark 1993)

As mining progresses, the gob gradually consolidates sufficiently to support large loads

resulting from the overburden weight (Pappas and 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

gob gas ventholes. Thus, an understanding of compaction phenomena and the resultant reservoir properties of gob material is very important for developing adequate methane control

strategies.

Accuracy of models for predicting methane flow in the vicinity of the longwall face

and methane control in longwall mines depends on an accurate representation of the gob

porosity and permeability. Despite the importance of these two reservoir properties, their

predictions are difficult have been conducted by only a few studies (Brunner 1985; Ren

et al. 1997; Wendt and Balusu 2002; Esterhuizen and Karacan 2005, 2007; Whittles et al.

2006; Karacan et al. 2007). This may be due to challenges and unknowns related to the

gob environment. These challenges are even more difficult due to the inaccessibility of the

gob environment for conducting direct measurements of porosity and permeability using

conventional tests.

This work presents a fractal porous medium model to predict gob porosity and permeability for controlling methane flow and to predict air leakage into the gob. For this approach,

the particle size data for simulated gob material given by Pappas and Mark (1993) were used

to calculate fragmentation fractal dimensions and particle size distributions before and after

controlled loading tests. The size distributions were scaled up to real dimensions to describe

the geometry and size of flow channels. 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 and 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 due to loading.

2 Experimental Methods to Obtain Gob Material Size Distribution

In this paper, rock fragment size distribution data were taken from Pappas and Mark (1993).

They conducted a laboratory study to estimate the gradation of actual gob material, to

evaluate the stress–strain behavior of simulated gob material using load deformation tests. In

this paper, the particle size gradation data were further analyzed with the developed models
to predict porosity and permeability. This section gives a brief description of the laboratory

tests conducted by Pappas and Mark (1993).

2.1 Obtaining Rock Particle Gradation Curves for Tests

Pappas and Mark (1993) conducted laboratory tests on rock materials with properties similar

to those of actual gob material. The characteristics 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. They reduced the rock size and gradation to a laboratory scale by proportionally shifting

the particle size distribution curve of the actual gob material determined from pictures taken

from the headgate entries (Fig. 1) where portions of the gob could be viewed. This shift was

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

particle size of the laboratory sample was reached. The detailed procedure for this can be

found in Pappas and Mark (1993).

2.2 Load Deformation Tests on Fragmented Rocks

For load deformation tests on fragmented rock material, a test chamber was built in which

the simulated gob material was subjected to different loads. The weight of the rock in the test

chamber was determined by subtracting the weight of the chamber from the total weight. This

value was used to calculate initial void volume, or porosity, of the simulated gob material.

Then, the loading on the rock material wasincreased at a specified rate and load and displacement were monitored. The maximum load applied to the simulated gob material increased

as the test series progressed.

Pappas and Mark (1993) conducted 20 tests by varying the maximum particle size in the

simulated gob material and the maximum applied stress. In this paper, the size distributions

obtained with maximum particle sizes of 5.1 cm (2 inches) for shale and weak sandstone

were used for modeling and comparison purposes. These tests will be referred to as Test-A

(shale-maximum size 5.1 cm), Test-B (shale-maximum size 8.9 cm), Test-C (weak sandstone-maximum size 5.1 cm), and Test-D (weak sandstone-maximum size 8.9 cm). Table 1

gives the test summary of each simulated gob material. Figures 2 and 3 show particle size

distributions of shale and weak sandstone before and after loading tests.

Table 1 Summary of simulated materials (from Pappas and Mark 1993) analyzed in this paper and their

loading tests

Test Rock type Maximum size Maximum Initial Porosity @ Porosity @

(cm) stress (MPa) porosity 5.4MPa max. stress

A Shale 5.1 19.0 0.802 0.365 0.160

B Shale 8.9 21.3 0.679 0.269 0.078

C Weak sst 5.1 19.3 0.790 0.313 0.152

D Weak sst 8.9 21.1 0.719 0.337 0.162

“sst” in tests C and D stands for “sandstone”

The data are for four

tests conducted on two different initialsize distributions with different initial maximum sizes.

Figure 4 shows the stress–strain data obtained during the loading tests. The data

showed that the stress–strain behavior of simulated materials was non-linear, indicating a
strain-hardening process during compression. Again, it was observed that the stress–strain

behavior was not dramatically affected by the rock type. However, it was observed that as

the maximum size of initial material increased from 5.1 cm to 8.9 cm, the granular medium

experienced somewhat less compaction at higher stresses.

In the following sections, particle size distribution data of simulated gob material before

and after loading tests (Figs. 2 and 3), as well as stress–strain data (Fig. 4), will be used to

calculate various fractal dimensions of flow channels and rock skeleton. These techniques

will develop a method for predicting porosity and permeability of a “fractal” gob during

compaction from overburden stress.

0 

20 

40 

60 

80 

100 

120 

% Finer by Weight 

5.1 cm max. diameter-initial distribution 

Shale-5.1 cm max. diameter-after 19 MPa load 

Weak sst-5.1 cm max. diameter-after 19.3 MPa load 

100 10 1 0.1 

Grain Size (mm) 

Fig. 2 Particle size distribution ofshale before and after compressing the simulated gob material in the loading

chamber

0 

20 

40 

60 

80 

100 

120 

% Finer by Weight 

8.9 cm max. dimater-initial distribution 

Shale-8.9 cm mx. diameter-after 21.3 MPa load 

Weak sst-8.9 cm max. diameter-after 21.1 MPa load 

100 10 1 0.1 

Grain Size (mm) 

Fig. 3 Particle size distribution of weak sandstone before and after compressing the simulated gob material

in the loading chamber
0.00 

2.00 

4.00 

6.00 

8.00 

10.00 

12.00 

14.00 

16.00 

18.00 

20.00 

Vertical Stress (MPa) 

Shale-5.1 cm max. diameter 

Shale-8.9 cm max. diameter 

Weak sst-5.1 cm max. diameter 

Weak sst-8.9 cm max. diameter 

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

Strain (%) 

Fig. 4 Stress–strain curves for shale and weak sandstone of different initial particle distributions (Figs. 2

and 3) subjected to various vertical loadings in the test chamber

3 Theory and Model Development

3.1 Fractal Fragmentation

Fragmentation is a 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.

Fragment size distribution and fragment density may influence the mechanics of fluids

in the fragmented medium. However, an explicit formulation of the relationship between

pore texture and the hydraulic properties of the porous media remains a challenge due to the

complex pore-particle geometry. Yu and Liu (2004) presented fractal analyzes of permeabilities for both saturated and unsaturated porous media based on the fractal nature of pores.

Schlueter et al. (1997) defined the pores of sedimentary rocks from thin sections and scanning electron microscope (SEM) images using a perimeter–area relationship. They discussed

how the pore-area fractal dimensions and the grid sizes used to measure fractal dimensions

can affect permeability predictions. Hydraulic properties in a fragmented porous medium

can also be determined either from simple sieve analyzes or from a laser-diffraction analysis of the particles (Karacan and Halleck 2003). Bitelli et al. (1999) and Perfect and Kay

(1995) used particle size data to calculate the fractal dimensions of different soil samples

exposed to various operational and weathering conditions using the power law scaling of the

number–particle size relation given by Turcotte (1986):

N D (ω > Q) = B F FQ− , (1)

where N is the cumulative number of particles of size ω greater than a characteristic size

Q, 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 aggregates, and

BF is a coefficient related to number of particles of unit diameter.

The first attempt to directly relate the particle size distribution to a soil water characteristic was made by Arya and Paris (1981). In their model, pore lengths (based on spherical

particles) were scaled to natural pore lengths using a scaling parameter α with an average
value of 1.38. This fitting parameter α was derived from least squares regression of the predicted water content to the measured water content and its use was justified by accounting

for the non-spherical nature of the particles. Tyler and Wheatcraft (1989) interpreted α as the

fractal dimension of a tortuous pore. Since, then, there has been growing interest in the use

of fractals to predict hydraulic properties from particle size distributions (Rieu and Sposito

1991; Tyler and Wheatcraft 1992; Comegna et al. 1988).

3.2 Model Development for Predicting Permeability in a Fragmented Porous Medium

Under Successive Compaction

3.2.1 Flow in Fractal Porous Media

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 and 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. For a straight capillary, the tortuous

length equals the representative length. Following Wheatcraft and Tyler (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 and Cheng 2002):

L 1 t(η) = D D η − T L T o . (2)

In this equation, DT is the tortuosity fractal dimension, which can take values between 1

and 2. DT represents the convolutedness 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 and Tyler 1988). However, Wheatcraft and 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:

 D

η P max N(L ≥ η) = . (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 −(D η η P 1) P max + d(η). (4)

In Eq. 4, DP is the pore-area fractal dimension (1 < DP < 2). A DP-value of 2 represents

a regular pore area (circle, square etc.) 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 Eqs. 3 and 4 as:
= max D

≥

 η P

Nt(L ηmin) . (5)

ηmin

Dividing Eq. 4 by Eq. 5 obtains an expression for the probability density function:

dN − Nt

= D D η P −(DP+1) P η d(η) = f (η)d(η), min (6)

where

f (η) = D D P −(DP Pη η min

+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-sectionalshape, with no-slip boundary conditions, the Hagen–Poiseuille (H–P) flow equation can be used. For these flow conditions, the

Navier-Stokes equation takes the form:

2 2 IP

(∂x + ∂y )ν(x, y) = . (8)

μL

H–P flow is often characterized by the hydraulic resistance,Ihyd = IP/Q, in which IPis

the pressure drop along the channel and Q is the flow rate. By including the drag and the

effect of surface area, hydraulic resistance can alternatively be given as Ihyd

∗ = 2 μL/A .

In this equation, μ is the dynamic viscosity of the fluid and L is the channel length. For an

irregular shape on the x–y plane, the cross-sectional area “A” is defined as dxdy. Using '

the same approach,

 

 dl is the perimeter of the object. Then, the hydraulic radius can be 'p

defined as: dxdy R '

H = (9) dl 'p

Figure 5 shows an arbitrary cross section ' with

 

perimeter 'p in the x- and y-plane for a

straight channel placed along the z-axis. In H–P flow, there is a no-slip boundary condition

at the walls of 'p, and thus the actual hydraulic resistance will depend on the perimeter

as well as the cross-sectional area. This dependence can therefore be characterized by a

dimensionless correction factor α for shape:

Ihyd α = . (10) I∗

hyd

The relation between the pressure drop IP, the velocity ν(x, y), and the geometrical correction factor becomes:

IP = IhydQ = αI∗

hydQ = αI ν( , ∗hyd

'

 

dxdy x y). (11)

If this equation is solved for a circular area (circular pipe) of ', we obtain:

πIP c η Qi = , (12) 16αμLt

where η is the hydraulic diameter and Lt is the total length of the channel. In circular pipes,

α = 8 and the exponent “c” of the equivalent channel diameter is 4. For some known shapes,
such as ellipse, rectangle, and triangle, the shape correction factor (α) is given by Mortensen

et al. (2005):

For ellipse:α = 4π

 a b

+

 

, (13) b a

where a is the half-length of the major axis and b is the half-length of the minor axis.

⎛ −1

3 π γ 2

 ∞ γ 2 nπγ For rectangular cross sections:α = ⎝ n4 − tanh , (14) 8 5n π π 5

3

(

2 n=1, ,5...

⎞

)

where γ is the width-to-height (w/h) ratio of the rectangular channel (Mortensen et

⎠

al. 2005).

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, Eq. 14 becomes:

12 5 π γ 2

α = , (15)

π γ5 − 186ζ (5)

where ξ (5) results from truncating the summation after fifth iteration in Eq. 14.

For a triangular crosssection: α √ = 20 3, if it is an equilateral triangle. If the crosssection

is another triangular shape, then:

√ 25 40 3 α = C + , (16) 17 17

where C = 8(a

1

+b+c 2)

J . In this equation, a, b, and c are the side lengths of the

b 2 (a2+ 2+c2 4 4 4 ) −(a 2 +b +c )

triangular channel.

Arya et al. (1999) attempted to determine shape factor “α” and exponent “c” empirically

using experimental flow data in different soils. As a result of their analysis performed on

different soil textures and materials (clays, sands, and sandy loams), they found that the

exponent of diameter “c”, in Eq.12, varied between 2.808 and 4.714. The average of all

tests indicated that this number is 3.531. For the shape correction factor α, however, they

determined that the logarithm of α for different textures changed from 2.759 to 6.742, with

an average of 5.507 (Arya et al. 1999). In this paper, the average diameter exponent, from

Arya et al.’s study (3.531), will be used in the H–P equation (Eq. 12) to describe flow in

single channels. The variable “α”, will be retained in the equation to run parametric analyzes

related to channel shape.

x
 

Ψp

Ψ

y

z

Fig. 5 An arbitrary cross section of a pore with cross sectional area ' and perimeter 'p
 

 

 

 

 

η3.531 π IP

Qi(η) = (17) 16 Lt(η)α μ

If the total flow rate is calculated for each channel between the minimum and maximum sizes

over the entire range of pore sizes, Eqs. 4 and 17 can be integrated to give the total flow rate:

ηmax

QT = − Qi(η)dN(η)

ηmin

 

L−DT 2.531+DT −DP π IP o DP 2.531−DT ηmin = ηmax 1 − . (18) 16 μ α 2.531 + DT − DP ηmax

Equation 18 honor 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 Eq. 18 with Darcy’s Law:

k AT IP

QT = , (19)

μ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:

 

L1−DT 2.531+DT −Dp π o DP 2.531+DT ηmin k = ηmax 1 − . (20) 16AT α 2.531 + DT − DP ηmax

Although Eq. 20 is a description of permeability in a general porous medium, the structural

parameters embedded in AT , for example, need to be defined.

3.2.2 Establishing Total Area “AT ” for Flow and Porosity for Completely Fragmented

Porous Media

Total area of flow is the combination of areas for both pores and the rock surface exposed

to the flow. This area gives the superficial velocity in the porous medium. If only the areas

of pores are taken account, this will be the interstitial velocity, which is the average velocity

in the actual pores or channels of the medium. 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 as(Rieu and Sposito 1991):

N

AT = AP = APi . (21)

i=1

In Eq. 21, APi is the pore area fractions.

Tyler and Wheatcraft (1989) conceptualized the flow channels as tubes in a close-packed

arrangement of grains and expressed the pore-size/grain-size relation with the following equation by replacing the scaling term α, in the Arya and Paris model (1981), by the tortuosity

fractal dimension:

1/2 2

φi N1−DT ηi = Qi i . (22) 3
In this equation, Ni is the number of particles of a certain size (Qi) lining up to create the

tortuous channel with diameter ηi and tortuosity fractal dimension DT . In Eq. 22, φi is the

partial porosity in each particle size class. The summation of partial porosities gives the total

(or effective) porosity of the fragmented porous medium φT (Rieu and Sposito 1991).

When Eq. 22 is inserted into Eq. 21, the total area is obtained in terms of particle size and

partial porosities:

s 

2 πQ 1 D A = A = i φ T T p i Ni

−

 

. (23) 6

i=1

Rieu and Sposito (1991) defined the virtual pore-size fraction concept for porosity. If a fractal

porous medium, whose porosity results from a broad range of pore sizes, decreasing in mean

diameter from ηo to ηs−1(s ≥ 1), is divided into “s” virtual pore-size fractions, then the total

porosity is the summation of all partial porosities. 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 Qi is given by Eq. 22, then the number of pore-size classes

can assumed to be equal to the number of particle size classes.

By using the virtual pore-size fraction concept and introducing the pore coefficient r, the

ratio of successive increments of pore volume to the corresponding partial volumes, partial

porosity can be written as (Rieu and Sposito 1991):

φi = 1 − (1 − s r) −i (24)

Pore coefficient r is closely related to the fractal dimension DF as:

3−D r = 1 − r F (r < 1,r < 1), (25)

where “r” is a similarity ratio that scales the pore sizes in a fractal porous medium. Using

these relations, the porosity of a completely fragmented porous medium is given as:

φ = 1 − 3−DF s (r ) . (26)

This equation, in terms of particle sizes, gives total porosity as:

φT = 1 −

 3 Q −DF min , (27) Qmax

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 Eq. 1.

3.2.3 Plastic Reduction in Porosity with Applied Stress Increment in Uniaxial Compression

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). 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, which can be measured by

various methods, such as sieving or laser-light scattering (McDowell et al. 1996).

Using the principles of grain crushing and a particle size distributionmodeled as a self-similar pattern, the plastic reduction in porosity with a one-dimensional applied stress increment

is given by (McDowell et al. 1996):
m ( 	 (1−D 1−D F ) F ) ε m(DF −1)

d 2 φT = − 1 β f 2 (2 DF )mσ σ 2 0 − dσ (28) (1 − κ)σ Q ̄ max−  ̄  ̄

In this equation, DF is the fragmentation fractal dimension defined by Eq. 1, ε is the surface

energy in linear-elastic fracture mechanics (Griffith 1920), m is the Weibull modulus based

on the observations that particle-survival 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 (A) and stress increment (dσ ̄ )

terms:

m(D 1

= −  ̄ F − )

d 2 φT Aσ 2 − dσ ,  ̄ 	 (29)

where A is defined as:

− ε m(1 D (1 D − ) F ) F

A = β f 2 2 (2 DF )mσ . (30) (1 max

− − κ)σ Q ̄ 0

Equation 29 can be integrated, assuming A is constant at each stress level, between two

successive stress increments (i to i + 1) to obtain:

⎡ m(DF −1)

− m(D 1 F 2 2

−1)

1

 ̄ 1

	 

σ −

− ⎢ i φTi φTi = −A m

+ σ 1

)

 ̄

⎤

− i

(

⎥ . (D 1 m D 1) (31) + F − −1 F − −1 2 2

In order to solve Eqs. 28–31, the plastic

⎣

compressibility index (A)should

⎦

be defined by using

either substituting typical values of physical parameters in Eq. 30 or a typical value of 0.1 for

sands and clays (McDowell et al. 1996).

4 Results and Discussion

4.1 Description of Fractal Model Parameters and their Calculations

In order to predict gob porosity and permeability under increasing compressive stresses,

Eqs. 20, 23 and 27 (Eqs. 28–31 for porosity change under compression) should be evaluated.

In effect, Eqs. 28–31 evaluate the change in porosity as a function of increasing vertical

stress and fragmentation fractal dimension. The new porosity replaces the porosity term (φi)

in Eq. 23 to modify the flow area (AT ). Finally, the updated flow area is inserted into Eq. 20

to estimate the permeability at a particular stress and fragmentation level.

The general solution procedure stated above requires determination of these parameters:

•	 Representative length (Lo)

•	 Maximum and minimum particle dimensions, Qmin and Qmax , in the range(s) thatshows

fractal behavior.

•	 Weibull modulus (m)

•	 Pore, tortuosity, and fragmentation fractal dimensions, DP, DT , and DF , respectively.

•	 Virtual pore-size classes (s)

•	 Pore coefficient (r) and similarity ratio (r).

•	 Minimum and maximum pore (channel) diameters, ηmin and ηmax .

Weibull modulus (m) is a measure of the variability in strength of the material. It increases

with decreasing the variability in strength in the particle. For chalk, brick, stone, and cement,
m is between 1 and 5 (McDowell et al. 1996). Similar values would be expected for carbonates. As “m” decreases, the impurities or structural variabilities in the grains increase and

the strength of the material decreases. Thus, a low m-value around 2 should be expected for

natural materials.

Most flow occurs in the horizontal directions within the gob due to the pressure differentials between points in the gob and the nearest gateroad or longwall face. Representative

length (Lo) thus 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. According to

Eq. 2, when there is no convoluteness (DT = 1), the length of a straight channel is equal to

the representative length in any flow direction. Therefore, the representative lengths change

based on the location of the evaluation points in the gob. In this paper, Lo was taken as 10m.

The porosity and permeability values are representative of 10m long segments in the flow

direction, experiencing an average macroscopic stress of σ ̄ .

Fragmentation fractal dimension (DF ) must be <3 (Kozak et al. 1996) in order for

the probability of grain fragmentation to be <1 (Comegna et al. 1988). A fractal dimension between 0 and 3 reflects a population dominated by relatively larger particles and by

a distribution dominated by smaller particles (Turcotte 1986) as the value approaches 3.

It has also been reported that confinement under compression impacts fractal dimension

with lower fractal dimensions caused by relatively unconfined conditions and higher fractal dimensions approaching 3 resulting from extensive grinding caused by high confinement.

In the study of Pappas and Mark (1993), the particle-size data was reported as a distribution of mass percent of total particles passing through a certain sieve size (Figs. 2

and 3), rather than the actual number of particles. Therefore, initial mass was calculated

using volume of the test chamber and assuming a constant rock density of 2.3 g/cc for

shale and 2.5 g/cc for sandstone. From this data, the number of particles in each size

class (sieve size) was calculated by assuming that the particles are spherical in shape.

Example linear fits to the number–size data are given in Fig. 6 (for tests A and B, before

loading and after loading). These data show that fragmentation fractal dimensions have

increased afterthe loading testsindicating the presence of a finer particle size produced during

crushing.

The cumulative number versus particle size relation (Fig. 6), which is used to determine

DF can also be used to find the maximum and minimum particle size (Hunt and Gee 2001).

Maximum and minimum particle diameters (Qmin and Qmax) are approximately the limits of

the linearity in the number–size relation. However, it should be noted that there is a deviation

from linearity in Fig. 6 (minimally in A and more in B before loading). This may be related

due to multi-fractal distribution of rock particle size classes before loading. In this case,

a piecewise analysis of the data is a possibility. However, toward larger particle diameters

where a large portion of porosity may reside, the slope becomes >3, which physically is not

possible. Thus, in order not to exclude that porosity in the initial material, it was assumed

that linearity in the log–log behavior in Fig. 6 is maintained for all cases through the entire

range of particle sizes.

The fragmentation fractal dimensions and maximum and minimum particle sizes before

and after loading tests are listed in Table 2. This table shows that after loading, fragmentation fractal dimensions of the initial particle size distributions have increased indicating the

shift in maximum and minimum particles sizes created during crushing, as shown in Figs. 2

and 3.
The number of virtual pore-size classes “s” is equal to the number of different particle size

fractions that are within the data set where fractal scaling is observed. Since a linear trend

is accepted between all particle size classes (Fig. 6), the number of particle size classes used

will be the same as the virtual pore-size fractions (Table2).

Fig. 6 Log–log plots for

number–particle size data (Eq. 1)

for Tests A and B before and after

loading. DF’s for Test A—before

compaction and after compaction

are 2.19 and 2.89, respectively.

DF’s for Test B—before

compaction and after compaction

is 2.25 and 2.91, respectively

N (>Omega) N (>Omega) 

1.E+05 

1.E+04 

1.E+03 

1.E+02 

1.E+01 

1.E+00 

Test A-After Compaction 

Test A-Before Compaction 

A 

0.1 

1.E+05 

1.E+04 

1.E+03 

1.E+02 

1.E+01 

1.E+00 

1 10 

Particle Size (cm) 

Test B-After Compaction 

Test B-Before Compaction 

B 

0.1 1 

Particle Size (cm) 

10 

Table 2 Results of DF and maximum (Qmax) and minimum (Qmin) particle sizes determined before and

after loading tests using number–size relationship

Test Before loading After loading

s DF Qmin(cm) Qmax (cm) s DF Qmin (cm) Qmax (cm)

A 9 2.19 1.5 5.1 10 2.89 0.5 5.1

B 10 2.25 1.7 8.9 10 2.91 0.6 8.9

C 9 2.19 1.5 5.1 10 2.81 0.4 5.1

D 10 2.25 1.7 8.9 10 2.87 0.8 8.9

The pore coefficient r, which is

the ratio of successive increments of pore volume to the corresponding partial volumes and

the similarity ratio r that scales the pore sizes in a fractal porous medium can be calculated

using Eqs. 25 and 26.

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. Mandelbrot et al. (1994) proposed

that “the difference between an object’s fractal dimension and its topological (or Euclidean)

dimension could be described asthe fractal increment or FI.” They furthersuggested that this

fractal increment could be used to estimate the fractal dimension of lower-order topological

units. For example, a line drawn on a fractal surface will have the same FI as the surface. If

the fractal dimension of the surface is 2.5 and its FI is 0.5 (2.5 − 2.0 = 0.5), then the fractal

dimension of a line drawn across this surface will be 1.5 (1.0 + 0.5 = 1.5). This approach

has been used in the literature (Sammis and 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 and Adler (1988) observed

fractal increments of 0.23 in pore structures of dolomitic limestone. Tyler and 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 and Halleck (2003) determined that an average

DT increment of 0.29 for incompletely fragmented porous media by 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 and Tyler (1988) as 1.08 using Monte-Carlo simulations,

values reported by Yu and Cheng (2002), and Yu and Liu (2004) for pores in 52% porosity medium between 1.10 and 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 and Drzymala 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 and Wheatcraft 1989), its Euclidean dimension is 3. Then, the calculated increments

(DF − 3) are added to the Euclidean dimension for area to obtain DP. The calculated values

are less than 2 (Table 3).

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 3. These values will be used in model calculations later in this paper.

4.2 Up-Scaling Particle Sizes used in Experiments to Actual Gob Material

Pappas and Mark (1993) simulated the gob material to determine the tensile and compressive

rock strengths, rock density, surface roughness, and rock shape. However, they scaled down

the grain-size distribution curves for determining the stress–strain properties of granular

materials representing a gob. The ratio of scaling between actual and laboratory gob material

was approximately 10.

In order to calculate porosity and permeability using actual gob-material sizes and their

distributions, the gradation curves (Figs. 2 and 3) were up-scaled by 10 times to their original

sizes. Up-scaling, 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 Tables 2 and 3 will still be valid.

Table 4 shows the up-scaled maximum and minimum particle sizes under the assumption of

scale-invariant fragmentation.

Table 3 DT and DP values determined for the tortuous and irregular flow channels and straight and smooth

channels

Test DT (straight) DT (tortuous) DP(uniform) DP(irregular)

Before loading

A 1

B 1

C 1

D 1

After loading

A 1

B 1

C 1

D 1

1.10

1.10

1.10

1.10

1.26

1.26

1.26

1.26

2 1.19

2 1.25

2 1.19

2 1.25

2 1.89

2 1.91

2 1.81

2 1.87

Table 4 Up-scaled particle sizes

for different tests before and after

loading

Before loading

Test Qmin (cm)

A 15

B 17

C 15

D 17

After loading

Qmax (cm) Test Qmin (cm) Qmax (cm)

51 A 5.0 51

89 B 6.0 89

51 C 4.0 51

89 D 8.0 89

4.3 Porosity of Initial and Compacted Gob Material

4.3.1 Calculation of Porosity Reduction from Stress–Strain Behavior and from Fractal

Model (Eq.27) Before and After Compaction

Before calculating the porosities using particle size distributions and the theory described

in Sect. 3, porosities that should be experienced in the gob were calculated using the experimental stress–strain curve (Fig. 4) and the initial porosities given in Table 1. The results of

this calculation are shown in Fig. 7.

The porosity evaluation during uniaxial compression ofsimulated gob materialshowsthat

at earlier loadings, the material could be compressed more when compared to later stages as

also suggested by their stress–strain behavior. At around a stress level of 19MPa (2,750 psi),

the final porosities were roughly 10% for Test B, and 19–20% for the other tests. These

numbers are close to the final porosities given in Table 1. It should be noted that the final

values reported in Table 1 were obtained at higher stress levels. Thus, a difference of 3–4%

between the values given in Fig. 7 and the values in Table 1 can be considered close for this

analysis. As an inner check point, the porosity values reported for 5.4MPa (800 psi) were

compared with the corresponding data shown in Fig. 8. These values are also within 1–2% of
each other. Thus, the data of Fig. 7 will be used to compare with the porosity values calculated

using fractal models (Eq. 24–31).

Fig. 7 Porosities calculated for 0.90 

compacting gob materials under 0.80 compressive loading

0.70 

0.60 

Porosity 

0.00 

0.10 

0.20 

0.30 

0.40 

0.50 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

0.00 

0.10 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 

0.80 

0.90 

Porosity 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

Fig. 8 Porosities predicted for compacting gob materials using Eqs. 28– 31

Initial and final porosities before and after compaction can also be calculated using fractal model (Eq. 27) and the particle size distributions of fragmented material. The porosity

values obtained using Eq. 27 and the values determined in the laboratory before loading

are compared in Table 5. The pore coefficients r and similarity values r calculated using

Eq. 25 and 26, respectively, for the initial gob material are also given in Table 5. The porosity

predictions given in this table are usually within 15% compared to laboratory measurements.

This difference can be attributed to the nature of independent measurements and calculations.

Table 5 also shows the pore coefficients r and similarity values r calculated by using

predicted porosity values in Eqs. 25 and 26. The values of these two numbers should be <1.

The pore coefficient r in a fractal porous medium can be interpreted as the partial porosity

contributed to each partial volume by pores of a certain size. Similarity number r, on the

other hand, is a linear similarity ratio that relates successive pore sizes by their particle sizes

(Rieu and Sposito 1991). Thus, having close and high r-valuesindicate close scaling relations

between each pore and size class.
Table 5 Comparison of porosities determined in the laboratory with the predictions using Eq. 27 for beforeloading cases

Test Porosity-ØL -(laboratory) Porosity-ØP -(Predicted) (Eq. 27) rØp (Eq. 26) rØp (Eq. 25)

A 0.802 0.629 0.873 0.104

B 0.679 0.711 0.847 0.117

C 0.790 0.629 0.873 0.104

D 0.719 0.711 0.847 0.117

The pore coefficients, r, and similarity ratios, r, for a fractal porous medium for both sets of porosity values

are also given

Table 6 The porosities, pore coefficients, r, and similarity ratios, r, for a fractal porous medium, calculated

using the crushed material size distribution data obtained after loading

Test Porosity @

max. stress

Final porosity

from strain

Porosity -ØP­

(Predicted)

rØp (Eq. 26) rØp (Eq. 25)

(Table 1) (Fig. 7) (Eq. 27)

A 0.160 0.19 0.225 0.793 0.025

B 0.078 0.09 0.216 0.764 0.024

C 0.152 0.19 0.383 0.775 0.047

D 0.162 0.19 0.269 0.786 0.031

The porosities of simulated gob material after loading tests were also calculated by using

the particle size distributions data after the tests and the fractal model given in Eq. 27. For

these calculations, lower bounds of the particle diameters after loading, where linear number–

size relations were observed, were determined as 0.5, 0.6, 0.4, and 0.8 cm (5, 6, 4, and 8 cm

for up-scaled diameters) for tests A, B, C, and D, respectively. The porosity values obtained

using fractal fragmentation data in Eq. 27 and the porosities determined in the laboratory

after loading tests and using the strain data are compared in Table 6. The pore coefficients r

and similarity r-values calculated using Eq. 25 and 26, respectively, are also given in Table 6.

This table shows that predicted porosities after loading are lower than the calculated initial

porosities (Table 5) due to crushing of grains. The comparison of predicted porosities in

Table 6 with the values calculated using strain data (Fig. 7) and with the reported porosities

at maximum stress show that the values calculated using non-fractal methods are somewhat

lower compared to the predicted ones.

4.3.2 Predicting Compaction “Path” of Gob Material and Evolving Porosities Using

Fractal Crushing Model

Progressively increasing stresses during mining compacts the initial gob material until it can

fully support the weight of the overburden. This stress reestablishment happens at a certain

distance from the longwall face. During the stress reestablishment period, the reservoir properties of the gob material continue changing as a result of compaction and crushing under

increasing stress. In thissection, the application and the results of the fractal crushing method

introduced in Sect. 3 are discussed for its potential to model gob compaction.

Porosities were calculated using fractal compression and crushing models presented in

Eqs. 28–31.The results ofthese calculations are shown in Fig. 8.In the calculations using these

equations, the plastic compressibility index A was taken as 0.1 as suggested by McDowell et
al. (1996). The Weibull modulus m, on the other hand, was taken as 2.0, considering the fact

that shale and weak sandstone are natural materials with many defects to lead to their high

probability of fragmentation. The initial porosities plotted in Fig. 8 are the predicted initial

porosities given in Table 5.

Figure 8 shows that as stress increases, the predicted porosities using the fractal model

(Eqs. 28–31) decrease in a manner similar to Fig. 7. Also, considering the use of predicted

initial porosities given in Table 5 and the use of Weibull and plastic compression constants

together with the fragmentation fractal dimensions given in Table 2 for the initial number–

size distribution of simulated gob material, the predicted porosity values are very close to

calculations shown in Fig. 7, as well as to the final porosities given in Table 6. This suggests

that 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.

In Fig. 8, the initial porosities and fragmentation fractal dimensions are the same for each

of the two test pairs (A–C) and (B–D) and their predicted porosities coincide as opposed to

being separate as in Fig. 7.

4.4 Calculation of Permeability of Gob Material

4.4.1 Determining Flow Area of Channels (AP) in the Fractal Porous Medium

In order to calculate gob permeability, Eq. 20 was used. The permeability equation requires

flow area of pores (AP, Eq. 23) and maximum and minimum channel sizes to be known. The

size of the grains in a fragmented porous media will directly affect the pore flow area and

pore diameters forming between the grains. Also, as compression increases on the initial gob

material, the porosity will decrease as discussed in the previous section. 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 23 is used to calculate the total flow area of pores. This equation can alternatively

be represented with an integral:

Qmax

π 

QφN1−DT dQ. (32) 3

Qmin

In this integral, maximum and minimum particle sizes are the up-scaled values and the number (N) is a function of grain size (Q). In order to find this function, the numbers of up-scaled

particles for each up-scaled size were plotted. By up-scaling the material mass (assuming

density stays constant) by 1,000 times due to a 10-fold increase in particle diameter, the

number–size relations are obtained. The functional relationships for these number–size plots

are given in Table 7. These equations show that the number–size data can be successfully

represented by power functions. The R2 values of the equations representing the data in

Table 7 are about 0.9927.

In Eq. 32, Ø-values are the partial porosities for each infinitesimal size interval. Integration between Qmin and Qmax will result in the total fractal porosity (ØT ), which will not be

dependent on Q. With these treatments, the fractal flow area (Eq. 32), for Test-A for example,
becomes:

Qmax

π −1.74 φ (10326 1 Q ) −DT T

 

QdQ. (33) 3

Qmin

Solving the integral gives the fractal pore area. The fractal flow areas calculated for the initial

gob material of all tests are given in Table 8. These equations can be directly used in permeability equation (Eq. 20) for flow area term. In these equations, ØT is the porosity which

changes as uniaxial stress progressively increases according to Eqs. 28–31.

Table 7 Functional relationships

of number of particles as a

function of up-scaled particle

diameters (Q is in cm)

Test N

A Q−1.74 10326.0

B Q−1.75 8176.7

C Q−1.74 9687.1

D Q−1.75 7158.0

Table 8 Flow areas calculated

for the initial gob material for

Tests A-D. Particle diameters are

the up-scaled values and their

units are in cm

Test

A

AP (Eq. 23)

103261−DT π

3 φT

⎡

⎣

Q0.26+1.74DT max

1.74DT + 0.26 − Q0.26+1.74DT min1.74DT + 0.26

⎤

⎦

B

8176.71−DT π

3 φT

⎡

⎣

Q0.25+1.75DT max

1.75DT + 0.25 − Q0.25+1.75DT min1.75DT + 0.25

⎤

⎦

C

9687.11−DT π

3 φT

⎡

⎣

Q0.26+1.74DT max

1.74DT + 0.26 − Q0.26+1.74DT min1.74DT + 0.26

⎤

⎦

D

71581−DT π

3 φT

⎡

⎣

Q0.25+1.75DT max

1.75DT + 0.25 − Q0.25+1.75DT min1.75DT + 0.25

⎤

⎦

4.4.2 Determining Size of Tortuous Flow Channels (η)

In order to determine maximum and minimum channel sizes, Eq. 22 is used. 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 Table 7. The partial

porosity term in the equation is, on the other hand, the amount of porosity in that partial

volume. 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” defined earlier in this paper (Table 2). With these additions, the maximum and minimum

channel diameters (ηmax and ηmin) as a function of Q can be calculated. These relations are

given in Table 9.
Table 9 Minimum and maximum sizes (ηmin and ηmax) of tortuous flow channels as a function of particle

size

Test ηmax (Eq. 22) ηmin (Eq. 22)

J J

A 0.82 φT

s (102361−DT )Q0.26+1.74DT max 0.82 φTs (102361−DT )Q0.26+1.74DT min

J J

B 0.82 φT

s (8176.71−DT )Q0.25+1.75DT max 0.82 φTs (8176.71−DT )Q0.25+1.75DT min

J J

C 0.82 φT

s (9687.11−DT )Q0.26+1.74DT max 0.82 φTs (9687.11−DT )Q0.25+1.74DT min

J J

D 0.82 φT

s (71581−DT )Q0.25+1.75DT max 0.82 φTs (71581−DT )Q0.25+1.75DT min

Particle diameters to be used are the up-scaled values. The units of both particle sizes and the resulting channel

sizes are in cm

4.4.3 Predicting Permeabilities for Circular Flow Channels with Different Tortuous Paths

and for Pore Area Irregularities

Figure 9 shows the predicted permeability of gob material composed of broken rocks and

circular flow channels as a function of uniaxial compressive load for shale and weak sandstone.

 

0 

0.01 

0.02 

0.03 

0.04 

0.05 

0.06 

0.07 

0.08 

Permeability (cm2) 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

DT =1.1 DP=2 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

Fig. 9 Predicted permeabilities for gob simulated in Tests A–D assuming circular channels with DT = 1.1

Figure 9 shows that the permeability is higher initially and decreases as a result of

crushing and compaction during loading. The predicted permeabilities are changing from

0.074 cm2(74 × 105 Darcy) to 0.026 cm2 5 (26 × 10 Darcy) for a gob initially composed of

bigger rocks and from 0 2 5 2 5 .026 cm (26 × 10 Darcy) to 0.009 cm (9 × 10 Darcy) for smaller

rock particles. The data shows that the percent reduction in permeability is close in both

cases, around 65%. Figure 9 also shows that permeabilities are mainly affected by the size of

the initial material, rather than the type of rock. This is partly due to using the same Weibull

moduli (m) and plastic compressibility index (A) in the porosity reduction equations. Nevertheless, this observation is consistent with the expectations since bigger particles are prone

to leave more pore space when they are casually dumped into a finite volume.

The predicted permeabilities in Fig. 9 with DT = 1.1 and DP = 2 conditions may be

considered as an optimistic approach for permeability calculation, since circular pores with

no irregularities and relatively low tortuosity may result in higher permeability than it actually is for a gob environment. Thus, the complexity of the model environment was increased
by changing the DT - and DP-values. DP-values were changed to 1.19 for shale and 1.25 for

weak sandstone as presented in Table 3 for the initial gob material. DT , on the other hand,

was increased to 1.26 based on literature data discussed in Sect. 4.1. It is also interesting to

note that a DT of 1.26 is almost the mid-point between 1 (straight channel) and 1.5 after

which tortuosity fractal dimension does not have a physical meaning according to Wheatcraft

and Tyler (1988). The gob permeabilities predicted for more tortuous and irregular circular

channels using the updated DT - and DP-values are shown in Fig. 10a, b.

 

 

Permeability (cm2) Permeability (cm2) 

0.04 

0.035 

0.03 

0.025 

0.02 

0.015 

0.01 

0.005 

0 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

DT =1.1 DP(Shale) =1.19 

DP(SST) = 1.25 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

A 

0.016 

0.014 

0.012 

0.01 

0.008 

0.006 

0.004 

0.002 

0 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

DT = 1.26 D P(Shale) = 1.19

 D P(SST) = 1.25 

B 

Fig. 10 Predicted gob permeabilities for different DP-values of 1.19 for shale and 1.25 for sandstone (a), and

with additional increasing flow channel tortuosity (DT = 1.26) (b)

Figure 10a, b shows that as the complexity of the flow path and the irregularity of pore

areas increase, the permeability of the fragmented porous medium decreases. This is in

agreement with the physical situation and with the practical experiences. For instance, for

the case presented in Fig. 10a, when the irregularity of circular channels is increased, the

permeability values decrease by almost 50% compared to the values shown in Fig. 9 that
represent perfectly circular channels. In addition to this, when we increase the tortuosity

of irregular circular channels in the porous medium by increasing DT from 1.1 to 1.26

(Fig. 10b), the permeability values decrease to even lower values. This is also in agreement with the expectations and the physics of flow in porous media. The data presented

in Fig. 10b show that the permeability for a gob composed of smaller rock fragments initially is 0.004 cm2(4 × 105 Darcy), which decreases to ∼ 0.001 cm2(1 × 105 Darcy) when

the gob is fully compacted under 19.4MPa stress. For gobs composed of bigger rock fragments initially (Tests B and D), permeabilities decrease from 0.015 cm2(15 × 105 Darcy) to

0 005 cm2 5×105 . ( Darcy). These values represent an almost 6-fold decrease as compared to

the values in Fig. 10.

4.4.4 Predicting Permeabilities for Gobs with Flow Channels in the Shape of Triangle

and Ellipse

In a gob environment, flow channels may vary in shape from circular to triangular to elliptical. Equation 20 predicts permeability based on different-shaped channels through the use

of a shape factor α, as discussed in Sect. 3.2.1. In this study, we have calculated the gob permeability for equilateral triangle- and ellipse-shaped flow channels as a function of loading.

In these calculations, DT and DP were 1.1 and 2, respectively. Since both of these shapes

are regular

√

2-D objects filling their area domain completely, their DP is also 2. Equation 13

and 20 3 were used in Eq. 20 as shape factors for an ellipse and a triangle, respectively.

Figure 11a, b presents the predicted values for these two situations.

Figure 11a, b shows that the permeability decreases about an order of magnitude in the

case of triangular pore channels(Fig. 11a) compared to the circular ones(Fig. 9) and about 20

times, when the channels are in the shape of an ellipse with a 5:1 aspect ratio (Fig. 11b). As

the aspect ratio is increased to 10:1 (shape factor 126.92) and tortuosity and pore-area fractal

dimensions are changed to values in Fig. 10b for more tortuous flow paths and irregularly

shaped pores, the permeabilities decrease even more. Permeabilities calculated for this situation are given Fig. 12. This may represent a situation where irregular, slit-like channels are

meandering in a fractal fragmented porous medium. This seems to be a closer approximation

to a gob where flatsandstone and shale rock blocks are piled (Fig. 1) and are compressed under

overburden loading. Calculations show that in a fractal porous medium as the complexity

and irregularity of pores increase, the permeabilities decrease accordingly.

4.4.5 Comparison of Fractal Model Prediction with Reported Values and the Predictions

of Carman–Kozeny (C–K) and Happel (H) Permeability Equations

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 −1×107Darcy). 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 andBalusu (2002) used

maximumvalues around 1×10−9m2 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 and Karacan (2007) reported calculated values of 1 × 10−9 m2(1 × 103 Darcy).

The differences in permeabilities may be related to different coal-seam geologies, different

panel layouts and to differences in caving characteristics. Nevertheless, it should be noted
0 

0.001 

0.002 

0.003 

0.004 

0.005 

0.006 

0.007 

0.008 

Permeability (cm2) 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

For flow channels in 

the shape of 

equilateral 

triangles a a 

a 

A 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

0 

0.0005 

0.001 

0.0015 

0.002 

0.0025 

0.003 

0.0035 

0.004 

Permeability (cm2) 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

For flow channels in the 

shape of an ellipse 

a=10 b=2 

B 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

Fig. 11 Predicted permeabilities for a gob simulated in Tests A–D assuming flow channels are in the shape

of an equilateral triangle (a) and an ellipse (b) with 5:1 aspect ratio between its major axes. DT = 1.1 was

used as tortuosity

 

0 

0.00005 

0.0001 

0.00015 

0.0002 

0.00025 

0.0003 

0.00035 

0.0004 

Permeability (cm2) 

Shale-Test A 

Shale-Test B 

Weak sst-Test C 

Weak sst-Test D 

For flow channels in the 

shape of an ellipse 

a=20 b=2 

DT = 1.26 D P(Shale) = 1.19

 D P(SST) = 1.25 

0.00 5.00 10.00 15.00 20.00 

Compressional Loading (MPa) 

Fig. 12 Predicted permeabilities for gob simulated in Tests A–D assuming flow channels are in the shape of

an ellipse with 10:1 aspect ratio between its major axes. DT = 1.26 and DP-values of 1.19 and 1.25 were

used for fractal dimensions
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 analyzes predict permeabilities close to the values reported in other studies.

In order to compare the permeability values obtained using the calculation methodology

developed in this paper with popular and simpler equations used for flow in aggregates, Carman–Kozeny (C–K) and Happel (H) equations (Li and Logan 2001) presented in Eqs. 34 and

35, respectively, were used. These equations were used to calculate permeability variations

only for Test-A. When using these equations, porosity reduction was accounted for as the

stress increased and characteristic aggregate size was taken as the average of maximum and

minimum particle sizes (33 cm) reported for this test. Figure 13 shows the calculations with

C–K and H equations.

0 

2 

4 

6 

8 

10 

12 

0.00 5.00 10.00 15.00 20.00 

Permeability (cm2)

Carman-Kozeny 

Happel 

Test - A 

Permeability without 

porosity reduction 

Compressional Loading (MPa) 

Fig. 13 Predicted permeabilities for gob using Carman–Kozeny and Happel equations (Eqs. 34 and 35)

As this figure indicates, they predict values close to each other. However, these values are orders of magnitude higher than the reported gob permeability values

in this and other studies. They are also far from being realistic for a compacted gob. Thus,

fractal models developed in this study predict gob permeabilities more realistically than C–K

and H equations and thus they are more adequate to use for permeability predictions of a gob

environment.

3

=  ̄ 2

k Q φ

C−K T	 

180	 1− 2 (34) ( φT )

k =	 

Q ̄ 2 3−4

 

5 .5φ+4.5φ

 

3 6 φ H 18 3 3+2 5 (35) φ ( φ )

−


 

√
 where φ = 3 (1 − φT )

( )

5 Summary and Conclusions

This study developed a predictive approach using the principles of scaling and fractal porous

medium combined with fluid flow. This approach allowed the calculation of porosity and

permeability from the size distribution of broken rock material in the gob using flow and

fractal crushing equations for granular materials. The specific conclusions made from this

study can be listed as follows.

(1)	 Particle size distributions of gob material can be determined from digital pictures or

from simulated tests. These data can be used to find the fragmentation fractal dimen­
sion of gob material based on a number–size relationship. This study showed that

the number–size relation for different rock materials and gradations can be described

with a linear relation. This indicates that the gob material before and after compaction

shows fractal scaling properties. Tortuosity and pore-area fractal dimensions can then

be determined using the “slit-island” theorem.

(2)	 Permeability in a fractal porous media 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, then which can be integrated with the Darcy equation to describe permeability of a fractal porous media in terms of textural and fractal

properties.

(3)	 Fractal crushing of granular materials can be used in the definition of flow area to simulate compaction of a gob. 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. Thus, it is possible using the number–size distribution of a

granular media and itsfragmentation fractal dimension to calculate changing porosities

as stress is progressively increased.

(4)	 It was observed that the developed fractal model predicts gob permeability at any loading level better than simple models,such asthe Carman–Kozeny and Happel equations,

since it includes complexities of pores and grains in a completely fragmented porous

medium. Predicted permeability values calculated using the fractal permeability model

are in better agreement with geomechanical model calculations and reported values in

the literature.

(5)	 It was observed that the size of the rock blocks in the gob make a more dramatic

difference in permeability compared to the type of rock material. However, it should

also be emphasized that particle size segregation is important as the size distribution

controls the fractal properties. Thus, if there is segregation in particle size distribution,

at different heights in the gob, this should affect the transport properties of porous

medium.

(6)	 Using the developed approach and the formulations, it is possible change the pore shapes

and irregularities determined from gob pictures or to change the tortuosity of flow channels based on expectations to obtain better results. Geometrically, the results showed

that circular and smooth channels give the highest permeability. As the irregularity and

aspect ratio of the channels increase, permeabilities decrease as expected.

(7)	 The ultimate benefit of using the approach presented in this paper will be the ability to

make better gob reservoir property predictions and to ensure better methane control in

the mine environment.

References

Ahmed, A.M.H., Drzymala, J.: Two-dimensional fractal linearization of distributed curves. Physicochem.

Probl. Miner. Process. 39, 129–139 (2005)

Arya, L.M., Paris, J.F.: A physicoemprical model to predict the soil moisture characteristics from particle-size

distribution and bulk density data. Soil Sci. Soc. Am. J. 45, 1023–1030 (1981)

Arya, L.M., Leij, J.F., Van Genucten, M.Th., Shouse, P.: Sealing parameter to predict the soil water characteristics from practicle size distribution data. Soil. Sci. Soc. Am. J. 63, 510–519 (1999)

Bitelli, M., Gaylon, S.C., Flury, M.: Characterization of particle-size distribution in soils with a fragmentation

model. Soil Sci. Soc. Am. J. 63, 782–787 (1999)

Brunner, D.J.: Ventilation models for longwall leakage simulation. 2nd US Mine Ventilation Symposium,

23–25 September, Reno, NV (1985)
Comegna, V., Damiani, P., Somella, A.: Use of fractal model for determining soil water retention curves. Geoderma 85, 307–323 (1988)

Esterhuizen, G.S., Karacan, C.Ö.: 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 (2005)

Esterhuizen, G.S., Karacan, C.Ö.: A methodology for determining gob permeability distributions and its application to reservoir modeling of coal mine longwalls. 2007 SME Annual Meeting, Denver, CO (2007)

Griffith, A.A.:The phenomena ofrupture and flow in solids. Philos.Trans.R. Soc.Lond. A221, 163–198 (1920)

Hunt, A.G., Gee, G.W.: Application of critical path analysisto fractal porous media: comparison with examples

from the Hanford site. Adv. Water Resour. 25, 129–146 (2001)

Jacquin, C.G., Adler, P.M.: Fractal porous media II: geometry of porous geologic structures. Transp. Porous

Media 2, 571–596 (1988)

Karacan, C.Ö., Halleck, P.M.: A fractal model for predicting permeability around perforation tunnels using

size distribution of fragmented grains. J. Pet. Sci. Eng. 40, 159–176 (2003)

Karacan, C.Ö., Esterhuizen, G.S., Schatzel, S.J., Diamond, W.P.: Reservoir simulation-based modeling

for characterizing longwall methane emissions and gob gas venthole production. Int. J. Coal Geol.

71(2–3), 225–245 (2007)

Kozak, E., Pachepsky, Y.A., Sokolowski, S., Sokolowska, Z., Stepniewski, W.: A modified number-based

method for estimating fragmentation fractal dimension ofsoils. Soil Sci. Soc. Am.J 60, 1291–1297 (1996)

Li, X.-Y., Logan, B.E.: Permeability of fractal aggregates. Water Resour. 35, 3373–3380 (2001)

Mandelbrot, B.B., Passaja, D.E., Paulley, A.J.: Fractal character of fracture surfaces of metals. Nature 21, 308 (1994)

McDowell, G.R., Bolton, M.D., Robertson, D.: The fractal crushing of granular materials. J. Mech. Phys.

Solids 44, 2079–2101 (1996)

Mortensen, N.A., Okkels, F., Bruss, H.: Reexamination of Hagen–Poiseuille flow: shape dependence of the

hydraulic resistance in microchannels. Phys. Rev. E 71(057301), 1–4 (2005)

Pappas, D.M., Mark, C.: Behavior of simulated longwall gob material. Report of Investigations RI No. 9458,

US Department of Interior, US Bureau of Mines (1993)

Perfect, E.: Fractal models for the fragmentation of rocks and soils: a review. Eng. Geol. 48, 185–198 (1997)

Perfect, E., Kay, B.D.: Applications of fractals in soil and tillage research: a review. Soil Tillage Res. 36,

1–20 (1995)

Ren, T.X., Edwards, J.S., Josefowicz, R.R.: CFD modeling of methane flow around longwall coal faces. 6th

International Mine Ventilation Congress, 17–22 May, Pittsburgh, PA (1997)

Rieu, M., Sposito, G.: Fractal fragmentation, soil porosity and soil water properties: theory. Soil Sci. Soc. Am.

J. 55, 1231–1238 (1991)

Sammis, C.G., Steacy, S.J. : Fractal fragmentation in crustal shear zones. In: Barton, C., La Pointe,

P.R. (eds.) Fractals in the Earth Sciences, pp. 265 Pleanum Publishin Co., New York, NY (1995)

Schlueter, E.M., Zimmerman, R.W., Witherspoon, P.A., Cook, N.G.W.: The fractal dimension of pores in

sedimentary rocks and its influence on permeability. Eng. Geol. 48, 199–215 (1997)

Turcotte, D.L.: Fractals and fragmentation. J. Geophys. Res. 91, 1921–1926 (1986)

Tyler, S.W., Wheatcraft, S.W.: Application of fractal mathematics to soil water retention estimation. Soil Sci.

Soc. Am. J. 53, 987–996 (1989)

Tyler, S.W., Wheatcraft, S.W.: Fractal scaling of soil particle size distributions: analysis and limitations. Soil

Sci. Soc. Am. J. 56, 363–369 (1992)

Weiss, J.: Fracture and fragmentation of ice: a fractal analysis of scale invariance. Eng. Fract. Mech. 68, 1975–

2012 (2001)

Wendt, M., Balusu, R.: CFD modeling of longwall goaf gas flow dynamics. Coal Saf. 20, 17–34 (2002)

Wheatcraft, S.W., Tyler, S.W.: An explanation to scale-dependent dispersivity in heterogeneous aquifers using

concepts of fractal geometry. Water Resour. Res. 24, 566–578 (1988)

Whittles, D.N., Lowndes, I.S., Kingman, S.W., Yates, C., Jobling, S.: Influence of geotechnical factors on gas

flow experienced in a UK longwall coal mine panel. Int. J. Rock Mech. Min. Sci. 43, 369–387 (2006)

Yavuz, H.: An estimation of cover pressure re-establishment distance and pressure distribution in the goaf

longwall coal mines. Int. J. Rock Mech. Min. Sci. 41, 193–205 (2004)

Yu, B.M., Cheng, P.: A fractal model for bi-dispersed porous media. Int. J. Heat Mass Transf. 45(11), 2983–

2993 (2002)

Yu, B.M., Liu, W.: Fractal analysis of permeabilities for porous media. AIChE J. 50(1), 46–57 (2004)