Abstract
Coal permeability is intrinsically anisotropic because of the cleat structure of coal. Therefore, coal permeability can be denoted by a second-order tensor under three-dimensional conditions. Our previous paper proposed an analytical model of the principal permeability tensor of coal during primary coalbed methane (CBM) recovery. Based on this model, 18 modeling cases were considered in the present study to evaluate how the principal permeabilities were influenced by representative coal properties (the areal porosity, the internal swelling ratio, and the Young modulus) during primary CBM recovery. The modeling results show that with regard to the influences of the areal porosity on the principal permeabilities, an increase in cleat porosity reduces the sensitivity of each principal permeability to pore pressure change. The magnitudes of the principal permeabilities are positively proportional to the internal swelling ratio. The principal permeabilities thus tend to monotonically increase with a depletion in the pore pressure when the internal swelling ratio increases. Because the internal swelling ratio represents the extent of gas-sorption-induced matrix deformation, an increase in the internal swelling ratio increases desorption-induced matrix shrinkage and thus induces an increase in permeability. The principal permeabilities are positively proportional to the isotropic principal Young moduli and the synchronously changing anisotropic principal Young moduli. On the other hand, the principal Young modulus within the plane of isotropy influences the principal permeabilities within this plane in diverse patterns depending on both the dip angle of the coalbed and the pitch angle of the cleat sets. The principal permeability perpendicular to the plane of isotropy is positively proportional to this principal Young modulus, and this correlation pattern is independent of both the dip angle and pitch angle.
1. Introduction
Coal permeability is a property that predominantly determines whether a coalbed methane (CBM) reservoir can be commercially recovered [1–3]. Therefore, the investigation on coal permeability has attracted worldwide attention, and various approaches such as experimental measurements, field observations, mathematical modeling, and numerical simulations have been used for the investigation [4–6].
CBM production triggers a series of coal-gas interactions such as changes in the effective stress and gas desorption from the coal matrix [3, 7, 8]. These interactions influence coal permeability during CBM production. When methane is extracted from a coalbed, the effective stress decreases and cleats are compressed; the compression causes a reduction in permeability. On the other hand, the extraction makes methane desorb from the coal matrix and causes the matrix block to shrink. This shrinkage increases permeability. Therefore, the manner in which coal permeability varies during CBM production is dependent on the net effects of the effective stress and gas desorption [9, 10].
The net effects of the effective stress and gas desorption on the variation pattern of coal permeability are influenced by coal properties such as porosity, Young’s modulus, Poisson’s ratio, and sorption strain. The influences of these properties have been evaluated on the basis of theoretical permeability models. Chen et al. [11] documented the influences of Young’s modulus, Poisson’s ratio, and sorption strain on the variation pattern of coal permeability. Cui and Bustin [12] evaluated how Poisson’s ratio influences the variation pattern. Guo et al. [13] evaluated the influences of sorption strain on the variation pattern. Liu et al. [14] discussed the influences of cleat porosity and sorption strain on the pattern. Liu and Harpalani [15] evaluated how porosity influences the pattern.
The abovementioned studies have considered coal permeability to be isotropic. However, coal permeability is intrinsically anisotropic because the fracture structure of coal is heterogeneously distributed [16]. Coalbed fractures are typically composed of three cleat sets: face cleats, butt cleats, and bedding planes [3, 16]. These three cleat sets are normally orthogonal to each other, and coal permeability is intrinsically anisotropic because of this orthogonal fracture structure [3]. Anisotropic permeabilities can be denoted by a second-order diagonal tensor under three-dimensional (3D) conditions. This tensor can be referred to as the “principal permeability tensor” [17, 18], and the diagonal elements of this tensor can be referred to as “principal permeabilities.” The principal permeabilities may have different responses to coal properties during primary CBM recovery. However, these responses have not been well documented in the literature.
The present study uses 18 modeling cases to evaluate how coal properties influence the principal permeabilities during primary CBM recovery. The modeling cases are considered on the basis of a permeability model developed in one of our previous studies [19].
2. Modeling Influences of Coal Properties on Principal Permeability Tensor during Primary CBM Recovery
2.1. Representations of Principal Permeability Tensor
Because the components of the uniaxial strain conditions are aligned along the vertical and horizontal directions, these components are noncoaxial to the cleat sets when the coalbed is inclined. Two coordinate systems can be formulated to align the components of the uniaxial strain conditions with the cleat sets. One is the uniaxial coordinate system (UCS), which is used to accommodate the uniaxial strain conditions, and the other is the cleat coordinate system (CCS), which is used to accommodate the cleat sets. Figure 1 shows the schematic of the two coordinate systems, which can be correlated to two featured properties. One is the dip angle of the inclined coalbed, and the other is the pitch angle of the cleat sets [19].

