Gas Transport Model in Organic Shale Nanopores Considering Langmuir Slip Conditions and Diffusion: Pore Confinement, Real Gas, and Geomechanical Effects

Abstract

Nanopores are extremely developed and randomly distributed in shale gas reservoirs. Due to the rarefied conditions in shale strata, multiple gas transport mechanisms coexist and need further understanding. The commonly used slip models are mostly based on Maxwell slip boundary condition, which assumes elastic collisions between gas molecules and solid surfaces. However, gas molecules do not rebound from solid surfaces elastically, but rather are adsorbed on them and then re-emitted after some time lag. A Langmuir slip permeability model was established by introducing Langmuir slip BC. Knudsen diffusion of bulk phase gas and surface diffusion of adsorbed gas were also coupled into our nanopore transport model. Considering the effects of real gas, stress dependence, thermodynamic phase changes due to pore confinement, surface roughness, gas molecular volume, and pore enlargement due to gas desorption during depressurization, a unified gas transport model in organic shale nanopores was established, which was then upscaled by coupling effective porosity and tortuosity to describe practical SGR properties. The bulk phase transport model, single capillary model, and upscaled porous media model were validated by data from experimental data, lattice Boltzmann method or model comparisons. Based on the new gas transport model, the equivalent permeability of different flow mechanisms as well as the flux proportion of each mechanism to total flow rate was investigated in different pore radius and pressure conditions. The study in this paper revealed special gas transport characteristics in shale nonopores and provided a robust foundation for accurate simulation of shale gas production.

1. Introduction

In recent decades, shale gas has drawn great attention globally, especially with its successful development in North America [1,2,3]. As a typical unconventional energy resource, shale gas transport mechanisms in organic shale nanopores display a number of unique characteristics; for example, the slippage effect, Knudsen diffusion, surface diffusion, and adsorption/desorption behaviors [4,5]. These unique transport mechanisms enhance the flow capability of shale gas and promote commercial production in shale gas reservoirs. Generally, Shale gas reservoir (SGR) is tight or ultratight formations, with pore sizes ranging from 1−200 nm [6] and matrix intrinsic permeability from 10⁻⁹–10⁻³ mD [7,8]. The stress dependence of shale properties, for example, porosity and permeability, complicates the study even further [9]. Thus, researching gas transport mechanisms in organic shale nanopores is a very complex and challenging task for petroleum engineers [10,11].

Gas adsorption is a special storage form in SGRs, which can account for up to 85% of the original gas in place due to the very large surface area of a shale formation [4]. On one hand, adsorbed gas increases the recoverable shale gas reserves, because it will desorb from pore surfaces when pressure decreases [12,13,14,15]; on the other, adsorbed gas occupies a certain amount of flow paths and changes the effective hydraulic radius, which cannot be ignored for nanometer-scale transport studies [10,11]. Surface diffusion of adsorbed gas occurs under a chemical potential energy gradient in shale nanopores [16] and facilitates shale gas production. Considering the high reservoir pressure and surface heterogeneity, a surface diffusion model of adsorbed gas in nanopores of SGRs was established and analyzed by Wu et al. [17]. A new formulation of apparent permeability (AP) for gas transport in SGRs was proposed by Zhang et al. [18], which took surface diffusion into account. Wang et al. [19] established an AP model for gas transport in shale nanopores considering multiple effects, such as non-Darcy flow and surface diffusion.

Viscous flow and Knudsen diffusion of free gas are two important bulk phase gas transport mechanisms. Due to the breakdown of the continuum hypothesis in nanopores, the viscous seepage theory must be modified to satisfy slip boundaries. Klinkenberg [20] devised an empirical formulation to account for the slippage effect, which linked slippage permeability to intrinsic permeability, average pressure, and an empirical parameter b_k. After that, much effort has been made to obtain this parameter b_k theoretically or empirically [21,22,23,24,25,26]. Basing their work on the dynamic slippage concept, Ertekin et al. [23] established a new AP model considering the slippage effect. Skjetne and Auriault [27] made a further study on the gas slippage flow phenomenon in porous media (PM) on the basis of the Navier–Stokes equation.

Another type of slip model is based on slip boundary conditions (BCs), including the first and second-order modifications. The Maxwell slip model [28] is a first-order approximation from kinetic theory, based on which many other permeability models [25,26,29] have been established for tight and shale gas reservoirs. Second-order slip velocity models have been successively proposed to improve calculation precision. Hisa et al. [30] investigated molecular rarefaction effects through experiments and proposed a new slip BC that added a second-order term to the Maxwell model [28]. Basing their work on high-order slip BCs and asymptotic analysis, Beskok and Karniadakis [31] proposed another widely used second-order slip velocity model. Considering the location difference between slip surface and solid wall, Zhang et al. [32] improved the Beskok–Karniadakis (BK) slip model by introducing a new parameter γ to describe the distance between the slip surface and the solid wall. Niu et al. [33] made a new second-order gas-permeability correlation for shale gas slip flow based on the study of Zhang et al. [32]. In summary, slip velocity models can uniformly be expressed as:

$$u_{\mathrm{slip}}=C_1\lambda {(\frac{\partial u}{\partial n})}_{\mathrm{s}}+C_2{\lambda }^2{(\frac{{\partial }^{2}u}{\partial n^2})}_{\mathrm{s}}$$

(1)

