Numerical simulation of a coupled gas flow and geomechanics process in fractured coalbed methane reservoirs

The pores and throats of coal seam are narrow, and the permeability is very low. It is a porous medium with complicated deformation mechanism. In the process of gas production, the change of stress, pore pressure, gas adsorption and desorption will lead to the change of grains, pore volume, and the seepage capacity of coal seam. Based on the theory of poroelastic and seepage mechanics, a physical model of fractured horizontal well and a coupled mathematic model of coal deformation and gas flow are established. The dynamic change model of permeability is deduced too. The development characteristics and laws of coalbed methane are numerically simulated by using numerical analysis software Comsol Multiphysics. Finally, the sensitivity of the key parameters in the model and the influence of the fracture parameters are analyzed. The results show that the application of fractured horizontal well is an effective measure to increase production of coalbed reservoirs; with the increase of Langmuir volumetric strain constant, fracture half‐length, and the numbers of fracture, the production of coalbed methane is increasing. In the case of constant total fracture length, number, and spacing, when fracture half‐lengths are long and short interlaced distribution, the production is the highest and it is the best fracture distribution pattern. The study provides a theoretical basis for the further study of coalbed deformation, gas migration, efficient, and rational development of coalbed methane reservoirs.

and pore volume. That will affect the change of permeability and ultimately affect the productivity of CBM. 4,5 At the same time, the pore throat of coal seam is narrow and the permeability is very low. 6 So, hydraulic fracturing and horizontal well are the key approaches to improve the reservoir seepage ability which increases the complexity of gas seepage in reservoir and greatly increases the difficulty of research. 7,8 In previous studies, many scholars have made a lot of researches on the effect of stress and adsorption on seepage capacity of coal seam. [9][10][11][12][13] However, there are relatively few studies on dynamic coupling and productivity of fractured horizontal well development in coal seam. Therefore, to study the coupled model of coal seam deformation and gas flow in fractured horizontal well is very important to the seepage mechanism, influencing factors and productivity analysis of CBM development. Therefore, it is necessary to further study the production performance and productivity of fractured horizontal wells with the effect of coal deformation and gas flow.
On the basis of multi-physics coupling research, a fully coupled mathematical model of coal deformation and gas flow is established in this paper, and the dynamic models of porosity and permeability are derived. Finally, the numerical analysis software Comsol Multiphysics is used to solve the model by finite element method. The sensitivity of the key physical parameters and the complex mechanism of gas flow in the fractured horizontal well are analyzed, and the productivity under different fracture parameters and geometries is studied.

| Physical model
In this paper, the poroelastic theory is used to define and describe the mathematical coupled model of coal seam deformation, gas flow, and their interaction in coal seam system. 14,15 Based on the following assumptions: 1. Coal seam is a homogeneous, isotropic, and elastic continuum.

2.
Strains are much smaller than the scale of coal seam. 3. Under the isothermal condition, the gas in the pore is ideal, and its viscosity is constant. 4. Coal is saturated by gas. 5. The rate of gas flow through the coal is defined by Darcy's law and neglects the influence of gravity.

| Governing equation of coal seam deformation
For a homogeneous, isotropic and elastic medium, the total strain of coal contains the strain induced by effective stress and sorption of coal seam. The total strain of coal can be expressed as (tensile is positive) 14 : where ɛ ij is the component of the total strain tensor; G is the shear modulus of coal, G = E∕2 (1 + v); σ kk = σ 11 + σ 22 + σ 33 ; K is the bulk modulus of coal, K = E∕3 (1 − 2v); E is the Young's modulus of coal; v is Passion's ratio of coal; and p is the gas pressure within the pores. The component of effective stress is defined as � ij = ij + p ij ; α = 1 − K/K s is the Biot coefficient; K s is the bulk modulus of the coal grains δ ij ; is the Kronecker delta.
The sorption-induced volumetric strain ɛ s is fitted onto Langmuir-type equation 16,17 : where the Langmuir volumetric strain, ɛ L , is a constant representing the volumetric strain at infinite pore pressure and the Langmuir pressure constant, p L , representing the pore pressure at which the measured volumetric strain is equal to 0.5 ɛ L . 18 Combining with the equilibrium equation σ ij,j + f i = 0 and Cauchy's equation ij = u i,j + u j,i ∕2, the governing equation for coal deformation can be expressed as: where u i is the component of the displacement; σ ij denotes the component of the total stress tensor; and f i denotes the component of the body force.

