Influences of Coal Properties on the Principal Permeability Tensor during Primary Coalbed Methane Recovery: A Parametric Study

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 desorptioninduced 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.


Introduction
Coal permeability is a property that predominantly determines whether a coalbed methane (CBM) reservoir can be commercially recovered [1][2][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][5][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].

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].

Geofluids
Δσ eXX = Δσ exx cos 2 θ d +Δσ ezz sin 2 θ d À Á cos 2 θ p +Δσ eyy sin 2 θ p , Δσ eYY = Δσ exx cos 2 θ d +Δσ ezz sin 2 θ d À Á sin 2 θ p +Δσ eyy cos 2 θ p , where E X , E Y , and E Z are the principal Young moduli in CCS; F in is the internal swelling ratio, which represents the ratio of the internal swelling to the known sorptioninduced strain of coal block; k X , k Y , and k Z are the principal permeabilities in CCS; k X0 , k Y0 , and k Z0 are the initial principal permeabilities in CCS; p is the pore pressure; p 0 is the initial pore pressure; p LX , p LY , and p LZ are the principal Langmuir pressure constants in CCS; X, Y, and Z are the axes of CCS; x, y, and z are the axes of UCS; ε LX , ε LY , and ε LZ are the principal Langmuir strain constants in CCS; θ d is the dip angle of coalbeds; θ p is the pitch angle representing the acute angle between the strike direction of the inclined coalbed and butt cleat; ν is the Poisson's ratio; σ eXX , σ eYY , and σ eZZ are the normal effective stresses in CCS; σ exx , σ eyy , and σ ezz are the normal effective stresses in UCS; ϕ aX0 , ϕ aY0 , and ϕ aZ0 are the initial principal areal porosities in CCS. Principal permeability k X Principal permeability k Y Principal permeability k Z Internal swelling ratio Figure 6: 2D contours showing distribution of rebound pressure within the 2DDP domain for Cases 4, 5, and 6. Domain I represents the region where the rebound pressure is greater than 5.0 MPa and the principal permeabilities increase monotonically with a depletion in the pore pressure. Domain II represents the region where the rebound pressure is 0.0-5.0 MPa and the principal permeabilities vary nonmonotonically with a depletion in the pore pressure. Domain III represents the region where the rebound pressure is lower than 0.0 MPa and the principal permeabilities decrease monotonically with a depletion in the pore pressure. 6 Geofluids Equation (1) shows how the principal permeabilities vary during gas sorption and effective stress change. The normal effective stresses in the CCS (Δσ eXX , Δσ eYY , and Δσ eZZ ) can be represented by the normal effective stresses in the UCS, as expressed in equation (2). The normal effective stresses Δσ exx , Δσ eyy , and Δσ exx 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. (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.

Modeling Configuration. Equation
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 E X and E Z , 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.

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   III   III   III   II   II   II   II   II  II   I   I   I  I  I  II   I   I    Principal permeability k Y Principal permeability k Z Figure 9: 3D isosurfaces showing principal permeabilities within the 3DDPP domain for Cases 10, 11, and 12. The color legends in these images have a uniform measurement unit of "millidarcy (mD)." 8 Geofluids 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 k Z increases monotonically with a depletion in the pore pressure, the acceleration of k Z reduces with an increase in the areal porosity. In this case, k Z is negatively proportional to the areal porosity, as shown by the red curves in Figure 4. When k Z decreases monotonically with a depletion in the pore pressure, the deceleration of k Z reduces with an increase in the areal porosity. In this case, k Z is positively proportional to the areal porosity, as shown by the green curves in Figure 4. When k Z evolves nonmonotonically with a depletion in the pore pressure, the increasing areal porosity reduces the deceleration of k Z when the pore pressure is greater than the rebound pressure but increases the acceleration of k Z when the pore pressure is lower than the rebound pressure. In this case, k Z 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, k Z can be either positively proportional or negatively proportional to the areal porosity depending on the acceleration of k Z 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.