where u_slip is the velocity at the slip surface; λ is average mean free path of gas molecules; and C₁ and C₂ are two slip coefficients, the values of which are summarized in

Table 1. Slip coefficients in different slip boundary conditions.

Reference	C1	C2
Maxwell [28]	(2 − σv)/σv	0
Hsia et al. [30]	(2 − σv)/σv	−1/2
Hwang et al. [34]	0.01807	−0.676775(4/π)0.58734 Kn−0.82532
Beskok et al. [31]	1/(1 − b Kn)	0
Aubert et al. [35]	(2 − σv)/σv	−9/8
Bahukudumbi et al. [36]	[1.2977 + 0.71851tan−1 (−1.17488 Kn0.58642)] (2 − σv)/σv	0
Hadjiconstantinou [37]	1.1466	−0.647
Ng et al. [38]	1.15	−0.25 Kn−0.65
Zhang HW et al. [32]	[1 − (1 − γ)σv]/σv	(1 − γ)[1 − (1 − γ)σv]/(2σv)

Note: σ_v is tangential momentum accommodation coefficient; Kn is Knudsen number; b is slip coefficient in B–K model [31]; γ is a parameter in Zhang’s model to describe slip surface location.

.

The values of some parameters in those slip velocity models, for example, the tangential momentum accommodation coefficient σ_v, can be very difficult to determine [39,40,41]. Although that coefficient is defined as the percentage of molecules reflected diffusively from solid walls, its physical meaning is still not very straightforward. The value has to be obtained by experiments or molecular dynamics (MD), which are also not easy tasks [41]. Besides, the slip velocity derived from those models may not be bounded, depending on specific situations [39,40], which may cause severe problems in numerical applications [42]. Except for Zhang et al. [32] and Niu et al. [33], all the studies in Table 1 ignored the distance between the slip surface and the solid wall. Meanwhile, in those studies, gas molecules through a slip surface are assumed to be comprised of molecules from solid wall and the surface λ away from the solid wall, which account for the same percentage of 50% (that is, n_λ = n_w = 0.5 n_s). This assumption is also invalid when gas molecules are not uniformly distributed across the pipe section. The uneven molecule distribution, which can be attributed to the existence of gas molecular accumulation or adsorption near solid walls, has been proved by MD [32,33].

In those slip velocity models, the content of molecules diffusively reflected from solid walls is artificially assumed to be σ_v, whereas (1 − σ_v) is the ratio of molecules specularly reflected to incident molecules. However, in reality gas molecules do not rebound elastically from solid walls after colliding [39,40,41,43,44,45], due to the existence of attractive forces between gas and solid molecules. Gas molecules can be adsorbed on solid surfaces and then desorbed when the pressure decreases, causing an inelastic rebound. Therefore, a Langmuir slip BC more accurately describes the real slip phenomenon from the perspective of a physical meaning. Knudsen diffusion and surface diffusion are important gas transport mechanisms in shale nanopores. Although the Langmuir slip permeability (LSP) model for SGRs was first studied by Singh et al. [41], the results are not very convenient to apply. Besides, Knudsen/surface diffusion, real gas effects, geomechanical effects, and thermodynamics phase changes were also not considered in their research. Song et al. [46] pointed out that slippage flow, Knudsen diffusion, and surface diffusion were three gas transfer mechanisms in nanometer pores. However, the effects of slippage flow and Knudsen diffusion were broadly coupled together by the BK model [31] in their research, which ignored gas accumulation near pore walls and failed to distinguish the relative intensity of slippage flow and Knudsen diffusion in different conditions. Besides, the surface diffusion model in their study did not consider the effects of PM properties.

In this paper, considering gas-solid interaction, the slip phenomenon is described by the Langmuir slip BC [39,40,41], and a new slip model for bulk phase gas transport was established. Based our work on the LSP model, we also propose a new gas transport model coupling viscous flow, slippage effect, Knudsen diffusion, surface diffusion, and adsorption/desorption. Also taken into consideration were the real gas effect, stress dependence of physical properties, thermodynamics phase changes due to pore confinement, and dynamic adsorption layer thickness. The advantages of our LSP model over previous Maxwell-based ones are also illustrated. To make our gas transport model more applicable to practical situations, we re-expressed it in an AP form and upscaled it from a single capillary model to a porous media form, which coupled porosity and tortuosity into AP model and considered the effects of adsorbed gas volume on effective porosity.

2. Flow Mechanisms and Their Description in Nanopores

The three main gas storage forms in SGRs are free gas, adsorbed gas, and dissolved gas. Due to the nanometer-scale effect, bulk phase free gas can be transferred by slippage viscous flow under pressure gradient and by Knudsen diffusion under a concentration field. Meanwhile, adsorbed gas can move from one site to another to achieve mass transport on pore surface by surface diffusion [47]. The typical gas storage and transport form in shale nanopores is shown in Figure 1.

2.1. Stress-Dependent Darcy Permeability

As is known, the intrinsic permeability of a nanometer capillary with a circular section can be given [19,46,48] as:

$$k_{\mathrm{ins}}=\frac{r^2}{8}$$

(2)

where r is the radius of the capillary tube.

To take stress dependence into account, we adopt Gangi’s model [46,49,50] to describe the permeability change under effective stress:

