Apparent permeability in tight gas reservoirs combining rarefied gas flow in a microtube

In tight gas reservoirs, the major flow channels are composed of micro/nanopores in which the rarefaction effect is prominent and the traditional Darcy law is not appropriate for gas flow. By combining the Maxwell first-order slip boundary condition and Navier–Stokes equations, a three-dimensional (3D) analysis of compressible gas slip flow in a microtube was presented, and the flux rate and pressure variation in the flow direction were discussed. Subsequently, by superimposing the Knudsen diffusion, a gas flux formula applicable to a larger Knudsen number was further proposed and satisfactorily verified by two groups of published experimental data in microtubes or microchannels in the membrane. The results indicate that slip flow and Knudsen diffusion make an important contribution to the total gas flow in the microtube, and their weight increases with an increase in the Knudsen number. By substituting the gas flux formula into Darcy’s law for compressible gas, a new apparent permeability model for tight gas reservoirs was proposed, in which the slippage effect and Knudsen diffusion were synthetically considered. The results indicate that the apparent permeability of tight reservoirs strongly depends on the reservoir pressure and pore-throat radius, and an underestimation value may be predicted by the previously published models. This study provides a case study for evaluating these apparent permeability models, which remains a challenging task in the laboratory.


Introduction
Tight sandstone gas is an important energy source. A sound understanding of reservoir properties and gas flow mechanisms is required to effectively develop these natural gas resources (Gensterblum et al., 1995;Sander et al., 2017;Rutter et al., 2022). Gas molecules generally move without rules; this movement is commonly described by the molecule mean free path. The Knudsen number (K n ), defined as the ratio of the mean free path and the characteristic length of the microchannels, is applied to determine whether the gas flow meets the continuity assumption (Freeman et al., 2011). As the characteristic length of microchannels approaches the mean free path, the rarefaction effect is significant, resulting in the breakdown of the continuity assumption (Javadpour et al., 2007;Civan, 2010). Based on K n , gas flow regimes can be divided into continuum flow (K n <10 −3 ), slip flow (10 −3 <K n <10 −1 ), transition regime (10 −1 <K n <10), and free molecular regime (K n >10) (Roy et al., 2003;Civan, 2010;Sakhaee-Pour and Bryant, 2011;Swami et al., 2012;Kim et al., 2016).
The Navier-Stokes equations and Darcy's law are suitable for continuum flow. For slip flow, two approaches have been proposed to describe the gas flow. One approach combines the Navier-Stokes equations with different slip boundary conditions (Shen et al., 2007), and the other uses a molecular-based model (Fukui and Kaneko, 1988). Compared with molecular-based models, slip boundary models are simpler and more efficient because they can be applied with a continuum description (Shen et al., 2007). Maxwell (1878) first proposed the slip model and expected two kinds of reflections existing as gas molecules to collide with the wall, a part of (σ v ) the diffuse reflection, and the remaining part of (1-σ v ) the specular reflection. σ v depends on the surface material, surface roughness, incident angles of gas molecules, gas types, temperature, and pressure (Arkilic and Schmidt, 1997). Although the derivation process of the slip model was not strictly physical, it correctly reflected the dependence of the slip velocity on the reflection and sectional velocity gradient. Subsequently, the 1.5order and second-order slip boundaries were proposed based on kinetic theory (Kennard, 1938;Mitsuya, 1993;Wu and Bogy, 2001;Hadjiconstantinou, 2003). These models usually act for K n < 0.1, and the second-order slip model is likely to overestimate the flow rate at large K n . Meanwhile, the legitimacy of 1.5-order and second-order slip models has also been debated (Shen et al., 2007). A general slip boundary condition declared to be appropriate for the entire Knudsen range was developed by Beskok and Karniadakis (1999). However, the slip coefficient needs to be determined using either experimental or DSMC data. Simultaneously, a similar slippage effect was obtained by employing the effective viscosity near the wall or applying the extended Navier-Stokes equations with no-slip boundary conditions (Beskok and Karniadakis, 1999;Chakraborty and Durst, 2007;Roy and Chakraborty, 2007;Arlemark et al., 2010;Michalis et al., 2010;Agrawal, 2011). For the transition and free molecular regimes, the collision between the gas molecules and the wall is more prominent. The Navier-Stokes equations have to be substituted by Boltzmann or Burnett equations (Agarwal et al., 2001;Freeman et al., 2011;Rahmanian et al., 2012), but it is difficult to obtain analytical solutions. Subsequently, some authors have attempted to estimate the total gas flux in microchannels by superimposing the slip flow with Knudsen diffusion (Javadpour, 2009;Zhang et al., 2015), which shows a reasonable match with experimental data (Javadpour, 2009). In addition, many experiments focusing on gas flow at different scales and shapes of microchannels have been conducted (Arkilic and Schmidt, 1997;Maurer et al., 2003;Hsieh et al., 2004;Ewart et al., 2006;Velasco et al., 2012), which have been reviewed in detail by Morini et al. (2011).
Unlike the Darcy permeability, the apparent permeability for gas transport in tight reservoirs results from the joint effects of rock properties and gas flow regimes, which makes a big difference in the production characteristics of natural gas. Pores and throats are small in tight gas reservoirs. As the gas reservoir pressure decreases, the apparent permeability becomes as much as 10 times larger than the Darcy permeability (Darabi et al., 2012;Kalarakis et al., 2012). Klinkenberg (1941) first developed a famous formula that considers the slippage effect based on the kinetic theory, but it is not appropriate for flows at high K n (Sakhaee-Pour and Bryant, 2011).
By adopting the lattice Boltzmann method, Tang et al. (2005) confirmed the Klinkenberg equation and concluded that the secondorder term of K n needed to be added for gas flow at high K n . Meanwhile, experimental results have shown that the gas slippage factor is a power function of the Darcy permeability versus porosity (Jones and Owens, 1980;Florence et al., 2007), but it might not hold for tight gas reservoirs (Swami et al., 2012). Michel et al. (2011) proposed a dynamic slippage factor considering the Knudsen diffusion coefficients, but the slippage effect was not accounted for (Swami et al., 2012). Moreover, by combining the Beskok-Karniadakis model (1999), Civan et al. (2011), andSakhaee-Pour andBryant et al. (2012) introduced apparent permeability as a function of K n . As the total gas flux in microand nanopores consists of viscous flow driven by the pressure gradient and Knudsen diffusion under the action of a gas concentration gradient, an accurate apparent permeability model may need to consider the aforementioned two parts in the meantime (Harley et al., 1995;Rushing et al., 2003;Javadpour, 2009). It is closely related to rock properties, gas types, and environmental conditions (Javadpour, 2009;Chen et al., 2015;Zhang et al., 2015).
This study theoretically analyzed the gas slip flow in a 3D microtube by combining the Navier-Stokes equations and firstorder slip boundary condition. A gas flux formula in the microtube was proposed by superimposing the Knudsen diffusion, which was suitable for a higher K n . Subsequently, a new apparent permeability model for tight gas reservoirs, including the slippage effect and Knudsen diffusion, was proposed, which contributes to a further understanding of the gas flow mechanism and production characteristics in tight reservoirs.
2 Theoretical gas flow formula in a microtube 2.1 The derivation of the theoretical gas flow formula As gas flows in a 3D microtube (Figure 1), ignoring the mass force and temperature variation and assuming that viscosity is Frontiers in Earth Science frontiersin.org 02 constant and the second viscosity coefficient μ′ equals 0, the Navier-Stokes equations for compressible gas under steady-state conditions are written as follows (Chen, 2017): where fluid pressure is denoted as p, Pa; gas viscosity μ, Pa s; and velocity V, m/s. By ignoring the changes in radial velocity (v) and angular velocity (w) over angle (θ) (Eq. 1), the continuity equation for steady-state and compressible gas flow was expanded in cylindrical coordinates: For ideal gases, p = ρR a T/M, where M represents the molar mass (kg/kmol); R a refers to the gas constant, 0.008314 (Pam 3 )/(kmolK); ρ is the density, kg/m 3 ; and T is the temperature (K). The nondimensionalization of these variables is as follows: streamwise velocity (u) and radial velocity (v) are normalized by the areaaveraged velocity at the outlet, the streamwise coordinate (z) by the length of the microtube (L), radial coordinate (r) by the radius of the microtube (R), and p and ρ are the pressure and density at the outlet, respectively. Meanwhile, ε = R/L <<1 and δp = p i -p o . With these assumptions, Eqs. 2-4 were expressed in a non-dimensional form as follows: It is clear that the radial velocity equals 0 at the wall of the microtube, where the streamwise velocity is expressed by the Maxwell first-order slip velocity (Eq. 8). Meanwhile, K n was replaced by the ratio of K n at the outlet (denoted by K o ) and the non-dimensional pressure to obtain the non-dimensional boundary condition:ũ After expandingp,ũ, andṽ in powers of ε Eqs. 9 10, we found that the 0th-order radial velocity was not included according to (Eq. 7):pr ur,z ( ) ũ 0r ,z ( )+ εũ 1r ,z ( )+ ε 2 . . .
Subsequently, based on Eqs 9, 11 and solving Eq. 6 with the boundary conditions (Eq. 8), the non-dimensional 0th-order streamwise velocity is written as (Eq. 12) By substituting Eq. 12 into Eq. 7 and integrating it once in the non-dimensional microtube radius, the expression of the nondimensional first-order radial velocity (Eq. 13) can be obtained: Combining Eq. 8, Eq. 13 was simplified to It was expanded in a second-order polynomial of the non-dimensional pressure (Eq. 14), which contributes to solving the non-dimensional pressure (Eq. 15). Note that the non-dimensional pressure is positive and equal to 1 at the outlet of the microtube: The outlet flow rate is Then, by substituting Eq. 14 into it, the gas rate at the outlet is written as Eq. 16, and the outlet mass flux is derived as shown in Eq. 17: where J o denotes the mass flux at the outlet, kg/s/m 2 ; and q refers to the flow rate at the outlet, m 3 /s. In addition, the mass flux brought about by the Knudsen diffusion can be described by the Knudsen diffusion model (Eq. 18) (Javadpour, 2009). The total mass flux at the outlet of the microtube, including viscous flow, slip flow, and the Knudsen diffusion, is presented by a combination of Eq. 17 and Eq. 18, which is theoretically valid for slip and transition flow regimes (P = p i /p o ) (Javadpour, 2009;Swami et al., 2012;Zhang et al., 2015). R b is the gas constant [8.314 J/(mol K)]: Ewart et al. (2006) conducted an experimental investigation of the gas mass rate for an isothermal steady flow in cylindrical  Figure 2 (black dots). The gas mass flux model proposed in this study (Eq. 19) was used to model the experiment (σ v = 0.8) (Javadpour, 2009), which showed good agreement with the experimental results, and the relative error was less than 5%. In addition, Roy et al. (2003) used a membrane with a pore size of 200 nm and thickness of 60 μm to study the gas micro-flow. Figure 3 shows the mass flux versus pressure drop during this experiment, with an exit K n of 7.36. The results calculated by Eq. 19 were simultaneously depicted in Figure 3 with an exit K n of 14.66. K n was expressed as the ratio of the mean free path to the microtube diameter in Roy's study. Based on Figure 3, the model results based on Eq. 19 approach the experimental results well at small pressure drops but deviate from the experimental results as the pressure drop increases. This is mainly because of the irregular shape of the microchannels in the membrane, which may lead to a smaller transport capacity than that of the regular microtube used in the model.

