An internal swelling factor model to examine the influence of permeability anisotropy on coalbed methane extraction

Owing to the characteristics of coal reservoirs, the gas flow capacity and permeability exhibit strong anisotropy. The anisotropy in terms of the magnitude, which corresponds to the permeability in the horizontal direction being several orders of magnitude larger than that in the vertical direction, has been investigated. However, the anisotropy in terms of the mechanical boundaries, specifically, the presence of constant volume and stress boundaries in the horizontal and vertical directions, respectively, has not been examined. Therefore, a coupling model of the gas flow and coal seam deformation was developed, and a model to reflect the permeabilities in the horizontal and vertical directions was established considering the internal swelling factor. The models were verified considering field production data, and a numerical analysis was performed. In the early and later stages of gas production, the gas appeared from the fracture and matrix systems, and it was extracted in different regions in the timescale and synchronously, respectively. Owing to the constant stress boundary condition in the vertical direction of the coal seam, the reduction in the gas pressure in the fractures decreased the horizontal fracture opening, thereby decreasing the permeability in the horizontal direction. Because the horizontal direction exhibited a constant volume boundary condition, the desorption of the adsorbed gas resulted in volumetric shrinkage of the matrix, thereby increasing the permeability in the vertical direction. Non‐Darcy effects reduced the gas flow rate and exerted considerable influence in the early stage of the extraction. Moreover, this effect exhibited anisotropy, which was more pronounced in the horizontal direction. The surface diffusion coefficient continuously increased, which promoted the flow of the adsorbed gas in the matrix. The proposed model can be used to estimate the impact of the permeability anisotropy on the coalbed methane and underground gas extraction.


| INTRODUCTION
Coal and the associated coalbed methane (CBM) are valuable energy resources and have attracted considerable research attention as CBM is a clean and environmentally friendly resource. 1 Furthermore, when the CBM in a reservoir has been exhausted, CO 2 , which is a greenhouse gas, can be injected in the reservoir to realize geological sequestration. 2,3 Therefore, several researchers have conducted in-depth research into the gas flow behavior and its influencing factors. [4][5][6] Permeability represents a key parameter that influences the gas flow capacity, 7,8 and its evolution model can facilitate the realization of linked multifield coupling, thereby attracting widespread research attention. 9 In the early models, the permeability was simply assumed to be a function of the effective stress. 10,11 However, this approach cannot explain the U-shape variation of the permeability change during the gas injection process; in particular, the permeability first decreases and then increases. Many scholars have attributed the complex evolution of the reservoir permeability to the competitive deformation induced by the effective stress and gas adsorption. Palmer and Mansoori 12 proposed the widely adopted P-M model that takes into account the effective stress and matrix shrinkage/swelling. Shi and Durucan 13 proposed a permeability model related to the gas pressure, in which the permeability change was determined considering the competition between the fracture compression and matrix shrinkage terms during the reduction in the gas pressure. However, the aforementioned models are only applicable in the case of uniaxial strain conditions. Robertson and Christiansen 14 proposed a permeability model applicable in the case of biaxial and triaxial stress conditions, in which the permeability was considered as a function of the fracture compression strain, matrix elastic strain, and adsorption strain. Zhang et al 15 proposed a permeability and porosity model considering the pore variation caused by the effective stress and gas adsorption. Based on the work of Zhang et al, Liu et al 16 considered the interaction between the matrix and fractures. In the dual-porosity medium model developed subsequently, the interaction between the fractures and matrix was focused on. 17 In the effective strain model, which has become popular in recent years, 18,19 this interaction is categorized in terms of the gas pressure and gas adsorption effects, both of which are controlled by the mass and effective stress transfer. 18 In summary, the permeability evolution during the gas production process is a complex process and not completely understood.
In general, to develop the permeability model, the following two aspects must be considered. The first aspect pertains to the anisotropy of coal. 20,21 Coal is a typical sedimentary resource, and its flow and mechanical properties vary considerably in the horizontal and vertical directions. 18 For most sedimentary resources, the permeability in the parallel direction of the bedding plane is considerably larger than that in the perpendicular direction. 22,23 At present, several anisotropic permeability models are available. Liu et al 24 investigated the anisotropic permeabilities by considering the directional strains and the reduction ratio of the elastic moduli. Wang et al 25 analyzed the effects of the gas adsorption-induced directional strains and internal swelling on the anisotropic permeability. An et al 26 established an anisotropic permeability model by considering the stress change and gas adsorption. Yang et al 27 proposed an anisotropic permeability model and discussed the effect of the anisotropy on the gas-solid coupling. The difference in the flow performance and anisotropy in the mechanical properties in the horizontal and vertical directions have been extensively examined in the existing reservoir simulations and permeability anisotropy models. 28 However, in an actual reservoir environment, the mechanical boundary conditions vary in different directions. 18 According to recent research, 18 the horizontal and vertical directions reflect a constant volume and constant stress boundary condition, respectively. Nevertheless, the anisotropy in the permeability owing to the different boundary conditions has not been considered in existing work.
In addition to the anisotropy, the internal swelling of the coal must be considered. 29,30 Internal expansion is usually measured considering the internal swelling factor f, 31,32 which indicates the influence of the gas adsorption strain on the permeability. 31,33 Specifically, when coal swells owing to gas adsorption, the coal body expands as a whole, simultaneously producing inward-expanding stress extrusion gaps, which reduce the coal permeability. 32 In the existing models, f is generally assumed to be constant, 34 which can overestimate 34 or underestimate 35 the effect of adsorption on the fracture permeability. In the former case, the gas sorption-induced strain only occurs in the fracture system, 34 and in the latter case, the adsorption leads to only the coal bulk strain. 35 The range of f thus remains controversial. 32,36 The related experiments 37 and theoretical analysis 35 show that f is a process quantity that varies through the gas transfer process; however, this aspect is not considered in most existing permeability models.
Considering these aspects, in this work, f, defined as the ratio of the fracture strain to the bulk strain, was investigated, and a permeability anisotropy model was established based on the different mechanical boundary conditions in the horizontal and vertical directions. Furthermore, a coupled model was developed considering the non-Darcy effect in the fractures and surface diffusion in the matrix system. The model was verified using field data collected from a coal mine, and the sensitivity of the gas production rate to different parameters was analyzed.