| Governing equation of coal gas flow
The equation for mass balance of the gas is defined as: where ρ g is the gas density, kg/m 3 ; v g is the Darcy velocity vector, m/s; Q s is the gas source or sink, kg/(m 3 · s); t is the time, s; and m is the gas content, kg/m 3 , including free-phase gas and absorbed gas 19 : where ϕ is porosity; ρ ga is the gas density at standard condition, kg/m 3 ; ρ c is the coal density, kg/m 3 . According to the ideal gas law and Darcy's law, we obtain: where p a is one atmosphere of pressure (101.325 kPa); k is the permeability, m 2 ; and μ is the dynamic viscosity of the gas, (N · s)/m 2 .
The permeability k is dependent on the porosity, ϕ, which is a function of the volumetric strain, ɛ v , and the sorptioninduced strain, ɛ s . These relationships will be discussed in the next section.

| Porosity and permeability model General porosity model
We assume that the sorption-induced strain for the coal is the same as that for the pore space. If the initial porosity is ϕ 0 at initial pressure p 0 and the initial volumetric strain is zero. According to the definition of coal volumetric strain and porosity, the expression of porosity of coal seam under the interactions of coal deformation and gas flow is defined as 14 :

Porosity model under uniaxial strain condition
In the process of depleting development mode, the coal seam is restrained laterally, and the lateral strains can be ignored. The coal seam is under the conditions of uniaxial strain and constant overburden. They are the most widely applied assumptions used to simplify the geomechanical description. 20 According to the actual situation in the field, assuming that ɛ 33 is the direction of uniaxial strain and overburden load which is basically unchanged, so ɛ 11 and ɛ 22 are the lateral strains that both equal to 0. 11,[20][21][22][23][24] At this time, the volumetric strain of the coal seam is as follows: where M (constrained axial modulus) is related to Young's modulus, E, and Passion's ratio, v, Combining with the equation (7), the porosity model under the conditions of uniaxial strain and constant overburden load is derived as: Considering the cubic law, we can obtain the expression of permeability under the conditions of uniaxial strain and constant overburden load: The expression of the fully coupled model is similar to the classical P-M permeability model. But the P-M model only considers the grains bulk modulus as infinity, neglecting the effect of the grains deformation on the pores space.

Rebound pressure
By observing the equations (9) and (10), we can find that the changes of coal seam porosity and permeability are affected by the compression of coal and coal grains, and affected by matrix shrinkage. Palmer and Mansoori define the pressure which the porosity and permeability begin to increase as rebound pressure. That is the pressure matrix shrinkage begins to dominate the processes of porosity and permeability changes. Its expression is 21 Similarly, the expression of the P-M rebound pressure does not reflect the effect of the grains compression, by using the similar method, we can deduce the rebound pressure of the fully coupled porosity model considering the influence of the grains compression. Under the conditions of uniaxial strain and constant overburden load, the expression of rebound pressure is It is worth noting that the rebound pressure does not depend on the initial pressure p 0 of the coal seam.

Coupled governing equations
The coupled governing equation (3) for coal seam deformation has been obtained in 2.1. For the coupled governing equation of gas flow, the partial derivative of porosity ϕ with respect to time t can be obtained by equation (7): Putting equation (13) into (6), the fully coupled governing equation for the gas flow affected by the coal seam deformation can be obtained: Therefore, Equations (3), (10), (13), and (14) define the fully coupled model of coal seam deformation and gas flow. The coupling relationships between the coal seam deformation and gas flow are shown in Figure 1.

RESULT ANALYSIS
The above governing equations constitute the mathematical model describing the coal deformation and gas flow. In this paper, the finite element method is applied to solve the numerical model by Comsol Multiphysics software. It is assumed that the length of the rectangular coal seam is the same as the horizontal well length, and the coal seam is symmetrical about the horizontal well. All the properties and laws of the two sides of the horizontal well are the same, so this paper only takes the analysis of the 1/2 rectangular coal seam. The model here is a 500 m × 200 m rectangle. The bottom boundary pressure is specified as bottom hole pressure (p w = p a ), and zero fluxes on the other three sides are specified. The initial gas pore pressure in the coal is set at 6.2 MPa. Nine fractured fractures with infinite conductivity are symmetrical about horizontal wells. For reservoir simulations with the 2D geometry, the hydraulic fractures are treated as one-dimensional line by representing the fractures as the interior boundaries of matrix. 25,26 The geometric model of reservoir and the main parameter used in calculation are shown in Figure 2 and Table 1, respectively.

