Apparent Permeability Prediction of Coal Matrix with Generalized Lattice Boltzmann Model considering Non- Darcy Effect

College of Civil Engineering, Hunan City University, Yiyang 413000, China Key Laboratory of Key Technologies of Digital Urban-Rural Spatial Planning of Hunan Province, Hunan City University, Yiyang 413000, China School of Mechanics & Civil Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China


Introduction
Permeability is a key factor in the exploitation of coalbed methane, which directly determines whether the target coal reservoir is valuable for mining. For a long time, people evaluated the permeability of the coal reservoir based on Darcy's law. However, recent extensive studies [1][2][3][4] show that the non-Darcy effect at the pore scale has a great contribution to the gas flow flux. If the influence is ignored, the permeability obtained will have a great deviation.
Coal is a typical tight porous medium, with a wide distribution of pore size, which is distributed in the range of 1-100 nm [5,6]. Therefore, gas flow in coal in the form of a multiscale flow process includes viscous flow, slip flow, transition flow, and the Knudsen diffusion (free molecular flow). Under the exploitation of coalbed methane, the dominant flow mechanism also changes with the reservoir pressure [7]. Specifically, the adsorbed gas in the nanopore of the matrix is desorbed as free gas when the pressure drops in the process of drilling and drainage, and the resulting pressure difference then leads the free gas to diffuse into the fracture system. After the gas pressure tends to be stable, gas flows in the form of viscous flow in the wellbore, macropore, or fracture network under a pressure gradient. However, due to the small pore throat of the coal matrix, when the average free path of gas molecules is close to the pore size, the gas molecules have a nonzero velocity on the wall, which will result in the occurrence of gas slippage [8]. In addition to the discontinuous effect causing the apparent permeability of the gas to be higher than the intrinsic permeability, the concentration difference of the adsorbed gas on the pore surface will also make the gas diffuse to the low concentration area; this phenomenon is called surface diffusion [9,10]. Previous studies [11,12] have confirmed that the total amount of adsorbed gas is more than 85% in coal reservoirs; therefore, the surface diffusion of adsorbed gas is also an important transfer mechanism [9,[13][14][15].
The coal seam is a kind of low permeability reservoir. According to well test data [16], the permeability of most CBM reservoirs is less than 1 mD, and permeability is one of the main control factors limiting the productivity of CBM. Villazon et al. [17] established a permeability model that took into consideration the non-Darcy effect and employed a lognormal distribution as pore size distribution. However, the model only considered the effect of gas slippage, and other involved flow mechanisms were not taken into account. Based on this model, Xiong et al. [18] proposed a fully coupled, free and adsorptive phase transport permeability model taking into consideration the impact of the effective pore radius and the non-Darcy gas flow. This model linearly superimposed the flow flux contributed by gas slippage and surface diffusion, and took into account the effect of adsorbed gas molecules on the pore space, but did not regard the influence of the Knudsen diffusion. Wu et al. [19] proposed the weighting factor for automatically adjusting the contribution of viscous flow and the Knudsen diffusion to the total flow flux by the Kn within different pore size ranges and reservoir pressures.
The maceral composition of coal consists of organic matter (OM) and inorganic matter (IOM) [11,12]. The basic organic components of coal that can be identified under an optical microscope are called organic macerals, which are macerals converted from plant residues. Inorganic macerals are the minerals in coal observed under a microscope. The IOM in coal is usually finely dispersed in OM. Zhao et al. [6] tested six groups of coal samples based on the American Society for Testing and Materials (ASTM 2009-D5142) and found that the contents of mineral components were 2.6%, 3.0%, 5.7%, 7.7%, 13.9%, and 17.5%, respectively. Therefore, the process of methane migration in the coal matrix involved complicated multiscale gas flows in porous media, which can be assumed to divide into three parts: (1) gas flow with non-Darcy effect and gas adsorption in OM, (2) gas flow with non-Darcy effect in IOM, and (3) gas free flow in microfractures.
In recent years, the Lattice Boltzmann Method (LBM) on fluid flow in porous media attracted wide attention [20]. LBM was naturally suitable for pore-scale modeling, since it can easily deal with fluid flow in a channel [21,22]. However, understanding the fluid dynamics in porous media and predicting effective transport properties (permeability, effective diffusivity, etc.) is of paramount importance for practical applications. It is essential to establish a model capable of dealing with the fluid flow at a larger scale. Guo and Zhao [20] developed a generalized Lattice Boltzmann model (GLBM) for fluid flow through porous media, then the properties such as porosity and permeability were introduced into each site of the porous media. GLBM is suitable for REV-scale modeling, which is capable of handling fluid flow in porous media with permeable and impermeable matter and channel. Chen et al. [23] first established a model considering Klinkenberg's effect and introduced the apparent permeability as a parameter into the GLBM. Wang et al. [10] proposed a DGM-GMS (i.e., the dusty gas model and the generalized Maxwell-Stefan model) coupled permeability model based on the GLBM, which takes into consideration the permeability of OM and IOM at the REV scale, respectively.
In this paper, based on previous studies, we proposed a GLBM model for predicting coalbed methane permeability. This model included the following sections: (1) The generalized Lattice Boltzmann method was employed to solve the generalized Navier-Stokes equations, which was applied to the calculation and prediction of local permeability at representative elementary volume (REV) scale. (2) The fully coupled viscous flow, slip flow, transition flow, and the Knudsen diffusion were considered in the simulation, and the contribution of these flow mechanisms to the total flow flux was automatically adjusted by the weight factor. (3) The impact of adsorbed gas molecules on the pore space and the contribution of surface diffusion to the total flow flux were taken into consideration. (4) The separate evolution of the local permeability of OM and IOM in the coal matrix was taken into consideration.

Generalized Model for Fluid Flow in Porous Media
For isothermal flows of incompressible fluids in porous media at the REV scale, the generalized Navier-Stokes equations which were proposed by Nithiarasu [24] are capable of the simulation.
2.1. Generalized Navier-Stokes Equations. The governing equations of mass and momentum for the generalized Navier-Stokes equations can be given by where ε is the porosity, ν e = νJ is the effective viscosity, J e is the viscosity ratio, and F represents the total body force including both medium resistance and external forces and can be given by where ν is the kinematic viscosity of the fluid and G is the external body force. F ε is the geometric function and K is the permeability of porous media; both these two parameters are related to the porosity. For a porous medium 2 Geofluids composed of solid particles with diameter d p , the Ergun correlation [25] gives Note that equation (4) was the original permeability of the GLBM; as for three-component porous media, there are three sets of permeability that should be decided: OM, IOM, and fractures. The permeability of fractures can be calculated directly by equation (4) with a porosity of 1, then the GLBM recovered to the LBM and the region of the fracture naturally becomes a channel. The apparent permeability of OM and IOM will be proposed in Section 3, which can be substituted into equation (3).

LB Model for Generalized Navier-Stokes Equations.
Guo and Zhao [20] constructed a LB model which can be used to solve the generalized Navier-Stokes equations, and the corresponding evolution equation of the particle distribution function is where where f i is the discrete density distribution function, f eq i is the local equilibrium equation, e i is the discrete velocity of particles, δ t is the time step, τ is the relaxation time, and ω i is the weighting coefficient.
As in the standard LB model, the density and velocity of flow are defined as Because the force F also contains the flow velocity u, the velocity can be explicitly given by where v is a temporary velocity defined by and the parameters c 0 and c 1 are given by By using the Chapman-Enskog technique with the pressure p = c 2 s ρ/ε, and the effective viscosity ν e = c 2 s ðτ − 0:5Þδ t , the generalized LB model can recover to equation (1) and equation (2) in the incompressible limit.

The Apparent Permeability Model of CBM Taking Multiple Flow Mechanisms into Consideration
The pores of OM in coal are usually less than 100 nm in diameter and have a large amount of gas adsorbed on their surfaces [26,27]. In nanopores, the probability of gas molecules colliding with the pore wall is significantly increased, which provides a nonnegligible momentum source for gas transfer near the pore wall. The flux contribution of gas slip to the total flow increases, which is a significant deviation from Darcy's law. For the nanopores in OM, the presence of the adsorbed gas molecules adds another complication. Firstly, the gas molecules adsorbed on the pore surface will transfer along the wall due to the concentration difference. Secondly, the gas adsorption amount increases with the pressure; thus, the adsorption efficiency of the pore space constituted by the reshaped channel will be related to the pressure. According to Poiseuille's law, for capillary tubes with a radius r, the intrinsic permeability (i.e., the absolute permeability, which is the viscous flow permeability of an unreactive ideal fluid) can be expressed as The volume flux of gas viscous flow is where μ is the dynamic viscosity coefficient of the gas (Pa ⋅ s), p is the pressure (Pa), and Q ∞ is the volume flow flux (m 3 /s).
As gas molecules adsorbed on the inner surface of a capillary, the loss of the cross-sectional area for free gas transmission may be large [18]. Let us assume that the gas adsorption conforms to the Langmuir monolayer adsorption theory. When the gas adsorption amount of specific sites in the capillary tube reaches saturation, the section of the capillary tube should subtract the molecular diameter of 3 Geofluids the adsorbed gas molecules. Therefore, the effective radius of the capillary tube allowing free gas to pass through can be written as The gas adsorption amount is a function of pressure. When the reservoir pressure is p, the coverage of adsorbed molecules on the inner surface of the capillary tube should be considered. Therefore, based on the Langmuir singlelayer adsorption theory [28], the effective radius can be modified to Then, the relationship between effective porosity and initial porosity is After considering the reduction effect of gas adsorption to the pore space, the volume flow flux can be rewritten as Consequently, the permeability involving gas adsorption becomes a function of pressure: The Navier-Stokes equations with appropriate slip boundary conditions are sufficient for modeling gas flow in the slip flow regime [29]. However, since the mean free path of the gas molecules is comparable to the pore size in the transitional flow regime, momentum transfer between gas molecules and between the gas molecules and the wall is significant. Therefore, the combined effect of viscous flow and the Knudsen diffusion should be considered simultaneously in the transitional flow regime [14,15].
Klinkenberg [30] proposed the apparent permeability equation, which gave the linear relationship between the apparent permeability and the reciprocal of the pore pressure, as follows: where f ðKnÞ is the correction coefficient and Klinkenberg's correction is the first-order correction; Beskok and Karniadakis [8] give a second-order correction, which is suitable for describing gas flow in four flow regimes (includes viscous flow, slip flow, transition flow, and free molecular flow): where b is the gas slippage coefficient, generally b = −1; α is a gas rarefied coefficient. Civan [31] gives The Knudsen number is defined as the ratio of the mean free path of molecules to the characteristic pore size of porous media: The mean free path is defined as Considering the impact of the effective pore size and gas slippage on apparent permeability, the effective radius r eff should be employed instead of the pore radius in equation (19), and the corresponding permeability should be corrected to The mass flow flux with the gas slippage taken into consideration is As the Knudsen number is in the range of 10 −3 < Kn < 10 −1 , the collision between gas molecules is dominant in the process of gas transport, and the probability of collision between gas molecules and the wall is relatively small, but it cannot be ignored. When the velocity of gas molecules near the wall is nonzero, gas slippage will occur, which can be described by the ideal gas slip flow equation [14]. If Kn > 10, the collision between gas molecules and the nanopore wall is dominant, and the behavior of gas transfer can be described as the Knudsen diffusion. According to Fick's law, the flux contributed by the Knudsen diffusion can be written as The corresponding volume flow flux is where D Kn is the Knudsen diffusion coefficient for gas flow in a single capillary, which can be defined as [2].

Geofluids
Then, the mass flow flux contributed by the Knudsen diffusion is Due to the different concentrations of adsorbed gas in each site on the inner surface of a nanopore, surface diffusion makes a significate role in gas transport. The adsorbed gas and the bulk flow gas can be connected by the relationship of the Langmuir isotherm adsorption equation: The mass flow flux contributed by surface diffusion is In the viscous flow regime, the probability of intermolecular collisions dominates the total collisions, which occur under high-pressure conditions or in systems without boundaries. For the Knudsen diffusion, only the collision between the gas molecules and the solid walls is considered, which occurs in systems with near-vacuum pressure or extreme restriction [10]. Due to the restriction of gas molecules in coal pores, both viscous flow and the Knudsen diffusion cannot be neglected under high pressure. Therefore, the combined effect of viscous flow and the Knudsen diffusion should be considered in gas transport. Besides, both viscous flow and the Knudsen diffusion are considered in the same computational region; it is necessary to determine the reasonable weighting coefficients of each and their respective contributions to the overall phase transport [19]. The combination of gas slippage and the Knudsen diffusion has been established by Darabi et al. [2,32] using different weighting coefficients and linear summing. Based on the analysis of collision frequencies between molecules and between molecules and pore walls, the weighting factors that distinguish the contribution of these two flow mechanisms to the total flow can be obtained [9,15] ω ν = 1 1 + Kn , Considering the influence of porosity and tortuosity of porous media, the total flux can be expressed as The total mass flow of gas transport in the coal matrix can be considered as a combination of viscous flow (includes the slip flow and transition flow), the Knudsen diffusion, and surface diffusion. In this paper, based on the combined effects  5 Geofluids of viscous flow, the Knudsen diffusion, and surface diffusion, the apparent permeability of OM is obtained: Nanopores are also distributed in the IOM of coal, which may have a certain adsorption effect on gas molecules, but they can be ignored compared to the strong adsorption effect of OM. Therefore, it is assumed that the pores of IOM only serve as channels for fluid percolation. According to equation (34), we can get the apparent permeability of IOM:  Table 1, and the simulation results are shown in Figure 1. It can be found that the simulation results are in good agreement with the experimental data. It is noteworthy that the pore size was approximately 200 nm; thus, the effect of surface diffusion was weak and ignored in the simulation. Therefore, the model proposed in this paper can accurately describe the coalbed methane in the nanopores in the real state.

The Impact of Component Content and Fracture on
Permeability. Figure 2(a) shows the microstructure of the coal reproduced from the previous study [35]. It can be observed that coal is a three-component porous medium composed of OM, IOM, and fracture. Among them, OM accounts for the majority of the total volume and is randomly divided into irregular shapes by the fracture, and IOM is finely dispersed in the OM. There are many methods for porous medium random reconstruction. Chen et al. [23] reconstructed shale matrix reconstruction images based on the Quartet Structure Generation Set (QSGS). However, from Figure 2(a), we can find that the dispersed IOM is roughly circular, and there are overlaps in many positions, which makes the distribution form more complex and increases the difficulty of reconstruction. According to the morphology and distribution characteristics of the three components, we reconstructed a three-component porous medium, as shown in Figure 2(b), which roughly represents the microstructure of the coal composed of microfracture (red), OM (blue), and IOM (green). Besides, due to the limitation of observation technology, both OM and IOM contain a large number of pores that cannot be shown in the figure. These pores play an important role in methane migration, which should be considered in the model.
Zhao et al. [6] tested the macerals of six groups of coal samples following ASTM 2009-d5142 and found that the contents of IOM were 2.6%, 3.0%, 5.7%, 7.7%, 13.9%, and 17.5%, respectively. The results show that the content of OM and IOM in coal is markedly different, and IOM has a wide range of contents.
To explore the influence of the volume fraction of OM and IOM on permeability, we reconstructed 8 groups of porous media with the random distribution of OM and  Figure 3. To ensure the generality of the research, the content of IOM in the reconstructed porous media covers from 7% to 52%. Other parameters used in the simulation are shown in Table 2. A global permeability is employed to evaluate the impact of the content of each component on the permeability. The global permeability can be obtained according to Darcy's law, where u is the mean fluid velocity (m/s) and ∇p is the pressure gradient (Pa/m). Figure 4, with the variation of the volume fraction of OM, global permeability K global shows a nearlinear change despite the random distribution of OM and IOM. The overall trend follows that K global decreases with an increase in OM content. This indicates that IOM may have a great contribution to gas transport due to its wide distribution of mesopores and macropores. It is worth mentioning that the results based on IOM are finely dispersed in the OM assumption. The porosity/permeability of IOM itself is not certainly greater than that of OM, but in the state of dispersion, the interface between inorganic particles or between inorganic particles and organic particles will increase the probability of macropores. As the coal matrix is characterized by a strong heterogeneity, it may also contain a large volume   8 Geofluids of inorganic crystals, and the porosity and permeability will be considerably lower than that of OM. In fact, in the process of coalbed methane mining, we are more concerned about the data of the whole production cycle. However, due to the nature of LBM and the limited range of pressure (density) fluctuation in the simulation, the results obtained directly from the real parameters of the coal reservoir may not be correct. In this part, the density and viscosity coefficient of coalbed methane under different pressures are used to calculate the global permeability, and the results are summarized to indirectly obtain the data of the whole pressure drop process. It should be noted that the density and viscosity coefficients of methane are discrete under different conditions. We employed an online calculation software (http://www.peacesoftware.de/einigewerte/ methan.html) to calculate the corresponding parameters 9 Geofluids according to different pressures (the temperature is maintained at 303.15 K). Through observation, when the pressure is below 20 MPa, the relationship between the pressure and density of methane is close to linear, as shown in Figure 5(a). Similarly, the relationship between pressure and dynamic viscosity can be fitted by a quadratic polynomial, as shown in Figure 5(b). Figure 5(c) shows the fitting relationship between kinematic viscosity and pressure, which is a vital parameter in LBM. The pressure-matched parameters are used in the following simulation. Reservoir temperature may have a great influence on the physical properties of the gas in the process of exploitation; however, the temperatures do have some connections with other sides, e.g., gas adsorption. Therefore, we only consider the isothermal conditions at the present research.

As shown in
There is a kind of double-pore structure composed of pores and fractures, in which fractures are common in coal.
As can be seen from Figure 2(a), the morphology of fractures is "tree-like," with many branches, which connect all parts of the coal.
To further explore the impact of fractures on permeability, four groups of three-component porous media are reconstructed in this part, as shown in Figure 6. The distribution and content of OM and IOM remain unchanged in the following research. In the four groups of porous media, the fracture distributions are as follows: (a) no fracture, (b) fractal dimension D f = 2:36 (the volume fraction of the fracture is 3.5%), (c) D f = 2:43 (the volume fraction of the fracture is 5.9%), and (d) D f = 2:50 (the volume fraction of the fracture is 8.4%). According to the relationship between density and the viscosity coefficient shown in Figure 6, the simulation is carried out under different pressures. The calculation input parameters are shown in Table 2, and the simulation results are shown in Figures 7 and 8. 10 Geofluids Figure 7(a) shows that the fluid velocity is very slow when no fracture exists in the matrix, and the fluid velocity in IOM is significantly higher than that in OM. With the existence of a fracture in the matrix, as shown in Figures 7(b)-7(d), the velocity in the fracture is several orders of magnitude higher than that in the matrix, which indicates that the contribution of fractures to the gas transport capacity is large; besides, the fluid velocities in IOM and OM are also significantly increased, which dramatically enhances the gas flow. Moreover, the trend in Figures 7(b)-7(d) suggests that with the increase of the content of the fracture (or the fractal dimension), the magnitude of the fluid velocity grows up, effectively improving the gas flow. Figure 8 summarizes global permeability under different pressures and fracture development. It can be found that with the increase of the content of the fracture, the trend of global permeability variation is similar, which increases with the decrease of pressure. When no fracture exists completely in the matrix, the global permeability is in the order of 10 -4 mD in the pressure drop range, which is extremely low permeability; the results were verified by the Finite Element Method as can be seen in Appendix B. In practice, the reservoir with such a small permeability has no value for exploita-tion. However, when the matrix contains a fracture with a width of about 10 μm, the permeability is greatly stimulated, and the permeability reaches 0.4 mD as the reservoir pressure is 5 MPa. The global permeability increases with the enlargement of the content of the fracture, and when the volume fraction is at 8.4%, the permeability exceeds 0.85 mD. Figure 8(e) collects the results of Figures 8(a)-8(d) and shows that the existence of fractures is crucial for improving the permeability of a coal reservoir, and it also indicates the necessity of reservoir stimulation techniques such as hydraulic fracturing. According to the statistical results of the permeability tests of previous studies [5,6,12,16], the permeability of coal reservoirs in China is basically in the range of 1 × 10 −3~1 mD, which shows that the permeability prediction results are reliable.

Conclusion
In this paper, we established a GLBM model for predicting coalbed methane permeability in the REV scale, which was fully coupled with the viscous flow, slip flow, transition flow, and the Knudsen diffusion. The model also considered the impact of the effective pore radius and surface adsorption  With fracture D f = 2.50 0 5 1 0 Pressure (MPa) 15 20 With fracture D f = 2.43 With fracture D f = 2.36 Without fracture (e) Figure 8: The impact of fracture development on permeability. and adopted the weighting factors to automatically adjust the contribution of viscous flow and the Knudsen diffusion to the total gas flow. Besides, based on the morphology and distribution characteristics of OM, IOM, and fractures in the coal matrix, a random reconstruction algorithm of three-component porous media is developed. For the accuracy of the simulation, we adopted a series of pressurematched parameters.
The results are as follows: (1) Due to the wide range of pore size distribution, the distribution of mesopores and macropores is relatively large, which makes the contribution of IOM to the gas transport capacity also relatively large. (2) The permeability of the coal matrix itself is extremely low, and the contribution of a fracture to gas permeability is very large. (3) The greater the content of the fracture (or the fractal dimension), the larger the permeability; it is necessary to employ reservoir stimulation techniques such as hydraulic fracturing in the development process of coalbed methane.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.