1 Introduction

Tectonically deformed coal is the coal whose chemical and mechanical properties have been dramatically altered by tectonic movements (Cheng and Pan 2020; Ju et al. 2014). Different from intact coal, it is usually characterized by a higher gas content, a higher gas desorption rate, a lower mechanical strength and a lower permeability (Black 2019; Wang et al. 2015; Wu et al. 2019; Zhao et al. 2015, 2019; Zhi and Elsworth 2016). Tectonically deformed coal is widely distributed in China, specifically, in the west of Henan Province, the south of Shanxi Province and the middle and north of Anhui Province, etc., where intense tectonic movements occurred during coal formation. According to the statistics, there are about 457 billion tons of tectonically deformed coal resources in China and up to 23.5% of coal production is from tectonically deformed coal seams (Cheng and Pan 2020).

The layered structure of coal seams comprising tectonically deformed sub-layers and intact sub-layers is often found in actual geological occurrence of coal (Lu et al. 2019a, b). Tectonically deformed coal seams caused by regional bedding sliding generally possess the combined vertical “intact-tectonically deformed-intact” layered structure (e.g., Daning Coal Mine in Shanxi Province) (Lu et al. 2019a, b) or “intact-tectonically deformed-intact” layered structure (e.g., Guhanshan Coal Mine in Henan Province) (Geological Report of Guhanshan Mine 2015; Zhang et al. 2019). In contrast, those caused by local faults tend to own an integrated continuous region of tectonically deformed coal around the fault areas (Fig. 1). The different combinations of tectonically deformed coal and intact coal results in their varying difficulty degrees of gas drainage. Therefore, commonly used original unified gas drainage methods like boreholes along the seam should be modified accordingly by amending parameters such as gas drainage time, gas drainage pressure and borehole spacing.

Fig. 1
figure 1

Typical combination forms of intact-tectonically deformed coal (Zhang et al. 2019; Lu et al. 2019a, b; Geological Report of Guhanshan Mine 2015)

In this study, layered coal seams comprising tectonically deformed sub-layers and intact sub-layers were taken as the research object. Three combined permeability forms, "high–low–high", "equal" and "low–high–low" representing "intact-tectonically deformed-intact", "intact" and "intact-tectonically deformed-intact" coal combinations, respectively, were considered. Besides, based on the proposed constitutive equations, the gas pressure distribution and gas drainage curve were numerically analyzed and verified by the actual drainage data. The results obtained are helpful to the design of gas drainage boreholes and the prevention of outbursts in coal seams with tectonically deformed sub-layers.

2 Constitutive equations

2.1 Equation of gas diffusion in matrixes

Coal is generally regarded as a heterogeneous dual-porosity medium consisting of matrixes and fractures (Fan and Liu 2018, 2019; Pan et al. 2010), and it is often simplified into an aggregation of sheets, match sticks and cubes in the simulation of gas transport in coal (Liu and Harpalani 2013). According to the sequence of gas transport, gas transports in coal by three steps: gas desorption from matrix surface to matrix pores, gas diffusion from matrix pores to fractures and gas flow in fractures.

In matrixes, gas exists in two states, namely the adsorbed state and the free state, of which the former obeys the Langmuir law, as written in Eq. (1):

$$X_{\text{x}} = \frac{{abp_{\text{m}} }}{{1 + bp_{\text{m}} }},$$
(1)

where, \(a\) is the gas adsorption capacity of coal, i.e. the Langmuir volume; \(b\) is the reciprocal of Langmuir pressure; \(p_{\text{m}}\) is the gas pressure in matrixes; and \(X_{\text{x}}\) is the gas adsorption amount.

The amount of free gas conforms to the following relationship:

$$X_{\text{y}} = \frac{{\phi_{\text{m}} p_{\text{m}} }}{{\rho p_{0} }},$$
(2)

where, Xy is the free gas content; \(\phi_{\text{m}}\) is the matrix porosity; \(p_{0}\) is the ambient gas pressure; and \(\rho\) is the gas density at standard temperature and pressure (STP).

