Gas Multiple Flow Mechanisms and Apparent Permeability Evaluation in Shale Reservoirs

: Gas ﬂow mechanisms and apparent permeability are important factors for predicating gas production in shale reservoirs. In this study, an apparent permeability model for describing gas multiple ﬂow mechanisms in nanopores is developed and incorporated into the COMSOL solver. In addition, a dynamic permeability equation is proposed to analyze the effects of matrix shrinkage and stress sensitivity. The results indicate that pore size enlargement increases gas seepage capacity of a shale reservoir. Compared to conventional reservoirs, the ratio of apparent permeability to Darcy permeability is higher by about 1–2 orders of magnitude in small pores (1–10 nm) and at low pressures (0–5 MPa) due to multiple ﬂow mechanisms. Flow mechanisms mainly include surface diffusion, Knudsen diffusion, and skip ﬂow. Its weight is affected by pore size, reservoir pressure, and temperature, especially pore size ranging from 1 nm to 5 nm and reservoir pressures below 5 MPa. The combined effects of matrix shrinkage and stress sensitivity induce nanopores closure. Therefore, permeability declines about 1 order of magnitude compare to initial apparent permeability. The results also show that permeability should be adjusted during gas production to ensure a better accuracy.


Introduction
The stimulated shale reservoir is referred to as an inherently heterogeneous media benefited from the successful implementation of hydraulic fracturing [1][2][3][4][5]. Shale reservoirs contain organic matter (kerogen, bitumen) and inorganic matter (quartz, clay) [6,7]. Most of the shale gas exists in the adsorbed state on the nanopore surface of organic matrix due to a large specific surface area. However, there is little flow in pores and fractures in a free state ( Figure 1) and only an insignificant amount of shale gas is dissolved in organic kerogen [5][6][7][8]. Atomic force microscopy (AFM) and scanning electron microscopy (SEM) reveal complex multi-scale pore structure of shale [8][9][10]. Shale nanopores are generally smaller than 10 nm, which is similar to the gas molecule free path [11,12]. Considering the collision between gas molecules and the collisions between gas molecules and the surface of nanopore walls, gas flow in shale nanopores exhibits the non-continuum phenomenon [13]. Therefore, shale reservoirs have complex structures and flow mechanisms.
The research method of the transport mechanisms was based on the Navier-Stokes equation with an idealized channel flow model or tube flow model [14,15]. Other modeling approaches, such as Lattice-Boltzmann (LBM) [16], direct simulation Monte Carlo (DSMC) [16], and dynamic Sustainability 2019, 11, 2114 2 of 21 molecular (MD) [17,18] can be utilized to model gas flow in nanopores. However, considering the large scale of the hydraulic fracturing area and the macroscopic pores and fractures, gas flow mechanisms in a shale reservoir are generally written in a format based on an extended N-S equation [19], and apparent permeability model has been widely applied [20,21]. In addition, the characterization of gas flow mechanisms frequently depend on the Knudsen number (Kn is the ratio of the gas molecule free path to the pore diameter) [22][23][24][25][26][27][28][29][30] to evaluate whether Knudsen diffusion, slip flow, transitional flow, or continuous flow exists in nanopores ( Figure 2). Researches on gas flow in narrow and closed spaces dated back to the works of Maxwell [21] and Knudsen [22], which described the slippage effect near the solid surface as a consequence of slip flow, and the models were widely modified and applied. [23]. Klinkenberg [24] was the first to consider the gas slip phenomena and derived an apparent permeability expression. Ertekin [25] developed a dual-mechanism transport model under the influence of a pressure field and a concentration field. The transport process includes the Knudsen flow and the slip flow. Similarly, Javadpour et al. [12,26] proposed an apparent permeability model that describe Knudsen diffusion and slip flow on the nanopore surfaces. Beskok and Karniadakis [15] deduced a unified Hagen-Poiseuille-type equation covering the fundamental flow regimes including free molecule flow, transition flow, slip flow, and viscous flow conditions in channels at microand nano-scales for the first time. The model describes four flow regimes in one equation based on the Knudsen number. Afterwards, Civan [27] developed a permeability equation same to the Beskok and Karniadakis model. In addition, the approach rigorously accounts for the effect of several empirical parameters in porous media including intrinsic permeability, rarefaction coefficient, and a slippage factor [27,28]. Surface diffusion was added to a portion of the gas flow in the nanopore in more recent studies and the driving force was a chemical potential gradient [29,30]. In addition, Wu et al. [30,31] contributed reasonable weighting coefficients to value complex transport mechanisms including surface diffusion. Previous studies have gradually analyzed gas flow mechanisms in nanopores. However, few investigations have been done on the separate characterizations of gas multiple flow mechanisms and their models have failed to consider the influence of pore size changes on the flow mechanisms and permeability during gas production.
Shale nanopore size is influenced by combined effects of matrix shrinkage and stress sensitivity [32][33][34][35][36][37]. Thus, apparent permeability and gas multiple flow mechanisms are dynamically changing along with gas production. The adsorbed gas desorbs from nanopore surfaces leading matrix shrinkage and resulting in an increase in the pore size and apparent permeability [32][33][34]. In contrast, the effective stress of the matrix skeleton increases as gas in nanopores flow out, causing the pore to be compressed [35,36]. Dong [35] and Mckee's [37] research results have shown that a decline in permeability has a negative exponential relationship with the effective stress. Therefore, the effects of matrix shrinkage and stress sensitivity are in need for the accurate characterization of shale apparent permeability.
In this article, gas multiple flow mechanisms are characterized by the separate equations in an apparent permeability model. Specifically, the combined effects of matrix shrinkage and stress sensitivity on pore size and permeability evaluation are also incorporated into COMSOL solver. For better understanding of gas multiple flow mechanisms, surface diffusion, slip flow, transitional diffusion, and Knudsen diffusion are described in the apparent permeability model. In addition, kerogen diffusion can be negligible in terms of the extremely low diffusion coefficient. Stress sensitivity and matrix shrinkage act in opposite directions in the changes of pore size. Thus a dynamic permeability equation associated with pore size is developed. As an extension of our preliminary study, comprehensive simulations including more complete physical process characterization and influenced factors can be performed to access more accurate reservoir information.
The rest of this article is organized as follows: in Section 2, the concept of gas multiple flow mechanisms is presented. In Section 3, the apparent permeability equation and dynamic permeability equation are proposed. In Section 4, the simulation model in COMSOL is constructed. In Section 5, the results are displayed and discussed. The conclusion is presented at last.  [31] and Chong et al. [38]).

Shale Gas Multiple Flow Mechanisms
The pore size is in the order of nanometers and the rate of gas flow is slow in shale reservoirs [39], resulting in an increase in the frequency of collisions between gas molecules and pore walls. Thus, gas flow characteristics in a low permeability shale reservoir are different from those of liquids in a conventional reservoir. Thus, gas flow gas multiple flow mechanisms in nanopores need to be characterized.

Kerogen Diffusion
Kerogen diffusion involves the flow of shale gas dissolved in the organic kerogen nanopore under gas concentration [40] (Figure 1). Unfortunately, the gas diffusion rate is relatively low in kerogen in terms of diffusion coefficient and total dissolved gas content. Fick's Second Diffusion Law [41] is usually applied to model gas diffusion in the kerogen matrix:  [31] and Chong et al. [38]).