| CONCEPTUAL MODEL
CBM reservoirs are usually characterized by an ultralow permeability, and therefore, hydraulic fracturing must be applied to enhance the gas production rate. 38 The gas in a coal seam exhibits both the free and adsorption states; specifically, the gas exists in the free and adsorbed state in the fracture system and fracture wall, respectively. In the matrix, the gas is present in an adsorbed state. The CBM extraction process is shown in Figure 1. Once the production commences, the gas flow in the coal reservoir occurs in two steps. In the first step, the depleted gas can be attributed to the present fractures. The gas rate is usually high because of the high conductivity of the fractures; however, the gas rate reduces rapidly. In the second step, the adsorbed gas in the matrix diffuses to the matrix-fracture interface because of the concentration gradient between two systems.

| Mechanical state of coal reservoir
The stress state of a coal seam during CBM extraction has been extensively researched but remains debatable. 39 CBM wells are usually small, with dimensions of 100-200 mm after drilling; moreover, the stress-influenced range is considerably small compared to the entire reservoir and can thus be neglected. 38 Consequently, most reservoirs can be considered to be in their original stress state. The gas reservoir is assumed to be horizontal and under the action of an invariant overburden stress. 40 Therefore, in gas reservoir simulations, a state involving a uniaxial strain and constant overburden stress is generally assumed. 9,41 According to the assumption of the uniaxial strain, the strain within the horizontal plane is zero, although vertical strain may be present. 42 According to the constant overburden stress assumption, the invariant stress overlies a reservoir and can be calculated considering the weight of the overburden. 9 In this case, the mechanical conditions in the horizontal and vertical directions can be assumed to be constant volume and constant stress conditions, respectively. 43 The permeability in one specific direction is defined as a function of the fracture opening change in the other directions. 44 For example, when the vertical stress decreases, the fracturing opening in the horizontal direction increases, and the horizontal permeability is enhanced. Specifically, the evolution of the horizontal and vertical permeabilities likely follows this trend under a constant stress and constant volume, respectively.

