Modeling of Gas Production from Shale Reservoirs Considering Multiple Transport Mechanisms

Gas transport in unconventional shale strata is a multi-mechanism-coupling process that is different from the process observed in conventional reservoirs. In micro fractures which are inborn or induced by hydraulic stimulation, viscous flow dominates. And gas surface diffusion and gas desorption should be further considered in organic nano pores. Also, the Klinkenberg effect should be considered when dealing with the gas transport problem. In addition, following two factors can play significant roles under certain circumstances but have not received enough attention in previous models. During pressure depletion, gas viscosity will change with Knudsen number; and pore radius will increase when the adsorption gas desorbs from the pore wall. In this paper, a comprehensive mathematical model that incorporates all known mechanisms for simulating gas flow in shale strata is presented. The objective of this study was to provide a more accurate reservoir model for simulation based on the flow mechanisms in the pore scale and formation geometry. Complex mechanisms, including viscous flow, Knudsen diffusion, slip flow, and desorption, are optionally integrated into different continua in the model. Sensitivity analysis was conducted to evaluate the effect of different mechanisms on the gas production. The results showed that adsorption and gas viscosity change will have a great impact on gas production. Ignoring one of following scenarios, such as adsorption, gas permeability change, gas viscosity change, or pore radius change, will underestimate gas production.


Introduction
Due to the increasing energy shortage in recent years, gas production from shale strata has played an increasingly important role in the volatile energy industry of North America and is gradually becoming a key component in the world's energy supply [1,2]. Shale strata or shale gas reservoirs (SGR) are typical unconventional resources with a critically low transport capability in matrix and numerous "natural" fractures. Core experiments have shown that 90% of shale bedrock permeability are less than 150 nd, and diameters of the pore throat are 4~200 nm [3]. Natural gas has been stored in SGRs in three forms [2]: (a) as free gas in the micro fractures and nano pores; (b) as dissolved gas in the kerogen [4]; and (c) as adsorption gas in the surface of the bedrock. About 20%~85% of gas in SGRs is stored in form [5].
Gas transport in extremely low-permeability shale formations is a complex process with many co-existing mechanisms, such as viscous flow, Knudsen diffusion, slip flow [6,7], and gas adsorption-desorption. Some scholars have studied gas transport behavior in SGRs. E. Ozkan used a dual-mechanism approach to consider transient Darcy and diffusive flows in the matrix and stress-dependent permeability in the fractures for naturally fractured SGRs [8]. However, they ignored the effects of adsorption and desorption. Moridis considered Darcy's flow, non-Darcy flow, stress-sensitivity, gas slippage, non-isothermal effects, and Langmuir isotherm desorption [9]. They found that the production data from tight-sand reservoirs can be adequately represented without accounting for gas adsorption whereas there will be significant deviations if gas adsorption were omitted in SGRs. But they did not consider gas diffusion in the Korogen. Bustin studied the effect of fabric on gas production from shale strata. However, it was assumed the matrix does not have viscous flow or diffusion mechanisms. It was also assumed that gas transport in the fracture system obeys Darcy's law, which is not sufficiently accurate [10]. Yu-Shu Wu proposed an improved methodology to simulate shale gas production, but the gas adsorption-desorption were ignored in the model [11]. It is necessary to develop a mathematical model that incorporates all known mechanisms to describe the gas transport behavior in tight shale formations.
The simulation work was greatly simplified when the apparent permeability was first introduced by Javadpour. In 2009, the concept of apparent permeability, considering Knudsen diffusion, slippage flow, and advection flow, was proposed [4]. Through this method, the flux vector term can be simply expressed in the form of Darcy's equation, which greatly reduces the computation complexity. Then the concept of apparent permeability was further applied to pore scale modeling for shale gas [12,13]. Civan and Ziarani derived the expression for apparent permeability in the form of the Knudsen number [14,15] on a unified Hagen-Poiseuille equation [16].
In this paper, a general numerical model for SGRs was proposed which was constructed based on the dual porosity model (DPM). Mass balance equations were constructed for both matrix and fracture systems using the dusty gas model. In the matrix, Knudsen diffusion, gas desorption, and viscous flow were considered. Gas desorption was characterized by the Langmuir isothermal equation. In the fracture, viscous flow and non-darcy slip permeability were considered. The increase of the pore radius due to gas desorption was calculated from the gas desorption volume from the pore wall. Gas viscosity was characterized as a function of the Knudsen number. Weak form equations were derived for the system based on the constant pressure boundary and got solved using COMSOL. A comprehensive sensitivity analysis was conducted, and a detailed investigation was done to determine their impact on the shale gas production performance. The results showed that considering gas viscosity change greatly increased gas production under given reservoir conditions and slowed down the production decline curve. Considering pore radius increase due to gas desorption from the pore wall resulted in a higher production, but the effect was not very significant under the given reservoir conditions. In SGRs, both the matrix and fracture permeability changed during production. Ignoring one of these factors such as Knudsen diffusion, slippage effect, desorption, viscosity decrease, and porosity increase would lead to a lower cumulative production. Therefore, it is crucial to incorporate these factors into the existing models to obtain a more accurate shale gas production prediction. authors' adherence to all the PLOS ONE policies on sharing data and materials, as detailed online in the guide for authors.

