A formalism to compute permeability changes in anisotropic fractured rocks due to arbitrary deformations

International Journal of Rock Mechanics & Mining Sciences 125 (2020) 104159

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences journal homepage: http://www.elsevier.com/locate/ijrmms

A formalism to compute permeability changes in anisotropic fractured rocks due to arbitrary deformations Andy Wilkins *, Qingdong Qu Coal Mining Research Program, Queensland Centre for Advanced Technologies, CSIRO, PO Box 883, Kenmore, Q4069, Australia

A R T I C L E I N F O

A B S T R A C T

Keywords: Permeability Fractured rock Deformation Permeability enhancement

The Liu-Elsworth (1997) formalism that predicts permeability in deformed fractured rock has been highly useful and successful in numerous engineering situations involving mining, coal-bed methane and shale gas extraction. The formalism assumes the rock’s fractures are aligned with three orthogonal axes and that the rock is subjected to strains along these three directions only. Here we relax these assumptions: our proposed formula yields the permeability tensor for an anisotropic arbitrarily-fractured rock subjected to arbitrary deformations. Unlike the original Liu-Elsworth approach, our formula may therefore be used confidently in situations involving any stress and strain states, including those that are non-diagonal, and with rocks with any distribution of fracture orientations, including those that are nonorthogonal. This is advantageous in numerical simulation and in mining situations that frequently involve such situations. Our formula adds no discernible computational load to numerical simulations. Our proposal agrees with lab and in-situ direct and indirect tests of permeability.

1. Introduction Permeability changes due to redistribution of stresses or strains1 are important in many engineering fields. For example, in mining,2–4 during tunnel drivage or resource exca­ vation, stress relaxation and strata deformation and fracture cause changes in hydraulic conductivity of several orders of magnitude.5–20 This may result in large groundwater flows into the working area.21–25 Operationally, this often reduces the rate of advance. It may also affect the groundwater system: drawdown depletes aquifers for domestic, agricultural or commercial use; clean aquifers can be mixed with dirty ones; surface water features such as swamps and ponds can be elimi­ nated; and baseflow to rivers may be reduced or otherwise impacted.26–31 Similarly, in coal mining permeability changes caused by stress/ strain redistribution can lead to large methane gas flows to the working face,8,11,32–34 and can effect outbursts,35 resulting in operational haz­ ards. Permeability changes can cause direct methane emissions to the biosphere. Conversely, the permeability changes can substantially improve gas drainage using boreholes drilled into the highly-permeable regions.34,36–47 In coalbed methane and shale gas extraction, the effect of mechanical

deformations on permeability can substantially alter the rate of gas flow to boreholes.48–55,56 Understanding these permeability changes is therefore clearly important. One of the most successful formulations was proposed in 1997 by Liu and Elsworth.1 It was used in many of the aforementioned citations. Liu and Elsworth conceptualize rock as a fractured medium. Fluid moves along these fractures, so the fracture apertures and the density of fractures control the ease of fluid flow (the permeability of the rock). When the rock is deformed, these fractures open or close, and perhaps new fractures are created, and this changes the permeability. There are a number of alternate models that describe permeability in rocks, and a huge number of articles concerning these models. Review articles52,57–62 can guide an interested reader in this field. The theory behind these models is frequently based on the “cubic law” (Eqn (5), below) which also lies at the heart of Liu-Elsworth’s model. Liu-Els­ worth’s model expresses the permeability as a function of strain, while other models use porosity (e.g., the Kozeny-Carman model63–65) or linear-elastic stress66 or elastic strain67 perhaps complimented by a consideration of damage.68 One facet of Liu-Elsworth’s model is that it is successful outside the linear-elastic regime, without having to assume porosity is unrealistically high. Alternately, permeability may be measured experimentally as a function of stress or strain, and empirical

* Corresponding author. E-mail address: [email protected] (A. Wilkins). https://doi.org/10.1016/j.ijrmms.2019.104159 Received 16 January 2019; Received in revised form 31 October 2019; Accepted 13 November 2019 Available online 27 November 2019 1365-1609/Crown Copyright © 2019 Published by Elsevier Ltd. All rights reserved.

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

relationships formed, perhaps guided by theory. For instance, perme­ ability changes caused by induced microcracking in highly compressive regimes has been measured,69,70 and diagonal stress states were observed to induce non-zero off-diagonal permeability-tensor compo­ nents,69 as also occurs in our formalism. It is difficult to upscale these laboratory observations to mining scenarios that involve large exten­ sional deformations. The Liu-Elsworth formulation is conceptually pleasing, which ac­ counts for part of its popularity. On the practical side, it contains only one free parameter (the modulus reduction ratio) making it desirable for numerical modellers and for straightforward application in complicated mining scenarios. Importantly, the Liu-Elsworth formulation appears to provide good agreement with experimental observations. It is based on the cubic law, which follows from the Navier-Skokes equations assuming laminar flow between smooth parallel plates, and is fairly well-established.59,71 The Liu-Elsworth model adds a “modulus-reduction” parameter to the cubic law. The model has been compared with lab tests using small stresses and strains.53,72,73 It has also been compared with in-situ measurements of permeability in mining scenarios.5,13,74 It has also been compared with indirect observations of rock-mass permeability in large-strain settings.3,6,7,9–11,21,26,27,32 In these works, rock-mass permeability is inferred by observations of water or gas flows through reservoirs, into boreholes or mining excavations, or by using piezometer responses, and then comparing with numerical modelling. This paper seeks to extend the Liu-Elsworth formulation so that it is applicable to more general situations. The following three extensions are considered. Firstly, the original formalism conceptualised the fractures as being infinite planes aligned with the x, y and z axes. Liu and Elsworth were particularly interested in coal, where the cleating system is largely mutually orthogonal, but their formalism has been used extensively for other rocks where the natural fracture sets are not always so. Relaxing this assumption allows consideration of rocks with arbitrarily-aligned fracture planes, or with a distribution of fracture-plane alignments. To do this, we build on Snow’s 1969 work.75 Secondly, the original formalism assumed the applied strain was diagonal. This was extended by Liu et al.73 to small shear strains by introducing an extra modulus-reduction parameter and making the as­ sumptions that the stress-to-strength ratio was low and the fractures are aligned with the x, y and z axes. Here we relax these assumptions and employ just one modulus-reduction parameter, which allows the permeability tensor of arbitrarily deformed rock to be computed. The result transforms correctly under rotations. Thirdly, the original formalism used infinitesimal strain. However in mining applications, the permeability is often enhanced by several or­ ders of magnitude. This suggests a finite strain approach may be pref­ erable, so our approach does not assume infinitesimal strain. We explore a situation involving shear strains in which considering large strains is necessary to produce the expected result. Related to this, the original formulation has the potential to produce negative permeabilities (if βε < 1, which occurs frequently in finite-element simulations during nonlinear convergence to a final solution): we propose a simple fix. Our proposal agrees with experimental data collected in the smallstrain regime,53,72,73 because these experiments recorded only the di­ agonal permeability components when samples were subjected to di­ agonal stresses, and in this case our proposal produces the same result as the original Liu-Elsworth formulation. For the same reason, it agrees with in-situ measurements of permeability in mining scenarios.5,13,74 We have performed 3D finite-element simulations (similar to Section 6) and have found no discrepancies with the indirect observations of perme­ ability reported in the literature.3,6,7,9,10,21,26,27 We could find no experimental data that can directly support or refute the strain-dependence of the off-diagonal permeability components in our formalism. However, because any applied stress or strain may be diag­ onalized by rotating the frame of reference, and our proposal agrees

with experiment for diagonal stress or strain, it is likely the off-diagonal permeability components are correct because our tensor transforms correctly under rotations. The organization of this paper is as follows. In Section 2 we sum­ marize Darcy flow, permeability, and flow through fractures. In Section 3 we summarize the original Liu-Elsworth formalism. In Section 4 we explain how to include distributions of arbitrarily-aligned fracture planes and write Snow’s75 formula for permeability using distributions. In Section 5 we explain our proposal for computing permeability changes resulting from general, finite deformations. In Section 6 we provide two case studies comparing our new formulation with the original. A number of illustrative examples are presented in Appendices. 2. Permeability A fluid’s Darcy velocity quantifies the flow of fluid through the pores or fractures of a medium (such as rock) through which it is flowing. Denote the components of the Darcy velocity by vi , with i ¼ ðx; yÞ in two dimensions and i ¼ ðx; y; zÞ in three dimensions. The Darcy velocity has SI units m.s 1. Darcy’s law76 states that the Darcy velocity is vi ¼

kij

μ

(1)

rj p :

(For instance vx ¼ ðkxx =μÞrx p ðkxy =μÞry p ðkxz =μÞrz p.) In this equation rj p is the gradient of the fluid porepressure p, μ is the fluid dynamic viscosity, and the components of the permeability tensor are kij . The SI units of rj p are Pa.m 1, μ has units Pa.s, and kij has units m2. In this equation there is an implied sum over the repeated index j (j is summed over x and y in 2D, or x, y and z in 3D). Gravity is not included in this formula9: it makes no difference to the arguments below. Henceforth, we use the Einstein convention that there is an implied sum over repeated indices, unless explicitly stated. Now consider fluid flowing through a planar fracture with unit normal vector b n i . Irrespective of the pressure gradient, the Darcy ve­ locity must be oriented along the fracture plane, that is b n i vi ¼ 0 (recall there is an implied sum over the repeated index, i). Put another way, the Darcy velocity should depend only on the tangential component of the porepressure gradient, because that is the component that is driving the fluid along the fracture plane. This tangential component is � n j rj p ¼ δij b n j rj p : ri p b ni b ni b (2) (Here δij ¼ 1 if i ¼ j, and is zero otherwise.) This means that the permeability of a material containing just such a fracture must be kij ∝ δij