| Internal swelling factor
Coal is a typical porous medium consisting of both matrix and fracture systems. When coal swells owing to gas adsorption, the coal body expands as a whole and produces inwardexpanding stress extrusion gaps, which reduce its permeability. [45][46][47] Herein, three adsorption-induced strains pertaining to the coal bulk (b), matrix (m), and fracture (f) are defined as s b , s m , and s f , respectively. In the initial state, where the subscript 0 represents the initial state. In the initial state at p 0 , the volumes of the coal bulk, matrix, and fracture can be expressed as Similarly at a given p, the corresponding volumes can be expressed as V b0 + V b0 . Equation (1) can be applied under both the initial and current conditions. The following expression indicates the relation among the matrix, fracture, and coal bulk components 48 : where Δ represents the increments.
The strain increment in bulk is a function of the matrix strain and fracture strain, and the porosity can be applied to link the aforementioned strains as follows: The Langmuir curve is used to calculate the adsorption strain 49,50 : where ε S is the gas adsorption-induced strain; ε L is the Langmuir strain; P L (MPa) is the Langmuir pressure constant; and the subscript i = b, m, represents the coal bulk and matrix, respectively.
f is defined as the ratio of the fracture strain to the bulk strain: From Equations (1)-(5), the following expression can be derived 48 : where φ 0 is the initial porosity of the coal, and p is the equilibrium pressure (Pa). It can be noted that f depends on the adsorption parameters of the coal matrix and coal bulk. The former and latter terms can be determined by performing adsorption experiments using coal powder or briquette coal (coal sample made of coal powder) and raw coal, respectively. Usually, ε Lm is larger than ε Lb ; however, the relationship between p Lm and p Lb is unknown. In the following section, the relationship between f and p is derived, with the parameters set as follows: φ 0 = 0.01, ε Lm = 0.04, ε Lb = 0.02, and p 0 = 0.1 MPa. The values of p Lm and p Lb are presented in Figure 2.
As shown in Figure 2, the values of f are less than zero, indicating that the internal swelling leads to a reduction in the permeability. Moreover, the absolute value is significantly larger than the unit value, because, in this study, the pore strain was considered instead of the bulk strain owing to the permeability being directly related to the pore strain. 51 The relative values of p Lm and p Lb considerably influence the evolution of f. Specifically, when p Lm > p Lb and p Lm < p Lb , f increases and decreases with the gas pressure, respectively.

| MATHEMATICAL MODEL
This section describes the formulation of the gas flow and reservoir deformation coupled model, taking into account gas flow in the fracture system, gas diffusion in the matrix system, and permeability in the horizontal and vertical directions.

| Gas flow in the fracture system
The Darcy equation is used to describe the gas flow in the fracture system. The mass conservation law can be expressed as 52 with where subscript f represents the fracture system, p represents the gas pressure (Pa), k denotes the permeability (m 2 ), and u is the velocity vector (Pa s) modified owing to the non-Darcy effect because of the large permeability. 53 The Forchheimer equation is applied to describe the non-Darcy effect as 53 where δ f is a correction coefficient defined as The parameter β f is an intrinsic parameter of the porous media 54 and can be expressed as a power law function of the permeability k f as 55 where γ and η are reservoir-specific coefficients.
ρ is the gas density defined as where p ga is the standard atmospheric pressure (Pa), and ρ ga is the gas density (kg/m 3 ) at this pressure. Both free and adsorbed gas exist in the fracture system: where φ f is the porosity of the fracture system, and ρ c (kg/m 3 ) is the coal density.