Shale Gas Multiple Flow Mechanisms
The pore size is in the order of nanometers and the rate of gas flow is slow in shale reservoirs [39], resulting in an increase in the frequency of collisions between gas molecules and pore walls. Thus, gas flow characteristics in a low permeability shale reservoir are different from those of liquids in a conventional reservoir. Thus, gas flow gas multiple flow mechanisms in nanopores need to be characterized.

Kerogen Diffusion
Kerogen diffusion involves the flow of shale gas dissolved in the organic kerogen nanopore under gas concentration [40] (Figure 1). Unfortunately, the gas diffusion rate is relatively low in kerogen in terms of diffusion coefficient and total dissolved gas content. Fick's Second Diffusion Law [41] is usually applied to model gas diffusion in the kerogen matrix:

Shale Gas Multiple Flow Mechanisms
The pore size is in the order of nanometers and the rate of gas flow is slow in shale reservoirs [39], resulting in an increase in the frequency of collisions between gas molecules and pore walls. Thus, gas flow characteristics in a low permeability shale reservoir are different from those of liquids in a conventional reservoir. Thus, gas flow gas multiple flow mechanisms in nanopores need to be characterized.
where λ is diffusion coefficient in kerogen in m 2 /s, and c is the gas concentration in mol/m 3 .

Surface Diffusion
Surface diffusion is the phenomenon of adsorbed gas flowing along the nanopore surface ( Figure 1). The driving force is the concentration gradient of adsorbed gas [30,31]. Shale matrix has a large specific surface area and a high concentration gradient. Therefore, surface diffusion is a non-neglected transport mechanism in nanopores. Previous work indicates that in the presence of surface diffusion, gas apparent permeability can be several orders of magnitude higher than that of conventional reservoirs [30,31]. In addition, the driving force of surface diffusion can be described by the Maxwell−Stefan method [42]: The surface diffusion coefficient D s can be obtained in combination with experimental data from gas adsorption on nanotubes [43]: where ∆H is the isostatic adsorption heat in J/mol, R is a molar gas constant (8.314 J/(mol · K)), and T is the reservoir temperature in K. The correction factor of surface diffusion ξ ms is related to the porosity φ and tortuosity τ [44]: where d M is the diameter of methane molecular in nm, r is the pore radius in m, φ is the matrix porosity (dimensionless), and τ is the pore tortuosity (dimensionless). The Langmuir monolayer adsorption concentration c s can be defined as [30]: where θ is the extent of the gas coverage on the pore wall for real gas (dimensionless), M is the molar mass in kg/m 3 , d M is the diameter of methane molecular in m, N A is Avogadro's constant (6.02 × 10 23 /mol).

Gas Flow Regimes in Nanopores
The diameters of shale pores below 2 nm are known as micropores and greater than 50 nm are macropores. The diameters in mesopores ranges from 2 to 50 nm [8,28]. The pore size in the shale ranges from nanometer to micrometer. However, organic and inorganic micro-nano pores provide a large specific surface area for gas accumulation, and micro-nano channels increase the frequency of collisions between gas molecules and pore walls. Hence, the gas flow regimes are different from traditional subsurface fluid viscous flow or molecular flow. The Knudsen number Kn is widely used to describe the gas flow regimes in shale nanopores. Continuum fluid flow (Kn ≤ 0.01), slip flow (0.01 < Kn < 0.1), transition flow (0.1 < Kn < 10), and free molecular flow (Kn ≥ 10) are presented in Figure 2. As shown, the gas flow regime is mainly dominated by Knudsen diffusion and transition diffusion in shale  Figures 1 and 2). With an increase in pore size or reservoir pressure, the collision frequency of gas molecules is enhanced. The gas flow regime is mainly transformed into slip flow and continuum flow. The Knudsen number Kn is defined as the ratio of molecular free path λ to pore diameter d [12,15,31]: The mean-free-path of gas molecules λ is given as where µ is the gas viscosity in Pa·s, p is the reservoir pressure in MPa, R is a molar gas constant (8.314 J/(mol· K) ), T is the reservoir temperature in K, and M is the molar mass in kg/m 3 , d is the pore diameter in nm. When gas flow in macropores (pore diameter > 50 nm) and the molecule free path is short under the high pressure, the frequency of intermolecular collisions is strongly high (Kn < 0.01) ( Figure 2). The flow mechanism is governed by the continuum flow and this transport is driven by a pressure gradient. Therefore, the traditional fluid dynamics equations can be used to describe fluid flow in porous media [45]: where k is the Darcy permeability in m 2 , ρ is the gas density in kg/m 3 , µ is the gas viscosity in Pa·s. When gas flow in macropores and mesopore (pore diameter > 2 nm) and the frequency of intermolecular collisions is high (Kn < 0.1) (Figure 2). This transport is still dominated by the collisions between gas molecules and is referred to as slip flow. Though the traditional hydrodynamic equations can be applied, the correction coefficient can more accurately describe the gas slip on the surface of the wall. Klinkenberg [24] published a slip factor to modify the permeability of mass flux: where ρ is the gas density in kg/m 3 , µ is the gas viscosity in Pa·s, λ is the molecular free path in m, b is a constant (dimensionless), and p is the reservoir pressure in MPa. The slip effect on the surface of nanopores accelerates gas flow [46]. Brown [47] also provided a theoretical dimensionless coefficient F to correct the slip velocity based on Poiseuille's Equation: where r is the pore radius in m, ρ is the gas density in kg/m 3 , µ is the gas viscosity in Pa·s, α is a constant (0.8), R is a molar gas constant (8.314 J/(mol· K) ), T is the reservoir temperature in K, and M is the molar mass in kg/m 3 . When 0.1 < Kn < 10, slip flow and Knudsen diffusion occur simultaneously and traditional hydrodynamic equations cannot describe this flow regime. Fortunately, this transport is usually corrected by Knudsen diffusion and skip flow, which is called transitional flow. In addition, transitional flow generally exists in tight gas reservoirs and it can be proposed as the combination of Knudsen diffusion and skip flow.
Knudsen flow or free molecular flow occurs under a low-pressure gradient when gas flow in micropores (pore diameter below 2 nm). Thus, the gas molecule free path is high and this transport is dominated by collisions between gas molecules and pore walls. The diffusion equation is established by gas pressure gradient and diffusion constant [12,48].
where D K is the Knudsen diffusion constant, which is defined as where r is the pore radius in m, R is a molar gas constant (8.314 J/(mol· K) ), T is the reservoir temperature in K, and M is the molar mass in kg/m 3 .

Model Establishment and Analysis
In this section, an apparent permeability model for describing gas multiple flow mechanisms is developed. In addition, a dynamic permeability equation that considers the matrix shrinkage and stress sensitivity effects is derived.