Pore Distribution and Gas Transport Process in SGRs
Understanding the pore distribution and geometry of shale strata is of great importance in understanding the gas transport process in the media. Fig 1 shows the gas distribution and geometry of shale strata from micro to macro scales. It informs us that in the fracture, only free gas exists. And in the matrix which is full of kerogen, free gas and adsorption gas co-exist. In such a complicated system, gas desorbs from the pore wall and transports in the matrix system. Due to pressure difference between the matrix and fracture system, gas transfers between the matrix and fracture. Then, gas flows to the wellbore and is produced to the surface [3]. This model constructed in this study was inspired according to this gas transport process in shale strata.

Model Construction
This theoretical study was based on the following basic assumptions: 1. Single component, one phase flow in SGRs; 2. Ignore the effect of gravity and heterogeneity on the gas flow; 3. Ideal gas behavior for the natural gas in shales, and gas deviation factor z = 1; 4. Formation rocks were incompressible, and the porosity change due to rock deformation was ignored; 5. Isothermal flow process was present in the whole reservoir life; 6. Gas adsorption-desorption kinetics obeyed the Langmuir curve, and can achieve equilibrium state quickly at any reservoir pressure (diffuse quickly).
Consider a generalized mass balance equation [16] in every grid block is as follows: where M is mass accumulation term; u is velocity; ρ is density; Q is source term; and t denotes time. For shale strata which are full of micro-fractures, two general methods can be used to deal with the fracture-matrix interaction. One is the dual porosity continuum method (DPCM), including dual porosity and dual permeability (DPDM) [17], or MINC (multiple interacting continua) [18]. The other one is the explicit discrete fracture model [19]. Though the later model is more rigorous, it is still limited in applications to field study due to its computational intensity and lack of detailed knowledge of fracture distribution [14]. In this study, the DPCM, which is shown in the Fig 2, was used in this study to describe the complex fracture distribution.
For the DPCM, there are two mass balance equations that correspond to fracture and matrix systems, respectively, as indicated by Warren and Root [17].
With the subscripts f and m representing fracture and matrix system, respectively, the two sets of equations are illustrated as follows: The first term on the left side is the mass accumulation term; the second term on the left side is the flow vector term; and the right side is the source/sink term.These terms are addressed correspondingly in the following sub-sections.