| Gas diffusion in the matrix system
The mass transfer between the matrix and fracture system can be defined in terms of the diffusion time as follows 53,56 : where subscript m refers to the matrix system. The diffusion time for the matrix system can be expressed as where D (m 2 /s) and ω denote the diffusion coefficient of the matrix and the shape factor (m −2 ), respectively. It is assumed that only adsorbed gas exists in the matrix system, and the diffusion term corresponds to the surface diffusion that leads to a chemical potential gradient in the adsorption phase and is related to the gas pressure and temperature. Based on the Langmuir adsorption theory, the adsorbed gas mass can be defined as where θ is the surface coverage of the adsorption layer, and P L is the Langmuir pressure. Chen and Yang 57 developed a kinetic model based on the random walk theory, which can be expressed as 56 The adsorbed gas content in equilibrium, m e (p wall ), is determined using the Langmuir isotherm: In addition, the adsorbed gas content in the matrix system is calculated using the Langmuir isotherm:

| Permeability model based on f
Based on the definition of the coal porosity and pore-elastic mechanics, the following equations can be obtained 6,10,58,59 : where = −σ kk /3 is the mean compressive stress (Pa). K is the bulk modulus (Pa), α and β are the Biot coefficients, and p is the equilibrium pressure (Pa).
Furthermore, By assuming that 1 Substituting the definition of f indication in Equation (15) into Equation (24) yields Then, the cubic function 60,61 is applied, and the permeability can be expressed as 48

| Constant confining pressure
As discussed, the mechanical conditions in the horizontal and vertical directions can be considered to be constant volume and constant stress conditions, respectively, as shown in Figure 3.
The confining pressure is invariable, and σ c − σ c0 = 0 ( Figure 3A). Under this condition, the permeability model can be expressed as By substituting f, the permeability model in the horizontal direction can be expressed as

| Constant volume
In the constant volume condition, the total strain is considered to be constant, which means that the increment in the coal bulk strain is zero (dε b = dV b /V b0 = 0). Through a similar formulation as in Equations (20)-(26), the governing equation for the permeability change under constant volume can be written as By substituting f, the permeability model for the vertical direction can be expressed as

| MODEL VERIFICATION
This section describes the verification of the proposed model by using field data collected from a coal mine in China. Specifically, this section introduces the field conditions, describes the model establishment process, and presents the verification results.

| Field conditions
The Gaojiazhuang Coal Mine in Shanxi Province is shown in Figure 4, and it is located in the south of the middle section of the Hedong Coalfield and on the SN-trending tectonic belt on the east bank of the Yellow River and the west slope of the Lvliang Mountain. The mining area is characterized with nearly NS-trending folds (the Lishi syncline and Wangjiahui anticline), NE-or NW-trending, west-dipping fractures (Tuanshuitou Fault, Zhujiadian Fault, etc), and EW-trending Lishi (Shanxi Province)-Wubu (Shaanxi Province) nose-shaped structures. The main coal seams are #2, #4, #8, #9, and #10, which contain coking and lean coal, and have gentle dip angles of 0° to 6°. The characteristics of the coal seam are listed in Table 1. The mine has a production capacity of 3.00 Mt/a with the maximum absolute and relative methane emission rates of 82.55 m 3 /min and 29.04 m 3 /t, respectively, and therefore, this mine is classified as a highly  the working face in the #8 drilling field in the belt roadway of the working face. The boreholes with a diameter of 135 mm were distributed in a fan-shaped arrangement, and different azimuths were set for their construction. The borehole parameters including the borehole lengths are summarized in Table 2. After drilling the boreholes, the extraction effects from the boreholes were continuously recorded using the extraction system.

