# Nonlinear Diffusion-Based Multifield Coupling Model for Thermally Enhanced CBM Recovery

Authors: Rui Yang; Weiqun Liu; Jiakun Lv
Journal: Geofluids
 (2022)

---

## Abstract

Thermal stimulation is a supplementary technique for enhancing gas recovery from coalbed methane (CBM) reservoirs that has received considerable attention worldwide. Investigating gas and heat transfer in coal seams during thermally enhanced CBM recovery is of great significance for predicting gas production and optimizing the extraction method. Gas diffusion in the coal matrix and coupled multiphysics are two of the most important aspects when analyzing gas migration and heat transport. However, previous studies either neglected the nonlinear diffusion process of gas or only assumed that the diffusion coefficient varies with production time. However, considerable experimental results indicate that the gas diffusion coefficient is not only determined by time but also by temperature, which strongly impacts multifield interaction during gas recovery. Thus, prior diffusion models and coupled models must be modified. In this paper, a time-and-temperature-dependent gas diffusion model is established based on fractal theory and experimental data. The proposed diffusion model is embedded into the coupled thermal-hydro-mechanical model to comprehensively describe the behavior of coal deformation, gas migration, and heat transport during CBM recovery. Additionally, both new diffusion model and coupled model were validated with experimental results or field test data, showing that these developed models are applicable for modeling long-term gas diffusion and gas production. Finally, the coupled model was implemented into COMSOL Multiphysics software, and a series of numerical simulations were conducted. The calculation results showed that heat injection could promote gas desorption and diffusion in the matrix while inhabiting gas flow in fractures near the injection well. Gas dynamic diffusion could inhabit gas migration in both matrix system and fracture system at a later production stage. This also means that ignoring the gas nonlinear diffusion process leads to a severe overestimation of coal permeability and gas production.


---

## Body

## 1. Introduction

Coalbed methane (CBM), as a type of high effective and clean gas resource, plays an increasingly important role in China energy structure [1, 2]. And methane outburst is also one of the most typical coal mine accidents [3, 4]. Thus, CBM recovery can not only provide a secure environment for coal mining but also meet the enormous energy demands of China. However, with complicated geological conditions, Chinese CBM reservoirs are generally in the state of ultralow permeability and high stress [5]. To enhance the CBM production and reservoir permeability, a number of enhancement methods have been investigated such as thermal stimulation techniques, which has gained great attention in past few decades [6]. But the temperature change can further complex the gas transport in coal seam and coal-gas interactions, which have a significant effect on CBM production and coal permeability [7, 8]. Therefore, an investigation on heat and gas transfer behavior during thermally enhanced coalbed methane (ECBM) recovery is urgently needed.The mechanical property of rock is extremely complex [9–15]. And abundant studies have found that temperature has a strong impact on coal permeability and adsorptivity. Ju et al. [16] investigated the influence of temperature on permeability of fractured coals by experiments. The results showed that as the temperature increases, the permeability decreased first and then increased. Wei et al. [17] studied the permeability evolution of gas-bearing coal with different temperature, and from the experimental data, the researchers proposed an equation to obtain the coal permeability under the combined effects of coal deformation and temperature. Levy et al. [18] found a linear decrease in the methane adsorption capacity of 0.12 m3/ton per one-degree increase in the temperature. Additionally, thermally enhanced coalbed methane recovery is a typical coupled thermo-hydro-mechanical process [19]. Zhu et al. [20] established a multifield coupling model to reveal the effects of thermal expansion, thermal convection, and thermal diffusion on gas transport and coal deformation during long-term CBM production. On the basis of their results, the temperature variation cannot be neglected in analyzing the coal-gas interactions. Teng et al. [21] built a more comprehensive model by considering the moisture volatilization as a result of high temperature, and the permeability changes of fracture and matrix induced by thermal stimulation is thoroughly discussed in this case. To further understand the rule of gas and heat transfer in coal seam, they evolved their investigation to a fully coupled heat-coal-gas model [22], which was well validated by both experiment and field data.This brief review above shows that certain scholars have studied the heat and gas transfer with multiphysical coupling for ECBM recovery. However, these models were proposed under the same assumption that gas diffusion is neglected or that the gas diffusion coefficient is a constant. Coal is a porous medium, and it is usually considered to be composed of matrix system and fracture system [23]. In gas recovery, the dual-porosity/single-permeability model is widely used to describe the process. In this model, the adsorbed and free gas in the matrix system firstly diffuses into the fracture system and then further flows into the production wells [24, 25]. Therefore, diffusion is a significant factor for controlling CBM production. But as previously discussed, the gas diffusion was generally seen as a linear process, which is defined by Fick’s law [26]. However, as a matter of fact, multiple experiments have shown that the traditional Fick’s law cannot describe the gas diffusion process over entire timescale [27–29], and many dynamic diffusion models have been derived (see Table 1). Li et al. [30] proposed a nonlinear diffusion model according to experimental results and fractal theory. The model illustrated that the gas diffusion coefficient decreases exponentially with time. Zhao et al. [31] established a time-dependent diffusion model based on force analysis of methane during desorption process, indicating that diffusion coefficient had an initial sudden reduce and a subsequent gentle decrease. Meanwhile, Yue et al. [32] tentatively built a power-law equation to obtain the real diffusion coefficient at any specific time, and the proposed model had been verified by a series of tests. Additionally, considering the complex pore structure of coal, Liu and Lin [33] modeled the both dynamic diffusivity of free gas and adsorbed gas. The investigation also showed that the diffusion coefficient is a strongly time-dependent parameter. Figure 1 presents the fitting results of different diffusion models with experimental data. It is clearly that nonlinear diffusion models can reflect the real diffusion process more accurately than traditional unipore diffusion model. Therefore, embedding nonlinear diffusion model into existing coupling model for CBM extraction is essential. Liu et al. [34] developed a coupled hydromechanical model taking the dynamic diffusion process into account firstly. They found that disregarding the attenuation of diffusion coefficient could overestimate gas production significantly. Then, the nonlinear diffusion model is wildly used in different coupling conditions [35–37], and all of the investigations proved that dynamic diffusion is an indispensable factor in analyzing fluid transfer and predicting gas production during CBM recovery.Table 1
Dynamic diffusion models from different studies.
ReferenceMathematical equationParametersDescriptionLi et al. [30]Dt=D0e−λtD0, λA semiempirical model based on experimental results and fractal characteristics of coalZhao et al. [31]Dt=D0ξtD0, ξA theoretical model on the basis of force analysisYue et al. [32]Dt=D0bt+1D0, bAn empirical power-law equation by analyzing experimental dataLiu and Lin [33]Dt=Db0e−ηt+Ds0e−γt+DrDb0, Ds0, η, γ, DrA theoretical model considering diffusion coefficient of free gas and adsorbed gas, respectively, based on pore structure of coalFigure 1
Fitting results of different diffusion models with the experimental dada from Kang et al. [41] (note: Qt is cumulative gas diffusion amount at time t, and Qinf is total gas diffusion amount).The above-mentioned diffusion models and coupling models are obtained under isothermal condition for primary CBM recovery. However, considerable evidences have shown that temperature have a marked impact on diffusion process [38, 39], which also has a strong feedback on desorption-induced temperature changes [40]. Thus, the previous dynamic diffusion models are inapplicable for thermally enhanced CBM recovery. For all we know, there is no study revealing the complex interactions between heat transfer and dynamic diffusion for CBM or shale gas production. Therefore, developing a coupled thermo-hydro-mechanical model considering the time-and-temperature-dependent diffusion coefficient and analyzing gas and heat transfer coupled the effect of nonlinear diffusion during thermally enhanced CBM recovery are essential and helpful.In this paper, a modified dynamic diffusion model is proposed based on previous studies and laboratory test results. In this model, the diffusion coefficient varies with gas drainage time, gas pressure, and temperature. Then, the new diffusion model is embedded into a coupled thermo-hydro-mechanical model for thermally enhanced CBM recovery. Both the new diffusion model and new coupled model are validated with experimental data or field test data in this study. Subsequently, COMSOL multiphysics software is adopted to numerically analyze the gas and heat transfer during long-term CBM recovery with heat injection. The impacts of extraction methods and diffusion model on spatial and temporal evolution of gas adsorption amount, temperature, and permeability and are also thoroughly discussed. Finally, the residual gas content in different systems (matrix system and fracture system) after twenty-year production is estimated to further analyze the mechanism of gas and heat transfer.Dt is dynamic diffusion coefficient, D0 is initial diffusion coefficient, t is diffusion time, λ is attenuation coefficient, ξ is a parameter related to the time when half of the total adsorbed gas desorbs, b is fitting coefficient, Db0 and Ds0 are initial diffusion coefficient of free gas and adsorbed gas, respectively, η and γ are attenuation coefficient of diffusivity of free and adsorbed gas, respectively, and Dr is residual diffusivity.


## 2. Modeling

In this section, we introduced the thermal governing equation and thermal effect on gas diffusion into the existing governing equation for DPSP geological model and multistage pore model to further clarify the mechanisms of gas and heat transfer during long-term thermal ECBM recovery.


### 2.1. Nonlinear Diffusion Model