Model results
According to the models in Section 2.1, the gas flow characteristics, including the pressure distribution, gas mass flux, and mass flow ratio for different flow regimes in the microtube, will be further discussed. Nitrogen was selected as the gas source, the temperature was set to 300 K, and σ v = 0.8 in the models.
Pressure distribution in the microtube for slip flow. Based on Eq. 15, the pressure distribution in the microtube at different inlet pressures is shown in Figure 4. The inlet pressures were 1, 2, 3, and 4 MPa (the outlet pressure was 0.7 MPa), and the average K n values were 0.083, 0.052, 0.038, and 0.03, respectively. At the same time, the length and diameter (d) of the microtube were 1 cm and 200 nm, respectively. When gas flows into the microtube, the flow rate increases, whereas the pressure decreases, which results in a change in the gas density in the flow direction. As shown in Figure 4, there was a non-linear variation in the pressure along the flow direction. The non-linearity of the pressure distribution in this microtube should be caused by gas compression, whereas the rarefaction effect should be the opposite. With a larger ratio of the inlet and outlet pressures, the curvature of the pressure curve is   Roy et al. (2003).

FIGURE 4
Pressure variation versus microtube coordinate at different inlet pressure.

Frontiers in Earth Science
frontiersin.org smaller, which can also be speculated from Eq. 15. In addition, the curvature was slightly smaller for the larger microtube at the same pressure ratio (Figure 5), and the influence of the gas type was almost negligible. Gas mass flux versus K n for different microtubes. Figure 6 shows the gas mass flux in microtubes of different scales with variational K n . The microchannel length was 1 cm, and the diameter of the microtube ranged from 50 to 1,500 nm. Simultaneously, the outlet pressure was set to 0.002 MPa. Figure 6 shows that the gas mass flux in the microtubes gradually decreases with an increase in K n . The gas mass flux for different microtubes is in accordance with K n smaller than 1, and the absolute slope of the coincident curves continues to decrease with an increase in K n . However, for K n >1, a slower decrease occurred in the smaller microtubes, resulting in the mass flux curves beginning to diverge between the different microtubes. Additionally, the slippage effect and Knudsen diffusion were more significant in smaller microtubes. However, the average pressure was higher for smaller microtubes to achieve the same K n .
Mass flow ratio versus K n . During gas flow in a microtube, the flow regimes include viscous flow, slip flow, and Knudsen diffusion, which are mixed in different proportions . According to Eqs 17-19, the gas mass flux ratios in different regimes with variations in K n are shown in Figure 7. The length and diameter of the microtubes were 1 cm and 200 nm, respectively. As shown in Figure 7, the viscous flow is dominant, and the Knudsen diffusion is almost neglected when K n is smaller than 0.1. As K n increases, the slip flow and Knudsen diffusion play a more important role in the gas flow. At the later stage of the transition flow, the Knudsen diffusion accounts for almost 50% of the total gas mass flux, which approaches the sum of viscous and slip flows. The ratio of Knudsen diffusion is constant in the free molecular flow regime and shows a slight increase with an increase in the microtube diameter. In brief, slip flow and Knudsen diffusion must be considered for the rarefied gas flow in the microchannels.