$$k_{\mathrm{ins}\mathsf{\sigma }}=k_0\frac{{[1-{(\frac{\sigma }{p_1})}^m]}^3}{{[1-{(\frac{{\sigma }_0}{p_1})}^m]}^3},$$

(3)

where k₀ is the formation permeability at original state, m²; p₁ and m are parameters related to the pore stiffness and roughness respectively; σ₀ is effective stress of initial state, Pa; σ is the effective stress during production, which can be obtained by:

$$\sigma =p_{\mathrm{s}}-\alpha p$$

(4)

where p_s is overburden formation pressure, Pa; α is Biot constant, dimensionless; p is reservoir pressure, Pa.

Considering the geomechanical impacts of formation compaction during reservoir depletion, we have the intrinsic permeability expression as follows by coupling Equations (2)–(4):

$$k_{\mathrm{ins}\sigma }=\frac{r_0^2}{8}\frac{{[1-{(\frac{p_{\mathrm{c}}-\alpha p}{p_1})}^m]}^3}{{[1-{(\frac{p_{\mathrm{c}}-\alpha p_0}{p_1})}^m]}^3}.$$

(5)

2.2. LSP Modeling

The adsorption and desorption behaviors of shale gas have been studied in petroleum industry for a long time. However, the interfacial interaction between gas molecules and solid surfaces has failed to attract wide attention in the industry. Based our work on the Langmuir slip BC [39,40,41], we will establish a new LSP model considering gas molecule residence on solid surfaces when the gas-solid collision happens.

The Langmuir adsorption isotherm can be expressed in the following form [44]:

$$\theta =\frac{bp}{1+bp}$$

(6)

where θ is the gas coverage ratio on a solid surface, dimensionless, and b is the Langmuir equilibrium adsorption constant, Pa⁻¹, which is usually expressed as the ratio of equilibrium adsorption rate to equilibrium desorption rate.

To unify our model with previous ones, the following expression of parameter b can be adopted [44]:

$$b=\frac{1}{4\omega Kn}\frac{1}{p_{\mathrm{r}}}$$

(7)

where ω is a parameter related to gas-solid interaction; p_r is the reference pressure, Pa; and Kn is the Knudsen number, dimensionless, which can be written as the ratio of gas mean free path to characteristic length [11,51] as:

$$Kn=\frac{\lambda }{R_{\mathrm{h}}}$$

(8)

The slip velocity BC could be expressed in the following form [39,40,41,52] once the surface coverage ratio is calculated:

$$u_{\mathrm{s}}=\theta u_{\mathrm{w}}+(1-\theta )u_{\mathsf{\lambda }}$$

(9)

where u_s is the slip surface velocity, m²/s; u_w is the velocity of solid wall, m²/s; and u_λ is gas velocity at the surface a mean free path away from the wall, m²/s.

In SGRs, the velocity of the wall can be presumed as zero. Moreover, the reference velocity can also be chosen as the midstream value, the reason for which can be found in the study of Eu [53].

Therefore, the Langmuir slip BC can be obtained by coupling Equations (6) and (9):

$$u_{\mathrm{s}}=\frac{1}{1+bp}u(x,r=0)$$

(10)

The Navier–Stokes equation can be adopted to describe the viscous incompressible Newtonian fluid flow [32,33] as:

$$\rho (\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}\cdot \nabla \overrightarrow{u})=-\nabla p+\mu {\nabla }^2\overrightarrow{u}+\rho \overrightarrow{g}$$

(11)

where ρ is the fluid density, kg/m³; $\overrightarrow{u}$ is the velocity field; μ is gas viscosity, Pa·s; and $\rho \overrightarrow{g}$ is the gravity force.

Usually, for fluid flow in a nanometer circular capillary, the influence of $\rho \overrightarrow{g}$ is negligible. Besides, the value of $\partial \overrightarrow{u}/\partial t$ can be set to zero when the flow is under steady state. Taking the cylindrical coordinate system (see Figure 2) and ignoring the term of inertia force and mass force, the following Stokes equation is adopted to describe steady state Hagen–Poiseuille fluid flow [32,33,54]:

$$\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}(r\frac{\mathrm{d}u}{\mathrm{d}r})=\frac{1}{\mu }\frac{\mathrm{d}p}{\mathrm{d}x}$$

(12)

According to practical situations, the fluid velocity in nanopores should be a finite value. Thus, the following definite solution condition could be given:

$$u{|}_{r=0}<+\infty$$

(13)

By solving the Hagen–Poiseuille flow Equation (12), with the definite solution condition of Equations (10) and (13), the fluid velocity distribution in a circular tube could be obtained:

$$u(r)=-\frac{R^2}{4\mu }\frac{\mathrm{dp}}{\mathrm{d}x}[1+\frac{4\omega p_{\mathrm{r}}}{p}Kn-{(\frac{r}{R})}^2]$$

(14)

Integrating the velocity distribution u over the radius r, the flow rate can be calculated

$$Q={\displaystyle \underset{0}{\overset{R}{\int }}2\mathsf{\pi }r\cdot u(r)\mathrm{d}r}$$

(15)

Based on the expression of flow rate Q, we can get the permeability modification coefficient F according to the research of Sakhaee-Pour and Bryant [54] as:

$$F=\frac{Q(Kn)}{Q(Kn\to 0)}=1+\frac{8\omega p_{\mathrm{r}}}{p}Kn.$$

(16)

