Stress-Dependent Deformation and Permeability of a Fractured Coal Subject to Excavation-Related Loading Paths

Abstract

The deformation and permeability of coal are largely affected by the presence and distribution of natural fractures such as cleats and bedding planes with orthogonal and abutting characteristics, resulting in distinct hydromechanical responses to stress loading during coal mining processes. In this research, a two-dimensional (2D) fracture network is constructed based on a real coal cleat trace data collected from the Fukang mine area, China. Realistic multi-stage stress loading is designed to sequentially mimic an initial equilibrium phase and a mining-induced perturbation phase involving an increase of axial stress and a decrease of confining stress. The geomechanical and hydrological behaviour of the fractured coal under various stress loading conditions is modelled using a finite element model, which can simulate the deformation of coal matrix, the shearing and dilatancy of coal cleats, the variation of cleat aperture induced by combined effects of closure/opening, and shear and tensile-induced damage. The influence of different excavation stress paths and directions of mining is further investigated. The simulation results illustrate correlated variations among the shear-induced cleat dilation, damage in coal matrix, and equivalent permeability of the fractured coal. Model results are compared with results of previous work based on conventional approaches in which natural fracture networks are not explicitly represented. In particular, the numerical model reproduces the evolution of equivalent permeability under the competing influence of the effective stress perpendicular to cleats and shear-induced cleat dilation and associated damage. Model results also indicate that coal mining at low stress rates is conducive to the stability of surrounding coal seams, and that coal mining in parallel to cleat directions is desirable. The research findings of this paper have important implications for efficient and safe exploitation of coal and coalbed methane resources.

1 Introduction

Coal is a highly heterogeneous geological medium embedded with both endogenic and exogenous fractures. The distribution of these fractures on all scales exerts fundamental controls on coal mining-induced fracture behaviour and gas flow dynamics (Wang et al. 2013). The understanding of the role of interior fracture network within coal is crucial for efficient and safe recovery of coal and coalbed methane resources, as well as prevention of underground hazards such as coal and gas outbursts frequently occurring in low-permeability seams (Guo et al. 2012; Cao et al. 2020).

Orthogonal to the orientation of bedding planes, coal fracture patterns often contain two sets of mutually orthogonal endogenic fractures, named face and butt cleats, both being progressively formed by the uniform shrinkage of coal matrix due to loss of water and other volatiles during the coal maturation process (Dawson and Esterle 2010). Face and butt cleats are generally perpendicular to the direction of minimum and maximum principal stresses, respectively (Laubach et al. 1998; Dawson and Esterle 2010). While face cleats are through-going, butt cleats generally terminate at intersections with face cleats. Cleats are opening-mode fractures and have no cohesion along the frictional sliding direction as well as little shear offset (Weniger et al. 2016).

Both geomechanical and hydrological behaviour of coal are controlled by properties of the coal matrix and intrinsic coal cleat patterns (Wang et al. 2013). Extensive experiments have been carried out to characterise mechanical and hydraulic properties of coal under applied loading (Durucan and Edwards 1986; Wang et al. 2011, 2013; Xie et al. 2012; Zhang et al. 2016). The most prevailing experimental design employed is the triaxial gas flow permeability experiment. Typical loading paths applied in experiments include fixed pore pressure with increasing deviatoric stress (either by increasing axial stress or unloading confining stress), and fixed deviatoric stress with increasing pore pressure (Wang et al. 2011, 2013; Chen et al. 2013). To mimic the impact of different mining layouts in laboratory, loading paths and rates were used in triaxial experiments corresponding to the evolution of mining-induced stress acting upon in-situ coal masses, i.e., a simultaneous increase of axial stress and decrease of confining stress (Xie et al. 2012; Zhang et al. 2016; Ju et al. 2017). In the recent decade, the focus of experimental investigation has gradually shifted from the permeability evolution of intact coal to that of fractured coal (Izadi et al. 2011; Wang et al. 2011). It was found that fracture geometry changes permeability by several orders of magnitude (Wang et al. 2011). In particular, a through-going fracture oriented aligned with the flow direction would control the timing of flow and thus the amount of time-dependent sorption and swelling. The investigation has also been extended to the evolution of the interior coal fracture network and associated methane permeability variation under various excavation stress paths (Ju et al. 2017).

In continuum-based modelling of coal mining or coalbed methane drainage, the coal cleat system is not directly modelled due to unaffordable computational requirements; instead, model elements are assigned with effective permeability and porosity values representative of coal (Guo et al. 2009; Si et al. 2015). The underlying assumption is that fractured coal consisting of interconnected factures can be idealised as an equivalent porous continuum, with hydraulic properties being represented by an equivalent anisotropic hydraulic conductivity tensor in relation to mean fracture aperture and spacing. However, although anisotropic features in hydrological behaviour are preserved, the spatial correlation and hierarchical organisation of interior coal cleat network (Laubach et al. 1998; Weniger et al. 2016) cannot be well represented.

Numerical modelling of rock/coal behaviour on laboratory scale could reveal fundamental coupling physics of the fracture-matrix system by explicitly representing small-scale discontinuities and their responses. In previous research, emphasis has been placed on progressive failure of fractured rocks and associated variation in hydromechanical properties. Considering small-scale defects as local heterogeneities or ubiquitous microcracks, a series of continuum-based coupled hydromechanical models have been developed by incorporating degradation of mechanical properties, confining pressure-dependent dilatancy, and deformation-dependent permeability (Yuan and Harrison 2005; Zhu et al. 2006; Lu et al. 2013). Discrete fracture network (DFN) models have also been introduced to model hydrological behaviour under geomechanical effects using continuum/extended-continuum methods, discontinuum methods and hybrid finite-discrete element methods (Lei et al. 2017a).

Recently the rapid development of non-destructive characterisation methods such as stereomicroscopy and CT scanning and relevant image processing techniques allowed acquisition of high-quality images of interior cleat system within coal (Jing et al. 2017; Wang et al. 2018a). Statistical analysis was also made possible for extracting some important parameters and attributes of coal cleat patterns, such as cleat spacing, length, aperture and connectivity. These developments allowed the coal cleat system reconstructed from core samples to be used as conduits in fluid flow modelling through coal (Jing et al. 2017). However, to the authors’ best knowledge, the coal cleat system as a fracture network has not been explicitly represented in modelling hydromechanical behaviour of coal.

This work aims to explicitly model coal cleat behaviour and associated gas flow dynamics subject to excavation-related stress loading paths involving an increase of axial stress and a decrease of confining stress. A 2D numerical model is constructed incorporating a real coal cleat network. A coupled hydromechanical model that incorporates fracture deformation, damage evolution and stress-dependent permeability is used for the simulation. The evolution of cleat displacement, damage distribution, aperture variation, and flow velocity structure reorganisation are examined at representative stress states to demonstrate the significant role of coal cleat network. Influences of excavation stress paths and principal stress orientations are further investigated, with implications for efficient and safe exploitation of coal and coalbed methane resources discussed.

2 Governing Equations

The coupled hydromechanical model used for simulating coal cleat behaviour and gas flow dynamics subject to excavation-related stress paths considers solid deformation, fracture deformation, damage evolution and stress-dependent permeability as well as fluid flow in both porous matrix and natural fractures. This section presents mathematical formulations of each component of this model.