Mass accumulation term
The general form of the mass accumulation term is: where β denotes the fluid phase, ϕ is porosity, S β is the faction of pore volume occupied by the phase β, ρ β is the density of the phase β. Specifically for a gas reservoir, S β = 1. However, in matrix system of shale strata, there are adsorption gas and free gas in the system which is shown in Fig 1. The free gas, which occupies the pores of the matrix, can be represented by M free ¼ ϕ P β S β ρ β , and the adsorbed gas in the surface of the bulk matrix can be represented as M adsorp = ∑(1 − ϕ)q a , where q a is the adsorption gas volume per unit bulk volume.
Gas desorption occurs when the pressure difference exists in organic grids (kerogen). With free gas production, the pressure in the pores decreases, which results in the pressure difference between the bulk matrix and the pores. Due to this pressure drop, the gas desorbs from the surface of the bulk matrix. The most commonly used empirical model that describes the sorption and desorption of gas in shale and provides a reasonable fit to most experimental data is the Langmuir single-layer isotherm model [12,14,20,21], which is expressed in Eq (5): where v L is the Langmuir volume, denoting the amount of gas sorbed at infinite pressure p 1 , and p L is the Langmuir pressure, corresponding to the pressure at which half of the Langmuir volume v L is reached. Normally, v L and p L are measured under standard temperature and pressure (STP) conditions. So, the adsorption gas volume per unit of bulk volume can be expressed in Eq (6): where q a is the adsorption gas per unit area of the shale surface, kg/m 3 ; V std is the mole volume under standard condition (0°C, 1atm), m 3 /mol; q std is the adsorption volume per unit mass of shale, m 3 /kg; V L is the Langmuir volume, m 3 /kg; p L is the Langmuir pressure, P a ; and ρ s is the density of shale core, kg/m 3 . The mass accumulation considering free gas and adsorption gas can be represented as M = ∑(ϕρ β + (1 − ϕ)q a ), and the corresponding partial differential form is Based on the equation of state (EOS): As there is only free gas in the fracture system, the mass accumulation term can be described as follows:

Flow vector term
This paper distinguishes the gas flow in shale gas reservoirs by flowing media: matrix and fracture systems. Due to the differences of pore sizes in those two media, gas flow mechanisms are different. Fig 3 shows the general flow process from matrix to fracture, then to wellbore in the shale strata.

