Mathematical Modeling and Simulation on the Stimulation Interactions in Coalbed Methane Thermal Recovery

: Heat stimulation of coalbed methane (CBM) reservoirs has remarkable promotion to gas desorption that enhances gas recovery. However, coalbed deformation, methane delivery and heat transport interplay each other during the stimulation process. This paper experimentally validated the evolutions of gas sorption and coal permeability under variable temperature. Then, a completely coupled heat-gas-coal model was theoretically developed and applied to a computational simulation of CBM thermal recovery based on a ﬁnite element approach of COMSOL with MATLAB. Modeling and simulation results show that: Although di ﬀ erent heat-gas-coal interactions have di ﬀ erent e ﬀ ects on CBM recovery, thermal stimulation of coalbed can promote methane production e ﬀ ectively. However, CBM thermal recovery needs a forerunner heating time before the apparent enhancement of production. The modeling and simulation results may improve the current cognitions of CBM thermal recovery.


Introduction
Coalbed methane (CBM) is a cleaner and cheaper resource among the fossil fuels [1][2][3]. It is reported that the accumulated reserve of Chinese CBM is about 10 billion cubic meters, and the recoverable resource represents about 47% [4]. In earlier years, CBM is treated as a nerve-wracking hazard in coal mining. However, people's attitudes are changing that CBM becomes an efficiently and environmentally friendly fuel now. However, the production of CBM meets a great challenge of lower reservoir permeability and higher methane sorption capacity in China. A survey research shows that the permeability of Chinese coal ranges from 1 × 10 −4 mD to 1 × 10 −3 mD while it ranges from 0.1 mD to 1 mD of American or Australian coal [5]. As a result, the traditional direct-recovery method of CBM cannot satisfy the demands of effective production. Therefore, researchers and engineers in this field tried some unconventional methods by using a series of manual treatments to the methane reservoirs, such as the fracturing by water and gas, the displacement of adsorption by carbon dioxide and the stimulation of increasing temperature.
Thermal stimulation to deep CBM reservoirs is an effective method to promote gas recovery. Researchers have explored many different basic theories and mining technologies [6][7][8][9][10]. Before the production of methane, water with high temperature of 80 • C is assumed to be injected into a hypothetical CBM reservoir to causes helpful coal-gas interactions [11]. The simulation results indicate results are mainly qualitative, it is necessary to propose a theoretically permeability model that can better connect the experimental rules with CBM thermal production.
CBM thermal recovery donates a multi-physical issue of coal, gas and heat [32,33]. Abed et al. [34] described a thermo-hydro-mechanical (THM) framework that suitable for modeling the behavior of unsaturated soils and rock. Due to the CBM thermal recovery method of microwave heating, Gao et al. [35] illustrated the interactions among temperature field, coalbed compaction, and methane transfer. According to their thermo-hydro-mechanical model developed in a finite element environment, higher stimulation temperature results in larger recovery radius. Fan et al. [36] discussed the competitive sorption of carbon dioxide with methane and its interactions with water under a thermo-hydro-mechanical-chemical condition. A two-phase-flow model was developed for CO 2 enhanced CBM recovery. Xia et al. [37] established a coupled hydro-thermo-mechanical model for the spontaneous combustion of underground coal and quantitatively predicted the spontaneous combustion locations of Dongtan coal mine. These precursor works have enlightening meanings for the study of multi-physical interactions in CBM thermal recovery. However, a fully coupled heat-gas-coal model is still necessary.
Following on our previous establishment of coal permeability model, this paper developed a completely coupled heat-gas-coal model for deep CBM thermal recovery by considering the interactions among coalbed deformation, methane delivery and heat transport. To evaluate these interactions among three physical fields and the production efficiency of CBM recovery, a numerical simulation using the finite element approach was validated. Finally, a series of analysis work based on the modeling and simulation are carried out.