Gas diffusion is regarded as the mass exchange between coal matrix and fracture in dual-porosity single-permeability model for CBM recovery. In general, the diffusion amount per unit timeQd obeys the Fick’s law, which reads:
(1)Qd=DχMgRTpm−pf,where D is diffusion coefficient, Mg is molar mass of gas, R is gas molar constant, pm and pf are gas pressure in matrix and fracture, respectively, and χ is the shape factor.Based on traditional unipore model, previous studies usually assumed that diffusion coefficient is a constant during long-time gas production (see on Figure2(a)). But in recent years, the multistage pore structure of coal has been investigated, which is shown in Figure 2(b). It can be seen that the gas diffusion path is consisted of series-connected multiple different-sized pores in multistage pore model. At the initial time (t=t0), the gas diffuses from the larger pore of matrix system into fracture system, when the diffusion coefficient is D0. But with the drainage time increases, the gas diffusivity decreases for longer diffusion path, smaller pore size, and greater diffusion resistant. And from Liu et al.’s study [34], the following equation can describe the gas dynamic diffusion coefficient under multistage pore model:
(2)D=D0exp−λt,where λ is the attenuation coefficient of the dynamic diffusion coefficient, t is diffusion time, and D0 is initial diffusion coefficient, which can be defined as follows according to Arrhenius equation [42]:
(3)D0=D∞exp−qRT,where D∞ is the value of D0 at infinite temperature and q represents the activation energy.Figure 2
Gas diffusion in (a) unipore model and (b) multistage pore model during long-term gas production.In addition, the values ofD∞ and q change with gas pressure based on Li et al.’s experimental result [42]. Figure 3 depicts the relationship between initial diffusion coefficient at infinite temperature D∞ and activation energy q and gas pressure. It can be clearly observed that the values of D∞ and q decrease with the increase of gas pressure. In this investigation, we fitted the laboratory data of D∞ with an exponential equation, which is expressed as:
(4)D∞=D∞aexp−Apm−pa,where D∞a is the value of D∞ at atmospheric pressure, pm is gas pressure of coal matrix, and pa is atmospheric pressure. For activation energy, the laboratory data is fitted to a linear relation written as:
(5)q=qa−Bpm−pa,where qa is the value of q at atmospheric pressure; the parameters A and B are coefficients determined by fitting curves. A larger values of A and B would result in more decrease of D∞ and q with gas pressure, respectively. By substituting Equations (3)–(5) into Equation (2), we obtain the new nonlinear diffusion model:
(6)D=D∞aexp−Apm−pa−qa−Bpm−paRT−λt,where D is dynamic diffusion coefficient. In this case, the value of D is not a constant, but varies with temperature, gas pressure, and diffusion time.Figure 3
Relationship between initial diffusion coefficient at infinite temperature (D∞) and activation energy q and gas pressure.


### 2.2. Coupled Thermal-Hydro-Mechanical Model

In the following, a series of governing equations are established to further clarify the gas and heat transfer during thermally enhanced coalbed methane recovery, which includes coal deformation equation, gas dynamic diffusion equation, gas seepage equation, and heat transfer equation. The new model is derived based on the following hypotheses:(1)
Coal is homogeneous and isotropic(2)
The plasticity of coal seam is not considered(3)
The dual-porosity single-permeability geological model (DPSP) is adopted for gas flow(4)
The porosity of coal matrix is invariant(5)
The water distribution in coal is disregarded(6)
The chemical reaction during thermal ECBM recovery is neglected


#### 2.2.1. Coal Deformation

Based on the constitutive relations of poroelasticity, the governing equation for coal deformation that considers desorption-induced strain and temperature-induced strain can be expressed as follows:(7)Gui,jj+G1−2νuj,ji−αmpm,i−αfpf,i−KαTT,i−Kεs,i+fi=0,where G is the shear modulus, ν is the Poisson’s ratio, K is the bulk modulus, αT is the thermal expansion coefficient, ui and fi are displacement and force in i-direction, respectively, αm and αf are Biot effective stress coefficients, T is temperature, pf is fracture pressure, and εs is desorption-induced strain, which can be calculated by a Langmuir-type equation:
(8)εs=εLpmpm+pLexp−b2T−T01+b1pm,where εL and pL represent Langmuir strain constant and Langmuir pressure constant, respectively, b1 and b2 are the pressure and temperature coefficients, respectively, and T0 is reference temperature.


#### 2.2.2. Gas Diffusion in Matrix

In DPSP model, the mass conservation law of coal matrix system can be described as the following:(9)∂mm∂t=−Q1−φf,where mm represents the methane content per unit volume of the matrix, φf is the porosity of fracture, and Q is the mass exchange between coal matrix and fracture, which can be defined as:
(10)Q=DχMgRTpm−pf,where χ is the shape factor and Mg is molar mass of gas. In general, the diffusion coefficient D is regarded as a constant, but in this study, we embed the new proposed nonlinear diffusion model into Equation (10), which is dependent on diffusion time, gas pressure, and temperature. In addition, the total methane content in coal matrix consists of two components: free gas and adsorbed gas. Thus the value of mm can be given as follows:
(11)mm=φmρ+ρaρcVsg1−φm−φf,where φm is the porosity of matrix, ρ is gas density, ρa is the gas density under standard condition, ρc is coal density, and Vsg represents the adsorbed gas content in coal matrix, which is defined by Langmuir volume equation:
(12)Vsg=VLpmpm+pLexp−b2T−T01+b1pm,where VL is Langmuir volume constant. By substituting Equations (10)–(12) and Equation (6) into Equation (9), we yield the governing equation for gas diffusion:
(13)ρaρc∂Vsg∂T−φmMgRT2pm∂T∂t+ρaρc∂Vsg∂pm+φmMRT∂pm∂t=MgχD∞aexp−Apm−pa−qa−Bpm−pa/RT−λtpf−pmRT.


#### 2.2.3. Gas Flow in Fracture

Gas transport in fracture obeys the Darcy’s law and mass conservation law, so the governing equation for gas seepage in fracture is expressed as:(14)∂φfρ∂t+∇∙−kμρ∇pf=Q,where k is permeability, μ is velocity of gas, and Q is the source or sink of gas. In DPSP model, the value of Q equals to diffusion amount of gas in matrix. Additionally, the coal deformation has a significant impact on fracture porosity and permeability. Based on our pervious study [43] and cubic law, the porosity and permeability varying with mean effective stress are shown as Equation (15) and Equation (16), respectively:
(15)φf=αf+φf0−αfexp−Δσ′K,(16)k=k0αfφf0+φf0−αfφf0exp−Δσ′K3,where σ′ is effective stress and φf0 and k0 are initial porosity of fracture and permeability, respectively.Substituting Equation (6), Equation (10), Equation (15), and Equation (16) into Equation (14), we obtain:
(17)−ρSK∂σ′∂t+MgφfRT∂pf∂t−ρφfT∂T∂t+∇∙−kμρ∇pf=MgχD∞aexp−Apm−pa−qa−Bpm−pa/RT−λtpm−pfRT,where
(18)S=φf−αfexp−Δσ′K.


#### 2.2.4. Heat Transfer

Considering thermal dilatation of gas and coal, thermal convection, thermal conduction, and gas adsorption energy, the governing equation for heat transfer can be described as:(19)∂ρCMT∂t⏟Energychange+TKgαg∇∙−kμ∇pf⏟Thermaldilatationofgas+TKαT∂εv∂t⏟Thermaldilatationofcoal+qstρcρaMg∂Vsg∂t⏟Gasadsorptionenergy=−∇∙−1−φf−φmλs+φf+φmλg∇T⏟Thermalconduction−ρCgkμ∇pfT⏟Thermalconvection,where Kg is bulk modulus of gas, αgαg=1/T is thermal expansion coefficient of gas, εv is volumetric strain, qst is the isosteric heat of adsorption, Cg is gas specific heat constant, λs and λg are thermal conductivities of coal and gas, respectively, and ρCM represents specific heat capacity of gas-filled coal, which can be expressed as:
(20)ρCM=φf+φmρCg+1−φf−φmρcCs,where Cs is coal specific heat constant.Assuming1−φf−φm≈1, 1−φf−φmλs≫φf+φmλg and Kg=pf [44], Equation (19) evolves:
(21)ρCM+ρcρaqstM∂Vsg∂T∂T∂t+ρcρaqstM∂Vsg∂pm∂pm∂t+TKαT∂εv∂t+pf∇∙−kμ∇pf=λs∇2T+ρCgkμ∇pf∇T.


#### 2.2.5. Cross Coupling

Equations (7), (13), (17), and (21) are four governing equations for thermally enhanced coalbed methane recovery. The cross-couplings among these physical fields are shown in Figure 4. Although the interactions among coal deformation, gas flow in fracture, and heat transfer have been comprehensively investigated, the impacts of nonlinear gas diffusion on the three mentioned physical fields and the corresponding feedback mechanisms are introduced for the first time in this study (see on Figure 4), which involve five aspects:
(a)
The impact of gas nonlinear diffusion on coal deformation is reflected by the third term and sixth term on the left side of Equation (7). They refer to the variation of gas pressure in matrix and adsorption-induced strain, respectively. In particular, both of them have strong correlation with the dynamic diffusion coefficient(b)
The impact of gas nonlinear diffusion on gas flow in fracture is reflected by the right side of Equation (17). It refers to the mass exchange of gas between coal matrix and fracture. In particular, it not only depends on the pressure difference between matrix and fracture but also diffusion time and temperature(c)
The feedback mechanism of gas flow in fracture to gas nonlinear diffusion is reflected by the right side of Equation (13). It equal and opposite with the right side of Equation (17). That ensures the mass of gas within a representative elementary volume is conserved(d)
The impact of gas nonlinear diffusion on heat transfer is reflected by the first and second term on the left side of Equation (21). They refer to the energy change induced by gas adsorption. In particular, the decrease rate of gas adsorption amount obviously decreases with prolonged production period(e)
The feedback mechanism of heat transfer to gas nonlinear diffusion is reflected by the first term on the left side and the right side of Equation (13). The former refers to the change of gas adsorption amount induced by temperature, and the latter refers to the time-and-temperature-dependent diffusion coefficientFigure 4
Cross-couplings among multiple physical fields.From aforementioned discussion, there are four variablespm,pf,T,εvorσ′ for the coupled mathematic model. All the variables are interactive, and all the governing equations are nonlinear partial differential equations (PDE). To solve the complex problem, COMSOL Multiphysics software is adopted which provides a powerful PDE-based modeling environment. In this study, the coal deformation equation is calculated by the solid mechanics module, and three general PDE modules are utilized to address gas dynamic diffusion, gas flow in fracture, and heat transfer, respectively.


## 2.1. Nonlinear Diffusion Model