Flow vector term in the matrix
The general Darcy's law informs us the relationship between velocity and pressure drop, as presented in the Eq (10): For low density fluids such as gases, it is assumed that the effect of gravity can be ignored. Therefore, a simplified empirical form of Darcy's law can be used for the flow vector term in pores ranging from tens to hundreds of microns: where K * is the permeability tensor. In the matrix system, where pores are in the range of nanometers, the conventional Darcy's law cannot be used to describe the flow process. Bird et al. concluded that gas transportation in nano pores is a multi-mechanism-coupling process that includes Knudsen diffusion, viscous convection, and slip flow.   mechanisms of gas transport in the nano pores of shale strata [22]. The following subsections describe those flow mechanisms further.
(1). Viscous flow. When the mean free path of the gas is smaller than the pore diameter (Knudsen number is far less than 1), then the motion of the gas molecules is mainly affected by the collision between molecules. There exists viscous flow caused by the pressure gradient between the single component gas, which can be described by Darcy's law [23]: where J v is the mass flow (kg/(m 2 Á s)), ρ m is gas density in bedrock (kg/ m 3 ), k mi is the intrinsic permeability  of bedrock (m 2 ), μ g is the gas viscosity (Pa Á s), and p m is the pore pressure in the bedrock (Pa).
(2). Knudsen Diffusion. When the pore diameter is small enough so that the mean free path of the gas is close to the pore diameter (i.e., Knudsen number > 1), the collision between the gas molecules and the wall surface dominates. The gas mass flow can be expressed by the Knudsen diffusion [23]: where J k is the mass flow caused by Knudsen diffusion (kg/(m2 Á s)), M g is the gas molar mass (kg/mol), D km is the diffusion coefficient of the bedrock (m2/s), and C m is the gas mole concentration (mol/ m 3 ).
(3). Slip flow. In low-permeability formations (less than 0.001 md) or when the pressure is very low, gas slip flow cannot be omitted when studying gas transport in tight reservoirs [6,24]. Under such kind of flow conditions, gas absolute permeability depends on gas pressure, which can be expressed as follows: where k a is the apparent permeability, k i * is the intrinsic permeability, and b is the slip coefficient. Some empirical models have been developed to account for slip-flow and Knudsen diffusion in the form of apparent permeability, including correlations developed from the Darcy matrix permeability [8], correlations developed based on flow mechanisms [3,4,25], and semiempirical analytical models by Moridis et al. in 2010 [9]. This research adopted Javadpour's method, which considered slip flow, viscous flow, and Knudsen diffusion [4]. The apparent permeability is given as follows: where α is the tangential momentum accommodation coefficient (TMAC), which characterizes the slip effect. α is a function of wall surface smoothness, gas type, temperature, and pressure, which varies from 0 to 1 [4,26,27]. In the shale matrix with pores in nano to micro scale, it is necessary to consider the effect of slip effect and Knudsen diffusion. Cao et al. (2005) also pointed out that when the gas free mean path is close to the pore radius, the flow falls in the slip regime [28]. The TMAC α in the Eq (16) can be used to characterize the slip effect or rarefaction effect. TMAC represents the part of gas molecules reflected diffusely from the tube wall relative to specular reflection [4]. The slip length can be related to the TMAC using equation l ¼ 2Àa a Kn [28]. From this we can find that under a certain Kn, the smaller the TMAC α, the larger the slip length. More gas molecules will under slip flow regime. Arkilic et al. obtained the α about 0.8 for light gases flowing in silicon microchannel [28,29]. Javadpour has also employed 0.8 as the value for tangential momentum accommodation coefficient in deriving gas permeability in shale matrix [4]. For simplification and consistency, we also employ 0.8 as the TMAC value in this paper. The accurate TMAC value for a specific gas in the mudrock can be obtained using lab experiments or numerical methods, such as Molecular Dynamics method, Monte Carlo method, and direct simulation Monte Carlo method, etc.
So, the flow vector term in the matrix is: For the matrix system, it should be noted that with the gas desorption from the pore wall, the pore radius increases. The volume of gas desorbed from the wall can be characterized by the Langmuir isotherm curve. Considering the gas desorption, the effective pore radius can be expressed as follows [18]: Here, r eff is the initial pore radius, d CH 4 is the molecular diameter of methane, P L is the Langmuir pressure, and r max is the maximum pore radius, which is the initial pore radius plus the molecular diameter. With the matrix pressure decreasing during the production process, gas desorbs from the pore wall. So, from Eq (17), we can find that the pore radius is increasing with production. If the matrix pressure can approach 0, which is the ideal scenario, the final pore diameter will be equal to initial radius r i . The process is shown in Fig 5.It can be found that after the gas desorbed from the pore wall,the gas transport channel will become larger.Due to the increase of the pore radius,the collision of gas moleculars with the pore wall will become less dominant.Eq (16) has shown this change.With pore radius increase,the slip effect coefficient b m will decrease. The contribution of Knudsen diffusion and slip flow will contribute less to the gas flux with pore radius increases. At the same time, with the pore radius increase, the viscous flow which is denoted using the blue dots in Fig 4 will become more dominant and contribute more to the gas flux. The total gas flux is the combined effect of these effects. Flow vector term in the fracture In this study, Knudsen diffusion and viscous flow are considered for the gas flow in the frature system. The mass flux can be expressed as the summation of the two mechanisms. Therefore, Or Using the form of conventional Darcy's flow equation, the fracture apparent permeability can be expressed as: where p f is the fracture pressure, k fi is the initial fracture permeability, k f is the apparent fracture permeability, b f is the Klinkenberg coefficient for the fracture system, D kf is the Knudsen diffusion coefficient for the fracture system, and ∅ f is the fracture porosity. Combining all of the above, the flow vector term in the fracture is:

Source and sink term
For DPCM, it is of great importance to consider the interaction between the matrix and the fracture in simulation. There are different methods for handling this issue in reservoir simulation. The boundary condition approach explicitly calculates the amount of the matrix-fracture petroleum fluid transfer by imposing boundary condition at each time step. It is very useful in well testing [8,30]. However, this method is impractical in full field simulation due to its high computational cost. Another method is the Warren-Root method [17,30]. The Warren-Root method calculates the cross flux between the fracture and matrix systems by assuming a pseudo-steady state flow between the matrix and fracture when the no-flow boundary is confined. The gas transfer between the matrix and fracture systems is represented in Eq (24): where s ¼ 4 1 is the crossflow coefficient between the matrix and fracture systems, and L x , L y , and L z are the fracture spacings in the x,y, z directions, respectively. In the fracture system, there exists a source term that flows into fracture from the matrix and a sink term that flows out of the fracture into the wellbore. Therefore, for the matrix system, the sink term that flows out of the matrix to the fracture can be described as follows: The sink term followed the model developed by Aronofsky and Jenkins [31], which considers gas production from vertical well and can be expressed using Eq (26): In Eq (26), when the production well is in the corner, θ = π/2. When the production well is in the center, θ = 2π [32]. In addition, p wf is the bottomhole flowing pressure; p f À is the average pressure in the fracture system; r w is the well radius; and r e is the drainage radius, which can be expressed as follows: where Δx, Δy, k x , and k y are the length of grid and permeability in the x and y directions. For the fracture sytem, Q f = T − q p .