The average gas concentration in matrixes is:

$$\overline{c}_{\text{m}} = \frac{{X_{\text{x}} + X_{\text{y}} }}{{V_{\text{m}} }} = \left( {\frac{{abp_{\text{m}} }}{{1 + bp_{\text{m}} }} + \frac{{\phi_{\text{m}} p_{\text{m}} }}{{\rho p_{0} }}} \right) \cdot \frac{\rho \text{M}}{{V_{\text{m}} }},$$
(3)

where, \(\overline{c}_{\text{m}}\) is the average gas concentration in matrixes; \(V_{\text{m}}\) is the matrix volume; and M is the molar mass of gas.

According to the Fick’s description, diffusion is a mass transfer process driven by gas concentration difference:

$$q_{\text{m}} = D\tau V_{\text{m}} \left( {\overline{c}_{\text{m}} - c_{\text{f}} } \right)$$
(4)

where, \(q_{\text{m}}\) is the mass flux; \(D\) is the diffusion coefficient; \(\tau\) is the shape factor of matrixes; and \(c_{\text{f}}\) is the gas concentration in fractures.

Gas pressure and gas concentration can be interconverted using the equation of state of real gas:

$$\left\{ \begin{gathered} c_{\text{m}} = \frac{M}{{Z_{\text{m}} RT}}p_{\text{m}} \hfill \\ c_{\text{f}} = \frac{M}{{Z_{\text{f}} RT}}p_{\text{f}} \hfill \\ \end{gathered} \right.$$
(5)

where, \(p_{\text{f}}\) is the gas pressure in fractures; \(Z_{\text{m}}\) and \(Z_{\text{f}}\) whose values can be assumed to be 1 are the coefficients of compressibility for gas in matrixes and fractures, respectively; \(R\) is universal gas constant, 8.314 J/(mol K); and \(T\) is the gas temperature, K.

Hence, Eq. (6) can be obtained by integrating Eq. (4) with Eq. (5):

$$ q_{\text{m}} = \frac{{D\tau MV_{\text{m}} }}{RT}\left( {p_{\text{m}} - p_{\text{f}} } \right) $$
(6)

Gas in matrixes diffuses into fractures, and the mass conservation equation of gas in matrixes is:

$$q_{\text{m}} = - \frac{{\partial \overline{c}_{\text{m}} V_{\text{m}} }}{\partial t}$$
(7)

Equation (8) can be obtained by substituting Eqs. (3) and (6) into Eq. (7):

$$\frac{{\partial p_{\text{m}} }}{\partial t} = - \frac{{D\tau \left( {p_{\text{m}} - p_{\text{f}} } \right)}}{{\left( {\frac{ab}{{1 + bp_{\text{m}} }} - \frac{{ab^{2} p_{\text{m}} }}{{\left( {1 + bp_{\text{m}} } \right)^{2} }} + \frac{{\phi_{\text{m}} }}{{\rho p_{0} }}} \right)\frac{\rho RT}{{V_{\text{m}} }}}}$$
(8)

2.2 Equation of gas flow in fractures

Gas flow in fractures is widely believed to obey the Darcy law (Pan et al. 2010), which can be expressed as:

$$v_{\text{g}} = \frac{{k_{\text{f}} }}{\mu }\nabla p_{\text{f}}$$
(9)

where, \(v_{\text{g}}\) is the gas flow rate; kf is the gas permeability; and \(\mu\) is the gas dynamic viscosity.

On the one hand, gas in fractures will leave fractures due to the pressure difference; on the other hand, gas in matrixes will supply free gas to fractures due to gas desorption and diffusion. Thus, the mass conservation equation of gas in fractures can be written as:

$$\frac{{\partial m_{\text{f}} }}{\partial t} - \nabla \cdot \left( {\rho v_{\text{g}} } \right) = q_{\text{m}}$$
(10)

where, \(m_{\text{f}}\) is the mass of free gas in fractures.

Substituting Eqs. (6) and (9) into Eq. (10), the governing equation for gas flow in fractures can be obtained:

$$\phi_{\text{f}} \frac{{\partial \left( {p_{\text{f}} } \right)}}{\partial t}{ + }p_{\text{f}} \frac{{\partial \left( {\phi_{\text{f}} } \right)}}{\partial t} - \frac{k}{\mu }\nabla \cdot \left( {p_{\text{f}} \cdot \nabla p_{\text{f}} } \right) - D\tau \cdot \left( {1 - \phi_{\text{f}} } \right)\left( {p_{\text{m}} - p_{\text{f}} } \right) = 0$$
(11)

where, \(\phi_{\text{f}}\) is the fracture porosity.

The above constitutive equations can be interconnected with the aid of the PDE module in COMSOL Multiphysics (Fig. 2).

Fig. 2
figure 2

Constitutive equations

3 Geometric model and boundary setting

Based on the three cases listed in Fig. 1, a simplified geometric model was established (Fig. 3). The model contains two sections, Zone I and Zone II, representing areas with drainage boreholes and areas without drainage boreholes, respectively. The combination type of intact coal and tectonically deformed coal is realized by changing coal permeability in Zones I and II. The permeability ratios of coal in Zones II and I (kII/kI) are 0.1, 1 and 10 for Cases 1, 2 and 3, respectively. The geometric model is a 3-m-long and 1-m-high rectangle whose boundaries were set to be zero flux. In addition, a borehole with the diameter of 100 mm was arranged in Zone II, and its boundaries were set to be −13 kPa pressure. After the simulation, the drainage curve of gas mass was plotted by taking the integral of gas pressure difference in the studied area. Moreover, Line OO’ (0, 0.25; 3, 0.25) and Point A (1, 0.5) were set to monitor the local change of gas pressure. It is noteworthy that the mechanical module was not introduced into the simulation because the permeability ratio (kII/kI) needs to be constant for comparison. All the physical parameters used are given in Table 1.

Fig. 3
figure 3

Simplified geometric model and boundary settings

Table 1 Parameter settings (An et al. 2013; Gong et al. 2019)

4 Result and discussion

4.1 Overall gas pressure distributions with different permeability ratios kII/kI

Figure 4 displays the simulation results of overall gas pressure with different permeability ratios under the same drainage pressure. In Case 1 where kII/kI = 0.1, gas pressure mainly decreases around the borehole in Zone II, exhibiting a clear lamp-cover-shaped distribution, while no obvious gas pressure gradient is obtained in the high-permeability coal in Zone I (with the same color). Regarding Case 2 where kII/kI = 1, the area of pressure decay expands from Zone I to Zone II. An obvious gas pressure gradient appears in Zone II, which is similar to results of gas pressure distribution in homogeneous coal with a single borehole. In contrast, Case 3 where kII/kI = 10 demonstrates a more obvious partition property than Cases 1 and 2. In this case, the blue area expands dramatically in Zone II, indicating the formation of a larger range of gas pressure decay compared with the previous cases. However, because of the low permeability in Zone I, the rate of gas pressure decay decelerates in Zone I, forming a clear strip-shaped gas pressure distribution. Considering the range of gas pressure decay, it can be concluded that the "low–high–low" permeability combination described in Case 3 brings about the best gas drainage performance, whereas the "high–low–high" permeability combination described in Case 1 leads to the greatest difficulty in gas drainage.

Fig. 4
figure 4

Cloud map of gas distribution in the whole area

4.2 Local gas pressure distributions with different permeability ratios kII/kI

Figure 5 presents the simulation results of gas pressure distribution at Line OO′ for 50, 100, 150 and 200 days, respectively. It can be found that the local pressure distributions in the three cases resemble the overall results displayed in Fig. 4. In Case 1 where kII/kI = 0.1, the gas pressure decay area mainly concentrates around the borehole. Specifically, the shapes of gas pressure decay in Zone I are nearly parallel straight lines, while those in Zone II are curved lines with high tangency slopes. The reason is as follows: Due to the much lower permeability of coal in Zone II, gas can hardly be desorbed from the coal and remains at a high pressure for a long time. On the contrary, gas in the low-permeability coal in Zone I can easily move away from fractures, resulting in a small pressure gradient and straight lines of gas pressure decay along the x axis. In Case 2 where kII/kI = 1, no obvious partition signs are observed, and the gas pressure distribution curves are relatively smooth in a single-butterfly shape along the x axis. In Case 3 where kII/kI = 10, a clear partition line exists between Zone I and Zone II. In each zone, the curves are distributed in a dual-butterfly shape along the x axis. The reason is as follows: Gas in the high-permeability coal in Zone II can easily move away, but Zone I fails to supply gas to Zone II in time owing to its lower permeability. Eventually, a discontinuous pressure distribution different from the one in Case 2 is formed. Hypothetically, Zone II acts as a borehole with an expanded diameter of 1 m which is much larger than the initially set diameter (100 mm) of the simulated borehole. Hence, the pressure distribution is supposed to feature dual stages, consisting of pressure decay curves for a smaller borehole and a hypothetical larger borehole.

Fig. 5
figure 5

Gas pressure distribution curves at Line OO

The roles of Zone I and II in gas pressure distribution curves can be obtained in light of the theory of water transport through pipes. In Case 1, Zone II acts as a valve that limits the "water flux" from Zone I. It determines the maximum flow rate of gas and thus is the main controller in this situation (Fig. 5d). On the contrary, if the "water flux" cannot reach the maximum flux that pipes can provide, which may result in an under-voltage phenomenon, the main controller would be the water flux itself. As a result, in Case 3, Zone I should be the main controller of gas flow (Fig. 5f). In Case 2, the roles of Zone I and Zone II are ideally matched, which means the actual water flux basically equals the maximum water flux that the valve can provide. In this situation, Zone 1 and Zone II equally control the gas flow (Fig. 5e).

4.3 Gas drainage curves with different permeability ratios kII/kI

To better illustrate the variation of gas drainage volume with time, the pressure at Point A (1, 0.5) was monitored during simulation (Fig. 6a). It can be found that the decreasing rate and range of residual pressure at Point A both rise with the increase of kII/kI. The relationship between the gas drainage volume and the time varies from a straight line to a dual-stage curve consisting of a fast pressure-decay curve and a slow one.

Fig. 6
figure 6

Residual gas pressure at Point A and cumulative gas drainage volume for whole area with different permeability ratios

The first-order kinetic model is extensively employed to describe the gas desorption behavior of coal (Airey 1968):

$$M\left( t \right) = M_{\infty } \exp \left( { - k t} \right),$$
(12)

where, \(M\left( t \right)\) and \(M_{\infty }\) are the amounts of desorbed gas mass at times t and t, respectively; k is the first-order kinetic model constant; and t is the time.

Similar to the mass-based first-order kinetic model, a pressure-based first-order kinetic model Eq. (13) is established based on the equation of state of gas (Busch et al. 2004; Zhao et al. 2014):

$$P\left( t \right) - P_{\infty } = \left( {P_{0} - P_{\infty } } \right)\exp \left( { - k \cdot t} \right),$$
(13)

where, \(P\left( t \right)\), \(P_{0}\) and \(P_{\infty }\) are the gas pressures at times t, t0 and t, respectively. Equations (11) and (12) are appropriate for describing gas desorption in a single uniform system. If the gas pressure variation is characterized by dual systems, i.e. tectonically deformed coal and intact coal, Eq. (13) can be evolved into Eq. (14):

$$\begin{gathered} Q_{r} \left( t \right) = \frac{{P\left( t \right) - P_{\infty } }}{{P_{0} - P_{\infty } }} = \frac{{P_{1} \left( t \right) + P_{2} \left( t \right) - P_{1\infty } - P_{2\infty } }}{{P_{10} + P_{20} - P_{1\infty } - P_{2\infty } }} \hfill \\ \begin{array}{ll} {} & {} \\ \end{array} = \frac{{\left[ {P_{1} \left( t \right) - P_{1\infty } } \right]/\left( {P_{10} - P_{1\infty } } \right)}}{{1 + \left( {P_{20} - P_{2\infty } } \right)/\left( {P_{10} - P_{1\infty } } \right)}} + \frac{{\left[ {P_{2} \left( t \right) - P_{2\infty } } \right]/\left( {P_{20} - P_{2\infty } } \right)}}{{\left( {P_{10} - P_{1\infty } } \right)/\left( {P_{20} - P_{2\infty } } \right) + 1}} \hfill \\ \begin{array}{*{20}c} {} & {} \\ \end{array} = \frac{\zeta }{1 + \zeta }\frac{{P_{1} \left( t \right) - P_{1\infty } }}{{P_{10} - P_{1\infty } }} + \frac{1}{1 + \zeta }\frac{{P_{2} \left( t \right) - P_{2\infty } }}{{P_{20} - P_{2\infty } }} = \lambda \frac{{P_{1} \left( t \right) - P_{1\infty } }}{{P_{0} - P_{1\infty } }} + \left( {1 - \lambda } \right)\frac{{P_{2} \left( t \right) - P_{2\infty } }}{{P_{0} - P_{2\infty } }} \hfill \\ \begin{array}{*{20}c} {} & {} \\ \end{array} { = }\lambda \exp \left( { - k_{1} \cdot t} \right) + \left( {1 - \lambda } \right)\exp ( - k_{2} \cdot t) \hfill \\ \end{gathered}$$
(14)

where, \(Q_{r}\) is the decreasing percentage of gas pressure; \(P_{1} \left( t \right)\), \(P_{2} \left( t \right)\), \(P_{10}\), \(P_{20}\), \(P_{1\infty }\) and \(P_{2\infty }\) are the gas pressures in intact coal and tectonically deformed coal at times t, t0 and t, respectively; \(\lambda\) is the contribution from intact coal to the decrease of gas pressures,\(\lambda { = }\zeta /\left( {1 + \zeta } \right)\); \(\zeta\) is a constant, \(\zeta { = }\left( {P_{20} - P_{2\infty } } \right)/\left( {P_{10} - P_{1\infty } } \right)\); and \(k_{1}\) and \(k_{2}\) are the first-order kinetic constants for intact coal and tectonically deformed coal, respectively.

Table 2 displays the fitting results of data obtained in Fig. 6a using Eq. (14). In Case 1, the fitting results show a single-phase characteristic of a zero contribution from intact coal (Zone II with a high permeability). Moreover, the first-order kinetic constants fitted for the two systems are of the same value (4.9 × 10–4), with an equally slow rate of mass transport. In Case 2, the contribution from intact coal rises to 47%, almost equal to the contribution from tectonically deformed coal (53%). With respect to the kinetic constant, \(k_{1}\) is about two orders of magnitude larger than \(k_{2}\). In Case 3, the fitting results exhibit a clear dual-phase characteristic. The contribution from intact coal is about 4 times larger than that from tectonically deformed coal. Besides, \(k_{1}\) (2.5 × 10–2) is much larger than \(k_{2}\) (1.0 × 10–4).

Table 2 Fitting results using Eq. (13)

Figure 6b shows the cumulative gas drainage volume by taking the surface integral of gas pressure. When kII/kI= 10, the curve is of the largest value and the greatest initial slope. As kII/kI decreases to 0.1, the curve, which approximates a straight line, is of the smallest value and the smallest initial slope. It can be inferred that, if the initial gas pressure is the same, completing gas drainage in the studied area consumes the longest time for Case 1 and the shortest time for Case 3. Thereby, it is suggested to design proper drainage methods for different gas occurrence conditions. For the situation in Case 1, longer drainage time, lower drainage pressure and smaller borehole spacing should be arranged in hope of achieving a satisfactory drainage effect.

5 Drainage method optimization for coal with tectonically deformed sub-layers

As revealed by the above analysis on gas drainage difficulties in coal seams with different combinations of tectonically deformed sub-layers and intact sub-layers, proper drainage methods should be adopted to ensure a satisfactory gas drainage volume in a fixed period, so that underground mining can be carried out at the scheduled rate. The drainage methods can be divided into two categories: commonly used drainage methods and enhanced drainage methods. The former, such as boreholes along the seam and boreholes crossing the seam, require longer drainage time, lower drainage pressure, smaller borehole spacing, larger borehole radius, etc., to achieve a better drainage performance. However, if these methods are not effective or economically allowable, enhanced drainage methods, such as long boreholes along the seam and hydraulic flushing boreholes, should be adopted to drain gas.

Figure 7 illustrates the schematic diagram of drainage methods used in Daning Coal Mine in Shanxi Province and Guhanshan Coal Mine in Henan Province, China. In Daning Coal Mine that possesses a typical "intact-tectonically deformed-intact" coal combination, long boreholes along the seam were constructed to drain gas in coal seams. The long boreholes can cover a huge area of coal seam and enhance coal permeability by means of initial water injection, conducing to efficient gas drainage in soft coal. Moreover, for the thickened tectonic sub-layers, branch-holes were supplemented to promote the gas drainage performance.

Fig. 7
figure 7

Schematic diagram of drainage methods used in Daning Coal Mine and Guhanshan Coal Mine, China (Lu et al. 2019a, b; Wang 2016; Zhang et al. 2017)

In Guhanshan Coal Mine which has a typical "intact-tectonically deformed" coal combination ("Tectonically deformed-intact-tectonically deformed" coal combination also exists locally.), a newly developed method of hydraulic flushing boreholes was applied to gas drainage in tectonically deformed sub-layers. These boreholes can enhance the coal permeability and the drainage radius in targeted places by flushing the soft tectonically deformed coal with high-pressure water. After applying the hydraulic flushing technology, the permeability difference between tectonic sub-layers and intact sub-layers is reduced, gradually forming a combination type described in Case 2. Accordingly, the pressure decay curves change from dual-stage ones to nearly one-stage ones. Figure 8 discloses a comparison between gas concentrations in hydraulic flushing boreholes and ordinary boreholes in literature. The gas concentration in ordinary boreholes falls at distinct rates in two stages (the fast stage and the slow stage), while that in hydraulic flushing boreholes declines as a smooth single-stage curve instead of showing dual-stage characteristics. Meanwhile, compared with the gas concentration in ordinary boreholes, that in hydraulic flushing boreholes corresponds to a reduced tangent curve slope and dramatically increased values, suggesting that shorter drainage time is needed if the gas volume to be drained is fixed. The result verifies correctness of the abovementioned numerical analysis.

Fig. 8
figure 8

Curves of gas concentration in hydraulic flushing boreholes and ordinary boreholes (Zhang et al. 2019)

6 Conclusions

In this study, the gas pressure distributions in layered coal seams consisting of tectonically deformed sub-layers and intact sub-layers were analyzed. The main conclusions are:

  1. (1)

    Different combinations of tectonically deformed coal and intact coal determine the varying pressure distributions in coal seams. In Case 1 where kII/kI= 0.1, gas pressure mainly decreases around the borehole. The shapes of gas pressure decay in the high-permeability zone are nearly parallel straight lines, while those in the low-permeability zone are curves with large tangency slopes. In Case 2 where kII/kI= 1, the area of pressure decay expands, and the gas pressure distribution curves are distributed in a single-butterfly shape. In Case 3 where kII/kI= 10, a clear partition line can be observed, and the curves are distributed in a dual-butterfly shape.

  2. (2)

    The simulated gas drainage curve varies from a single-stage line to a dual-stage curve as the permeability ratio kII/kI grows. For the smallest kII/kI, it takes the longest time to drain all the gas in the studied area. Therefore, longer drainage time, lower drainage pressure and smaller borehole spacing or other enhanced drainage methods should be employed for achieving a satisfactory drainage effect.

  3. (3)

    In-situ gas drainage data for hydraulic flushing boreholes and ordinary boreholes exhibit a clear dual-stage characteristic for the coal seam with tectonically deformed sub-layers and a single-stage characteristic for the permeability-enhanced coal seam, which verifies correctness of the simulated pressure decay curves.