| Numerical model
The geometric model was established based on the layout shown in Figure 5, with a length, width, and height of 300 m, 200 m, and 6 m, respectively, as shown in Figure 6A. A free tetrahedral was applied to mesh the geometry model with 93 146 elements and a quantity of 0.9, as shown in Figure 6B. The aforementioned coupled gas flow model was applied and solved using COMSOL Multiphysics (version 5.4), which is a commercially available solver for partial differential equations (PDE). For the gas flow in the fracture system, a variable bottom-hole pressure (Equation 31) was applied to the borehole to simulate the extraction pressure, without considering the flow boundary conditions at other boundaries. For the gas diffusion in the matrix, no-flow boundary conditions were used because of the absence of any direct contact with hydraulically formed fractures. The employed parameters are listed in Table 3, wherein the general parameters are those pertaining to the gas properties, and the assumed parameters are those considered in the calculation, as reported in the literature. 28,62,63 In Table 3, the subscript f represents the fracture system, and subscripts x, y, and z denote the permeability in the three mutually orthogonal coordinate directions.
A time-dependent pressure was applied to the borehole to simulate the pressure change during production and realize better convergence 64 : where p b is the final borehole pressure (0.1 MPa), p 0 is the initial gas pressure, and t d is the characteristic time to control the pressure drop.
The presence of the non-Darcy flow renders the model simulation challenging as the correction coefficient δα and velocity vα are interrelated. To address this issue, four subequations must be solved simultaneously: three subequations for the velocity vector in the three directions (Equation 9 and Equation 7), in contrast to the one equation to be solved in the case of Darcy flows. The aforementioned relationship can be solved through the iterations of the four subequation. 56 To this end, the General Form PDE Interface in the COMSOL with four variables (velocity vα in the three directions and the liquid pressure) was used. A desktop computer, with an Intel(R) Core (TM) i7-6700 @ 3.40 GHz as the CPU and 24.0 GB of memory, was used for the computation. The calculation time was approximately 15 minutes. Figure 7 shows the comparison between the numerical fitting results and the results of the gas production in the field. Field data for approximately 1 year, collected from the coal mine, were used for the comparison. The fit goodness of the model results with the field data was 0.9, indicating reasonable agreement. The initial gas flow rate from the field data was approximately 0.45 m 3 /min, which rapidly decreased to 0.1 m 3 /min in approximately 3 months and then remained constant. The numerically simulated rate of gas production in the initial stage was higher than that in situ, likely because the model did not consider the water-gas two-phase flow. The change in the proportion of the gas flow in the fractures and matrix with time is shown in Figure 7. It can be noted that the gas in the early stage comes from the fracture system, whereas it comes from the matrix system in the later stage, wherein it accounts for almost 100% of the gas.

| RESULTS AND DISCUSSION
This section describes the use of the verified numerical model as a benchmark in a numerical analysis conducted to investigate the impact of the permeability evolution, non-Darcy effect, and gas diffusion in the matrix on the gas flow.

| Gas pressure distribution
The numerical model was used to obtain the distribution and evolution of the gas pressure in the fracture system after different durations of gas extraction (Figure 8). After gas extraction for 3 days, the gas pressure in the fracture system in the region near the borehole decreases. An evident pressure reduction zone occurs in a fan-shaped region with a radius of approximately 80 m; however, the gas pressure in the fractures does not decrease significantly beyond the areal coverage of the borehole group. With the increase in the extraction time, the fan-shaped region gradually enlarges and propagates to the bottom of the borehole. As the gas is extracted over 3 months, the fan-shaped region expands to the bottom  of the borehole; in contrast, the gas pressure in fractures does not decrease significantly in the remote triangular region (lower right corner of the pressure distribution), although this value begins to decrease after 3 years. When gas extraction is conducted for 30 years, the gas pressure in the fractures decreases across the entire range of the geometric model. In summary, in the early stage of extraction, the gas pressure in the fractures decreases from the vicinity to the end of the borehole and from the high-density region to the low-density region of the boreholes. In the first 3 months, the gas pressure in the fractures within the areal coverage of the borehole group decreases rapidly, and the gas can be extracted rapidly during this stage. After 3 months, the fractures gradually propagate from the drill site to the far end of the borehole, and the gas pressure in the fractures outside the areal coverage of the borehole group starts to decrease. In particular, in the first 3 months of gas extraction, the gas extracted from the fractures comes from within the areal coverage of the borehole group. After this period, the gas is slowly extracted from the fractures outside this region.
Similarly, the distribution and evolution of the gas pressure in the matrix system after different durations of gas extraction were determined ( Figure 9). As the gas is extracted over 3 days, the gas pressure in the matrix throughout the coal seam is the same as the initial gas pressure. When gas extraction occurs over 3 months, the gas pressure in the matrix changes slightly with an average value of 1.07 MPa because of the low gas transport capacity of the matrix. When the gas is extracted over 3 years, the gas pressure in the matrix decreases: The gas pressure in the matrix within the areal coverage of the borehole group is slightly smaller than that outside the region. The gas pressure in the matrix differs slightly across the whole model. When gas extraction is conducted over 30 years, the gas pressure in the matrix decreases significantly. Specifically, the gas pressure in the matrix during the extraction differs only slightly spatially and exhibits an overall synchronous downward trend. In the early period, the gas pressure in the matrix, as determined by the model, is almost constant, and the gas in the region is not extracted in this stage. After 3 months, the gas pressure in the matrix begins to decrease, and the gas in the matrix is extracted.
The distribution and evolution of the gas pressure in both the fractures and the matrix indicate that the difference in the maximum and minimum gas pressures in the fracture model is larger than that in matrix in the corresponding stages of the gas extraction. In particular, the gas in the fractures is extracted from different regions, whereas the gas in the matrix is synchronously extracted. In the early stage of the gas extraction (first 3 months), the extracted gas is derived from the fractures within the areal coverage of the borehole group. From 3 to 30 years, the extracted gas is derived from the fractures outside the areal coverage and from the matrix.