Gas diffusion is regarded as the mass exchange between coal matrix and fracture in dual-porosity single-permeability model for CBM recovery. In general, the diffusion amount per unit timeQd obeys the Fick’s law, which reads:
(1)Qd=DχMgRTpm−pf,where D is diffusion coefficient, Mg is molar mass of gas, R is gas molar constant, pm and pf are gas pressure in matrix and fracture, respectively, and χ is the shape factor.Based on traditional unipore model, previous studies usually assumed that diffusion coefficient is a constant during long-time gas production (see on Figure2(a)). But in recent years, the multistage pore structure of coal has been investigated, which is shown in Figure 2(b). It can be seen that the gas diffusion path is consisted of series-connected multiple different-sized pores in multistage pore model. At the initial time (t=t0), the gas diffuses from the larger pore of matrix system into fracture system, when the diffusion coefficient is D0. But with the drainage time increases, the gas diffusivity decreases for longer diffusion path, smaller pore size, and greater diffusion resistant. And from Liu et al.’s study [34], the following equation can describe the gas dynamic diffusion coefficient under multistage pore model:
(2)D=D0exp−λt,where λ is the attenuation coefficient of the dynamic diffusion coefficient, t is diffusion time, and D0 is initial diffusion coefficient, which can be defined as follows according to Arrhenius equation [42]:
(3)D0=D∞exp−qRT,where D∞ is the value of D0 at infinite temperature and q represents the activation energy.Figure 2
Gas diffusion in (a) unipore model and (b) multistage pore model during long-term gas production.In addition, the values ofD∞ and q change with gas pressure based on Li et al.’s experimental result [42]. Figure 3 depicts the relationship between initial diffusion coefficient at infinite temperature D∞ and activation energy q and gas pressure. It can be clearly observed that the values of D∞ and q decrease with the increase of gas pressure. In this investigation, we fitted the laboratory data of D∞ with an exponential equation, which is expressed as:
(4)D∞=D∞aexp−Apm−pa,where D∞a is the value of D∞ at atmospheric pressure, pm is gas pressure of coal matrix, and pa is atmospheric pressure. For activation energy, the laboratory data is fitted to a linear relation written as:
(5)q=qa−Bpm−pa,where qa is the value of q at atmospheric pressure; the parameters A and B are coefficients determined by fitting curves. A larger values of A and B would result in more decrease of D∞ and q with gas pressure, respectively. By substituting Equations (3)–(5) into Equation (2), we obtain the new nonlinear diffusion model:
(6)D=D∞aexp−Apm−pa−qa−Bpm−paRT−λt,where D is dynamic diffusion coefficient. In this case, the value of D is not a constant, but varies with temperature, gas pressure, and diffusion time.Figure 3
Relationship between initial diffusion coefficient at infinite temperature (D∞) and activation energy q and gas pressure.


## 2.2. Coupled Thermal-Hydro-Mechanical Model

In the following, a series of governing equations are established to further clarify the gas and heat transfer during thermally enhanced coalbed methane recovery, which includes coal deformation equation, gas dynamic diffusion equation, gas seepage equation, and heat transfer equation. The new model is derived based on the following hypotheses:(1)
Coal is homogeneous and isotropic(2)
The plasticity of coal seam is not considered(3)
The dual-porosity single-permeability geological model (DPSP) is adopted for gas flow(4)
The porosity of coal matrix is invariant(5)
The water distribution in coal is disregarded(6)
The chemical reaction during thermal ECBM recovery is neglected


### 2.2.1. Coal Deformation

Based on the constitutive relations of poroelasticity, the governing equation for coal deformation that considers desorption-induced strain and temperature-induced strain can be expressed as follows:(7)Gui,jj+G1−2νuj,ji−αmpm,i−αfpf,i−KαTT,i−Kεs,i+fi=0,where G is the shear modulus, ν is the Poisson’s ratio, K is the bulk modulus, αT is the thermal expansion coefficient, ui and fi are displacement and force in i-direction, respectively, αm and αf are Biot effective stress coefficients, T is temperature, pf is fracture pressure, and εs is desorption-induced strain, which can be calculated by a Langmuir-type equation:
(8)εs=εLpmpm+pLexp−b2T−T01+b1pm,where εL and pL represent Langmuir strain constant and Langmuir pressure constant, respectively, b1 and b2 are the pressure and temperature coefficients, respectively, and T0 is reference temperature.


### 2.2.2. Gas Diffusion in Matrix

In DPSP model, the mass conservation law of coal matrix system can be described as the following:(9)∂mm∂t=−Q1−φf,where mm represents the methane content per unit volume of the matrix, φf is the porosity of fracture, and Q is the mass exchange between coal matrix and fracture, which can be defined as:
(10)Q=DχMgRTpm−pf,where χ is the shape factor and Mg is molar mass of gas. In general, the diffusion coefficient D is regarded as a constant, but in this study, we embed the new proposed nonlinear diffusion model into Equation (10), which is dependent on diffusion time, gas pressure, and temperature. In addition, the total methane content in coal matrix consists of two components: free gas and adsorbed gas. Thus the value of mm can be given as follows:
(11)mm=φmρ+ρaρcVsg1−φm−φf,where φm is the porosity of matrix, ρ is gas density, ρa is the gas density under standard condition, ρc is coal density, and Vsg represents the adsorbed gas content in coal matrix, which is defined by Langmuir volume equation:
(12)Vsg=VLpmpm+pLexp−b2T−T01+b1pm,where VL is Langmuir volume constant. By substituting Equations (10)–(12) and Equation (6) into Equation (9), we yield the governing equation for gas diffusion:
(13)ρaρc∂Vsg∂T−φmMgRT2pm∂T∂t+ρaρc∂Vsg∂pm+φmMRT∂pm∂t=MgχD∞aexp−Apm−pa−qa−Bpm−pa/RT−λtpf−pmRT.


### 2.2.3. Gas Flow in Fracture

Gas transport in fracture obeys the Darcy’s law and mass conservation law, so the governing equation for gas seepage in fracture is expressed as:(14)∂φfρ∂t+∇∙−kμρ∇pf=Q,where k is permeability, μ is velocity of gas, and Q is the source or sink of gas. In DPSP model, the value of Q equals to diffusion amount of gas in matrix. Additionally, the coal deformation has a significant impact on fracture porosity and permeability. Based on our pervious study [43] and cubic law, the porosity and permeability varying with mean effective stress are shown as Equation (15) and Equation (16), respectively:
(15)φf=αf+φf0−αfexp−Δσ′K,(16)k=k0αfφf0+φf0−αfφf0exp−Δσ′K3,where σ′ is effective stress and φf0 and k0 are initial porosity of fracture and permeability, respectively.Substituting Equation (6), Equation (10), Equation (15), and Equation (16) into Equation (14), we obtain:
(17)−ρSK∂σ′∂t+MgφfRT∂pf∂t−ρφfT∂T∂t+∇∙−kμρ∇pf=MgχD∞aexp−Apm−pa−qa−Bpm−pa/RT−λtpm−pfRT,where
(18)S=φf−αfexp−Δσ′K.


### 2.2.4. Heat Transfer

Considering thermal dilatation of gas and coal, thermal convection, thermal conduction, and gas adsorption energy, the governing equation for heat transfer can be described as:(19)∂ρCMT∂t⏟Energychange+TKgαg∇∙−kμ∇pf⏟Thermaldilatationofgas+TKαT∂εv∂t⏟Thermaldilatationofcoal+qstρcρaMg∂Vsg∂t⏟Gasadsorptionenergy=−∇∙−1−φf−φmλs+φf+φmλg∇T⏟Thermalconduction−ρCgkμ∇pfT⏟Thermalconvection,where Kg is bulk modulus of gas, αgαg=1/T is thermal expansion coefficient of gas, εv is volumetric strain, qst is the isosteric heat of adsorption, Cg is gas specific heat constant, λs and λg are thermal conductivities of coal and gas, respectively, and ρCM represents specific heat capacity of gas-filled coal, which can be expressed as:
(20)ρCM=φf+φmρCg+1−φf−φmρcCs,where Cs is coal specific heat constant.Assuming1−φf−φm≈1, 1−φf−φmλs≫φf+φmλg and Kg=pf [44], Equation (19) evolves:
(21)ρCM+ρcρaqstM∂Vsg∂T∂T∂t+ρcρaqstM∂Vsg∂pm∂pm∂t+TKαT∂εv∂t+pf∇∙−kμ∇pf=λs∇2T+ρCgkμ∇pf∇T.


### 2.2.5. Cross Coupling

Equations (7), (13), (17), and (21) are four governing equations for thermally enhanced coalbed methane recovery. The cross-couplings among these physical fields are shown in Figure 4. Although the interactions among coal deformation, gas flow in fracture, and heat transfer have been comprehensively investigated, the impacts of nonlinear gas diffusion on the three mentioned physical fields and the corresponding feedback mechanisms are introduced for the first time in this study (see on Figure 4), which involve five aspects:
(a)
The impact of gas nonlinear diffusion on coal deformation is reflected by the third term and sixth term on the left side of Equation (7). They refer to the variation of gas pressure in matrix and adsorption-induced strain, respectively. In particular, both of them have strong correlation with the dynamic diffusion coefficient(b)
The impact of gas nonlinear diffusion on gas flow in fracture is reflected by the right side of Equation (17). It refers to the mass exchange of gas between coal matrix and fracture. In particular, it not only depends on the pressure difference between matrix and fracture but also diffusion time and temperature(c)
The feedback mechanism of gas flow in fracture to gas nonlinear diffusion is reflected by the right side of Equation (13). It equal and opposite with the right side of Equation (17). That ensures the mass of gas within a representative elementary volume is conserved(d)
The impact of gas nonlinear diffusion on heat transfer is reflected by the first and second term on the left side of Equation (21). They refer to the energy change induced by gas adsorption. In particular, the decrease rate of gas adsorption amount obviously decreases with prolonged production period(e)
The feedback mechanism of heat transfer to gas nonlinear diffusion is reflected by the first term on the left side and the right side of Equation (13). The former refers to the change of gas adsorption amount induced by temperature, and the latter refers to the time-and-temperature-dependent diffusion coefficientFigure 4
Cross-couplings among multiple physical fields.From aforementioned discussion, there are four variablespm,pf,T,εvorσ′ for the coupled mathematic model. All the variables are interactive, and all the governing equations are nonlinear partial differential equations (PDE). To solve the complex problem, COMSOL Multiphysics software is adopted which provides a powerful PDE-based modeling environment. In this study, the coal deformation equation is calculated by the solid mechanics module, and three general PDE modules are utilized to address gas dynamic diffusion, gas flow in fracture, and heat transfer, respectively.


