Numerical Investigation of an Anisotropic Permeability Model for Bedded Coal Based on the Equivalent Fracture Aperture Coefficient

Permeability evaluation of bedded coal is always a difficult problem, thereby an improved permeability model is established based on the equivalent fracture aperture coefficient. The equivalent fracture aperture coefficient can reflect the anisotropy of coal in mechanical and seepage properties. This model is executed into COMSOL Multiphysics to perform numerical simulations. The evolutions of permeability, gas pressure, and cumulative gas production on three typical boundary conditions are investigated. Our simulation results show that this model is competent to represent evolution of anisotropic permeability in bedded coal during gas production, which is governed by the competitive effects of desorption and effective stress. Under the same boundary conditions, different equivalent fracture aperture coefficients lead to significantly different trends in anisotropic permeability evolution. Anisotropic permeability changes with gas pressure and time during production, and the magnitude of change depends on both the equivalent fracture aperture coefficient and applied boundary conditions. Our model shows that the bedding characteristics of permeability have a significant effect on gas production forecasts.


Introduction
Coalbed methane (CBM) [1] has received increasing attention in China in recent years, but the scale of CBM commercialization in China is still too small compared with other countries (such as the US and Australia). e reason is that the reservoir conditions in China are relatively complex, especially many reservoirs have obvious bedding structure. In the process of CBM production, coal deformation and gas migration is a complex coupling process involving many factors, among which permeability is a key factor [2]. e typical bedding structure of coalbed methane reservoirs in China will lead to the anisotropy of mechanical and permeability characteristics, so a better permeability model needs to be advanced to re ect the characteristics of this bedding structure for predicting CBM production.
Numerous studies have shown that permeability evolution during gas migration is controlled by many factors such as coal structure, sorptive capacity, e ective stress, and temperature . For the in uence of e ective stress, these studies generally show that increase of permeability is accompanied by decrease of e ective stress and increase of gas pressure. For bedding structure, the permeability of coal seams with horizontal beddings is generally higher than that with bedding planes perpendicular to the horizontal direction. Some scholars have paid attention to the combined e ects of multiple factors [23]. For example, Seidle and Huitt [25] studied permeability evolution under combined e ect of sorptive e ect and stress, assuming that coal is isotropic to simplify the problem. However, coal is a typical sedimentary rock in nature, and the existence of bedding planes generates isotropic mechanical properties and corresponding anisotropic permeability. Pomeroy and Robinson [26] rst found that the permeability of coal was di erent along di erent bedding directions. Koenig and Stubbs [27] found an order of magnitude di erence in permeability along bedding direction. Gash et al. [28] claimed that the presence of bedding and joints led to directional compression or displacement of coal, resulting in anisotropic permeability. In our previous studies, we conducted anisotropic permeability tests and found that the sensitivity of permeability to stress is similar to that of previous scholars. In addition, we identified an exponential relationship between permeability and the bedding angle [29].
Despite these aforementioned experimental studies on the coal anisotropic permeability, how to quantify the influence of the bedding on coal permeability remains largely unsolved. erefore, in the following work, we propose a permeability evolution model that can reflect the bedding characteristics and predict gas production.

Permeability Evolution Model for Bedded Coal
Korsnes et al. [30] and Wang et al. [31] have shown that the susceptibility of temperature and water on permeability assessment cannot be neglected. When exposed to water and varying temperature condition, effective porosity can be reduced and swelling of coal matrix can be induced. To simplify the problem, we assume an as-received water content and isothermal condition.