| Permeability evolution
Based on the developed numerical model, the evolution of the permeability in different directions (horizontal and vertical) was obtained (Figure 10). It can be noted that the permeability in the horizontal direction begins to decrease after the gas extraction. When the extraction time approaches 3 years, k/k 0 is minimized. However, the permeability in the vertical direction begins to increase after the gas extraction and reaches a  Figure 10 shows the evolution of f. It can be noted that the absolute value of f decreases with time during the gas depletion process, which is consistent with the observations from Figure 2.
At the beginning of the gas extraction, the extracted gas is mostly derived from the fractures, and the gas pressure in the fractures begins to decrease. Because of the constant stress boundary condition in the vertical direction, the reduced gas pressure in the fractures leads to a reduction in the fracture volume, and the channels for the gas migration in the horizontal direction become narrow; consequently, the permeability in the horizontal direction decreases. Furthermore, as the gas in the fractures is continuously extracted, because of the constant volume boundary condition in the horizontal direction, the gas pressure in the fractures decreases, and the adsorbed gas in the matrix starts to be desorbed. The continuous desorption and release of gas in the matrix cause the matrix to shrink, thereby increasing the fracture volume and permeability in the vertical direction. When the matrix volume stabilizes, the permeability exhibits its maximum value and stabilizes. Subsequently, the model was used to calculate the gas flow rates without considering the permeability variation in both the horizontal and vertical directions, as shown in Figure 11. It can be noted that when the permeability variation in the horizontal direction is not considered, the gas flow in the early stage of the gas extraction is larger than that calculated using the benchmark model. After 1 year of extraction, the values are consistent with the calculation results obtained using the benchmark model, because the gas in the later stage mainly comes from the matrix system. Furthermore, when the permeability variation in the vertical direction is not considered, the gas flow at the beginning of the gas extraction is smaller than that calculated using the benchmark model, indicating reasonable agreement. The comparative analysis also indicates that influence of the permeability in the horizontal direction is more significant than that of the permeability in the vertical direction. This phenomenon likely occurs because (a) the permeability in the horizontal direction is considerably larger than that in the vertical direction, and (b) the coal reservoir primarily expands in the horizontal direction and is fairly narrow in the vertical direction.

