SPE Hanyi Wang

SPE, Schlumberger Matteo Marongiu-Porcu

Impact of Shale-Gas Apparent Permeability on Production Combined Effects of Non-Darcy Flow/Gas Slippage, Desorption, and Geomechanics*

Many shale gas and ultralow permeability tight gas reservoirs can have matrix permeability values in the range of tens to hundreds of nanodarcies. The ultrafine pore structure of these rocks can cause violation of the basic assumptions behind Darcy's law. Depending on a combination of pressure-temperature (P-T) conditions, pore structure and gas properties, Non-Darcy flow mechanisms such as Knudsen diffusion and/or Gas-Slippage effects will impact the matrix apparent permeability. Even though numerous theoretical and empirical models have been proposed to describe the increasing apparent permeability due to Non-Darcy flow/Gas-Slippage behavior in nano-pore space, few literature have investigated the impact of formation compaction and the release of adsorption layer on shale matrix apparent permeability, which can be coupled together with the Non-Darcy flow/Gas-Slippage mechanism, during reservoir depletion.

In this article, we first presented a thoroughly review on gas flow in shale nano-pore space and discussed what factors can impact shale matrix apparent permeability, besides the well-studied Non-Darcy flow/Gas-Slippage behavior. Then, a unified shale matrix apparent permeability model is proposed to bridge the effects of Non-Darcy flow/Gas-Slippage, geomechanics (formation compaction) and the release of adsorption layer into a single, coherent equation. In addition, a mathematical framework for an unconventional reservoir simulator that developed in this study is also presented.

Different matrix apparent permeability models are implemented in the numerical simulator to examine what factors impact the matrix apparent permeability most. Furthermore, how these factors influence the average apparent permeability in the Simulated Reservoir Volume (SRV) and long-term production are also discussed. Finally, the impact of natural fracture networks on matrix apparent permeability evolution is investigated. The results indicate that even though the conductive fracture networks play a vital role in shale gas production, the matrix apparent permeability evolution during pressure depletion can still make a difference.

Introduction

In classic fluid flow mechanics where continuum theory holds, fluid velocity is assumed to be zero at the pore wall [START_REF] Sherman | The Transition From Continuum to Molecular Flow[END_REF]. This is a valid assumption for conventional reservoirs having pore radii in the range of 1 to 100 micrometer, because fluids flow as a continuous medium. Correspondingly, Darcy's equation, which models pressure-driven viscous flow, works properly for such reservoirs. However, in shale reservoirs, which have pore throat radii in the range of 1 to 200 nanometers, fluid continuum theory breaks down and gas molecules follow a somewhat random path while still maintaining a general flow direction governed by pressure gradient. Molecules strike against the pore walls and tend to slip at pore walls instead of having zero velocity [START_REF] Sherman | The Transition From Continuum to Molecular Flow[END_REF]. Although these phenomena can be modeled accurately using molecular physics, this is not practical for modeling flow of natural gas through shale reservoirs, which requires intensive computation that falls beyond current computational capabilities. It is thus desirable to integrate molecular flow behavior with correlations based on the standard Darcy's equation, in order to capture these mechanisms and maintain a relatively light load of computations. [START_REF] Klinkenberg | The Permeability of Porous Media to Liquid and Gases[END_REF] was the first contributor to identify gas slippage in porous media by studying rarified gas flowing at various pressures. He observed that the actual gas flow rate was larger than the expected that predicted with Darcy's equation. He attributed this higher flow to gas slippage along the pore walls and proposed the apparent gas permeability adjusted by a slippage factor.

In order to propose corrections for Non-Darcy flow over different flow regions in nanopore space, numerous authors [START_REF] Civan | Effective Correlation of Apparent Gas Permeability in Tight Porous Media[END_REF][START_REF] Darabi | Gas flow in ultra-tight shale strata[END_REF][START_REF] Florence | Improved Permeability Prediction Relations for Low-Permeability Sands[END_REF][START_REF] Fathi | Shale Gas Correction to Klinkenberg Slip Theory[END_REF][START_REF] Javadpour | Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone)[END_REF][START_REF] Javadpour | Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone)[END_REF][START_REF] Sakhaee-Pour | Gas Permeability of Shale[END_REF][START_REF] Singh | Non-empirical apparent permeability of shale[END_REF] have quantified these effects by modifying the slippage factor or determining the apparent permeability as a function of Knudsen number. [START_REF] Swami | Non Darcy Flow in Shale Nanopores: Do We Have a Final Answer[END_REF] investigated and compared 10 theoretical and empirical models for quantifying Non-Darcy flow/Gas-Slippage effects in unconventional gas reservoirs; although these different approaches provided different results, they all indicated that the flow rate at nanoscale level is greater than that predicted from Darcy's equation. Specifically, compared to other theoretical models investigated in the study, the correlations proposed by [START_REF] Sakhaee-Pour | Gas Permeability of Shale[END_REF] are considered to be most reliable because they are supported by well-designed laboratory experiment.

Even though numerous correlations have been developed for apparent permeability enhancement in nano-pore structures, to account for the effects of Non-Darcy flow/Gas-Slippage behavior, they ignore the changes in pore-structure geometry induced by stress variations and the release of attached gas molecules during production. The assumption of constant pore radius is valid under steady-state flow, however, within these ranges of ultralow permeability, transient and pseudo steady-state flow dominate the entire production life of shale formaitons. With constant decreasing pore pressure and increasing effective stress, the formation microstructure will be affected, which may lead to decrease in matrix permeability. In fact, many studies [START_REF] Pedrosa | Pressure Transient Response in Stress-Sensitive Formation[END_REF][START_REF] Soeder | Porosity and Permeability of Eastern Devonian Gas Shale[END_REF]Vera and Ehlig-Economides, 2013;[START_REF] Xiao | Dynamic Characteristic Evaluation Methods of Stress Sensitive Abnormal High Pressure Gas Reservoir[END_REF] revealed that shale matrix permeability can also be stress dependent. Specifically, the increase in the effective stress as a result of pore pressure depletion during production often results in a loss of productivity associated to shale matrix permeability impairment, especially for overpressured shale gas reservoirs that exhibit stress-sensitive permeability characteristics. Laboratory experiment [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF] on shale samples even indicates that shale porosity and permeability decrease with increasing effective stress, following a power law relationship. However, whether the decreasing permeability is induced by the shrinkage of nano-pore space or the closure of micro-fractures is not well studied.

It is very important to realize that, while all these studies focusing on Non-Darcy flow/Gas-Slippage behavior in porous media do not account for nano-pore structure changes during production, the other set of studies that focus on stress-dependent permeability do not consider nano-flow mechanism and adsorption layer thickness variation in the shale matrix. Because these two sets of studies lead to opposite conclusions on shale matrix permeability evolution as pore pressure decreases and effective stress increases, neglecting either aspect may be an unacceptable oversimplification and result in unrealistic permeability prediction. In order to address this paradox and bridge the contradictory impacts of nano-flow mechanisms and nano-pore structure alteration on shale matrix permeability during production, a unified permeability model is thus particularly suited.