(3)

nj : b ni b

This formula was explored in detail by Snow.75 For the special case of sets of orthogonal fractures, Eqn (3) has been explored in detail by Liu et al.5 These authors present explicit examples along with beautiful di­ agrams illustrating the angular dependence of permeability. Throughout this paper we use the projection tensor Pij Pij � δij

(4)

nj : b ni b

n (Pij b n j ¼ 0) and is idem­ It is symmetric (P ¼ P), orthogonal to b potent (Pij Pjk ¼ Pik ). An example in 2D is given in Eqn (A.2) and in 3D in Eqn (D.1). In Eqn (4) we have used the � symbol: we mean “defined to be equal to”. In summary, for a material containing a single fracture with unit normal vector, kij ∝Pij . The objective of this paper is to generalize this to a material containing multiple fractures and state the coefficient(s) of proportionality. T

2

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

3. The Liu-Elsworth formalism Liu and Elsworth1 conceptualize fractured rock as consisting of solid rock containing sets of fractures whose normals are aligned with the x, y and z axes. The fractures have apertures bx , by and bz , respectively. The spacing between successive fracture planes normal to the x axis is sx ; while sy and sz are defined similarly. All these parameters are assumed to be measured in the undeformed rock, and bi and si are just 6 numbers, not the components of covariant vectors. The permeability tensor, k, is anisotropic in this situation, but di­ agonal. If the fractures are infinite planes, thin, and not filled with sand or other obstructions, the non-zero components of k are71 kxx ¼

b3y b3 þ z ; 12sy 12sz

kyy ¼

b3 b3x þ z ; 12sx 12sz

kzz ¼

b3y b3x þ : 12sx 12sy

Fig. 1. Conventions used for the azimuthal and polar angles, θ and φ.

(5)

(The so-called “cubic law”.) Liu and Elsworth do not demand that these formulae hold: instead they simply assume that the permeability, spacings and apertures have been measured in the undeformed rock. Nevertheless, Eqn (5) motivate the Liu-Elsworth formulae for the deformed case, below. Suppose that the rock is deformed via � i ~x ¼ δij þ εij xj : (6) ~i are the final positions after Here xi are the original positions and x application of the deformation. In this paper, a tilde above a quantity designates “the quantity after deformation has been applied”. Liu and Elsworth assume the strain tensor ε, is diagonal and infini­ tesimal: 1 0 εxx 0 0 ε ¼ @ 0 εyy 0 A ​ and ​ ε2 � 0 ; (7) 0 0 εzz

Fig. 2. A fracture (shaded grey) of aperture b is initially oriented along the y axis. It is rotated and deformed by the application of a general, linear strain of Eqn (18), which results in the blue-shaded fracture. The tangent vector, t, transforms contravariantly as specified in Eqn (19). This means the normal vector, n, must transform covariantly, as in Eqn (22). A point on the original fracture right-hand surface, bnT , is transformed to bMnT on the deformed fracture right-hand surface, as in Eqn (32). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

damage.5,73,74,77 Having a formalism with just one tunable parameter is quite advantageous in practice and is one of the reasons the Liu-Elsworth formalism has been so popular and successful. By way of motivation and explanation of the Liu Elsworth formula, consider the case Rm ¼ 0 where the deformation is entirely due to changes in fracture aperture (all the strain is accommodated by the fractures and the solid rock matrix does not deform). Then

Liu and Elsworth postulate that the permeability tensor becomes � � ~kxx ¼ 1kxx 1 þ βy εyy 3 þ ð1 þ βz εzz Þ3 ; 2 � ~kyy ¼ 1kyy ð1 þ βx εxx Þ3 þ ð1 þ βz εzz Þ3 ; 2 �� ~kzz ¼ 1kzz ð1 þ βx εxx Þ3 þ 1 þ βy εyy 3 ; 2

(8)

~ bi þ ~si ¼ ðbi þ si Þð1 þ εii Þ ¼ ~ bi þ si :

where βi ¼ 1 þ ð1

Rm Þ

si : bi

(10)

Here i is either x, y or z (there is no sum over i). The first equality is the definition of strain, while the second holds because Rm ¼ 0 so ~si ¼ si . This yields � � � � si ~ bi ¼ bi 1 þ 1 þ εii ¼ bi ð1 þ βi εii Þ ; (11) bi

(9)

Here i is either x, y or z (there is no sum over i), and the βi are not components of a covariant vector: they are just three numbers in the Liu Elsworth formalism. In Eqn (9) Rm is the so-called “modulus reduction ratio”. It is the only unknown input parameter to the formalism. The modulus reduction ratio satisfies 0 � Rm � 1 and is physically thought of as the ratio of rock-mass modulus to rock-matrix modulus. It apportions applied strain between the fracture and matrix material, and controls the magnitude of the permeability changes resulting from application of a strain: if Rm ¼ 0 the permeability changes are maximized, while if Rm ¼ 1 the perme­ ability changes are minimized. It can be quite difficult to choose an appropriate Rm in a finite-element simulation of real-life phenomena, although it may be correlated with other measures of rock

where the last equality holds by the definition of β in Eqn (9) with Rm ¼ 0 (again, there is no sum over i). Using the “deformed” versions of Eqn (5) yields equations such as ~kxx ¼

�3 b3y b3 1 þ βy εyy þ z ð1 þ βz εzz Þ3 ; 12sy 12sz

(12)

which are similar to the Liu-Elsworth forms given in Eqn (8). Before moving on, note that the strain applied to the rock is ε, but Liu and Elsworth postulate that the permeability “sees” an effective strain of βε. Because β≫1 (usually) the permeability changes can be very large: 3

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

~ k≫k.

(B.2), (C.2) and so on).

4. Arbitrarily-aligned fracture planes

5. Application of a general, finite deformation

4.1. Notation used

Liu and Elsworth postulated that applying an infinitesimal and di­ agonal strain to the fractured rock enhances the fractures and so en­ hances the permeability, and their reasoning led to Eqn (8), above. In this section, we extend Eqn (8) to include general, finite deformations. In contrast to the proposal by Liu et al.73 where infinitesimal shear strains were also considered, we do not introduce an extra modulus-reduction parameter and do not assume small stress-to-strength ratios and do not assume 3 orthogonal sets of fractures. Our result transforms correctly under rotations of the global coordinate system. It should be stated that our proposal may need modification in the case of cyclic loading or shear failure, both of which may alter the fracture-plane morphology, and thereby potentially also impact fracture aperture. We conceptualize the deformation as simply rotating and translating the fracture-plane surfaces, which opens, closes and rotates the fractures without changing their physical characteristics. A natural way of incorporating morphological changes would be to include damage dependence in the function f.

The alignment of fracture planes in two dimensions can be param­ eterized by the angle, θ, to the x axis, while in three dimensions the conventions used in this paper for the azimuthal angle θ and the polar angle φ are shown in Fig. 1. The unit normal vector to a fracture plane, b n, determines its orientation: � ðcos θ; sin θÞ in 2D b ni ¼ (13) ðsin φ cos θ; sin φ sin θ; cos φÞ in 3D For instance, Fig. 2 shows a fracture (shaded grey) that has normal lying along the x axis. In the following presentation it is notationally convenient to use Ω to represent either θ (in 2D) or both ðθ; φÞ (in 3D). Our convention for integration over the angular coordinate(s) is 8Z π > > dθ in 2D > Z < 0 dΩ ¼ Z θ¼π Z ϕ¼π (14) > > > sin ϕ dϕdθ in 3D : θ¼0

ϕ¼0

5.1. Finite deformations

R 2π The standard integration over θ is 0 dθ, but in the case of fracture planes, a fracture plane oriented at angle θ is indistinguishable from a fracture plane oriented at angle θ þ π, so we integrate θ from 0 to π only.

The theory of infinitesimal and finite strains is found in many text books, for instance see.79 Denote the positions in the undeformed rock by xi (which are ðx; yÞ in two dimensions and ðx; y; zÞ in three di­ mensions). Deform the fractured rock using the vector ui ¼ ui ðxÞ, so that the new position vector is

4.2. Fracture planes The distribution of fracture plane orientations, apertures and spac­ ings may be parameterized by two functions of Ω. �

Infinitesimal strain is when ∂u =∂x ≪1 for all i and j. As an example, consider the general linear 2D deformation of the form i

b3 12s ðΩÞ,

which is the quotient of aperture-cubed and spacing (divided by 12 for convenience). Introduce the distribution K ðΩÞ defined by

K dΩ �

b3 : 12s

(17)

~xi ¼ xi þ ui : j

ux ¼ εxx x þ εxy y ​ and ​ uy ¼ εyx x þ εyy y :

(15)

(18)

Here the εij are components of a strain tensor (but ε need not be infin­ itesimal in the following) but a rigid-body rotation has also been included since ε is not symmetric. The effect of this deformation on a fracture aligned along the y axis is shown in Fig. 2. This example is referred to throughout the text below. The tangent vector(s) to a fracture plane transform as contravariant vector(s):

� K is a distribution, meaning that its units are m2.rad 1 in 2D and m2. steradian 1 in 3D. Examples are given in Eqn (A.3), (B.1), (C.1), (D.2) and (E.1). � bs ðΩÞ, which is the quotient of aperture and spacing. Recall that the normal vector to the fracture plane, Eqn (13), de­ termines its orientation.

~ti ¼

∂~xi j t � Mij tj ; ∂xj

(19)

A pictorial example is shown in Fig. 2. Here we have introduced

4.3. Permeability Our proposal is that the permeability of the undeformed rock is Z kij ¼ dΩ K ðΩÞf ðΩÞPij ðΩÞ: (16)

Mij �

∂~xi : ∂xj

(20)

For the general 2D linear deformation of Eqn (18), M has components � � εxy 1 þ εxx Mij ¼ : (21) εyx 1 þ εyy