Up to now, the establishment of the LSP model has been completed, which considers slippage effect, molecular residence on solid wall after collisions, and geomechanical effects of stress dependence during depressurization:

$$k_{\mathrm{slip}}=F(p,Kn)k_{\mathrm{ins}\sigma }$$

(17)

2.3. Knudsen Diffusion

Knudsen diffusion is another important transport mechanism in nanometer pores [55,56,57,58]. Considering solid surface roughness [18,59], we can express Knudsen diffusion in the following form:

$$J_{\mathrm{k}}=-{\delta }^{D_{\mathrm{f}}-2}M{D}_{\mathrm{k}}\cdot \nabla c_{\mathrm{g}}$$

(18)

where δ is the pore surface roughness, dimensionless; D_f is fractal dimension of pore surface, dimensionless; M is molar mass, kg/mol; J_k is Knudsen mass flux, kg/(m²·s); c_g is gas concentration, mol/m³; and D_k is the Knudsen diffusion coefficient, m²/s, the value of which is shown in

Table 2. Value of Knudsen diffusion coefficient in different models.

Number	Dk Expression	Reference
1	Dk = 2r[8RT/(πM)]0.5/3	Igwe [56]
2	Dk = 31.54(kins)0.67/M0.5	Ertekin et al. [23]
3	Dk = 9700r(T/M)0.5	Ruthven [60]
4	Dk = 4σ[πkinsφRT/(2M)]0.5	Florenceet al. [24]

Note: T is temperature, K; φ is porosity.

.

According to the gas state equation, gas concentration C can be expressed with the consideration of real gas effect as:

$$c_{\mathrm{g}}=\frac{p}{Z{R}_{\mathrm{cst}}T}$$

(19)

where R_cst is the gas constant, 8.314 J/(mol·K).

Finally, the Knudsen diffusion flux can be expressed in terms of pressure gradient as:

$$J_{\mathrm{k}}=-{\delta }^{D_{\mathrm{f}}-2}\frac{M{D}_{\mathrm{k}}}{RT}\cdot \frac{p{c}_{\mathrm{g}}}{Z}\nabla p.$$

(20)

2.4. Surface Diffusion

As we have introduced in Section 2.3, gas molecules could adsorb on solid surfaces and move under a chemical potential gradient [61]. Surface diffusion mass flux can be expressed [10,11,17,59] as:

$$J_{\mathrm{s}}=-\frac{D_{\mathrm{s}}{C}_{\mathrm{s}}M}{p}\nabla p$$

(21)

where D_s is the surface diffusion coefficient, m²/s, and C_s is the concentration of adsorbed gas, mol/m³.

Based on the research of Hwang and Kammermeyer [34] and combing the methane/carbon adsorption experiment data, the surface diffusion coefficient for methane could be given [11,16,52] as:

$$D_{\mathrm{s}}=8.29\times {10}^{-7}{T}^{0.5}\exp (-\frac{\Delta H^{0.8}}{RT})$$

(22)

where ∆H is the equivalent adsorption heat, J/mol.

According to the Langmuir adsorption theory, the adsorbed gas concentration can be expressed as:

$$C_{\mathrm{s}}=\frac{{\rho }_{\mathrm{gsc}}{G}_L{\theta }_{\mathrm{real}}}{M}$$

(23)

where G_L is the Langmuir volume, m³/m³, and ρ_gsc is gas density at the standard condition, kg/m³.

Note that the expression of gas concentration C_s in Equation (23) includes the influence of the gas coverage ratio. Thus, there is no need to modify diffusion coefficient D_s when pressure changes, which differs from previous studies [10,11,17,59] and is very convenient in application.

In Equation (23), θ_real is the gas coverage ratio on solid surfaces considering the real gas effect and thermodynamics phase changes, which can be expressed as:

$${\theta }_{\mathrm{real}}=\frac{bp/Z}{1+bp/Z}$$

(24)

2.5. Pore Confinement and the Real Gas Effect

A gas deviation factor Z is adopted to account for the real gas effect and thermodynamic phase changes of nanopore fluids due to the pore confinement effect. To calculate this parameter, critical pressure and temperature changes in nanometer-scale pores should be calculated first, which can be obtained by a modified van der Waals equation proposed by Islam et al. [46] in the following form [62]:

$$T_{\mathrm{c}}=\frac{8}{27b{R}_{\mathrm{cst}}}[a-2{\zeta }^3\epsilon N^2\frac{\zeta }{r_{\sigma }}(2.6275-0.6743\frac{\zeta }{r_{\sigma }})]$$

(25)

$$p_{\mathrm{c}}=\frac{8}{27b^2}[a-2{\zeta }^3\epsilon N^2\frac{\zeta }{r_{\sigma }}(2.6275-0.6743\frac{\zeta }{r_{\sigma }})]$$

(26)

where T_c is pseudocritical temperature, K; p_c is pseudocritical pressure, MPa; a is the van der Waals energy parameter, 10⁻⁶ Pa·m⁶/mol²; b is the van der Waals volume parameter, 10⁻³ m³/mol; ζ is the Lennard–Jones size parameter, m; ε is the Lennard–Jones energy parameter, dimensionless; and N is Avogadro’s constant, 6.02 × 10²³.

Then, according to the corresponding state principle, pseudoreduced pressure and temperature are expressed as:

$$p_{\mathrm{pr}}=\frac{p}{p_{\mathrm{c}}}$$

(27)

and:

$$T_{\mathrm{pr}}=\frac{T}{T_{\mathrm{c}}}$$

(28)

where p_pr is the pseudoreduced pressure, dimensionless, the and T_pr is pseudoreduced temperature, dimensionless.

The gas deviation Z-factor can be solved out by the following correlation formula for high-pressure gases:

$$Z=(0.702\times e^{-2.5T_{\mathrm{pr}}})p_{\mathrm{pr}}^2-5.524\times e^{-2.5T_{\mathrm{pr}}}{p}_{\mathrm{pr}}+0.044\times T_{\mathrm{pr}}^2-0.164\times T_{\mathrm{pr}}+1.15$$

(29)

Gas viscosity can be obtained by the following formula proposed by Izadmehr et al. [63]:

$$\mu =a\times T_{\mathrm{pr}}+b\times p_{\mathrm{pr}}+c\times \sqrt{p_{\mathrm{pr}}}+d\times T_{\mathrm{pr}}^2+e\times \frac{p_{\mathrm{pr}}}{T_{\mathrm{pr}}}+f$$

(30)

The values of coefficients a through f are shown in

Table 3. Values of correlation coefficients for viscosity calculation.

Coefficient	Value
a	0.00850507486545010
b	−0.00104065426590739
c	−0.00217777225933512
d	−0.000510724061609292
e	0.00595154429253907
f	−0.000548942531453252

.

2.6. Effect of Adsorbed Layer on Pore Radii

The effect of gas molecule adsorption on nanometer-scale pore radii is not negligible, because the molecules are large enough to change the effective hydraulic radius. During reservoir production, adsorption has two effects on permeability: (1) the hindering effects due to its occupation of flow paths and (2) enhancement due to pore enlargement after desorption during depressurization. Coupling Equations (2) and (3), we have the expression of pore radius under the effective stress σ:

$$r_{\mathrm{ins}\sigma }=r_0{[\frac{1-{(\frac{\sigma }{p_1})}^m}{1-{(\frac{{\sigma }_0}{p_1})}^m}]}^{1.5}$$

(31)

Assuming a monolayer adsorption on solid surfaces, the effect of adsorbed gas on pore radii can be given [46,59] by:

$$r_{\mathrm{eff}}=r_{\mathrm{ins}\sigma }-d_{\mathrm{m}}{\theta }_{\mathrm{real}}$$

(32)

where r_eff is the effective hydraulic radius considering the effect of stress dependence and gas adsorption on solid surface, m.

In the end, the permeability considering pore stress dependence and adsorbed gas occupation on flow paths as well as slippage effect can be written as:

$$k_{\mathrm{eff}}=\frac{r_{\mathrm{eff}}^2}{8}$$

(33)

where k_eff represents the effective permeability considering stress dependence and adsorption layer thickness, m².

2.7. AP Modeling

To unify our model with continuum media flow equation of Darcy’s law, Knudsen and surface diffusion flux were also transformed into the form of equivalent permeability. Comparing corresponding diffusion flux expression with Darcy’s law, different types of equivalent permeability expressions could be obtained.

For slippage viscous flow, considering the effects of stress dependence and dynamic adsorption/desorption, the equivalent permeability could be expressed from the LSP model as:

$$k_{\mathrm{slip}}=F(Kn,p)\cdot \frac{r_{\mathrm{eff}}^2}{8}$$

(34)

In contrast to previous slip models, the slip coefficient F is not only a function of Knudsen number Kn, but also dependent on pore pressure.

For Knudsen diffusion, the equivalent permeability has the following form, where the Knudsen diffusion coefficient is determined by the method of Igwe [56]:

$$k_{\mathrm{knudsen}}={\delta }^{D_{\mathrm{f}}-2}{D}_{\mathrm{k}}{c}_{\mathrm{g}}\mu$$

(35)

Similarly, the equivalent permeability of surface diffusion could be expressed as:

$$k_{\mathrm{surface}}=D_{\mathrm{s}}{G}_{\mathrm{L}}{\theta }_{\mathrm{real}}\frac{\mu Z}{p^2}\frac{T{p}_{\mathrm{sc}}}{T_{\mathrm{sc}}}$$

(36)

According to previous research by Xiong et al. [64], Wasaki and Akkutlu [51], and Javadpour [55], the AP can be written as the sum of bulk transport and adsorbed gas transport. Therefore, the AP of a single nanoscale capillary can be written as:

$$k_{\mathrm{app}}=F(Kn,p)\cdot \frac{r_{\mathrm{eff}}^2}{8}+{\delta }^{D_{\mathrm{f}}-2}{D}_{\mathrm{k}}{c}_{\mathrm{g}}\mu +D_{\mathrm{s}}{G}_{\mathrm{L}}{\theta }_{\mathrm{real}}\frac{\mu Z}{p^2}\frac{T{p}_{\mathrm{sc}}}{T_{\mathrm{sc}}}$$

(37)

2.8. Model Correlation Considering PM Properties

The AP expression of Equation (37) represents a single capillary (Figure 3b), without considering PM properties and their geomechanical effects during depressurization. Thus, it needs to be upscaled to incorporate the effective porosity and tortuosity [65] (Figure 3a) to represent numerous nanopores in SGRs.