Experimental Program
In this section, the sorption characteristics of coal to pure CH 4 under variable temperature is observed in a self-developed gas adsorption apparatus, see Figure 1. Coal samples that acquired from the exploratory borehole of Jinjia lignite mine China are powdered into size of 150 µm. The experimental temperature varies from 25 • C to 85 • C, concretely 25 • C, 45 • C, 65 • C and 85 • C as representatives. At each temperature, the adsorption content was tested in a pressure vessel at the pressure that increases from the atmospheric pressure to a maximum of 10 MPa.
Processes 2018, 6, x FOR PEER REVIEW 3 of 17 above experimental results are mainly qualitative, it is necessary to propose a theoretically permeability model that can better connect the experimental rules with CBM thermal production. CBM thermal recovery donates a multi-physical issue of coal, gas and heat [32,33]. Abed et al. [34] described a thermo-hydro-mechanical (THM) framework that suitable for modeling the behavior of unsaturated soils and rock. Due to the CBM thermal recovery method of microwave heating, Gao et al. [35] illustrated the interactions among temperature field, coalbed compaction, and methane transfer. According to their thermo-hydro-mechanical model developed in a finite element environment, higher stimulation temperature results in larger recovery radius. Fan et al. [36] discussed the competitive sorption of carbon dioxide with methane and its interactions with water under a thermo-hydro-mechanical-chemical condition. A two-phase-flow model was developed for CO2 enhanced CBM recovery. Xia et al. [37] established a coupled hydro-thermo-mechanical model for the spontaneous combustion of underground coal and quantitatively predicted the spontaneous combustion locations of Dongtan coal mine. These precursor works have enlightening meanings for the study of multi-physical interactions in CBM thermal recovery. However, a fully coupled heatgas-coal model is still necessary.
Following on our previous establishment of coal permeability model, this paper developed a completely coupled heat-gas-coal model for deep CBM thermal recovery by considering the interactions among coalbed deformation, methane delivery and heat transport. To evaluate these interactions among three physical fields and the production efficiency of CBM recovery, a numerical simulation using the finite element approach was validated. Finally, a series of analysis work based on the modeling and simulation are carried out.

Experimental Program
In this section, the sorption characteristics of coal to pure CH4 under variable temperature is observed in a self-developed gas adsorption apparatus, see Figure 1. Coal samples that acquired from the exploratory borehole of Jinjia lignite mine China are powdered into size of 150 μm. The experimental temperature varies from 25 °C to 85 °C, concretely 25 °C, 45 °C, 65 °C and 85 °C as representatives. At each temperature, the adsorption content was tested in a pressure vessel at the pressure that increases from the atmospheric pressure to a maximum of 10 MPa.

Effects of Temperature on Methane Adsorption Content
The directly observed methane adsorption content under variable temperature is represented by the colorful spots in Figure 2. One can obtain that temperature has a significant influence on the

Effects of Temperature on Methane Adsorption Content
The directly observed methane adsorption content under variable temperature is represented by the colorful spots in Figure 2. One can obtain that temperature has a significant influence on the adsorption capacity of Jinjia lignite. For example, when the temperature increases from 25 • C to  A Langmuir equation is often used to describe the methane adsorption content under constant temperature [38][39][40] as: where, sg V is the methane adsorption content, m 3 /kg, L V and b donate the Langmuir volume and pressure constants. However, Equation (1) must be modified if the temperature changes. Here, an additional exponential term is introduced to revise the equation as: where, L P donates the Langmuir pressure constant, Pa. ref T is the reference temperature for methane sorption, K. 1 c and 2 c are the coefficients for pressure and temperature, Pa −1 and K −1 , respectively. Figure 2 also shows the validation results of Equation (2) by the obtained experimental data. The well match of the fitting curves with the experimental results indicates that the modified Langmuir Equation (2) can effectively describe the evolution of methane sorption under variable temperature. This domino effect of gas sorption with changeable temperature indicates useful implication for CBM thermal recovery. A Langmuir equation is often used to describe the methane adsorption content under constant temperature [38][39][40] as:

Volumetric Stain Induced by Methane Desorption
where, V sg is the methane adsorption content, m 3 /kg, V L and b donate the Langmuir volume and pressure constants. However, Equation (1) must be modified if the temperature changes. Here, an additional exponential term is introduced to revise the equation as: where, P L donates the Langmuir pressure constant, Pa. T re f is the reference temperature for methane sorption, K. c 1 and c 2 are the coefficients for pressure and temperature, Pa −1 and K −1 , respectively.

Volumetric Stain Induced by Methane Desorption
Coal matrix expands when it adsorbs gas. According to [41], the gas sorption-induced volumetric strain has a linear relationship with the gas adsorption content as: where α sg is the expansion coefficient. Thus, the volume deformation of coal caused by methane adsorption at different temperature can be expressed as: where ∆T = T − T 0 shows the change of temperature.

Mathematical Model
Thermal recovery of CBM is a coupled process of coalbed deformation, methane delivery and heat transport. Among these three fields, one affects another. For example, the compaction of coalbed block causes the reduction of coal porosity and permeability for gas flow while the decrease of methane pressure changes the effective stress to promote coalbed deformation. Moreover, thermal expansion that caused by temperature change affects both coal deformation and gas flow. Following section is to establish a completely coupled heat-gas-coal model for CBM thermal recovery by considering the interactions among coalbed deformation, methane flow and heat transport. Before any derivation, coalbed is assumed as one kind of homogeneous, isotropic continuum while methane is treated as ideal gas. Methane transport in porous coalbed is treated as a Darcy's flow.

Coalbed Deformation Equation
Based on rock elasticity theory, the deformation of coal reservoir can be governed by a Navier-type equilibrium equation. In this paper, the equilibrium equation is optimized by the variable temperature [42] as: in which F i and p represent the body force of coal and the pressure of methane, MPa. G = 0.5E/(1 + ν) and K = E/(3(1 − 2ν)) are the shear and bulk modulus, respectively. E and ν represent the elasticity modulus and Poisson's ratio, respectively. α = 1 − K/K s , (α < 1) means the Biot's coefficients, where K s donates the modulus of coal grains. The volumetric stain that induced by temperature change is defined [43,44] as: where, α T is the thermal expansion coefficient, K −1 .

Methane Flow Equation
The CBM flow in reservoirs obeys a mass conservation equation [37]: Processes 2019, 7, 526 6 of 16 in which Q sκ represents the methane source. v g is the velocity vector of flow, m/s; m is the content of methane in coalbed that can be expressed as: where, the first and second terms represent the free and adsorbed components, respectively. φ is the porosity of coal, ρ c represents coal density while ρ ga represents the methane density at standard conditions. The real gas density can be expressed as: in which, R donates the universal gas constant and M g means the molar mass.
Darcy's velocity vector v g in Equation (7) is proportional to the pressure gradient ∇p and permeability k: where µ represents the viscosity coefficient of methane. Following on our previous work [31], the permeability of coalbed can be expressed as: in which, K f donates the equivalent modulus of coal fracture. Volume strain ε v that induced by effective stress is expressed as: in which, σ is the mean stress of coal. Substituting Equations (4)-(6) into Equation (12), the permeability model reads as In the experimental testing of coal permeability, Equation (13) can be modified as: where, Figure 3 shows the fitting results of the experimental data [30] by the proposed permeability model. Table 1 lists the fitting parameters of Equation (14). The matching result shows that the proposed coal permeability model can be well used to describe the evolution of coal permeability under variable temperature.  Substituting Equations (8)- (13) into Equation (7), we can obtain the methane flow equation:  Substituting Equations (8)- (13) into Equation (7), we can obtain the methane flow equation: in which,