Establishment of Apparent Permeability Model
In the typical shale reservoirs, the Knudsen number is in the range of 0.01-6 [30]. Thus, the multiple flow mechanisms in shale nanopores include kerogen diffusion, surface diffusion, Knudsen diffusion, transmit flow, and slip flow, which are described in Section 2. These transport mechanisms should be presented in the model so that the evolution law and influencing factors of shale permeability can be more accurately analyzed. However, kerogen diffusion can be negligible at the stimulated reservoir scale in terms of the extremely low diffusion coefficient and total dissolved gas content. In addition, transmit flow is the mixing of Knudsen diffusion and slip flow. Therefore, all flow mechanisms in nanopores can be described by the combination of surface diffusion (Equation (2)), slip flow (Equation (9)), and Knudsen diffusion (Equation (12)). The total gas flux through the nanopore is developed as To model the gas permeability, the total gas flux needs to be converted to a form of permeability. Thus, the relationship between the gas flow rate and permeability in shale reservoir can be described as [12,45] The Darcy flow rate equation [45] is usually applied for characterizing gas flow rate in conventional reservoirs and partially tight reservoirs. Similarly, considering the meso-and macropores are the main flow channel of shale gas, the gas flow in shale reservoir can be determined in Equation (16): By substituting Equations (16) and (14) into Equation (15), the characterized term of gas seepage capacity in shale reservoirs called apparent permeability k app is obtained as Note that the characterization of the gas flow capacity is shown in Equation (17). The apparent permeability k app is a more accurate expression and interpretation of shale permeability, similar to the Darcy permeability in conventional reservoirs. In addition, the apparent permeability depends on the properties of shale gas at the different reservoir pressure, temperature, and pore size. This dependence is expected because the interaction of gas molecules with the pore surface and multiple flow mechanisms in nanopores are presented in this model.
The gas flow channel can be conceptualized as a circular tube for porous media [45]. The permeability relies on the size of the flow space. Mathematically, fluid flow in a circular tube or pore can be empirically given by Hagen-Poiseuille's equation [49]. Thus, the Darcy permeability k d of porous media can be assumed as [49] Previous research results shown that the production of tight gas or shale gas is higher than expected for gas production from these reservoirs [12,26]. Furthermore, laboratory tests indicated that shale permeability is higher than expected from these reservoirs [49]. The reason is the fact that the Darcy permeability equation of conventional reservoir can hardly evaluate the production rate of tight gas. Comparing to conventional reservoir, the apparent permeability in shale reservoirs is affected by geological conditions and multiple flow mechanisms. To address this problem, the apparent permeability is divided by Darcy permeability type [45]: where D s is the surface diffusion coefficient in m 2 /s, C s is the Langmuir monolayer adsorption concentration in kg/m 3 , ξ ms is the correction factor of surface diffusion (dimensionless), r is the pore radius in m, µ is the gas viscosity in Pa·s, R is a molar gas constant (8.314 J/(mol· K) ), T is the reservoir temperature in K, and M is the molar mass in kg/m 3 , b is the gas slip factor (dimensionless), and p is the reservoir pressure in MPa. Note that the Equation (19) shows that the apparent permeability in shale reservoirs is larger than the Darcy permeability expected from these reservoirs. In addition, Equation (19) suggests that the smaller the size of the pores, the larger the difference between apparent permeability and Darcy permeability. This equation also shows that the deviation between the apparent permeability and Darcy permeability is larger at lower pressures. Therefore, apparent permeability is similar to Darcy permeability by increasing the pore size or increasing reservoir pressure.

Permeability Dynamic Changes under Matrix Shrinkage and Stress Sensitivity
The following introduces the effects of matrix shrinkage and stress sensitivity on the pore size to address the apparent permeability dynamic evaluation rather than remain constant. Shale apparent permeability shows a high sensitivity to pore size and the change in pore size is highly related to the deformation of the matrix skeleton [50][51][52]. In shale reservoirs, during pressure drawdown of the reservoir by primary production, the nanopore is compressed and the permeability is reduced due to effective stress increases. In contrast, pressure drawdown results in gas desorption, and this is accompanied by matrix shrinkage, which enlarges the nanopore and leads permeability to increase ( Figure 3). Thus, the change in nanopore opening and shale permeability could be positive (increase) or negative (reduction). Sustainability 2019, 11, x FOR PEER REVIEW 8 of 21 The research results from shale reservoirs show that the content of adsorbed gas in shale accounts for 20-85% [53,54]. Furthermore, nanopores in shale reservoirs are mostly oil-wet, resulting in a significant adsorption capacity. The adsorption of shale gas satisfies the Langmuir isotherm adsorption curve [55]. The content of gas adsorption is associated with the reservoir pressure p , the Langmuir limiting adsorption pressure L P , and the Langmuir limiting adsorption volume L V [55], which can be expressed as where V is the gas adsorption amount in 3 cm /g , L V is the Langmuir volume in 3 cm /g , L P is the Langmuir pressure in MPa , and p is the adsorption-desorption dynamic equilibrium pressure in MPa .
Adsorption gas gradually desorbs from the matrix pore surface as a decrease in reservoir pressure, resulting in matrix shrinkage and pore size enlargement, this is simply called matrix positive strain effect. Fortunately, the strain of matrix shrinkage ε Δ caused by the changes in the adsorbed gas content V can be described by the Bangham solid deformation theory [56].
V dp EV p (21) where ρ is the shale density in 3 kg/m , R is a molar gas constant ( 8.314 J/(mol K) ⋅ ), T is the reservoir temperature in K , p is the reservoir pressure in MPa , E is the Young's modulus in GPa , V is the gas adsorption amount in 3 cm /g , and L P is the Langmuir pressure in MPa . Seidle [57,58] presented a model describing the effect of matrix shrinkage deformation on the reservoir porosity. Therefore, the relationship between the change of porosity φ and matrix strain ε Δ can be proposed as A large number of nanopores in shale can be assumed to be a series of capillary tubes [58]. The mass flow of ideal gas in a circular tube is derived from the Hagen-Poiseuille law [49,58]. Therefore, the dynamic changes in permeability and porosity is proposed. The research results from shale reservoirs show that the content of adsorbed gas in shale accounts for 20-85% [53,54]. Furthermore, nanopores in shale reservoirs are mostly oil-wet, resulting in a significant adsorption capacity. The adsorption of shale gas satisfies the Langmuir isotherm adsorption curve [55]. The content of gas adsorption is associated with the reservoir pressure p, the Langmuir limiting adsorption pressure P L , and the Langmuir limiting adsorption volume V L [55], which can be expressed as where V is the gas adsorption amount in cm 3 /g, V L is the Langmuir volume in cm 3 /g, P L is the Langmuir pressure in MPa, and p is the adsorption-desorption dynamic equilibrium pressure in MPa. Adsorption gas gradually desorbs from the matrix pore surface as a decrease in reservoir pressure, resulting in matrix shrinkage and pore size enlargement, this is simply called matrix positive strain effect. Fortunately, the strain of matrix shrinkage ∆ε caused by the changes in the adsorbed gas content V can be described by the Bangham solid deformation theory [56].
where ρ is the shale density in kg/m 3 , R is a molar gas constant (8.314 J/(mol· K) ), T is the reservoir temperature in K, p is the reservoir pressure in MPa, E is the Young's modulus in GPa, V is the gas adsorption amount in cm 3 /g, and P L is the Langmuir pressure in MPa. Seidle [57,58] presented a model describing the effect of matrix shrinkage deformation on the reservoir porosity. Therefore, the relationship between the change of porosity φ and matrix strain ∆ε can be proposed as A large number of nanopores in shale can be assumed to be a series of capillary tubes [58]. The mass flow of ideal gas in a circular tube is derived from the Hagen-Poiseuille law [49,58]. Therefore, the dynamic changes in permeability and porosity is proposed. where k 0 is the initial permeability in m 2 , φ 0 is the initial porosity (dimensionless), r 0 is the initial pore size in nm, and µ is the gas viscosity in Pa·s. Substituting Equation (21) and Equation (22) into Equation (24), gas desorption from matrix shrinkage initiates the changes in pore size can be described as where ρ is the shale density in kg/m 3 , R is a molar gas constant (8.314 J/mol·K), T is the reservoir temperature in K, p is the reservoir pressure in MPa, E is the Young's modulus in GPa, V is the gas adsorption amount in cm 3 /g, and P L is the Langmuir pressure in MPa, k 0 is the initial permeability in m 2 , φ 0 is the initial porosity (dimensionless), and r 0 is the initial pore size in nm. Shale gas recovery is accompanied by a decrease in the reservoir pressure, resulting in pore closure and a rapid decline in shale permeability. Unfortunately, this is simply called matrix negative strain effect. In addition, the ratio of shale permeability to initial permeability has a negative exponential relationship with the effective stress [35][36][37]: where k 0 is the initial permeability in m 2 , p s is the change of pressure in MPa, and a is a modified coefficient (dimensionless). Combining Equations (24) and (26), effective stress triggers changes in pore size is established as where p s is the pressure change in MPa, a is a modified coefficient (dimensionless), and r 0 is the initial pore radius in nm. Matrix shrinkage and effective stress act in opposite directions on pore radius and permeability. The initial pore radius is r 0 , the pore radius decrease (r 2 − r 0 ) relies on pressure drawdown, and the pore radius increment (r 1 − r 0 ), associated with gas desorption, are considered in dynamic changes of permeability. Therefore, substituting Equations (25) and (27) into Equation (24), the relationship between dynamic permeability k and apparent permeability k app is obtained as