Apparent permeability
At present, gas production in tight gas reservoirs presents a tendency for rapid growth. To a great extent, it is determined by the apparent permeability of gas transport in porous media. The gas flow mechanism in the microtubes can be used to promote the understanding of gas flow in tight reservoirs. Based on Darcy's law for compressible gas, the apparent permeability for tight porous media (Eq. 20) can be calculated (Javadpour, 2009;Zhang et al., 2015): Pressure variation versus microtube coordinate at different microtube diameters.

FIGURE 6
Gas mass flux versus K n at different diameters of microtubes.

FIGURE 7
Dynamic mass flux ratio by different flow parts with K n variations.
where k g denotes the permeability, m 2 ; and q m refers to the gas mass rate (kg/s). By combining the new gas flux formula in the microtube presented in Section 2.1 (Eq. 19), a new apparent permeability model (Eq. 21) was derived, which includes viscous flow, slip flow, and Knudsen diffusion. Meanwhile, the ratio of apparent permeability and Darcy's permeability in porous media was written as Eq. 22. This indicates that the difference between the apparent permeability and the Darcy permeability depends primarily on the pore radius and pressure. The apparent permeability of tight reservoirs significantly depends on pressure and pore radius (Roy et al., 2003;Freeman et al., 2011): The correlation of the k app /k D ratio with pressure and pore diameter is shown in Figure 8. The outlet pressure was 0.002 MPa, and the pore diameter ranged from 1 to 6,000 nm. As shown in Figure 8, the ratio is largest at the smallest inlet pressure of 0.01 MPa. In addition, the ratio decreased with increasing pressure, which indicates that in tight reservoirs, the viscous flow was dominant at high pressure, whereas the production gas was mainly driven by slip and Knudsen diffusion at low pressure.
A comparison of the new apparent permeability model and the previously published models is shown in Figure 9. The k app / k D ratio based on Klinkenberg's (1941) model is the smallest, especially for K n >0.1 because only the slippage effect was considered. In the model of Chen et al. (2015), the Knudsen diffusion was considered, but slippage was excluded. The k app /k D ratio of Chen et al.'s model is in good agreement with that of Civan's model (2010), which introduced Beskok and Kaniadakis's (1999) correction factor and is slightly larger than that in Klinkenberg's model. Zhang et al. (2015) combined the Klinkenberg model (1941) with the Knudsen diffusion. Javadpour (2009) combined slip flow using the theoretical dimensionless coefficient (F) introduced by Brown et al. (1946) with the Knudsen diffusion. The ratio predicted by the models of Zhang et al. and Javadpour is approximate while being larger than that of the models proposed by Klinkenberg, Chen et al., and Civan. The results of the new model most approximate those of Javadpour's model and are slightly larger than those of Chen et al.'s model because of the different slippage factors considered. It can be concluded that these apparent permeability models are all more or less equivalent in the slip flow regime and differ in the transition and free molecular regimes because the different gas flow mechanisms and theoretical models were accounted for.