| Model verification and parameter analysis
The mathematical model established in this paper is applied to the physical model and solved numerically. The results are shown in the figures. Figures 3 and 4 are the pore pressure distribution of coal seam and the pore pressure of two points at different times, respectively. The point 1 (225, 50) is near the fractures and the point 2 (225, 150) is far away from the fractures area. It can be found that the minimum value of gas pressure is distributed around hydraulic fractures. With the increase of time, the pore pressure is gradually decreased, and the pressure around the fractures decreases faster than that of the area far away from the fractures. Figure 5 is the permeability ratio distribution of coal seam at different times. Through observation, it is found that the overall permeability of the coal seam has increased, and the permeability ratio around the fractures is larger than that in the area far from the fractures. This phenomenon is due to the rapid decline of gas pressure around the fractures, gas first desorbs from the matrix around the fractures, and then the matrix around the fractures starts to shrink. Finally, when the gas pressure is reduced to the bottom hole pressure, the gas pressure and permeability all tend to be equal.  Figure 6 compares the relationship between permeability ratio (k/k 0 ) and gas pressure under the three conditions of fully coupled model, P-M model and the model without desorption-induced deformation. It is found that the permeability of the fully coupled model and the P-M model increases with the decrease of pore pressure. And the permeability of the model without adsorption deformation decreases with the decrease of pore pressure. With the decrease of gas pressure, the increase of effective stress will lead to the compression of the whole coal seam and reduce the permeability. At the same time, if the adsorption deformation is considered, the gas desorption will lead to the shrinkage of the coal matrix, which will increase the permeability. The two factors jointly control the change of permeability. The increase of permeability is due to the mechanism that the effect of desorption deformation is greater than the effect of effective stress on coal compression. [27][28][29][30] Because the P-M model ignores the deformation of coal grains, the permeability calculated by this model is a little larger. When the absorption deformation is not taken into account, the change of permeability is controlled only by effective stress. The increase of effective stress compresses the coal seam and reduces the permeability. As shown in Figure 7, when the Langmuir volumetric strain constant ɛ L is 0.01 and 0.015, respectively, the permeability decreases first and then increases with the decrease of pore pressure. At first, the shrinkage deformation caused by desorption is small, and the compression of coal seam plays a leading role, so the permeability is less than the initial permeability. When the pore pressure continues to decrease, the influence of the matrix shrinkage increases gradually. After the turning point, the advantage of the matrix shrinkage is gradually larger than the compression of the coal seam, so the permeability increases gradually. The turning point in Figure 7 is the point where the permeability begins to increase after the decrease. The turning point pressure can be calculated by the rebound pressure formula given by formula (12). By calculating the pressure of ɛ L at 0.01 and 0.015, the corresponding rebound pressures are 2.36 and 4.27 MPa, respectively, which are consistent with the results obtained in numerical simulation in Figure 7. When the ɛ L is >0.2, the effect of the matrix shrinkage is greater than that of the effective stress compression at the beginning, so the calculated rebound pressure is higher than the initial pore pressure p 0 . When the pore pressure is gradually on decrease, the larger the ɛ L is, the greater the effect of matrix shrinkage is, and the larger the permeability ratio (k/k 0 ) is. The results of Figure 8 show that the improvement of final permeability of coal seam with lower initial porosity is better. Therefore, it is more necessary to consider the effect of coal seam deformation on the seepage capacity of coal seam with lower porosity.  Figure 9 compares the production of CBM in three conditions:

| The effect of different conditions and parameters on production
1. Fully coupled condition: coupled the process of coal deformation and gas flow which is derived in this study; 2. Uncoupled condition: ignoring coal deformation that means porosity and permeability are constant in the process of gas production; 3. Without adsorption deformation condition: considering ɛ s = 0 which means ignoring the effect of sorption-induced strain on porosity and permeability.
The results show that the fully coupled model predicts the highest daily and cumulative productions of CBM. And when the adsorption deformation is not considered, the predicted daily and cumulative productions are the least. It is because with the gas production, the permeability of the fully coupled model is gradually increasing, which effectively improves the seepage capacity of the coal seam. Without considering the adsorption deformation of coal seam, the whole coal seam has only compression effect during the process of gas production, and the permeability is gradually decreasing, so the production is the lowest. It can be found from Figure 10 that the larger the Langmuir volumetric strain constant is, the stronger the matrix shrinkage is and the better the improvement of permeability is, so the cumulative production of coalbed methane is higher.