Our previous paper [19] proposed a model of the principal permeability tensor during primary CBM recovery. The representations of this model are shown in equation (1), equation (2), and equation (3). The present paper will only summarize the basic assumptions and representations of this model. More details can be obtained from our previous paper. The model was developed for a coalbed that is subjected to uniaxial strain conditions. The model assumes that only elastic deformation occurs in the coalbed. The coalbed contains three orthogonal cleat sets (Figure 1). The orientations of the principal permeability tensor are coaxial to those of the cleat sets. Although stress change may induce fracture initiation or propagation in fractured rocks like coal [20, 21], our proposed model assumes that no fracture initiation or propagation occurs in coal during CBM recovery. These assumptions are seemingly simple, but they are widely adopted in the literature and have been validated against both experimental and field data [4, 22, 23]. where , , and are the principal Young moduli in CCS; is the internal swelling ratio, which represents the ratio of the internal swelling to the known sorption-induced strain of coal block; , , and are the principal permeabilities in CCS; , , and are the initial principal permeabilities in CCS; is the pore pressure; is the initial pore pressure; , , and are the principal Langmuir pressure constants in CCS; , , and are the axes of CCS; , , and are the axes of UCS; , , and are the principal Langmuir strain constants in CCS; is the dip angle of coalbeds; is the pitch angle representing the acute angle between the strike direction of the inclined coalbed and butt cleat; is the Poisson’s ratio; , , and are the normal effective stresses in CCS; , , and are the normal effective stresses in UCS; , , and are the initial principal areal porosities in CCS.
Equation (1) shows how the principal permeabilities vary during gas sorption and effective stress change. The normal effective stresses in the CCS (, , and ) can be represented by the normal effective stresses in the UCS, as expressed in equation (2). The normal effective stresses , , and under uniaxial strain conditions can be represented by equation (3). Note that in equation (3), the principal Young moduli are assumed to be transversely isotropic. The combination of equation (1), equation (2), and equation (3) can represent the variation in the principal permeabilities during primary CBM recovery. In our previous papers [3, 19], the formulations of equation (1), equation (2), and equation (3) have been described and these equations have been validated against experimental and field data. Hence, these details are not repeated in this paper.
2.2. Modeling Configuration
Equation (1), equation (2), and equation (3) show that the principal permeability tensor can be influenced by multiple properties. These properties include the dip angle, the pitch angle, the areal porosity, the internal swelling ratio, and the Young modulus. The influences of both the dip angle and pitch angle have been evaluated in our previous study. This study will use 18 modeling cases (divided into six groups) to evaluate how the areal porosity, the internal swelling ratio, and the Young modulus influence the principal permeabilities during primary CBM recovery. The parameters used in each modeling case are presented in Table 1.
Group 1 evaluates how the isotropic areal porosity influences the principal permeabilities. Group 2 focuses on the influences of the isotropic internal swelling ratio. Because the internal swelling ratio is the rate-controlling parameter of sorption strain, this group represents the influences of gas sorption. Because the Young modulus is anisotropic, the influences of the Young modulus will be discussed in four groups. Group 3 evaluates the influences of the isotropic Young modulus. Group 4 evaluates the influences of the anisotropic and synchronously changing principal the Young moduli. Groups 5 and 6 evaluate the influences of the principal Young moduli and , respectively.
In each modeling case, the principal permeabilities were plotted within a 3D domain formulated by the dip angle, pitch angle, and pore pressure. This domain is referred to as the “3DDPP domain” in this paper. The workflow for the 3D visualization of the principal permeabilities is shown in Figure 2. Equation (1), equation (2), and equation (3) and the related properties were implemented in the 3DDPP domain to compute each principal permeability in this domain. This step was achieved using a self-writing C++ code. The 3DDPP domain was surrounded by three pairs of two-dimensional (2D) boundaries. The computed data were then implemented in the Tecplot 360 Package for 3D visualization.