Numerical Model and Module
The gas flow process in a shale reservoir mode was simulated in COMSOL solver. COMSOL software is a multi-physics coupling solver based on commercial PDE (Partial Differential Equation). To simplify the simulation, the shale reservoir is assumed to be blocks, which are divided by network fractures (shown in Figure 1). The size of the simulated model is a square with a length of 1 dm as a characterization unit of shale based on the average spacing of fractures in stimulated shale reservoirs [39,51] (shown in Figure 4).

Physical Boundary Conditions and Input Parameters
The 2D time-dependent Rock Mechanic module and Subsurface Flow module in COMSOL are applied to simulate the gas flow process (Figure 4). Appropriate boundary conditions must be employed for this coupled matrix deformation and gas flow model. For shale deformation during gas production, the vertical pressure (50 MPa) and horizontal pressure (53 MPa) were applied to the top and the right side [39,51], and the left and bottom boundaries of model were both rollered. For gas flow due to the pressure gradient, the initial pore pressure of the reservoir was assumed to be 30 MPa [53,59]. A constant boundary pressure of 10 MPa was applied to the left and right boundaries of the model [60], which were considered the interface between the fracture and the matrix. No flow conditions were set for other boundaries.

Physical Boundary Conditions and Input Parameters
The 2D time-dependent Rock Mechanic module and Subsurface Flow module in COMSOL are applied to simulate the gas flow process (Figure 4). Appropriate boundary conditions must be employed for this coupled matrix deformation and gas flow model. For shale deformation during gas production, the vertical pressure (50 MPa) and horizontal pressure (53 MPa) were applied to the top and the right side [39,51], and the left and bottom boundaries of model were both rollered. For gas flow due to the pressure gradient, the initial pore pressure of the reservoir was assumed to be 30 MPa [53,59]. A constant boundary pressure of 10 MPa was applied to the left and right boundaries of the model [60], which were considered the interface between the fracture and the matrix. No flow conditions were set for other boundaries. The permeability equation was set as a function of variables and parameters, and it was calculated and updated at each time step in COMSOL. The input properties are listed in Table 1. The values of these properties were chosen from the shale reservoirs of Fuling and Barnet [39,51]. Reservoir initial pressure, P ( MPa ) 30 Overburden pressure, y σ ( MPa ) 50 Confining pressure, x σ ( MPa ) 53 Boundary pressure, 0 P ( MPa ) 10 Langmuir pressure constant, L P ( MPa ) 4 Langmuir volume constant, L V ( 3 m kg ) 3 3 10 − × Initial porosity of matrix, 0 φ 0.05 The permeability equation was set as a function of variables and parameters, and it was calculated and updated at each time step in COMSOL. The input properties are listed in Table 1. The values of these properties were chosen from the shale reservoirs of Fuling and Barnet [39,51].

Apparent Permeability Evaluation during Production
In this section, the shale block pressure distribution is shown in Figure 5. With gas production, the pressure of the matrix block near the outflow boundary slowly becomes equal to the fracture pressure (10 MPa). The pore size significantly affects the reservoir pressure distribution. An increased pore size of the shale promotes rapid penetration of gas and the shale block pressure rapidly decreases to the fracture pressure. When the nanopores are micropores and mesopores (pore radius below 2 nm), after 2000 days of production, the pressure in the middle of the model is still higher than the well pressure. This suggests that more intensive fractures should be developed in shale reservoirs to obtain more gas.
The monitoring line is used to show the pressure distribution from the center of the reservoir block to the boundary at different pore sizes and production times (shown in Figure 6). Note that the gas pressure drops significantly at the boundary of the fracture. However, when the nanopores are micropores and mesopores (pore radius below 2 nm), the gas pressure in the block still remains unchanged even the gas seepage into the fracture boundary for a long time. The reason is the fact that gas migrates at a very slow rate due to the low permeability in micropores and mesopores. The pressure distribution at the pore radii ranges of 5 nm and 25 nm are shown in Figure 6c,d, respectively. As shown, the pore size enlargement increases rapidly the reservoir seepage capacity, and the gas flows into the low-pressure region of the fracture in a very short period of time. Thus, it should be noted that the pore size plays an important part in gas flow process and permeability evaluation. In the following section, the evolution of apparent permeability should be analyzed in more detail.