In this study, we propose a unified matrix apparent permeability model for shale gas formations, which unifies the mechanisms of Non-Darcy flow/Gas-Slippage, gas adsorption layer and geomechanical effects into a coherent global workflow, as shown in Fig. 1. The impact of different permeability models upon long-term production is assessed and compared by developing an unconventional shale gas reservoir simulator, which involves fluid flow within the shale formation matrix and fractures, gas adsorption/desorption as a function of pressure and real gas properties. The results of this study provide us with a univocal understanding of the net outcomes derived from the complex mechanisms that dominate shale matrix permeability and offer a solid foundation for further investigations and analyses on shale gas production forecasting and reservoir modeling. Transport properties of the shale matrix control the transport of (desorbed) gas toward the hydraulic fractures. Therefore, matrix flow can directly impact the long-term productivity of fractured shale gas reservoirs, even with the presence of natural fractures, fissures, cleats, etc.

Fig. 1-Mechanisms that alter shale matrix apparent permeability during production.

The structure of this paper is as follows. First, Non-Darcy flow/Gas-Slippage mechanisms are reviewed and the effects of geomechanics and adsorption layer are incorporated into the basic formula for matrix apparent permeability. Next, the general mathematical framework for numerical simulation and reservoir modeling (a unit SRV pattern) are presented. Then, the simulation results are analyzed and compared between different permeability models. Finally, conclusion remarks and discussions are presented.

Formulation of Unified Apparent Permeability Model in Shale Matrix

Basic Formula for Non-Darcy Flow/Gas-Slippage Effects Darcy's law cannot describe the actual gas behavior and transport phenomena in nano-porous media. In such nano-pore structure, fluid flow departs from the well-understood continuum regime, in favor of other mechanisms such as slip, transition, and free molecular conditions. The Knudsen number is a dimensionless parameter that can be used to differentiate flow regimes in conduits at micro and nanoscale, and it is defined as the ratio of the molecular mean free path length, 𝜆, to characteristic length of the channel, which in the case of shale matrix is the effective pore radius, 𝑟 [START_REF] Knudsen | Die Gesetze der Molekularstromung und der inneren Reibungsstromung der Gase durch Roḧren (The laws of molecular and viscous flow of gases through tubes)[END_REF]:

𝐾 𝑛 = 𝜆 𝑟 (1)
The mean free path can be calculated by [START_REF] Civan | Shale-Gas Permeability and Diffusivity Inferred by Improved Formulation of Relevant Retention and Transport Mechanisms[END_REF]:

𝜆 = 𝜇 𝑔 𝑃 √ 𝜋R𝑇 2M (2)
where 𝜇 𝑔 is the gas viscosity, T is the reservoir temperature, 𝑃 is the reservoir pressure, M is the gas average molecular weight and R is the universal gas constant. Substituting the mean free path by including real gas, Z-factor yields the following relationship:

𝐾 𝑛 = 𝜇 𝑔 𝑍 𝑃𝑟 √ 𝜋𝑅𝑇 2𝑀 (3)
Table 1 shows how these different flow regimes, which correspond to specific flow mechanisms, can be classified within different ranges of 𝐾 𝑛 . The assumptions underlying Darcy's Law, assuming that momentum transfer by means of bulk phase viscosity and fluid velocity matches solid velocity at walls, is only valid within continuum regime. However, the Knudsen number 𝐾 𝑛 in most shale reservoirs lies between 10 -3 to 1 [START_REF] Ziarani | Knudsen's Permeability Correction for Tight Porous Media[END_REF], which indicates that the slip and transition flow regimes are most likely to be encountered in most shale gas reservoirs. The apparent permeability of shale matrix can be represented by the following general form:

𝑘 𝑎 = 𝑘 ∞ 𝑓(𝐾 𝑛 ) (4)
where 𝑘 ∞ is the intrinsic permeability of the porous medium, which is defined as the permeability for a viscous, non-reacting ideal liquid, and 𝑓(𝐾 𝑛 ) is the correlation term that relates the matrix apparent permeability and intrinsic permeability. Different models [START_REF] Aguilera | Incorporating Capillary Pressure, Pore Throat Aperture Radii, Height above Free Water Table, and Winland r35 Values on Pickett Plots[END_REF][START_REF] Beskok | A Model for Flows in Channels, Pipes, and Ducts at Micro and Nano Scales[END_REF] have been developed, based on experiments, to quantify the relationship between intrinsic permeability and nano-pore structure in porous media. For a capillary tube of radius, 𝑟, the intrinsic permeability can be derived from Eq. ( 23) in [START_REF] Beskok | A Model for Flows in Channels, Pipes, and Ducts at Micro and Nano Scales[END_REF]:

𝑘 ∞ = 𝑟 2 8 (5)
From Eq. ( 5), it can be noticed that the matrix intrinsic permeability 𝑘 ∞ is an intrinsic property of the porous medium itself, and is only related the nano-pore geometry that allows fluid flow effectively.

Sakhaee-Pour and Bryant (2012) developed correlations for 𝑓(𝐾 𝑛 ) based on extensive lab experiments. They proposed a firstorder permeability model in the slip regime and a polynomial form for the permeability enhancement in the transition regime:

𝑓(𝐾 𝑛 ) = { 1 + 5𝐾 𝑛 Slip Regime 0.8453 + 5.4576𝐾 𝑛 + 0.1633𝐾 𝑛 2 Transition Regime 0.1 < 𝐾 𝑛 < 0.8 (6)
Intuitively, Eqs. ( 3) and ( 6) predict net increases of 𝐾 𝑛 and 𝑓(𝐾 𝑛 ) with decreasing pore pressure, respectively. However, decreasing pore pressure can also lead to reduction in pore radii and, in turn, reduce the intrinsic permeability 𝑘 ∞ . Eq. ( 4) shows that the evolution of apparent permeability is determined by two combined mechanisms (i.e., the reduction in 𝑘 ∞ and the increase of 𝑓(𝐾 𝑛 ) during production). Consequently, if geomechanical effects dominate the system ( 𝑘 ∞ is very sensitive to effective stress variation), then 𝑘 𝑎 will decrease over production time. On the other hand, if nano-flow mechanisms dominate the system (𝑓(𝐾 𝑛 ) increases dramatically with decreasing pore pressure), 𝑘 𝑎 will increase over production time.

Formula for Coupled Non-Darcy Flow and Geomechanical Effects

The next step is to incorporate the impact of formation effective stress upon pore structure, intrinsic permeability, and Non-Darcy flow. Laboratory measurements by [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF] shows that the relationship between porosity and stress can be expressed by a power law:

∅ = ∅ 0 ( 𝑃 𝑒 𝑃 0 ) -𝐶 ∅ (7)
where ∅ 0 represents the formation porosity under reference pressure 𝑃 0 , 𝑃 𝑒 is the effective confined pressure and 𝐶 ∅ is a dimensionless material-specific constant that can be determined by lab experiments. For silty-shale samples, the values of 𝐶 ∅ range from 0.014 to 0.056. It should be mentioned that all the presented variables with a subscript 0 corresponds to their value at the reference state, which can be laboratory or initial reservoir conditions. For simplicity, in this paper, the reference state is chosen to be standard conditions under laboratory environment, then 𝑃 0 can be set to be 0.1MPa and Eq. ( 7) can be modified as the following:

∅ = ∅ 0 ( 𝜎 𝑚 𝑃 0 ) -𝐶 ∅ (8)
where 𝜎 𝑚 is the mean in-situ effective stress. In a general porous media, the loss in cross-section area is equivalent to the loss in porosity,

∅ ∅ 0 = 𝑟 2 𝑟 0 2 (9)
Combining Eq. ( 8) and Eq. ( 9) leads to

r = 𝑟 0 ( 𝜎 𝑚 𝑃 0 ) -0.5𝐶 ∅ (10)
According to Eq. ( 5), we can have

𝑘 ∞ 𝑘 ∞0 = 𝑟 2 𝑟 0 2 (11)
Combining Eq. ( 10) and Eq. ( 11) leads to

𝑘 ∞ = 𝑘 ∞0 ( 𝜎 𝑚 𝑃 0 ) -𝐶 ∅ ( 12
)
where 𝑟 0 and 𝑘 ∞0 are pore radius and matrix intrinsic permeability under reference pressure 𝑃 0 (e.g., 0.1 MPa at standard conditions). Eq. ( 12) finally allows incorporating the effects of geomechanics into nano-flow behavior as depicted in Eq. ( 13)

𝑘 𝑎 = 𝑘 ∞0 ( 𝜎 𝑚 𝑃 0 ) -𝐶 ∅ 𝑓(𝐾 𝑛 ) (13)
Compare to Eq. ( 4), which only describes the effect of Non-Darcy flow/Gas-Slippage on matrix apparent permeability, Eq. ( 13) reflects the coupled effects of both Non-Darcy flow and reservoir compaction on matrix apparent permeability Formula for Coupled Non-Darcy Flow, Geomechanical and Adsorption Layer Effects The Langmuir isotherm (1916), which characterizes monolayer adsorption, is the most-used gas adsorption theory for shale gas reservoir modeling, and the thickness of the gas adsorption layer can be interpolated based on a Langmuir type functional relationship:

δ = 𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 (14)
where 𝑑 𝑚 is the average diameter of gas molecules residing on the pore surface, and 𝑃 𝐿 is the Langmuir pressure. When the effects of the absorption layer are considered, the effective pore radius in Eq. ( 10) should then be modified as:

r = 𝑟 0 ( 𝜎 𝑚 𝑃 0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 (15)
Thus, the corresponding porosity and intrinsic permeability in Eq. ( 9) and ( 12) can be rewritten as:

∅ = ∅ 0 (𝑟 0 ( 𝜎 𝑚 𝑃 0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 ) 2 𝑟 0 2 (16) 𝑘 ∞ = 𝑘 ∞0 (𝑟 0 ( 𝜎 𝑚 𝑃 0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 ) 2 𝑟 0 2 (17)
Thus, combine Eq. ( 4) and Eq. ( 17), the in-situ matrix apparent permeability can be calculated by:

𝑘 𝑎 = 𝑘 ∞0 (𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 ) 2 𝑟 0 2 𝑓(𝐾 𝑛 ) (18)
where 𝑓(𝐾 𝑛 ) represents the correction Non-Darcy flow/Gas-Slippage behavior, and the term in front of it (also reflected in Eq.( 17)) represents the overall correction for matrix intrinsic permeability, which is solely determined by the nano-pore geometry that can be altered by formation compaction and release of adsorption layer.

In summary, Eq. ( 4), ( 13) and ( 18) are three permeability models that will be investigated thoroughly in this study to examine the impact of the coupled effects on shale matrix apparent permeability and long term shale gas production. These variables that do not contain the subscript 0 refer to the values at in-situ conditions. From the derivation of these equations, it can be noticed that some assumptions have been made in the process in order to combine the mechanisms of Non-Darcy flow, geomechanics and gas adsorption layer into a single, coherent equation. However, it should be mentioned that all the assumptions used in the derivations are based on empirical relationships that obtained from laboratory experiment published in literature, and thus, can be deemed as reasonable assumptions.

Even though in this study, we use Sakhaee-Pour and Bryant (2012) laboratory derived correlation (Eq. ( 6)) to calculate 𝑓(𝐾 𝑛 ) for Non-Darcy flow/Gas Slippage correction, other theoretical or empirical models may also be applied to determine 𝑓(𝐾 𝑛 ).

Because the matrix apparent permeability model developed in this section is a general formula, and the Non-Darcy flow/ Gas-Slippage correction factor 𝑓(𝐾 𝑛 ) is independent from the correction of intrinsic permeability (Eq. ( 17)). Compare the simulation results of various Non-Darcy flow/ Gas Slippage models on 𝑓(𝐾 𝑛 ) and the associated coupled effects is out of the scope of this paper, but a comprehensive comparison between most widely used Non-Darcy flow/ Gas Slippage models can be found in Swami et al. ( 2012)

Coupled Real Gas Flow and Induced Stresses in a Porous Medium

As the pore pressure decreases during production, formation porosity, intrinsic permeability, pore radius and Knudsen number are subject to time and space variation. Consequently, a fully coupled unconventional reservoir simulator is particularly needed in order to investigate the impact of these different permeability models on prediction of long-term shale gas production. Previous studies [START_REF] Shabro | Numerical Simulation of Shale-Gas Production: From Pore-Scale Modeling of Slip-Flow, Knudsen Diffusion, and Langmuir Desorption to Reservoir Modeling of Compressible Fluid[END_REF]2012) have developed theoretical models for gas flow in organic shales with the effects of gas diffusion in kerogen and advection in nano-pores, but how to quantify these effects in complex shale systems are still not well understood within scientific and industry community, and the validation of these theoretical models by tailored laboratory experiment will be needed for practical reservoir simulation. In this paper, apparent permeability model derived from previous section is used to model the gas flow rate within the shale matrix:

𝒒 𝑔 = - 𝐤 a 𝜇 𝑔 • ∇𝑃 (19)
where 𝒒 𝑔 is the gas velocity vector, 𝐤 𝐚 is the gas matrix apparent permeability vector. The continuity equation within shale gas formation can be written as:

𝜕𝑚 𝜕𝑡 + ∇ • (𝜌 𝑔 𝒒 𝑔 ) = 𝑄 𝑚 ( 20
)
where 𝑚 is the gas mass content per unit volume, 𝜌 𝑔 is gas density, 𝑄 𝑚 is the source term and 𝑡 is the generic time. The gas mass content 𝑚 is obtained from two contributions:

𝑚 = 𝜌 𝑔 ∅ + 𝑚 𝑎 (21)
𝜌 𝑔 ∅ is the free gas mass in the shale pore space per unit volume of formation, while 𝑚 𝑎 is the adsorbed gas mass per unit volume of formation, which can be determined from the Langmuir isotherm [START_REF] Langmuir | The Constitution and Fundamental Properties of Solids and Liquids[END_REF]:

𝑚 𝑎 = 𝜌 𝑚 𝜌 𝑔𝑠𝑡 𝑉 𝐿 𝑃 𝑃 + 𝑃 𝐿 (22)
where 𝜌 𝑚 is the shale matrix density, 𝜌 𝑔𝑠𝑡 is the gas density at standard conditions and 𝑉 𝐿 is the Langmuir volume.

The in-situ gas density is calculated according to the real gas law:

𝜌 𝑔 = 𝑃𝑀 𝑍𝑅𝑇 (23)
The Z-factor can be calculated with an explicit correlation [START_REF] Mahmoud | Development of a New Correlation of Gas Compressibility Factor (Z-Factor) for High Pressure Gas Reservoirs[END_REF] based on the pseudo-reduced pressure (p pr ) and pseudo-reduced temperature (T pr ):

𝑍 = (0.702𝑒 -2.5𝑇 𝑝𝑟 )(𝑝 𝑝𝑟 2 ) -(5.524𝑒 -2.5𝑇 𝑝𝑟 )(𝑃 𝑝𝑟 ) + (0.044𝑇 𝑝𝑟 2 -0.164𝑇 𝑝𝑟 + 1.15) ( 24
)
Gas viscosity is an intrinsic property of gas itself, which can be determined based on its compositions, pressure and temperature. Lee et al. Correlation (1966) is used to estimate the gas viscosity: The porous medium is assumed to be perfectly elastic so that no plastic deformation occurs. The constitutive equation can be expressed in terms of effective stress ( 𝜎 𝑖𝑗 , ), strain ( 𝜀 𝑖𝑗 ), and pore pressure (P):

𝜇 𝑔 =
𝜎 𝑖𝑗 = 2G𝜀 𝑖𝑗 + 2G𝜈 1 -2𝜈 𝜀 𝑘𝑘 𝛿 𝑖,𝑗 -𝛼𝑃𝛿 𝑖,𝑗 (26)
where G is the shear modulus, 𝜈 is the Poisson's ratio, 𝜀 𝑘𝑘 represents the volumetric strain, 𝛿 𝑖,𝑗 is the Kronecker delta defined as 1 for 𝑖 = 𝑗 and 0 for 𝑖 ≠ 𝑗, and 𝛼 is the Biot's effective stress coefficient, which is assumed to be 1 in this study. The straindisplacement relationship and equation of equilibrium are defined as:

𝜀 𝑖𝑗 = 1 2 (𝑢 𝑖,𝑗 + 𝑢 𝑗,𝑖 ) (27) 𝜎 𝑖𝑗,𝑗 + 𝐹 𝑖 = 0 ( 28
)
where 𝑢 𝑖 and 𝐹 𝑖 are the components of displacement and net body force in the i-direction. Combining Eq. ( 26)-( 28), we have a modified Navier equation in terms of displacement under a combination of applied stress and pore pressure variations:

𝐺∇𝑢 𝑖 + 𝐺 1 -2𝜈 𝑢 𝑗,𝑖𝑖 -𝛼𝑃𝛿 𝑖,𝑗 + 𝐹 𝑖 = 0 (29)
Flow along the Hydraulic Fracture Darcy's law is used to describe the flow behavior inside the hydraulic fracture:

𝒒 𝑓 = - 𝒌 𝑓 𝜇 𝑔 • 𝑑 𝑓 ∇ 𝑇 𝑃 (30)
where 𝐪 𝐟 is the gas volumetric flow rate vector per unit length in the fracture, 𝐤 𝐟 is the fracture permeability vector, and ∇ T P is the pressure gradient tangent to the fracture surface. The continuity equation along the fracture reflects the generic material balance within the fracture:

𝑑 𝑓 𝜕∅ 𝑓 𝜌 𝑔 𝜕𝑡 + ∇ 𝑇 • (𝜌 𝑔 𝒒 𝑓 ) = 𝑑 𝑓 𝑄 𝑓 (31)
where 𝑑 𝑓 is the fracture width, ∅ f is the fracture porosity, and Q f is the mass source term, which can be calculated by adding the mass flow rate per unit volume from two fracture walls (left and right):

𝑄 𝑓 = 𝑄 𝑙𝑒𝑓𝑡 𝑓 + 𝑄 𝑟𝑖𝑔ℎ𝑡 𝑓 (32 -𝑎) 𝑄 𝑙𝑒𝑓𝑡 𝑓 = - 𝒌 𝑎 𝜇 𝑔 𝜕𝑃 𝑙𝑒𝑓𝑡 𝜕𝒏 𝒍𝑒𝑓𝑡 (32 -𝑏) 𝑄 𝑟𝑖𝑔ℎ𝑡 𝑓 = - 𝒌 𝑎 𝜇 𝑔 𝜕𝑃 𝑟𝑖𝑔ℎ𝑡 𝜕𝒏 𝒓𝒊𝒈𝒉𝒕 (32 -𝑐)
where 𝒏 is the vector perpendicular to fracture surface. These general equations for modeling flow in hydraulic fracture can be also applied to natural fractures.

Model Description

Based on the principles introduced in the previous section, an unconventional reservoir simulator was constructed and all the derived coupled equations were solved numerically. This simulator makes use of a finite element analysis [START_REF] Zienkiewicz | The Finite Element Method[END_REF] scheme based on the Newton-Raphson method [START_REF] Ypma | Historical development of the Newton-Raphson method[END_REF], which requires evaluating a Jacobian matrix starting from initial guess values in order to approximate the solution throughout successive iterations. All the variables in the system are updated at the end of each time increment and input as initial values at the start of the next time increment. In order to avoid extreme mesh refinement inside and around fracture and guaranteeing solution convergence, a lower-dimension interior boundary, using tangential derivative conditions based on Eqs. ( 30) through ( 32), is implemented to represent the fluid flow behavior along the fracture. The numerical simulator program developed in this study is based on previous numerical framework by [START_REF] Wang | Conceptual Study of Thermal Stimulation in Shale Gas Formations[END_REF]. Even though it is possible to simulate an entire section of a horizontal well with multiple fractures, it is more efficient to simulate a unit SRV pattern (shown in Fig. 2) and apply symmetric conditions along the boundaries.

Fig.2-A plane view of half wing transverse hydraulic fracture and discretized mesh.

Table 2 shows all the input parameters used for the following simulations, which includes laboratory data, reservoir conditions, drainage geometry, formation properties, real gas, and gas adsorption parameters. In order to make a consistent assessment on how the Non-Darcy flow, adsorption layer, and geomechanical effects impact the shale matrix permeability during production, the reference input values for porosity, intrinsic permeability, and pore radius for all models are considered at standard conditions under laboratory environment, while their corresponding values at initial reservoir conditions are calculated based on in-situ reservoir pressure-temperature and stress conditions. The rock material constant (C ∅ ), is set to be 0.035, which is the median value of laboratory measurements [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF]; the diameter of adsorption gas molecules ( d m ) is assumed to be 0.414 nm, which is the diameter of a methane molecule.

The permeability of shale rocks is often measured on crushed core samples and core plugs by a variety of methods [START_REF] Mallon | How should permeability be measured in fine-grained lithologies? Evidence from the chalk[END_REF]. The crushed core permeability measurements yield true matrix values, eliminating unrepresentative contributions of pre-existing microfractures. It is also often used when good quality core plugs cannot be obtained. For the specific computations presented in this study, the permeability measured at laboratory conditions is converted to apparent permeability at initial reservoir conditions at the beginning of the simulation through different permeability models. e.g., when formation compaction is considered, the intrinsic permeability that corresponds to the initial reservoir effective stress is calculated using Eq. ( 12), starting from the 100-nD intrinsic permeability at laboratory conditions; the non-Darcy/gas-slippage correction term 𝑓(𝐾 𝑛 ) corresponding to the initial Knudsen number is calculated using Eq. ( 6). Then, the final apparent permeability at initial reservoir conditions is determined by multiplying these two variables using Eq. ( 13), providing the value of 120 nD for the given initial reservoir conditions.

For the purpose of simplification and easy interpretation of simulation results, we assume that the formation is isotropic, so the permeability tensor vector in the fluid flow equations can be treated as scalar and pore size is a representative average value.

For a more realistic simulation, stochastic permeability model with characterization of laboratory derived pore size distributions (Naraghi and Javadpour, 2015) can be applied.

Results and Analysis

In this section, we present our analysis results and demonstrate how the Non-Darcy flow, adsorption layer, and geomechanical effects will impact the shale matrix permeability during production (including both transient linear flow and pseudo-steady state flow period), and how these mechanisms can influence long-term shale gas production. The hydraulic fracture conductivity is set to remain constant as its initial value over time, in order to obtain a direct assessment of the abovementioned effects on the shale matrix apparent permeability.

First, the base case simulation results are presented and the impact of the absorption layer is assessed. Then, the simulation results for different combinations of permeability variation mechanisms are compared for 30 years of production history. In addition, sensitivity studies for the different permeability variation mechanisms, effective stress and pore radius are examined. Finally, the effects of the presence of conductive fracture networks are investigated.

Fig. 3 shows the pressure distribution and the adsorbed gas per unit volume in the simulated drainage area at the end of 5 years of production history. It can be observed that the pressure disturbance has reached the boundary and pseudo steady-state flow has already ensued from initial transient flow, to dominate the rest of the production life. In turn, the adsorbed gas layer has already begun to desorb as pressure declines locally. It appears that almost 60% of such adsorbed gas has already been produced near the pressure sink along the hydraulic fracture walls, but it still remains untapped within the region beyond the fracture tip. Fig. 4 shows porosity distribution within the simulated SRV. The results demonstrate that the adsorption layer affects considerably the formation porosity during production. The left-hand side of Fig. 4 does not include the gas adsorption layer on the pore walls, and we can observe that the porosity in the regions adjacent to the hydraulic fracture (i.e., at lower pressure) can decrease to as low as 5% of the that at initial reservoir conditions as a result of the pore size reduction due to increased effective stress. However, when the presence of the attached adsorbed gas molecules on the pore surface is considered (righthand side of Fig. 4), a reversed trend is observed. The porosity appears to increase toward the regions adjacent to the hydraulic fracture (i.e., as pore pressure decreases); this implies that the increase in pore cross section due to the desorption of the adsorbed gas is able to offset the adverse effect of matrix compaction, although at this specific level of simulated depletion the overall porosity values are still lower than the ones simulated without the gas adsorption layer. Fig. 5 shows the Knudsen number distribution within the simulated SRV; left-hand side does not include the gas adsorption layer on the pore walls, while the right-hand side does. We can observe that the values of Knudsen number with the gas adsorption layer are slightly larger than the ones without the gas adsorption layer, and this is due to the reduced pore radii caused by the presence of the adsorbed molecules attached on the pore surface. Nevertheless, in both cases, the values of Knudsen number fall within the slip and transition regimes, according to the classification in Table 1. Furthermore, the slip regime dominates the entire simulated drainage area at initial reservoir conditions, but while the pore pressure decreases, the transition regime manifests itself in the regions adjacent to the hydraulic fracture (i.e., at lower pressure) and gradually spreads out along with declining reservoir pressure.

Fig. 5-Knudsen number distribution with and without adsorption layer after 5 years of production.

Fig. 6 shows the shale matrix apparent permeability distribution within the simulated SRV. The left-hand side does not include the gas adsorption layer on the pore walls, while the right-hand side does. For both cases, we can observe that after 5 years of production, the apparent shale matrix permeability in the regions adjacent to the hydraulic fracture (i.e., at lower pressure) is approximately 2 to 3 times greater than the initial reservoir matrix permeability. This indicates that in the more depleted areas Non-Darcy flow dominates the apparent permeability evolution mechanisms. However, we can also observe the presence of a dark contour that defines an "envelope" in which the shale matrix permeability is actually lower than the ones at initial conditions. The extension of this permeability envelope is larger when considering the gas adsorption layer on the pore walls. This phenomenon is due to the fact that formation compaction can temporarily dominate the matrix apparent permeability evolution, where the reduction in intrinsic permeability due to shrinkage of pore size cannot be compensated by Non-Darcy flow/Gas-Slippage behavior. However, as pressure continues to decline, the impact of formation compaction will be alleviated and the matrix apparent permeability will start to increase.

Impact of Different Permeability Models

Next, we compare and investigate three different matrix apparent permeability models (i.e., Non-Darcy flow; Non-Darcy flow with geomechanical effects; Non-Darcy flow with both geomechanical and adsorption layer effects) and their composite impact on the average matrix apparent permeability in the SRV during long term production history. The average matrix apparent permeability is calculated by integrating the simulated 𝑘 𝑎 profiles over the entire domain and divided by its volume.

The input parameters are the same for all different models, as depicted in Table 2; the only difference is the shale matrix apparent permeability is calculated with different equations (i.e., Eq. ( 4), ( 13) and ( 18)).

Fig. 7 shows the evolution of the average matrix apparent permeability over time; notice that the value at initial reservoir conditions for the Non-Darcy flow curve and the Non-Darcy flow with geomechanical effects curve is 120 nD, which is estimated (using Eq.( 4) and ( 13)) from the reference value at laboratory conditions from Table 2. It can also be noted that the initial apparent permeability is lower than 120 nD when the adsorption layer is included due to the reduced intrinsic permeability from the shrinkage of effective pore radius (using Eq.( 18)), and the whole curve is shifted downward, but the general increasing trend is preserved. This result is extremely significant, because it provides a robust theoretical validation to the acclaimed experimental work presented in Sakhaee-Pour and Bryant (2012) when the adsorption layer is considered. Our simulated average matrix apparent permeability in the SRV increases during production due to the dominant effects of Non-Darcy flow and the gradual release of the adsorption layer from the pore surface, over the negative geomechanical effect (formation compaction). It can be noted that the three curves in Fig. 7 exhibit a consistent discontinuity (i.e., around 5 to 7 years of simulated production) due to the relief of the local "dark envelope", where the matrix permeability is lower than its initial permeability, as seen in Fig. 6.

Fig. 7-Average matrix apparent permeability in SRV with different permeability models.

Even though Fig. 7 shows three distinct permeability trends for the three permeability models considered, the long-term simulated production performance (shown in Fig. 8) shows that the geomechanical effects (formation compaction) have negligible impact on the long-term production. This behavior is interpreted as an offset between the negative compaction effect on the shale matrix intrinsic permeability and the positive effects of matrix porosity reduction (more free gas expelled from pore space). Nevertheless, and more importantly, when the presence of the adsorption layer on the matrix pore surface is considered, the long-term production performance is visibly impaired, consistent with the matrix apparent permeability downward shift in Fig. 7. It should be mentioned that, even though different apparent permeability models are implemented throughout this study, all the simulations are based on the same input parameters that provided in Table 2, unless otherwise specified. For example, the switch on and off of the adsorption layer effects only impact the matrix apparent permeability through Eq. ( 18), and it does not impact the initial adsorption gas in-place, which is determined by the same input parameter for the Langmuir Isotherm (based on Eq.( 22)) in all simulations. For clarification, an additional curve (i.e., blue dash line) that only considers free gas in-place with initial reservoir permeability is shown in Fig. 8. It can be observed that the cumulative production is much lower if adsorption gas is not included in the simulation model for material balance calculation. Next, we increase the rock material constant, C ∅ , from 0.035 to 0.056, which is the highest value reported by the laboratory measurements in shale samples [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF]. The higher the value of C ∅ , the more sensitive the porosity will be to the effective stress, such that the geomechanical effects can have larger influence on shale matrix pore size and intrinsic permeability.

Fig. 9 shows the porosity and matrix apparent permeability distribution within the simulated SRV. Similarly to Fig. 6, we can observe the presence of a dark contour that defines an "envelope" in which the porosity and the shale matrix permeability are actually lower than the ones at initial reservoir conditions. This is because of the shrinkage of pore size cannot be compensated by the release of gas adsorption layer and Non-Darcy flow behavior in some regions. However, in the more depleted region inside the envelope (i.e., closer to the fracture wall) the actual porosity is greater than the one at initial reservoir conditions, due to the beneficial effect of gas desorption from the pore surface. Interestingly, under the considered increased porosity sensitivity to the effective stress, the matrix apparent permeability will still increase as much as three times with respect to the one at initial reservoir conditions within the low-pressure region adjacent to the pressure sink along the hydraulic fracture. Similarly to Fig. 7, Fig. 10 shows the evolution of the average matrix apparent permeability over time within the simulated SRV, for different values of C ∅ . The illustrated results indicate that even for more stress-sensitive shale formations, the Non-Darcy flow mechanisms still dominate the overall increasing trend of matrix apparent permeability, regardless the presence of the gas adsorption layer. However, for more stress-sensitive shale formations, the average matrix apparent permeability curves appear to be shifted down because of the larger reduction in porosity and pore radius when converted from laboratory values to initial reservoir conditions.

Impact of Pore Radius

While all the simulated results presented so far have been obtained considering an average reference pore radius, 𝑟 0 , of 5 nm (as reported in Table 2), we now study the effect of smaller and larger 𝑟 0 (i.e., 3 nm and 8 nm, respectively). We also use the ratio of average matrix apparent permeability to initial reservoir permeability (𝑘 𝑎 /𝑘 𝑎𝑖 ) as an indicator to assess the impact of these different radii for the three different matrix apparent permeability models (i. geomechanical effects; Non-Darcy flow with both geomechanical and adsorption layer effects) and their composite impact on the SRV-average matrix apparent permeability during our simulated long-term production history.

Fig. 11 shows that when the 3-nm pore radius is considered, the average matrix apparent permeability in the SRV increases quite significantly with time, while the 8-nm pore radius offers a minimal increase. This happens because larger values of Knudsen number, 𝐾 𝑛 , can be achieved under low-pressure regimes for smaller-pore radii, indicating that Non-Darcy flow/ Gas-Slippage phenomena in nano-pores are more severe for the smaller-pore radii. Fig. 11 also shows that when the 3-nm pore radius is considered, the gas adsorption layer has a much more significant impact on the matrix apparent permeability evolution than for the case at 8-nm pore radius. This happens because the gas adsorption layer accounts for a larger percentage of the pore throat cross section area for smaller pore radii. Thus, the gradual desorption of the gas molecules form pore walls will lead to a more pronounced net increase in the matrix apparent permeability.

Fig. 12 demonstrates once more that the geomechanical effects provide a very limited impact on the net increase trend of the matrix apparent permeability over time, for all pore radii considered. Even though the pore radius is likely to shrink to some degree during production due to formation compaction, this reduction will intensify the Non-Darcy flow/Gas-Slippage effects, which will offset the loss of intrinsic matrix permeability. This confirms our previous observations that Non-Darcy flow and gas adsorption layer impact the matrix apparent permeability evolution during production in a much more severe way than formation compaction.

Impact of Natural Fractures

It is a well-known fact that most shale formations have massive pre-existing natural fracture networks that are generally sealed by precipitated materials weakly bonded with mineralization. Such poorly sealed natural fractures are generally reported to interact heavily with the hydraulic fractures during the injection treatments, serving as preferential paths for the growth of complex fracture network. The reactivation of such pre-existing planes of weakness (i.e., natural fractures, micro-fractures, fissures) is well documented and observed in microseismic monitoring [START_REF] Cipolla | A Practical Guide to Interpreting Microseismic Measurements[END_REF][START_REF] Peyret | Subsurface To Surface Microseismic Monitoring for Hydraulic Fracturing[END_REF][START_REF] Zakhour | Real-Time Use of Microseismic Monitoring for Horizontal Completion Optimization Across a Major Fault in the Eagle Ford Formation[END_REF].

Because of the low viscosity fracturing fluid generally used for shale gas fracturing (i.e., slickwater) and the narrow, staggered natural fracture paths, the proppant material cannot reach and fill the created fracture network during fracturing stimulation, but the slippage/shear is still considered the main contribution to the increased fracture network conductivity within the SRV. This final section presents our investigation on how the matrix apparent permeability evolution affects the long-term production performance within this complex, realistic scenario. A discrete fracture network (DFN) is randomly generated with 0.1 mm fracture width and 500-md fracture permeability. First, we simulate an ideal case in which fracture network conductivity doesn't change with time. In a second, more realistic case, the expected reduction in permeability (pressuredependent permeability, PDP) of the simulated fracture network caused by the increasing effective stress during production is simulated as [START_REF] Cho | Pressure-Dependent Natural-Fracture Permeability in Shale and Its Effect on Shale-Gas Well Production[END_REF]Raghavan et al., 2004):

𝑘 𝑛𝑓 = 𝑘 𝑛𝑓,𝑖 𝑒 -𝑏𝜎 𝑚 (32)
where 𝑘 𝑛𝑓 is the natural fractures permeability, 𝑘 𝑛𝑓,𝑖 is the natural fracture permeability at initial reservoir conditions, and b is a parameter that can be determined from experimental data for specific reservoir rocks.

Fig. 13 and Fig. 14 show the distribution of pressure and matrix apparent permeability after 1 year and 10 years of production, respectively. It can be clearly observed that with the presence of the conductive fracture network, the reservoir pressure depletes more rapidly (comparison with Fig. 3). The matrix apparent permeability manifests the same behaviors portrayed in Fig. 6 and Fig. 9, at the expected, much-faster pace.

Fig. 15 shows the simulated cumulative production for the three matrix permeability composite models under the presence of the complex fracture network (constant permeability and PDP). Obviously, the highest simulated cumulative production occurs for the constant fracture network conductivity scenario, while the lowest simulated cumulative production occurs when the fracture network is not considered. The more realistic scenario with PDP fracture network falls in between the two theoretical scenarios, for any of the three matrix apparent permeability models.

These results are not intended to replace or diminish the extraordinary current research effort focusing on "unconventional fracture modeling" and unconventional production unstructured grid simulators for non-planar complex fracture networks, but rather provide a novel perspective on the almost-always neglected importance of the matrix apparent permeability evolution on long-term production. The matrix apparent permeability evolution in the nano-pore domain can-and will-indeed make a difference in the long-term production, even within the dominating presence of complex fracture networks. The results also reveal that if the overall effective permeability (consider both fracture network and matrix conductivity) in the SRV region decreases during production, the cause may stem from the reduction in fracture conductivity (i.e., due to proppant crash or embedment), rather than the decrease in matrix permeability.

Conclusions and Discussions

Permeability is one of the most fundamental reservoir rock properties required for modeling hydrocarbon production. Even though the increased apparent permeability in shale nano-pore space due to Non-Darcy flow/Gas-Slippage has been extensively studied in recent years, the impact of increasing effective stress and release of adsorption layer on the nano-pore geometry and the associated intrinsic permeability has been largely ignored.

In this article, we first presented a comprehensive literature review on gas flow behavior in shale nano-pore space, then a unified matrix apparent permeability model is proposed, which not only accounts for the Non-Darcy flow/Gas-Slippage correction, but also introduces the correction for intrinsic permeability that capture the changes in nano-pore structure due to formation compaction and gas adsorption layer. Next, we presented a mathematical framework for developing a shale gas reservoir simulator and three matrix apparent permeability models are implemented in the following simulations. Finally, factors that influence the average matrix apparent permeability in the SRV, as well as long-term shale gas production are investigated. The most relevant conclusions emerging from this research work are:

1. The matrix apparent permeability in shale formation is not only determined by the Non-Darcy flow/Gas-Slippage behavior, but also by the intrinsic permeability within the nano-pore structure. The intrinsic permeability can be altered by the competitive mechanisms of the increasing effective stress and the release of the adsorbed molecules during production, while the Non-Darcy flow/Gas-Slippage effects are mostly controlled by the pressure depletion and effective pore radius.

2. The shale matrix permeability measured at laboratory conditions needs to be correlated to reservoir conditions cautiously, because the Non-Darcy flow/Gas-Slippage behavior, adsorption layer, and geomechanical effects are extremely sensitive to the pore radius and can cause a severe alteration of the in-situ matrix apparent permeability and simulation results.

3. The geomechanical effects (formation compaction) can temporarily dominate the matrix apparent permeability evolution during the early production stages (where the reduction in intrinsic permeability due to the shrinkage of pore size cannot be compensated by Non-Darcy flow/Gas-Slippage effects). However, as the formation pressure keeps declining, the compaction effects will be offset by Non-Darcy flow/Gas-Slippage behavior and gas desorption effects, and finally the matrix apparent permeability will benefit from pressure depletion and start to increase.

4. The most likely cause for the impaired productivity in fractured shale formations during depletion is the reduction in fracture conductivity, rather than the reduction in matrix apparent permeability.

5. The gradual release of the molecular adsorption layer (desorption) has a significant impact on the matrix apparent permeability evolution when the pore radii are small, while this impact is limited at large pore radii.

6. The matrix apparent permeability evolution during production can make a difference in well performance and long term shale gas production-even when a conductive natural fracture network is present.

Besides the proposed unified shale matrix permeability model, the fully coupled unconventional reservoir simulator developed in this study provides a very comprehensive tool to investigate how the combined real gas nano-flow mechanisms and rock deformations impact the matrix apparent permeability evolution, as well as the long-term productivity of hydraulically fractured shale gas formations. Even though geomechanical effects on the evolution of shale matrix apparent permeability during production is negligible, but it has substantial impact on well performance and long-term gas production, due to the sensitivity of fracture conductivity to in-situ effective stress. Because the evolution of both matrix apparent permeability and fracture conductivity in the SRV region during pressure depletion can inevitably impact well productivity, both comprehensive reservoir characterization and good understanding of propped/unpropped fracture conductivity under various effective stresses are critical for effective reservoir simulation and shale gas production prediction.

Nomenclature

Fig. 3 -

 3 Fig.3-Pressure distribution and adsorbed gas per unit volume after 5 years of production.

Fig. 4 -

 4 Fig.4-Porosity distribution with and without adsorption layer after 5 years of production.

Fig. 6 -

 6 Fig. 6-Matrix apparent permeability distribution with and without adsorption layer after 5 years of production.

Fig. 8 -

 8 Fig. 8-Cumulative production prediction with different permeability models.Impact of Porosity Sensitivity to Stress

Fig. 9 -

 9 Fig.9-Porosity and matrix apparent permeability distribution after 5 years of production with 𝐂 ∅ =0.056.

Fig. 10 -

 10 Fig. 10-Average matrix apparent permeability in SRV with different C ∅ .

  e., Non-Darcy flow; Non-Darcy flow with Flow, Geomechanics, Adsorption Layer,C∅=0.052 Non-Darcy Flow, Geomechanics,C∅=0.036 Non-Darcy Flow, Geomechanics, Adsorption Layer,C∅=0.036

Fig. 11 -

 11 Fig. 11-Impact of Non-Darcy flow, geomechanical effects and adsorption layer with different pore radii.

Fig. 13 -Fig. 14 -

 1314 Fig. 13-Pressure and matrix apparent permeability distribution after 1 year of production.

Table 1 . Fluid Flow Regimes Defined by Ranges of 𝑲 𝒏 (Roy et al., 2003)

 1

	𝐾 𝑛	0 -10 -3	10 -3 -10 -1	10 -1 -10 1	> 10 1
	Flow Regime Continuum	Slip	Transition	Free Molecular

Table 2 . Input Parameters for Base Case Simulation Input Parameters Value Initial

 2 reservoir pressure, 𝑝 𝑖 34.447[MPa] Wellbore pressure, 𝑝 𝑤𝑓 3.447[MPa] Intrinsic permeability measured at laboratory conditions, 𝑘 ∞0 100[nD] Hydraulic fracture permeability, 𝑘 𝑓 10 4 [mD] Formation porosity measured at laboratory conditions, ∅ 0 0.08 Pore radius measured at laboratory conditions, 𝑟 0 5[nm] Material constant for porosity change, 𝐶 ∅

	0.035

SRV 25 Reservoir temperature, T 366[K] Density of formation rock, 𝜌 𝑚 2600 [kg /𝑚 3 ] Langmuir volume constant, 𝑉 𝐿 0.0113[m 3 /kg] Langmuir pressure constant, 𝑃 𝐿 10.34[MPa]

  𝑏 = Parameter for fracture pressure dependent permeability, L𝑡 2 / m, 1/𝑃𝑎 𝐶 ∅ = Material constant 𝑑 𝑓 = Fracture width, L, 𝑚 𝑑 𝑚 = Diameter of absorbed gas molecules, L, 𝑚 𝐹 = Net body force, m/L𝑡 2 , 𝑃𝑎 𝐺 = Shear Modulus, m/L𝑡 2 , 𝑃𝑎 𝑘 ∞ = Intrinsic permeability, 𝐿 2 , 𝑚𝐷 𝑘 ∞0 = Intrinsic permeability at reference condition, 𝐿 2 , 𝑚𝐷 𝑘 𝑎 = Matrix apparent permeability, 𝐿 2 , 𝑚𝐷 𝒌 𝒂 = Matrix apparent permeability vector,𝐿 2 , 𝑚𝐷 𝑘 𝑎𝑖 = Matrix apparent permeability at initial reservoir condition,𝐿 2 , 𝑚𝐷 𝑘 𝑛𝑓,𝑖 = Natural fracture permeability at initial reservoir condition, 𝐿 2 , 𝑚𝐷 𝑘 𝑛𝑓 = Natural fracture permeability, 𝐿 2 , 𝑚𝐷 𝑘 𝑓 = Hydraulic fracture permeability, 𝐿 2 , 𝑚𝐷 𝒌 𝒇 = Hydraulic fracture permeability vector, 𝐿 2 , 𝑚𝐷 𝐾 𝑛 = Knudsen number 𝑚 = Total gas content, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝑚 𝑎 = Gas adsorption mass per unit volume, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝑀 = Molecular weight, m/n, 𝑘𝑔/mol 𝑃 𝑝𝑟 = Pseudo-reduced pressure 𝑃 = Reservoir Pressure, m/L𝑡 2 , 𝑃𝑎 𝑃 0 = Pressure at reference condition, m/L𝑡 2 , 𝑃𝑎 𝑃 𝑒 = Effective confined pressure, m/L𝑡 2 , 𝑃𝑎 𝑃 𝐿 = Langmuir pressure, m/L𝑡 2 , 𝑃𝑎 𝒒 𝒇 = Flow rate vector in the fracture per unit length, 𝐿 2 /t, 𝑚 3 /𝑠/𝑚 𝒒 𝒈 = Velocity vector of gas phase, L/t, 𝑚/s 𝑄 𝑓 = Mass source term in fracture, m/𝐿 3 /𝑡, 𝑘𝑔/𝑚 3 /𝑠 𝑄 𝑚 = Mass source term in matrix, m/𝐿 3 /𝑡, 𝑘𝑔/𝑚 3 /𝑠 𝑟 = Effective pore radius, L, 𝑚 𝑟 0 = Effective pore radius at reference condition, L, 𝑚 𝑅 = Universal gas constant, m𝐿 2 /𝑡 2 /𝑛/𝑇, 8.3145 𝐽/𝑚𝑜𝑙/𝐾 𝑇 = Reservoir temperature, T, 𝐾 𝑇 𝑝𝑟 = Pseudo-reduced temperature 𝑢 𝑖,𝑗 = Component of displacement, L, m 𝑉 𝐿 = Langmuir volume, 𝐿 3 /m, 𝑚 3 /𝑘𝑔 𝑍 = Gas deviation factor 𝛼 = Biot's coefficient 𝛿 = Thickness of gas adsorption layer, L, 𝑚 𝛿 𝑖,𝑗 = Kronecker delta 𝜀 𝑖𝑗 = Elastic strain 𝜀 𝑘𝑘 = Volumetric strain 𝜆 = Characteristic length of flow path, L, 𝑚 𝜇 𝑔 = Gas viscosity, m/Lt, Pa• s 𝜈 = Poisson's ratio 𝜌 𝑔 = Gas density, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝜌 𝑚 = Matrix density, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝜌 𝑔𝑠𝑡 = Gas density at standard condition, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝜎 𝑖𝑗 = Effective stress, m/L𝑡 2 , 𝑃𝑎 𝜎 𝑚 = Mean effective stress, m/L𝑡 2 , 𝑃𝑎 ∅ = Matrix in situ porosity ∅ 0 = Matrix porosity at reference condition ∅ 𝑓 = Fracture porosity