Taking the thermal expansion/contraction and sorption-induced swelling/shrinkage of the coal matrix into consideration, the component of the total strain tensor for a non-isothermal dual-porosity coal can be represented by (Wu et al., 2010; Liu et al., 2011; Ge et al., 2019) Eq. 1: ε i j = 1 2 G σ i j − ( 1 6 G − 1 9 K ) σ k k δ i j + α p m δ i j + β p f δ i j + K ( α T T δ i j + ε s δ i j ) 3 K (1) Where G = D / 2 ( 1 + υ ) ; D = [ 1 / E + 1 / a K n ] − 1 K = D / 3 ( 1 − 2 υ ) ; α = 1 − K / K s ; β = 1 − K / a ⋅ K n ; δ _ij = 1 (i = j) or δ _ij = 0 (i≠j) is the Kronecker tensor; G is the shear modulus of coal [MPa]; K is the bulk modulus of coal [MPa]; α and β are the Biot coefficients [-] for fracture and coal matrix, respectively; p _f and p _m are methane pressure [MPa] in the fracture and the coal matrix, respectively; α _T is the coefficient of thermal expansion [1/K]; T is the reservoir temperature [K]; ε s is matrix deformation induced by adsorption/desorption, [Note: this variable could be represented by the typical Langmuir-type equation ε s = ε L p m / ( p m + p L ) ]; ε L is the Langmuir strain constant for the matrix [-]; p L the Langmuir pressure constant [MPa]; E is the elastic modulus of coal [MPa]; υ is the Poisson’s ratio [-]; K _s is the bulk modulus of coal grains [MPa]; K _n is the normal stiffness of fracture [MPa/m].