Apparent Permeability Evaluation during Production
In this section, the shale block pressure distribution is shown in Figure 5. With gas production, the pressure of the matrix block near the outflow boundary slowly becomes equal to the fracture pressure (10 MPa). The pore size significantly affects the reservoir pressure distribution. An increased pore size of the shale promotes rapid penetration of gas and the shale block pressure rapidly decreases to the fracture pressure. When the nanopores are micropores and mesopores (pore radius below 2 nm), after 2000 days of production, the pressure in the middle of the model is still higher than the well pressure. This suggests that more intensive fractures should be developed in shale reservoirs to obtain more gas.
The monitoring line is used to show the pressure distribution from the center of the reservoir block to the boundary at different pore sizes and production times (shown in Figure 6). Note that the gas pressure drops significantly at the boundary of the fracture. However, when the nanopores are micropores and mesopores (pore radius below 2 nm), the gas pressure in the block still remains unchanged even the gas seepage into the fracture boundary for a long time. The reason is the fact that gas migrates at a very slow rate due to the low permeability in micropores and mesopores. The pressure distribution at the pore radii ranges of 5 nm and 25 nm are shown in Figure 6c,d, respectively. As shown, the pore size enlargement increases rapidly the reservoir seepage capacity, and the gas flows into the low-pressure region of the fracture in a very short period of time. Thus, it should be noted that the pore size plays an important part in gas flow process and permeability evaluation. In the following section, the evolution of apparent permeability should be analyzed in more detail.   Apparent permeability evaluation positively correlated with changes in pore size is shown in Figure 7. A major difference of the permeability evolution for Figure 7a,b can be observed. As shown in Figure 7a, the pore radius increases from 20 nm to 100 nm, and the permeability is increased by about 20 times. However, the pore size increases from 1 nm to 5 nm, and the permeability is only increased by eight times in Figure 7b. It is found that the change in permeability is inconsistent with the pore size, which means that gas multiple flow mechanisms and geological conditions such as reservoir pressure, temperature, and pore size could affect the permeability evaluation in nanopores. The ratio of apparent permeability to Darcy permeability at different pore sizes and reservoir pressures is shown in Figure 8. As displayed, the  Apparent permeability evaluation positively correlated with changes in pore size is shown in Figure 7. A major difference of the permeability evolution for Figure 7a,b can be observed. As shown in Figure 7a, the pore radius increases from 20 nm to 100 nm, and the permeability is increased by about 20 times. However, the pore size increases from 1 nm to 5 nm, and the permeability is only increased by eight times in Figure 7b. It is found that the change in permeability is inconsistent with the pore size, which means that gas multiple flow mechanisms and geological conditions such as reservoir pressure, temperature, and pore size could affect the permeability evaluation in nanopores.  Apparent permeability evaluation positively correlated with changes in pore size is shown in Figure 7. A major difference of the permeability evolution for Figure 7a,b can be observed. As shown in Figure 7a, the pore radius increases from 20 nm to 100 nm, and the permeability is increased by about 20 times. However, the pore size increases from 1 nm to 5 nm, and the permeability is only increased by eight times in Figure 7b. It is found that the change in permeability is inconsistent with the pore size, which means that gas multiple flow mechanisms and geological conditions such as reservoir pressure, temperature, and pore size could affect the permeability evaluation in nanopores. The ratio of apparent permeability to Darcy permeability at different pore sizes and reservoir pressures is shown in Figure 8. As displayed, the  The ratio of apparent permeability to Darcy permeability at different pore sizes and reservoir pressures is shown in Figure 8. As displayed, the k app : k d ratio is higher in a small radius, especially a pore radius below 5 nm. These results show one to two orders of magnitude change in permeability and explain the unusual gas production from these tight reservoir systems. Thus, this implies that the changes in pore size should be fully considered and shale permeability should be adjusted in the shale reservoirs. Furthermore, with the pore size and reservoir pressure gradually increasing, the shale pore characteristics are similar to those of traditional high permeability reservoirs, and the apparent permeability is gradually reduced to the Darcy permeability. Pore radius as small as 100 nm show no differences between the apparent permeability and the Darcy permeability. This also fits well with the actual situation that the pore size in conventional reservoirs is macroscopic. permeability and explain the unusual gas production from these tight reservoir systems. Thus, this implies that the changes in pore size should be fully considered and shale permeability should be adjusted in the shale reservoirs. Furthermore, with the pore size and reservoir pressure gradually increasing, the shale pore characteristics are similar to those of traditional high permeability reservoirs, and the apparent permeability is gradually reduced to the Darcy permeability. Pore radius as small as 100 nm show no differences between the apparent permeability and the Darcy permeability. This also fits well with the actual situation that the pore size in conventional reservoirs is macroscopic. In addition to the pore size, the reservoir pressure and temperature also influence the ratio of apparent permeability to Darcy permeability ( Figure 9). As shown in Figure 9, increasing temperature results in higher apparent permeability at any pressure step because the collisions between molecules increase. For low reservoir pressure, the mean free path of gas molecules increases, leading to a higher collision frequency between the nanopore surface and gas molecules. Hence the difference between the apparent permeability and the initial Darcy permeability is expanded. The ratio of apparent permeability to Darcy permeability gradually becomes equal to 1 along with the pressure increases. The reason for this is the fact that gas molecules move and collide rapidly and then the gas flow regime is in close proximity to the continuum flow. As a result, the apparent permeability should be carefully adjusted along with the production of shale reservoirs. Simultaneously, lab measurements must be performed under reservoir geological conditions as well as at different pressures and temperatures to generate a data bank of apparent permeability.  In addition to the pore size, the reservoir pressure and temperature also influence the ratio of apparent permeability to Darcy permeability ( Figure 9). As shown in Figure 9, increasing temperature results in higher apparent permeability at any pressure step because the collisions between molecules increase. For low reservoir pressure, the mean free path of gas molecules increases, leading to a higher collision frequency between the nanopore surface and gas molecules. Hence the difference between the apparent permeability and the initial Darcy permeability is expanded. The ratio of apparent permeability to Darcy permeability gradually becomes equal to 1 along with the pressure increases. The reason for this is the fact that gas molecules move and collide rapidly and then the gas flow regime is in close proximity to the continuum flow. As a result, the apparent permeability should be carefully adjusted along with the production of shale reservoirs. Simultaneously, lab measurements must be performed under reservoir geological conditions as well as at different pressures and temperatures to generate a data bank of apparent permeability. permeability and explain the unusual gas production from these tight reservoir systems. Thus, this implies that the changes in pore size should be fully considered and shale permeability should be adjusted in the shale reservoirs. Furthermore, with the pore size and reservoir pressure gradually increasing, the shale pore characteristics are similar to those of traditional high permeability reservoirs, and the apparent permeability is gradually reduced to the Darcy permeability. Pore radius as small as 100 nm show no differences between the apparent permeability and the Darcy permeability. This also fits well with the actual situation that the pore size in conventional reservoirs is macroscopic. In addition to the pore size, the reservoir pressure and temperature also influence the ratio of apparent permeability to Darcy permeability ( Figure 9). As shown in Figure 9, increasing temperature results in higher apparent permeability at any pressure step because the collisions between molecules increase. For low reservoir pressure, the mean free path of gas molecules increases, leading to a higher collision frequency between the nanopore surface and gas molecules. Hence the difference between the apparent permeability and the initial Darcy permeability is expanded. The ratio of apparent permeability to Darcy permeability gradually becomes equal to 1 along with the pressure increases. The reason for this is the fact that gas molecules move and collide rapidly and then the gas flow regime is in close proximity to the continuum flow. As a result, the apparent permeability should be carefully adjusted along with the production of shale reservoirs. Simultaneously, lab measurements must be performed under reservoir geological conditions as well as at different pressures and temperatures to generate a data bank of apparent permeability.