Heat Transport Equation
The conservation of energy for CBM thermal recovery obeys an equilibrium [45] as: in which, the specific heat capacity of coalbed reads: K eq is the effective thermal conductivity. C g and C c represent the specific heat coefficients for methane and coal. Q T is the heat source.
Combining Equation (18) with Equation (17), one gets the modified conservation equation of energy as: The Equations (5), (16) and (19) make up a completely coupled heat-gas-coal model. The mathematical model is used in simulating CBM thermal production. Figure 4 chooses a modeling and simulation domain to represent a half area of CBM thermal recovery from the research of Shahtalebi et al. [15]. Based on a partial-differential-equation solver of COMSOL with MATLAB that runs on Windows 10 environment with Intel Core i7-6500U and RAM 16.0 GB hardware, the simulator is implemented in 1049 s. In Figure 4, the simulation domain consists of 4208 elements. The computation time is 10 9 s that are derived into 100 steps.

Modeling on CBM thermal Recovery
The length and width of the simulation domain are one hundred meters and forty meters, respectively. The depth of methane recovery well into coalbed is seventy meters. A thermal stimulation well in size of thirty meters locates in the center. For the CBM reservoir, the initial methane pressure and temperature are 3.5 MPa and 298 K, respectively. For coal deformation, the boundary AB and BC are confined to normal force of 15 and 8 MPa, while the boundary AD and CD are constrained by normal displacement. For gas flow, the boundary AB, BC, CD and DE are symmetric boundary with no flow, while AE is common pressure boundary of 0.1 MPa. For heat transfer, the boundary AB, BC, CD and DA are thermal insulation, while FG is common temperature boundary of 373 K. To contrast simulation results, comparison sites, no. 1, 2, 3 and 4, are pre-set. Other simulation parameters are listed in Table 2.  Figure 4 chooses a modeling and simulation domain to represent a half area of CBM thermal recovery from the research of Shahtalebi et al. [15]. Based on a partial-differential-equation solver of COMSOL with MATLAB that runs on Windows 10 environment with Intel Core i7-6500U and RAM 16.0 GB hardware, the simulator is implemented in 1049 s. In Figure 4, the simulation domain consists of 4208 elements. The computation time is 10 9 s that are derived into 100 steps.    From the figures, one can obviously find that the larger distances from thermal stimulation well correspond to lower increment of reservoir temperature. For earlier years, temperature change is only confined to the vicinity of thermal well. It indicates that CBM thermal recovery needs a long period of heating time to get beneficial output. For example, after one year (about 3 × 10 7 s), the reservoir temperature keeps almost it initial temperature, the temperature at point 4 is 320 K while the   From the figures, one can obviously find that the larger distances from thermal stimulation well correspond to lower increment of reservoir temperature. For earlier years, temperature change is only confined to the vicinity of thermal well. It indicates that CBM thermal recovery needs a long period of heating time to get beneficial output. For example, after one year (about 3 × 10 7 s), the reservoir temperature keeps almost it initial temperature, the temperature at point 4 is 320 K while the From the figures, one can obviously find that the larger distances from thermal stimulation well correspond to lower increment of reservoir temperature. For earlier years, temperature change is only confined to the vicinity of thermal well. It indicates that CBM thermal recovery needs a long period of heating time to get beneficial output. For example, after one year (about 3 × 10 7 s), the reservoir temperature keeps almost it initial temperature, the temperature at point 4 is 320 K while the temperatures at point 3, 2 and 1 are lower than 300 K. However, when the heating time is longer than 10 years, obvious increment of temperature is observed in each observation point. After a heating time of 30 years, the reservoir temperature is higher than 350 K.

