A new desorption-induced permeability dynamic model and its application in primary coalbed methane recovery

In the process of dewatering and recovery of coalbed methane, coal permeability exhibits a quite unique feature due to the interference of matrix shrinkage and stress effects. A new theoretical dynamic model was proposed for coal permeability based on the assumptions of matchstick geometry of the coal and uniaxial strain condition. Distinct from previous models such as P&M and S&D models, our model relates the gas-sorption-reduced strain to the change of surface energy of coal solids. One of the advantages of this model is that it does not require the sorption-reduced strain as an essential input, and therefore eliminates the related laborious and expensive laboratory measurement. The model was validated by fitting it to two sets of public data and shows an excellent match with the observed data. The results also indicate that our model has a better performance in predicting the permeability dynamics than P&M and S&D models. Additionally, a sensitivity analysis of the effect of input parameters on permeability dynamics was conducted by gray-relation theory, and the initial porosity and reservoir temperature are demonstrated to exert a most distinguished effect on the permeability dynamics. Finally, the proposed model was incorporated into a numerical simulator and successfully applied to conduct a history match of the gas and water production rate in a developed territory.


Introduction
Coalbed methane, as a high-quality new energy, has the characteristics of high calorific value, no pollution, and good safety (Xue 2015). It has developed rapidly in China in the past ten years. Unlike the conventional natural gas reservoirs, coal seams are both source rock and reservoir rock which exhibits a conspicuous dual-porosity property , especially for medium-high-rank coals. Methane is mainly stored on the surface of micropores (matrix system) in an adsorption state under initial reservoir conditions and desorbs and subsequently transfers through macropores (cleat system) once reservoir pressure drops below a certain level (desorption pressure) (Fan et al. 2018). Practically, the permeability of matrix system is negligible while the cleat system is the primary migration channel of coalbed methane (Dong et al. 2019). Stress-induced permeability changes in conventional rocks have been in-depth investigated previously (Wang et al. 2010(Wang et al. , 2012. However, permeability dynamics exhibit unique different characteristics in coals due to gas adsorption/desorption phenomenon. In the process of coalbed methane exploitation, the 1 3 permeability is mainly affected by effective stress and matrix shrinkage. Many studies and production experiences show that permeability is the key to coalbed methane production (Durucan et al. 1986;Tao et al. 2014). Permeability is sensitive to effective stress and matrix shrinkage during CBM exploitation. With the decrease of reservoir pressure, the importance of these two effects on the permeability of coal reservoir changes dynamically . During the depletion of CBM reservoir, a reduction in permeability will occur due to the increasing effective stress imposed to the coal seam ). However, this reduction will be opposed or even offset by the permeability enhancement resulting from the matrix shrinkage associated with methane desorption (Yan et al. 2020).
For in-situ CBM reservoir, it is generally considered to be in uniaxial strain state, and the lateral boundary in this state is fixed (zero lateral strain) (Chen et al. 2015;Liu et al. 2014). Because the overburden is constant when CBM is exhausted, the vertical stress remains constant. Uniaxial strain condition was first proposed by Geertsma in 1966 to describe the deformation of reservoir during depletion (Geertsma 1966). In his work, it was assumed that with continuous production, the reservoir deformation was mainly in the vertical direction because its lateral size was larger than the vertical one (Mitra et al. 2012;Liu et al. 2020). In 1998, Palmer and Mansoori extended this boundary to coalbed methane reservoirs to simulate the permeability changes during depletion (Palmer et al. 1998). A theoretical permeability change prediction model (P&M model) is proposed from the rock mechanics point of view, which considers both effective stress change and matrix shrinkage effect. Since then, based on the assumption of uniaxial strain, a series of permeability dynamics models have been proposed (Table 1). Among which Shi-Durucan (S&D) model and Cui-Bustin (C&B) model gained much attention and popularity (Feng et al. 2016;Shi et al. 2004;Cui et al. 2005). The common characteristics among P&M, S&D, and C&B models lie in that all are developed for uniaxial strain conditions. Liu et al. presented a model that explicitly considers fracture-matrix interaction during coal deformation processes based on a newly proposed internal swelling stress concept (Liu et al. 2010). In this concept, the effect on permeability of partial separation of matrix blocks by fractures that do not completely cut through the whole matrix is first examined. Liu's model is also developed for uniaxial conditions.
Apart from the hypothesis of uniaxial stress condition, the assumption of constant volume of underground formation was also used for developing analytical permeability model. Massarotto et al. suggested that coalbed methane reservoirs are actually under a constant volume condition and stated Empirical formula, constants have no real physical meaning Levine (1996) Lack of clear boundary conditions and strict derivation process Palmer and Mansoori (1996) Overburden rock pressure is constant, uniaxial strain, effective stress coefficient is 1 Overburden rock pressure is constant, uniaxial strain, effective stress coefficient is 1 Physical simulation experiment of coal adsorption expansion is suitable for high rank coal reservoir Cui and Bustin (2005) Overburden rock pressure is constant, uniaxial strain, effective stress coefficient is 1 Fu et al. (2005) Constant temperature, uniaxial strain, regardless of particle compressibility Zhou et al. (2009) Desorption and adsorption only result in volumetric strains, which are equal in all three principal stress directions 1 3 that an essentially constant volume of the coal reservoir was maintained by the methan and tensile strength of the overburden across the reservoir and over time. A simplified model based on this philosophy was lately proposed that is complex in equation form but calls for less parameters (Manssarotto et al. 2009). However, the application of constant volume theory in permeability prediction till now is far behind of that of uniaxial strain conditions. In general, the permeability evolution models can be classified into two categories: porosity-based and stress-strainbased permeability models (Ma et al. 2011). For the porosity-based models, the permeability was correlated with porosity through cubic laws (Palmer and Mansoori 1998;Palmer 2009;Wei et al. 2015), For the stress-strain-based model, the permeability evolution was related to the effective stress variation through the exponential function Durucan 2004, 2005), On this basis, a dynamic model of permeability evolution in CBM well development is established considering the effective stress effect and matrix shrinkage effect (Cui and Bustin 2005;Zhang and Tong 2008;Zhang and Wang 2008;Zhou et al. 2009;Liu and Rutqvist 2010).
However, most of the rock mechanics parameters used in the existing models are measured from the core in the laboratory, which ignores the scale effect of the change of mechanical parameters from the core scale to the reservoir scale, which is easy to cause large errors. In this paper, a new method is introduced based on the surface physicochemical principle (linking the reduced strain of gas adsorption with the change of coal solid surface energy). This method does not need to take the reduced strain of adsorption as the necessary input, thus eliminating the associated laborious and expensive laboratory measurements. Proposed a new theoretical dynamic model for coal permeability during CBM production based on the assumptions of matchstick geometry of the coal and uniaxial strain condition (Meng et al. 2018). A mathematical model of strain permeability under the coupling action of effective stress and matrix shrinkage is established. The gas-sorption-induced strain is treated as directly analogous to the change of surface energy of coal solids. Model validation was conducted by fitting the public data, followed by sensitivity analysis. The model was then incorporated into a CBM reservoir simulator and succeeded in history matching gas and water rate for a CBM well in Ordos basin.