Weight Ratios of Gas Flow Mechanisms and Its Influenced Factors
Based on Equations (14) and (17), the weight ratios for different flow mechanisms as the functions of the pore size, reservoir pressure, and temperature are obtained in this part. These ratios represent the percentage of the contribution of different flow regime to the total flux. As shown in Figure 10, surface diffusion and Knudsen diffusion account for a large proportion of total flux. In particular, surface diffusion dominates when the pore size is below 2 nm. The flow in the nanopores is almost entirely slip flow when the pores are macropores (r > 25 nm). Due to the pore radius being less than 2 nm and a higher concentration of gas molecules on the surface of nanopores, the weight coefficient of surface diffusion increases. A larger pore size leads to an increase in the proportion of slip flow due to the increased frequency of collisions between gas molecules. This also confirms that increased pore size causes the gas flow to transform into a viscous continuous flow. The results in the figure show that the proportion of different flow regimes varies greatly, as the pore size ranges from 1 to 5 nm, indicating the flow mechanism is more sensitive to micropores and mesopores. In addition, the results suggest that the apparent permeability is carefully influenced by gas multiple flow mechanisms in micro-and meso-pores.
The effect of the reservoir pressure on the weight of different flow mechanisms is presented in Figure 11. Under a low pressure condition, the velocity of gas flow slows down, gas molecules accumulate on the surface of the pore walls, thus surface diffusion accounts for a larger proportion. With an increase in reservoir pressure, more gas molecules participate in the collisions between the gas molecules and slip flow dominates. However, Knudsen diffusion is related to the collisions between the gas molecules and pore walls, and its weight increases first and then decreases as the reservoir pressure rises. This is probably a consequence of the pressure transition region between low to high pressure. The results show that when the reservoir pressure is lower than 5 MPa, the flow mechanism changes significantly, which indicates that the apparent permeability is more sensitive to low pressure. Based on Equations (14) and (17), the weight ratios for different flow mechanisms as the functions of the pore size, reservoir pressure, and temperature are obtained in this part. These ratios represent the percentage of the contribution of different flow regime to the total flux.
As shown in Figure 10, surface diffusion and Knudsen diffusion account for a large proportion of total flux. In particular, surface diffusion dominates when the pore size is below 2 nm. The flow in the nanopores is almost entirely slip flow when the pores are macropores (r > 25 nm). Due to the pore radius being less than 2 nm and a higher concentration of gas molecules on the surface of nanopores, the weight coefficient of surface diffusion increases. A larger pore size leads to an increase in the proportion of slip flow due to the increased frequency of collisions between gas molecules. This also confirms that increased pore size causes the gas flow to transform into a viscous continuous flow. The results in the figure show that the proportion of different flow regimes varies greatly, as the pore size ranges from 1 to 5 nm, indicating the flow mechanism is more sensitive to micropores and mesopores. In addition, the results suggest that the apparent permeability is carefully influenced by gas multiple flow mechanisms in micro-and meso-pores.
The effect of the reservoir pressure on the weight of different flow mechanisms is presented in Figure 11. Under a low pressure condition, the velocity of gas flow slows down, gas molecules accumulate on the surface of the pore walls, thus surface diffusion accounts for a larger proportion. With an increase in reservoir pressure, more gas molecules participate in the collisions between the gas molecules and slip flow dominates. However, Knudsen diffusion is related to the collisions between the gas molecules and pore walls, and its weight increases first and then decreases as the reservoir pressure rises. This is probably a consequence of the pressure transition region between low to high pressure. The results show that when the reservoir pressure is lower than 5 MPa, the flow mechanism changes significantly, which indicates that the apparent permeability is more sensitive to low pressure.  Based on Equations (14) and (17), the weight ratios for different flow mechanisms as the functions of the pore size, reservoir pressure, and temperature are obtained in this part. These ratios represent the percentage of the contribution of different flow regime to the total flux.
As shown in Figure 10, surface diffusion and Knudsen diffusion account for a large proportion of total flux. In particular, surface diffusion dominates when the pore size is below 2 nm. The flow in the nanopores is almost entirely slip flow when the pores are macropores (r > 25 nm). Due to the pore radius being less than 2 nm and a higher concentration of gas molecules on the surface of nanopores, the weight coefficient of surface diffusion increases. A larger pore size leads to an increase in the proportion of slip flow due to the increased frequency of collisions between gas molecules. This also confirms that increased pore size causes the gas flow to transform into a viscous continuous flow. The results in the figure show that the proportion of different flow regimes varies greatly, as the pore size ranges from 1 to 5 nm, indicating the flow mechanism is more sensitive to micropores and mesopores. In addition, the results suggest that the apparent permeability is carefully influenced by gas multiple flow mechanisms in micro-and meso-pores.
The effect of the reservoir pressure on the weight of different flow mechanisms is presented in Figure 11. Under a low pressure condition, the velocity of gas flow slows down, gas molecules accumulate on the surface of the pore walls, thus surface diffusion accounts for a larger proportion. With an increase in reservoir pressure, more gas molecules participate in the collisions between the gas molecules and slip flow dominates. However, Knudsen diffusion is related to the collisions between the gas molecules and pore walls, and its weight increases first and then decreases as the reservoir pressure rises. This is probably a consequence of the pressure transition region between low to high pressure. The results show that when the reservoir pressure is lower than 5 MPa, the flow mechanism changes significantly, which indicates that the apparent permeability is more sensitive to low pressure.  The effect of temperature on the contribution of the flow mechanism to the total flux is presented in Figure 12. The result shows that temperature has the lowest effect on the flow regime. According to Equations (6) and (7), a higher reservoir temperature increases the free path of gas molecules and reduces the Knudsen number. The gas flow regime changes from Knudsen flow to slip flow (shown in Figure 2). Thus, the weight coefficient of slip flow increases under the higher temperature in Figure 12. However, surface diffusion shows the least sensitivity to temperature mainly because the diffusion coefficient is only associated with gas concentration on the nanopore surface.
Sustainability 2019, 11, x FOR PEER REVIEW 15 of 21 Figure 11. Contribution of the flow mechanism to the total flux at different pressures. Pore radius = 5 nm, temperature = 350 K.
The effect of temperature on the contribution of the flow mechanism to the total flux is presented in Figure 12. The result shows that temperature has the lowest effect on the flow regime. According to Equations (6) and (7), a higher reservoir temperature increases the free path of gas molecules and reduces the Knudsen number. The gas flow regime changes from Knudsen flow to slip flow (shown in Figure 2). Thus, the weight coefficient of slip flow increases under the higher temperature in Figure 12. However, surface diffusion shows the least sensitivity to temperature mainly because the diffusion coefficient is only associated with gas concentration on the nanopore surface. Figure 12. Contribution of the flow mechanism to the total flux at different temperatures. Pore radius = 5 nm, reservoir pressure = 10 MPa.

Permeability Dynamic Change
In this part, the changes of pore size and permeability caused by matrix shrinkage and stress sensitivity were calculated based on Equation (28). The influence of matrix shrinkage on pore size changes is shown in Figure 13, gas desorbed from the surface causes matrix skeleton shrinkage, resulting in pore size enlargement. The reservoir pressure decreases, leading more gas to desorb from the nanopore wall and a greater increment in the pore size. Thus, gas adsorption significantly causes matrix shrinkage and increases shale permeability. The effect of stress sensitivity on the decrease in the pore radius is presented in Figure 14. As displayed, during depressurization development of shale reservoirs, the effective stress of matrix