Gas viscosity change
Currently, some reservoir simulators have considered gas viscosity change, which is a function of pressure. However, this study considered the gas viscosity change to be a function of the Knudsen number (the ratio of the free molecular length to the pore diameter), which is dependent on the gas transport period, as discussed previously. Gas viscosity changes with the Knudsen number in the production process [33]. Fig 6 shows that the ratio of the effective viscosity to the initial viscosity changes with the Knudsen number, which is a function of gas pressure. Eq (28) shows the relationship between gas viscosity and the Knudsen number. Beskok and Karniadakis (1999) have derived the expression of the Knudsen number using a more general form, as shown in Eq (29), which is easier to be considered in the numerical simulation [16]. Also, the results of this analytical equation has good agreement with DSMC results and with the linearized Boltzmann solution.
where μ 0 is the initial viscosity at the initial reservoir pressure (p m = p f = p i ). This suggested a Knudsen dependence of the rarefaction parameter β and provided an analytical expression for this dependence. Eq (29) expresses the definition of the Knudsen number, which can be related to be the basic parameters that change with pressure: The final mathematical model to characterize gas transport in SGRs can be expressed as follows: With following boundary conditions considered in this paper: Boundaryconditionforfracture : In this paper, the model Eqs (30)- (34) and four versions of its simplification have been considered: (1) considering the basic DPM, which has considered none of adsorption, non-Darcy flow, gas viscosity, and pore radius change. (2) considering the basic DPM with adsorption which corresponds to the first term of the left part of Eq (30); (3) considering the basic DPM with adsorption and non-Darcy permeability change. That is, b m 6 ¼ 0 and b f 6 ¼ 0; (4) considering the basic DPM with adsorption, non-Darcy permeability change, and gas viscosity change which is represented in the Eq (28); and (5) considering adsorption, non-Darcy permeability change, gas viscosity change, and pore radius change. Pore radius change is represented in Eq (17), in which r eff represents the changing radius as shown in b m (r = r eff in Eq (16)).
The finite element method was used to solve the PDE equations from Eqs (30)- (34), which are listed above. First, multiply a test function v(x,y) on both sides of original Eq (30) where O is the domain; Γ 1 is outside boundary; and Γ 2 is the inner boundary.
Using Green's Theorem (Divergence Theory and Integration by parts in multi-dimension), we obtain that: If we just consider the solution on the domain O, the weak form for the governing Eq (30) is: Taking the similar procedures on governing Eq (31), we obtain:

Model Verification
As there is no real field data which is same with this therotical case, a common method to verify the righteousness is to compare the results of the numerical simulation against the analytical results. Wu et al. derived the analytical solution for 1D steady-state gas transport [34], which is shown in Eq (41). Thus, we can verify our model by comparing the results from simpilification version of the model in this paper with the analytical results. The reservoir properties and Klinkenberg properties used in this verification were obtained from the experimental study of the welded tuff at Yucca Mountain.
where b is the Klinkenberg factor, 7.6×10 5 Pa; p(L) is the gas pressure at the outlet (x = L), 1.0 × 10 5 Pa; q m is the air injection rate, 1.0 × 10 −6 kg/s; μ g is the gas dynamic viscosity, 1.84 × 10 −5 Pa.s; L is the length from the inlet (x = 0) to the outlet, 10 m; x is the task location along the gas transport path, m; k 1 is the intrisic permeability for welded tuff atYucca Mountain, 5.0 × 10 −19 m 2 ; and β is the compressibility factor, 1.8 × 10 −5 Pa -1 m -1 . Fig 7 presents the match results between the model presented in this paper and the analytical solution. As shown in Fig 7, the results obtained using the model presented in this paper agrees well with the analytical solution [34] provided by Wu et al., which validated our mathematical model.