Considering the existence of adsorbed gas on solid surfaces, the effective radius and area for bulk phase gas transport are r_eff and A_eff respectively; the radius and cross-section area of the single capillary are r₁ and A₁ respectively. Meanwhile, the length of a single capillary is l₁ and the permeability is k₁. Correspondingly, the PM is assumed to be composed of n capillaries, with permeability k_t, cross section area A_t, and length l_t.

For bulk phase gas transport in n single capillaries without considering adsorbed gas, the volume flux can be given according to Darcy’s law using the equivalent permeability conception as:

$$q_{\mathrm{n}}=-n{A}_1\frac{k_1}{\mu }\frac{dp}{d{l}_1}.$$

(38)

The volume flux of UPM (upscaled porous media) model is:

$$q_{\mathrm{t}}=-A_{\mathrm{t}}\frac{k_{\mathrm{t}}}{\mu }\frac{dp}{d{l}_{\mathrm{t}}}$$

(39)

According to equivalent state, the volume flux calculated from Equations (38) and (39) should be the same, so we have:

$$k_{\mathrm{t}}=\frac{n{A}_1}{A_{\mathrm{t}}}\frac{d{l}_{\mathrm{t}}}{d{l}_1}{k}_1$$

(40)

From a physical meaning perspective, we have the following two expressions [66]

$$n{A}_1=A_{\mathrm{t}}\phi$$

(41)

$$\tau =\frac{l_1}{l_{\mathrm{t}}}$$

(42)

Coupling Equations (40)–(42), the relation between k₁ and k_t is obtained as:

$$k_{\mathrm{t}}=\frac{\phi }{\tau }{k}_1$$

(43)

Considering the effect of adsorbed gas on porosity, we have:

$${\phi }_{\mathrm{eff}}=\phi {(\frac{r_{\mathrm{eff}}}{r_1})}^2$$

(44)

Thus, the correction factor from single capillary model to UPM model for bulk phase gas is:

$${\gamma }_{\mathrm{b}}=\frac{\phi }{\tau }{(1-\frac{r_{\mathrm{ad}}}{r_1})}^2$$

(45)

In Equation (45), r_ad is the thickness of the adsorption layer, which can account for the occupation of adsorbed gas in flow paths and the pore enlargement due to gas desorption during pressure-depletion production:

$$r_{\mathrm{ad}}=d_{\mathrm{m}}{\theta }_{\mathrm{real}}$$

(46)

where d_m is the gas molecular diameter, m.

Assuming the cross section area for surface diffusion is A_s and taking similarly procedures as bulk gas transport for surface diffusion, the volume flux in n single capillaries is calculated as:

$$q_{\mathrm{ns}}=-n{A}_{\mathrm{s}}\frac{k_{\mathrm{s}}}{\mu }\frac{dp}{d{l}_1}$$

(47)

For surface diffusion in UPM model, the volume flux is:

$$q_{\mathrm{st}}=-A_{\mathrm{t}}\frac{k_{\mathrm{ts}}}{\mu }\frac{dp}{d{l}_{\mathrm{t}}}$$

(48)

According to the single capillary model (Figure 3b), the ratio of A_s to A₁ is:

$$\frac{A_{\mathrm{s}}}{A_1}=1-{(\frac{r_{\mathrm{eff}}}{r_1})}^2$$

(49)

Coupling Equations (41) and (42) and Equations (47)–(49), we have the correction factor from single capillary to UPM model for surface diffusion as:

$${\gamma }_{\mathrm{s}}=\frac{\phi }{\tau }(1-{(1-\frac{r_{\mathrm{ad}}}{r_1})}^2)$$

(50)

In Equations (45) and (50), the tortuosity τ can be estimated by the empirical formula [18,67,68]

$$\tau =1+0.63\ln (1/\varphi )$$

(51)

Finally, the AP considering PM properties and adsorbed gas occupation can be expressed as:

$$k_{\mathrm{app}}={\gamma }_{\mathrm{b}}F(Kn,p)\cdot \frac{r_{\mathrm{eff}}^2}{8}+{\gamma }_{\mathrm{b}}{\delta }^{D_{\mathrm{f}}-2}{D}_{\mathrm{k}}{c}_{\mathrm{g}}\mu +{\gamma }_{\mathrm{s}}{D}_{\mathrm{sur}}{G}_{\mathrm{L}}{\theta }_{\mathrm{real}}\frac{\mu Z}{p^2}\frac{T{p}_{\mathrm{sc}}}{T_{\mathrm{sc}}}$$

(52)

3. Model Validation

Currently, the use of Langmuir velocity slip BC for gas transport studies in SGRs is relatively rare and new in the petroleum industry. Thus, the accuracy of our LSP model, bulk phase gas transport model, single capillary model, and UPM model are all validated in this section.

First, the experimental data and lattice Boltzmann method (LBM) simulation results from Fathi et al. [69] are adopted in this section to validate our LSP model. The experimental data was collected from a routine gas uptake measurements using 5.68 mm size fragments. As can be seen from Figure 4, the matching results from three different sources are reasonable, which can prove the accuracy of our LSP model. Corresponding fitting parameters are presented in

Table 4. Input parameters for LSP model validation.