2.1 Solid Deformation

The mechanical equilibrium of the fractured porous rock is governed by:

$$\sigma_{ij,j} + f_{i} = 0,$$

(1)

where σ_ij is the stress and f_i is the body force. The stress–strain relation of intact rocks follows the linear poroelasticity law:

$$\sigma^{\prime}_{ij} = \sigma_{ij} - \alpha p\delta_{ij} = \frac{E\upsilon }{{\left( {1 + \upsilon } \right)\left( {1 - 2\upsilon } \right)}}\varepsilon_{kk} \delta_{ij} + \frac{E}{1 + \upsilon }\varepsilon_{ij} ,$$

(2)

where σ′_ij is the effective stress, ε_ij is the strain, E is the Young’s modulus, ν is the Poisson’s ratio, α is the Biot’s coefficient, p is the fluid pressure and δ_ij is the Kronecker delta. Compressive stress and strain are denoted to be positive, and tensile stress and strain to be negative. The Biot’s coefficient for coal is typically less than 1, varying between 0.76 and 0.95 (Durucan et al. 2009; Guo et al. 2014).

When the rock is coal as a reservoir to sorbing gases, the total strain contains the strain induced by effective stress and gas sorption-induced volumetric strain. The sorption-induced strain, governed by a Langmuir-type equation (Durucan et al. 2009), leads to matrix swelling proportional to the volume of gas sorbed. The sorption-induced strain may dominate the absolute permeability of coal at low gas pressures in coalbed reservoirs, where the geomechanical process could be simplified to be uniaxial strain and constant vertical stress conditions (Palmer and Mansoori 1996; Shi and Durucan 2005). The competing influence of effective stress and gas sorption is crucial to interpreting coal permeability for applications such as primary depletion, enhanced recovery or CO₂ sequestration operations. However, effective stress change induced by coal mining has an overriding effect on the deformation of coal matrix and variation in cleat aperture (Liu and Elsworth 1997; Adhikary and Guo 2015; Si et al. 2015); therefore, the sorption-induced strain is neglected in this model.

2.2 Fracture Deformation

The normal deformation of pre-existing fractures under compression is based on a hyperbolic relation (Fig. 1a) (Bandis et al. 1983):

$$v_{{\text{n}}} = \frac{{\sigma^{\prime}_{{\text{n}}} v_{{\text{m}}} }}{{k_{{{\text{n0}}}} v_{{\text{m}}} + \sigma^{\prime}_{{\text{n}}} }},$$

(3)

where v_n is the normal closure, σ′_n is the effective normal compressive stress, k_n0 is the initial normal stiffness, and v_m is the maximum allowable closure. Thus, the normal stiffness of the fracture also exhibits a nonlinear behaviour:

$$k_{{\text{n}}} = \frac{{\left( {\sigma^{\prime}_{{\text{n}}} + k_{{{\text{n0}}}} v_{{\text{m}}} } \right)^{2} }}{{k_{{{\text{n0}}}} v_{{\text{m}}}^{2} }}.$$

(4)

Fig. 1

Fracture deformation model: a the normal closure v_n as a function of the effective normal stress σ_n; b the shear stress τ as a function of the shear displacement u; c the dilational displacement v_s as a function of the shear displacement u

The shear stress of rock fractures is based on an elasto-plastic constitutive model with strain-softening (Fig. 1b) (Goodman 1976):