According to continuum mechanics theory, the stress equilibrium equation and the strain–displacement relationships for coal deformation can be calculated from Eq. 2: (Zhao Y. et al., 2020): { σ i j , j + f i = 0 ε i j = 1 2 ( u i , j + u j , i ) (2) where f i is the body force component and u i is the displacement component.

By substituting Eq. 1 into Eq. 2, the Navier-type equation that governs coal deformation, Eq. 3, is obtained: G u i , j j + G 1 − 2 ν u j , i i − α p m , i − β p f , i − K α T T , i − K ε L P L ( p m + P L ) 2 p m , i + f i = 0 (3)

Methane migration in a coal seam during drainage obeys the law of mass conservation the equation for which can be expressed as (Zhang et al., 2008): ∂ m ∂ t + ∇ ⋅ ( ρ g q g ) = Q S (4) where t is the time of methane flow [s]; ρ g is the methane density [kg/m³]; q _g is the seepage velocity vector [m/s]; Q _S is the methane flow rate from its source [kg/(m³∙s)]; m is the total methane mass including both free gaseous methane and adsorbed methane [kg/m³]. Assuming that methane sorption only occurs in the matrix material, the methane stored in the matrix and fracture can be written as (Wu et al., 2010; Zhang et al., 2021): { m m = ϕ m p m ρ g m p a + ρ g m ρ c V L p m p m + P L ( 1 − A − W ) m f = ϕ f ρ g f (5) Where m _m and m _f are the methane mass [kg/m³] in the coal matrix and fractures, respectively; ϕ m is the porosity of the matrix system [-]; ϕ f is the porosity of the cleat system [-]; ρ c is the coal density [kg/m³];V _L is the Langmuir pressure constant [m³/kg]; A is coal ash [%]; W is the water content in the coal [%]; p a is standard atmospheric pressure [MPa]; ρ g m and ρ g f are the methane density in the coal matrix and fractures, respectively. These densities can be expressed by Eq. 6: { ρ g m = M g R T p m ρ g f = M g R T p f (6) where M _g is the molecular mass of CH₄ [kg/mol], and R is the universal gas constant [J/(mol•K)].

Taking into consideration the equations for fluid exchange between matrix and fractures (Lim and Aziz 1995) and with reference to previous research results (Ge et al., 2019), the governing equations for methane flow in the matrix–fracture system could be calculated from Eq. 7: { ∂ m m ∂ t + ∇ ( − k m μ ρ g m ∇ p m ) = − ρ g m k m ψ μ ( p m − p f ) ϕ f ∂ ρ g f ∂ t + ρ g f ∂ ϕ f ∂ t + ∇ ( − k f μ ρ g f ∇ p f ) = ρ g m k m ψ μ ( p m − p f ) (7) Where k _m is the matrix permeability [m²]; k _f is the fracture permeability [m²]; μ is the dynamic viscosity of methane [Pas]; ψ = 4 ( 1 / a x 2 + 1 / a y 2 ) is the shape factor [1/m²]; a x and a y are the matrix spacing in the x- and y-directions [m].

Coal seams contain fluid and solid phases whose thermodynamic parameters are completely different. Any pre-existing geothermal equilibrium in a coal seam is disrupted during CMM extraction. By applying the law of energy conversation, the heat transfer in a dual-porosity coal seam under non-isothermal conditions has been described in the literature (Zhou et al., 1998; Zhu et al., 2011; Ge et al., 2019) and can be expressed by Eq. 8: ( p C ) M ∂ T ∂ t + T K g α g ( ∇ ( − k m μ ∇ p m ) + ∇ ( − k f μ ∇ p f ) ) + T K α T ∂ ε v ∂ t = λ m ∇ 2 T + ρ g m k m C g μ ∇ p m ∇ T + ρ g f k f C g μ ∇ p f ∇ T (8) Where ( p C ) M = ϕ m ( ρ g m C g ) + ϕ f ( ρ g f C g ) + ( 1 − ϕ m − ϕ f ) ( ρ c C s ) is the effective heat capacity of methane-bearing coal; λ M = ( ϕ m + ϕ f ) λ g + ( 1 − ϕ m − ϕ f ) λ s is the heat conductivity coefficient for methane-bearing coal [J/(m∙s∙K)]; λ _s is the heat conductivity coefficient for the coal’s skeleton [J/(m∙s∙K)]; λ _g is the heat conductivity coefficient for methane [J/(m∙s∙K)]; C _g and C _s are the specific heat capacity [J/(kg∙K)] of the methane and the coal skeleton, respectively; K _g is the bulk modulus of methane [MPa], and α g = 1 / T is the thermal expansion coefficient of methane [1/K].

Coal is a fractured rock with permeability arising from both pores in the coal matrix and fractures that cut through the coal (Zhang et al., 2018; Cheng et al., 2021; Du et al., 2021; Lv et al., 2021). Since the impact of thermal changes effects on the coal’s permeability cannot be ignored (Perera et al., 2012; Li et al., 2020), we extended the matrix porosity and permeability model proposed by Zhang et al. (2008) to encompass the effect of the temperature on permeability. Using this extended model, the dynamic evolution of matrix porosity and permeability evolution during temperature changes can be calculated from the following Eq. 9: { ϕ m = 1 1 + S [ ( 1 + S 0 ) ϕ m 0 + α ( S − S 0 ) ] k m = k m 0 ( 1 1 + S [ ( 1 + S 0 ) + α ϕ m 0 ( S − S 0 ) ] ) 3 (9) where S = ε v + p m / K s − ε s − α T T ; S 0 = ε v 0 − p m 0 / K s − ε L p m 0 / p m 0 + P L − α T T 0 ; k _m0 is the initial coal matrix permeability at the initial methane pressure p _m0 and matrix porosity ϕ _m0 .

Expanding on the results of the previous research (Liu et al., 2010), the governing equation for the coal’s fracture permeability, an equation that incorporates sorption-induced swelling/shrinkage, thermal expansion, and effective stress, can be written as: k f i k f 0 = [ 1 + 2 ( 1 − R m ) ϕ f 0 ( Δ ε j − 1 3 α T Δ T − 1 3 Δ ε s ) ] 3 , i ≠ j (10) Where R _m = E/E _S is the reduction ratio of the elastic modulus [-]; ϕ _f0 and k _f0 are the initial fracture porosity and initial fracture permeability, respectively.

Equations 3, 7–10 constitute the THM coupling model used to simulate methane transport in a coal seam during CMM drainage. Figure 2 depicts the interactions and relationships simulated by the THM model. When CMM suction drainage begins using negative pressure, the free methane and heat in the fractures around the boreholes will rush out first. After that, the methane in the matrix will gradually be drained. The coal will be deformed by the change in effective stress and temperature and the coal matrix and fracture permeabilities, and thus methane flow will respond to these changes. The THM model governing equations are interdependent and interrelated, as shown in Figure 2. For this study, the whole set of coupled equations were solved discretely with the aid of Comsol Multiphysics (COMSOL Inc., Burlington, MA, United States), a finite element solver.

To select the optimum layout for multiple tree-type CMM drainage boreholes, the methane drainage efficiency of several possible tree-type sub-borehole configurations were compared. Figure 3 shows possible sub-borehole configurations for simulated methane drainage from nine tree-type boreholes. This study modeled four sub-boreholes from each main borehole. The sub-borehole configurations are a square layout, a rhombic layout, a parallel-mixed layout and a cross-mixed layout. The numerical models for the borehole layouts in Figure 3 were 50 by 50 m. The main tree-type boreholes were spaced 6 m apart and each main tree-type borehole had four 3 m long sub-boreholes. The main boreholes and the sub-boreholes were 94 and 25 mm in diameter, respectively. The boundary conditions for the models shown in Figure 3 were defined in accordance with the stresses on the Ji_15-17 seams in the Shoushan mine in Henan province, China. To reduce the number of main pre-drainage boreholes required in the Ji_15-17 seams and improve methane drainage, the Shoushan mine implemented a full field test of the TTBD method in one of its floor-level methane drainage roadways. The details of the geology of the Ji_15-17 seams are introduced in a previous study (Zhang et al., 2021). In the models constructed for the present study, the maximum horizontal stress, 19.5 MPa, was imposed on the models’ upper surfaces and the minimum horizontal stress, 12.5 MPa, was imposed on their right sides. The other sides were defined as roller boundaries. The coal was initially saturated with methane at a pressure of 1.38 MPa. A constant methane pressure of 116 mm Hg (∼15.5 kPa) was applied to all the tree-type boreholes; no flow conditions were applied to the other boundaries. The initial reservoir temperature was set to 323 K and the tree-type borehole walls were defined as temperature boundaries with a temperature of 293 K. The numerical models’ input parameters are given in Table 1 (Zhu et al., 2011; Xia et al., 2014; Li et al., 2016; Zheng et al., 2016; Zhang et al., 2021).

Figure 4 shows methane pressure in a coal seam at different times during modeled CMM drainage for different tree-type sub-borehole configurations. The letters S, R, P, and C designate the square, rhombic, parallel mixed and cross mixed sub-borehole layouts, respectively, and the “xxd” labels show the number of days since gas extraction was initiated. For example, the diagram labeled “R215d” shows methane pressure for a rhombic sub-borehole layout after 215 days of drainage. From Figure 4, it can be seen that methane pressure near the tree-type boreholes gradually decreases and the drainage zone expands as the gas extraction time increases. When the drainage time is less than the completion time, there are zones between tree-type sub-boreholes where the methane has not yet been drained (the red zones in Figure 5). As the drainage time lengthens, the methane pressure in these zones will gradually decrease and eventually it will fall below 0.74 MPa. When this occurs, CMM extraction has been completed. For different sub-borehole configurations, although the number of sub-boreholes is the same, the times required to complete extraction differ. The times required to drain a zone completely can differ by up to 44.4%. The differences in completion times can be clearly observed from Figures 4, 6. It is worth noting that if the drainage completion time is used as the only assessment criterion for drainage effectiveness, the rhombic tree-type borehole layout is better than the square, parallel mixed, and cross mixed layouts. These model results suggest that optimizing the tree-type sub-borehole layout can improve drainage performance. To check the accuracy of the THM model results, multiple tree-type boreholes with rhombic sub-borehole layouts were used to drain CMM from the Ji_15–17 seams in the Shoushan mine. After 198 days of gas extraction, the methane pressure between tree-type boreholes dropped to less than 0.74 MPa. The relative error in our simulated drainage completion time is 8.72%. Figure 7 shows field production data for methane extraction from the rhombic-configured tree-type boreholes in the Shoushan mine and simulated methane production results from the rhombic THM model. It can be seen that the simulated results and the field methane data match very well. These results show that the THM model is valid and it can be used to simulate tree-type borehole methane drainage accurately. The model can also be used to determine the most appropriate tree-type sub-borehole layout.

The design and optimization of tree-type borehole spacing for methane pre-drainage from coal reservoirs is an important issue that needs to be resolved. Drilling many closely spaced boreholes from an underground roadway into a coal seam is expensive and requires long construction times. A practical spacing would minimize the number boreholes while still achieving excellent methane drainage effect. Figure 8 displays the sizes of the effective drainage zones around multiple tree-type boreholes at different spacings for extraction times from zero to 2,000 days. The gas extraction completion times for each borehole spacing are also shown. From Figure 8, it is apparent that the size of the effective drainage zone for any borehole spacing first rises dramatically as the drainage time increases but then the zone size increases much more slowly after the pre-drainage completion time is reached. In a certain period in the early drainage stages, the shorter the tree-type borehole spacing, the larger the drainage zone. However, at longer drainage times, this relationship is reversed. The reason for this reversal is that the methane in the coal between closely spaced boreholes can be rapidly drained due to the strong interactions between the boreholes but with the passage of drainage time, drainage zone growth rates are restricted by the tree-type borehole spacing and the coal’s low permeability. To explore the relationship between tree-type borehole spacing and drainage efficiency, tree-type borehole spacings are shown plotted against drainage completion times in Figure 9. The figure shows that the time required to complete drainage differs considerably for different tree-type borehole spacings; the tree-type borehole spacing is related to the pre-drainage completion time by a power law. When the distance between the tree-type boreholes is increased from 6 to 10 m, the required drainage time increases from 215 to 1,522 days, an increase of more than six times. Thus, an appropriate distance between tree-type boreholes should be determined to obtain the most efficient CMM extractions possible. Doing so will improve the efficiency of mining operations.

Actually, to promote pre-drainage construction when a coal face is not being worked and to eliminate the risk of coal and gas outbursts, coal mining companies will have to allocate sufficient lead time to allow for adequate methane drainage. The lead time is usually longer than the drainage completion time. How borehole spacing is designed is significantly important because the required lead time at the Shoushan mine is approximately 400 days. The distribution of methane pressure on the measuring line from point A (0, 25 m) to point B (50, 25 m) after 400 days of extraction is discussed. As shown on Figure 10, the larger the tree-type borehole spacing, the higher methane pressure in the coal seam. For the tree-type borehole spacing of 6 or 7 m, the methane pressure between the tree-type boreholes has been successfully reduced to below 0.74 MPa after 400 days of drainage. Thus, taking the drainage effect and the lead time allowed into account, a diamond-shaped sub-borehole layout in tree-type boreholes spaced 7 m apart could be recommended for extracting the CMM from Ji_15–17 coal seams in the Shoushan coal mine. These THM model simulation results can provide reliable support for optimizing the layout of tree-type boreholes to improve the efficiency of methane drainage programs in the field.

Figure 11 gives the influence of different Langmuir volume constants, V _L , on methane drainage performance of tree-type boreholes. After draining for the similar amount of time, the effective drainage zone areas decrease as the Langmuir volume constant increases (Figure 11A). Meanwhile, it can be found in Figure 11 that coal with different Langmuir volume constants reach drainage completion after different periods of time. The greater the Langmuir volume constant, the more time is needed to complete CMM extraction. The relationship between drainage completion time and the Langmuir volume constant is linear (Figure 11B). As the Langmuir volume constant is changed from 0.01 to 0.05 m³/kg, the time required for pre-drainage completion increases by a factor of approximately 3.9. The reason is that a higher Langmuir volume constant indicates that more methane is adsorbed and stored in the coal and this means there is more methane to be extracted.

The Langmuir pressure constant P _L , is another parameter used to characterize adsorption. It is the pore pressure corresponding to one-half of the Langmuir volume (V _L /2). A graph showing how different P _L ’s change effective drainage zone sizes and drainage completion times is shown as Figure 12. For the same drainage time, as the Langmuir pressure constant increases, the effective drainage zones expand (Figure 12A), but the final drainage completion times show a downward trend (Figure 12B). The reduction in the drainage time–P _L relationship follows a power law (Figure 12B). The time required for pre-drainage is reduced by 30.4 and 41.1% when the Langmuir pressure constant increases from 0.5 to 1.729 and 2.5 MPa, respectively. This occurs because the higher Langmuir pressure constants increase the gas desorption rate during drainage to some extent and this shortens the time required for pre-drainage. To summarize, the drainage performance of multiple boreholes is affected by the adsorption factors, therefore, it is fairly important that simulations consider the adsorption and that the coal’s adsorption parameters be measured accurately.

Figure 13 displays the changes in the drainage performance for different initial fracture permeabilities, k _f0’s. As show in the figure, for the same amount of drainage time, an increase of the fracture permeability leads, as expected, to a larger effective drainage zone and the time required to complete the entire pre-drainage process decreases. The drainage completion times fall rapidly as k _f0 increases in the lower ranges of the initial fracture permeability considered whereas the drainage times decrease much more slowly in the higher initial fracture permeability ranges. The time required for pre-drainage is reduced by 92.2 and 58.8% when the initial fracture permeability is changed from 2 × 10⁻¹⁸ m² to 4 × 10⁻¹⁷ m² and from 4 × 10⁻¹⁷ m² to 1 × 10⁻¹⁶ m², respectively. A power function has been fitted to the drainage completion time–initial fracture permeability data (Figure 13B). Compared with fracture permeability, the effect of matrix permeability on the time required for pre-drainage is relatively minor. As Figure 14A shows, higher matrix permeabilities, k _m0 ’s, make the effective drainage zone for any given drainage time slightly larger. This results in a slow and linear downward trend in drainage completion times on the drainage completion time–initial matrix permeability plot (Figure 14B). The time needed to complete methane extraction decreases from 226 d for a matrix permeability of 2 × 10⁻¹⁸ m² to only 142 d for a matrix permeability of 1 × 10⁻¹⁷ m². Drainage performance is less sensitive to changes in matrix permeability than to changes in fracture permeability. The reason for this difference is that the fractures are the dominant pathways for methane migration. Indirectly, this also demonstrates that the TTBD method for enhancing fracture permeability through stress relief and the addition of sub-boreholes to connect the main boreholes to the fractures in the seams should be effective for CMM drainage.