Apparent Permeability Model for Gas Transport in Multiscale Shale Matrix Coupling Multiple Mechanisms

Abstract

Apparent gas permeability (AGP) is a significantly important parameter for productivity prediction and reservoir simulation. However, the influence of multiscale effect and irreducible water distribution on gas transport is neglected in most of the existing AGP models, which will overestimate gas transport capacity. Therefore, an AGP model coupling multiple mechanisms is established to investigate gas transport in multiscale shale matrix. First, AGP models of organic matrix (ORM) and inorganic matrix (IOM) have been developed respectively, and the AGP model for shale matrix is derived by coupling AGP models for two types of matrix. Multiple effects such as real gas effect, multiscale effect, porous deformation, irreducible water saturation and gas ab-/de-sorption are considered in the proposed model. Second, sensitive analysis indicates that pore size, pressure, porous deformation and irreducible water have significant impact on AGP. Finally, effective pore size distribution (PSD) and AGP under different water saturation of Balic shale sample are obtained based on proposed AGP model. Under comprehensive impact of multiple mechanisms, AGP of shale matrix exhibits shape of approximate “V” as pressure decrease. The presence of irreducible water leads to decrease of AGP. At low water saturation, irreducible water occupies small inorganic pores preferentially, and AGP decreases with small amplitude. The proposed model considers the impact of multiple mechanisms comprehensively, which is more suitable to the actual shale reservoir.

1. Introduction

With the gradual depletion of conventional energy, unconventional resource such as shale gas, methane hydrate and geothermal energy has drawn extensive attention [1,2]. In 2019, American shale gas production estimated by EIA was about 25.28 trillion cubic feet (715.85 billion cubic meters), which amounts to 17.42% of total world natural gas production [3]. With the huge reserve and great potential exploitation, shale gas will become a significant part of the energy supply [4,5]. Production of shale gas wells quickly enters a long-term stage of low production after the initial stage of rapid production decline. The ultimate production is determined by the long-term stage of low production, which indicates the importance of gas flow behavior in shale matrix [6]. Therefore, accurate evaluation of transport capacity for shale matrix is significant for improving the recovery and economic benefits of shale gas reservoirs.

Shale matrix can be classified as organic matrix (ORM) and inorganic matrix (IOM) based on component, and the pore sizes in shale matrix are about 1–1000 nm [4,7]. Gas storage and transport characteristics in two categories of matrix are distinct, due to the different component, surface wettability and pore structures [8]. Core scanning images shows that the majority of pores in ORM are circular, and pores in IOM mainly exhibit slit shape [9]. Gas mainly occurs in ORM as free and adsorbed states owing to hydrophobicity and huge specific surface area. Bulk gas flow and surface diffusion coexist in ORM. IOM is strongly hydrophilic due to its high component of clay and quartz, and irreducible water mainly occurs in IOM [10]. Free gas occupies the center of inorganic pore, and irreducible water is absorbed in the pore wall as water film. Bulk gas flow is dominated transport pattern in IOM. As irreducible water saturation ranges from 0.15 to 0.6 in shale reservoir, influence of irreducible water on gas transport behavior in nanopores cannot be neglected. Furthermore, real gas effect, gas ab-/de-sorption effect and porous deformation also influence gas transport behavior in shale matrix [11]. Thus, accurate evaluation of transport capacity for shale matrix must consider coupling impact of multiple mechanisms [12].

Capacity of gas transport throw shale matrix is usually represented by apparent gas permeability (AGP), which can be determined by laboratory experiment, numerical simulation and analytical model. With extremely low-permeability, experimental methods make it difficult to obtain the AGP of shale sample accurately [13]. Simulation methods can accurately investigate the fluid transport behavior and storage at nanoscale [14]. However, simulation methods are not suitable for reservoir field due to its expensive and time-consuming process in nowadays [15]. Analytical model under some reasonable assumption can describe gas transport behavior under the coupling influence of some important mechanisms, which is suitable for multiscale medium [16]. Analytical models can be classified into two types, one of which is modifying nonslip boundary to describe gas transport characteristics in nanopore, such as 1st order, 1.5 order, 2nd order and Langmuir slip condition [17,18]. Beskok and Karniadakis [19] investigated gas transport behavior in micro channels and established a second order slip boundary model (B&K model). Considering phase behavior change and multiple gas transport mechanisms, Song et al. [20] proposed AGP models for ORM and IOM. Considering the influence of stress effect, real gas and of water film on gas flow, Sun et al. [21] improved B&K model to describe gas transport in organic and inorganic pore. With multiple mechanisms coupling, Li et al. [22] developed an AGP model for shale matrix based on the Langmuir slip condition. However, multiscale effect of shale matrix is neglected in Li’s model. The accuracy of slip boundary model depends on empirical coefficients, so more experimental or molecular simulation data is needed due to complex composition of shale matrix [23]. Another type of analytical model is using weighting coefficient to couple multiple transport mechanisms. Darabi et al. [13] and Javadpour et al. [24] coupled multiple mechanisms through liner superposition method, but this method overestimates the AGP of shale matrix in transform regime. Using collision frequency as weighting coefficient, Wu et al. [12,16,25] proposed a series of AGP models to describe gas transport within different cross section shapes’ nanopores. Wang et al. [26] improved Wu’s model to investigate impact of water film and PSD on gas flow in shale matrix. In general, gas transport behavior is still controversial and unclear due to rarefaction effect at nanoscale. At present, multiple mechanisms such as porous deformation, gas ad/desorption, complicated flow patterns and real gas effect are considered in published AGP models. However, effective hydraulic aperture is usually used in most models and few models consider pore size distribution. Meanwhile, irreducible water distribution in shale matrix is not addressed in previous models. Gas transport behavior in nanopores is strongly related to aperture. Ignoring the impact of multiscale effect and irreducible water distribution will lead to overestimating AGP of shale matrix.