Parameter	Value	Parameter	Value
Overburden pressure, pc	90 MPa	Pore radius, r	10 nm
Pore pressure, p	0.69 to 40 MPa	Pressure at original state, p0	80 MPa
Pore roughness parameter, m	0.5	Gas constant, Rcst	8.314 J/(mol∙K)
Temperature, T	353 K	Biot constant, α	1
Pore stiffness parameter, p1	150 MPa	Reference pressure for LSP, pr	1 MPa
Adsorption parameter, ω	0.1

. Then, model comparison was done among first-order correlation [28], second-order correlation [30], corrected second-order correlation [33], and our LSP model (defining ω = 0.6, p_D = 1, σ_v = 0.6, and C = 0.55), as can be seen in Figure 5. During the slip flow regime, the first and second-order correlation results are much higher than corrected second-order correlations and our LSP results, whereas corrected second-order correlation results overlap with the LSP results. This is because the first and second-order correlations assumed elastic gas—solid collisions and uniform distribution of gas molecules across the cross-section area of nanopores.

The corrected second-order correlation method considered the gas accumulation effects near solid surfaces according to MD, whereas our LSP model considered gas adsorption on solid surfaces for slip flow study. Gas accumulation and adsorption both narrow flow paths and hinder gas transport in small pores. The first- and second-order correlation methods ignored the hindering effect and may overestimate the slip permeability in formation evaluation.

Then, to validate the bulk phase gas transport model, we compared the results from our model to those published data from linearized Boltzmann solutions [70] and experimental data [71]. Note that bulk phase gas is defined as that transported by slip viscous flow and Knudsen diffusion, and bulk phase permeability is defined as the sum of slip permeability k_slip and the equivalent Knudsen permeability k_nudsen. As can be seen from Figure 6, three different data sources adequately match with each other.

Further, the AP (considering bulk phase and adsorbed gas transporting) and bulk phase permeability (without considering surface diffusion) from our model (defining pore pressure p = 2 MPa, formation temperature T = 383 K, ω = 0.2, reference pressure p_r = 0.1 MPa, pore stiffness parameter m = 0.5, pore fractal dimension D_f = 2.6, and pore roughness δ = 0.25) is compared with those from the generalized lattice Boltzmann model in Wang et al. [72]. Clearly, both AP and bulk phase gas permeability match the Lattice Boltzmann method (LBM) results very well in different pore size conditions, as shown in Figure 7. Surface diffusion plays an important role in small pores and throats, but it can be ignored when pore radius is larger than 20 nm in this condition.

Finally, model comparison was done for AP models to further validate the proposed model’s correctness. The calculated results from different models are shown in Figure 8 (corresponding parameters used for comparison are presented in

Table 5. Parameters used for model comparison.

Parameter	Value	Parameter	Value
Gas type	CH4	Pore radius, r	30 nm
Pore pressure, p	10–70 MPa	Porosity, ϕ	0.060
Overburden pressure, pc	90 MPa	Reference pressure, pr	1 MPa
Adsorption parameter, ω	2	Temperature, T	414 K
Pore roughness parameter, m	0.5	Pore roughness, δ	0.7
Fractal dimension of pore surface, Df	2.5

). The model of Michel et al. [73] describes only bulk phase gas transport and ignores adsorbed gas. Thus, the results are quite lower than those of other models and basically constant with pressure changes. Because both adsorbed gas and bulk phase gas are considered, the results from Xiong et al. [64], Wasaki et al. [9], and our UPM model display a similar changing trend; that is, an increase of k_app/k_ins as pressure decreases. However, our model results are similar to Xiong et al. [64] at low pressures and tend to match those of Wasaki et al. [9] as pressure increases. This is because inelastic gas–solid collisions, thermodynamic phase changes due to pore confinement, the dynamic hindering effect of adsorbed gas, real gas effects, and geomechanical effects are simultaneously considered in our model. Besides, Knudsen diffusion and the slippage effect were discussed separately in our model, whereas in the other two models they are treated by a modification function of Knudsen number Kn. Furthermore, surface diffusion is influenced in real time by pressure and temperature in our model, whereas constant diffusion coefficients were chosen for whole production time in the other two models. All these factors led to the different decreasing trend of AP models in Figure 8.

4. Results and Discussion

In this section, the parameters sensitivity analysis will be done in different conditions, including the influences of pore radius, pressure, temperature, effective stress, and other parameters on the flux of each fluid flow mechanism; that is, viscous flow, Knudsen diffusion, and surface diffusion. Corresponding parameters used in our model are listed in

Table 6. Input parameters of gas and formation for model analysis.

Parameter	Value	Parameters	Value
Gas type	CH4	Pore radius, r	0.5–100 nm
Pore pressure, p	1–70 MPa	Porosity, ϕ	0.05
Langmuir pressure, pL	13.79 MPa	Gas molar weight, M	0.016 kg/mol
Temperature, T	373 K	Biot constant, α	1
Fractal dimension, Df	2.7	Reference pressure, pr	1 MPa
Surface roughness, δ	0.4	Gas constant, Rcst	8.314 J/(mol∙K)
Equivalent adsorption heat, ∆H	16,000 J/mol	vdW energy parameter, a	0.22998 m6 Pa/mol2
Pore stiffness parameter, p1	150 MPa	vdW energy parameter, b	4.28 × 10−5 m3/mol
Pore roughness parameter, m	4	Lennard–Jones size parameter, ζ	3.73 × 10−9 m
Adsorption parameter, ω	2	Lennard–Jones size parameter, ε	2.0434 × 10−21
Avogadro’s constant, N	6.02 × 1023	Langmuir volume, GL	20 m3/m3