## 2.2.1. Coal Deformation

Based on the constitutive relations of poroelasticity, the governing equation for coal deformation that considers desorption-induced strain and temperature-induced strain can be expressed as follows:(7)Gui,jj+G1−2νuj,ji−αmpm,i−αfpf,i−KαTT,i−Kεs,i+fi=0,where G is the shear modulus, ν is the Poisson’s ratio, K is the bulk modulus, αT is the thermal expansion coefficient, ui and fi are displacement and force in i-direction, respectively, αm and αf are Biot effective stress coefficients, T is temperature, pf is fracture pressure, and εs is desorption-induced strain, which can be calculated by a Langmuir-type equation:
(8)εs=εLpmpm+pLexp−b2T−T01+b1pm,where εL and pL represent Langmuir strain constant and Langmuir pressure constant, respectively, b1 and b2 are the pressure and temperature coefficients, respectively, and T0 is reference temperature.


## 2.2.2. Gas Diffusion in Matrix

In DPSP model, the mass conservation law of coal matrix system can be described as the following:(9)∂mm∂t=−Q1−φf,where mm represents the methane content per unit volume of the matrix, φf is the porosity of fracture, and Q is the mass exchange between coal matrix and fracture, which can be defined as:
(10)Q=DχMgRTpm−pf,where χ is the shape factor and Mg is molar mass of gas. In general, the diffusion coefficient D is regarded as a constant, but in this study, we embed the new proposed nonlinear diffusion model into Equation (10), which is dependent on diffusion time, gas pressure, and temperature. In addition, the total methane content in coal matrix consists of two components: free gas and adsorbed gas. Thus the value of mm can be given as follows:
(11)mm=φmρ+ρaρcVsg1−φm−φf,where φm is the porosity of matrix, ρ is gas density, ρa is the gas density under standard condition, ρc is coal density, and Vsg represents the adsorbed gas content in coal matrix, which is defined by Langmuir volume equation:
(12)Vsg=VLpmpm+pLexp−b2T−T01+b1pm,where VL is Langmuir volume constant. By substituting Equations (10)–(12) and Equation (6) into Equation (9), we yield the governing equation for gas diffusion:
(13)ρaρc∂Vsg∂T−φmMgRT2pm∂T∂t+ρaρc∂Vsg∂pm+φmMRT∂pm∂t=MgχD∞aexp−Apm−pa−qa−Bpm−pa/RT−λtpf−pmRT.


## 2.2.3. Gas Flow in Fracture

Gas transport in fracture obeys the Darcy’s law and mass conservation law, so the governing equation for gas seepage in fracture is expressed as:(14)∂φfρ∂t+∇∙−kμρ∇pf=Q,where k is permeability, μ is velocity of gas, and Q is the source or sink of gas. In DPSP model, the value of Q equals to diffusion amount of gas in matrix. Additionally, the coal deformation has a significant impact on fracture porosity and permeability. Based on our pervious study [43] and cubic law, the porosity and permeability varying with mean effective stress are shown as Equation (15) and Equation (16), respectively:
(15)φf=αf+φf0−αfexp−Δσ′K,(16)k=k0αfφf0+φf0−αfφf0exp−Δσ′K3,where σ′ is effective stress and φf0 and k0 are initial porosity of fracture and permeability, respectively.Substituting Equation (6), Equation (10), Equation (15), and Equation (16) into Equation (14), we obtain:
(17)−ρSK∂σ′∂t+MgφfRT∂pf∂t−ρφfT∂T∂t+∇∙−kμρ∇pf=MgχD∞aexp−Apm−pa−qa−Bpm−pa/RT−λtpm−pfRT,where
(18)S=φf−αfexp−Δσ′K.


## 2.2.4. Heat Transfer

Considering thermal dilatation of gas and coal, thermal convection, thermal conduction, and gas adsorption energy, the governing equation for heat transfer can be described as:(19)∂ρCMT∂t⏟Energychange+TKgαg∇∙−kμ∇pf⏟Thermaldilatationofgas+TKαT∂εv∂t⏟Thermaldilatationofcoal+qstρcρaMg∂Vsg∂t⏟Gasadsorptionenergy=−∇∙−1−φf−φmλs+φf+φmλg∇T⏟Thermalconduction−ρCgkμ∇pfT⏟Thermalconvection,where Kg is bulk modulus of gas, αgαg=1/T is thermal expansion coefficient of gas, εv is volumetric strain, qst is the isosteric heat of adsorption, Cg is gas specific heat constant, λs and λg are thermal conductivities of coal and gas, respectively, and ρCM represents specific heat capacity of gas-filled coal, which can be expressed as:
(20)ρCM=φf+φmρCg+1−φf−φmρcCs,where Cs is coal specific heat constant.Assuming1−φf−φm≈1, 1−φf−φmλs≫φf+φmλg and Kg=pf [44], Equation (19) evolves:
(21)ρCM+ρcρaqstM∂Vsg∂T∂T∂t+ρcρaqstM∂Vsg∂pm∂pm∂t+TKαT∂εv∂t+pf∇∙−kμ∇pf=λs∇2T+ρCgkμ∇pf∇T.


## 2.2.5. Cross Coupling

Equations (7), (13), (17), and (21) are four governing equations for thermally enhanced coalbed methane recovery. The cross-couplings among these physical fields are shown in Figure 4. Although the interactions among coal deformation, gas flow in fracture, and heat transfer have been comprehensively investigated, the impacts of nonlinear gas diffusion on the three mentioned physical fields and the corresponding feedback mechanisms are introduced for the first time in this study (see on Figure 4), which involve five aspects:
(a)
The impact of gas nonlinear diffusion on coal deformation is reflected by the third term and sixth term on the left side of Equation (7). They refer to the variation of gas pressure in matrix and adsorption-induced strain, respectively. In particular, both of them have strong correlation with the dynamic diffusion coefficient(b)
The impact of gas nonlinear diffusion on gas flow in fracture is reflected by the right side of Equation (17). It refers to the mass exchange of gas between coal matrix and fracture. In particular, it not only depends on the pressure difference between matrix and fracture but also diffusion time and temperature(c)
The feedback mechanism of gas flow in fracture to gas nonlinear diffusion is reflected by the right side of Equation (13). It equal and opposite with the right side of Equation (17). That ensures the mass of gas within a representative elementary volume is conserved(d)
The impact of gas nonlinear diffusion on heat transfer is reflected by the first and second term on the left side of Equation (21). They refer to the energy change induced by gas adsorption. In particular, the decrease rate of gas adsorption amount obviously decreases with prolonged production period(e)
The feedback mechanism of heat transfer to gas nonlinear diffusion is reflected by the first term on the left side and the right side of Equation (13). The former refers to the change of gas adsorption amount induced by temperature, and the latter refers to the time-and-temperature-dependent diffusion coefficientFigure 4
Cross-couplings among multiple physical fields.From aforementioned discussion, there are four variablespm,pf,T,εvorσ′ for the coupled mathematic model. All the variables are interactive, and all the governing equations are nonlinear partial differential equations (PDE). To solve the complex problem, COMSOL Multiphysics software is adopted which provides a powerful PDE-based modeling environment. In this study, the coal deformation equation is calculated by the solid mechanics module, and three general PDE modules are utilized to address gas dynamic diffusion, gas flow in fracture, and heat transfer, respectively.


## 3. Model Validation

In order to ensure the reliability of the proposed diffusion model and coupled model, both of them are verified by experimental results or field test data. In the verification of nonlinear diffusion model, two groups of accurate laboratorial data derived by Ma et al. [45] and Charriere et al. [46] are used to compare with the fitting curves. Additionally, field history data of methane is adopted to achieve the calibration of the coupled model and its implementation into COMSOL Multiphysics.


### 3.1. Validation of the Nonlinear Diffusion Model

Ma et al. [45] conducted a series of experiments to investigate the gas diffusion behavior under different conditions in coal. Coal samples used in the paper originates from the Lorraine basin, France. These samples were ground mechanically into coal particles with size of 0.5-1.0 mm. The experiment was carried out adopting gravimetric sorption system (Rubotherm, Bochum, Germany). During the experiment, the whole diffusion processes of pure CH4 were measured at the temperature of 298.15 K (25°C) and at pressures of 0.1 MPa and 5.0 MPa, respectively. Based on our proposed nonlinear diffusion model, the fitting results between theoretical value of diffusion coefficient and experimental data from Ma et al. [45] are shown in Figure 5(a).Figure 5
Fitting results of our diffusion model with the experimental data from (a) Ma et al. [45] and (b) Charriere et al. [46].
(a)(b)Furthermore, Charriere et al. [46] also analyzed the temperature-dependent gas diffusion process by experiments. The coal samples originating from Zhaogu coal mine (China) were selected for diffusion test in the paper. Before the experiment, the authors sieved the coal particles using sieves with pore size of 0.25-0.5 mm firstly and then dehydrated the prepared particles to ensure the completely dry state. In experimental process, the gas diffusion rates were obtained by a self-assembled sorption-diffusion system. To fully verify our model, we select the experimental results of diffusion coefficient with different temperatures (30°C and 80°C) at constant pressure (0.5 MPa) to compare with the fitting curves, which is shown in Figure 5(b).From Figure5, it is clear that the theoretical results of diffusion coefficient calculated by our model match the experimental data accurately. And this figure also shows that the diffusion coefficient decreases with diffusion time, and this tendency is more obviously in the early stage of gas diffusion. The conditions of high gas pressure and high temperature can promote the diffusivity of gas significantly. All of the above demonstrate that the proposed nonlinear diffusion model can well describe the gas diffusion process in various kinds of coal under different gas pressure and different temperature. In addition, as shown in Figure 5, for gas diffusion in coal originating from the same producer, the values of parameters in Equation (6) are only slightly different under various diffusion conditions, but for gas diffusion in different types of coal, these values are quite different, indicating that the above parameters are dominated by coal category.