In this work, circle nanopores in ORM and slit nanopores in IOM are combined to characterize complicated pore structure in shale matrix, and AGP model for multiscale shale matrix is established. The reliability of the proposed model is validated by result from experiment and published models. Further, the impacts of physical effects on AGP are serious discussed, such as pressure, aperture, pore geometry, porous deformation and irreducible water saturation. The effective PSD, AGP and relative permeability of Balic shale sample at different water saturation are obtained. Distinguishing features of AGP and relative permeability calculated by effective hydraulic aperture and PSD are discussed at last. The proposed model reveals the impact of irreducible water saturation on the gas transport characteristics in multiscale shale matrix and provides reference for engineering application.

2. Theoretical Model

Shale matrix is composed of ORM and IOM, and gas storage and transport mechanisms are distinctly different in two types of matrixes (Figure 1). Expression of effective pore aperture under comprehensive impact of multiple mechanisms is derived firstly in this section. Then, AGP models of ORM and IOM are derived respectively. Finally, the AGP model of actual shale matrix coupling multiple mechanisms is established. The basic assumptions for AGP models are as follows:

(1)

Nanopores in ORM and IOM are represented by circle tubes and slit tubes separately.

(2)

The aperture of nanopores in shale matrix is not uniform.

(3)

Gas occurs in organic pores as free and absorbed state, reservoir water can be neglected.

(4)

Free gas occupies the center of inorganic pore, and irreducible water is absorbed on the pore wall.

(5)

Shale matrix is compressible.

2.1. Real Gas Effect

During reservoir depletion, the shale gas properties change with change of pressure. Gas properties including gas deviation factor Z and viscosity are the important parameters, which affect gas transport behavior at nanoscale. Equation of State (EOS) can predict Z accurately, but this method involves iterative computation and many parameters [27,28]. Explicit empirical corrections can directly calculate Z, and which is reasonable under some conditions [29]. Ehsan et al. [30] proposed an explicit method by multiple regression analysis, and the expression is shown as follows:

$$Z=\frac{A_{\mathrm{s}1}+A_{\mathrm{s}2}\ln p_{\mathrm{pr}}+A_{\mathrm{s}3}{\left(\ln p_{\mathrm{pr}}\right)}^2+A_{\mathrm{s}4}{\left(\ln p_{\mathrm{pr}}\right)}^3+\frac{A_{\mathrm{s}5}}{T_{\mathrm{pr}}}+\frac{A_{\mathrm{s}6}}{T_{\mathrm{pr}}^2}}{1+A_{\mathrm{s}7}\ln p_{\mathrm{pr}}+A_{\mathrm{s}8}{\left(\ln p_{\mathrm{pr}}\right)}^2+\frac{A_{\mathrm{s}9}}{T_{\mathrm{pr}}}+\frac{A_{\mathrm{s}10}}{T_{\mathrm{pr}}^2}}$$

(1)

$$p_{\mathrm{pr}}=p/p_{\mathrm{pc}},T_{\mathrm{pr}}=T/T_{\mathrm{pc}}$$

(2)

where A_S1-A_s10 are coefficients; p is the pore pressure, Pa; p_pc is the critical pressure of gas, Pa; T is the temperature, K; T_pc is the critical temperature, K.

Equation (1) is applicable for predicting Z in the range of 0.2 < p_pr < 15 and 1.6 < T_pr < 3 [17]. Shale gas (mainly composition with methane) reservoir pressure and temperature condition is in the range of 5 MPa < p < 60 MPa and 310 K < T < 400 K, i.e., 1.02 < p_pr < 13.04 and 1.62 < T_pr < 2.09. Therefore, Equation (1) can be used for prediction Z in shale reservoir.

Expressions of viscosity and compressibility for shale gas in reservoir are given as follows [31]:

$$\mu ={\mu }_{1\mathrm{atm}}\left[1+\frac{7.9}{T_{\mathrm{r}}^5}\left(\frac{p_{\mathrm{r}}^4}{T_{\mathrm{r}}^{20}+p_{\mathrm{r}}^4}\right)+9\times {10}^{-6}{\left(\frac{p_{\mathrm{r}}}{T_{\mathrm{r}}}\right)}^2+0.28\frac{p_{\mathrm{r}}}{T_{\mathrm{r}}}\right]$$

(3)

$$c_{\mathrm{g}}=\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial p}\right)}_{\mathrm{T}}=\frac{1}{p}-\frac{1}{Z\left(p\right)}{\left(\frac{dZ\left(p\right)}{dp}\right)}_{\mathrm{T}}$$

(4)

where c_g is the compressibility factor, 1/Pa; μ is gas viscosity, Pa·s; μ_1atm is ideal gas viscosity when pressure equals 1 atm, Pa·s.

The molecular mean free path λ is defined as [32]:

$$\lambda =\frac{{\mu }_{\mathrm{g}}}{p}\sqrt{\frac{\pi ZRT}{2M}}$$

(5)

where M is gas molecular weight, kg/mol; λ is the molecular mean free path, m; R is ideal gas constant, 8.314 J/(mol·K).

2.2. Flow Channel of Divers Matrix

During shale reservoir depletion, aperture of nanopores in ORM and IOM is affected by porous deformation, gas ab/de-sorption effect and absorbed water film. The expression of effective pore size under the influence of multiple mechanisms is derived in this section.

2.2.1. Flow Channel for ORM