| Non-Darcy effect
The changes in the correction coefficient of the non-Darcy effect in different directions with time were calculated, as shown in Figure 12. It can be noted that the correction coefficient is the smallest in the X-direction where the non-Darcy f ISF,f effect is most obvious, followed by that in the Y-direction, and the correction coefficient in the Z-direction is the largest. In accordance with Equation (8), the correction coefficient of the non-Darcy effect is inversely proportional to the permeability and gas flow rate. The permeabilities and gas flow rates in the X-, Y-, and Z-directions decrease successively, as do the correction coefficients therein. With an increase in the extraction time, the gas flow and gas flow rate constantly decrease and tend to zero; consequently, the correction coefficients in the X-, Y-, and Z-directions tend to 1. The influence of the correction coefficient of the non-Darcy effect on the gas production is illustrated in Figure 13. The gas flows obtained considering and not considering the non-Darcy effect differ significantly only in the early stage of extraction, and they appear to be consistent at other extraction times. The main reason is that the non-Darcy effect primarily pertains to the influence of the gas pressure and gas flow rate in the fractures on the permeability; however, the gas is mainly extraction from the fractures only in the early stage. In this stage, owing to the large gas pressure and gas flow rate in the fractures, the correction coefficient is less than one, which leads to the above-mentioned phenomenon. In the later stage of the extraction, the extracted gas is derived from the matrix. In this stage, the correction coefficient is approximately 1, thereby leading to results that are consistent with those obtained without considering the non-Darcy effect.

| Gas diffusion in the matrix system
The changes in the surface diffusion coefficient (D S ) with the duration of gas extraction when the Langmuir pressure constant −P Lm values were 0.5, 1.5, and 3 MPa are shown in Figure 14. The calculation results (Equations 15 and 16) demonstrate that D S constantly increases with P Lm : These changes can be divided into two stages pertaining to a rapid and gradual increase, which are consistent with the findings of the Langmuir adsorption equation.
The changes in the gas flows with the gas extraction time under P Lm of 0.5, 1.5, and 3 MPa are shown in Figure 15. The calculation results indicate that the gas flows at different P Lm values are consistent in the early stage of the gas extraction and later differ. The main reason is that P Lm characterizes the matrix flow regime. The gas in fractures is mainly extracted in the early stage, whereas the gas extraction in the matrix mainly occurs thereafter. P Lm does not influence the gas flow in the fractures; however, it influences the gas extraction from the matrix after the early stages of the extraction. Moreover, a larger P Lm corresponds to a smaller D S and therefore a smaller gas flow. Therefore, in the later stage of extraction, the gas flow decreases with the increase in P Lm .

| CONCLUSION
The influence of the permeability anisotropy in CBM reservoirs on the gas extraction efficiency was investigated. In contrast to the existing research, the effects of the anisotropy of the mechanical boundary conditions on the permeability and gas production of the coal seam were examined. An anisotropy permeability model was developed considering the mechanical boundary conditions in different directions by using f, and a coupling model was established thereafter. The model was verified based on field production data, and the sensitivity of the gas production rate to different parameters was analyzed. The following conclusions were derived: 1. During gas extraction, the gas is derived from the fracture system in the early stages and is derived from the matrix system thereafter. In terms of the timescale, the gas is extracted from different regions of the fracture system, that is, the gas pressure of the fracture system in the affected and nonaffected regions differs considerably; however, the matrix system corresponds to synchronous extraction, which means that the gas pressure of the matrix system in the affected and nonaffected regions is nearly the same. 2. The permeability evolves differently in the horizontal and vertical directions. Because of the constant stress boundary conditions in the vertical direction, the reduction of the gas pressure in the fractures reduces the fracture volume in the horizontal direction, and the channels for gas migration become more narrow; therefore, the permeability in the horizontal direction decreases. Furthermore, owing to the constant volume boundary conditions in the horizontal direction, the reduced gas pressure in the fractures and the induced gas desorption in the matrix lead to shrinkage of the matrix volume, thereby increasing the fracture volume and permeability in the vertical direction. 3. The non-Darcy effect can reduce the rate of gas flow, albeit only in the early stage of extraction, in which it exhibits anisotropy. The correction coefficient for the non-Darcy effect is minimized in the strike direction, in which the non-Darcy effect is the most obvious, followed by that in the dip direction. This coefficient is maximized in the vertical direction, in which the non-Darcy effect is the weakest. 4. During gas extraction, the surface diffusion coefficient increases constantly, which promotes the flow of the adsorbed gas in the matrix; furthermore, the coefficient increases with a decrease in the Langmuir pressure constant. These changes can be divided into two stages, that is, a rapid and gradual increase stage, which is consistent with the results obtained using the Langmuir adsorption equation.