### 3.2. Validation of the New Coupled Model

In this section, a field test data of direct recovery [47] is utilized to compare with the numerical simulation result based on the coupled model. The coalbed methane gas field is located in the south of Shanxi Province, China, and the data originates from October 4, 2007, to September 3, 2008. The site situation, initial state, and boundary conditions of the gas field are shown in Figure 6, and the key parameters used for simulation are listed in Table 2. All the listed parameters are form real CBM field.Figure 6
The site situation, initial state, and boundary conditions of the gas field utilized for validation.Table 2
The key parameters used for the simulation in Section3 are listed as follows.
ParameterValueUnitParameterValueUnitDensity of coal1650kg m-3Young’s modulus3GPaPoisson’s ratio0.3—Initial porosity of fracture0.0302—Permeability of coal1.63×10−13m2Langmuir strain constant0.03—Langmuir pressure constant2.7×106PaLangmuir volume constant0.045—Figure7 depicts the actual data and simulation result of gas production rate. It is obviously that the predictive value of gas production rate obtained by the proposed coupled model has an appropriate agreement with the actual values. Both Figures 5 and 7 indicate that our models are applicable in modeling the gas and heat transfer during long-term direct and thermally enhanced CBM recovery.Figure 7
Comparison of our model and the field test data.


## 3.1. Validation of the Nonlinear Diffusion Model

Ma et al. [45] conducted a series of experiments to investigate the gas diffusion behavior under different conditions in coal. Coal samples used in the paper originates from the Lorraine basin, France. These samples were ground mechanically into coal particles with size of 0.5-1.0 mm. The experiment was carried out adopting gravimetric sorption system (Rubotherm, Bochum, Germany). During the experiment, the whole diffusion processes of pure CH4 were measured at the temperature of 298.15 K (25°C) and at pressures of 0.1 MPa and 5.0 MPa, respectively. Based on our proposed nonlinear diffusion model, the fitting results between theoretical value of diffusion coefficient and experimental data from Ma et al. [45] are shown in Figure 5(a).Figure 5
Fitting results of our diffusion model with the experimental data from (a) Ma et al. [45] and (b) Charriere et al. [46].
(a)(b)Furthermore, Charriere et al. [46] also analyzed the temperature-dependent gas diffusion process by experiments. The coal samples originating from Zhaogu coal mine (China) were selected for diffusion test in the paper. Before the experiment, the authors sieved the coal particles using sieves with pore size of 0.25-0.5 mm firstly and then dehydrated the prepared particles to ensure the completely dry state. In experimental process, the gas diffusion rates were obtained by a self-assembled sorption-diffusion system. To fully verify our model, we select the experimental results of diffusion coefficient with different temperatures (30°C and 80°C) at constant pressure (0.5 MPa) to compare with the fitting curves, which is shown in Figure 5(b).From Figure5, it is clear that the theoretical results of diffusion coefficient calculated by our model match the experimental data accurately. And this figure also shows that the diffusion coefficient decreases with diffusion time, and this tendency is more obviously in the early stage of gas diffusion. The conditions of high gas pressure and high temperature can promote the diffusivity of gas significantly. All of the above demonstrate that the proposed nonlinear diffusion model can well describe the gas diffusion process in various kinds of coal under different gas pressure and different temperature. In addition, as shown in Figure 5, for gas diffusion in coal originating from the same producer, the values of parameters in Equation (6) are only slightly different under various diffusion conditions, but for gas diffusion in different types of coal, these values are quite different, indicating that the above parameters are dominated by coal category.


## 3.2. Validation of the New Coupled Model

In this section, a field test data of direct recovery [47] is utilized to compare with the numerical simulation result based on the coupled model. The coalbed methane gas field is located in the south of Shanxi Province, China, and the data originates from October 4, 2007, to September 3, 2008. The site situation, initial state, and boundary conditions of the gas field are shown in Figure 6, and the key parameters used for simulation are listed in Table 2. All the listed parameters are form real CBM field.Figure 6
The site situation, initial state, and boundary conditions of the gas field utilized for validation.Table 2
The key parameters used for the simulation in Section3 are listed as follows.
ParameterValueUnitParameterValueUnitDensity of coal1650kg m-3Young’s modulus3GPaPoisson’s ratio0.3—Initial porosity of fracture0.0302—Permeability of coal1.63×10−13m2Langmuir strain constant0.03—Langmuir pressure constant2.7×106PaLangmuir volume constant0.045—Figure7 depicts the actual data and simulation result of gas production rate. It is obviously that the predictive value of gas production rate obtained by the proposed coupled model has an appropriate agreement with the actual values. Both Figures 5 and 7 indicate that our models are applicable in modeling the gas and heat transfer during long-term direct and thermally enhanced CBM recovery.Figure 7
Comparison of our model and the field test data.


## 4. Numerical Model and Solving Environment

Before the numerical simulation and analysis, the geometric model is established (Figure8) based on the five-spot pattern well for thermally enhanced CBM recovery (Figure 9). In Figure 9(a), the production system consists of one injection well and four production wells, and the distance between injection and production well is 80 m. Generally, the quadrant of the system is extracted to model the process of recovery, as shown in Figure 9(b). To analyze the gas and heat transfer between the two wells efficiently and accurately, in this paper, we simplify the three-dimensional model into a two-dimensional model (ABCD), which is illustrated in Figure 8. In the computational model, four monitoring points (MP) and a monitoring line is set, which are depicted in purple (see on Figure 8). And the boundary conditions of the established geometric model are also shown in Table 3. In addition, by fitting the nonlinear diffusion model to the experimental data of Li et al. [42], in this study, the values of D∞a, qa, A, and B are estimated to 23.35×10−5(cm2/s), 23.04 (kJ/mol), 0.847, and 2.316, respectively, as shown in Figure 3. Some other parameters inputted in the computational model are listed in Table 4.Figure 8
Computational model for simulating.Figure 9
Configuration of the five-spot pattern well for thermally enhanced CBM recovery.Table 3
Boundary conditions of the computational model.
ConditionCoal deformationGas flow in fractureGas diffusionHeat transferBoundary ABRoller boundaryNo-flow boundary—Dirichlet boundary (363.15 K)Boundary BCFixed boundaryNo-flow boundary—Insulated boundaryBoundary CDRoller boundaryDirichlet boundary (0.1 MPa)—Dirichlet boundary (298.15 K)Boundary DAConstant stress boundary (15 MPa)No-flow boundary—Insulated boundaryTable 4
Key parameters of the numerical model.
ParameterValueUnitParameterValueUnitYoung’s modulus (E)3GPaPoisson’s ratio (ν)0.3—Thermal expansion coefficientαT2.4×10−5K-1Langmuir strain constant (εL)0.005—Langmuir pressure constant (pL)2.7MPaLangmuir volume constantVL0.045m3/kgInitial porosity of fracture (φf0)0.03—Porosity of matrix (φm)0.01—Shape factor (χ)1.18×106m2Attenuation coefficient (λ)2×10−8—Velocity of gas (μ)1.84×10−5Pa·sInitial permeability (k0)1×10−16m/sGas specific heat constant (Cg)2160J/(kg·K)Coal-specific heat constant (Cs)1250J/(kg·K)Isosteric heat of adsorption (qst)15000J/MolThermal conductivity of coalλs0.2W/(m·K)


## 5. Results and Discussions

To comprehensively analyze the impacts of dynamic diffusion and heat injection on gas and heat transfer during long-term CBM recovery, three cases of simulations are installed according to different extraction methods and different diffusion models (note: when simulating long-term direct CBM recovery, the boundary condition of AB in Figure8 for heat transfer is insulated boundary, instead of Dirichlet boundary). The representation and comparison of the three cases are illustrated in Table 5.Table 5
Three cases for simulation and analysis.
NameDiffusion modelExtraction methodRepresentationCase 1Traditional diffusion (TD) modelDirect recoveryControl caseCase 2Traditional diffusion (TD) modelThermally enhanced recoveryTo analyze the impact of heat injectionCase 3New nonlinear diffusion (ND) modelThermally enhanced recoveryTo analyze the coupled impact of dynamic diffusion and heat injection


### 5.1. The Variation of Temperature

In this section, the temporal evolution and spatial distribution of temperature under the three cases are revealed to investigate the rule of heat transfer in coal seam.


#### 5.1.1. The Temporal Evolution

Figure10 shows the evolution of temperature with production time at different monitoring points (MP). It can be seen that because of heat injection, the temperature of case 2 and case 3 is obviously higher than case 1. At MP 1, the temperature increases with time after a short decrease. It can be explained that at the early stage of production, the temperature change is dominated by gas desorption, which absorbs energy and lowers the temperature in coal seam as a result of the existence of gas adsorption energy, while at the later stage of production, the temperature change is dominated by thermal conduction and thermal convection. At MP2, the temperature of case 2 decreases with time firstly and then generally levels off. That because the thermal effect of production well is weakened, and the temperature increase caused by thermal transfer and temperature decrease caused by gas desorption are in equilibrium. At MP3 and MP4, compared with case 1, the temperature of case 2 has a huge raise with time due to the strong thermal effect of injection well. Additionally, through the companion between temperature evolution of case 2 and case 3, the model ignoring the gas nonlinear diffusion can underestimate the value of temperature in coal seam, especially in the later production stage. But at MP4, the above phenomenon is unapparent. The reason can be explained as follows: in ND model, the gas diffusion coefficient decreases exponentially with time, whereas in TD model, the diffusion coefficient is constant over time; thus, with time increases, it is more restricted for gas in coal matrix to diffuse into the fracture considering dynamic diffusion. That indicates in the later stage, the gas desorption amount is much less for case 3, resulting in less energy absorption and higher temperature in coal seam.Figure 10
The temporal evolution of temperature for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