Permeability Dynamic Change
In this part, the changes of pore size and permeability caused by matrix shrinkage and stress sensitivity were calculated based on Equation (28). The influence of matrix shrinkage on pore size changes is shown in Figure 13, gas desorbed from the surface causes matrix skeleton shrinkage, resulting in pore size enlargement. The reservoir pressure decreases, leading more gas to desorb from the nanopore wall and a greater increment in the pore size. Thus, gas adsorption significantly causes matrix shrinkage and increases shale permeability.
Sustainability 2019, 11, x FOR PEER REVIEW 15 of 21 Figure 11. Contribution of the flow mechanism to the total flux at different pressures. Pore radius = 5 nm, temperature = 350 K.
The effect of temperature on the contribution of the flow mechanism to the total flux is presented in Figure 12. The result shows that temperature has the lowest effect on the flow regime. According to Equations (6) and (7), a higher reservoir temperature increases the free path of gas molecules and reduces the Knudsen number. The gas flow regime changes from Knudsen flow to slip flow (shown in Figure 2). Thus, the weight coefficient of slip flow increases under the higher temperature in Figure 12. However, surface diffusion shows the least sensitivity to temperature mainly because the diffusion coefficient is only associated with gas concentration on the nanopore surface. Figure 12. Contribution of the flow mechanism to the total flux at different temperatures. Pore radius = 5 nm, reservoir pressure = 10 MPa.

Permeability Dynamic Change
In this part, the changes of pore size and permeability caused by matrix shrinkage and stress sensitivity were calculated based on Equation (28). The influence of matrix shrinkage on pore size changes is shown in Figure 13, gas desorbed from the surface causes matrix skeleton shrinkage, resulting in pore size enlargement. The reservoir pressure decreases, leading more gas to desorb from the nanopore wall and a greater increment in the pore size. Thus, gas adsorption significantly causes matrix shrinkage and increases shale permeability. The effect of stress sensitivity on the decrease in the pore radius is presented in Figure 14. As displayed, during depressurization development of shale reservoirs, the effective stress of matrix The effect of stress sensitivity on the decrease in the pore radius is presented in Figure 14.
As displayed, during depressurization development of shale reservoirs, the effective stress of matrix skeleton increases, resulting in matrix skeleton expansion and pore closure. It is found that the changes of the pore size are insignificant in the period of depressurization, as the pore is micropore (pore radius below 2 nm). The reason is the fact that a smaller nanopore size leads to more contacts among the matrix particles, which decreases the compression of the matrix framework. skeleton increases, resulting in matrix skeleton expansion and pore closure. It is found that the changes of the pore size are insignificant in the period of depressurization, as the pore is micropore (pore radius below 2 nm). The reason is the fact that a smaller nanopore size leads to more contacts among the matrix particles, which decreases the compression of the matrix framework. The apparent permeability and dynamic permeability evaluation are shown in Figure 15. Permeability is positively correlated with the pore size and is more sensitive to large pore sizes, which is also shown in Figure 7. The complex flow mechanisms in micro-and meso-pores lead to inconspicuous changes in permeability. In addition, the changes of permeability are influenced by stress sensitivity and matrix shrinkage. With the decrease in reservoir pressure, gas desorption leads to matrix shrinkage and pore size expansion. In contrast, the pore size decreases due to an increase in the effective stress of matrix skeleton. However, the combined actions of gas desorption and effective stress lead a reduction in the pore diameter. As a result, dynamic permeability is lower than apparent permeability due to pore size reduction. The simulation results show that the dynamic permeability declines monotonically by about 1 order of magnitude compared to the apparent permeability. This also clearly confirms the fact that the permeability of the reservoir is declining during the production process. Thus, the changes in dynamic permeability could more accurately reflect the permeability evaluation during gas production.  The apparent permeability and dynamic permeability evaluation are shown in Figure 15. Permeability is positively correlated with the pore size and is more sensitive to large pore sizes, which is also shown in Figure 7. The complex flow mechanisms in micro-and meso-pores lead to inconspicuous changes in permeability. In addition, the changes of permeability are influenced by stress sensitivity and matrix shrinkage. With the decrease in reservoir pressure, gas desorption leads to matrix shrinkage and pore size expansion. In contrast, the pore size decreases due to an increase in the effective stress of matrix skeleton. However, the combined actions of gas desorption and effective stress lead a reduction in the pore diameter. As a result, dynamic permeability is lower than apparent permeability due to pore size reduction. The simulation results show that the dynamic permeability declines monotonically by about 1 order of magnitude compared to the apparent permeability. This also clearly confirms the fact that the permeability of the reservoir is declining during the production process. Thus, the changes in dynamic permeability could more accurately reflect the permeability evaluation during gas production. skeleton increases, resulting in matrix skeleton expansion and pore closure. It is found that the changes of the pore size are insignificant in the period of depressurization, as the pore is micropore (pore radius below 2 nm). The reason is the fact that a smaller nanopore size leads to more contacts among the matrix particles, which decreases the compression of the matrix framework. The apparent permeability and dynamic permeability evaluation are shown in Figure 15. Permeability is positively correlated with the pore size and is more sensitive to large pore sizes, which is also shown in Figure 7. The complex flow mechanisms in micro-and meso-pores lead to inconspicuous changes in permeability. In addition, the changes of permeability are influenced by stress sensitivity and matrix shrinkage. With the decrease in reservoir pressure, gas desorption leads to matrix shrinkage and pore size expansion. In contrast, the pore size decreases due to an increase in the effective stress of matrix skeleton. However, the combined actions of gas desorption and effective stress lead a reduction in the pore diameter. As a result, dynamic permeability is lower than apparent permeability due to pore size reduction. The simulation results show that the dynamic permeability declines monotonically by about 1 order of magnitude compared to the apparent permeability. This also clearly confirms the fact that the permeability of the reservoir is declining during the production process. Thus, the changes in dynamic permeability could more accurately reflect the permeability evaluation during gas production.  Figure 15. Apparent permeability k app and dynamic permeability k as a function of the pore radius: (a) micro-and meso-pore, (b) meso-and macro-pore.