Permeability model
As in Fig. 1, the coal seams are assumed to be ideal bundledmatchstick geometry (Shi et al. 2004).
The matrix shrinkage associated with methane desorption in coal seams can be treated as directly analogous to the volumetric thermal contraction of thermoelastic porous medium (Palmer et al. 1998). According to the derivation of the stress-strain relationship for a homogenous, isotropic thermoelastic porous medium, the stress-strain relationship for isothermal coal seams can be written as where α is the Boit coefficient; σ c is the confining pressure, MPa; ε b is the bulk volumetric strain, MPa; G is the shear modulus, MPa; ν is the drained Poisson's ratio, K is the macroscopic bulk modulus, MPa; K s is the solid grain bulk modulus, MPa; p is the pore pressure, MP; ε V is the macroscopic volumetric matrix shrinkage strain induced by coalbed methane desorption; δ ij is the Kronecker delta, is zero when i ≠ j and one when i = j; subscript i, j = x、y、z.
Note that the shear modulus G and the bulk modulus K are related to Young's modulus E and Poisson's ratio ν, respectively, through isotropic elasticity theory (Palmer et al. 1998).
Substituting Eqs. (5) and (6) into Eq. (9), and solving for ε b yields Considering the coal sample as a porous medium containing solid volume V s and pore volume V p , the bulk volume can be defined as V b = V p + V s and the porosity is given byφ = V p /V b . According to Eq. (12), the volumetric evolution of coal-bearing porous medium can be described in terms of the bulk volume strain and the pore volume strain, respectively. The bulk volume strain increment and the pore volume strain increment can be expressed as (Zimmerman et al. 2000) where β is a dimensionless effective stress coefficient. According to the bulk modulus characteristic of the porous frame defined by Brown and Korringa, the pore compressibility can be expressed as (Zimmerman et al. 2000;Berryman 1992) The pore bulk modulus K p can be shown to be dependent of K, K s , and the porosityφ for linear poroelasticity (Berryman 1992; Zimmerman et al. 1986). The relation is known in general that According to Eq. (2), we can write Eq. (17) as Using the definition of porosity, the porosity change of deforming coal seams can be described in terms of the bulk volume strain and the pore volume strain Substituting Eqs. (2), (13), (14), and (18)