#### 5.1.2. The Spatial Distribution

Figure11 shows the distribution of temperature at the monitoring line after different production times. It is clearly that because of insufficient heat supplement, the temperature is always lower than initial temperature for case 1. However, for case 2, the temperature is obviously higher than initial value near the injection well, and the range of heat injection influence expands with the production time, which is 18, 20, 26.5, and 35 meters for 1 year, 4 years, 10 years, and 20 years of production, respectively. It also can be found in Figure 11 that the distribution of temperature near the production well for case 2 has the similar trend and feature with that for case 1 in whole process of recovery, indicating the injected heat still has no significant effect on temperature change near the production well after long-term production. This result may be attributed to the low thermal conductivity of coal. Additionally, the nonlinear diffusion of gas has almost no impact on the spatial distribution of temperature at the early production stage, as shown in Figures 11(a)–11(c). But for the later production stage (Figure 11(d)), the gas diffusion and desorption in case 3 drops off severely based on Equation (6), thus considering the nonlinear diffusion process can cause a lower energy absorption and a higher value of temperature. This phenomenon is more obvious in the middle region between the injection and production wells, where the thermal effects of the two wellheads are weak.Figure 11
The spatial distribution of temperature for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


### 5.2. The Variation of Gas Adsorption Amount

In this section, the temporal evolution and spatial distribution of gas adsorption amount under the three cases are revealed to investigate the rule of gas transfer in coal matrix.


#### 5.2.1. The Temporal Evolution

Figure12 shows the evolution of gas adsorption amount with production time at different monitoring points. For case 1, it can be found that the gas adsorption amount decreases with time, and the value is much higher near the injection well. But for case 2 and case 3, the results are quite different. In Figures 12(a) and 12(b), because MP1 and MP2 are not affected by the injected heat, the gas adsorption amount of case 2 has the similar trend with that of case 1. At MP3 and MP4, the calculated value by case 2 is far below the value obtained by case 1. After 20 years of production, the remaining gas adsorption amount at MP3 and MP4 of case 2 are just 65.7% and 35.2% of these values of case 1. That means the injected heat and the increase of temperature promote the recovery of adsorbed gas, and the enhanced gas production by thermal stimulation are mainly due to the fully extraction of gas in coal matrix (adsorbed gas). In addition, the impact of nonlinear diffusion on gas adsorption content is apparent at MP1, MP2, and MP3 in the later production stage as a result of the exponential attenuation in diffusion coefficient with time. However, as shown in Figure 12(d), there is no obvious difference between the gas adsorption content of case 2 and case 3 at MP4. The reason can be explained as follows. At the beginning of production, the temperature at MP4 has a sudden increase, resulting in most of the adsorbed gas desorbs into the matrix pore at the early production stage, thereafter, the free gas in coal matrix further diffuses into the fracture with a high diffusion coefficient, when the impact of nonlinear diffusion is weak, and at the later stage, in spite of the strong impact of nonlinear diffusion on gas migration from matrix to fracture, the remaining gas content in matrix is relatively less; thus, the real impact of the ND model on gas adsorption amount is also inapparent. It also can be deduced that neglecting the process of nonlinear diffusion can overestimate the gas production in thermally enhanced CBM recovery.Figure 12
The temporal evolution of gas adsorption amount for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


#### 5.2.2. The Spatial Distribution

Figure13 shows the distribution of gas adsorption amount at the monitoring line after different production times. There is no doubt that the gas adsorption amount of case 1 increases with the distance from the production well. But the situations are quite distinct for case 2 and case 3. All the four figures illustrate that the gas adsorption amount of case 2 and case 3 decreases sharply near the injection well as a result of the temperature increase, and this conclusion is consistent with the previous study [21, 48, 49]. At the middle of the reservoir, the gas adsorption content of case 2 is slightly higher than that value of case 1. This phenomenon may be due to the higher gas matrix pressure of case 2 and needs to be further investigated in the future. In addition, because the gas in matrix becomes more difficult to diffuse into the fracture, the impact of ND model on the distribution of gas adsorption amount gradually appears over production time. In Figure 13(d), after 20 years of production, at a distance of 17 m, the gas adsorption amount of case 3 has the same value with that of case 1, namely, critical distance. It also can be concluded that when the distance is less than the critical distance (17 m), the change of gas adsorption amount is dominated by temperature, while when the distance is further than the critical distance (17 m), the change of gas adsorption amount is dominated by dynamic diffusion.Figure 13
The spatial distribution of gas adsorption amount for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


### 5.3. The Variation of Permeability

In this section, the temporal evolution and spatial distribution of permeability under the three cases are revealed to investigate the rule of gas transfer in fracture.


#### 5.3.1. The Temporal Evolution

Figure14 shows the evolution of permeability ratio with production time at different monitoring points. Based on the previous study, the permeability of coal is mainly determined by gas desorption, temperature increase, and gas pressure drop in fracture during thermally enhanced CBM recovery. It can be seen that for case 1, the permeability increases with time at all four monitoring points, indicating the coal shrinkage induced by gas desorption dominates the permeability change during direct CBM recovery. For case 2, the permeability evolution has the similar trend with that of case 1 at MP1, MP2, and MP3, but because of the matrix swelling induced by temperature increase, the value obtained by case 2 is lower. At MP4, the permeability of case 2 also increases with production time, but there is a sudden raise at the early stage as a result of the mass gas desorption induced by temperature increase. In addition, from the results obtained in case 3, the impact of ND model on permeability evolution is obvious. With an increase of time, the result of case 3 has the similar trend with that of case 2. This is because the diffusion coefficient calculated by ND model is not too different from the value derived by TD model at the early production stage, and both the permeability of case 2 and case 3 are dominated by gas desorption. With the further increase of time, the permeability of case 3 has different degrees of decline at the four monitoring points. This phenomenon can be explained that with the rapid decline of diffusion coefficient, the adsorption gas desorbs and diffuses slowly, which results in the reduction of matrix shrinkage and further leads to permeability decrease. That also means the permeability evolution of case 3 at the later production stage is dominated by temperature increase and decrease of fracture pressure.Figure 14
The temporal evolution of permeability ratio for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


#### 5.3.2. The Spatial Distribution

Figure15 shows the distribution of permeability ratio at the monitoring line after different production times. Compared with the result of case 1, the permeability of case 2 is higher in the early production stage but is lower in the later production stage. The former is as a result of mass desorption of gas induced by temperature increase, and the latter is caused by smaller matrix shrinkage. In addition, as the discussion above, the ND model has almost no impact on permeability at the early production stage. But after 10 years of production, considering nonlinear diffusion can obtain a lower value of permeability, which is attributed to the restriction of gas diffusion. After 20 years of production, the curve of case 3 is complex. The permeability increases firstly along with the distance from the injection well, then decreases at the distance of 6 m and finally increases at the location far away from the heat injection well. The minimum permeability ratio is 1.03 at the distance of 25 m, which is induced by the elevated temperature. At the production well, the permeability of case 3 is also higher and is still dominated by the gas desorption.Figure 15
The spatial distribution of permeability ratio for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


### 5.4. The Residual Gas Content in Different Systems of Coal after 20 Years of Production

In this section, the residual gas content in matrix system and fracture system after twenty-year production is estimated to further investigate the heat and gas transfer during thermally enhanced CBM recovery coupled the effect of nonlinear diffusion. And the gas in matrix system is further divided into two parts to comprehensively describe, which are adsorbed gas and free gas. Figure16 shows the residual gas content for different cases after 20 years of production. From the figure, the difference of residual gas in fracture, residual free gas in matrix and residual adsorbed gas between case 1 and case 2 is 0.52%, 0.03%, and 8.26%, respectively. These values prove once again that the fully extraction of adsorbed gas is the most primary reason why heat injection can enhance the CBM recovery. However, the above result is obtained under the assumption that the diffusion coefficient is constant during the recovery process. In fact, the gas diffusion is varied with production time, temperature, and pressure. It is clearly in Figure 16 that considering nonlinear gas diffusion can derive a much higher value of residual gas in matrix (34.77% for free gas and 44.09% for adsorbed gas) but a much lower value of residual gas in fracture (0.99%). This can be explained that with time increases, the gas diffusion coefficient decreases exponentially based on Equation (6), and more gas remains in coal matrix, which results in the less source of gas in fracture and further leads to the lower value of residual gas in fracture. Therefore, how to further extract the adsorbed gas in actual CBM recovery engineering is still an urgent problem and needs to be investigated in the further. It also should be noted that although the residual gas content in matrix of case 3 is higher than that of case 1 (see on Figure 16), it not means the heat injection blocks the gas migration and production, because case 1 also neglects the dynamic gas diffusion process. Additionally, the impact of dynamic diffusion on the process of direct CBM recovery has been studied in our previous work [36].Figure 16
The residual gas content in different systems of coal for the three cases after twenty-year production.


## 5.1. The Variation of Temperature

In this section, the temporal evolution and spatial distribution of temperature under the three cases are revealed to investigate the rule of heat transfer in coal seam.


### 5.1.1. The Temporal Evolution

Figure10 shows the evolution of temperature with production time at different monitoring points (MP). It can be seen that because of heat injection, the temperature of case 2 and case 3 is obviously higher than case 1. At MP 1, the temperature increases with time after a short decrease. It can be explained that at the early stage of production, the temperature change is dominated by gas desorption, which absorbs energy and lowers the temperature in coal seam as a result of the existence of gas adsorption energy, while at the later stage of production, the temperature change is dominated by thermal conduction and thermal convection. At MP2, the temperature of case 2 decreases with time firstly and then generally levels off. That because the thermal effect of production well is weakened, and the temperature increase caused by thermal transfer and temperature decrease caused by gas desorption are in equilibrium. At MP3 and MP4, compared with case 1, the temperature of case 2 has a huge raise with time due to the strong thermal effect of injection well. Additionally, through the companion between temperature evolution of case 2 and case 3, the model ignoring the gas nonlinear diffusion can underestimate the value of temperature in coal seam, especially in the later production stage. But at MP4, the above phenomenon is unapparent. The reason can be explained as follows: in ND model, the gas diffusion coefficient decreases exponentially with time, whereas in TD model, the diffusion coefficient is constant over time; thus, with time increases, it is more restricted for gas in coal matrix to diffuse into the fracture considering dynamic diffusion. That indicates in the later stage, the gas desorption amount is much less for case 3, resulting in less energy absorption and higher temperature in coal seam.Figure 10
The temporal evolution of temperature for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