Geofluids
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. (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 "initialvalue lower" (IVL) isosurfaces. Because the principal permeabilities are positively proportional to the internal swell-ing 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.

Influences of Internal Swelling Ratio. Equation
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][29][30][31][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    Figure 13: Division of the 2DDP domain in terms of variation patterns of principal permeabilities with respect to increasing E X . In Domain i, the principal permeabilities increase monotonically with an increase in E X . In Domain ii, the principal permeabilities vary nonmonotonically with an increase in E X , and the curve of each principal permeability has one maximum point. In Domain iii, the principal permeabilities evolve nonmonotonically with an increase in E X , and the curve of each principal permeability has one minimum point. In Domain iv, the principal permeabilities decrease monotonically with an increase in E X .

Geofluids
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 E 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 E X and E Z can be correlated by a ratio constant (R c ), 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 E X and E Z but are dependent on their ratio R c . When E X and E Z 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 E X and E Z 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 E X and E Z 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 E X increases and E Z remains invariant, the normal effective stresses in the CCS are dependent on both E X and E Z . In this case, the influences of E X 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 E X for each principal permeability. The IVG isosurfaces gradually encroach on the IVL isosurfaces with an increase in E X for k Z . The maximum and minimum values of k Z are positively proportional to E X . This correlation indicates that k Z is positively proportional to E X within the entire 3DDPP domain.
However, the influences of E X on k X and k Y 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 E X for k X and k Y . 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 k X , as illustrated in 12 Geofluids Figures 11(a), 11(d), and 11(g). Similarly, for k Y , 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 k X and k Y are not positively proportional to E X . When E X increases from 500 to 1700 MPa, the minimum values of k X and k Y increase, whereas their maximum values decrease. When E X increases from 1700 to 3000 MPa, the minimum values of k X and k Y decrease, whereas their maximum values do not change. Finally, the normalized curves of k X are plotted (Figure 12), showing the diverse patterns formed by the influences of E X with a variation in the dip angle and pitch angle. The curves in Figure 12 indicate that the patterns showing the influences of E X on k X are independent of the pore pressure. The red curves show that k X increases monotonically with an increase in E X . The green curves show that k X decreases monotonically with an increase in E X . The yellow curves indicate that k X evolves nonmonotonically with an increase in E X . Because the patterns showing the influences of E X on k X and k Y 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 k X and k Y evolve differently, showing diverse patterns, with an increase in E X . 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 k Z is positively proportional to E X , and this pattern is independent of both the dip angle and pitch angle.

Geofluids
The variation in the principal permeabilities with a depletion in the pore pressure also shows a complex dependence on E X because of the effects of the dip angle and pitch angle. Figure 14 shows how the rebound pressure varies with an increase in E X . The rebound pressures of k X and k Y are dependent on E X in diverse ways. Domains I and II encroach on each other with an increase in E X for both k X and k Y . Domain I encroaches on Domain II along the right edge of the 2DDP domain. In this encroached domain, the k X -p and k Y -p curves change from nonmonotonically evolving ones to monotonically increasing ones with an increase in E X . Domain II encroaches on Domain I along the bottom edge of the 2DDP domain for k X and along the top edge of the 2DDP domain for k Y . In this encroached domain, the k X -p and k Y -p curves change from monotonically increasing ones to nonmonotonically evolving ones with an increase in E X . The images in the right column of Figure 14 show that the rebound pressure of k Z is positively proportional to E X in the entire 2DDP domain. Accordingly, Domain I always tends to encroach on Domains II and III with an increase in E X for k Z . This encroachment indicates that the k Z -p curves tend to become positively proportional to E X within the 2DDP domain.
When E Z increases and E X remains invariant, the normal effective stresses in the CCS are dependent on both E X and E Z . In this case, the influences of E Z 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 E Z for k X and k Y . This indicates that the magnitudes of k X and k Y seem to be positively proportional to E Z within the entire 3DDPP domain. However, for k Z , the IVL isosurfaces encroach on the IVG isosurfaces with an increase in E Z . This indicates that the magnitudes of k Z seem to be negatively proportional to E Z within the 3DDPP domain. Figure 16 shows that the rebound pressures of k X and k Y are positively proportional to E Z . Accordingly, Domain I tends to encroach on Domains II and III with an increase in E Z . This encroachment indicates that the k X -p and k Y -p curves tend to become positively proportional to E Z . The images in the right column of Figure 16 show that the rebound pressure of k Z is negatively proportional to E Z in the entire 2DDP domain. Accordingly, Domain II tends to encroach on Domain I with an increase in E Z . This indicates that the k Z -p curves tend to become negatively proportional to E Z .
A comparison of the modeling results for Groups 6 and 7 indicates that an increase in E X or E Z causes an increase in their respective normal principal permeabilities. However, the influences of E X or E Z on their respective parallel principal permeabilities are distinct. The variation in k X and k Y shows diverse patterns with a depletion in the pore pressure and an increase in E X , and this variation depends on both the dip angle and pitch angle. The principal permeability k Z is negatively proportional to E Z , and this pattern of correlation is independent of both the dip angle and pitch angle.

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 E X or E Z causes an increase in their respective normal principal permeabilities. However, the influences of E X or E Z on their respective parallel principal permeabilities are distinct

Data Availability
There is no experimental data in this manuscript.