into Eq. (19) yields
Generalized Terzaghi's effective stress law as follow (Nur and Byerlee 1971) where σ e is the volumetric effective stress and σ ij e is the effective stress tensor.
Using this definition, Eq. (20) can be written as Assuming constant K and K φ , integrating Eq. (23) yields where 0 is the initial fracture porosity, e 0 is the initial effective stress, and p 0 is the initial reservoir pressure.
Furthermore, Eq. (24) can be approximated using Taylor's expansion as: On the assumption that reservoirs are under uniaxial strain conditions, that is, ε xx = ε yy = 0, the effective horizontal stress σ xx e or σ yy e is gotten as follows: According to Eqs. (3), (21), and (22) Substituting Eq. (26) and Eq. (28) into Eq. (27), and solving for e xx yields Transposing and reducing, one obtains and the effective stress increment becomes The substitution of Eq. (30) into Eq. (25) leads to In the limit of incompressible grains, that is, → 1 , thus, a model similar to the P&M model is yielded: To derive the coal matrix shrinkage strain, we use the Bangham equation, which relates the change of strain of a porous body with a change of surface energy. This thought of equivalence was first mentioned and verified by Bangham and his associates (Bangham 1952). Their experimental results showed that upon the reduction of the free surface energy by the adsorption of gases and vapors, as calculated from the Gibbs adsorption isotherm, the solid increases in size by amounts of the order of 0.5% in one dimension. It has been shown in the case of charcoal and coal that, within the limits of accuracy of the experiments, the percentage linear expansion of the solid is directly proportional to a change in the free surface energy per unit area. The form of the equation is: where ε is the strain change of coal solid; ρ c is the density of coal, t/m 3 ; S is specific surface area, m 2 /t; △γ is the change of surface energy of coal solid, J/m 2 , which is caused by the adsorption/desorption of gas to/from the surface of coal solid. According to Gibbs equation, △γ can be written as: where γ 0 and γ is the surface energy of coal solid under vacuum condition and after adsorbing gases, respectively, Γ = V∕V0S , the difference between surface concentration and bulk phase concentration; combining Eq. (33) with Eq.
(34) results in: V in Eq. (35) can be gained from equation of Langmuir isotherm adsorption: For the case of CBM production, during the process of depressurization when pressure falls from critical desorption pressure p r to p, the change of strain is obtained from Eq. (37): where V L is Langmuir volume constant, m 3 /t; b is Langmuir pressure constant, MPa −1 ; T is absolute temperature, K; R is universal gas constant, 8.3143 J/(mol⋅K); V 0 is gas molar volume under standard condition, 22.4L/mol; p r is critical desorption pressure, MPa. The following relationship for the effective stress change is obtained: Assume that the fracture permeability varies as the porosity as follows (Pan et al. 2012) Our model bears a similarity with P&M (Eq.41) and S&D model (Eq.42) in the composition of equation, i.e., all models include two parts controlling permeability dynamics, the stress effect sector, and matrix shrinkage sector. However, the theoretical basis of our model is reliable, the form is simple and the parameters involved in the model are easy to obtain, and hence, eliminates the elaborate as well as expensive work of determination of lithologic parameters including constrained modulus M bulk modulus K and isotherm parameter sorption-induced strain ε l (ε g in S&D model); On the other hand, the commonly used prediction models such as P&M model is established based on the assumption that porosity is much less than 1, for coal reservoirs of the high porosity, the calculation of the model will appear bigger error, while our model to make up for the defects, also can be applied to coal reservoir of the high porosity, which has a much wider application scope.