### 5.1.2. The Spatial Distribution

Figure11 shows the distribution of temperature at the monitoring line after different production times. It is clearly that because of insufficient heat supplement, the temperature is always lower than initial temperature for case 1. However, for case 2, the temperature is obviously higher than initial value near the injection well, and the range of heat injection influence expands with the production time, which is 18, 20, 26.5, and 35 meters for 1 year, 4 years, 10 years, and 20 years of production, respectively. It also can be found in Figure 11 that the distribution of temperature near the production well for case 2 has the similar trend and feature with that for case 1 in whole process of recovery, indicating the injected heat still has no significant effect on temperature change near the production well after long-term production. This result may be attributed to the low thermal conductivity of coal. Additionally, the nonlinear diffusion of gas has almost no impact on the spatial distribution of temperature at the early production stage, as shown in Figures 11(a)–11(c). But for the later production stage (Figure 11(d)), the gas diffusion and desorption in case 3 drops off severely based on Equation (6), thus considering the nonlinear diffusion process can cause a lower energy absorption and a higher value of temperature. This phenomenon is more obvious in the middle region between the injection and production wells, where the thermal effects of the two wellheads are weak.Figure 11
The spatial distribution of temperature for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


## 5.1.1. The Temporal Evolution

Figure10 shows the evolution of temperature with production time at different monitoring points (MP). It can be seen that because of heat injection, the temperature of case 2 and case 3 is obviously higher than case 1. At MP 1, the temperature increases with time after a short decrease. It can be explained that at the early stage of production, the temperature change is dominated by gas desorption, which absorbs energy and lowers the temperature in coal seam as a result of the existence of gas adsorption energy, while at the later stage of production, the temperature change is dominated by thermal conduction and thermal convection. At MP2, the temperature of case 2 decreases with time firstly and then generally levels off. That because the thermal effect of production well is weakened, and the temperature increase caused by thermal transfer and temperature decrease caused by gas desorption are in equilibrium. At MP3 and MP4, compared with case 1, the temperature of case 2 has a huge raise with time due to the strong thermal effect of injection well. Additionally, through the companion between temperature evolution of case 2 and case 3, the model ignoring the gas nonlinear diffusion can underestimate the value of temperature in coal seam, especially in the later production stage. But at MP4, the above phenomenon is unapparent. The reason can be explained as follows: in ND model, the gas diffusion coefficient decreases exponentially with time, whereas in TD model, the diffusion coefficient is constant over time; thus, with time increases, it is more restricted for gas in coal matrix to diffuse into the fracture considering dynamic diffusion. That indicates in the later stage, the gas desorption amount is much less for case 3, resulting in less energy absorption and higher temperature in coal seam.Figure 10
The temporal evolution of temperature for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


## 5.1.2. The Spatial Distribution

Figure11 shows the distribution of temperature at the monitoring line after different production times. It is clearly that because of insufficient heat supplement, the temperature is always lower than initial temperature for case 1. However, for case 2, the temperature is obviously higher than initial value near the injection well, and the range of heat injection influence expands with the production time, which is 18, 20, 26.5, and 35 meters for 1 year, 4 years, 10 years, and 20 years of production, respectively. It also can be found in Figure 11 that the distribution of temperature near the production well for case 2 has the similar trend and feature with that for case 1 in whole process of recovery, indicating the injected heat still has no significant effect on temperature change near the production well after long-term production. This result may be attributed to the low thermal conductivity of coal. Additionally, the nonlinear diffusion of gas has almost no impact on the spatial distribution of temperature at the early production stage, as shown in Figures 11(a)–11(c). But for the later production stage (Figure 11(d)), the gas diffusion and desorption in case 3 drops off severely based on Equation (6), thus considering the nonlinear diffusion process can cause a lower energy absorption and a higher value of temperature. This phenomenon is more obvious in the middle region between the injection and production wells, where the thermal effects of the two wellheads are weak.Figure 11
The spatial distribution of temperature for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


## 5.2. The Variation of Gas Adsorption Amount

In this section, the temporal evolution and spatial distribution of gas adsorption amount under the three cases are revealed to investigate the rule of gas transfer in coal matrix.


### 5.2.1. The Temporal Evolution

Figure12 shows the evolution of gas adsorption amount with production time at different monitoring points. For case 1, it can be found that the gas adsorption amount decreases with time, and the value is much higher near the injection well. But for case 2 and case 3, the results are quite different. In Figures 12(a) and 12(b), because MP1 and MP2 are not affected by the injected heat, the gas adsorption amount of case 2 has the similar trend with that of case 1. At MP3 and MP4, the calculated value by case 2 is far below the value obtained by case 1. After 20 years of production, the remaining gas adsorption amount at MP3 and MP4 of case 2 are just 65.7% and 35.2% of these values of case 1. That means the injected heat and the increase of temperature promote the recovery of adsorbed gas, and the enhanced gas production by thermal stimulation are mainly due to the fully extraction of gas in coal matrix (adsorbed gas). In addition, the impact of nonlinear diffusion on gas adsorption content is apparent at MP1, MP2, and MP3 in the later production stage as a result of the exponential attenuation in diffusion coefficient with time. However, as shown in Figure 12(d), there is no obvious difference between the gas adsorption content of case 2 and case 3 at MP4. The reason can be explained as follows. At the beginning of production, the temperature at MP4 has a sudden increase, resulting in most of the adsorbed gas desorbs into the matrix pore at the early production stage, thereafter, the free gas in coal matrix further diffuses into the fracture with a high diffusion coefficient, when the impact of nonlinear diffusion is weak, and at the later stage, in spite of the strong impact of nonlinear diffusion on gas migration from matrix to fracture, the remaining gas content in matrix is relatively less; thus, the real impact of the ND model on gas adsorption amount is also inapparent. It also can be deduced that neglecting the process of nonlinear diffusion can overestimate the gas production in thermally enhanced CBM recovery.Figure 12
The temporal evolution of gas adsorption amount for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


### 5.2.2. The Spatial Distribution

Figure13 shows the distribution of gas adsorption amount at the monitoring line after different production times. There is no doubt that the gas adsorption amount of case 1 increases with the distance from the production well. But the situations are quite distinct for case 2 and case 3. All the four figures illustrate that the gas adsorption amount of case 2 and case 3 decreases sharply near the injection well as a result of the temperature increase, and this conclusion is consistent with the previous study [21, 48, 49]. At the middle of the reservoir, the gas adsorption content of case 2 is slightly higher than that value of case 1. This phenomenon may be due to the higher gas matrix pressure of case 2 and needs to be further investigated in the future. In addition, because the gas in matrix becomes more difficult to diffuse into the fracture, the impact of ND model on the distribution of gas adsorption amount gradually appears over production time. In Figure 13(d), after 20 years of production, at a distance of 17 m, the gas adsorption amount of case 3 has the same value with that of case 1, namely, critical distance. It also can be concluded that when the distance is less than the critical distance (17 m), the change of gas adsorption amount is dominated by temperature, while when the distance is further than the critical distance (17 m), the change of gas adsorption amount is dominated by dynamic diffusion.Figure 13
The spatial distribution of gas adsorption amount for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


## 5.2.1. The Temporal Evolution

Figure12 shows the evolution of gas adsorption amount with production time at different monitoring points. For case 1, it can be found that the gas adsorption amount decreases with time, and the value is much higher near the injection well. But for case 2 and case 3, the results are quite different. In Figures 12(a) and 12(b), because MP1 and MP2 are not affected by the injected heat, the gas adsorption amount of case 2 has the similar trend with that of case 1. At MP3 and MP4, the calculated value by case 2 is far below the value obtained by case 1. After 20 years of production, the remaining gas adsorption amount at MP3 and MP4 of case 2 are just 65.7% and 35.2% of these values of case 1. That means the injected heat and the increase of temperature promote the recovery of adsorbed gas, and the enhanced gas production by thermal stimulation are mainly due to the fully extraction of gas in coal matrix (adsorbed gas). In addition, the impact of nonlinear diffusion on gas adsorption content is apparent at MP1, MP2, and MP3 in the later production stage as a result of the exponential attenuation in diffusion coefficient with time. However, as shown in Figure 12(d), there is no obvious difference between the gas adsorption content of case 2 and case 3 at MP4. The reason can be explained as follows. At the beginning of production, the temperature at MP4 has a sudden increase, resulting in most of the adsorbed gas desorbs into the matrix pore at the early production stage, thereafter, the free gas in coal matrix further diffuses into the fracture with a high diffusion coefficient, when the impact of nonlinear diffusion is weak, and at the later stage, in spite of the strong impact of nonlinear diffusion on gas migration from matrix to fracture, the remaining gas content in matrix is relatively less; thus, the real impact of the ND model on gas adsorption amount is also inapparent. It also can be deduced that neglecting the process of nonlinear diffusion can overestimate the gas production in thermally enhanced CBM recovery.Figure 12
The temporal evolution of gas adsorption amount for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


## 5.2.2. The Spatial Distribution

Figure13 shows the distribution of gas adsorption amount at the monitoring line after different production times. There is no doubt that the gas adsorption amount of case 1 increases with the distance from the production well. But the situations are quite distinct for case 2 and case 3. All the four figures illustrate that the gas adsorption amount of case 2 and case 3 decreases sharply near the injection well as a result of the temperature increase, and this conclusion is consistent with the previous study [21, 48, 49]. At the middle of the reservoir, the gas adsorption content of case 2 is slightly higher than that value of case 1. This phenomenon may be due to the higher gas matrix pressure of case 2 and needs to be further investigated in the future. In addition, because the gas in matrix becomes more difficult to diffuse into the fracture, the impact of ND model on the distribution of gas adsorption amount gradually appears over production time. In Figure 13(d), after 20 years of production, at a distance of 17 m, the gas adsorption amount of case 3 has the same value with that of case 1, namely, critical distance. It also can be concluded that when the distance is less than the critical distance (17 m), the change of gas adsorption amount is dominated by temperature, while when the distance is further than the critical distance (17 m), the change of gas adsorption amount is dominated by dynamic diffusion.Figure 13
The spatial distribution of gas adsorption amount for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