The complex pore structure in shale matrix can be represented by capillary bundle model. Based on Kozeny–Carman equation, permeability and porosity of ORM can be expressed by hydraulic radius; the expressions are shown as below [33,34]:

$$k_{\mathrm{t}0}=\frac{{\eta }_{\mathrm{org}}{\varphi }_0\cdot r_0^2}{8{\tau }^2}=\frac{{\varphi }_{\mathrm{t}0}\cdot r_0^2}{8{\tau }^2}$$

(6)

$${\varphi }_{\mathrm{t}0}=n\pi \tau r_0^2$$

(7)

where n is the number of organic pores per unit area, 1/m²; k_t0 is the intrinsic permeability of ORM, m²; ϕ_t0 is the porosity of ORM at the initial condition, dimensionless; r₀ is the hydraulic radius of ORM at initial conditions, m; τ is the tortuosity, dimensionless; η_org is volume content of organic pores, dimensionless; ϕ₀ is the porosity of shale matrix.

During the reservoir depletion, porosity, radius and permeability of shale matrix will decrease by increase of effective pressure. Permeability has exponent relation to effective stress; permeability of ORM considering porous deformation can be expressed as below [35,36]:

$$k_{\mathrm{stress}}=k_0{e}^{-{\psi }_{\mathrm{t}}\left(p_{\mathrm{c}}-p\right)}$$

(8)

where ψ_t is the compressibility coefficient, pa⁻¹; k_stress is the permeability of ORM considering porous deformation, m²; p_c is the overlying formation pressure, Pa.

Substituting Equations (6) and (7) into (8), and ratio of permeability at different effective tress can be expressed as follows:

$$\frac{{\varphi }_{\mathrm{stress}}{r}_{\mathrm{stress}}^2}{{\varphi }_0{r}_0^2}=\frac{r_{\mathrm{stress}}^4}{r_0^4}=e^{-{\psi }_{\mathrm{t}}\left(p_{\mathrm{c}}-p\right)}$$

(9)

Rearranging Equation (9), the hydraulic radius of organic pores can be obtained as:

$$r_{\mathrm{stress}}=r_0{e}^{-0.25{\psi }_{\mathrm{t}}\left(p_{\mathrm{c}}-p\right)}$$

(10)

where r_stress is the hydraulic radius of organic pore considering porous deformation, m; ϕ_stress is the porosity of ORM considering porous deformation, dimensionless.

The presence of absorbed gas will decrease the effective radius of organic pores (Figure 2). Influence of absorbed gas layer on hydraulic radius cannot be ignored, due to gas molecular diameter becoming comparable to organic pore size [37]. Gas coverage is commonly used to qualify thickness of absorbed gas layer for monolayer adsorption, and which is defined as follows [17]:

$$\theta =\frac{p/Z}{p_{\mathrm{L}}+p/Z}$$

(11)

$$r_{\mathrm{ad}}=d_{\mathrm{m}}\cdot \theta$$

(12)

where θ is the Gas coverage in surface of organic pores, dimensionless; p_L is the Langmuir pressure, Pa; d_m is the molecular diameter of gas, 0.38 nm; r_ad is the thickness of absorbed gas layer, m.

Effective radius of organic pore accounting for coupling influence of porous deformation and adsorption can be expressed as follows:

$$r_{\mathrm{eff}}=r_{\mathrm{stress}}-r_{\mathrm{ad}}$$

(13)

Porosity of ORM is directly proportional to the square of equivalent radius

$${\varphi }_{\mathrm{eff}-\mathrm{or}}={\varphi }_{\mathrm{t}0}\cdot {\left(\frac{r_{\mathrm{ef}}{}_{\mathrm{f}}}{r_0}\right)}^2={\varphi }_{\mathrm{t}0}\cdot e^{-0.5{\psi }_{\mathrm{t}}\left(p_{\mathrm{c}}-p\right)}$$

(14)

where r_eff is effective radius, m; ϕ_eff-or is the effective porosity of ORM, dimensionless.

2.2.2. Flow Channel for IOM

The pores in IOM mostly present slit type, rectangular pore with the same height and cross-sectional area is used to represent inorganic pores in this work (Figure 3).

$$S=W_{0}H=\xi H_0^2$$

(15)

where S is the cross-sectional area, m²; H₀ is the width of inorganic pores, m; W₀ is the height of inorganic pore, m; ξ is the aspect ratio, dimensionless.

The aspect ratio is defined as follows:

$$\xi =\frac{W_0}{H_0}$$

(16)

Reservoir water molecules are absorbed on the surface of inorganic pore as water film and they are immobile under low-pressure gradient [38]. Irreducible water and porous deformation have great impact on effective hydraulic aperture and porosity of IOM. The thickness of absorbed water film in slit nanopores of IOM can be quantified by Li’s model [39].

The initial irreducible water saturation of IOM can be expressed as:

$$S_{\mathrm{w}}\approx \frac{2h}{H_0}$$

(17)

where S_w is the irreducible water saturation of nanopore, dimensionless; h is the thickness of water film, m.

Similar to organic pores, effective width of inorganic pore and effective porosity of IOM accounting for coupling absorbed water film and porous deformation can be obtained by:

$$H_{\mathrm{eff}}=H_0{e}^{-0.25}{}^{{\psi }_{\mathrm{t}}\left(p_{\mathrm{c}}-p\right)}-2h$$

(18)

$${\varphi }_{\mathrm{eff}-\mathrm{IOM}}=\left(1-{\eta }_{\mathrm{org}}\right){\varphi }_{\mathrm{s}0}{\left(e^{-{0.25}^{-{\psi }_{\mathrm{t}}\left(p_{\mathrm{c}}-p\right)}}-S_{\mathrm{w}}\right)}^2$$