CBM reservoir model
The following set of partial differential equations are given to formulate the water and gas flows (Zhang et al. 2009): where λ is the unit conversion factor; p l is the fracture pressure for the l-phase; k f is the fracture permeability; k rl is the l-phase relative permeability; μ l is the l-phase viscosity; B l is the l-phase formation volume factor; ρ l is the l-phase density; S l is the l-phase saturation for the fracture system; t is time, ϕ f is the fracture porosity; q vl is the l-phase volumetric flow rate per unit volume reservoir; l = g, w; H is the elevation; V m is the average gas concentration; τ is adsorption time; q vmf is the desorbed gas flow rate through the coal matrix per unit volume reservoir; ρ c is the coal density; and V E is the equilibrium methane concentration, which can be described by Langmuir isotherm as follows:V E V L , p L are the Langmuir volume and pressure constants, respectively.
In this paper, we use the closed outer boundary and constant pressure inter boundary conditions. The auxiliary equations are given as: where p c is the capillary pressure between gas and water phases.
In time, we use a fully implicit formulation solved at every time step with an Orthomin method preconditioned with incomplete LU decomposition (Zhang et al. 2009), in which pressure and saturation are both solved implicitly. At each iteration step, the porosity and permeability (along three dimensions of x, y, z) is updated with Eqs. (39) and (40), respectively.

Validation with field data
MCGovern once presented a general pressure-dependent permeability multiplier curve for the San Juan basin coalbeds (Fig. 2), which turns out to be cut out for model validation in this paper. Figure 1 demonstrates a permeability ratio k/k 0 up to a factor of 7 as reservoir pressure drops from 5.516 to 0.069 MPa after the reservoir is largely dewatered (MCGovern 2005). It's noted that the data by McGovern was once used by Shi and Durucan for verifying their model and the value of reservoir parameters used in this study are primarily determined based on Shi's work (Table 2).
Though coal seams are practically under-saturated at initial conditions (the volume of adsorbed methane is less than that corresponding to initial pressure p 0 on Langmuir curve), in this study, we found that an assumption of the critical pressure p r to be 8.0 MPa yields a good agreement with the field data. This assumption is to, some extent feasible due to large degree of depressurization. As a matter of fact, validation of Shi's model was achieved by varying initial reservoir pressure from 8.27 to 11.38 MPa.
To derive the Langmuir volume, we directly borrowed the α s (CH 4 ) = 0.00038 that was used for history matching the flue-gas-injection micropilot in Shi. This leads to a V L of 34 m 3 /t. It is worth pointing out that laboratory researches on coal sorption properties have demonstrated a phenomenon referred to as adsorption/desorption hysteresis, i.e., (Bush et al. 2003), the desorption curve does not follow exactly the adsorption curve under decreasing and increasing pressure, respectively. This hysteresis results in a decrease in Langmuir pressure in depressurization process compared with that in pressurization process. Since the primary recovery of CBM from underground reservoir is a depressurization  process, it is recommended that desorption curve be applied when simulating permeability performance during primary recovery. Therefore, the value of P L used by Shi (4.31 MPa) was induced to 3 MPa to better simulate the permeability performance (Shi et al. 2004).

Comparison with previous models
Palmer and Mansoori reported a good history matching of the bottom-hole pressure and gas production rate for a CBM well in San Juan basin with a permeability curve with strong rebound at lower pressure drawdown, which conforms to curve 3 in Fig. 3. Later, Shi and Durucan gave the permeability responses with S&D and P&M models based on the history matched data. In their research, the initial porosity and pressure were set to 0.00085 and 10.4 MPa, respectively, and the simulation results are plotted in Fig. 3. We can see that the S&D model generates a much stronger rebound at lower drawdown pressure than the P&M model, but both failed to yield a good match.
Mavor and Vaughn reported in-situ permeability changes from drill-stem and shut-in tests in three wells in San Juan Fruitland formation coal seams (Mavor et al. 1998). As is noted previously, isotherm parameters exhibit some kind of difference due to adsorption/desorption hysteresis, and hence, we modified Mavor's V L to 31m 3 /t and P L to 2.8 MPa. Assume a normal pressure and temperature gradient, an insitu pressure of 10.4 MPa corresponds to a temperature of about 330 K. The simulation result was demonstrated in Fig. 3. It is indicated that the newly developed model yields a much better match compared with the above two models.

