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

: Apparent gas permeability (AGP) is a signiﬁcantly important parameter for productivity prediction and reservoir simulation. However, the inﬂuence of multiscale e ﬀ ect 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 e ﬀ ects such as real gas e ﬀ ect, multiscale e ﬀ ect, 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 signiﬁcant impact on AGP. Finally, e ﬀ ective pore size distribution (PSD) and AGP under di ﬀ erent 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.


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 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.

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.
Energies 2020, 13, x 3 of 24 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.

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.

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:

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: p pr = p/p pc , T pr = T/T pc (2) Energies 2020, 13 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]: 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]: 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).

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.

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]: where n is the number of organic pores per unit area, 1/m 2 ; k t0 is the intrinsic permeability of ORM, m 2 ; φ t0 is the porosity of ORM at the initial condition, dimensionless; r 0 is the hydraulic radius of ORM at initial conditions, m; τ is the tortuosity, dimensionless; η org is volume content of organic pores, dimensionless; φ 0 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]: where ψ t is the compressibility coefficient, pa −1 ; k stress is the permeability of ORM considering porous deformation, m 2 ; 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: Rearranging Equation (9), the hydraulic radius of organic pores can be obtained as: 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]: 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.
Energies 2020, 13, x 5 of 24 Substituting Equations (6) and (7) into (8), and ratio of permeability at different effective tress can be expressed as follows: Rearranging Equation (9), the hydraulic radius of organic pores can be obtained as: where rstress 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]: where θ is the Gas coverage in surface of organic pores, dimensionless; pL is the Langmuir pressure, Pa; dm is the molecular diameter of gas, 0.38 nm; rad 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: Porosity of ORM is directly proportional to the square of equivalent radius where reff is effective radius, m; ϕeff-or is the effective porosity of ORM, dimensionless.

Flow Channel for IOM
The pores in IOM mostly present slit type, rectangular pore with the same height and crosssectional area is used to represent inorganic pores in this work ( Figure 3).  Effective radius of organic pore accounting for coupling influence of porous deformation and adsorption can be expressed as follows: Porosity of ORM is directly proportional to the square of equivalent radius where r eff is effective radius, m; φ eff-or is the effective porosity of ORM, dimensionless.

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).
where S is the cross-sectional area, m 2 ; H 0 is the width of inorganic pores, m; W 0 is the height of inorganic pore, m; ξ is the aspect ratio, dimensionless. where S is the cross-sectional area, m 2 ; H0 is the width of inorganic pores, m; W0 is the height of inorganic pore, m; ξ is the aspect ratio, dimensionless.
The aspect ratio is defined as follows: 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: where Sw 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: where ϕs0 is the porosity of IOM at the initial condition, dimensionless; Heff is effective width of inorganic pore, m; ϕeff-IOM is the effective porosity of IOM, dimensionless.

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. The aspect ratio is defined as follows:

Surface Diffusion
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: 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: 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.

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.

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]: where J s is mass flux of surface diffusion in organic pore Kg/s; A s is annulus area of absorbed gas layer, m 2 ; D s is the surface diffusion coefficient, m 2 /s; C smax is the maximum absorbed gas concentration of shale sample, mol/m 3 . Expression of annulus area of absorbed gas layer and the maximum concentration of absorbed gas in shale sample can be written: where ρ s is density of shale, kg/m 3 ; V L is the Langmuir volume, m 3 /m 3 ; V std is molar volume of gas, m 3 /mol. According to the equation, surface diffusion flux is given in the form of pressure gradient: where ∆p is the pressure gradient of along the capillary, Pa. The surface diffusion coefficient D s is defined as [40]: where H is the equivalent adsorption heat, J/mol; D s0 is surface diffusion coefficient at zero gas coverage, m 2 /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: where ρ b is the density of free gas, Kg/m 3 .

Bulk Gas Flow
Knudsen number (Kn) of gas within organic pores is written as: 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]: where k slip-tube is AGP of slip flow within organic pore, m 2 . 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]: where δ is the surface roughness of pore, dimensionless; k Kn-tube is AGP of Knudsen diffusion within organic pore, m 2 ; 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]: 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: (32) where k bt is AGP of bulk gas within organic pore, m 2 .

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: where Q is the gas flux in ORM, m 3 /s; L is the characteristic length of ORM, m.
The flux in ORM can be also expressed as the form of Darcy's law: where k app-tube is AGP for ORM, m 2 . Substituting Equations (7) and (33) into (34), the AGP of ORM coupling multiple transport mechanisms and physical effects can be derived as: Equation (35) encompasses multiple gas transport mechanisms in organic matrix, which can describe gas transport behavior for all flow regime.

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: 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]: where k slip-slit is AGP of slip flow within inorganic pore, m 2 ; A(ξ) is the section shape factor for slip flow, dimensionless. Expression of A(ξ) is shown as follows: Considering impact of surface roughness, the AGP of Knudsen diffusion for real gas transports in inorganic pores is given by [16]: where k kn-slit is AGP of Knudsen diffusion within inorganic pore, m 2 ; B(ξ) is the section shape factor for Knudsen diffusion dimensionless; Its expression is shown as follows: Based on molecular collision frequency, the expressions of weighting coefficient for diverse transport mechanisms in inorganic pore are given as: 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: (42) where k app-slit is AGP for IOM, m 2 .

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: where k app-org is the AGP for ORM considering multiscale, m 2 ; k app-iom is the AGP for IOM considering multiscale, m 2 ; 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: where k app is the AGP of shale matrix, m 2 . The initial irreducible water saturation of shale matrix can be expressed as follows: 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.

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.

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.

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.

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.

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

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.  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.

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.  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 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. 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.  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. 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.

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.  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.

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 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.

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 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.

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 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.
Energies 2020, 13, x 16 of 24 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.

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

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. 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 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 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 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. Energies 2020, 13, x 18 of 24 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 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 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 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. Energies 2020, 13, x 19 of 24 Figure 18. The contribution of ORM to total AGP varies water saturation.
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%. 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 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%. 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.

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.  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.

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.
Author Contributions: Conceptualization, X.L., S.L. and X.T.; methodology and investigation Y.L. and J.L.; writing-original draft preparation, S.L. and X.T.; writing-review and editing, X.L., X.T. and F.W. All authors have read and agreed to the published version of the manuscript.