Isotropic Permeability Model.
Many permeability models considering fractured coal consist of matrix and cheats, as shown in Figure 1. Isotropic deformation is typically assumed in the matrix and anisotropic deformation in orthogonal fractures due to sorptive and mechanical behaviors. In such a discontinuous medium, the key to understanding permeability evolution is to quantify the change of fracture aperture. Liu et al. [32] indirectly obtained the quantitative relationship for fracture aperture by subtracting the individual matrix deformation from the total deformation of both the matrix and the cheats. Coal porosity can be defined by a uniform cheat spacing s and cleat width b as where ϕ is the porosity of the coal. According to Liu et al. [32], effective stress affects both the matrix and the fractures. From Hooke's law, the elastic modulus is needed to solve the deformation. e elastic modulus of matrix in this model is denoted by E m , but the elastic modulus of cheats is difficult to obtain. In view of this situation, elastic modulus reduction ratio as R m � E/E m is introduced.
ereby, by assuming that the discontinuous medium as a whole has an elastic modulus of E, we can indirectly obtain the fracture aperture change Δb by subtracting the matrix deformation from the total deformation of the entire medium, which can be written as Because b << s, equation (2) can be simplified to the effective strain increment can then be defined as Δε e � Δσ e /E, and equation (3) can be rewritten as where b 0 represents the initial fracture aperture of the coal. From equation (1), the initial coal porosity can be calculated by Many field and laboratory experiments [8,9,[12][13][14] have shown that mechanical behavior of coal is affected by two types of competitive behaviors as follows: (i) Desorption of gases such as CH 4 and CO 2 reduces the gas pressure and increases the effective stress, resulting in coal compression and fracture closure (ii) Gas desorption in the matrix causes the matrix to expand and the cheats to compress erefore, when considering sorptive gas, the total effective strain change Δε e can be expressed as where Δε t represents the whole strain increment and Δε s represents shirking/swelling strain increment. Substituting equations (5) and (6) into (4) yields Many experimental studies [10,33,34] point that there is a cubic relationship between porosity and permeability.  Shock and Vibration erefore, combining the cubic relationship with (1), we obtain the isotropic permeability equation as where k is the dynamic permeability and k 0 is the initial coal permeability. Substituting equation (7) into (8), equation (8) can be simplified to 2.2. Equivalent Fracture Aperture. Many simplified models have been proposed to clarify mechanical and transport properties of raw coal, of which the Warren-Root model [35] is one of the most used. is model is composed of block units consisting of a matrix system as the gas storage source and fractures providing gas flow. However, gas production predicted by this type of models often shows a large deviation compared with field data. e reason could be that this simplified structure of coal does not realistically reflect the actual bedding structure of coal. While the Warren-Root model is applicable to coal seams with bedding direction parallel or perpendicular to the horizontal plane, it cannot effectively describe coal seams with bedding direction oblique to the horizontal plane. Yang et al. [29] identified an exponential function based on natural constant e to quantity permeability and bedding direction and found that that the fracture aperture is related to the bedding direction. erefore, we propose an improved Warren-Root model considering the bedding structure to evaluate mechanical and migration characteristics of bedded coal. In this simplified coal model shown in Figure 2, the original fracture aperture in different bedding directions is b θ , while an equivalent fracture aperture b e is proposed to replace the original fracture aperture.
Here, we introduce an equivalent fracture aperture coefficient η, based on which the fracture aperture in any bedding direction may be written as follows: When the bedding direction is parallel to the horizontal plane, i.e., θ � 0°, η � 1, the value range is 0 < η ≤ 1. Understandably, replacing the original fracture aperture with the equivalent fracture aperture will inescapably cause differences of mechanical properties, which are mainly reflected in the modulus of coal and are also claimed by other scholars [36,37]. In their experiments, authors found that modulus of coal is very different depending on whether the bedding direction is parallel or perpendicular to the horizontal plane, and experiments displayed that coal's elastic modulus exhibits a trend close to an exponential function as the bedding direction relative to the sample changes. erefore, we introduce a modulus correction coefficient δ for bedded coal to obtain the equivalent elastic modulus, expressed as where E is coal's elastic modulus with bedding direction parallel to the horizontal plane and E θ is coal's elastic modulus in different bedding directions.

Improved Anisotropic Permeability
Model. From Section 2.2, substituting equations (10) and (11) into (9), we obtain an improved permeability model based on the equivalent fracture aperture coefficient, which is written as Subsequently, the change of permeability takes the form of which can be further simplified as where k 0 is the coal permeability at a certain reservoir pressure as θ � 0°.

Model
Validation. e improved permeability model based on the equivalent fracture aperture coefficient proposed in this study is expected to reflect the influence of bedding planes on mechanical and transport characteristics of coal. Next, we use the laboratory experimental data conducted on raw coal samples under isothermal conditions to verify the fidelity of the model [29]. It should be noted that after the sample was recovered from the coal seam, deformation had occurred due to the release of the in situ stress. Certainly, it is hard to reappear coal permeability under the real original in situ stress state. However, our experimental results demonstrated that the internal structure of the coal samples tends to stabilize during multiple loading and unloading cycles under different hydrostatic stresses, which to some degree avoided the test error caused by the compressive damage. erefore, as a reference, the permeability data from the third and fourth loading cycles are used to verify the improved model.
Many experimental studies have shown the dependence between coal permeability and effective stress. For example, Seidle and Huitt [26] and Yang et al. [29] experimentally demonstrated that gas permeability gradually decreases with increasing stress. For instance, In the HN-01 sample [29], permeability was reduced from 6e −18 m 2 to 0.7e −18 m 2 , as the confining stress increased from 3 MPa to 17 MPa. It is worth mentioning that this behavior occurs in multiple loading and unloading stress cycles. Moreover, permeability can also exhibit different anisotropic characteristics based on the angle between bedding direction and load direction, as shown in Figures 3 and 4.
As shown in Figure 3, the quantitative relationship between the bedding direction and permeability can be identified by our new model proposed in this study. Specifically, we assume that both the equivalent fracture aperture coefficient and modulus correction coefficient are 1 when the bedding angle is horizontal (θ � 0°). Our model is very close to the experimental data with increase in effective stress and decrease in permeability. e equivalent fracture aperture coefficient is 0.82 when the bedding angle is 30°, which means that the permeability is reduced due to the equivalent fracture aperture and the porosity reduction. It should be noted that equivalent elastic modulus is changed with the bedding direction as well, and the corresponding modulus correction coefficient is 1.22. By considering the influence of structural and mechanical parameters, this improved model is also representative of describing permeable change caused by effective stress. Similar fitting curve is also obtained when the bedding direction is perpendicular to the horizontal plane (θ � 90°). As illustrated in Figure 4, similar trends can be seen from the fourth loading period as well. Importantly, this model is able to reasonably match all experimental data during several loading and unloading cycles for the coal samples with different bedding angles.
To further validate our model, we also consider experiments conducted by other scholars, particularly, the permeability tests on coal samples with different bedding angles by Pan et al. [13]. In their tests, the greater effective stress and bedding angle, the smaller permeability, consistent with permeability evolution observed in Yang et al. [29]. Results in Figure 5 show that when the equivalent fracture aperture coefficient is 0.85 and the modulus correction coefficient is 1.16, the model has the best fitting curve for the coal permeability results of test X1 with oblique bedding planes.
Similarly, when the equivalent fracture aperture coefficient is 0.57 and the modulus correction coefficient is 1.1, the model can best match the coal permeability results of test P1 with vertical bedding planes.

Coupled Gas-Solid Model for Bedded Coal
On the desorption process, permeability is not only a single variable but also closely related to coal deformation and gas flow, which need to be fully coupled and examined simultaneously. To simplify the analysis, the following assumptions are considered in this work: (i) Coal is a transversely isotropic continuous medium (ii) Strain is infinitesimal and linear elastic (iii) e desorption process is isothermal (iv) Gas migration is laminar flow and follows Darcy's law

Control Equation of Coal Deformation.
Based on theory of linear elastic porous media [38], defined the effective stress σ e ij , the total strain ε ij , and the sorptive strain ε s , the bedded coal constitutive model can be written as where δ ij is the Kronecker symbol and C ijkl is the matrix stiffness. We specify that the compression is negative. e effective stress σ e ij can be written as Li et al. [10].
where σ ij is effective stress tensor and p is gas pressure in matrix and cheats. e sorptive strain ε s can be written as Harpalani and Chen [39].
where the sorptive strain follows the Langmuir type and ε l and p L are sorption constants.  Shock and Vibration e differential equilibrium equations of continuum can be written as where f i is body force in i direction and x j is j th coordinate. According to geometric equations ε ij � (u j,i + u i,j )/2, taking equations (15)-(17) into (18), the control equation of coal deformation can be written as follows (see Liu et al. [40,41]): Equation (19) represents the basic control equation for the coal's deformation. e characteristics of bedding structure are mainly embodied in the change of elastic modulus of coal seam at different bedding angles. [42] claimed that coalbed methane exists in coal seams in two forms: free gas and adsorbed gas. Cheats are filled with free gas as migration channel and matrices are filled with sorptive gas as a container. So, the gas density C can be expressed as

Control Equation of Gas Flow. George and Barakat
where ρϕ is the free gas in cheats and (1-ϕ)ρ g ρ s [V Lp /(p L + p)] is the sorptive gas in coal matrix. ρ g and ρ s represent the gas density at atmospheric pressure and the skeletal density of coal, respectively. According to gaseous equation of state, density of gas is rewritten as pressure as follows: Substituting equation (21) into (20) yields where M is the ideal gas molar mass, R is the gas constant, and T is the temperature of the coal reservoir. Gas migration is laminar flow and follows Darcy's law, and it can be expressed as where μ is the gas viscosity. Finally, the gas flow governing equation can be written as Equation (24) represents the basic control equation for gas migration. k θ is the key parameter reflecting seepage characteristics of coal seams with different bedding angles.   Figure 5: Coal permeability as a function of effective stress. Experimental data are from Pan et al. [13]. e coupled gas-solid equations for coal contain equations (14), (19), and (24), which will be implemented in the COMSOL Multiphysics for numerical simulations in the next section.

Model Applications
In this section, we use the gas-solid coupled model to predict field production. As illustrated in Figure 6, a two-dimensional (2D) model with a size of 4 m × 4 m is built in COMSOL Multiphysics and a small hole with a radius of 0.1 m representing a well is located in the middle (top figure in Figure 6). e model boundary conditions are shown in Figure 6. Noflow boundary conditions are applied on all model boundaries, and the atmospheric pressure is applied to the well. e initial gas pressure of the model is 2 MPa. e monitoring point for properties and parameters (e.g., gas pressure and permeability) is located at C.

Shock and Vibration
Typically, there are three types of mechanical boundary conditions [14]: constant stress condition, constant volume condition, and uniaxial strain condition. However, these three types of boundary conditions are not universally applicable to any stress environment in the field; therefore, it is important to choose the appropriate boundary conditions based on the actual situation of the field. In this study, different boundary conditions are used to quantify the effects of coal structure, gas adsorption, effective stress, and other factors on the evolution of permeability. e parameters of this model are listed in Table 1 [2]. Figure 7, constant stress boundary condition is applied to the top and right boundaries of this model.

Constant Stress Condition. As shown in
It can be observed that when gas pressure is greater than 1.2 MPa, increase of permeability is accompanied by decrease of effective stress; when gas pressure is 1.2 MPa, permeability reaches the lowest value; when gas pressure is less than 1.2 MPa, permeability rebounds to increase with increasing effective stress. e explanation of permeability rebounding is that effective stress effect and sorptive effect compete with each other during gas production process. In    Shock and Vibration the initial stage of production, gas pressure in the coal seam is high, change behavior of effective stress on permeability is stronger than that of desorption behavior; in other words, effective stress is dominant in permeability evolution. Conversely, change behavior of effective stress on permeability is weaker than that of desorption at low gas pressure during production, where gas desorption causes matrix shrinkage and exceeds change of fracture apertures. As a result, the influence of desorption on permeability is dominant at low gas pressure. is rebounding behavior is consistent with earlier studies, such as Wang et al. [43]. For the next case, the equivalent fracture aperture coefficient is 0.8 and the modulus correction coefficient is 1.2, which means that the bedding direction in the model is oblique to the horizontal plane. As shown in the blue line in Figure 8(a), with reduction of gas pressure and increment of effective stress, permeability ratio of coal increases from 0.5 to nearly 1. e reason is that the introduction of equivalent fracture aperture coefficient and modulus correction coefficient leads to the increase in the equivalent elastic modulus of coal. At the same effective stress, the impact of deformation associated with the permeability is smaller. erefore, the influence of desorption has a dominant effect on permeability evolution.
Notice that the smaller the equivalent fracture aperture coefficient is, the more perpendicular the bedding orientation is to the horizontal plane. Concurrently, the less permeable the coal is, the slower the gas pressure changes. For instance, as shown in Figure 9(a), at production time t � 100D, the gas pressure of the model with bedding parallel to the horizontal plane drops to 0.28 MPa (black line); the gas pressure of the model with bedding oblique to the horizontal plane drops to  ). is is because, among the three cases, the last one has the smallest permeable path, which slows the gas desorption rate and the gas pressure decline rate. Figure 9(b) shows the comparison of cumulative gas production volume for these three cases. We find that gas production reaches the peak at approximately 200 days when both the equivalent fracture aperture coefficient and the modulus correction coefficient are 1 while gas production reaches the peaks at about 250 days and 400 days for the other two cases, respectively (θ � 45°and θ � 90°). Figure 10, the uniaxial strain condition is a state between the constant stress and the constant volume boundary conditions, where all vertical boundaries have zero displacement, and the stresses on the upper and lower boundaries remain constant [6,9]. e effective stress only affects the permeability in vertical direction of the model, while desorption affects permeability in all directions under this condition. It can be predicted that desorption effect of the model would dominate evolution of permeability in this case. During gas desorption, the permeability increases with decreasing gas pressure due to the dominant effect of coal matrix shrinkage on the net change of fracture aperture. Figure 11 shows that coal permeability ratio increases from 1.3 to 6 with decreasing gas pressure and increasing effect stress (black line) when both the equivalent fracture aperture coefficient η and the modulus correction coefficient δ are 1, which means that the bedding orientation is parallel to the horizontal plane. When η � 0.8 and δ � 1.2, coal permeability ratio increases from 1.5 to 2.5 during gas production. Lastly, when η � 0.6 and δ � 1.2, coal permeability ratio increases from 0.2 to 0.3. For all three cases with three different bedding orientations, coal permeability always increases with decreasing gas pressure during production, indicating the dominant role desorption plays in the evolution of permeability under uniaxial strain boundary condition. As shown in Figure 12(a), at production time t � 100D, the gas pressure of the model with horizontal bedding planes reduces to 0.12 MPa; the gas pressure of the model with bedding orientation oblique to the horizontal plane declines to 0.23 MPa; and the gas pressure of the model with vertical bedding planes drops to 0.74 MPa. Identical to the reason mentioned before, this is because the vertical bedding direction of the model has the smallest permeable flow path that slows the desorption rate and in turn slows down the pressure decline. Figure 12(b) illustrates that cumulative gas production reaches its peak at about 100 days when both the equivalent fracture aperture coefficient η and the modulus correction coefficient δ are 1, while cumulative gas production reaches their peaks at about 200 days and 600 days for the other two cases, respectively. Compared with cases     under constant stress boundary condition, the cumulative CBM production reaches the peak earlier, which is also indicative of the strong effect of desorption on permeability.