Evolution of Methane Pressure
The distribution of methane pressure in reservoirs changes a lot during CBM thermal recovery. It can be seen from the contour distribution of methane pressure after 1, 10, 20 and 30 years in Figure 7 that the methane pressure decreases in the first 10 years. From the figure, we can easily find that the mean gas pressure is about 2.1 MPa after 10 years recovery, whereas 0.65, 0.33 and 0.24 MPa after 20, 30 and 40 years.   According to the state equation of ideal gas, higher temperature leads to a higher gas pressure in a confined space. Figure 7 shows that higher coalbed temperature doesn't have prominent domination to methane pressure in reservoirs, even in the area near the thermal stimulation source. This is because the CBM is not confined in coal blocks seriously, but strongly transported out from the narrow fracture network. Figure 8 shows the evolution of coalbed permeability at different comparison sites. From the figure, one can draw two conclusions. The first is that coal permeability enlargers firstly and then decreases a little tiny bit. The second is that the permeability ratio increases with the decreasing distance from thermal stimulation well. Further, the normalized permeability ratio enlargers to the maximal value of about 1.53 firstly and then decreases. Although the evolution trends of normalized permeability ratio at 4 different observation points are similar, the corresponding production time node for the maximum permeability ratio delays with the increasing distance from the heat injection well, circumstantiate 8 × 10 8 , 7 × 10 8 , 5 × 10 8 and 2 × 10 8 s, respectively. From Equation (8), and the numerical results of Figures 6 and 8, one can explain the permeability evolution trend as that the accelerative effect of thermal desorption due to increased temperature plays a dominant role in enlarging coal permeability in the increasing stage, however the effective stress induced coal compaction gradually takes over the dominant role to reduce coal permeability in the decreasing stage where the reservoir temperature is stabilized.

Evolution of Coalbed Permeability
Processes 2018, 6, x FOR PEER REVIEW 12 of 17 According to the state equation of ideal gas, higher temperature leads to a higher gas pressure in a confined space. Figure 7 shows that higher coalbed temperature doesn't have prominent domination to methane pressure in reservoirs, even in the area near the thermal stimulation source. This is because the CBM is not confined in coal blocks seriously, but strongly transported out from the narrow fracture network. Figure 8 shows the evolution of coalbed permeability at different comparison sites. From the figure, one can draw two conclusions. The first is that coal permeability enlargers firstly and then decreases a little tiny bit. The second is that the permeability ratio increases with the decreasing distance from thermal stimulation well. Further, the normalized permeability ratio enlargers to the maximal value of about 1.53 firstly and then decreases. Although the evolution trends of normalized permeability ratio at 4 different observation points are similar, the corresponding production time node for the maximum permeability ratio delays with the increasing distance from the heat injection well, circumstantiate 8 × 10 8 , 7 × 10 8 , 5 × 10 8 and 2 × 10 8 s, respectively. From Equation (8), and the numerical results of Fig. 6 and Figure 8, one can explain the permeability evolution trend as that the accelerative effect of thermal desorption due to increased temperature plays a dominant role in enlarging coal permeability in the increasing stage, however the effective stress induced coal compaction gradually takes over the dominant role to reduce coal permeability in the decreasing stage where the reservoir temperature is stabilized.  Figure 9a and 9b show the cumulative methane production and recovery efficiency in CBM thermal recovery cases with different stimulation temperatures. From the cumulative production, one can conclude that the thermal recovery with higher stimulation temperature has greater promotion to the final production of methane. However, the equal increment of coalbed temperature generates a larger promoting efficiency at lower temperature level due to the great enhancement of temperature to methane desorption. Figure 9b shows the efficiency of the increasing yield of coalbed methane thermal production. From the figure, one can see that the maximum production efficiencies of methane with stimulation temperatures of 50 °C, 75 °C and 100 °C are enlarged by 20%, 30% and 33% respectively compared with the conversional production method. However, Figure 9(b) also indicates that thermal recovery of coalbed methane needs a forerunner heating time to obtain apparent enhancement of production.   Figure 9a,b show the cumulative methane production and recovery efficiency in CBM thermal recovery cases with different stimulation temperatures. From the cumulative production, one can conclude that the thermal recovery with higher stimulation temperature has greater promotion to the final production of methane. However, the equal increment of coalbed temperature generates a larger promoting efficiency at lower temperature level due to the great enhancement of temperature to methane desorption. Figure 9b shows the efficiency of the increasing yield of coalbed methane thermal production. From the figure, one can see that the maximum production efficiencies of methane with stimulation temperatures of 50 • C, 75 • C and 100 • C are enlarged by 20%, 30% and 33% respectively compared with the conversional production method. However, Figure 9b also indicates that thermal recovery of coalbed methane needs a forerunner heating time to obtain apparent enhancement of production. Processes 2018, 6, x FOR PEER REVIEW 13 of 17 (a) Cumulative production.