Results and Analysis
The 2-D reservoir model is shown in Fig 8. For simplification, we only studyed the ¼ area of the whole reservoir. The reservoir parameters are shown in Table 1. The parameters used in this paper were selected based on the literature review regarding shale gas simulation [10]. The simulated reservoir is located at a depth of 5463 ft with a pressure gradient of 0.53 psi/ft and a temperature gradient of 0.065 K / ft, which corresponded to the initial reservoir pressure and  Table 1.
In this study, firstly we investigated the impact of parameters such as the initial reservoir pressure, matrix permeability, fracture permeability, matrix porosity, and fracture porosity on shale gas production. Then, the mechanisms discussed earlier were gradually incorporated into the simulation model to investigate their impact on gas transport in the shale strata. These  mechanisms are adsorption, permeability change due to comprehensive slip effect, viscosity change, and radius change due to gas adsorption.

Effect of Gas adsorption
In this study, we compared scenarios in which adsorption was considered and was not considered. When adsorption was considered, there is an adsorption term on the left hand side of Eq (30). Only free gas is considered in the matrix system if adsorption is ignored. The following analysis are all based on the scenario in which adsorption was considered. First, we need to know the effect of adsorption on the gas production. Fig 9 illustrates the effect of adsorption on the production rate and cumulative production. It is obvious that adsorption has a great effect on the production rate and cumulative production. Ignoring the effect of adsorption gas leads to great underestimations in the gas production in such situations.

Effect of non-Darcy flow
Effect of non-Darcy flow can be considered by adjusting the slip coefficient which appereas in Eq (30). If non-Darcy flow is considered, the slip coefficient b m = b f = 0. Fig 10 illustrates the matrix permeability and fracture permeability change with time. From the Fig 10, it is clear that the gas permeabilities in the matrix and fracture gradually increase during pressure depetion. The matrix permeability has a greater increase then the fracture permeability. This result has also validated the righteous of our model. The fluid flow channels are larger in the fractures than those in matrix. So, the slip effect coefficient in the fracture (b f ) will be smaller compared with those in matrix (b m ). So, the change of gas permeability will also be small. As shown in Fig  11, the impact of non-Darcy on the production rate and the cumulative production are plotted on the log-log scales. It is clear that considering non-Darcy flow increases the production rate and cumulative production; however, there is not a significant difference between these two scenarios. This is because though matrix permeability increases a lot as shown in Fig 10(A), however, fracture permeability is critical in determining fliud flow rate from the fracture system to the wellbore. Therefore, it is important to increase the fracture permeability rather than matrix permeability by including hydraulic fractures in using stimulations.

Effect of gas viscosity change
Effect of gas viscosity change was evalutaed on the basis of the previous study, which considered adsorption and non-Darcy flow. That is, if the gas viscosity change is considered, then b m 6 ¼ 0, b f 6 ¼ 0. The gas viscosity change is represented in Eq (28). If the gas viscosity change is not considered, then b m 6 ¼ 0, b f 6 ¼ 0, and the gas viscosity is a constant. As shown before, the gas viscosity decreases with production and pressure depletion. Fig 12 shows the comparison of the production rate and cumulative production between these two scenarios. It can be seen that gas viscosity has a great impact on the production rate and cumulative production under the production condition of this paper. From Eq (28), we can find that gas viscosity is a nonlinear function of the pressure. During the gas production process, pressure decrease leads to a decrease in the gas viscosity, which results in an increase in the production. Therefore, ignoring the effect of gas viscosity change decreases the estimated production.

Effect of pore radius change
Pore radius change is represented in Eq (17), which will change the radius shown in the part of b m (r = r eff ) in Eq (16). All the factors that have been considered before are still considered when analyzing the effect of pore radius change: adsorption, non-Darcy flow, and gas viscosity change. From Eq (17), we can find that the maximum pore radius in the production will be intrinsic radius plus the diameter of CH 4 . Therefore, the pore radius actually does not change too much. However, r eff is affected by the matrix pressure, which changes during reservoir depletion; therefore, it is still important to consider this effect. The effect of the pore radius change on the production rate and cumulative production are shown in Fig 13. From Fig 13, we can find that although there is a slight increase in the production rate and cumulative production. The effect of the radius change is not very significant. Fig 14 shows the comparison of the above 5 models: (1) basic DPM which has not considered adsorption, non-Darcy flow, gas viscosity, and pore radius change; (2) Considering adsorption; (3) Considering adsorption and non-Darcy permeability change; (4) Considering adsorption, non-Darcy permeability change, and gas viscosity change; (5) Considering adsorption, non-Darcy permeability change, gas viscosity change, and pore radius change. It can be found that the mechanisms that have considered will increase the production. Ignoring any one of these factors will lead to an underestimated gas production.