Discussion
Pores and throat spaces in tight gas reservoirs are complex pore network systems. According to Zou et al. (2012), the pore diameter of tight gas reservoirs is in the range of 5-700 nm. When the pore diameter is less than 10 nm, surface diffusion occurs on the surface of the pores (Yang et al., 2016). Therefore, surface diffusion was not considered in this study because the vast majority of pores in reservoirs are larger than 10 nm. The range of K n in tight gas reservoirs as the reservoir pressure decreases from 35 MPa to 1 MPa, as shown in Figure 10. Figure 10 suggests that K n decreases with an increase in reservoir pressure and pore diameter, whereas the main flow Effect of pressure and microchannel diameter on the k app /k D ratio.

FIGURE 9
Comparison of the new apparent permeability model and the previously published models.
Frontiers in Earth Science frontiersin.org 06 regimes are the slip flow and transition regime (K n <10). Darcy's law is not appropriate, and the new apparent permeability model is theoretically valid for tight gas reservoirs. However, it may be impractical to test the validity of these models in a laboratory study of reservoir rocks. As gas flows in true reservoir rocks, the variation in apparent permeability comes from two parts: the changes in gas flow regimes and the stress sensitivity of rocks (i.e., the pore space changes in rocks). For a high K n , it is very difficult to weigh the contribution of the rock stress sensitivity because of the complicated variation in the gas flow regimes. Some authors excluded the slippage effect in tight sandstones by conducting experiments at high pressure (Li et al., 2009) or directly neglected the effect of effective stress on apparent permeability (Yuan et al., 2016). Therefore, it may be better to directly include the effective stress in the apparent permeability model together with the transport mechanism. Some researchers have attempted to do some work and have obtained apparent permeability models based on the linear effective stress law (Cao et al., 2016;Xiao et al., 2019). However, effective non-linear stress is common in tight rocks (Li et al., 2014;Xiao et al., 2019), and the apparent model needs to be improved in future work.
Moreover, based on network simulations, Bernabé et al. (2010) found that the matrix permeability in 3D simple cubic, FCC, and BCC networks obeyed the "universal" power laws, k∝(z-z c )β, where β is a function of the standard deviation of the pore radius distribution and z c is the percolation threshold in terms of the coordination number (z). Permeability can be further expressed as a combination of scale-invariant parameters (Bernabé et al., 2010): where r H is the mean pore radius, l is the mean pore separation distance, and w is a function of the standard deviation of the pore radius distribution. Our results showed that rock permeability is related to the pore structure. Thus, it may be beneficial to considering the pore structure parameters in different apparent permeability models.

Conclusion
A 3D analysis of the Navier-Stokes equations for compressible gas flow in a microtube was presented by combining the first-order slip boundary condition. The non-linear pressure variation was analyzed for the gas slip flow. Subsequently, a new gas flux formula in the microtube, including viscous flow, slip flow, and Knudsen diffusion, was proposed, which is in good agreement with the published experimental data.
According to the calculated Knudsen number range, it can be concluded that slip and transition flow are the main flow regimes in tight gas reservoirs. The gas mass flux driven by slip and Knudsen diffusion significantly contributes to the total gas mass flux in the micropores and nanopores, and their weight increases with an increase in the Knudsen number.
A new apparent permeability model for tight gas reservoirs, considering the slippage effect and Knudsen diffusion, is presented based on the new gas flux formula and Darcy's law. By comparing the new model with previous models, the results of the new model were close to those of Javadpour's model, and other previous models may underestimate the apparent permeability of tight reservoirs. The results show that the apparent permeability strongly depends on the reservoir pressure and pore-throat radius.

FIGURE 10
Range of K n in tight gas reservoirs as p = 1-35 MPa.
Frontiers in Earth Science frontiersin.org