## 5.3. The Variation of Permeability

In this section, the temporal evolution and spatial distribution of permeability under the three cases are revealed to investigate the rule of gas transfer in fracture.


### 5.3.1. The Temporal Evolution

Figure14 shows the evolution of permeability ratio with production time at different monitoring points. Based on the previous study, the permeability of coal is mainly determined by gas desorption, temperature increase, and gas pressure drop in fracture during thermally enhanced CBM recovery. It can be seen that for case 1, the permeability increases with time at all four monitoring points, indicating the coal shrinkage induced by gas desorption dominates the permeability change during direct CBM recovery. For case 2, the permeability evolution has the similar trend with that of case 1 at MP1, MP2, and MP3, but because of the matrix swelling induced by temperature increase, the value obtained by case 2 is lower. At MP4, the permeability of case 2 also increases with production time, but there is a sudden raise at the early stage as a result of the mass gas desorption induced by temperature increase. In addition, from the results obtained in case 3, the impact of ND model on permeability evolution is obvious. With an increase of time, the result of case 3 has the similar trend with that of case 2. This is because the diffusion coefficient calculated by ND model is not too different from the value derived by TD model at the early production stage, and both the permeability of case 2 and case 3 are dominated by gas desorption. With the further increase of time, the permeability of case 3 has different degrees of decline at the four monitoring points. This phenomenon can be explained that with the rapid decline of diffusion coefficient, the adsorption gas desorbs and diffuses slowly, which results in the reduction of matrix shrinkage and further leads to permeability decrease. That also means the permeability evolution of case 3 at the later production stage is dominated by temperature increase and decrease of fracture pressure.Figure 14
The temporal evolution of permeability ratio for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


### 5.3.2. The Spatial Distribution

Figure15 shows the distribution of permeability ratio at the monitoring line after different production times. Compared with the result of case 1, the permeability of case 2 is higher in the early production stage but is lower in the later production stage. The former is as a result of mass desorption of gas induced by temperature increase, and the latter is caused by smaller matrix shrinkage. In addition, as the discussion above, the ND model has almost no impact on permeability at the early production stage. But after 10 years of production, considering nonlinear diffusion can obtain a lower value of permeability, which is attributed to the restriction of gas diffusion. After 20 years of production, the curve of case 3 is complex. The permeability increases firstly along with the distance from the injection well, then decreases at the distance of 6 m and finally increases at the location far away from the heat injection well. The minimum permeability ratio is 1.03 at the distance of 25 m, which is induced by the elevated temperature. At the production well, the permeability of case 3 is also higher and is still dominated by the gas desorption.Figure 15
The spatial distribution of permeability ratio for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


## 5.3.1. The Temporal Evolution

Figure14 shows the evolution of permeability ratio with production time at different monitoring points. Based on the previous study, the permeability of coal is mainly determined by gas desorption, temperature increase, and gas pressure drop in fracture during thermally enhanced CBM recovery. It can be seen that for case 1, the permeability increases with time at all four monitoring points, indicating the coal shrinkage induced by gas desorption dominates the permeability change during direct CBM recovery. For case 2, the permeability evolution has the similar trend with that of case 1 at MP1, MP2, and MP3, but because of the matrix swelling induced by temperature increase, the value obtained by case 2 is lower. At MP4, the permeability of case 2 also increases with production time, but there is a sudden raise at the early stage as a result of the mass gas desorption induced by temperature increase. In addition, from the results obtained in case 3, the impact of ND model on permeability evolution is obvious. With an increase of time, the result of case 3 has the similar trend with that of case 2. This is because the diffusion coefficient calculated by ND model is not too different from the value derived by TD model at the early production stage, and both the permeability of case 2 and case 3 are dominated by gas desorption. With the further increase of time, the permeability of case 3 has different degrees of decline at the four monitoring points. This phenomenon can be explained that with the rapid decline of diffusion coefficient, the adsorption gas desorbs and diffuses slowly, which results in the reduction of matrix shrinkage and further leads to permeability decrease. That also means the permeability evolution of case 3 at the later production stage is dominated by temperature increase and decrease of fracture pressure.Figure 14
The temporal evolution of permeability ratio for three different cases at (a) MP1, (b) MP2, (c) MP3, and (d) MP4.
(a)(b)(c)(d)


## 5.3.2. The Spatial Distribution

Figure15 shows the distribution of permeability ratio at the monitoring line after different production times. Compared with the result of case 1, the permeability of case 2 is higher in the early production stage but is lower in the later production stage. The former is as a result of mass desorption of gas induced by temperature increase, and the latter is caused by smaller matrix shrinkage. In addition, as the discussion above, the ND model has almost no impact on permeability at the early production stage. But after 10 years of production, considering nonlinear diffusion can obtain a lower value of permeability, which is attributed to the restriction of gas diffusion. After 20 years of production, the curve of case 3 is complex. The permeability increases firstly along with the distance from the injection well, then decreases at the distance of 6 m and finally increases at the location far away from the heat injection well. The minimum permeability ratio is 1.03 at the distance of 25 m, which is induced by the elevated temperature. At the production well, the permeability of case 3 is also higher and is still dominated by the gas desorption.Figure 15
The spatial distribution of permeability ratio for three different cases at monitoring line after (a) one-year, (b) four-year, (c) ten-year, and (d) twenty-year production.
(a)(b)(c)(d)


## 5.4. The Residual Gas Content in Different Systems of Coal after 20 Years of Production

In this section, the residual gas content in matrix system and fracture system after twenty-year production is estimated to further investigate the heat and gas transfer during thermally enhanced CBM recovery coupled the effect of nonlinear diffusion. And the gas in matrix system is further divided into two parts to comprehensively describe, which are adsorbed gas and free gas. Figure16 shows the residual gas content for different cases after 20 years of production. From the figure, the difference of residual gas in fracture, residual free gas in matrix and residual adsorbed gas between case 1 and case 2 is 0.52%, 0.03%, and 8.26%, respectively. These values prove once again that the fully extraction of adsorbed gas is the most primary reason why heat injection can enhance the CBM recovery. However, the above result is obtained under the assumption that the diffusion coefficient is constant during the recovery process. In fact, the gas diffusion is varied with production time, temperature, and pressure. It is clearly in Figure 16 that considering nonlinear gas diffusion can derive a much higher value of residual gas in matrix (34.77% for free gas and 44.09% for adsorbed gas) but a much lower value of residual gas in fracture (0.99%). This can be explained that with time increases, the gas diffusion coefficient decreases exponentially based on Equation (6), and more gas remains in coal matrix, which results in the less source of gas in fracture and further leads to the lower value of residual gas in fracture. Therefore, how to further extract the adsorbed gas in actual CBM recovery engineering is still an urgent problem and needs to be investigated in the further. It also should be noted that although the residual gas content in matrix of case 3 is higher than that of case 1 (see on Figure 16), it not means the heat injection blocks the gas migration and production, because case 1 also neglects the dynamic gas diffusion process. Additionally, the impact of dynamic diffusion on the process of direct CBM recovery has been studied in our previous work [36].Figure 16
The residual gas content in different systems of coal for the three cases after twenty-year production.


## 6. Conclusions

To investigate the gas and heat transfer during thermally enhanced CBM recovery coupled the effect of nonlinear diffusion, we developed a nonlinear diffusion model and a coupled thermal-hydro-mechanical model in this study. In the new diffusion model, the multistage pore structure of coal and Arrhenius equation are taken into account, and the gas diffusion coefficient are varied with production time, temperature, and pressure. Additionally, the coupled model incorporates coal deformation, heat transport, gas flow in fracture, and the proposed gas nonlinear diffusion process. Both the proposed diffusion model and coupled model are validated with experimental results or field test data and shows great applicability in modeling long-term CBM recovery. Through the validations and the results of a series of numerical simulations under different cases, the following conclusions are drawn:(1)
The nonlinear diffusion model can well describe the gas diffusion process in various kinds of coal. For gas diffusion in coal originating from the same producer, the same parameters (D∞a, qa, λ, A, B) in Equation (6) have just slight difference under different diffusion conditions; but for gas diffusion in different types of coal, the values of the same parameters are quite different. These results indicate that the parameters describing the behavior of gas diffusion are dominated by coal category(2)
The simulation results show that the impact of heat injection on gas and heat transfer is mainly embodied near the injection well. The specific performance is that the injected heat could increase the temperature, promote the extraction of adsorbed gas, and inhabit the gas flow in fracture as a result of the temperature-induced coal swelling. In addition, the impact of dynamic diffusion on gas and heat transfer is mainly reflected in the later production stage. The specific performance is that considering gas nonlinear diffusion could advance the increase of temperature, inhabit the gas diffusion in coal matrix, and decrease the coal permeability(3)
The fully extraction of adsorbed gas is the most primary reason why heat injection can enhance the CBM production. However, disregarding the gas nonlinear diffusion process could seriously underestimate the residual gas content in matrix system after long-term production, which might further lead to the overestimate of gas production during thermally enhanced CBM recoveryBased the conclusion of our study, there are two important problems that need to be solved urgently for actual CBM recovery engineering. The first one is how to widen the influence scope of heat injection, in other words, how to enhance the heat transfer efficiency in coal seam. And the other one is how to address the inadequate extraction of adsorbed gas induced by the attenuation of gas diffusivity. Additionally, it also should be noted that the proposed model still neglected the two-phase flow process and the discreteness of coal fracture. These drawbacks remain to be further overcame in the future.



---

*Source: 4013305-2022-05-12.xml*