3. Modeling Results and Discussion
3.1. Influences of Areal Porosity
Figure 3 shows the principal permeabilities within the 3DDPP domain for Cases 1, 2, and 3. For each principal permeability, the structure of the isosurfaces in the 3DDPP domain remains invariant with an increase in the areal porosity. The color legends in Figure 3 show that the increasing areal porosity increases the lower bound but decreases the upper bound of each principal permeability.

Figure 4 provides further details regarding the influences of the areal porosity on the evolution of each principal permeability during primary CBM recovery. When the principal permeability increases monotonically with a depletion in the pore pressure, the acceleration of reduces with an increase in the areal porosity. In this case, is negatively proportional to the areal porosity, as shown by the red curves in Figure 4. When decreases monotonically with a depletion in the pore pressure, the deceleration of reduces with an increase in the areal porosity. In this case, is positively proportional to the areal porosity, as shown by the green curves in Figure 4. When evolves nonmonotonically with a depletion in the pore pressure, the increasing areal porosity reduces the deceleration of when the pore pressure is greater than the rebound pressure but increases the acceleration of when the pore pressure is lower than the rebound pressure. In this case, is always positively proportional to the areal porosity when the pore pressure is higher than the rebound pressure, as indicated by the yellow curves in Figure 4. When the pore pressure is lower than the rebound pressure, can be either positively proportional or negatively proportional to the areal porosity depending on the acceleration of in different cases. Based on the curves shown in Figure 4, the influences of the areal porosity on the principal permeabilities can be summarized as follows: an increase in cleat porosity reduces the sensitivity of each principal permeability to pore pressure change.

The modeling cases described in this subsection are based on the assumption of isotropic cleat porosity. The isotropic nature of cleat porosity is consistent with the experimental observations of Gash et al. [26], who reported that the anisotropy between the cleat porosities in different directions is minor for Fruitland Formation coals. However, the literature also contains reports on cases where the anisotropy of cleat porosity has been found to be substantial. Wang et al. [27] observed that the cleat porosity parallel to the bedding plane could be four times that normal to the bedding plane for coals from the Junggar Basin. When cleat porosity has substantial anisotropy, the principal permeability along the direction of the minimum cleat porosity will have a greater sensitivity to pore pressure or stress change than the principal permeabilities along the directions of the medium and maximum cleat porosities.
3.2. Influences of Internal Swelling Ratio
Equation (1) indicates that all the three principal permeabilities are positively proportional to the internal swelling ratio. The 3D isosurfaces shown in Figure 5 confirm this correlation. The isosurfaces can be classified into two categories for each principal permeability: isosurfaces having values higher than the initial value of each principal permeability, and isosurfaces having values lower than the initial value of each principal permeability. In this study, the first category is referred to as “initial-value greater” (IVG) isosurfaces, whereas the second category is referred to as “initial-value lower” (IVL) isosurfaces. Because the principal permeabilities are positively proportional to the internal swelling ratio, the IVG isosurfaces encroach on the IVL isosurfaces with an increase in the internal swelling ratio. In addition to this encroachment, the increasing internal swelling ratio also increases the maximum and minimum values of each principal permeability, as shown by the color legends in Figure 5.

Because the magnitudes of the principal permeabilities change with the internal swelling ratio, their variation pattern changes according to changes in the pore pressure. The changing variation pattern affects the rebound pressure of the principal permeabilities. The increasing internal swelling ratio changes the distribution of the rebound pressure in the 2D domain formed by the dip angle and pitch angle (2DDP domain), as shown in Figure 6. Figure 6 shows that the rebound pressure increases with an increase in the internal swelling ratio in the entire 2DDP domain for each principal permeability. Because of the increasing rebound pressure, the gray-yellow Domain I tends to encroach on the other two domains. This encroachment indicates that all the three principal permeabilities tend to increase monotonically during pore pressure depletion with an increase in the internal swelling ratio.

