Analysis of gas transport behavior in organic and inorganic nanopores based on a unified apparent gas permeability model

Different from the conventional gas reservoirs, gas transport in nanoporous shales is complicated due to multiple transport mechanisms and reservoir characteristics. In this work, we presented a unified apparent gas permeability model for real gas transport in organic and inorganic nanopores, considering real gas effect, organic matter (OM) porosity, Knudsen diffusion, surface diffusion, and stress dependence. Meanwhile, the effects of monolayer and multilayer adsorption on gas transport are included. Then, we validated the model by experimental results. The influences of pore radius, pore pressure, OM porosity, temperature, and stress dependence on gas transport behavior and their contributions to the total apparent gas permeability (AGP) were analyzed. The results show that the adsorption effect causes Kn(OM) > Kn(IM) when the pore pressure is larger than 1 MPa and the pore radius is less than 100 nm. The ratio of the AGP over the intrinsic permeability decreases with an increase in pore radius or pore pressure. For nanopores with a radius of less than 10 nm, the effects of the OM porosity, surface diffusion coefficient, and temperature on gas transport cannot be negligible. Moreover, the surface diffusion almost dominates in nanopores with a radius less than 2 nm under high OM porosity conditions. For the small-radius and low-pressure conditions, gas transport is governed by the Knudsen diffusion in nanopores. This study focuses on revealing gas transport behavior in nanoporous shales.