Comparions of Apparment Permeability Model to Previous Models
In this section, the ratio of apparent permeability to Darcy permeability in this article was compared to the results obtained from previous models with the same settings [15,24]. The ratio was calculated by Equation (19) at the reservoir pressure of 5 MPa and the temperature of 350 K. The comparison results were shown in Figure 16. It is found that the ratios in these models are gradually expanding along the decrease of pore size, which is caused by gas multiple flow mechanisms in micro-and meso-pore (below 20 nm). However, some major differences can be found between these models. As shown, the ration in this paper higher than the results from Klinkenberg's model [24] and Beskok and Karniadaki's model [15] within the pore size range presented in Figure 16. The reason is the fact that gas multiple flow mechanisms are characterized by different equations in this paper, whereas only skip flow is considered in Klinkenberg's model and a single equation describes flow regimes in Beskok and Karniadaki's model. In addition, the ratio in this paper has a more significantly declining trend and the deviation of the ratio is higher when the pore is micropore (below 2 nm). This is caused by the fact that the surface diffusion is added in the part of gas flow mechanisms in this paper compared to other models. Thus, the ration of apparent permeability to Darcy permeability presented in Figure 16 suggests that surface diffusion plays a non-negligible part in the gas flow process, especially if the nanopore is a micropore. Moreover, gas multiple flow mechanisms should be more accurately characterized by different methods in the apparent permeability model.
Compared to previous studies, permeability evaluation caused by stress sensitivity and matrix shrinkage was added and calculated by Equation (28). The ratio of dynamic permeability to apparent permeability that shows the changes of permeability caused by stress sensitivity and matrix shrinkage at different reservoir pressures and pore sizes is presented in Figure 17. As displayed, along with production, the ratio of dynamic permeability to apparent permeability has been consistently below 1 as a consequence of matrix shrinkage and effective stress. However, the ratio drops along with the reservoir pressure and then increases gradually when the nanopore is a micropore (below 2 nm). This is caused by the fact that more shale gas is adsorbed on the pore surface due to the larger the specific surface area in micropores. In addition, the gas flow is slow at low pressure and the differences of stress are low. Thus, the matrix expansion caused by the effective stress should be weakened and matrix shrinkage accounts for the more changes of pore size under low reservoir pressure. As shown, the ratio increases gradually along with reservoir pressure as the pore is meso-and macro-pore (above 5 nm). The reason is the fact that there is less adsorbed gas and higher gas flow velocity in nanopores as the pore size enlarges. Moreover, the particle connection of matrix is weakened and more susceptible to effective stress, resulting in a significant decrease in pore size. The results in Figure 17 suggest that the effect of stress sensitivity and matrix shrinkage should be carefully included in the shale permeability evaluation models. More accurate knowledge of the factors affecting the pore size and permeability may help in the optimization of shale gas production.
was calculated by Equation (19) at the reservoir pressure of 5 MPa and the temperature of 350 K. The comparison results were shown in Figure 16. It is found that the ratios in these models are gradually expanding along the decrease of pore size, which is caused by gas multiple flow mechanisms in micro-and meso-pore (below 20 nm). However, some major differences can be found between these models. As shown, the ration in this paper higher than the results from Klinkenberg's model [24] and Beskok and Karniadaki's model [15] within the pore size range presented in Figure 16. The reason is the fact that gas multiple flow mechanisms are characterized by different equations in this paper, whereas only skip flow is considered in Klinkenberg's model and a single equation describes flow regimes in Beskok and Karniadaki's model. In addition, the ratio in this paper has a more significantly declining trend and the deviation of the ratio is higher when the pore is micropore (below 2 nm). This is caused by the fact that the surface diffusion is added in the part of gas flow mechanisms in this paper compared to other models. Thus, the ration of apparent permeability to Darcy permeability presented in Figure 16 suggests that surface diffusion plays a non-negligible part in the gas flow process, especially if the nanopore is a micropore. Moreover, gas multiple flow mechanisms should be more accurately characterized by different methods in the apparent permeability model.
Compared to previous studies, permeability evaluation caused by stress sensitivity and matrix shrinkage was added and calculated by Equation (28). The ratio of dynamic permeability to apparent permeability that shows the changes of permeability caused by stress sensitivity and matrix shrinkage at different reservoir pressures and pore sizes is presented in Figure 17. As displayed, along with production, the ratio of dynamic permeability to apparent permeability has been consistently below 1 as a consequence of matrix shrinkage and effective stress. However, the ratio drops along with the reservoir pressure and then increases gradually when the nanopore is a micropore (below 2 nm). This is caused by the fact that more shale gas is adsorbed on the pore surface due to the larger the specific surface area in micropores. In addition, the gas flow is slow at low pressure and the differences of stress are low. Thus, the matrix expansion caused by the effective stress should be weakened and matrix shrinkage accounts for the more changes of pore size under low reservoir pressure. As shown, the ratio increases gradually along with reservoir pressure as the pore is meso-and macro-pore (above 5 nm). The reason is the fact that there is less adsorbed gas and higher gas flow velocity in nanopores as the pore size enlarges. Moreover, the particle connection of matrix is weakened and more susceptible to effective stress, resulting in a significant decrease in pore size. The results in Figure 17 suggest that the effect of stress sensitivity and matrix shrinkage should be carefully included in the shale permeability evaluation models. More accurate knowledge of the factors affecting the pore size and permeability may help in the optimization of shale gas production. Figure 16. The ratio of apparent permeability to Dracy permeability (k qpp : k d ) as a function of the pore radius under different models (modified from Klinkenberg [24] and Beskok and Karniadakis [15]). Reservoir pressure = 5 MPa, temperature = 350 K.  [24] and Beskok and Karniadakis [15]). Reservoir pressure = 5 MPa, temperature = 350 K.

Conclusions
Gas flow mechanisms and permeability evolution are particularly important to evaluate the capacity of shale gas production. In this paper, gas multiple mechanisms in shale nanopores were studied. Specifically, an apparent permeability model and dynamic permeability equations were developed and incorporated into a numerical model by the COMSOL solver. According to theoretical research and numerical analysis, the following conclusions were obtained.
Gas complex flow mechanisms and the collisions between gas molecules and pore walls increase the apparent permeability. As calculated, the ratio of the apparent permeability to the Darcy permeability presents 1-2 orders of magnitude change and explains the unusual gas production from shale reservoirs. The evaluation of apparent permeability shows more sensitivity to small pore sizes (1-10 nm) and low reservoir pressures (below 5 MPa). The surface diffusion, slip flow, and Knudsen flow exists in shale nanopores, and their weights vary widely when the pore size ranges from 1 to 5 nm and the pressure is below 5 MPa. This is why the shale apparent permeability is affected by geological factors and multiple flow mechanisms during gas production. In addition, the results infer that gas flow mechanisms should be characterized by separate equations in the apparent model.
Compared to the apparent permeability, the dynamic permeability declines monotonically by about 1 order of magnitude due to the pore size reduction caused by the combined effects of stress sensitivity and matrix shrinkage. This verifies that shale permeability gradually decreases over a long period of production. Moreover, matrix shrinkage and stress sensitivity act opposite directions in shale permeability, and the degree of action is associated with pore size and reservoir pressure. Thus, the comparison result suggests that permeability should be carefully adjusted along with the production of shale reservoirs. Grateful appreciation is expressed for these supports.

Conflicts of Interest:
The authors declare no conflict of interest. Figure 17. The ratio of dynamic permeability to apparent permeability (k : k qpp ) as a function of the reservoir pressures and the pore size considering stress sensitivity and matrix shrinkage. Temperature = 350 K.

Conclusions
Gas flow mechanisms and permeability evolution are particularly important to evaluate the capacity of shale gas production. In this paper, gas multiple mechanisms in shale nanopores were studied. Specifically, an apparent permeability model and dynamic permeability equations were developed and incorporated into a numerical model by the COMSOL solver. According to theoretical research and numerical analysis, the following conclusions were obtained.
Gas complex flow mechanisms and the collisions between gas molecules and pore walls increase the apparent permeability. As calculated, the ratio of the apparent permeability to the Darcy permeability presents 1-2 orders of magnitude change and explains the unusual gas production from shale reservoirs. The evaluation of apparent permeability shows more sensitivity to small pore sizes (1-10 nm) and low reservoir pressures (below 5 MPa). The surface diffusion, slip flow, and Knudsen flow exists in shale nanopores, and their weights vary widely when the pore size ranges from 1 to 5 nm and the pressure is below 5 MPa. This is why the shale apparent permeability is affected by geological factors and multiple flow mechanisms during gas production. In addition, the results infer that gas flow mechanisms should be characterized by separate equations in the apparent model.
Compared to the apparent permeability, the dynamic permeability declines monotonically by about 1 order of magnitude due to the pore size reduction caused by the combined effects of stress sensitivity and matrix shrinkage. This verifies that shale permeability gradually decreases over a long period of production. Moreover, matrix shrinkage and stress sensitivity act opposite directions in shale permeability, and the degree of action is associated with pore size and reservoir pressure. Thus, the comparison result suggests that permeability should be carefully adjusted along with the production of shale reservoirs. Grateful appreciation is expressed for these supports.