| The effect of hydraulic fracture parameters on production
As shown in Figure 11, when the half-length of the hydraulic fracture increases from 0 m (without hydraulic fracturing) to 120 m, the fractures contact the further area of matrix, which enlarges the area of permeability improvement of coal seam and increases the cumulative production of coalbed methane. When the coal seam is not fractured, the production of coalbed methane is very low. However, after fracturing the coal seam, the daily and cumulative productions of coalbed methane have greatly improved, and when the hydraulic fracture half-length is 120 m, the cumulative production is 17.56 times higher than that of 0 meters of the hydraulic fracture half-length. It is obvious that hydraulic fracturing plays an important role in increasing the production of CBM reservoir. The result of Figure 12 shows that when the fracture spacing and half-length are constant, with the numbers of fracture (N f ) increasing from 3 to 9, the area of permeability improvement has expanded along the horizontal well direction. And it has a significant effect on the cumulative production of CBM.

| The effect of fracture geometry on production
For any gas reservoir, considering the production problem after hydraulic fracturing, there is optimal fracture geometry. In this paper, the four kinds of fracture geometry are compared and analyzed in the case of constant fracture length, number, and spacing. In order to study the influence of fracture geometry on coalbed methane production, and select the best fractures distribution forms. Assuming that the nine fractures are symmetrical about horizontal well, the four kinds of fracture geometry are as follows (Figure 13): (a) Fracture half-length decreases along the middle to both ends; (b) Fracture half-length increases progressively along the middle to both ends; (c) Fracture half-lengths are long and short interlaced distribution; and (d) Fracture half-length is equilong and uniform distribution. As can be seen from Figure 13, the forms and areas of pressure transmission of the four geometries are different. Among them, the pressure spread area of geometry 3 is the largest; and the other three geometries are smaller in the area of pressure propagation and the pressure drop, so there are more areas that the pressure drop has not spread. As shown in Figure 14, geometry 3 is the best way to contact with the matrix, and the area of permeability improvement is also the largest. However, the other three geometries have larger areas of poor permeability improvement far away from the horizontal well, and the improvement is worse than that in geometry 3. Figure 15 shows that the daily and cumulative productions of geometry 3 are higher than those of the other three geometries, and the geometry 4 has the minimum daily and cumulative productions. The longer fractures of geometry 3 guarantee the good communication of fractures to the areas far away from horizontal well. At the same time, the shorter fractures distributed between the longer fractures make up for the lack of communication caused by too much spacing between the two longer fractures. Therefore, we should choose the fractures distribution of geometry 3 to carry out coal seam hydraulic fracturing to achieve the best productivity.

| CONCLUSION
Based on the theory of poroelastic, seepage mechanics, and considering adsorption and desorption deformation, a fully coupled mathematical model of coal deformation and gas flow is established, and the dynamic models of porosity and permeability are derived. The models are applied to the numerical simulation of the development of fractured horizontal well in coal seam.

1.
The permeability of coal seam increases gradually with the decrease of pore pressure, and the increase area expands gradually from the fractures to the matrix. Eventually, the permeability tends to be equal. Without considering the grains deformation, the permeability ratio (k/k 0 ) obtained by the simulation is a bit higher. If the adsorption deformation is not considered, the permeability ratio (k/k 0 ) will be relatively lower. The larger the Langmuir volumetric strain constant ɛ L is, the higher the permeability ratio (k/k 0 ) is, and the higher the pore pressure is when permeability begins to increase. The lower the initial porosity is, the better the permeability of coal seam will eventually be improved. Therefore, the lower porosity is, the greater the effect of coal deformation on seepage capacity is. 2. The fully coupled model has the highest daily and cumulative productions, and the daily and cumulative productions are the least when the model is not taken the adsorption deformation into account. The bigger the Langmuir volumetric strain constant ɛ L is, the better the permeability of coal seam will be. Therefore, the more the adsorption deformation of coal seam is, the higher the cumulative production of coalbed methane will be. 3. When the fracture half-length increases from 0 meters (without hydraulic fracturing) to 120 meters, the production of coalbed methane will also increases. If the fracture spacing and half-length are constant, the cumulative production of CBM will increase with the increase