Methane Production with Different Thermal Stimulation Temperature
(b) production efficiency. Figure 9. Gas production with different stimulation temperatures.

Methane Production with Different Initial Permeability
To evaluate the effects of initial reservoirs permeability on CBM thermal recovery, five numerical cases with different permeability and consistent stimulation temperature of 100 °C are carried out. Evolutions of total methane production are shown in Figure 10a. We can conclude that the methane production enlargers exponentially with production time. To measure the production efficiency, we define an index of output proportion as a percentage of cumulative production to total content of reservoir. Figure 10b shows the evolution of production time with the initial permeability when the output proportion percentage takes value of 25%, 50%, 75% and 90%, respectively. From Figure 10, one can obtain that larger value of initial permeability represents higher methane production in earlier time. For example, it takes 17 and 0.5 years respectively to recover 75% content of the total reserve when the reservoir permeability is 5×10 −19 m 2 and 5×10 −17 m 2 , respectively.

Methane Production with Different Initial Permeability
To evaluate the effects of initial reservoirs permeability on CBM thermal recovery, five numerical cases with different permeability and consistent stimulation temperature of 100 • C are carried out. Evolutions of total methane production are shown in Figure 10a. We can conclude that the methane production enlargers exponentially with production time. To measure the production efficiency, we define an index of output proportion as a percentage of cumulative production to total content of reservoir. Figure 10b shows the evolution of production time with the initial permeability when the output proportion percentage takes value of 25%, 50%, 75% and 90%, respectively. From Figure 10, one can obtain that larger value of initial permeability represents higher methane production in earlier time. For example, it takes 17 and 0.5 years respectively to recover 75% content of the total reserve when the reservoir permeability is 5×10 −19 m 2 and 5×10 −17 m 2 , respectively. Processes 2018, 6, x FOR PEER REVIEW 14 of 17 (a) Cumulative production.
(b) production time for different output proportion Figure 10. Gas production with different initial reservoir permeability.

Conclusions
Methane adsorption and coal permeability decrease with increasing temperature. This paper experimentally validated the evolutions of gas sorption and coal permeability under variable temperature. It also established a completely coupled heat-gas-coal model including multi-physics of coalbed deformation, methane delivery and heat transport. The mathematical model was applied to a computational simulation of CBM thermal recovery and solved by COMSOL with MATLAB in a finite element approach environment. To evaluate the thermal recovery process and the production efficiency, a series of analysis work were carried out. The results show that: (1) Thermal stimulation with higher temperature contributes to greater promotion to CBM recovery production. However, it needs a forerunner heating time before an apparent enhancement of production. (2) The normalized coal permeability ratio increases to the maximal value of about 1.53 firstly, and then decreases slightly. The corresponding production time node for the maximum permeability ratio delays with the

Conclusions
Methane adsorption and coal permeability decrease with increasing temperature. This paper experimentally validated the evolutions of gas sorption and coal permeability under variable temperature. It also established a completely coupled heat-gas-coal model including multi-physics of coalbed deformation, methane delivery and heat transport. The mathematical model was applied to a computational simulation of CBM thermal recovery and solved by COMSOL with MATLAB in a finite element approach environment. To evaluate the thermal recovery process and the production efficiency, a series of analysis work were carried out. The results show that: (1) Thermal stimulation with higher temperature contributes to greater promotion to CBM recovery production. However, it needs a forerunner heating time before an apparent enhancement of production. (2) The normalized coal permeability ratio increases to the maximal value of about 1.53 firstly, and then decreases slightly.