(19)

where ϕ_s0 is the porosity of IOM at the initial condition, dimensionless; H_eff is effective width of inorganic pore, m; ϕ_eff-IOM is the effective porosity of IOM, dimensionless.

2.3. Gas Transport in ORM

Gas mainly occurs in ORM as absorbed and free state. During shale gas depletion, absorbed gas not only desorbs into free gas but transport in the form of surface diffusion. Mass flux of gas in ORM is the sum of bulk gas flow and surface diffusion.

2.3.1. Surface Diffusion

Driven by concentration gradient, absorbed gas could be transported in organic pores via surface diffusion. The jumped model is applied to describe surface diffusion; mass flux of surface diffusion in organic pore can be expressed as [20]:

$$J_{\mathrm{s}}=-M{A}_{\mathrm{s}}{D}_{\mathrm{s}}\nabla \left(C_{\mathrm{smax}}\theta \right)$$

(20)

where J_s is mass flux of surface diffusion in organic pore Kg/s; A_s is annulus area of absorbed gas layer, m²; D_s is the surface diffusion coefficient, m²/s; C_smax is the maximum absorbed gas concentration of shale sample, mol/m³.

Expression of annulus area of absorbed gas layer and the maximum concentration of absorbed gas in shale sample can be written:

$$A_{\mathrm{s}}=\pi r_{\mathrm{ad}}\left(2r_{\mathrm{eff}}-r_{\mathrm{ad}}\right)$$

(21)

$$C_{\mathrm{smax}}=\frac{{\rho }_{\mathrm{s}}{V}_{\mathrm{L}}}{V_{\mathrm{std}}}$$

(22)

where ρ_s is density of shale, kg/m³; V_L is the Langmuir volume, m³/m³; V_std is molar volume of gas, m³/mol.

According to the equation, surface diffusion flux is given in the form of pressure gradient:

$$J_{\mathrm{s}}=-M{A}_{\mathrm{s}}{D}_{\mathrm{s}}{C}_{\mathrm{smax}}\frac{d\theta }{dp}\nabla p=-M{A}_{\mathrm{s}}{D}_{\mathrm{s}}{C}_{\mathrm{smax}}\frac{Z{p}_{\mathrm{L}}}{{\left(Z{p}_{\mathrm{L}}+p\right)}^2}\nabla p$$

(23)

where Δp is the pressure gradient of along the capillary, Pa.

The surface diffusion coefficient D_s is defined as [40]:

$$D_{\mathrm{s}}=D_{\mathrm{s}}^0\frac{\left(1-\theta \right)+\frac{\kappa }{2}\theta \left(2-\theta \right)+\left[H\left(1-\kappa \right)\right]\left(1-\kappa \right)\frac{\kappa }{2}{\theta }^2}{{\left(1-\theta +\frac{\kappa }{2}\theta \right)}^2}$$

(24)