This involves the distribution K of Eqn (15) (but not the function b= s: that will appear in the formula for the permeability in the deformed rock, below) and the projection operator of Eqn (4). The function fðΩÞ accounts for sand,71 etc, that may partially fill the fractures, which effectively decreases b3 =s in Eqn (5). In the following presentation we mostly use f ¼ 1 to reduce the complexity of the equations. Eqn (16) was written in 1969 by Snow,75 using sums rather than integrals and distributions, who proceeded to make a detailed Monte-Carlo study of the permeability of special cases of 1, 2 and 3 fracture sets. In 2D it was further studied for a variety of different fracture distributions by Panda and Kulatilake.78 Examples of computing kij in 2D and 3D are presented in Appendices (see Eqn (A.4),

The normal vector must transform covariantly to maintain a zero inner product with the tangent vector(s): � � ∂xj b ni ¼ i b ~ nj N ¼ b n j M ji 1 N : (22) ∂~x ~ j ¼ 1, viz Here N is a normalization to ensure that j b n � � n M 1 j2 : N 2 ðΩÞ ¼ �b

(23)

An example in 2D is Eqn (A.7), and a pictorial representation is 4

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

shown in Fig. 2.

(11) in the context of finite strain. Denote the normal vector to the ~ (Eqn (22)). fracture before deformation by b n , and after deformation by b n

Suppose that one surface of a fracture passes through the origin. The

5.2. Permeability

T other surface passes through the point b b n , where b is the fracture aperture (we consider the point to be a column vector for convenience below). For ease of notation, consider translation-free deformations, in which the origin is mapped to the origin. The generalization to include translations will be seen to be trivial. Such a situation is shown in Fig. 2. Under such deformations, the point on the other surface deforms as:

~ ¼ x þ uðxÞ, the permeability tensor After applying the deformation x components are Z ~kij ¼ ~ K~ ðΩÞ ~ ~f ðΩÞ ~ P ~ ~ ij ðΩÞ: dΩ (24) This follows directly from Eqn (16). The tildes indicate the quantity after the deformation has been applied. The question is, what exactly are ~ etc? The next sections derive and propose forms that allow Eqn K~ , P, (24) to be written in terms of the original K , etc, and the deformation u.

where M is the matrix defined in Eqn (20). The projection of this vector ~ is the new fracture aperture, ~b: onto b n

5.2.1. Rotations Before considering the effect of fracture opening and closing on the permeability, consider the simpler case of u describing a rotation. In this ~ under the action of u. For instance, in 2D, θ→ ~θ ¼ θþ α for case Ω→Ω some constant α. ~ ¼ dΩ. For instance, in the 2D example, d~θ ¼ dθ. Clearly dΩ ~ ¼ K ðΩÞ. That is, the The distribution K →K~ , and clearly K~ ðΩÞ ~ is just the value of K at the angle Ω new value of K at a certain angle Ω which transformed into this angle. For instance, with ~θ ¼ θ þ α, if K ðθÞ ¼ sin θ, then ~ ¼ K ðΩÞ ¼ sin θ ¼ sinð~ K~ ðΩÞ θ

� ~ n ⋅bMb ~ nT ¼ b N ; b¼b

All terms — β, b and N — are functions of Ω in general. 5.2.4. Proposed Liu-Elsworth extension

Combining Eqn (30) with (33) our proposal is that

This may be trivially extended to the f 6¼ 1 case yielding ~ ~fðΩÞ ~ ¼ K ðΩÞfðΩÞ. K~ ðΩÞ ~: The projection operator depends on b n ~ ¼ δij ~ ij ðΩÞ P

~ K~ ðΩÞf ~ ðΩÞ ~ ¼ dΩK ðΩÞf ðΩÞ 1 þ βðΩÞ N 1 ðΩÞ dΩ

Eqn (22) yields ~ ¼ δij ~ ij ðΩÞ P

M

b n k M ki1 b n l M lj 1 N2

~ ¼ Mik ðδkl ~ ij ðΩÞ In the case of a pure rotation, P

¼ MT and N ¼ 1. Putting all this together yields Z ~kij ¼ dΩK fMik Pkl M Tlj ¼ Mik kkl M Tlj : 1

(27)

β ¼ βðΩÞ ¼ 1 þ ð1

Rm Þ

s ; b

(36)

~ are found in Appendices. is a function of bs ðΩÞ. Examples of computing k Previously we mentioned that when βε < 1, the Liu-Elsworth formula produces negative permeabilities. In Eqn (35) if βðN 1 1Þ < 1 the same problem occurs. To overcome this problem,

(28)

our proposal is the replacement � �� 1 þ β N 1 1 →H 1 þ β N 1 1 ;

5.2.2. No change in fracture aperture or spacing Generalizing the preceding section slightly, suppose that u describes a transformation that does not change b or s. Suppose that b and s are ~ and functions of Ω. The transformation causes Ω→Ω,

(37)

where � HðxÞ ¼

(29)

x 0

for ​ x � 0 for ​ x < 0

(38)

This does not preclude permeability from decreasing in response to a compressional strain, but simply limits the reduction so permeability is never negative.

which is similar to Eqn (25): after transformation the value of b3 = 12s ~ is just the value of b3 f= 12s at the (modified by f) at a certain angle Ω angle that transformed into this angle. Unlike a pure rotation, dΩ, may change under the action of such a ~ However, without any changes in fracture transformation: dΩ→dΩ. aperture and spacing, clearly ~ K~ ðΩÞ ~ ~f ðΩÞ ~ ¼ dΩK ðΩÞf ðΩÞ : dΩ

(34)

In this expression, M is given by Eqn (20), N by Eqn (23) and

b nk b n l ÞMTlj , since

This is just a rotation of the permeability tensor, as expected.

3 ~ b ~ ~ ~ b3 ðΩÞf ðΩÞ ¼ ðΩÞf ðΩÞ ; 12~s 12s

��3 1 :

Substituting this into Eqn (24) and using the transformation of the projection tensor in Eqn (27) yields ! Z b n k M ki1 b n l M lj 1 ��3 ~kij ¼ δij : (35) dΩ K f � 1 þ β N 1 1 N2

(26)

b ~ b ~ : n i ðΩÞ ~ n j ðΩÞ ~

(32)

where N is the normalization in Eqn (23). Refer to Fig. 2. The preceding paragraph has assumed that the rock and fracture deform identically. Following Liu and Elsworth, it is reasonable to introduce a modulus-reduction factor, encoded in β, and assume that the fracture aperture after deformation is �� ~ b¼b 1þβ N 1 1 : (33)

(25)

αÞ :

(31)

bb n T →bMb nT ;

6. Case studies Before presenting the case studies, it is worth mentioning that our proposal adds negligible computational expense to numerical models. Although formulae such as Eqn (35) and (E.3) look substantially more complicated than their Liu-Elsworth cousins (8), their impact on com­ puter runtime has not been discernible in any of our simulations.

(30)

since both sides are equal to b3 f=12sðΩÞ by Eqn (29). 5.2.3. Changes in fracture aperture and spacing Now let us perform similar motivational arguments to Eqns (10) and 5

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

because the fracture planes are rotated by the shear strain, as in the top ~xy is positive, because a gradient of drawing in the schematic of Fig. 3. k porepressure that drives fluid towards the positive y direction will also drive fluid to the right in the positive x direction (fluid will flow along the green arrow in Fig. 3). If the strain is too large, then βðN 1 1Þ < 1 and Eqn (37) zeroes the permeability tensor. In stark contrast, the isotropic case develops very large and negative kxy . Physically this is because the application of shear strain opens fractures lying along the line y ¼ x, as depicted in the schematic of Fig. 3. Similar effects were also found by Liu et al.73 Hence, kxy scales as ε3 for ε larger than about 0.01. kxy is negative, since a gradient of por­ epressure that drives fluid towards the positive y direction will also drive fluid to the left (roughly along the green arrow in Fig. 3) along these open fractures lying along y ¼ x. For the isotropic case, the exact direction of fluid flow in response to ~yy =k ~xy . a porepressure gradient in the y direction is dictated by the ratio k ~yy is nonzero for ε ¼ 0, and scales as ε2 for large ε. Eqn (E.7) reveals that k Therefore, the flow direction will be predominantly in the y direction for small ε, but as ε becomes larger, the flow direction will rotate until it approaches a 45∘ angle to the y axis.

Fig. 3. Schematic views: (a) an orthogonal distribution of fractures; (b) the orthogonal distribution after applying a pure shear; (c) an isotropic distribution of fractures; (d) the isotropic distribution after applying a pure shear, where the red-colored fractures have enhanced aperture. In each case the green arrow indicates the fluid-flow direction in response to a porepressure gradient in the y direction. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

6.2. Permeability change around a longwall goaf Extraction of coal using the longwall mining method creates large voids underground. The weight of the overlying rock means it collapses into the voids, creating the “longwall goaf”. The permeability of the rock in and around the goaf area is of paramount importance in controlling fluid flow in the region. In this section we calculate the permeability in and around the goaf using the original Liu-Elsworth formula and our proposal. We demon­ strate that both formulations predict virtually identical vertical perme­ ~zz . However, the Liu-Elsworth formula always keeps the offability k

6.1. Shear strains in 2D Consider two 2D situations. Firstly, orthogonal fractures aligned with the x and y axes with bx ¼ by ¼ 10 4 m, sx ¼ sy ¼ ð1 =12Þ m. Sec­ ondly, fracture planes isotropically distributed with aperture b0 ¼ 10 4 m and spacing s0 ¼ ðπ =24Þ m. According to Eqn (A.4) and Eqn (C.2), these distributions both result in kxx ¼ kyy ¼ 10

12

m2 ;

(39)

diagonal permeability components zero, while our proposal yields offdiagonal components that are of similar magnitude to the on-diagonal components. In this section, results from a simple 3D finite-element simulation of this process are presented. The model geometry is shown in Fig. 5. A single longwall panel is modelled, which is 400 m deep. It is 1 km long, 300 m wide and 3 m thick. Because of symmetry along the panel’s long axis, only 150 m of the panel needs to be modelled, with appropriate boundary conditions on the symmetry axis. (The half-panel sits in 0 � x � 150, 0 � y � 1000 and 400 � z � 397.) Only the rock

(approximately 1 Darcy) and kxy ¼ 0. Suppose that Rm ¼ 0:5. Then β ¼ 418 in the orthogonal-fracture case, and β ¼ 655 in the isotropic case. Apply a pure shear deformation, viz Eqn (18) with εxx ¼ 0 ¼ εyy and εxy ¼ ε ¼ εyx , as shown in the schematic of Fig. 3. ~xy , which is given by Eqn (A.11) and Eqn (E.4). The LiuConsider k ~xy ¼ 0. The results are shown in Fig. 4. Elsworth formalism always sets k ~xy only develops In the orthogonal fracture-plane case, a non-zero k

Fig. 4. resulting from a pure shear in 2D. Left: orthogonal fracture planes originally aligned with the x and y axes. Right: isotropic distribution of fracture planes. Note that the negative permeability component is plotted in the right-hand graph. 6

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

Fig. 5. Finite element mesh of the coal-mine model. Left: plan view from above, where the green rectangle shows the area of coal that will be excavated from the underlying coal seam. Top-right: section view along the panel, where the thin green line at the bottom shows the coal to be excavated. Bottom-right: 3D view, as seen from a little below the model. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

strata lying above the mine are modelled in this simulation ( 400 � z � 0). The finite-element mesh consists of fairly coarse hexahedra. In the panel region there are 17 elements lying along the 1 km long axis, while across the panel the elements have size 10 m. The vertical mesh size varies from 3 m near the mining seam to 70 m near the top surface. The extraction of coal is modelled by successively removing finiteelements from the longwall panel volume, starting at y ¼ 0 and mov­ ing step-by-step towards y ¼ 1 km. (In practice, this is achieved by successively setting the element’s Young’s modulus to zero.) At each step, an equilibrium solution is sought. Because only the rock strata above the longwall panel are simulated, there is no floor material to prevent the roof from collapsing to z ¼ ∞, so a constraint is placed upon the model that limits uz � 3 m. The model is solved using MOOSE.80,81 The rock strata are modelled by a single rock mass with uniform properties (for simplicity, there are no layers, such as sandstone, claystone, etc, that are typical around real coal mines). Despite the simplistic geometry, rather complicated elasto-plasticity is needed to suitably model the caving process. The rock’s Young’s modulus is 8 GPa and the Poisson’s ratio is 0.25. Mohr-Coulomb plasticity82,83 is used with a cohesion of 3 MPa, a friction angle of 37∘ and a dilation angle of 15∘. This is truncated by a Ran­ kine-like84 tensile-yield cap, where the maximum extensional principle stress cannot exceed the tensile strength of 1 MPa. In addition to these plasticity models, the rock is assumed to contain horizontal weak planes that are governed by nonassociated Coulomb friction with a tensile cutoff. These have shear strength of 0.1 MPa, a friction angle of 20∘, a dilation angle of 10∘ and a tensile strength of 0. The implementation of plasticity in MOOSE is quite sophisticated85,86 and allows the combi­ nation of any number of yield surfaces and flow rules. It is common in more realistic models to employ softening laws for many of these moduli, particularly the cohesive strengths and friction angles, to cap­ ture rock-mass weakening as it deforms. In this simple study, only per­ fect plasticity is used.

The deformation and plastic strains resulting from excavating the longwall mine are shown in Fig. 6. In this model, there is significant rock breakage in the caved zone — via Mohr-Coulomb shear failure and Rankine tensile failure — as shown in (a) and (b). There is less failure of the weak planes — as shown in (c) and (d), although it is more uniformly distributed throughout the overburden. Readers should understand that quantitatively the results may be significantly different in real coal mines, with different rock properties and more realistic geology. The in-situ horizontal permeability is assumed to be 10 14 m2 and the in-situ vertical permeability is 10 15 m2. The horizontal fracture spacing is 1 m and the vertical fracture spacing 0.1 m (these quantities were chosen for convenience to maintain the horizontal-to-vertical ratio). The modulus reduction ratio is assumed to follow Rm ¼ 0:4 þ 0:6

z þ 400 ​ for ​ z < 200

200 ;

(40)

and Rm ¼ 1 for z � 200. Physically, this means that the rock close to the seam (ze 400) is more broken — it has a lower rock mass rating5,73,74,77 — than the rock further up in the overburden. Two cases are studied: the first with fracture planes oriented orthogonally to the axes using the Liu-Elsworth formula; the second with a transverse-isotropic distribution of fracture planes using our proposal (see Appendix E). ~zz on a vertical slice across the excavation. Evidently, Fig. 7 shows k there is little difference in the results using the Liu-Elsworth formula or the formula proposed in this paper. Qualitatively, the results are as ex­ pected from field observations13,15,21,32,87–89 and physical models.16,18,22 There is a region of enhanced permeability above the longwall panel (in the caving and fractured zone). The zone of enhanced permeability extends horizontally outside the region of the longwall panel due to horizontal strains. There is a “lobe” of enhanced perme­ ability near the edges of the longwall panel that extends into the roof. Readers should understand that quantitatively the results may be

7

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

Fig. 6. Deformation and plastic strain on a vertical slice across the excavation. The longwall goaf is shown in green. The white arrows show the displacement. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7. on a vertical slice perpendicular to the longwall’s long axis. Left: Using the Liu-Elsworth formula, Eqn (8). Right: Using a transverse-isotropic distribution of fracture planes and the formula proposed in this paper, Eqn (35). Both approaches give more-or-less identical results, as expected.

significantly different in real coal mines, with different rock properties and more realistic geology. It is expected that Liu-Elsworth and our proposal produce very ~zz are similar results. The diagonal components of permeability, such as k

diagonal infinitesimal strain the two formulae are identical. ~zx on the same vertical slice, according to the formula Fig. 8 shows k proposed in this paper. The Liu-Elsworth formula always maintains ~zx ¼ 0, however, Fig. 8 illustrates that the off-diagonal components of k

chiefly controlled by the diagonal components of strain, and with purely 8

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

~xz ¼ 0. Right: Using the formula proposed Fig. 8. on a vertical slice perpendicular to the longwall’s long axis. Left: Using the Liu-Elsworth formula that always sets k ~ in this paper, Eqn (35), where this permeability component is of similar magnitude to kzz of Fig. 7.

Fig. 9. Flow directions that would result from a uniform porepressure gradient in the z direction (displayed on a vertical slice perpendicular to the longwall’s long axis). The arrows are colored by the magnitude of permeability, which is proportional to flow speed in this case.

permeability can be quite large in real applications. Again, this is to be expected: Liu-Elsworth maintain the off-diagonal permeability components are zero, while in our case these off-diagonal components indicate fluid is flowing at an angle to the axes, which obviously occurs if dipping fractures are opened by shear strains. Fig. 9 shows the flow directions that would result from a por­ epressure gradient in the z direction. Fig. 7 has demonstrated that the ~zz ) is similar for both models. vertical flow velocity (dictated by k

orthogonal fracture planes subjected to infinitesimal diagonal strains along these three directions. In many applications our proposal yields very similar results to that predicted by Liu-Elsworth because the two formalisms are identical in the case of infinitesimal and diagonal strains. However, in situations where shear strains open fractures that are not aligned with the axes, the off-diagonal components of the permeability are important and our proposal allows computation of the complete permeability tensor. Our proposal may need modification in the case of cyclic loading or shear failure, both of which may alter the fracture-plane morphology, which we assume does not occur, but which might be incorporated using a damage-dependent f function. Our proposal agrees with lab and in-situ direct and indirect tests of permeability3,5–7,9,10,13,21,26,27,53,72,73,74 inasmuch as the original Liu-Elsworth formalism also agrees with these experimental data. Our permeability tensor transforms correctly under rotations, and hence describes permeability in any reference frame, including those where the stress and strain are non-diagonal.

However, in the zone immediately above the panel this is dwarfed by the ~zx for the transverselarge horizontal flows induced by the large k

isotropic case. In real mining scenarios, the porepressure gradient will not be purely in the z direction, especially around the rib areas, but Fig. 9 clearly shows the impact of the non-zero off-diagonal terms in the deformed permeability tensor. 7. Conclusions We have proposed a formula that enables the calculation of the permeability tensor of fractured rock subjected to arbitrary de­ formations. The rock may contain an arbitrary distribution of fractures, and our formula transforms correctly under rotations of the global co­ ordinate system. This generalizes the highly successful and useful LiuElsworth1 formula that is only applicable to rock containing three sets of

Acknowledgment This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

9

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

Appendix A. Example: Orthogonal planes oriented along the axes in 2D

Z

In the appendices δðxÞ is the standard Dirac delta distribution with property f ðxÞδðx

(A.1)

x0 Þdx ¼ f ðx0 Þ:

In 2D, where the unit normal is given by Eqn (13), the projection operator is � � 1 cos2 θ sin θ cos θ nj ¼ Pij � δij b ni b 2 sin θ cos θ 1 sin θ

(A.2)

Suppose that there exist two sets of fracture planes. The first has normal in the x direction with aperture bx and spacing sx . The second has normal in the y direction with aperture by and spacing sy . Then � � b3y b3 1 K ðθÞ ¼ x δðθÞ þ δ θ π : (A.3) 2 12sx 12sy Performing the integration in Eqn (16) with f ¼ 1 using the distributions of Eqn (A.3) yields Z π � b3y kxx ¼ dθK ðθÞ 1 cos2 θ ¼ ; 12sy 0 Z π kxy ¼ dθK ðθÞð cos θ sin θÞ ¼ 0;

(A.4)

0

Z

π

kyy ¼

dθK ðθÞ 1 0

� b3 sin2 θ ¼ x ; 12sx

which are the 2D versions of Eqn (5). A general linear 2D deformation of Eqn (18) yields the transformation matrix of Eqn (20) with components � � 1 þ εxx εxy Mij ¼ : εyx 1 þ εyy

(A.5)

This has inverse � 1 1 þ εyy M ij 1 ¼ εyx Δ

(A.6)

εxy 1 þ εxx

� ;

where Δ is the determinant of M. The transformed normal vector of Eqn (22) has components 1 b n~ i ¼ NΔ

� 1 þ εyy cos θ

εyx sin θ ;

� 1 þ εyy cos θ

εyx sin θ

�

(A.7)

εxy cos θ þ ð1 þ εxx Þsin θ ;

with ðNΔÞ2 ¼

�2

þ ð1 þ εxx Þsin θ

εxy cos θ

�2

(A.8)

;

to ensure unit normal. In the case at hand, only the results at θ ¼ 0 and θ ¼ π =2 are needed, because of the Delta functions in K . Evaluating at θ ¼ 0 yields � 1 b n i ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ εyy ; εxy ; ​ when ​ θ ¼ 0 ~ 2 1 þ εyy þ ε2xy

N � Nθ¼0

(A.9)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 þ εyy þ ε2xy � ¼ ; ​ when ​ θ ¼ 0 : ð1 þ εxx Þ 1 þ εyy εxy εyx

~ follows immediately from Eqn (A.7), but to explain this physically, note that the tangent vector when θ ¼ 0 is just t ¼ ð0; 1Þ, so the trans­ (The b n formation yields ~t ¼ ðεxy ; 1 þ εyy Þ. This is shown in Fig. 2.) Evaluating at θ ¼ π=2 yields 1 b n i ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ ð1 þ εxx Þ2 þ ε2yx

N � Nθ¼π=2

�

εyx ; 1 þ εxx : ​ when ​ θ ¼ π

, 2 (A.10)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , ð1 þ εxx Þ2 þ ε2yx � ¼ ; ​ when ​ θ ¼ π 2 : ð1 þ εxx Þ 1 þ εyy εxy εyx

This may be explained physically be considering the tangent vector t ¼ ð1; 0Þ. 10

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

Finally, inserting these results into Eqn (35) and performing the integration over θ yields � � ~kxx ¼ kxx 1 þ βy N 1 θ¼π =2

1 ​ þ kyy 1 þ βx N θ¼0

� � ~kxy ¼ kxx 1 þ βy N 1 θ¼π =2 1 ​ þ kyy 1 þ βx N θ¼0

� � ~kyy ¼ kxx 1 þ βy N 1 θ¼π =2

1 ​ þ kyy 1 þ βx N θ¼0

��3 1 ��3 1

ð1 þ εxx Þ2 ð1 þ εxx Þ2 þ ε2yx

ε2xy 1 þ εyy

�2

þ ε2xy

;

��3 ε ð1 þ ε Þ yx xx 1 ð1 þ εxx Þ2 þ ε2yx � ��3 εxy 1 þ εyy 1 ; �2 1 þ εyy þ ε2xy ��3 1 ��3 1

(A.11)

ε2yx ð1 þ εxx Þ2 þ ε2yx

�2 1 þ εyy : �2 1 þ εyy þ ε2xy

Here βx ¼ 1 þ ð1 Rm Þsx =bx and βy ¼ 1 þ ð1 Rm Þsy =by . This is the 2D version of the Liu-Elsworth Eqn (8) generalised to arbitrary finite strain. The nonzero off-diagonal components arise from rotation of the fracture planes. If εxy ¼ 0 ¼ εyx , the permeabilities read � ~kxx ¼ kxx 1 þ βy εyy 3 ; ~kxy ¼ 0 ; (A.12) ~kyy ¼ kyy ð1 þ βx εxx Þ3 ; which are exactly the Liu-Elsworth results. To first order, the shear strains εxy and εyx have no effect whatsoever on the diagonal components of permeability. Physically this is because the shear strain εxy may be thought of an expansion along the x þ y direction with a corresponding contraction along the x y direction. The expansion causes an increase in bx while the contraction causes a decrease, with a zero net change. For finite strain, the shear strain causes a change in permeability. Physically this is because a finite shear strain does cause lengths along the x and y axes to change, even though this is a second-order effect. Appendix B. Example: Orthogonal planes oriented at an angle to the axes in 2D Suppose that there exist two sets of fracture planes. The first has normal b n ¼ ðcosα; sinαÞ with aperture bα and spacing sα . The second has normal � � �� � � b with aperture bαþπ=2 and spacing sαþπ=2 . This is illustrated in Fig. 10. Then n ¼ cos α þ12 π ; sin α þ12 π K ðΩÞ ¼

b3α δðθ 12sα

αÞ þ

b3αþπ=2 12sαþπ=2

� δ θ

α

1 π 2

� (B.1)

:

Fig. 10. Two orthogonal sets of fracture planes oriented at an angle to the ðx; yÞ axes. The red set have aperture bα and spacing sα . The blue set have aperture bαþπ=2 and spacing sαþπ=2 .

Performing the integration in Eqn (16) using the distributions of Eqn (B.1) yields

11

A. Wilkins and Q. Qu

Z

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

� � b3αþπ=2 � b3 1 cos2 θ ¼ α sin2 α þ sin2 α þ π ; 2 12sα 12sαþπ=2

π

kxx ¼

dθK ðθÞ 1 0

Z

π

kxy ¼

dθK ðθÞð cos θ sin θÞ ¼ 0

Z

dθK ðθÞ 1 0

�

� � � � 1 1 sin α þ π cos α þ π ; 2 2 12sαþπ=2 b3αþπ=2

(B.2)

� � b3αþπ=2 � b3 1 sin2 θ ¼ α cos2 α þ cos2 α þ π : 2 12sα 12sαþπ=2

π

kyy ¼

b3α sin α cos α 12sα

A rotation to the ð~ x; ~ yÞ reference frame indicated in Figure 10 reveals that these forms are reasonable: � � �� �� � k~x~x k~x~y kxx kxy cos α sin α cos α sin α ¼ : k~x~y k~y~y kxy kyy sin α cos α sin α cos α

(B.3)

Explicitly: k~x~x ¼

b3αþπ=2 12sαþπ=2

; (B.4)

k~x~y ¼ 0 ; k~y~y ¼

b3α ; 12sα

which are the forms expected from Eqn (5) and (A.4). Appendix C. Example: A uniform distribution of planes Suppose that there exists a distribution of planes (in 2D or 3D), with the aperture, b0 , and spacing, s0 , both being independent of the fracture normal direction. Then K ðΩÞ ¼

b30 12s0

(C.1)

Performing the integration in Eqn (16) using the distributions of Eqn (C.1) yields, in 2D, Z π � πb3 kxx ¼ dθK 1 cos2 θ ¼ 0 ; 24s 0 0 Z π kxy ¼ dθK ð cos θ sin θÞ ¼ 0; 0

Z

π

kyy ¼

dθK 1 0

(C.2)

� π b3 sin2 θ ¼ 0 : 24s0

Apart from the trivial factor of π=2 resulting from the integration, this tensor is indistinguishable from Eqn (A.4) if bx ¼ by and sx ¼ sy , and Eqn (B.2) if bα ¼ bαþπ=2 and sα ¼ sαþπ=2 . The computation in 3D yields kxx ¼ kyy ¼ kzz ¼

πb30 9s0

(C.3)

;

with the off-diagonal components being zero. Appendix D. Example: Orthogonal planes oriented along the axes in 3D In 3D, where the unit normal is given by Eqn (13), the projection operator is 0 1 sin2 φ sin θ cos θ sin φ cos φ cos θ 1 sin2 cos2 θ 2 2 2 @ nj ¼ Pij � δij b ni b 1 sin sin θ sin φ cos φ sin θ A sin φ sin θ cos θ sin φ cos φ cos θ sin φ cos φ sin θ sin2 φ

(D.1)

Suppose that there exist three sets of fracture planes. The first set has normal in the x direction with aperture bx and spacing sx . The second set has normal in the y direction with aperture by and spacing sy . The third set has normal in the z direction with aperture bz and spacing sz . This is the case considered by Liu and Elsworth. Then , � � � � � � b3 b3 b3 1 1 1 K ðΩÞ ¼ x δðθÞδ φ π þ y δ θ π δ φ π þ z δðφÞ ðπ sin φÞ : (D.2) 2 2 2 12sx 12sy 12sz The π sin φ in the denominator is to correctly account for the measure at φ ¼ 0. 12

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

Performing the integration in Eqn (16) using the distributions of Eqn (D.2) yields kxx ¼

b3y b3 þ z ; 12sy 12sz

kyy ¼

b3 b3x þ z ; 12sx 12sz

kzz ¼

b3y b3x þ ; 12sx 12sy

(D.3)

with all other components being zero. These forms have already been written in Eqn (5). A general linear 3D deformation is ux ¼ εxx x þ εxy y þ εxz z ; uy ¼ εyx x þ εyy y þ εyz z ; uz ¼ εzx x þ εzy y þ εzz z :

(D.4)

Here ε is clearly the infinitesimal strain tensor, but because it is not symmetric, rigid-body rotations are included. The transformation matrix of Eqn (20) has components 1 0 1 þ εxx εxy εxz @ εyx 1 þ εyy εyz A ; Mij ¼ εzx εzy 1 þ εzz which has inverse � 0 1 þ εyy ð1 þ εzz Þ εyz εzy 1@ 1 M ij ¼ εyz εzx εyx ð1 þ �εzz Þ Δ εzy εyx 1 þ εyy εzx

�

(D.5)

1

εxz εzy εxy ð1 þ εzz Þ εxy εyz εxz 1 þ εyy ð1 þ εxx Þð1 þ εzz Þ εxz εzx εyx εxz ð1 þ ε�xx Þεyz A ; εzx εxy ð1 þ εxx Þεzy ð1 þ εxx Þ 1 þ εyy εxy εyx

(D.6)

where Δ is the determinant of M. The transformed normal vector of Eqn (22) is only needed for ðθ; φÞ being ð0; π =2Þ, ðπ =2; π =2Þ and ðθ; 0Þ, corre­ sponding to n ¼ ð1; 0; 0Þ, n ¼ ð0; 1; 0Þ and n ¼ ð0; 0; 1Þ, respectively. These are � b n i ¼ M xi1 Nx ​ for ​ n ¼ ð1; 0; 0Þ ; ~ . b n i ¼ M yi1 Ny ​ for ​ n ¼ ð0; 1; 0Þ ; ~ (D.7) � 1 b n i ¼ M zi Nz ​ for ​ n ¼ ð0; 0; 1Þ : ~ where N2i ¼

P�� 1 2 �Mij j ensures unit normalization. j

The arbitrary deformation given in Eqn (D.4) may be applied to the orthogonally-fractured rock of Eqn (D.2). The algebra is somewhat tedious but elementary and the result is ! X b3 ��3 M ai1 M aj1 a kij ¼ 1 þ βa N a 1 1 δij : (D.8) 12sa N 2a a2ðx;y;zÞ This is our finite-strain version of Eqn (12). Here βa ¼ 1 þ ð1 Rm Þsa =ba for a ¼ x, y or z. To leading order in ε, the Liu-Elsworth forms of Eqn (12) are recovered: kxx ¼

�3 b3y b3 1 þ βy εyy þ z ð1 þ βz εzz Þ3 ; 12sy 12sz

kyy ¼

b3 b3x ð1 þ βx εxx Þ3 þ z ð1 þ βz εzz Þ3 ; 12sx 12sz

kzz ¼

�3 b3y b3x 1 þ βy εyy ; ð1 þ βx εxx Þ3 þ 12sx 12sy

(D.9)

kxy ¼ 0 ; kxz ¼ 0 ; kyz ¼ 0 : In particular, the shear strains εxy , εxz and εyz have no effect whatsoever on the permeability, for the same reasons as in 2D (Section A).

13

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

Appendix E. Example: A transversely isotropic distribution in 3D Suppose that there exists a transversely isotropic distribution of planes with the following characteristics. There is a uniform distribution of fracture planes with normals in the ðx; yÞ plane has aperture bh and spacing sh . There is also a distribution of fracture planes with normals in the z direction with aperture bz and spacing sz . Then � � � b3 b3 1 K ðΩÞ ¼ h δ φ π þ z δðφÞ sin φ : (E.1) 2 sh sz Performing the integration in Eqn (16) using the distributions of Eqn (E.1) yields k0xx ¼ k0yy ¼ k0zz

¼

π b3h 12sh

π b3h 24sh

πb3z

þ

12sz

; (E.2)

;

with the off-diagonal components being zero. For the general linear 3D deformation of Eqn (D.4) the integration in Eqn (35) cannot be performed to produce a useful formula. Instead, assume that the strain is small and the expand to lowest order in strain. Specifically, the projection operator of Eqn (27) using 3D unit normal (13) and transformation inverse (D.6) may be easily expanded to first order. The coefficient of β, that is N 1 1 may be similarly expanded to first order. However, because β can be large, βðN 1 1Þ can be very large compared with 1 (which is the reason that the Liu-Elsworth formula produces order-ofmagnitude permeability changes), all orders of β must be retained. The result is kxx ¼

πb3z 12sz

ð1 þ βz εzz Þ3

� � � 3 � π b3 1 2 þ εxx εyy þ βh εxx þ ε2xx ε2xy þ ε2yx þ 3εyy ε2yy þ h 8 12sh 4 �� � � � � 3 þ β2h 5ε3xx þ 4ε2yx þ ε2xy 4 9εyy 2εxy εyx 4 þ εyy þ 7ε2yx εyy þ 20ε2yy 64 � � � �� 7ε2xy þ 2εxy εyx þ 9ε2yx 5ε3yy þ ε2xx 4 þ εyy þ εxx 8 þ εyy εyy � 1 þ β3h 7ε4xx 128

3ε4xy

þε3xx 5 þ 2εyy � þε2xx

�

6ε3xy εyx þ 3ε4yx þ 15ε2yx εyy þ 12ε2yx ε2yy þ 35ε3yy

3ε2xy εyy

� � 5 þ 6εyy þ 6εxy εyx ε2yx

� � 5 þ εyy εyy

πb3h 1 12sh 4

�

3 8

εxy þ εyx þ βh ð

3 � þ β2h ε3xy þ ε2xx 13εxy 64

� 3εyx þ 3ε2xy εyx þ 3εxy ε2yx þ ε3yx

� � 1 3 βh 3ε3xy þ ε3xx 22εxy 8εyx þ 3ε2xy εyx 128 � � � � � þε3yx 3 þ 6εyy εxy ε2yy 15 þ 8εyy þ εyx ε2yy ��

� � 1 þ 2εyy εyx

1 þ 2εxx Þεxy þ

3εxy ε2yy þ 13εyx ε2yy þ 2εxx εxy þ εyx

5 þ 4εyy

� 18εyy

� � ���i 2εyy ε2yy þ 9ε2yx 1 þ 2εyy

�

kxy ¼

þεxy

7ε4yy

�� � � � 12ε2xy þ 6εxy εyx þ 9 2ε2yx þ εyy þ εxx 18εxy εyx þ ε2xy 9 � þ 15

(E.3)

� þ 6εxx ε3xy þ 2ε2xy εyx þ εyx εyy

�

� 8εxy εyy

4 þ 3εyy �

8εyx εyy

��

2 xy yx

3 þ 2εyy þ 3ε ε

3 þ 4εyy

� � � 15 þ 22εyy þ 3ε2xx εyx � � 3 þ 2εyy þ εxy ε2yx þ

(E.4)

� �

5 þ 2εyy

� ���i 3 þ εyy εyy

14

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

kxz ¼ εzx

πb3z 12sz

ð1 þ βz εzz Þ3

� � � πb3 1 3 þ h εxz þ βh 3εxx εxz þ εxz εyy þ εxy þ εyx εyz 8 12sh 2 � 3 2� 2 β 5εxx εxz þ ε2xy εxz þ εxz ε2yx þ εxz ε2yy þ 2εyx εyy εyz 16 h � � �� þ2εxx εxz εyy þ εxy þ εyx εyz þ 2εxy εxz εyx þ εyy εyz � 1 þ β3h 35ε3xx εxz þ 9εxz ε2yx εyy þ 5εxz ε3yy þ 3ε3xy εyz þ 3ε3yx εyz þ 15εyx ε2yy εyz 128 � � � þ9ε2xy εxz εyy þ εyx εyz þ 15ε2xx εxz εyy þ εxy þ εyx εyz � � þ3εxx 5ε2xy εxz þ 10εxy εxz εyx þ 5εxz ε2yx þ 3εxz ε2yy þ 6εxy εyy εyz þ 6εyx εyy εyz þ

(E.5)

� ��i þ3εxy 6εxz εyx εyy þ 3ε2yx εyz þ 5ε2yy εyz kyx ¼ kxy

(E.6)

kyy ¼ kxx ​ but ​ with ​ x ↔ y :

(E.7)

kyz ¼ kxz ​ but ​ with ​ x ↔ y :

(E.8)

kzx ¼ kxz

(E.9) (E.10)

kzy ¼ kyz � � π b3 3 kzz ¼ h 1 þ βh εxx þ εyy 2 12sh � 3 � þ β2h 3ε2xx þ 2εxx εyy þ 3ε2yy þ ε2xy þ 2εxy εyx þ ε2yx 8 �� �� 1 3 þ βh εxx þ εyy 5ε2xx 2εxx εyy þ 5ε2yy þ 3ε2xy þ 6εxy εyx þ 3ε2yx 16

(E.11)

This paper is accompanied by a Mathematica notebook that derives these results.

References

13 Adhikary DP, Guo H. Modelling of longwall mining-induced strata permeability change. Rock Mech Rock Eng. 2015;48:345–359. https://doi.org/10.1007/s00603014-0551-7. 14 Adhikary D, Poulsen B. Estimating the height of mining induced connective fracturing. In: Smiley P, ed. 50th US Rock Mechanics/Geomechanics Symposium; June 27-30; Huston, Texas. ARMA; 2016:10. 15 Meng Z, Shi X, Li G. Deformation, failure and permeability of coal-bearing strata during longwall mining. Eng Geol. 2016;208:69–80. https://doi.org/10.1016/j. enggeo.2016.04.029. 16 Wang W, Jiang T, Faybishenko B, Wang Z, Hu W, Zhao Q. Closure of fracture due to cover stress Re-establishment after coal mining. Geotech Geol Eng. 2016;34: 1525–1537. https://doi.org/10.1007/s10706-016-0059-x. 17 Xu D, Peng S, Xiang S, He Y. A novel caving model of overburden strata movement induced by coal mining. Energies. 2017;10:476. https://doi.org/10.3390/ en10040476. 18 Zhang Y, Xu Y, Wang K, et al. The fracturing characteristics of rock mass of coal mining and its effect on overlying unconsolidated aquifer in Shanxi, China. Arabian J. Geosci. 2018;11:666. https://doi.org/10.1007/s12517-018-4034-0. 19 Poulsen B, Adhikary D, Guo H. Simulating mining-induced strata permeability changes. Int J Rock Mech Min Sci. 2018;237:208–216. https://doi.org/10.1016/j. enggeo.2018.03.001. 20 Khanal M, Guo H, Adhikary D. 3D numerical study of underground coal mining induced strata deformation and subsequent permeability change. Geotech Geol Eng. 2019;37:235–249. https://doi.org/10.1007/s10706-018-0605-9. 21 Guo H, Adhikary DP, Craig MS. Simulation of mine water inflow and gas emission during longwall mining. Rock Mech Rock Eng. 2009;42:25–51. https://doi.org/ 10.1007/s00603-008-0168-9. 22 Hu S, Jin Z, Zhang K, et al. Analyzing water resource protection techniques when mining shallow coal seam. Electron J Geotech Eng. 2016;21:1131–1140. 23 Zhang J, Peng S. Water inrush and environmental impact of shallow seam mining. Environ Geol. 2005;48:1068–1076. https://doi.org/10.1007/s00254-005-0045-8. 24 Wang WX, Sui WH, Faybishenko B, Stringfellow WT. Permeability variations within mining-induced fractured rock mass and its influence on groundwater inrush. Environ. Earth Sci. 2016;75:326. https://doi.org/10.1007/s12665-015-5064-5. 25 Adhikary D, Guo H. Assessment of mine water ingress and impact of mining on ground and surface water system. In: International Conference on NexGen Technologies

1 Liu J, Elsworth D. Three-dimensional effects of hydraulic conductivity enhancement and desaturation around mined panels. Int J Rock Mech Min Sci. 1997;34:1139–1152. https://doi.org/10.1016/S1365-1609(97)80067-6. 2 Shen B, Poulsen B, Qu Q. Overburden strata movement and stress change induced by longwall mining. In: 45th US Rock Mechanics/Geomechanics Symposium; 26-29 June, 2011. San Francisco, California: American Rock Mechanics Association; 2011:8. 3 Guo H, Qu Q, Adhikary D. Overburden response to longwall mining. In: Peng S, ed. Advance in Coal Mine Ground Control. Elsevier; 2017:111–155. 4 Adhikary D, Poulsen B, Khanal M. Assessment of longwall mining induced connective fracturing of overburden strata. ACARP. 2019. Report C24020. 5 Liu J, Elsworth D, Brady BH. Linking stress-dependent effective porosity and hydraulic conductivity to RMR. Int J Rock Mech Min Sci. 1999;36:581–596. https:// doi.org/10.1016/S0148-9062(99)00029-7. 6 Adhikary DP, Guo H. A coupled Cosserat two-phase double porosity flow model. MODSIM. 2005:1189–1195. 7 Adhikary DP, Guo H. A continuum model for simulating mine water inflow and gas emission. In: Proceedings of FED SM2008. ASME Fluids Engineering Conference; 2008:1–10. https://doi.org/10.1115/FEDSM2008-55044. 8 Yuan L, Guo H, Shen B, Qu Q, Xue J. Circular overlying zone at longwall panel for efficient methane capture of multiple coal seams with low permeability. J China Coal Soc. 2011;36(3):357–365. 9 Adhikary DP, Guo H. A coupled model of two-phase diffusion and flow through deforming porous cosserat media. Defect Diffusion Forum. 2010;297:281–293. https ://doi.org/10.4028/www.scientific.net/DDF.297-301.281. 10 Adhikary DP. Numerical simulation of mine water inflow. In: ISRM International Symposium EUROCK. 1–12. 2012. 11 Guo H, Yuan L, Shen B, Qu Q, Xue J. Mining-induced strata stress changes, fractures and gas flow dynamics in multi-seam long wall mining. Int J Rock Mech Min Sci. 2012;54:129–139. https://doi.org/10.1016/j.ijrmms.2012.05.023. 12 Adhikary D, Guo H. Measurement of longwall mining induced strata permeability. Geotech Geol Eng. 2014;32(3):617–626. https://doi.org/10.1007/s10706-014-97378.

15

A. Wilkins and Q. Qu

26 27 28

29

30

31

32 33 34 35 36 37

38

39 40 41 42 43 44 45 46

47

48 49 50

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

for Mining and Fuel Industries (NxGnMiFu-2017); 15-17 February 2017. New Delhi India: NxGnMiFu; 2017:1–5. Adhikary D, Wilkins A. Reducing the impact of longwall extraction on groundwater systems. ACARP. 2012. Report C18016. Adhikary D. Simulation of mining impact on water resources. In: EUROCK 2015 & 64th GEOMECHANICS COLLOQUIUM; 7-10 October 2015; Salzburg, Austria. 2015. Herron NF, Frery E, Crosbie R, et al. Observations Analysis, Statistical Analysis and Interpolation for the Hunter Subregion. Product 2.1-2.2 for the Hunter Subregion from the Northern Sydney Basin Bioregional Assessment. Australia: Department of the Environment and Energy, Bureau of Meteorology, CSIRO and Geoscience Australia; 2018. http://data.bioregionalassessments.gov.au/product/NSB/HUN/2.1-2.2. Herron NF, Peeters L, Crosbie R, et al. Groundwater Numerical Modelling for the Hunter Subregion. Product 2.6.2 for the Hunter Subregion from the Northern Sydney Basin Bioregional Assessment. Australia: Department of the Environment and Energy, Bureau of Meteorology, CSIRO and Geoscience Australia; 2018. http://data.bioregi onalassessments.gov.au/product/NSB/HUN/2.6.2. Herron NF, Macfarlane C, Henderson BL, et al. Impact and Risk Analysis for the Hunter Subregion. Product 3-4 for the Hunter Subregion from the Northern Sydney Basin Bioregional Assessment. Australia: Department of the Environment and Energy, Bureau of Meteorology, CSIRO and Geoscience Australia; 2018. http://data. bioregionalassessments.gov.au/product/NSB/HUN/3-4. Herron NF, Crosbie R, Peeters L, et al. Water Balance Assessment for the Hunter Subregion. Product 2.5 for the Hunter Subregion from the Northern Sydney Basin Bioregional Assessment. Australia: Department of the Environment and Energy, Bureau of Meteorology, CSIRO and Geoscience Australia; 2018. http://data. bioregionalassessments.gov.au/product/NSB/HUN/2.5. Guo H, Yuan L, Shen B, Qu Q, Xue J. Mining-induced strata stress changes, fractures and gas flow dynamics in multi-seam longwall mining. Int J Rock Mech Min Sci. 2012; 54:129–139. https://doi.org/10.1016/j.ijrmms.2012.05.023. Qu Q, Xu J, Wu R, Qin W, Hu G. Three-zone characterisation of coupled strata and gas behaviour in multi-seam mining. Int J Rock Mech Min Sci. 2015;78:91–98. Qu Q, Guo H, Todhunter C, Kerr H, Babic L, McGeachi R. Groundwater and gas behaviours. In: 24th World Mining Congress; 18-21 October 2016; Rio de Janeiro, Brazil. Rio de Janeiro Brazil World Mining Congress; 2016:17–26. Li Z. Defining Outburst-free Zones of a Protected Coal Seam in Multiple Seam Mining with Seam Gas Content Method. University of Wollongong Masters Thesis; 2016. Karacan CO. Reconciling longwall gob gas reservoirs and venthole production performances using multiple rate drawdown well test analysis. Int J Coal Geol. 2009; 80:181–195. https://doi.org/10.1016/j.coal.2009.09.006. Dougherty HN, Karacan CO, Goodman GVR. Reservoir diagnosis of longwall gobs through drawdown tests and decline curve analyses of gob gas venthole productions. Int J Rock Mech Min Sci. 2010;47:851–857. https://doi.org/10.1016/j. ijrmms.2010.04.010. Qu Q, Guo H, Shen B. Characterization of longwall strata behavior by 3D geotechnical modelling. In: Guo, Weijia, Shen, et al., eds. 3rd International Workshop on Mine Hazards Prevention and Control; 19-21 November 2013. Brisbane, Queensland: Atlantis Press; 2013:1–6. Xia T, Zhou F, Liu J, Hu S, Liu Y. A fully coupled coal deformation and compositional flow model for the control of the pre-mining coal seam gas extraction. Int J Rock Mech Min Sci. 2014;72:138–148. https://doi.org/10.1016/j.ijrmms.2014.08.012. Guo H, Todhunter C, Qu Q, Qin J. Longwall horizontal gas drainage through goaf pressure control. Int J Coal Geol. 2015;150:276–286. https://doi.org/10.1016/j. coal.2015.09.003. Guo H, Yuan L, Liang Y, Qu Q, Qin J, Xue S. Xie J Co-extraction of coal and methane. CIM Journal. 2015;6(11):5–14. Qin J, Yuan L, Guo H, Qingdong Q. Investigation of longwall goaf gas flows and borehole drainage performance by CFD simulation. Int J Coal Geol. 2015;150–151: 51–63. Qu Q, Guo H, Loney M. Analysis of longwall goaf gas drainage trials with surface directional boreholes. Int J Coal Geol. 2016;156:59–73. https://doi.org/10.1016/j. coal.2016.02.001. Qin J, Qu Q, Guo H. CFD simulations for longwall gas drainage design optimisation. In: The 9th International Symposium on Green Mining; 28-30 November 2016. University of Wollongong. University of Wollongong: Elsevier; 2016:1–6. Guo H, Todhunter C, Qu Q, Kerr H, Qin J. An integrated approach to optimisation of longwall gas drainage. In: 24th World Mining Congress - Mining in a world of innovation; 18-21 October 2016; Rio de Janeirovols 57–72. Guo H, Todhunter C, Qu Q, Kerr H, Qin J. An innovative drainage system for coal mine methane capture optimisation and abatement maximisation. In: Aziz N, ed. Coal Operators Conference; 11&12 Feb 2016; University of Wollongong. Wollongong: University of Wollongong; 2016:390–403. Qu Q, Guo H, Todhunter C, Kerr H, Babic L, RcGeachi R. Integrated field studies for optimal goaf gas drainage design in multi-seam longwall mining. In: Aziz N, ed. Coal Operators’ Conference; Feb 11-12 2016. the University of Wollongong. University of Wollongong; 2016:394–404. Rutqvist J, Stephansson O. The role of hydromechanical coupling in fractured rock engineering. Hydrogeol J. 2003;11:7–40. https://doi.org/10.1007/s10040-002-02415. Wu Y, Liu J, Elsworth D, Chen Z, Connell L, Pan Z. Dual poroelastic response of a coal seam to CO2 injection. Int. J. Greenh. Gas Control. 2010;4:668–678. https://doi.org/ 10.1016/j.ijggc.2010.02.004. Wu Y, Liu J, Elsworth D, Miao X, Mao X. Development of anisotropic permeability during coalbed methane production. J Nat Gas Sci Eng. 2010;2:197–210. https://doi. org/10.1016/j.jngse.2010.06.002.

51 Liu J, Chen Z, Elsworth D, Qu H, Chen D. Interactions of multiple processes during CBM extraction: a critical review. Int J Coal Geol. 2011;87:175–189. https://doi.org/ 10.1016/j.coal.2011.06.004. 52 Pan Z, Connell LD. Modelling permeability for coal reservoirs: a review of analytical models and testing data. Int J Coal Geol. 2012;92:1–44. https://doi.org/10.1016/j. coal.2011.12.009. 53 Wang S, Elsworth D, Liu J. A mechanistic model for permeability evolution in fractured sorbing media. J. Geophys. Res. B: Solid Earth. 2012;117:17pp. https://doi. org/10.1029/2011JB008855. 54 Chen Z, Liu J, Kabir A, Wang J, Pan Z. Impact of various parameters on the production of coalbed methane. SPE J. 2013;18:910–923. https://doi.org/10.2118/ 162722-PA. 55 Wang C, Wei MY, Jiang JD, Liu JS. Evolution of gas diffusivity in coals. In: IRSM International Symposium — 8th Asian Rock Mechanics Symposium. 2014. ISRM-ARMS82014-055. 56 Yang D, Qi X, Chen W, Wang S, Dai F. Numerical investigation on the coupled gassolid behavior of coal using an improved anisotropic permeability model. J Nat Gas Sci Eng. 2016;34:226–235. https://doi.org/10.1016/j.jngse.2016.06.058. 57 Peng S, Zhang J. Stress-dependent permeability. In: Engineering Geology for Underground Rocks. Berlin, Heidelberg: Springer; 2007. https://doi.org/10.1007/ 978-3-540-73295-2_8. 58 Zhang L. Aspects of rock permeability. Front Struct Civ Eng. 2013;7(2):102–116. https://doi.org/10.1007/s11709-013-0201-2. 59 Zou L, Tarasov BG, Dyskin AV, Adhikary DP, Pasternak E, Xu W. Physical modelling of stress-dependent permeability in fractured rocks. Rock Mech Rock Eng. 2013;46: 67–81. https://doi.org/10.1007/s00603-012-0254-x. 60 Zheng J, Zheng L, Liu HH, Ju Y. Relationships between permeability, porosity and effective stress for low-permeability sedimentary rock. Int J Rock Mech Min Sci. 2015; 78:304–318. https://doi.org/10.1016/j.ijrmms.2015.04.025. 61 Ma J. Review of permeability evolution model for fractured porous media. J. Rock Mech. Geotech. Eng. 2015;7(3):351–357. https://doi.org/10.1016/j. jrmge.2014.12.003. 62 Feng R, Chen S, Bryant S, Liu J. Stress-dependent permeability measurement techniques for unconventional gas reservoirs: review, evaluation, and application. Fuel. 2019;256:115987. https://doi.org/10.1016/j.fuel.2019.115987. 63 Kozeny J. Über kapillare Leitun des Wassers im Boden. Sitzungsber Akad. Wiss. Wien. 1927;12(2a):271–306. 64 Carman P. Fluid flow through a granular bed. Trans Inst Chem Eng. 1937;15, 150–17. 65 Oelkers EH. Physical and chemical properties of rocks and fluids for chemical mass transport calculations. Rev Mineral Geochem. 1996;34:131–191. 66 Liu HH, Rutqvist J. A new coal-permeability model, internal swelling stress and fracture-matrix interaction. Transp Porous Media. 2010;82(1):157–171. https://10 .1007/s11242-009-9442-x. 67 Bustin RM, Cui XJ, Chikatamarla L. Impacts of volumetric strain on CO2 sequestration in coals and enhanced CH4 recovery. AAPG (Am Assoc Pet Geol) Bull. 2008;92(1):15–29. https://doi.org/10.1306/08220706113. 68 Lyakhovsky V, Hamiel Y. Damage evolution and fluid flow in poroelastic rock. Izv Phys Solid Earth. 2007;43(1):13–23. https://doi.org/10.1134/S106935130701003X. 69 Oda M, Takemura T, Aoki T. Damage growth and permeability change in triaxial compression tests of Inada granite. Mech Mater. 2002;34:313–331. https://doi.org/ 10.1016/S0167-6636(02)00115-1. 70 Jiang T, Shao JF, Xu WY, Zhou CB. Experimental investigation and micromechanical analysis of damage and permeability variation in brittle rocks. Int J Rock Mech Min Sci. 2010;47:703–713. https://doi.org/10.1016/j.ijrmms.2010.05.003. 71 Witherspoon PA, Wang JSY, Iwai K, Gale JE. Validity of Cubic Law for fluid flow in a deformable rock fracture. Water Resour Res. 1980;16:1016–1024. https://doi.org/ 10.1029/WR016i006p01016. 72 Liu J, Chen Z, Elsworth D, Miao X, Mao X. Linking gas-sorption induced changes in coal permeability to directional strains through a modulus reduction ratio. Int J Coal Geol. 2010;83:21–30. https://doi.org/10.1016/j.coal.2010.04.006. 73 Liu J, Elsworth D, Brady BH, Muhlhaus HB. Strain-dependent fluid flow defined through rock mass classification schemes. Rock Mech Rock Eng. 2000;33:75–92. https://doi.org/10.1007/s006030050036. 74 Liu J, Elsworth D, Brady BH. A coupled hydromechanical system defined through rock mass classification schemes. Int J Numer Anal Methods Geomech. 1999;23: 1945–1960. https://doi.org/10.1002/(SICI)1096-9853(19991225)23:151945::AIDNAG423.0.CO;2-N. 75 Snow DT. Anisotropic permeability of fractured media. Water Resour Res. 1969;5: 1273–1289. https://doi.org/10.1029/WR005i006p01273. 76 Darcy H. Les fontaines publiques de la ville de Dijon. Paris: Dalmont; 1856. 77 Nicholson GA, Bieniawski ZT. A nonlinear deformation modulus based on rock mass classification. Int J Min Geol Eng. 1990;8:181–202. https://doi.org/10.1007/ BF01554041. 78 Panda BB, Kulatilake PHSW. Study of the effect of joint geometry parameters on the permeability of jointed rock. In: 35th U.S. Symposium on Rock Mechanics. USRMS; 1995:273–278. 79 Mase GE. Shaum’s Outline of Theory and Problems of Continuum Mechanics. New York: Schaum’s Outline Series McGraw-Hill; 1970. 80 Alger B, Andr�s D, Carlsen RW, et al. MOOSE Web page. https://mooseframework. org; 2019. Accessed 1 July 2019. 81 Gaston DR, Permann CJ, Peterson JW, et al. Physics-based multiscale coupling for full core nuclear reactor simulation. Ann Nucl Energy. 2015;84:45–54. 82 Doghri I. Mechanics of Deformable Solids: Linear and Nonlinear, Analytical and Computational Asepects. London: SpringerVerlag; 2000. https://doi.org/10.1007/ 978-3-662-04168-0.

16

A. Wilkins and Q. Qu

International Journal of Rock Mechanics and Mining Sciences 125 (2020) 104159

83 Abbo AJ, Lyamin AV, Sloan SW, Hambleton JP. A C2 continuous approximation to the Mohr–Coulomb yield surface. Int J Solids Struct. 2011;48:3001–3010. https://doi. org/10.1016/j.ijsolstr.2011.06.021. 84 Rankine W. On the Stability of Loose Earth. Philosophical Transactions of the Royal Society of London; 1857:147. 85 Adhikary DP, Jayasundara CT, Podgorney RK, Wilkins AH. A robust return-map algorithm for general multisurface plasticity. Int J Numer Methods Eng. 2017;109(2): 218–234. https://doi.org/10.1002/nme.5284. 86 Wilkins A, Spencer BW, Jain A, Gencturk B. A method for smoothing multiple yield functions. Int J Numer Methods Eng. 2019;1–16. https://doi.org/10.1002/nme.6215.

87 Karacan CO, Goodman G. Hydraulic conductivity changes and influencing factors in longwall overburden determined by slug tests in gob gas ventholes. Int J Rock Mech Min Sci. 2009;46:1162–1174. https://doi.org/10.1016/j.ijrmms.2009.02.005. 88 Wang S, Li X, Wang D. Void fraction distribution in overburden disturbed by longwall mining of coal. Environ. Earth Sci. 2016;75:151. https://doi.org/10.1007/ s12665-015-4958-6. 89 Li H, Chen Q, Shu Z, Li L, Zhang Y. On prevention and mechanism of bed separation water inrush for thick coal seams: a case study in China. Environ. Earth Sci. 2018;77: 759. https://doi.org/10.1007/s12665-018-7952-y.

17