Constant Volume Condition.
e constant volume condition is illustrated in Figure 13 where the displacement on all boundaries of the model is zero. In this case, the effective stress remains constant while the desorption affects the permeability in all directions; in other words, the desorption of the model will completely dominate the evolution of permeability and mechanical deformation [14,44]. Figure 14 shows that when η � 1 and δ � 1, the ratio of coal permeability increases from 1.7 to 12 with decreasing gas pressure. When η � 0.8 and δ � 1.2, coal permeability ratio increases from 0.8 to 4 during gas production. When η � 0.6 and δ � 1.2, coal permeability ratio increases from 0.2 to 0.3. Similar to the uniaxial strain boundary condition, independent of the bedding orientation relative to the horizontal plane, the permeability always increases during gas production, indicating the dominant effect of gas desorption on permeability evolution.
At production time t � 50D, as shown in Figure 15(a), gas pressure of this model with horizontal bedding planes declines to 0.23 MPa; the gas pressure of the model with oblique bedding planes reduces to 0.23 MPa; and the gas pressure of the model with vertical bedding planes drops to 0.54 MPa. e same reason mentioned earlier, smaller permeable flow path leading to slower desorption rate, applies to this case as well. Figure 15(b) illustrates that cumulative gas production reaches the maximum at about 50 days when η � 1 and δ � 1 and at about 100 days and 1000 days in the other two cases, respectively. Compared with the constant stress and uniaxial strain conditions, the CBM production reaches the production limit the earliest, indicating the strongest desorption effect on the evolution of permeability under constant volume condition.

Conclusions
A better understanding in permeability for bedded coal is of great necessity to accurate evaluation of coalbed methane production. By considering the influence of coal bedding structure on mechanical behavior and seepage characteristics, an improved permeability model for embodying anisotropy was established based on the equivalent fracture aperture coefficient. e permeability model was executed into the finite element software COMSOL and applied to three classical boundary conditions. Our simulation results show that evolutions of permeability, gas pressure, and cumulative production are capable of representing the anisotropic permeability of bedded coal. Different equivalent fracture aperture coefficients lead to the significant evolution tendency of anisotropic permeability under the same conditions. e magnitude of anisotropic permeability changes with gas pressure and time depends on both the equivalent fracture aperture coefficients. Our model shows that the bedding characteristics of permeability have a significant effect on gas production forecasts. However, the coal's nonlinear deformation behavior is not within the scope of this paper which actually exists in the process of coalbed methane extraction, aimed for simplifying the problem. We will continue to study these issues in the future.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest with respect to the research, authorship, and/or publication of this article.