Sensitivity analysis and discusses
The influence of different factors on the permeability ratio k/k 0 was evaluated by the gray relation theory (Liu et al. 2004), and the closeness of correlation between different sequences is determined by the level of similarity of their corresponding geometric configurations. In this research, the sequence of influence factor i (i = 1, 2,…, 9) and that of k 1 /  1 3 k 0 (k 1 corresponds to a reservoir pressure of p 1 = 0.5 MPa) are set as the behavior sequence X i and the feature sequence X 0 , respectively. The integrated incidence degree is adopted as a quantitative index to represent the closeness of each behavior and feature sequence. The influence factor with a greater integrated incidence degree exerts a more profound effect on permeability when reservoir pressure drops from p 0 to p 1 . The base value and range of each factor are listed in Table 3. Figure 4 shows the value variation of integrated incidence degree between k 1 /k 0 and influence factors. We can see that the initial porosity, coal density, and temperature have relatively high integrated incidence degrees (0.638 at the lowest), which indicates that these three parameters exhibit a dominant influence on the permeability model.
To examine the performance of permeability dynamics, we tried different combinations of parameters and found four sets (Table 4) are quite typical in generating four distinctive dynamic curves (Fig. 5). Coal seams were assumed to be saturated by pure methane, and hence, desorption pressure was equal to initial reservoir pressure of 6.1 MPa. Both adsorption and desorption processes are simulated. In the region p < 6.1 MPa, curve 1 displays a proximately exponential increase with decreasing pressure, and results in an increase of more than fourfold, while curve 2 displays a progressive decline. In case 3, permeability exhibits descending trend at initial pressure drawdown period and then rebounds strongly at lower pressure drawdown (about 1 MPa) and up to a factor of 2.16 at 0.15 MPa, more than doubling the initial permeability. Curve 4 shows a relatively "flat" trend in both adsorption and desorption processes, indicating little variation of permeability.
The occurrence of varied shapes of dynamic curves is probably attributed to the mutual interference of stress effect and matrix shrinkage. In the primary recovery of CBM reservoir, pore pressures decline gradually, resulting in an increase in effective stress and methane desorption. The former effect tends to shrink the pore volume while the latter can enhance that. If those two effects are evenly matched, permeability will exhibit negligible changes (Curve 4). Once the stress effect is much stronger than matrix shrinkage, permeability will decline without rebound. By contrast,

Field application
Hancheng area in the southeast margin of Ordos basin is one of regions possessing coalbed methane exploration and development potential. The 3 # coal seam in the Permian Shanxi Formation and the 5 # , 11 # coal seams in the Carboniferous Taiyuan Formation are the main target formations for coalbed methane recovery. History matching of gas and water production rate for a CBM well WL1 in Hancheng area was conducted by a coalbed methane simulator incorporated with the proposed permeability model. The 3 # , 5 # coal seams are penetrated and drained by well WL1. Reservoir parameters are listed in Table 5.
The simulated gas and water production curves, as shown in Fig. 6, are consistent with the field observations, indicating a fairly fine performance of the proposed model.
The dynamic change curves of coal permeability were illustrated in Fig. 7 during the coalbed methane development, which demonstrates an overall higher permeability in the near-wellbore region than in the far-field region. At the early stage of dewatering, permeability experiences a 10% decrease in near-wellbore region, which affects the drainage of water adversely. After about 150 days of drainage, permeability in this region begins to rebound. It is also the case with far-wellbore region except for the beginning time and amplitude of the rebound.
It is also interesting to point out that fluctuation of gas production rate has a close relation to near-wellbore permeability dynamic changes. Comparing Fig. 6(b) with Fig. 7, a fairly fine correspondence between those two curves can be found, i.e., the ascension/decline sections on gas rate curve have a corresponding congruent trend on permeability curve in near-wellbore region. After a drainage period of 150 days, permeability in near-wellbore region begins to rebound. We can also see that after about 5.5 years of drainage, permeability rebounds and then exceeds initial value (k/k 0 > 1) and average gas production rate climbs from less than 2000m 3 /d to 2300m 3 /d.

Conclusion
A new dynamic model of coal permeability was developed by the physical chemistry theory that relates desorption to changes of surface energy and the equivalence assumption between sorption-induced and surface-induced strains. The advantages of the new model over the previous include demand for less input parameter and allowing for application under variable temperatures. Model validation with published data indicates a fine performance of the new model. Following incorporating the model to a CBM numerical simulator, we succeeded in history match gas and water production rate for a CBM well in Ordos basin. Simulation results showed a close relation between the gas production rate changes and permeability dynamics.

Conflict of interest
The authors declare that they have no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.