$$H\left(1-\kappa \right)=\{\begin{array}{l}0,\kappa \ge 1\\ 1,0\le \kappa \le 1\end{array}$$

(25)

$$D_{\mathrm{s}}^0=-8.29\times {10}^{-7}\times T^{0.5}\exp \left[-\nabla H^{0.8}/\left(RT\right)\right]$$

(26)

where ▽H is the equivalent adsorption heat, J/mol; D_s0 is surface diffusion coefficient at zero gas coverage, m²/s; κ is the ratio of rate constant for blockage to the rate constant for forward migration, dimensionless.

Equation (26) is obtained from fitting experimental data of methane/active carbon at low pressure condition, and D_s is obtained by modifying D_s0 with gas coverage to be applicable to high pressure condition. The main component of shale gas is methane; therefore, D_s is suitable for temperature and pressure condition of shale gas reservoir.

According to Equation (23), the AGP of surface diffusion is obtained:

$$k_{\mathrm{s}}=\frac{J_{\mathrm{s}}\mu }{{\rho }_{\mathrm{b}}\pi r_{\mathrm{eff}}^2}=\mu M{D}_{\mathrm{s}}{C}_{\mathrm{smax}}\frac{Z{p}_{\mathrm{L}}}{{\rho }_{\mathrm{b}}{\left(Z{p}_{\mathrm{L}}+p\right)}^2}\left(1-\frac{r_{\mathrm{ad}}^2}{r_{\mathrm{eff}}^2}\right)$$

(27)

where ρ_b is the density of free gas, Kg/m³.

2.3.2. Bulk Gas Flow

Knudsen number (Kn) of gas within organic pores is written as:

$$K\mathrm{n}=\frac{\lambda }{r_{\mathrm{eff}}}$$

(28)

When Kn is in the range of 0.001~0.1, the dominated transport mechanism in nanopores is slip flow. According to B&K model, the AGP of slip flow in organic pores can be expressed by [20]:

$$k_{\mathrm{slip}-\mathrm{tube}}=f\left(K\mathrm{n}\right)k_{\infty }=\left(1+\frac{128}{15{\pi }^2}{\tan }^{-1}\left[4.0K{\mathrm{n}}_{}^{0.4}\right]K\mathrm{n}\right)\left(1+\frac{4K\mathrm{n}}{1+K\mathrm{n}}\right)\frac{r_{\mathrm{eff}}^2}{8}$$

(29)

where k_slip-tube is AGP of slip flow within organic pore, m². When Kn > 10, Knudsen diffusion is the dominated transport pattern in capillaries. With the solid surface roughness increase, the residence time for molecules staying around surface increases after collisions with the walls. Surface roughness will decrease diffusion efficiency. Considering the impact of solid surface roughness, the AGP of Knudsen diffusion is given by [13]:

$$k_{\mathrm{Kn}-\mathrm{tube}}={\delta }^{D_{\mathrm{f}}-2}\frac{2\mu C_{\mathrm{g}}{r}_{\mathrm{eff}}}{3}{\left(\frac{8Z\mathrm{R}T}{\pi M}\right)}^{0.5}$$

(30)

where δ is the surface roughness of pore, dimensionless; k_Kn-tube is AGP of Knudsen diffusion within organic pore, m²; D_f is the fractal dimension of shale sample, dimensionless. The Kn of shale gas under reservoir condition ranges from 0.0002 to 6, so the transport mechanisms of bulk gas flow in shale matrix encompassing slip flow and Knudsen diffusion [16]. Based on molecular collision frequency, the expressions of weighting coefficient of diverse mechanisms in organic pore are given as [26]:

$$f_{\mathrm{slip}-\mathrm{tube}}=\frac{1}{1+K\mathrm{n}},f_{\mathrm{Kn}-\mathrm{tube}}=\frac{K\mathrm{n}}{1+K\mathrm{n}}$$

(31)

where f_slip-tube is weighting coefficient of slip flow in organic pore, dimensionless; f_Kn-tube is weighting coefficient of Knudsen diffusion in organic pore, dimensionless.

Here, weighted superposition method is applied to couple slip flow and Knudsen diffusion to describe complicated gas transport behavior in organic pore reasonably. Therefore, AGP of bulk gas flow in a single organic pore is expressed as:

$$k_{\mathrm{bt}}=f_{\mathrm{slip}-\mathrm{tube}}\cdot k_{\mathrm{slip}-\mathrm{tube}}+f_{\mathrm{Kn}-\mathrm{tube}}\cdot k_{\mathrm{Kn}-\mathrm{tube}}$$

(32)

where k_bt is AGP of bulk gas within organic pore, m².

2.3.3. AGP Model of ORM

Gas transports pattern in organic pores encompassing bulk gas flow and surface diffusion, the total flux is the sum of the diverse mechanisms [25]. Assuming that ORM consists of capillaries with the same radius, combine Equations (27) and (32), gas flow flux in ORM can be calculated by:

$$Q=n\pi r_{\mathrm{eff}}^2\frac{\left(k_{\mathrm{s}}+k_{\mathrm{bt}}\right)\Delta p}{\mu \tau L}$$

(33)

where Q is the gas flux in ORM, m³/s; L is the characteristic length of ORM, m.

The flux in ORM can be also expressed as the form of Darcy’s law:

$$Q=\frac{k_{\mathrm{app}-\mathrm{tube}}\Delta p}{\mu L}$$

(34)

where k_app-tube is AGP for ORM, m². Substituting Equations (7) and (33) into (34), the AGP of ORM coupling multiple transport mechanisms and physical effects can be derived as:

$$k_{\mathrm{app}-\mathrm{tube}}=n\pi r_{\mathrm{eff}}^2\frac{L}{\tau L}\left(k_{\mathrm{s}}+k_{\mathrm{bt}}\right)=\frac{{\varphi }_{\mathrm{eff}-\mathrm{ORG}}}{{\tau }^2}\left(k_{\mathrm{s}}+k_{\mathrm{bt}}\right)$$

(35)

Equation (35) encompasses multiple gas transport mechanisms in organic matrix, which can describe gas transport behavior for all flow regime.

2.4. Gas Transport in IOM

For slit-like pores, the characteristic length is the width. Considering the influence of coupling porous deformation and water film, Kn of real gas within inorganic pore is defined as:

$$K\mathrm{n}=\frac{\lambda }{H_{\mathrm{eff}}}$$

(36)

Water film is immobile due to greater capillary force and adsorption force in inorganic pores. Only bulk gas flow exists in inorganic pore. Similar to organic pore, weighted superposition method is also applied to capture complex bulk gas transport behavior in inorganic pore. Considering correction of cross section shape, the AGP of slip flow in inorganic pore can be written as follows [41]:

$$k_{\mathrm{slip}-\mathrm{slit}}=A\left(\xi \right)\cdot \left(1+\frac{128}{15{\pi }^2}{\tan }^{-1}\left[4.0K{\mathrm{n}}_{}^{0.4}\right]K\mathrm{n}\right)\left(1+\frac{6K\mathrm{n}}{1+K\mathrm{n}}\right)\frac{H_{\mathrm{eff}}^2}{12}$$

(37)

where k_slip-slit is AGP of slip flow within inorganic pore, m²; A(ξ) is the section shape factor for slip flow, dimensionless. Expression of A(ξ) is shown as follows:

$$A\left(\xi \right)=1-\frac{192}{\xi {\pi }^2}{\displaystyle \sum _{\mathrm{i}=1,3,5\dots }^{\infty }\frac{\tanh \left(i\pi \xi /2\right)}{i^5}}$$

(38)

Considering impact of surface roughness, the AGP of Knudsen diffusion for real gas transports in inorganic pores is given by [16]:

$$k_{\mathrm{Kn}-\mathrm{slit}}={\delta }^{D_{\mathrm{f}}-2}\mu C_{\mathrm{g}}B\left(\xi \right)\cdot H_{\mathrm{eff}}{\left(\frac{\mathrm{ZR}T}{2\pi M}\right)}^{0.5}$$

(39)

where k_kn-slit is AGP of Knudsen diffusion within inorganic pore, m²; B(ξ) is the section shape factor for Knudsen diffusion dimensionless; Its expression is shown as follows:

$$B\left(\xi \right)=\frac{1}{\xi }\cdot \left[{\xi }^2\ln \left(\frac{1}{\xi }+\sqrt{1+\frac{1}{{\xi }^2}}\right)+\xi \ln \left(\xi +\sqrt{1+{\xi }^2}\right)-\frac{{\left(1+{\xi }^2\right)}^{3/2}}{3}+\frac{\left(1+{\xi }^3\right)}{3}\right]$$

(40)

Based on molecular collision frequency, the expressions of weighting coefficient for diverse transport mechanisms in inorganic pore are given as:

$$f_{\mathrm{slip}-\mathrm{slit}}={\left[1+\frac{K\mathrm{n}\left(1+1/\xi \right)}{2}\right]}^{-1},f_{\mathrm{Kn}-\mathrm{slit}}={\left[1+\frac{2}{K\mathrm{n}\left(1+1/\xi \right)}\right]}^{-1}$$

(41)

where f_slip-slit is weighting coefficient of slip flow for inorganic pore, dimensionless; f_Kn-slit is weighting coefficient of Knudsen diffusion for inorganic pore, dimensionless.

Thus, the AGP of IOM coupling multiple transport mechanism and physical effects is obtained as follows:

$$k_{\mathrm{app}-\mathrm{slit}}=\frac{{\varphi }_{\mathrm{eff}-\mathrm{IOM}}}{{\tau }^2}\left(f_{\mathrm{slip}-\mathrm{slit}}\cdot k_{\mathrm{slip}-\mathrm{slit}}+f_{\mathrm{Kn}-\mathrm{slit}}\cdot k_{\mathrm{Kn}-\mathrm{slit}}\right)$$

(42)

where k_app-slit is AGP for IOM, m².

2.5. AGP Model for Shale Matrix

Actually, shale matrix consists of a series of capillary with different pore size. Combining Equations (35) and (42) and pore size distribution, the AGP of homogeneous porous medium can be extended to real porous medium. The AGP of real ORM and IOM can be respectively expressed as:

$$k_{\mathrm{app}-\mathrm{org}}={\displaystyle \underset{r_{\min }}{\overset{r_{\max }}{\int }}\left[k_{\mathrm{app}-\mathrm{tube}}\left(r\right)\cdot f\left(r\right)\right]}\cdot d{r}_{}$$

(43)

$$k_{\mathrm{app}-\mathrm{iom}}={\displaystyle \underset{H_{\min }}{\overset{H_{\max }}{\int }}\left[k_{\mathrm{app}-\mathrm{slit}}\left(H\right)\cdot f\left(H\right)\right]}\cdot d{H}_{}$$

(44)

where k_app-org is the AGP for ORM considering multiscale, m²; k_app-iom is the AGP for IOM considering multiscale, m²; f (r) and f (H) is the volume fraction of organic pore and inorganic pore in diverse matrix, respectively, dimensionless.

Based on the capillary bundle model, organic and inorganic pores are parallel and without mass exchange. Therefore, the total flux in shale matrix is the sum of flux in two types of matrix, AGP of shale matrix can be written:

$$k_{\mathrm{app}}=k_{\mathrm{app}-\mathrm{ior}}+k_{\mathrm{app}-\mathrm{org}}$$

(45)

where k_app is the AGP of shale matrix, m².

The initial irreducible water saturation of shale matrix can be expressed as follows:

$$S_{\mathrm{w}-\mathrm{shale}}=\left(1-{\eta }_{\mathrm{org}}\right){\displaystyle \underset{H_{\min }}{\overset{H_{\max }}{\int }}{S}_{\mathrm{w}}\left(H\right)\cdot F\left(H\right)}dH$$

(46)

where S_w-shale is the irreducible water saturation of shale matrix, dimensionless.

The proposed AGP model of shale matrix coupling impacts multiscale, surface roughness, pore geometry, porous deformation, ad/desorption, real gas effect and water distribution and volume content of organic pores on multiple gas transport. Compared with the published AGP models, the proposed AGP model encompasses comprehensive impact of internal and external factors, especially the impact of multiscale and irreducible water saturation, which is a more suitable shale gas reservoir.

3. Model Validation

According to the above derivation, gas transport behavior in shale matrix is influenced by multiple mechanisms such as porous deformation, irreducible water saturation. Experimental results considering coupling impact of multiple mechanisms are sparse. Thus, molecular simulation and experimental data from literature and AGP calculated from published models is applied to validate the proposed AGP model of ORM and IOM in this paper, respectively. Thus, the accuracy of the proposed AGP model for actual shale matrix in this paper is confirmed indirectly. Parameters for model validation are exhibited in

Table 1. Basic parameters for model validation.

Parameters	Value
Temperature, T	300 K
Overlying formation pressure, pc	50 MPa
Critical pressure, ppc	4.599 MPa
Critical temperature, Tpc	190.56 K
Viscosity of ideal gas, μ1atm	1.49 × 10−5 Pa·s
Methane molar mass, M	0.016 kg/mol
Langmuir pressure, pL	13.7 MPa [22]

.

3.1. Validation of AGP Model for ORM

Three published AGP models [21,22,24] are compared with the proposed model to further demonstrate accuracy of the proposed model considering porous deformation. AGP of ORM calculated by four models versus pressure is plotted in Figure 4. The variation of gas density with pressure was neglected in Li’s model, so the AGP is smaller than results calculated by other models. AGP from Javadpour’s model is larger than other three models, especially at low pressure condition, this is because multiple transport mechanisms are reconsidered in his model. When the pressure is larger 10 MPa, the AGP from the proposed model are approximate to result from Sun’s model. Sun’s model is based on Civan’s model, which will overestimate AGP of nanopore in transition flow regime. Our model is more reasonable than Sun’s model at low reservoir pressure condition. Thus, the proposed AGP model is suitable to describe gas transport in ORM.

3.2. Validation of AGP Model for IOM

Experimental data from published literature is compared with the result calculated from the proposed AGP model of IOM (Equation (42)). Wu et al. [42] conducted gas displaying water experiment in a nanofluidic chips approach, which is consisted of 100 slit nano-channels with dimension of 100 nm (depth) × 5 μm(width). Rushing et al. [43] measured AGP at different pressure and water saturation in tight sandstone core, and the hydraulic pore width of sample equals 15 nm. Phenomena observed in their experiments is (gas core flows in the center of channels, and immobile water film attaches to pore wall) similar to assumption of the proposed model. From Figure 5, AGP results calculated from proposed AGP model of IOM are close to the experimental data from homogeneous porous medium and natural sample. The proposed model coupling multiple physical effects and multiple mechanisms is suitable for IOM.

4. Discussion

In this section, the influence of multiple effects, including water saturation, real gas effect, pressure, pore geometry, pore size, porous deformation on the AGP in ORM, IOM and shale matrix is analyzed based on newly developed models. It is worth noting that the sensitive analysis and discussion for ORM and IOM is under the assumption of uniform aperture. Numerical value of basic parameters used for analysis is exhibited in

Table 2. Basic parameters of fluid and formation for sensitivity analysis.

Parameters	Value
Pore radius of ORM, r0	2–50 nm [16]
Pore width of IOM, H0	2–50 nm [16]
Porosity of ORM, Φor0	0.018 [44]
Porosity of IOM, Φior0	0.05 [44]
Temperature, T	353 K
Pore pressure, p	0.5~50 MPa
Surface diffusion coefficient at zero gas coverage, DS0	5.69 × 10−6 m2/s [16]
Equivalent adsorption heat, ▽H	1.6 × 104 J
The ratio of rate constant for blockage to the rate Constant for forward migration, κ	0.5
Tortuosity, τ	2.2 [21]
Maximum concentration absorbed gas in shale core sample, Cmax	328.7 mol/m3
Compressibility coefficient of ORM and IOM, Ct	0.01 MPa−1 [36]
Fractal dimension of core sample, Df	2.6 [37]
Water saturation, Sw-shale	0–0.5

.

4.1. Gas Transport in ORM

Figure 6 presents the variation of AGP for multiple transport mechanisms in ORM. During reservoir depletion, surface diffusion increases monotonously and Knudsen diffusion and slip flow initially decrease and then increase. Nanoscale effect promotes gas flow, due to Kn increases with decrease of pressure. Porous deformation has disadvantage on effective pore radius, meanwhile gas desorption has advantage on effective radius. Thus, influence of porous deformation on AGP dominated in the high pressure, coupling influence of nanoscale effect and gas desorption prevail in low pressure.

From Figure 7, the gas conductance of ORM consists of triple rival transport mechanisms, and the contribution of multiple mechanisms are different in ORM with various radius. Surface diffusion is the dominated transport pattern in ORM consists of small pores (r < 2 nm). With the increase of radius in ORM, contribution of slip flow increases. Slip flow is the dominated flow pattern in ORM, when organic pores radius is larger than 15 nm. Generally, Knudsen diffusion has little contribution to AGP.

From Figure 8, AGP of ORM depends on pressure and radius, and the AGP of ORM with different radius exhibits different variation characteristics. Due to surface diffusion is the dominated transport mechanism (Figure 7a), AGP of ORM consists of small pores (r < 5 nm) increases monotonously. AGP of ORM exhibits shape of approximate “V” for radius larger than 5 nm. Real gas effect enhances AGP of ORM consisting of significant small pores. Deviation of AGP between two types of gas is less than 5% when radius in ORM equals 5 nm. AGP in actual shale matrix mainly contributes by pore with large radius, thus real gas effect can be ignored.

In order to simplify the AGP model, analysis of AGP affected by different effects coupling is conducted (Figure 9). Compared with total AGP, the deviation of AGP without considering surface diffusion and porous deformation are 25.93% and 27.78%, respectively. Thus, porous deformation and surface diffusion are the significant effects on the AGP.

4.2. Gas Transport in IOM

In this part, we focus on the coupling impact of irreducible water distribution, porous deformation, pore geometry and nanoscale effect on AGP. The aspect ratio of inorganic pores is assumed as four.

From Figure 10, AGP highly depends on aperture and pressure. AGP of IOM is positively correlated with width of inorganic pore. During reservoir depletion, the AGP of IOM decreases initially and then increases, and the pressure of turning point approximately equals 5. This is because the impact of nanoscale effects exceeds impact of porous deformation at low pressure condition. Impact of real gas effect on AGP of IOM can be ignored.

From Figure 11, AGP of IOM is negatively correlated with irreducible water saturation. High water saturation means thicker water film, which leads to decrease of effective pore width and AGP. When the irreducible water saturation reaches 0.5, the AGP reduces about 95%. It should be noted that the decrease rate of AGP significantly mitigates. This feature is caused by nanoscale effect, which enhances gas transport capacity of IOM at high water saturation condition.

From Figure 12, the AGP increases as the increase of aspect ratio, when the width of pore and irreducible water saturation in IOM is constant. The higher aspect ratio of inorganic pore means the greater cross-sectional area, which leads to higher gas transport capacity.

From Figure 13, both porous deformation and irreducible water saturation have negative influence on AGP in IOM. Deviation between total AGP and AGP without considering porous deformation increases with the decrease of pressure, which achieves an average of 44.47%. When the saturation in inorganic pore equals 0.3, the deviation compared with AGP without considering water saturation achieves an average of 333.77%. Thus, AGP in inorganic pores is affected by water saturation and porous deformation significantly.

4.3. Gas Transport in Shale Matrix

Actually, aperture of nanopore in shale matrix is not uniform. Gas transport capacity and irreducible water distribution in nanopore with diverse aperture is quite different. PSD of shale sample published by Kuila et al. [44] is substituted into Equation (46) to investigate characteristics of AGP for gas transport in shale matrix. The shale sample is from the Baltic Basin shales, the total organic carbon (TOC) is 5%, the total porosity is 6.8% and volume ratio of organic pore is 23.4%. Figure 14 exhibits the PSD of organic and inorganic pore, and PSD of the two type pores obey similar logarithmic normal distribution. Equivalent hydraulic diameter of nanopore in ORM is 10.2 nm, effective hydraulic width of nanopore in IOM is 75.94 nm and the aspect ratio equals 4.

From Figure 15, irreducible water preferentially occupies small inorganic pore, water content in small pores is greater than that in large pore. The irreducible water distribution in IOM is similar to the NMR T2 spectrum from literature [45], so the data calculated by theory model is reasonable. Pores with small width are blocked by irreducible water at a certain water saturation. For instance, nanopores with width smaller than 10 nm are fully blocked by water, when irreducible water saturation in shale matrix equals 0.3. For large pores, water film reduces effective width and porosity, so irreducible water has little impact on large pore.

From Figure 16, the effective PSD of inorganic pores with different irreducible water saturation also obeys approximate logarithmic normal distribution. The range and peak value of effective PSD decreases as water saturation. According to above analysis, pressure, aperture and irreducible water saturation are important effects which have a significant impact on AGP. In this part, we will calculate and analyze characteristics of AGP of actual shale matrix at different conditions.

From Figure 17a, AGP of shale matrix exhibits shape of approximate “V” during reservoir depletion. AGP decreases as irreducible water saturation increases; it should be noted that decrease rate significantly magnifies. This feature of decrease rate is contrary to the result calculated by mean aperture (Figure 11). From Figure 17b, relative permeability is negatively correlated with irreducible water saturation; the AGP decreases rapidly when the water saturation exceeds 0.2. Meanwhile, pressure has little impact on relative permeability at certain water saturation.

From Figure 18, the contribution of ORM to total AGP is positively correlated with irreducible water saturation, and which is negatively correlated with pore pressure. The contribution of ORM is much smaller than IOM due to the aperture and porosity in IOM is much larger than that of ORM. The contribution of ORM exceeds 0.1 only when the water saturation was greater than 0.6, so inorganic pores in IOM are the dominated flow channels.

From Figure 19, AGP calculated by pore size distribution is less than AGP from equivalent hydraulic aperture when the shale matrix does not contain water. The results calculated by the two apertures are similar at high pressure condition, while the deviation is obvious at low pressure condition. The AGP was mainly contributed by pores with large aperture in shale matrix. This characteristic is in agreement with the experimental results conducted by Moghaddam et al. [46]. Ignoring the influence of PSD in shale matrix will overestimate transport capacity of shale matrix.

From Figure 20, relative permeability calculated by equivalent aperture is far less than relative permeability calculated by PSD at the same water saturation. In actual shale matrix, irreducible water blocks small pores and has little impact on pores with large aperture. AGP of actual shale matrix decreases with small amplitude at low water saturation. Therefore, ignoring distribution of irreducible water in shale matrix will underestimate the gas transport capacity significantly. The deviation between two types aperture increases with the increase of water saturation. When irreducible water saturation equals 0.3, the underestimation of AGP calculated by equivalent aperture achieves 71.52%.

Above derivation, water saturation is a significant impact for AGP of shale matrix. Nowadays, the combination of horizontal well and hydraulic fracturing is the main exploration pattern. Massive fracture fluid is captured in the formation after hydraulic fracturing work, which may cause formation damage. Therefore, using liquid nitrogen and supercritical carbon hydroxide as fracture fluid can effectively avoid reservoir damage [47]. Compared with slick water, water-free fracture fluid could decrease environment pollution.

5. Conclusions

In this paper, a novel AGP model is derived to describe gas transport through shale matrix. Compared with previous models, the proposed model considers the impact of pore size distribution and irreducible water distribution on gas flow, which can better reveal gas nonlinear transport behavior in shale matrix. The proposed AGP model can provide help to a certain extent for numerical simulation, high-quality reservoir evaluation and production prediction. According to parameters sensitivity analysis, the main conclusions can be made:

(1)

During shale reservoir depletion, slip flow is the main transport mechanism in shale matrix. Surface diffusion and Knudsen diffusion mainly exist in small pores, especially at low pressure condition.

(2)

Under the comprehensive impact of porous deformation and nanoscale effect, AGP curve of shale matrix mainly exhibits approximate shape of “V” during pressure decreases.

(3)

Pore pressure, PSD, porous deformation and irreducible water saturation have significant impact on AGP. Real gas effect has little impact on AGP, which can be ignored in engineering application.

(4)

Irreducible water preferentially occupies small inorganic pore, and water content in small pores is higher than in large pores. At different irreducible water saturation, effective PSD of inorganic pores obeys approximate logarithmic normal distribution. The range and peak value of effective PSD decreases as water saturation increases.

(5)

For shale matrix, inorganic pores are the dominated gas transport channel. AGP decreases as irreducible water saturation increases, and the decrease rate significantly magnifies. Neglecting irreducible water distribution will lead to large error for AGP evaluation.