Sensitivity Analysis
To study the effect of reservoir parameters on gas production, we used the model that considered the effect of adsorption and non-Darcy flow. The effect of the following reservoir parameters was analyzed. (1) initial pressure; (2) matrix permeability; (3) fracture permeability; (4) matrix porosity; and (5) fracture porosity. For each scenario, we compared the production rate and cumulative production for three stages. In this study, when discussing the effect of certain parameters, the other parameters were kept the same as those listed in the Table 1.

Effect of initial pressure
Initial reservoir pressure is very important in the production, especially for gas reservoir which has no other driving force. Here, we analyzed the effect of the initial pressure on gas production. We considered the scenarios of the initial pressure being equal to 0.2 MPa, 2.0 MPa, and 20 MPa. Other parameters were kept the same as those listed in the Table 1. Fig 15 shows the comparison of the production rate and cumulative production under different initial reservoir pressures. From the plots, we can find that initial pressure has a great effect on the gas production; there is 3-to 5-fold increase in the cumulative production when the intial pressure increases 10-fold. With the increase of reservoir initial pressure, gas production rate and gas cumulative production increase.

Effect of matrix and fracture permeability
Permeability is the king in the reservoir production. For DPM, matrix permeability and fracture permeability are not the same and the effects of these two kinds of permeability are also different. In this study, the effects of matrix permeability and fracture permeability on gas production were evaluated together to show which effect will be more important. Fig 16 shows the effects of matrix permeability when the permeability is varied from 1.0×10 −4 mD to 1.0×10 −5 mD to 1.0×10 −6 mD. Fig 17 shows the effects of fracture permeability when the permeability changes from 0.1 mD to1 mD to 10 mD. From these figures, we can find that the production rate and cumulative production increase when the matrix permeability or fracture permeability increases. However, from the comparison between Fig 16(A) and Fig 16(B), and between Fig  17(A) and Fig 17(B), we can find that increasing fracture permeability has a more apparent effect on the production rate and cumulative production compared with the matrix permeability. This result has validated that the fracture system is the mian fluid flow channels which is the also the assumption for the dual porosity model. In the dual porosity model, the matrix system is the main storage space and the fracture system is the main fluid flow space.

Effect of matrix and fracture porosity
Porosity decides the gas amount which can be stored in the reservoir. Study on the effect of matrix and fracture porosity can help us choose target reservoir and speed history match. For studying the effect of porosity, three levels of matrix porosity and fracture porosity were considered and evaluated. We have compared the difference when we change the matrix porosity from 1% to 10% to 20% and the difference when we change the fracture porosity from 0.01% to 0.1% to 1%. From Figs 18 and 19, we can find that porosity increase leads to production increase. Also, from the comparison between Fig 18(A) and Fig 18(B), and between Fig 19(A) and Fig 19(B), it can be found that matrix porosity has a more siginificant effect on shale gas production compared with fracture porosity, which conforms to the assumption of dual porosity model.

Conclusion
This paper presents a theoretical model and mechanism study for shale gas production. Following conclusions were obtained from this study: A new mathematical model that considers adsorption, non-Darcy permeability change, gas viscosity change, and pore radius increase due to gas desorption was constructed; Results analysisshowed that considering one of following scenarios: adsorption, non-Darcy permeability change, gas viscosity change, and pore radius change increases the production estimate. Among these mechanisms, adsorption and gas viscosity change have a great impact on gas production. Ignoring one of these effects decreases gas production; Sensitivity analysis of the reservoir parameters showed that initial reservoir pressure has a great impact on gas production. Fracture permeability has a more important effect than the matrix permeability. Porosity increase leads to the increase of gas production. However, matrix porosity is more important than fracture porosity.