The influences of gas sorption on permeability evolution have been widely discussed in the literature [10, 22, 28–32]. Generally, pore pressure depletion induces gas desorption and matrix shrinkage. This process causes cleat dilation and results in increased permeability. On the other hand, an increase in the pore pressure causes gas adsorption and matrix swelling. This process compresses the cleats and reduces permeability. Because the internal swelling ratio represents the extent of gas-sorption-induced matrix deformation, an increase in the internal swelling ratio increases desorption-induced matrix shrinkage and thus induces an increase in permeability.
3.3. Influences of Young’s Modulus
When Young’s modulus is isotropic, the representations of the principal permeabilities are reduced to those expressed in equation (4). This equation set shows that the normal effective stresses in the CCS are independent of the isotropic Young modulus. When the isotropic Young modulus increases, the influences of the effective stress decrease and each principal permeability increases accordingly. All the three principal permeabilities are thus positively proportional to the isotropic Young’s modulus in the entire 3DDPP domain. The 3D isosurfaces shown in Figure 7 confirm this correlation. Figure 7 shows that the IVG isosurfaces gradually encroach on the IVL isosurfaces in the 3DDPP domain for each principal permeability. The maximum and minimum values of each principal permeability also increase with an increase in the isotropic Young modulus, as shown by the color legends in Figure 7. where is the isotropic Young modulus.

Figure 7 also shows that the increasing isotropic Young modulus changes the variation in the principal permeabilities with a depletion in the pore pressure. Generally, as the isotropic Young modulus increases, all the three principal permeabilities tend to increase monotonically when the pore pressure depletes. The rebound pressure also increases with an increase in the isotropic Young modulus, as illustrated in Figure 8. The increasing rebound pressure causes the gray-yellow Domain I to expand for all the three principal permeabilities.

When the Young modulus is anisotropic and the principal Young moduli increase synchronously, the principal Young moduli and can be correlated by a ratio constant (), as expressed in equation (5). The principal permeabilities can thus be expressed by equation (6). This equation set shows that the normal effective stresses in the CCS are independent of both and but are dependent on their ratio . When and increase synchronously, the effects of the effective stress reduce and the principal permeabilities increase accordingly. Therefore, the magnitudes of the principal permeabilities are positively proportional to the synchronously increasing and in the entire 3DDPP domain. The 3D isosurfaces in Figure 9 confirm this correlation.

Because the magnitudes of the principal permeabilities are positively proportional to the synchronously increasing principal Young moduli and within the entire 3DDPP domain, the principal permeabilities tend to increase monotonically during pore pressure depletion, as indicated by the 3D isosurfaces shown in Figure 9. The rebound pressure thus increases accordingly for each principal permeability, as illustrated in Figure 10. The increasing rebound pressure causes the gray-yellow Domain I to expand for all the three principal permeabilities.

When increases and remains invariant, the normal effective stresses in the CCS are dependent on both and . In this case, the influences of on the principal permeabilities cannot be seen visually from the representation of each principal permeability. Figure 11 shows the 3D isosurfaces for the principal permeabilities within the 3DDPP domain for Cases 13, 14, and 15. Figure 11 shows how the isosurface distribution changes with an increase in for each principal permeability. The IVG isosurfaces gradually encroach on the IVL isosurfaces with an increase in for . The maximum and minimum values of are positively proportional to . This correlation indicates that is positively proportional to within the entire 3DDPP domain.

However, the influences of on and seem to present diverse patterns depending on both the dip angle and pitch angle. First, the IVG and IVL isosurfaces encroach on each other with an increase in for and . The IVL isosurfaces encroach on the IVG isosurfaces along the lower pitch angle boundary, whereas the IVG isosurfaces invade the IVL isosurfaces along the upper dip angle boundary for , as illustrated in Figures 11(a), 11(d), and 11(g). Similarly, for , the IVG isosurfaces encroach on the IVL isosurfaces along the upper dip angle boundary, whereas the IVL isosurfaces invade the IVG isosurfaces along the lower pitch angle boundary, as shown in Figures 11(b), 11(e), and 11(h). Second, the maximum and minimum values of both and are not positively proportional to . When increases from 500 to 1700 MPa, the minimum values of and increase, whereas their maximum values decrease. When increases from 1700 to 3000 MPa, the minimum values of and decrease, whereas their maximum values do not change. Finally, the normalized curves of are plotted (Figure 12), showing the diverse patterns formed by the influences of with a variation in the dip angle and pitch angle. The curves in Figure 12 indicate that the patterns showing the influences of on are independent of the pore pressure. The red curves show that increases monotonically with an increase in . The green curves show that decreases monotonically with an increase in . The yellow curves indicate that evolves nonmonotonically with an increase in .