$$\tau = \left\{ {\begin{array}{*{20}l} {k_{{\text{s}}} u, \, } & {u \le u_{{\text{p}}} \, } \\ {\frac{{\tau_{{\text{p}}} - \tau_{{\text{r}}} }}{{u_{{\text{p}}} - u_{{\text{r}}} }}\left( {u - u_{{\text{r}}} } \right) + \tau_{{\text{r}}} ,} & {u_{{\text{p}}} < u \le u_{{\text{r}}} } \\ {\tau_{{\text{r}}} , \, } & {u \ge u_{{\text{r}}} \, } \\ \end{array} } \right.,$$

(5)

where τ is the shear stress, u is the shear displacement, k_s = ηk_n is the shear stiffness, η is the shear-to-normal stiffness ratio, u_p is the peak shear displacement, τ_p = k_su_p is the peak shear stress, u_r is the residual shear displacement, ϕ_f is the friction angle, and τ_r = σ’_ntanϕ_f is the residual shear stress. The dilational displacement v_s under a constant normal stress condition is related to the shear displacement in an incremental form as (Fig. 1c) (Saeb and Amadei 1992):

$${\text{d}}v_{{\text{s}}} = - \tan \phi_{{\text{d}}} {\text{d}}u,$$

(6)

where ϕ_d is the dilation angle given as:

$$\phi_{{\text{d}}} = \left\{ {\begin{array}{*{20}l} {\arctan \tau_{{\text{p}}} - \arctan \tau_{{\text{r}}} ,} & {{\text{for }}u \le u_{{\text{r}}} } \\ {0, \, } & {{\text{for }}u > u_{{\text{r}}} } \\ \end{array} } \right..$$

(7)

The fracture aperture is thus derived as (Cao et al. 2021, under review; Lei et al. 2016):

$$a = \left\{ {\begin{array}{*{20}l} {a_{0} + w{, }} & {\sigma^{\prime}_{{\text{n}}} < 0} \\ {a_{0} - v_{{\text{n}}} - v_{{\text{s}}} ,} & {\sigma^{\prime}_{{\text{n}}} \ge 0} \\ \end{array} } \right.,$$

(8)

where a₀ is the initial aperture, and w is the separation of opposing fracture walls if the fracture is under tension. The permeability of the fracture obeys the cubic law of the aperture a for laminar flow between parallel plates (Witherspoon et al. 1980):

$$k_{{\text{f}}} = \frac{{a^{2} }}{12}.$$

(9)

2.3 Damage Evolution

Rocks (e.g. coals) naturally contain pre-existing fractures and initial flaws (e.g., microcracks and cavities) (Kranz 1983; Anders et al. 2014), and damage in rocks usually initiates from the tips of pre-existing fractures/flaws. Damage evolution in rock material is modelled using an isotropic damage model (Fig. 2) based on the crack band approach (Bažant and Oh 1983). Damage in tension emerges in intact rocks if:

$$\sigma^{\prime}_{3} \le - \sigma_{{{\text{t0}}}} ,$$

(10)

where σ₃ is the minimum principal stress and σ_t0 is the tensile strength. The damage parameter D_t for tensile failure is then derived as (Jirásek and Bauer 2012):

$$D_{{\text{t}}} = \left\{ {\begin{array}{*{20}l} {0, \, } & {\varepsilon_{{\text{t}}} \ge \varepsilon_{{{\text{t}}0}} } \\ {1 - \frac{{\varepsilon_{{{\text{t}}0}} }}{{\varepsilon_{{\text{t}}} }}\exp \left( { - \frac{{\varepsilon_{{\text{t}}} - \varepsilon_{{{\text{t}}0}} }}{{\varepsilon_{{{\text{tf}}}} - \varepsilon_{{\text{t}}} }}} \right),} & {\varepsilon_{{\text{t}}} < \varepsilon_{{{\text{t}}0}} } \\ \end{array} } \right.,$$

(11)

where ε_t0 = − σ_t0/E₀ is the limit elastic tensile strain; E₀ is the initial Young’s modulus; ε_t is the equivalent tensile strain, defined as the maximum previously reached level of the tensile strain (Jirásek and Bauer 2012). The equivalent tensile strain in this work is derived as the norm of the positive part of the strain tensor (Mazars 1984):

$$\varepsilon_{{\text{t}}} = - \sqrt {\sum\limits_{i = 1}^{3} {\left( {\frac{{\left| {\varepsilon_{i} } \right| + \varepsilon_{i} }}{2}} \right)}^{2} } ,$$

(12)

where ε_i are the principal strains; the post-peak softening is controlled by the parameter ε_tf given as (Jirásek and Bauer 2012):

$$\varepsilon_{{{\text{tf}}}} = \frac{{G_{{\text{I}}} }}{{h_{{\text{b}}} E_{0} \varepsilon_{{{\text{t0}}}} }} + \frac{{\varepsilon_{{{\text{t0}}}} }}{2},$$

(13)

where G_I is the mode I energy release rate, and h_b is the width of crack band (at the scale of the size of finite elements where damage is localised; Jirásek and Bauer 2012).

Fig. 2

Damage model based on an elasto-plastic constitutive law with strain-softening: a stress–strain relation and b damage-strain relation

Damage in shear initiates when the Mohr–Coulomb failure criterion is satisfied (Tang et al. 2008):

$$\sigma^{\prime}_{1} - \frac{{1 + \sin \phi_{{\text{i}}} }}{{1 - \sin \phi_{{\text{i}}} }}\sigma^{\prime}_{3} \ge \sigma_{{{\text{c0}}}} ,$$

(14)

where σ′₁ is the effective maximum principal stress, ϕ_i is the internal friction angle, and σ_c0 is the uniaxial compressive strength. The damage parameter D_s for shear failure is then derived as:

$$D_{{\text{s}}} = \left\{ {\begin{array}{*{20}l} {0, \, } & {\varepsilon_{{\text{c}}} \le \varepsilon_{{{\text{c}}0}} } \\ {1 - \frac{{\varepsilon_{{{\text{c}}0}} }}{{\varepsilon_{{\text{c}}} }}\exp \left( { - \frac{{\varepsilon_{{\text{c}}} - \varepsilon_{{{\text{c}}0}} }}{{\varepsilon_{{{\text{cf}}}} - \varepsilon_{{\text{c}}} }}} \right),} & {\varepsilon_{{\text{c}}} > \varepsilon_{{{\text{c}}0}} } \\ \end{array} } \right.,$$

(15)

where ε_c0 is the limit elastic compressive strain given as:

$$\varepsilon_{{{\text{c0}}}} = \frac{1}{{E_{0} }}\left[ {f_{{{\text{c0}}}} + \frac{{1 + \sin \phi_{{\text{i}}} }}{{1 - \sin \phi_{{\text{i}}} }}\sigma^{\prime}_{3} - \tan \phi_{{\text{i}}} \left( {\sigma^{\prime}_{1} + \sigma^{\prime}_{2} } \right)} \right],$$

(16)

where ε_c = max(ε₁, 0) is the equivalent compressive strain; ε_cf controls the post-peak softening and is given as (Jirásek and Bauer 2012):

$$\varepsilon_{{{\text{cf}}}} = \frac{{G_{{{\text{II}}}} }}{{h_{{\text{b}}} E_{0} \varepsilon_{{{\text{c0}}}} }} + \frac{{\varepsilon_{{{\text{c0}}}} }}{2},$$

(17)

where G_II is the mode II energy release rate.

The damage parameter D takes the largest one of shear and tensile failure:

$$D = \max \left( {D_{{\text{t}}} ,D_{{\text{s}}} } \right).$$

(18)

During the damage evolution, the Young’s modulus of the rock matrix degrades as:

$$E = \left( {1 - D} \right)E_{0} .$$

(19)

The permeability of the rock matrix under stress and damage varies as (Zhu and Wei 2011):

$$k_{{\text{m}}} = k_{{0}} \left( {\frac{{\varphi_{{\text{m}}} }}{{\varphi_{{0}} }}} \right)^{\beta } \exp \left( {\gamma D} \right),$$

(20)

where k₀ is the initial matrix permeability, β = 3 is a power law exponent, γ = 5 is a damage-permeability coefficient, φ₀ is the initial porosity, and φ_m is the current porosity given as (Yuan and Harrison 2005):

$$\varphi_{{\text{m}}} = 1 - \left( {1 - \varphi_{{0}} } \right)\exp \left( { - \varepsilon_{{\text{v}}} } \right),$$

(21)

where ε_v is the volumetric strain.

2.4 Fluid Flow

The model considers steady-state fluid flow through both permeable matrix and multiple intersecting fractures in the rock mass (Lei et al. 2020). The continuity equation for fluid flow through the fractured porous media is given by:

$$\nabla \cdot q = 0,$$

(22)

where q is the flux. The equation of motion is given by Darcy’s law:

$$q = - \frac{k}{\mu }\nabla p,$$

(23)

where k is the permeability of rock matrix or fractures, and μ is the fluid viscosity.

3 Numerical Modelling

3.1 Coal Cleat Pattern and Properties

The aforementioned coupled hydromechanical model has been implemented in the finite element method (FEM)-based solver Comsol Multiphysics (Comsol 2019) to investigate the behaviour of fractured coals in response to excavation-induced stress variations. Before studying complex coal cleat systems, we have validated the model against analytical solutions (Renshaw and Schulson 2001) for the problem of wing crack growth emanating from the tips of a single pre-existing fracture in a solid under biaxial compressions (see Appendix A). The model has also been validated against analytical solutions (Snow 1969) for single-phase steady-state fluid flow through persistent fracture sets in porous media (see Appendix B). Close matches are obtained between the numerical and analytical results, justifying the accuracy and applicability of our model for analysing damage evolution and fluid flow behaviour in fractured media. Furthermore, a detailed validation of additional aspects of this hydromechanical coupling model can be found in Lei et al. (2021).

We then investigate the behaviour of fractured coals based on a laboratory-scale coal sample embedded with a clear cleat pattern (Fig. 3a). The coal cleat pattern was reconstructed from a 10 × 10 × 5 cm³ low-rank coal block collected in the Fukang mine area, China (Wang et al. 2018a). Surfaces of the coal block were hydraulically polished to eliminate all scratches and clearly expose the endogenic cleat geometries.

Fig. 3

Coal cleat pattern extracted from low-rank coal in the Fukang mine area, China (after Wang et al. 2018a, b): a polished coal surface, b binary image of the coal cleat pattern, and c digital coal cleat pattern imported into the numerical model

Surface topography images of the coal block were photographed by a stereomicroscope, with a resolution of 50 µm in length and 10 µm in width (Wang et al. 2018a). The images were then converted to binary images to distinguish cleat traces and coal matrix (Fig. 3b). The proportions of apertures between 10 µm and 50 µm are between 49.39 and 86.62% for coal samples collected from the same mine area.

The coal cleat system within a 5 × 5 cm² area at the bottom of the coal sample is extracted, digitalised and imported into the 2D numerical model (Fig. 3c). Individual coal cleats are represented by 1-D interface elements. The initial aperture of each cleat set (face or butt) is assumed to be constant, while a face to butt cleat aperture ratio of 1.3 is set to account for the anisotropy in aperture distribution (Wang et al. 2018b). The model domain containing distributed cleats is discretised using an unstructured mesh with an average element size of 1.5 mm.

To ensure the numerical model constructed to realistically represent the macroscopic mechanical behaviour of fractured coal, numerical calibration of mechanical properties of both coal matrix and cleats is performed before modelling the progressive coal mass failure and associated methane flow under excavation stress paths. Mechanical properties used for the calibration were characterised through a series of uniaxial compression tests (4 samples) and Brazilian tests (5 samples) for coal collected from the same area (Han et al. 2014). Physical and mechanical properties of the fractured coal are presented in Table 1.

Table 1 Physical and mechanical properties of the fractured coal (after Han et al. 2014)

The calibration is performed by matching against laboratory experiments the macroscopic response of deformation and failure behaviour of the fractured coal through numerical uniaxial compression tests and Brazilian tests (Potyondy and Cundall 2004; Li et al. 2014). Herein, the deformation of fractured coal is characterised by elastic properties such as Young’s modulus and Poisson’s ratio, while the failure behaviour is dominated by strength parameters such as uniaxial compressive strength and tensile strength. The calibration is based on a series of coal samples randomly cropped from Fig. 3c.

Figure 4 presents examples of calibration results with final damage patterns in uniaxial compression tests and Brazilian tests. The stress–strain curve of the uniaxial compression test reproduces four typical stages observed in laboratory experiments, namely initial fracture closure stage, elastic deformation stage, plastic deformation damage and post-failure stage. Of note is that the explicit modelling of closure and slippage of embedded fractures allows to successfully capture the fracture closure stage, which has rarely been reported in previous literature. This indicates that the macroscopic behaviour of fractured coal is influenced by the interplay of the coal matrix and cleat system. In addition, the constant increment of damaged area starting at the plastic deformation stage is in good consistency with the increasing number of acoustic emission (AE) events representing localised rock failure at the same stage in laboratory experiments. The failure of coal samples in either Fig. 4a or b is accompanied by a surge in the damaged area, corresponding to the macroscopic shear or tensile failure in laboratory experiments. In addition, damage in coal matrix generally originates from cleat tips in both cases. While the final damaged area runs through the coal sample in the numerical uniaxial compression test, only limited damaged area appears along the middle line of the Brazilian disk.

Fig. 4

Calibration results of Fukang coal: the stress–strain curve and the final failure pattern in a uniaxial compression test, and b Brazilian test (the red bars represent increments of proportion of damaged area to total domain area). Insets show the distribution of damage in the samples

Physical and mechanical properties of coal matrix and cleats used in following modelling work are presented in Table 2. While physical parameters (density, porosity, etc.) of coal are directly used, deformation (Young’s modulus, Poisson’s ratio, etc.) and strength parameters (friction angle, cohesion, etc.) of both coal matrix and cleats are obtained through calibration. In this model, calibrated values of G_II for coal matrix is smaller than G_I to better match the experimental measurements.

Table 2 Physical and mechanical properties of the fractured coal used in the model

3.2 Numerical Model Setup

The 5 cm × 5 cm numerical model incorporating the extracted coal cleat pattern (Fig. 3c) is used to simulate the geomechanical and hydrological response of fractured coals under excavation stress paths. The 2D model represents the horizontally placed fully-saturated fractured coal ahead of a coal mining face under an isotropic horizontal effective stress of 15 MPa (where the initial horizontal stress is 21 MPa, the pore pressure is 7.5 MPa, and the Biot’s coefficient for coal is 0.8). It is worth mentioning that unsaturated conditions may occur when the mining face reaches the target coal cell, which is not modelled in the current study. The advance of coal mining face leads to an increased deviatoric effective stress, and the prevailing stress field can be characterised by the strike-slip stress regime.

Three commonly practiced coal mining methods, namely protective seam mining, sub-level caving mining and mining without coal pillars, are considered in designing excavation stress paths in our numerical modelling. Figure 5a illustrates the typical stress path sequentially experienced by the model domain as the mining face approaches for the sub-level coal caving mining scenario. The coal mining perturbation marks a transition of the maximum principal stress σ₁ from initial stress to abutment stress and then to peak stress, followed by the post-failure residual stress. Note that the maximum peak stress 2.5–3.0σ₀ at a low confining stress 0.2–0.4σ₀ can only be achieved for strong coal, while weak coal is usually fractured and loses bearing capacity before reaching this stress.

Fig. 5

Typical stress paths and stress states used in the numerical modelling: a a typical stress path undergone by coal in sub-level coal caving mining (after Xie et al. 2012); b three loading paths representing different coal mining layouts (after Ju et al. 2017), and c schematic description of the far-field stress applied at different angles to the model domain

Stress paths applied in the numerical modelling sequentially mimic an initial equilibrium phase (OA in Fig. 5b) and a mining-induced perturbation phase (AB, AC or AD in Fig. 5b) involving an increase of axial stress σ₁ and a decrease of confining stress σ₃ until the coal fails. Excavation stress paths for these mining methods are characterised by the stress rate, defined as the ratio of axial loading to lateral unloading rates (Ju et al. 2017). Stress rates used for excavation stress paths of the three mining methods are 1, 2 and 3, respectively.

During the simulation, principal stresses following the three designated excavation-related loading paths are applied orthogonally to the model domain until failure, representing the advance of coal mining parallel to the direction of face cleats. While stresses are imposed to the top and lateral model boundaries, the bottom of the model is constrained by a roller boundary. To investigate the influence of the direction of mining, a set of modelling scenarios are further examined by applying the far-field stress at a range of different angles θ with an interval of 30° (Fig. 5c).

The equivalent permeability of the fractured coal over the simulated coal mining process is calculated simultaneously based on single-phase steady-state fluid flow simulations, where a prescribed macroscopic pressure differential is imposed at two opposite boundaries while the other boundaries are impervious. Assuming that gas migration through coal is an isothermal process and follows the ideal gas law, the gas permeability can be calculated by (Wang et al. 2013):

$$k = \frac{\mu L}{A}\frac{{2P_{0} Q}}{{P_{2}^{2} - P_{1}^{2} }},$$

(24)

where L and A are the length and cross-sectional area of the sample, μ is the viscosity of methane, Q is the flow rate, P₀ is the atmospheric pressure, and P₁ and P₂ denote downstream and upstream pressures, respectively.

Laboratory gas permeability measurement for coal can be influenced by gas slippage along the pore or fracture walls, known as the Klinkenberg effect (Klinkenberg 1941). When the pore radius approaches to the mean free path of gas molecules, the appreciable collision between gas molecules and the pore wall results in an additional flux and an apparent increase in permeability at low gas pressures. The Klinkenberg effect becomes negligible under very high gas pressures and for permeabilities above about 20 mD (Somerton et al. 1975). This gas slippage effect is not modelled in the gas permeability measurement in this numerical modelling.

The gas flowing through the deformed coal mass with stress-dependent apertures is methane, whose viscosity and compressibility are 1.10 × 10^–5 Pa s and 10^–5 Pa⁻¹, respectively. The permeability of coal matrix is assumed to be 1 × 10^–17 m², which gives a high contrast to cleat permeability so that methane flow through the domain is dominated by the cleat system (Seidle 2011; Weniger et al. 2016). A low pressure differential (P₂− P₁) is used for permeability measurement to minimise the influence of effective stress due to gas flow on the fracture behaviour and damage evolution in coal.

4 Model results and Analysis

4.1 Geomechanical and Hydrological Response at Representative Stress States

Figure 6 shows the geomechanical and hydrological behaviour of the fractured coal in response to sub-level caving coal mining (stress rate = 2). Representative stress states examined for the model include the unstressed state (σ₁ = σ₃ = 0), initial stress state (σ₁ = σ₃ = σ₀), abutment stress (σ₁ = 1.5σ₀, σ₃ = 0.75σ₀), and peak stress (σ₁ ≈ 2σ₀, σ₃ ≈ 0.5σ₀). The post-failure behaviour of the model largely depends on the extent of coal failure and the specific residual stress, and is thus not investigated here.

Fig. 6

Geomechanical and hydrological response of the fractured coal to sub-level caving coal mining (stress rate = 2): a damage distribution, b fracture displacement, c fracture aperture, d normalised velocity field for gas flow in the x direction, and e normalised velocity field for gas flow in the y direction (velocity values are normalised with the mean velocity under the unstressed state)

When the model domain is loaded by the initial stress, damage of coal matrix initiates between neighbouring cleat tips, where high stress concentrates (Fig. 6a). These local damaged areas gradually expand along cleat traces with the increase of axial stress. When reaching the peak stress, these damage areas coalesce into a macroscopic shear failure zone across the coal cleat system.

Figure 6b illustrates the variation of shear displacement (or shear offset) of cleats under the excavation stress path with stress rate = 2. While the cleat displacement field does not exhibit apparently heterogeneous distribution at initial stress, localised large offsets manifest at abutment stress along long cleats (including some through-going face cleats), which have little resistance against shear slippage. The nonuniform spatial distribution of cleat displacement further exacerbates at peak stress. Combining with Fig. 6a, it is worthwhile to note that damaged areas closely follow the spatial distribution of coal cleats experiencing relatively large slippage, indicating that damage in coal matrix is strongly affected by the geomechanical behaviour of the coal cleat system.

Figure 6c presents the evolution of cleat apertures under the prescribed excavation stress path. Both face and butt cleats are largely compressed at initial and abutment stresses. When reaching peak stress, the low confining stress as well as the macroscopic failure across the model domain leads to opening or shear dilation of a large proportion of cleats, in particular butt cleats. To quantitatively analyse the cleat aperture distribution, the probability density functions of cleat apertures at representative stress states are presented in Fig. 7. Without any stress loading, face and butt cleats have associated with constant apertures of 0.050 and 0.038 mm, respectively. The cleat apertures obey a quasi-Gaussian distribution with an average of 0.025 mm under the initial or abutment stress state. The distribution of cleat aperture under the abutment stress state has a slightly larger spread than that at initial stress, because higher axial stress leads to closure of face cleats (i.e., more small apertures), while lower confining stress causes opening of butt cleats (i.e., more large apertures). When reaching the peak stress state, the aperture of coal cleats experiencing opening or shear dilation increases tremendously up to several times of initially prescribed cleat apertures.

Fig. 7

Probability density functions of fracture aperture under different stress stages during sub-level caving coal mining (stress rate = 2)

As illustrated in Fig. 6d, e, coal cleats aligned with the gas flow direction form primary pathways through the domain. The gas flow velocity field is not only related to cleat aperture distribution, but also influenced by cleat orientation, density and connectivity. In both gas flow cases, a series of sub-parallel pathways, constituted by coal cleats directly connected to pervious boundaries, allow relatively large gas flow velocities. Streamlines in the x direction are more continuous as gas flow is along through-going face cleats, while those in the y direction are more tortuous along pathways involving non-persistent butt cleats and segments of face cleats. There are more pathways in the y direction than in the x direction due to a relatively small spacing for butt cleats. The maximum gas flow velocity could reach over 3 times of the mean velocity under the unstressed state. It is interesting to note that when a face cleat branches under the unstressed state, the gas flow velocity is reduced accordingly due to the mass conservation constraint. When reaching the initial or abutment stress state, the closure of coal cleats subjected to stress loading significantly suppresses the gas flow through the model domain. Gas flow breaks through the model domain more easily under the peak stress state as opposed to that under the initial or abutment stress state, and large gas flow velocities occur along pathways constituted by connected segments of coal cleats having large apertures.

Figure 8 presents the normalised velocity distribution for gas flow through the model domain in both x and y directions at different representative stress states. Although gas flow velocities at most coal cleats concentrate around the mean velocity under unstressed state, normalised gas flow velocities in both directions span over several orders of magnitude. It is interesting to correlate normalised velocities above 1 with dominant fluid pathways through the model domain. The coal cleat aperture field at unstressed state and initial stress has the largest and smallest proportion of cleats having normalised gas flow velocities above 1, respectively. As the axial stress increases from initial to abutment stress and finally to peak stress, the proportion of normalised gas flow velocities above 1 first gently and then remarkably increases, due to continuous increase in aperture by shear offset as well as opening of coal cleats. The maximum gas flow velocity in either flow direction is achieved at peak stress with macroscopic failure, which could be 50 times larger than the corresponding mean velocity under the unstressed state.

Fig. 8

Probability density functions of flow velocity at different stress stages during sub-level caving coal mining (stress rate = 2): a gas flow in the x direction, and b gas flow in the y direction (velocity values are normalised with the mean velocity under the unstressed state)

4.2 Effect of Excavation Stress Path

Figure 9 presents geomechanical and hydrological response of the fractured coal when the confining stress is relieved to 0.75 times of the initial stress (σ₃ = 0.75σ₀) for coal mining using different mining methods (i.e., characterised by different stress rates). Coal mining characterised by a larger stress rate causes more nonuniform cleat displacement distribution and associated damaged area, as shown in Fig. 9a, b. In contrast, the change in the face cleat aperture field is not sensitive to the axial stress, as governed by the hyperbolic relation between the cleat closure and effective normal stress (Eq. (3) and Fig. 1a). Nevertheless, a few butt cleats open (shown in red colour) for coal mining characterised by a stress rate of 3 (see Fig. 9c).

Fig. 9

Geomechanical and hydrological response of the fractured coal under the abutment stress state (σ₃ = 0.75σ₀) for coal mining at different stress rates: a damage distribution, b fracture displacement, c change in fracture aperture, d normalised velocity field for gas flow in the x direction, and e normalised velocity field for gas flow in the y direction (velocity values are normalised with the mean velocity under the initial stress state)

Largely dependent on the cleat aperture distribution, gas flow velocities through the domain in both x and y directions are not very sensitive to the stress rate, as presented in Fig. 9d, e. To quantify the gas flow velocity distribution in both flow directions, the probability density functions of velocities normalised with the mean velocity under the initial stress state is plotted in Fig. 10 for coal mining at different stress rates. Coal mining at a lower stress rate exerts less axial stress normal to face cleats, yielding higher normalised gas flow velocities in the x direction (Fig. 10a). In contrast, the difference in the stress rate (axial stress) has a minimal influence on butt cleat apertures and thus on gas flow velocities in the y direction (Fig. 10b).

Fig. 10

Probability density functions of normalised flow velocity through the model domain for coal mining at different stress rates: a gas flow in the x direction, and b gas flow in the y direction (velocity values are normalised with the mean velocity under the initial stress state)

Figure 11 presents the damage and equivalent permeability evolution under different stress paths. The fractured coal approached by coal mining at a larger stress rate fails at a higher confining stress, corresponding to a larger deviatoric stress. The variation of equivalent permeability of the system is dependent on the interior cleat aperture field, which is in turn under combined yet counteractive influences of effective normal stress and shear displacement (see Eq. (8)). As the deviatoric stress increases, the equivalent permeability in the x direction first decreases when the increasing effective normal stress plays a dominant role, and then begins to increase as shear dilation exerts more influence. The fluctuation of the equivalent permeability at the late stage is caused by the fact that the effective normal stress increases continuously driving a gradual permeability decline, while shear displacements of newly slipped cleats as well as formation of new damage in coal matrix occur consecutively leading to permeability increments. In comparison, equivalent permeability in the y direction first increases constantly as the confining stress is unloaded, and then at an accelerated rate as the damage provides more conductive paths. It is also noted that the initial equivalent permeability of the fractured coal in the x direction is 2.11 times of that in the y direction, which is attributed to the intrinsic anisotropy in connectivity and aperture of the coal cleat system.

Fig. 11

Global damage proportion and equivalent permeability associated with deviatoric stress for coal mining at different stress rates: a gas flow in the x direction, and b gas flow in the y direction

4.3 Effect of far-Field Stress Orientation

The model domain with cleats oriented obliquely to the direction of coal mining is susceptible to damage and fails before the axial stress reaches 1.5σ₀, therefore, geomechanical and hydrological response of the fractured coal under a lower abutment stress (σ₁ = 1.25σ₀, σ₃ = 0.875σ₀) is presented in Fig. 12 for different principal stress orientations (stress rate = 2).

Fig. 12

Geomechanical and hydrological response of the fractured coal at abutment stress (σ₁ = 1.25σ₀, σ₃ = 0.875σ₀) for different principal stress orientations: a damage distribution, b fracture displacement, c change in fracture aperture, d normalised velocity field for gas flow in the x direction, and e normalised velocity field for gas flow in the y direction (velocity values are normalised with the mean velocity under the initial stress state)

As shown in Fig. 12a, damage inflicted at orthogonal settings (θ = 0 or 90°) concentrates at the tips of neighbouring cleats having high stress concentration. In contrast, damage also concentrates at the left middle part and right bottom corner in addition to these vulnerable areas for θ = 30° and 60° model scenarios, and at the left bottom corner for θ = 120° and 150° model scenarios. The spatial distribution of shear displacement is closely associated with damage and resembles the damage distribution for the respective principal stress orientation (Fig. 12b). As the model domain is biaxially stressed over the excavation process, aperture enlargement of some cleats (Fig. 12c) is primarily related to shear dilation induced by large shear displacement (Fig. 12b).

As coal cleats oriented obliquely to the direction of axial loading are more susceptible to shear slippage, the Mohr circle diagram (Fig. 13) is applied to analyse the geomechanical response of the fractured coal for different principal stress orientations. When the principal stress is orientated parallel or orthogonal to the model domain, neither face nor butt cleats are aligned with the direction of shear slippage so they stay quite suppressed (Fig. 12b). In comparison, both butt and face cleats are susceptible to slippage along their respective directions in θ = 30° and 150° model scenarios, and stress states for butt cleats are closer to failure envelopes of fractures. The opposite is true for the slippage potential of face and butt cleats in θ = 60° and 120° model scenarios. Of note is that face and butt cleats are not perfectly oriented horizontally or vertically, so the geomechanical response of the model domain is not the same for θ = 30° and 150° model scenarios, and for θ = 60° and 120° model scenarios. For example, the slightly left inclined butt cleat at the top middle part of the model domain is activated in θ = 120° and 150° model scenarios, but remains suppressed in θ = 30° and 60° model scenarios (Fig. 12b).

Fig. 13

Mohr circle analysis of the cleat system subjected to far-field stress loading at different orientations

It could be seen from Fig. 12d that the θ = 90° scenario has larger normalised gas flow velocities in the x direction than the θ = 0° scenario, as face cleats aligned with the x direction are less stressed. Likewise, methane flows more easily in the y direction through butt cleats in the θ = 0° scenario than in the θ = 90° scenario (Fig. 12e). In oblique principal stress orientation scenarios, through-going primary fluid pathways are also formed with large gas flow velocities, as influenced by the shear dilation induced aperture increase (Fig. 12b, c). To further elucidate this phenomenon, normalised velocity distributions for gas flow through the domain in both x and y directions are plotted in Fig. 14 for different principal stress orientations. There are higher gas flow velocities in the x direction in oblique principal stress orientation scenarios than in the θ = 90° scenario (Fig. 14a), and higher gas flow velocities in the y direction in oblique principal stress orientation scenarios than in the θ = 0° scenario (Fig. 14b).

Fig. 14

Probability density functions of flow velocity for different far-field stress orientations: a gas flow in the x direction, and b gas flow in the y direction (velocity values are normalised with the mean velocity under the initial stress state)

Figure 15 presents the damage and equivalent permeability evolution for different far-field stress orientations. The fractured coal fails at relatively high deviatoric stresses (> 20 MPa) for orthogonal stress orientation scenarios (θ = 0° or 90°), and at relatively low deviatoric stresses (< 10 MPa) for oblique stress orientation scenarios (θ = 30° or 60°). As shown in Fig. 15a, the equivalent permeability in the x direction first decreases and then increases with the increased deviatoric stress in all scenarios except for the θ = 90° case, where face cleats are unloaded to allow a continuous permeability increase in the x direction. The model in the θ = 0° scenario has the lowest x direction equivalent permeability before failure due to the largest effective stress applied normal to face cleats. Likewise, the y direction equivalent permeability experiences an initial decrease followed by an increase in all scenarios except for the θ = 0° case (Fig. 15b). The lowest y direction equivalent permeability is achieved in the θ = 90° scenario, when butt cleats aligned with the gas flow direction are most significantly compressed.

Fig. 15

Global damage proportion and equivalent permeability associated with deviatoric stress for different cleat orientations with respect to the direction of coal mining: a gas flow in the x direction, and b gas flow in the y direction

5 Discussion

5.1 Controlling Factors of Coal Permeability

Permeability of coal in situ is usually dominated by bedding planes, with the highest permeability in general being parallel to bedding planes in the face cleat direction, and the lowest perpendicular to bedding planes (Li et al. 2014; Seidle 2011). During the coal mining process that is parallel to coal bedding planes, it is expected that significant unloading takes place in the horizontal direction for coals ahead of the mining face, which may cause enlargements of cleat apertures, but less variations in overburden stresses normal to bedding planes. Therefore, although the cleat system may not exert key controls on coal permeability in natural states, it is of great importance to investigate the stress-dependency of coal cleat permeability under excavation unloading conditions. Our model results have shown that the horizontal coal permeability perpendicular to the unloading direction could increase by one order of magnitude during the mining process (Figs. 11, 15). This suggests that our current 2D model focussing on cleat patterns may provide useful insights into the permeability variation of coals subject to mining operations. We would like to point out that it is essential to also consider the influence of bedding planes on coal permeability when solving real engineering problems, for which a 3D numerical model (like the ones developed in Lei et al. 2015,2017b) may be needed to realistically represent the hydromechanical behaviour of 3D layered coal systems.

The numerical model developed in this work provides a platform to analyse permeability evolution dominated by coal cleat patterns against effective stress and associated damage. The equivalent permeability of the fractured coal before failure measured in the numerical model is in the order of 10^–14 m² along either face or butt cleat orientation. The numerically calculated permeability is in general comparable to experimentally measured values for coal cores collected from the same mine area, which is at 10^–15-10^–14 m² level (Tao et al. 2017). The in-situ aperture of coal cleats governing the equivalent permeability in the model ranges from 9 to 54 µm with a median of 27 µm (Fig. 7), which fall in the typical estimated range varying between 0.1 and 100 µm (Close 1993; Weniger et al. 2016).

The parallel bedding direction fracture profile of the fractured coal modelled exhibits mesh-like and anisotropic features (Fig. 3). Face cleats are more persistent, while butt cleats are more tortuous and highly branched conduits. The anisotropic character in horizontal coal permeability could be quantified by the permeability anisotropy ratio (PAR), defined as the ratio of permeability in the face cleat direction over that in the butt cleat direction (Massarotto et al. 2003). The PAR is dependent on the face to butt cleat aperture ratio, the face to butt cleat ratio (or face to butt cleat spacing ratio), and tortuosity of percolation pathways. The PAR for the coal permeability under the initially isotropic stress in the model is 2.11, very close to the measured PAR (2.34) for coal samples collected from the same mine area (Wang et al. 2018b). This PAR also falls into ranges previously reported for both laboratory-measured and field-derived PARs, i.e., 1.35–19 (Massarotto et al. 2003) and 1.8–17 (Koenig and Stubbs 1986), respectively. For this specific coal cleat system, the PAR measured in the model seems to follow closely the cubic law of the face to butt cleat aperture ratio (1.3³ = 2.20). This similarity between PAR and face to butt cleat aperture ratio indicates that primary in-situ conduits for gas flow through coal are cleat apertures in the tens of micrometre range.

This work demonstrates the significant role of interior cleat patterns on damage and permeability evolution of coal by explicitly modelling the coal cleat behaviour under whole excavation stress paths. Both favourable directions of slippage in geomechanical behaviour and dominant gas flow directions in hydrological behaviour lie in coal cleat patterns characterised by the distinctive mesh-like feature. The anisotropic characteristics of coal is more influenced by face cleats, which have continuous run length and are sometimes through-going. Wang et al. (2011) investigated the role of fracture geometry in coal permeability by conducting permeability experiments for three different samples containing either multiple small fractures, a longitudinal through-going fracture or a radial through-going fracture, and found that only the sample with a longitudinal through-going fracture did not show sorption-induced swelling behaviour. This is believed to be because the sorption-induced swelling strain cannot be converted into fracture closing strains without the presence of connected bridges between opposite walls of through-going fractures (Izadi et al. 2011).

5.2 Implications for Excavation-Related Gas Flow Dynamics

The explicit representation of coal cleat pattern and consideration of fracture hydromechanical behaviour in our physics-based numerical model allow us to investigate the mechanisms that govern coal permeability dynamics subject to excavation-related stress paths. Previous numerical modelling of coal permeability (e.g., Guo et al. 2009; Si et al. 2015) was usually based on some experimentally-derived stress-permeability relationships, e.g. expressed as the cubic law of porosity or fracture aperture influenced by the normal effective stress. Such a relationship presumes that coals remain intact without newly formed porosity or fractures, which is in general the case in coalbed methane drainage practices. However, coals subject to excavation-induced stress changes are susceptible to failure, and the resulting shear dilation of fractures and evolution of damages would have a large impact on the coal permeability. Thus, the empirical coal stress-permeability model does not consider the full physics involved in the excavation process, and cannot accurately reflect the dynamic permeability variation subject to excavation-related stress paths, in particular the permeability enhancement due to shear dilation of coal cleats, and formation of new fractures. As opposed to presuming any stress-permeability relationship, we tackle this coupled hydromechanical problem directly solving the governing equations of fundamental physical processes, involving normal opening/closure and shear dilation of fractures, deformation and shear/tensile damage of coal matrix, and fluid flow in fractured porous media. Therefore, this physics-based model provides a way to investigate the underlying mechanistic processes that govern the variation of coal permeability, faithfully honouring the combined effects of coal cleat pattern, stress path, and fracture/damage behaviour.

The fractured coal simulated undergoes typical coal mining-induced stresses states until failure. The numerical model realistically reproduces the dynamic coal permeability evolution under the competing influence of the effective stress perpendicular to cleats, shear dilation of cleats, and associated damage. The equivalent permeability of the fractured coal in parallel face advance direction under excavation stress paths first decreases as pre-existing cleats close at low deviatoric stresses, followed by an increase when the influence of shear dilation of pre-existing cleats and associated damage takes over at high deviatoric stresses (Figs. 11, 15). Model results in terms of damage and permeability are in good consistency with the schematic process-based model for coal mining-induced stress, pore pressure and permeability changes (Wang et al. 2013), simulated stress and permeability profile at the coal mining horizon (Durucan 1981), and triaxial gas flow experiments measured permeability evolution against the deviatoric stress (Wang et al. 2013; Ju et al. 2017). Equivalently, the current model represents the behaviour of coal cells at a specific time but at different locations with respect to the mining face, characterised by variable stress conditions (initial stress, abutment stress and peak stress states).

The modelling work suggests that the protective seam mining method (stress rate = 1) is the safest amongst the three mining methods evaluated by achieving the least damage and highest coal permeability under the same confining stress condition. Coal mining parallel to directions of cleats, as practiced almost since the start of coal mining (McCulloch et al. 1974), is confirmed to be preferable owing to the minimum fracture slippage and damage to surrounding coals as well as maximum possible directional permeability towards mine workings. In addition, as coal tends to break along cleats, it was easier to mine along cleats than at an angle to them before the application of continuous miners (McCulloch 1975).

5.3 Representativeness of the Coal Cleat Pattern

The coal cleat system considered in this work is in general representative of typical coal cleat patterns. Geometrically, the coal cleat pattern is comprised of the orthogonal sets of face and butt cleats. Both face and butt cleat spacings are within the typical range from 0.1 to 5 cm (Seidle 2011). The range of apertures (mostly between 10 and 50 µm) is also comparable to that for other coals, such as those from San Juan Basin (ranging from 0.01 to 0.30 mm) (Close 1993). Values for other fracture attributes (fracture porosity, surface density and connectivity rate, etc.) of this cleat pattern, as measured by Wang et al. (2018a), are also similar to those from different coal mines. Nevertheless, features of coal cleat systems are closely related to the geological setting and history of the coal-bearing strata (Laubach et al. 1998). Specifically, coal cleats are spaced more sparsely when formed under higher tectonic stresses. The intensity and direction of tectonic movement determine the size and extension direction of coal cleats, and the stress anisotropy causes the difference in cleat apertures. This coal sample was taken from a coal seam located on the thrust nappe in the southern part of the Junggar Basin (Wang et al. 2018a, b). The coal-bearing strata has undergone strike-slip movement and intense tectonic deformations in the geological history, leading to extrusion nappe characteristics. Given the geological setting in the mining area, coal is subject to enormous horizontal compressive stresses. Therefore, the coal cleat pattern considered is more representative of coal cleat systems in tectonically disturbed and horizontal stress-dominated regions.

In the present work, the numerical model incorporating a 2D coal cleat system may not fully represent the complex 3D coal cleat structure in coalbeds. To inform numerical investigation with more precise interior structure of coal samples, great efforts have been devoted to extract fine fracture network within coal using X-ray CT tomography, and to quantify descriptive parameters of cleats such as cleat length, size, aperture, intensity, density, spacing and shape parameter aided by image processing methods (Weniger et al. 2016; Busse et al. 2017). Better image quality and resolution can be achieved by other promising measuring techniques, including the µCT scanning technique capable of capturing new fractures (Jing et al. 2016; Ju et al. 2017), and the scanning electron microscope (SEM) technique to obtain precise aperture sizes smaller than the image resolution (Jing et al. 2016). The use of these techniques greatly facilitates the justification of numerical models and comparison of model results against experimental observations. For example, observations using the µCT technique indicated that new fractures generated prior to failure under external loads have poor connectivity throughout coal samples barely contributing to coal permeability (Ju et al. 2017); therefore, exogeneous fractures (or mining-induced fractures) are justified to be mimicked by damaged areas around coal cleats.

It is worth mentioning that a constant initial aperture is assumed for either face or butt cleats in the present numerical model. However, aperture of fractures in geological media is found to follow a power law relationship with fracture length (Bonnet et al. 2001; Olson 2003). The cleat aperture as a function of cleat length in coal may exhibit two types of relationships: one with cleat aperture and length being positively correlated, and the other with very short cleats having larger apertures than average (Busse et al. 2017). In addition, spatial variability is also present in the initial aperture distribution. Fracture aperture is often characterised by a power law distribution (Barton and Zoback 1992; Bonnet et al. 2001). Coal cleat aperture was also found to follow a power law frequency distribution for both face and butt cleats (Busse et al. 2017). The aperture variability effect has recently been addressed by introducing a lognormal distribution for the variable initial aperture field in addition to stress effect, with applications to studying the anomalous transport (Kang et al. 2019) in subsurface fractured rocks. Consideration of such an aperture variability effect may also be included in the future coal cleats-related research.

Based upon outcomes from this cleat-scale modelling work, future work includes the use of the coupled hydromechanical model developed to carry out engineering-scale modelling of excavation-related gas flow behaviour using large-scale coal fracture networks. The capability of our model to study field-scale problems has already been demonstrated in our recent work of fluid injection into fractured crystalline rocks (Lei et al. 2021). Coal fractures patterns may be derived from surface outcrops of coal deposits reconstructed using laser scanning, or underground coal fracture systems captured by borehole logging/geophysical techniques such as 3D seismic, ground penetrating radar and electrical exploration (Pan et al. 2019, 2020). A number of typical coal fracture patterns will be considered in modelling the gas flow behaviour under excavation scenarios prevalent in underground mining practice.

6 Conclusions

In this work, a coupled hydromechanical model that explicitly represents coal cleat pattern and considers cleat deformation and matrix damage is developed and applied to simulate the gas flow dynamics of fractured coal under realistic excavation-related stress paths during mining operations. The influence of different excavation stress rates and directions of mining on damage and permeability evolution of coal embedded with a cleat network is investigated. This work demonstrates the significant role of interior cleat patterns by presenting correlated variations of the shear-induced cleat dilation, damage in coal matrix, and equivalent permeability of the fractured coal.

Model results are compared favourably with results of previous work based on conventional approaches in which natural fracture networks are not explicitly represented. In particular, the numerical model reproduces the evolution of equivalent permeability under the competing influence of the effective stress perpendicular to cleats and shear-induced cleat dilation and associated damage. Our results indicate that coal mining at low stress rates is conducive to the stability of surrounding coal seams, and that coal mining in parallel to cleat directions is desirable. The results have important implications for designing efficient and safe coal mining operations (in terms of mining method and layout, direction of mining, etc.) and guiding in-seam gas pre-drainage ahead of mine development to mitigate gas-related hazards in underground coal mines. The coupled hydromechanical model developed is promising to be extended to applications with more complex fracture patterns or those on a larger scale in coal beds.

Appendices

Appendix A: Model Validation for Damage Evolution Around a Pre-existing Fracture

Stable wing cracks are deemed to form at tips of pre-existing primary fractures under biaxial compression. These wing cracks tend to extend along the maximum principal stress under little to no confinement, while wing crack growth is restricted with localised damage along a narrow band under low to moderate confinement. The length ratio of wing cracks and the primary fracture asymptotically approaches the limiting value (Renshaw and Schulson 2001):

$$\frac{{l_{{\text{w}}} }}{{l_{{\text{p}}} }} = \frac{{\sigma_{1} /\sigma_{3} - 1 - \mu (\sigma_{1} /\sigma_{3} + 1)}}{4.3},$$

(25)

where l_w is the length of wing cracks, 2l_p is the length of the primary fracture, σ₁ and σ₃ are the maximum and minimum principal stresses, respectively, and µ is the friction coefficient along the sliding segment of the primary fracture.

A numerical model containing a single pre-existing fracture is constructed to validate the model for simulating damage evolution around the fracture under different confinement ratios. No friction is assumed between opposing walls of the primary fracture. Model results reproduce the wing crack propagation along the maximum principal stress, as shown in Fig. 

Fig. 16

Validation of the numerical model against the analytical solution (Renshaw and Schulson 2001) for the limiting ratio of wing crack length to primary fracture length as a function of the far-field stress ratio given no friction along the primary fracture

16. The length ratios of wing cracks and the primary fracture obtained from the simulation is consistent with those from analytical solutions.

Appendix B: Model Validation for Steady-State Fluid Flow Through Persistent Fractures

The effective permeability of fractured porous media can be calculated by superimposing the contributions of matrix permeability and fracture permeability (Snow 1969):

$$k_{{{\text{eff}}}} = \frac{{(L - n_{{\text{f}}} b_{{\text{f}}} )k_{{\text{m}}} + n_{{\text{f}}} b_{{\text{f}}} k_{{\text{f}}} }}{L},$$

(26)

where L is the model domain dimension, n_f is the number of fractures along x or y directions, b_f is the fracture aperture, and k_m and k_f ($= b_{{\text{f}}}^{2} /12$) are the matrix and fracture permeabilities, respectively.

A numerical model containing two orthogonal sets of persistent fractures is constructed to validate the model for calculating the effective permeability of fractured porous media. The model has a dimension of L = 100 m, and the model contains 50 equidistant persistent fractures in each direction with the fracture aperture b_f varying from 0.01 to 1 mm. The matrix permeability is kept at k_m = 5 × 10^–19 m². Steady-state fluid flow through persistent fractures in each direction is simulated by imposing a macroscopic pressure differential at two opposite boundaries while the other two boundaries have a no-flow condition. An excellent match is achieved between the numerical model results and the analytical solution in both directions (see Fig. 

Figure17

Validation of the numerical model against the analytical solution (Snow 1969) for single-phase steady-state fluid flow in fractured porous media with two orthogonal sets of persistent fractures

17 for comparison in the x direction).