.

With an increase in pore radius, slip flow permeability (Figure 9) as well as its contribution to total flux (Figure 10) also increases. In large pores (r > 20 nm) with pressures larger than 5 MPa, the effect of stress dependence is more obvious than that of the slippage effect. Thus, during the depressurization production, the slip permeability becomes less and less as effective stress increases. In small pores, the slippage effect dominates compared to stress dependence, and the slippage permeability increases as the pore pressure decreases. When the pore pressure is small enough, the permeability enhancement of the slippage effect is stronger than the permeability abatement of stress dependence in all pores, as shown in Figure 9.

Figure 10 shows that the contribution of slip flow increases with an increase in pore radius, and the increasing trend is faster in high-pressure and small pore radius conditions. In small pores, low pressure condition has a larger contribution than high pressure, while high pressure has a bigger contribution than small pressure in large pores. The slip flow contribution increases as pressure increases in pores with radii larger than 30 nm, although it displays complex changing rules in pores with radii smaller than 10 nm.

Knudsen diffusion becomes stronger with an increase in pore radius and a decrease in pressure (Figure 11). This is because the equivalent permeability of Knudsen diffusion is influenced mainly by the Knudsen diffusion coefficient in a certain pressure and temperature condition. For a constant pore radius, Knudsen diffusion increases with the decrease in pore pressure. Although gas viscosity increased as pressure increases, Knudsen diffusion permeability k_nudsen is more influenced by compressibility coefficient c_g, the two of which have a nearly inverse proportion relation. The contribution of Knudsen diffusion to total flux is very complex in small pores, as shown in Figure 12. Generally, that contribution decreases with pore radius in low-pressure conditions (p < 5 MPa), whereas it increases rapidly at first and slowly decreases later in high-pressure conditions (p > 10 MPa).

The change of surface diffusion with pore radius in high-pressure (p > 10 MPa) conditions is quite different from that in low-pressure (p < 5 MPa) conditions (Figure 13). It increases at first and decreases later as pore radius increases in low-pressure conditions, whereas it has a persistent decrease under high pressure. Actually, the detailed changing trend of surface diffusion intensity is not important, because its contribution is quite limited in large pore radii (r > 3 nm), no matter whether the value is low or high (Figure 14). However, its contribution to total flux is prominent in small pores. This is because it is advantageous for the occurrence of surface diffusion in small pores, when the porosity of PM is a constant. Meanwhile, the contribution of surface diffusion can be ignored in low-pressure conditions, although pore radii are small enough for the occurrence of surface diffusion. This is because adsorbed gas concentration is a function of the gas coverage ratio, which is very small in low-pressure conditions (Figure 15) and weakens the occurrence of surface diffusion.

Generally, non-Darcy flow mechanisms are dominant mainly in low-pressure and small pore radius conditions (Figure 16). With an increase in pore radii, the influence of the non-Darcy flow weakens and the flow tends to be viscous flow only.

For gas flow in pores with radii greater than 0.7 nm, the non-Darcy flow becomes more and more obvious during depressurization production. However, in contrast with previous studies, for gas flow in pores with radii less than 0.7 nm, with pressure decreasing, the non-Darcy flow effect is stronger at first, and then weakens abruptly. This is because surface diffusion dominates the flow in small pores when the pressure is larger than 10 MPa, but plays only a slight role when the pressure is less than 10 MPa (Figure 14).

As shown in Figure 17, the AP increases as pressure decreases in small pores, because the effect of the non-Darcy flow is stronger in low-pressure and small radius conditions. In large pores and throats, the AP is obviously influenced by both the non-Darcy flow effect and the geomechanical effect of stress dependence. The permeability decreases as pressure decreases at first due to stress sensitivity, and increases abruptly at low pressures, because the influence of the non-Darcy flow is more obvious than the geomechanical effects.

5. Conclusions

In this paper, a coupled AP model is established by introducing a Langmuir slip BC, simultaneously considering viscous flow, slippage effect, Knudsen diffusion, surface diffusion and adsorption/desorption. Both single capillary and capillary bundle models are established, which were validated by experimental data and LBM results. Coupled into the AP model were the effects of real gas, stress dependence, thermodynamic phase changes due to pore confinement, adsorbed layer thickness, pore enlargement due to gas desorption, surface roughness, and inelastic gas-solid collisions. Based on the AP model, the following conclusions can be made:

• The LSP model considers both gas slippage and gas adsorption effects on solid surfaces in nanopores. Gas slippage has a permeability enhancement effect, whereas gas adsorption effect is important in the estimation of permeability because adsorbed gas narrows flow paths.
• The non-Darcy flow is prominent in small pore size and low-pressure conditions and displays an abnormal increase with a decrease of pore radii in small pores due to the coexistence of surface diffusion and Knudsen diffusion.
• The geomechanical effect has a greater influence on AP of larger pores than the non-Darcy effect, whereas the AP of small pores is more influenced by the non-Darcy flow effect.
• The intensity of slip flow, Knudsen diffusion, and surface diffusion as well as their contribution to total flux was studied quantitatively; these intensities change with great complexity due to the combined effects of real gas, thermodynamics phase changes, and geomechanics.