Because the patterns showing the influences of on and are dependent on the dip angle and pitch angle, the 2DDP domain can be divided according to these diverse patterns, as shown in Figure 13. Figures 13(a) and 13(b) indicate that and evolve differently, showing diverse patterns, with an increase in . These patterns include a monotonically increasing pattern, nonmonotonically evolving pattern with a maximum point, nonmonotonically evolving pattern with a minimum point, and monotonically decreasing pattern. Figure 13(c) shows that is positively proportional to , and this pattern is independent of both the dip angle and pitch angle.

(a)

(b)

(c)
The variation in the principal permeabilities with a depletion in the pore pressure also shows a complex dependence on because of the effects of the dip angle and pitch angle. Figure 14 shows how the rebound pressure varies with an increase in . The rebound pressures of and are dependent on in diverse ways. Domains I and II encroach on each other with an increase in for both and . Domain I encroaches on Domain II along the right edge of the 2DDP domain. In this encroached domain, the and curves change from nonmonotonically evolving ones to monotonically increasing ones with an increase in . Domain II encroaches on Domain I along the bottom edge of the 2DDP domain for and along the top edge of the 2DDP domain for . In this encroached domain, the and curves change from monotonically increasing ones to nonmonotonically evolving ones with an increase in . The images in the right column of Figure 14 show that the rebound pressure of is positively proportional to in the entire 2DDP domain. Accordingly, Domain I always tends to encroach on Domains II and III with an increase in for . This encroachment indicates that the curves tend to become positively proportional to within the 2DDP domain.

When increases and remains invariant, the normal effective stresses in the CCS are dependent on both and . In this case, the influences of on the principal permeabilities cannot be seen visually from their representations. Figure 15 shows the principal permeabilities within the 3DDPP domain for Cases 16, 17, and 18. The IVG isosurfaces gradually encroach on the IVL isosurfaces with an increase in for and . This indicates that the magnitudes of and seem to be positively proportional to within the entire 3DDPP domain. However, for , the IVL isosurfaces encroach on the IVG isosurfaces with an increase in . This indicates that the magnitudes of seem to be negatively proportional to within the 3DDPP domain.

Figure 16 shows that the rebound pressures of and are positively proportional to . Accordingly, Domain I tends to encroach on Domains II and III with an increase in . This encroachment indicates that the and curves tend to become positively proportional to . The images in the right column of Figure 16 show that the rebound pressure of is negatively proportional to in the entire 2DDP domain. Accordingly, Domain II tends to encroach on Domain I with an increase in . This indicates that the curves tend to become negatively proportional to .

A comparison of the modeling results for Groups 6 and 7 indicates that an increase in or causes an increase in their respective normal principal permeabilities. However, the influences of or on their respective parallel principal permeabilities are distinct. The variation in and shows diverse patterns with a depletion in the pore pressure and an increase in , and this variation depends on both the dip angle and pitch angle. The principal permeability is negatively proportional to , and this pattern of correlation is independent of both the dip angle and pitch angle.
4. Conclusions
In this study, 18 modeling cases were used to evaluate how the variation patterns of principal permeabilities during primary CBM recovery were influenced by three representative coal properties: the areal porosity, the internal swelling ratio, and Young’s modulus. Based on the study results, the following conclusions are drawn: (1)The correlation between the principal permeabilities and the areal porosities is complex and depends on the dependence of the principal permeabilities on pore pressure. The influences of the areal porosity on the principal permeabilities can be summarized as follows: an increase in cleat porosity reduces the sensitivity of each principal permeability to pore pressure change(2)The principal permeabilities thus tend to increase monotonically with a depletion in the pore pressure when the internal swelling ratio increases. Because the internal swelling ratio represents the extent of gas-sorption-induced matrix deformation, an increase in the internal swelling ratio increases desorption-induced matrix shrinkage and thus induces an increase in permeability(3)The principal permeabilities are positively proportional to the isotropic principal Young moduli and the synchronously changing anisotropic principal Young moduli. An increase in the principal Young modulus or causes an increase in their respective normal principal permeabilities. However, the influences of or on their respective parallel principal permeabilities are distinct
Data Availability
There is no experimental data in this manuscript.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge the support of the National Natural Science Foundation of China (51804312, 51704184, and 51874314), the State Key Laboratory Cultivation Base for Gas Geology and Gas Control (WS2019A03), and the Fundamental Research Funds for the Central Universities (2021YQAQ01).