Introduction
In North America and China, shale gas with rich reserves and great potential has been developed efficiently (Su et al. 2015). The multiple gas transport mechanisms and reservoir features of the shales are different from those of other unconventional reservoirs. Therefore, it is essential to figure out the complicated gas transport behavior in nanoporous shales on the basis of physical experiments (Wang et al. 2015(Wang et al. , 2016a, numerical methods (Botan et al. 2013;Chen et al. 2015;Sun et al. 2017;Wu et al. 2015a, b, c;Yao et al. 2013), and theoretical apparent gas permeability (AGP) models. Experimental results (Ambrose et al. 2012;Zou et al. 2012) have shown that the pores in shale reservoirs range in nanoscale size (generally smaller than 5 nm), which causes the approximation between the molecular mean free path and shale nanopores size, and continuity assumption invalid (Song et al. 2016). Organic matter (OM) and inorganic matter (IM) have been observed in the typical shale samples, and the schematic diagram can be shown in Fig. 1.
The Knudsen number (Kn) is widely used to characterize the gas transport behavior in the microscale pores (Civan 2010;Civan et al. 2011Civan et al. , 2013. Different flow regimes can Qi Zhang and Wen-Dong Wang contributed equally to this work. Edited by Yan-Hua Sun be divided based on the Knudsen number Kn (the ratio of a molecular mean free path to the average pore diameter) (Wu et al. 2015c), including continuum flow ( Kn < 0.001 ), slip flow ( 0.001 < Kn < 0.1 ), transition flow ( 0.1 < Kn < 10 ), and molecular scale flow ( Kn > 10 ). Under typical shale gas reservoir conditions, Kn ranges from 0.0001 to 10 (Wu et al. 2014), which means continuum flow, slip flow, and transition flow coexist in shale nanopores. Moreover, Li et al. (2013) reported that the maximal proportion of adsorbed gas can reach around eighty percent of the total gas-in-place in shale gas reservoirs . Adsorbed gas decreases the organic nanopore size, and the surface diffusion of adsorbed gas molecules will occur (Akkutulu and Fathi 2012). Therefore, scholars have proposed various AGP models to characterize gas transport mechanisms. Civan et al. developed the Beskok and Karniadakis (B-K) model (Beskok and Karniadakis 1999) into the characterization of gas transport in the porous media and presented the AGP models based on the Knudsen number, considering the viscous flow and the rarefaction effect. Xiong et al. (2012) and Sigal (2013) improved Civan's models by taking the adsorbed gas into account. Based on the pore structure and Maxwell theory, Javadpour (2009) presented the AGP model with several empirical factors. Singh and Javadpour (2016) developed the AGP model by linear superposition of viscous flow and Knudsen diffusion. In addition, although Wang et al. (2016b) derived the AGP model for real gas transport in nanopores, the Knudsen diffusion is ignored. Song et al. (2016) introduced the phase behavior and presented the AGP models for IM and OM, which ignored the structural parameters in the surface diffusion equation. Weighting coefficients relevant to Kn of different mechanisms are introduced into the AGP models by Shi et al. (2013) and Wu et al. (2015aWu et al. ( , b, c, 2016Wu et al. ( , 2017. Few work (Shi et al. 2013;Zhang et al. 2017) focused on the effects of the OM porosity on the AGP of shale reservoirs. The comparison of different AGP models is shown in Table 1. In brief, an accurate AGP model is essential for the analysis of gas transport behavior in shale reservoirs.
Different from the previous work (Zhang et al. 2018), we modified a unified AGP model and focused on the comparison of different characterizations for the gas transport mechanisms, and more influences of various parameters on gas transport behavior are analyzed in this work. Viscous-slip flow, Knudsen diffusion, monolayer and multilayer gas adsorption, surface diffusion, the weighting coefficients of viscous-slip flow and Knudsen diffusion, and stress dependence are included in nanopores. The presented AGP model is validated by the experimental results. Then, on the one hand, the influences of Knudsen number, pore pressure, pore temperature, pore radius, OM porosity, and stress dependence on AGP are analyzed. On the other hand, the contributions of different mechanisms to the total AGP are conducted to reveal gas transport behavior. The proposed AGP model in this work is expected for macrosimulation for the development of shale reservoirs.

AGP mathematical model
The schematic diagram of gas transport in IM and OM nanopores is shown in Sects. 2.1 and 2.2. Some assumptions are as follows: single phase and isothermal methane transport  in nanopores; there is no mass transfer between organic and inorganic pores; gas adsorption on pore wall surface follows BET isotherm equation in OM; porosity and tortuosity, surface diffusion, and OM porosity are considered. Then, we derived the AGP models.

Characterization of gas transport in IM
The real gas effect, viscous-slip flow, and Knudsen diffusion are considered in IM nanopores as shown in Fig. 2. The Knudsen number Kn(IM) for gas transport in IM can be expressed as (Bird 1994;Michel et al. 2011) where r eI is the effective inorganic pore radius, m; κ B is the Boltzmann constant, MPa m 3 /K; T is the pore temperature, K; Z is the gas deviation factor, dimensionless; p is the pore pressure, MPa; and δ is molecular collision diameter, m. For real gas, Z is related to p pr = p p c and T pr = T T c can be expressed as (Mahmoud 2014) where p pr is the pseudo-reduced pressure, dimensionless; p c is the critical pressure, MPa; T pr is the pseudo-reduced temperature, K; and T c is the critical temperature, K. Lee et al. (1966) gave the expression of gas viscosity as where μ is the gas viscosity, MPa s; ρ is gas density, kg/m 3 ; M is the gas molar mass, kg/mol; and R is the gas universal constant, J/(mol K). Because of the collision between gas molecules, the viscous-slip flow for gas transport in IM nanopores can be written as (Karniadakis et al. 2005;Wu et al. 2016) (1) 5T pr − 5.524p pr e −2.5T pr + 0.044T 2 pr − 0.164T pr + 1.15 (3) = 10 −13 (9.379 + 0.01607M)T , where J I v+s is the viscous and slip flow mass flux in IM nanopores, kg/(m 2 s); I e is the effective porosity of IM, dimensionless; τ is the tortuosity, dimensionless; μ e is the effective gas viscosity, MPa s; b is the slip coefficient, dimensionless (here, b = −1 ); and ∇p is the gas pressure gradient, MPa/m. C o n s i d e r i n g t h e r a r e f i e d g a s e f fe c t , i f = 128 15π 2 tan −1 4.0Kn 0.4 , Eq. (4) becomes (Beskok and Karniadakis 1999) Considering the collision between gas molecules and pore walls, Knudsen diffusion for real gas is given by (Darabi et al. 2012) where J I K is the Knudsen diffusion mass flux in IM, kg/ (m 2 s); C g = 1 p − 1 Z dZ dp ; D m is the methane molecule diameter, m; for characterizing various transport mechanism weight on gas transport, weighting coefficients of viscous-slip flow, and Knudsen diffusion were presented (Shi et al. 2013;Wu et al. 2016). Thus, the total mass flux of free gas in IM pores can be expressed as where J I is the mass flux in IM nanopores kg/(m 2 s).
Two terms of the weighting coefficients based on the Knudsen number Kn reported (Shi et al. 2013;Wu et al. 2016) are listed as follows: where ω v+s is the weighting coefficient of viscous and slip flow, dimensionless; ω K is the weighting coefficient of the Knudsen number, dimensionless; and n = 5 , Kn 0.5 = 4.5 (Shi et al. 2013). Figure 3 illustrates the relationship between weighting coefficients and Knudsen number Kn.

Characterization of gas transport in OM
For adsorbed gas coexisting in pores, gas adsorption and surface diffusion are included as shown in Fig. 4.
The mass flux of free gas transport in OM can be expressed as where J O b is the bulk gas mass flux in OM, kg/(m 2 s); and O e is the effective porosity of OM, dimensionless. Both monolayer adsorption and multilayer adsorption have been found in shale reservoir conditions, and the adsorbed gas volume can be calculated by isotherm equation. Then, we can (10) obtain the effective pore radius of the nanocapillary tube for free gas transport in OM. Thus, the adsorption volume described by the BET isotherm equation is written as where V is the adsorption volume of real gas in OM nanopores, m 3 /kg; V L is the maximum adsorbed gas volume, m 3 /kg; C is a constant related to the net heat of adsorption, dimensionless; p 0 is the saturated adsorption pressure of the gas, MPa; and n is the adsorbed layer number, dimensionless. When n = 1 , Eq. (12) is simplified as where p L is the Langmuir pressure, MPa, Eq. (13) becomes the Langmuir isotherm. If V∕ V L = , the effective pore radius as shown in Fig. 4 can be expressed as (detailed derivation presented in our previous work (Zhang et al. 2018)): where r 0O is the original radius of organic nano-capillary tube, m; and for monolayer adsorption, the pore radius except first adsorption layer r ′ O is defined as (Xiong et al. 2012) where the gas coverage = p∕ Z p L +p∕ Z . The variation of the Knudsen number ratio Kn(OM)∕ Kn(IM) with pore radius and pore pressure is shown in Fig. 5. The related parameters are set as p L = 4 MPa, C = 20 , p 0 = 80 MPa, n = 2.
As shown in Fig. 5, in the organic nanopore, when the pore pressure is less than 1 MPa, Kn(OM)∕ Kn(IM) ≈ 1 . It suggests that the decrease in the pore radius caused by adsorbed gas has a slight effect on gas transport behavior. When the pore pressure is greater than 1 MPa, Kn(OM)∕ Kn(IM) increases with the increase in pore pressure and becomes larger than 1. Moreover, with the increase in pore radius, Kn(OM)∕ Kn(IM) decreases at the same pore pressure. Xiong et al. (2012) model just considered monolayer adsorption by using the pore radius r ′ O given by Eq. (16). Therefore, the results are greater than those calculated by the proposed model using the effective pore radius r eO when p < 20 MPa. Meanwhile, Kn(OM)∕ Kn(IM) of multilayer adsorption is smaller than that of monolayer adsorption when p < 20 MPa. When p > 20 MPa, Kn(OM)∕ Kn(IM) of multilayer adsorption increases greatly and is much larger than that of monolayer adsorption. Additionally, as the pore radius increases, Kn(OM)∕ Kn(IM) decreases, and especially, when the pore radius rises to 100 nm, Kn(OM)∕ Kn(IM) ≈ 1 , meaning the effects of adsorption on gas transport behaviors are negligible.
The surface diffusion flux for a real gas can be expressed as where J O s is the mass flux of surface diffusion in OM nanopores, kg/(m 2 s); D s is the surface diffusion coefficient, m 2 /s; and ρ r is the rock density, kg/m 3 .
The surface diffusion coefficient can be derived by the kinetic method under a high-pressure condition, and Chen and Yang (1991) presented the basic form as follows where D 0 s is surface diffusion coefficient when Γ = 0, m 2 /s. From Eqs. (18) and (19), it can be found that if b > m , > 1 and H(1 − ) = 0 causing surface diffusion stop. When b < m , < 1 , and H(1 − ) = 1 , surface diffusion occurs. Figure 6 shows the relationship among D s D 0 s , Γ, and κ. With the increase in gas coverage Γ, D s D 0 s increases at the same κ. When < 0.2 , κ influences D s D 0 s slightly. However, when > 0.2 , smaller κ causes a larger D s D 0 s at the same Γ, as shown in Fig. 6a. Besides, when > 0.4 , as κ increases, D s D 0 s gradually decreases, as shown in Fig. 6b. Particularly, when = 1 , D s D 0 s reduces sharply with the rising of κ. When < 0.4 , as κ increases, D s D 0 s increases slightly.
Combining Eqs. (11) and (16), the total mass flux of a real gas transport in OM can be expressed as where J O is the total mass flux of OM, kg/(m 2 s).
In the parallel connection of organic and inorganic pores, the total gas mass flux can be expressed as where J I + O is the total mass flux of a nanoporous medium, kg/(m 2 s).
We assume that the effective porosity of IM is I e , and the effective porosity of bulk gas in OM is O e . The gas mass flux in a porous medium can be obtained by Darcy's law as where k is the permeability, m 2 .
The AGP model in this work by comparing Eqs. (21) with (22) can be written as where k a is the apparent permeability of nanopores, m 2 ; k O a is the apparent permeability of OM nanopores, m 2 ; and k I a is the apparent permeability of IM nanopores, m 2 .
The AGP of OM and IM is given by, respectively, The pore radius considering the stress dependence can be given by (detailed derivation in "Appendix"): where r 0 is the original radius of the nanopore, m; p e is the effective stress, MPa; p′ is the atmospheric pressure MPa; q is the porosity coefficient, dimensionless; and s is the permeability coefficient, dimensionless.

Model verification
The experimental results (Tison 1993) and results obtained from the linearized Boltzmann (Loyalka and Hamoodi 1990) and Song et al. model (Song et al. 2016) were used to validate the presented model in this study. The validation indicates that the proposed AGP model in this work is reliable in describing gas transport behavior in nanopores, especially when the Knudsen number is less than 0.2 (Fig. 7).

Influences of pore radius and pressure on gas transport behavior
The influences of various parameters can be analyzed based on the parameters given in Table 2. The permeability ratio k I a k I 0 decreases as the pore radius increases when the pore pressure is 0.1-30 MPa, as shown in Fig. 8a. When r > 5 μm, k I a k I 0 ≈ 1 at any pore pressure    Fig. 8 The AGP of IM nanopores versus pore radius greater than 0.1 MPa. Then, the larger the pore pressure is, the smaller the pore radius is when k I a k I 0 ≈ 1 . When r < 10 nm, k I a k I 0 is larger than 1 at any pore pressure less than 30 MPa. Meanwhile, as shown in Fig. 8b, the AGP of IM increases with an increase in pore radius and then is equal to the intrinsic permeability at a certain pore radius. Moreover, a smaller pore pressure leads to a larger AGP at the same pore radius condition.
As shown in Fig. 9a, the surface diffusion coefficient D s affects the gas transport behavior when r < 10 nm and pore pressure is equal to 10 MPa. The larger surface diffusion coefficient D s , the greater k O a k O 0 is under the same pore radius condition. However, when the surface diffusion coefficient D s is less than 1 × 10 −5 m 2 /s, it has a little effect on the AGP. Figure 9b shows that when p = 0.1 MPa or r < 10 nm, the effect of D s can be ignored. In addition, under the same pore pressure condition, the AGP of OM nanopores increases with an increase in D s .

Influences of OM porosity, temperature, and stress dependence on gas transport behavior
The OM porosity influences gas transport behavior in nanoporous shale through the gas adsorption and surface diffusion. If the surface diffusion coefficient D s is set as 1 × 10 −4 m 2 /s, and the pore pressure is 1 MPa, the ratio of total AGP to intrinsic permeability of the nanoporous shale k a k 0 decreases slightly with the OM porosity increase when r < 5 nm. This suggests that the gas adsorption and small surface diffusion coefficient cause the AGP of OM smaller than that of IM, as shown in Fig. 10. The effects of the OM porosity and surface diffusion coefficient on gas transport behavior are analyzed in the next section.
As shown in Fig. 11a, there are slight differences k a k 0 under four different temperature conditions when r < 10 nm and p = 1 MPa. With an increase in temperature, k a k 0 increases. Furthermore, with an increase in pressure, the weighting coefficient of the Knudsen diffusion decreases, which leads to a decrease in k a k 0 . This suggests that the influence of temperature on gas transport behavior should be considered under small pore radius conditions. With the development of shale gas reservoirs, the stress dependence will occur in nanopores. Hence, different porosity and permeability coefficients of OM and IM nanopores are set to discuss the influence of stress dependence on the total AGP, as shown in Fig. 11b. When r < 40 nm, the total AGP is larger than the intrinsic permeability. This is because the total AGP is mainly affected by the Knudsen diffusion and surface diffusion. Moreover, considering the stress dependence, the total AGP is smaller than that without consideration of stress dependence. Additionally, when the porosity and permeability coefficients of IM nanopores q I > q O and s I > s O , the total AGP is smaller than that if q I < q O and s I < s O . This is because the stress dependence not only affects the gas viscous flow and slip flow in nanopores but also the Knudsen diffusion and surface diffusion. To a certain extent, the total AGP is affected by the surface diffusion coefficient and OM porosity.

Contributions of different mechanisms to the total AGP
Based on the above derivation, the ratio of each mechanism to the total AGP can be obtained. k O s is the apparent permeability of surface diffusion in OM, m 2 ; k I + O K is the apparent permeability of Knudsen diffusion in nanopores, m 2 ; and k I + O v + s is the apparent permeability of viscous and slip flow in nanopores, m 2 . The surface diffusion coefficient has a significant influence on the contribution of OM AGP to the total AGP. The influence of the OM porosity on gas transport behavior is discussed in Fig. 12. When D s = 1 × 10 −4 m 2 /s, k O a k a increases as the OM porosity increases, and the surface diffusion dominates the OM AGP. If r = 1 nm, OM porosity O e = 0.03 and IM porosity I e = 0.02 , k O a k a ≈ 85% when p = 5 MPa, and k O s k a is about 70%. Therefore, when the pore radius is small, the gas transport in OM nanopores dominates the total AGP under the large surface diffusion coefficient and high OM porosity. Figure 13 shows the contributions for three mechanisms and two pore types versus pore pressure under different pore radii (2 and 5 nm) and different pore pressures (1 and 20 MPa) when D s = 1 × 10 −4 m 2 /s, OM porosity O e = 0.01 , and IM OM porosity/IM porosity = 20% Fig. 12 Variation of the contributions of surface diffusion and OM to the total AGP with pore pressure under different r and D s porosity I e = 0.04 . In Fig. 13a, for pore radius of 2 nm, the increasing pore pressure decreases the Knudsen diffusion ratio k I + O K k a and increases the viscous-slip flow ratio k I + O v + s k a . The viscous-slip flow dominates gas transport when p > 20 MPa. For the pore radius of 5 nm, the viscous-slip flow is dominant at a smaller pore pressure (10 MPa) than that in the pore with 2 nm radius in Fig. 13b.
When the pore pressure is 1 MPa, the Knudsen diffusion still dominates in the pores with a radius smaller than 10 nm. k I + O K k a increases first and decreases later with an increase in pore radius, as shown in Fig. 13c. When the pore pressure continues to increase, as shown in Fig. 13d, the peak of k I + O K k a drops, and the surface diffusion, Knudsen diffusion, and viscous-slip diffusion collectively govern gas transport behavior in pores with a radius smaller than 10 nm. Moreover,   the pore radius range in which the viscous-slip flow dominates is larger with an increase in pore pressure. The influences of surface diffusion and Knudsen diffusion on gas transport can be neglected in the pores with a radius smaller than 10 nm.

Conclusions
In this study, a unified gas transport model in the parallel connection of inorganic and organic shale nanopores was developed. The model was validated by the experimental results. According to the discussion of the results, the following conclusions can be drawn: 1. The parallel connection of inorganic and organic nanopores is introduced into the AGP model to characterize the coexistence of different pore types in shale gas reservoirs. Otherwise, weighting coefficients based on Knudsen number are considered. The effect of adsorbed gas on Kn(OM)∕ Kn(IM) can be neglected when p < 1 MPa or r > 100 nm. When 1 MPa < p < 20 MPa , Kn(OM)∕ Kn(IM) caused by multilayer adsorption is lowest and then becomes largest if p > 40 MPa. 2. By validating the proposed model with experimental results, the weighting coefficients and surface diffusion equation can be accepted. The proposed model can be used to analyze different gas transport mechanisms in nanoporous shale and microsimulation in the shale gas reservoirs. 3. The pore radius and pressure have significant influences on the AGP. The larger pore radius or the lager pressure causes the smaller k a k 0 , meaning the viscous flow dominates gas transport. The surface diffusion has a large effect on gas transport behavior in nanopores with a radius less than 10 nm when D s > 1 × 10 −5 m 2 /s. The OM porosity affects k a k 0 slightly when the pore radius is larger than 5 nm. Additionally, the temperature has a small effect on gas transport behavior. 4. When the pore radius is the same, the larger the OM porosity is, the greater the value of k O s k a is. For a small pore and a low pressure, the surface diffusion and Knudsen diffusion are dominant. When the pore is larger and pressure is higher, gas transport behavior is controlled by viscous flow.
Considering the complex fracture network around the horizontal well in shale reservoirs (Zeng et al. 2016(Zeng et al. , 2018, the presented model can be coupled with the natural fracture system and hydraulic fracture network to simulate the field scale development.