Regular Article
Physical and mathematical modeling of gas production in shale matrix
Petrochina Research Institute of Petroleum Exploration & DevelopmentLangfang,
065007
Langfang  China
^{*} Corresponding author: yurz201169@petrochina.com.cn
Received:
12
June
2017
Accepted:
8
March
2018
Shale gas mainly stores in shale matrix, and gas production in shale matrix is very important during exploration. In order to clarify gas production and transport mechanism in shale matrix, an experimental modeling of gas production in shale matrix was designed and conducted with Longmaxi shale samples collected from South of Sichuan. The experimental results show that gas production decline curve displays a “L” pattern which indicates initial production is high and declines rapidly, while late production is low and declines moderately; meanwhile, pressure propagation in shale matrix is quite slow due to ultralow permeability. Based on the results, a mathematical model was derived to describe gas production in shale matrix. The comparison between numerical solution of mathematical model and experimental results shows that the mathematical model can well describe gas transport in shale matrix. In addition, factors affecting gas production were investigated on the basis of the mathematical model. Adsorbed gas can replenish gas pressure in pores by desorption and delay pressure propagation, and gas production decreases very quickly when there is no adsorbed gas. Other parameters (diffusion coefficient, permeability and porosity) also need to be considered in shale gas development.
© W. Guo et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Hydraulic fracturing technology is widely applied in shale reservoirs to significantly increase gas production. Due to extremely low permeability of shale matrix, shale gas production is strongly relied on gas supply of shale matrix. So it is very important to study gas production in shale matrix during exploration. Shale matrix is mainly composed of clay minerals and organic matters, which provide flow channels with different sizes (Curtis, 2002; Loucks et al., 2009; Passey et al., 2010; Liu et al., 2011; Pan et al., 2015; de Swaan, 2014; RomeroSarmiento et al., 2015).
In order to make clear gas flow mechanism in shale nanopores, many scholars made theoretical studies and obtained some useful information. Javadpour found that gas flow in shale nanopores can be modeled with a diffusive transport regime and negligible viscous effects (Javadpour et al., 2007); in 2009, presenting another formulation for gas flow in the nanopores of mudrocks based on Knudsen diffusion and slip flow, he drew a conclusion that the ratio of apparent permeability to Darcy permeability increases sharply as pore sizes reduce to smaller than 100 nm (Javadpour, 2009). Civan presented a model which can accommodate a wide range of fundamental flow mechanisms, such as continuum, slip, transition, and free molecular flow, depending on the prevailing flow conditions characterized by Knudsen number (Civan, 2010). After improving the description of shale matrix flow by considering diffusive flow in nanopores, Ozkan found that Knudsen flow takes over and contributes to fluids transfer from matrix to fracture network (Ozkan et al., 2010). Appling the dustygas model to the study of microscale flow behavior in shale gas system, Freeman showed that for rock with small pore throat diameters, the composition of flowing gas is significantly different than the situ stagnant conditions (Freeman et al., 2011). Villazon generalized a transport equation valid for all flow regimes for an ideal gas in a capillary tube and extended the formulation to quantify the effect of a distribution of pore sizes for different flow regimes (Villazon et al., 2011). Swami studied the contribution of Knudsen diffusion, gas slippage, gas desorption and gas diffusion from kerogen to total shale gas production by a proposed pore scale gas flow model (Swami et al., 2012). Shabro introduced a new numerical algorithm to forecast gas production in organic shale, simultaneously taking into account gas diffusion in kerogen, slip flow, Knudsen diffusion, and Langmuir desorption (Shabro et al., 2012). Alharty developed a dualporosity and a tripleporosity finitedifference model for singlephase flow in shale including advective, diffusive, and Knudsen flow mechanisms (Alharty et al., 2012). Using Lattice Boltzmann Method (LBM) to simulate shale gas transport in kerogen, Fathi drew a conclusion that mechanisms of slippage and surface transport would enhance gas transport in the organic nanopores (Fathi et al., 2012). Deng presented a seepage model for shale gas reservoir based on Beskok and Karniadakis equation (Beskok and Karniadakis, 1999; Deng et al., 2014). Cao investigated multiscale gas transport behavior in organic rich shale by theoretical modeling (Cao et al., 2016). Albinali constructed a 1D anomalous diffusion model to analyze oil and gas production in fractured nanoporous media (Albinali et al., 2016).
In summary, previous studies mainly focused on theoretical calculation in a micro scale. Those previous investigations did not obtain the experimental data to answer questions concerning gas production in shale matrix or pressure propagation in shale matrix. In this study, we conducted an experiment to model gas production in shale matrix and studied gas transport mechanism in shale matrix in a macro scale.
2 Materials and experiments
2.1 Sample preparation
Experimental sample was collected from the target layer of Longmaxi Formation in Southern Sichuan Basin because Longmaxi Shale is relatively stable and can better represent the formation. Core plug with diameter of 10.5 cm and length of 14.2 cm was collected in a cored well. The sample was put into experimental apparatus immediately when the sample was lifted from cored well after cleaning the drilling fluid on the sample surface. The cored depth of the sample is 1408 m; other samples along with the cored process were collected for further study of petrophysical characteristics.
2.2 Apparatus
The experimental apparatus is shown in Figure 1. The apparatus is an assemblage of seven components. They are an axis pressure pump, a confining pressure pump, a core holder, a gas injection system, a backpressure regulator, a flowmeter, and a data acquisition system. The core holder is made of highpressure resistant stainless steel with a rubber sleeve and sealing ring. Radial and confining pressures are exerted by axis pressure pump and confining pressure pump respectively. The gas inlet and outlet are equipped with pressure sensors. By using the ISCO pump, methane gas can be injected into the core. The data acquisition system records the overburden pressure, inlet, outlet pressures and backpressure regulator pressure.
Fig. 1 Experimental apparatus of gas production in shale matrix. 
2.3 Experiment procedure
Four steps are needed to model gas transport in shale matrix as follows:

measure diameter and length of the sample and sample weight, then put the sample into the core holder and adjust confining pressure and axis pressure to 40 MPa;

open valve 1 and valve 2 to let methane flow into core holder and saturate shale sample, both the inlet and outlet pressure are maintained 29 MPa by ISCO pump. The purpose of this step is to recover underground shale gas storage pressure condition. When the ISCO pump fluid does not change within 3 days, and it can be considered that the sample is fully saturated. The time for saturating in this experiment is 240 days;

close valve 1 and valve 2, set backpressure regulator pressure to 2 MPa, open valve 2 and record inlet pressure, outlet pressure, gas flow rate and time;

in order to study the effects of different working modes on gas production, backpressure valve is removed after 37 days and outlet pressure is changed to atmospheric pressure.
3 Results and discussion
3.1 Petrophysical characteristics of parallel samples
Table 1 presents compositional analysis of two parallel samples. Generally speaking, clay mineral content is the highest; the second highest mineral is quartz. Overall, the average values of clay, quartz and calcite content accounts for 87.5% of the whole composition minerals. Total organic content of the sample 1 is 1.88% and total organic content of the sample 2 is 2.02%.
Figure 2 shows adsorptiondesorption isotherms of the two parallel samples, according to isothermal experiments conducted at liquid nitrogen temperature using nitrogen as the adsorbent. Both adsorption and desorption isotherms show hysteresis of type IV (as per IUPAC classification). During liquid nitrogen adsorption process, when the relative pressure is greater than 0.4, the multilayer adsorption will occur and be accompanied by capillary condensation process. And desorption is the inverse process of adsorption. The hysteresis loops in the multilayer adsorption range with physical adsorption isotherms are usually related to capillary condensation and can be well reproduced. Due to the abundant 2 ∼ 10 nm mesopores with different shapes in shale matrix, the hysteresis loop between adsorption curve and desorption curve will appear and close at p/p_{0} around 0.45 ∼ 0.5. Due to the different pore size distribution, the hysteresis loops of two samples are different but both close at p/p_{0} around 0.45 ∼ 0.5. The hysteresis pattern, termed H3 as per IUPAC classification, indicates presence of slitlike pores. Figure 3 shows the pore size distribution obtained from desorption curves by BJH method. The distribution curves of two shale samples show a peak radius between 1 and 3 nm. The dominance of nanopores in shale matrix makes gas store and transport different from conventional gas sand.
Mineral composition of two parallel samples.
Fig. 2 Liquid nitrogen adsorption and desorption experimental results. 
Fig. 3 Poresize distribution of shale samples. 
3.2 Experimental result
Figure 4 shows the result of inlet pressure and outlet pressure as time goes on. As is seen from the figure, pressure propagation in shale matrix is very slow, and pressure propagates to the inlet after 100 days. The pressure drop at the inlet is very slow after the pressure propagates to the inlet. Figure 6 is daily gas production, showing an “L” pattern which displays a high initial gas production and a very long tail. At the beginning of the experiment, the gas production was high and decreased fast. At the later stage of the experiment, gas production was low but decreased as very slowly. This was consistent with production decline characteristic of real shale gas well.
Fig. 4 Inlet and outlet pressure curves. 
Fig. 6 Comparison of experimental result and numerical calculation result. 
3.3 Mathematical model
In order to describe gas production in shale matrix, we derived the following mathematical model. Suppose shale matrix as homogeneous, the gas flow in shale matrix includes the seepage term and the diffusion term. Therefore, its total flow velocity in shale matrix v^{T} is the sum of Darcy flow velocity v^{D} and diffusion flow velocity v^{K}: $${v}^{T}={v}^{D}+{v}^{K}\text{.}$$(1)
The seepage term can be described by Darcy's law: $${v}^{D}=\frac{k}{\mu}\nabla p\text{,}$$(2) where k is permeability of shale matrix, m^{2}; µ is viscosity of methane, Pa · s; p is gas pressure, Pa.
The diffusion term can be described by Knudsen diffusion equation: $${v}^{K}=\frac{D}{{\rho}_{g}}\nabla C\text{,}$$(3) where D is diffusion coefficient, m^{2}/s; C is mass concentration of gas in shale matrix, kg/m^{3}; ρ_{g} is gas density, kg/m^{3}.
Shale gas exists in adsorption gas and free gas state. Take onedimensional case for an example, Assume shale volume is V_{r}, m^{3}, the shale rock density is ρ_{r}, kg/m^{3}, and gas in the shale can be calculated by the real gas state equation, and the molecular weight of gas is M_{g}, mol/kg, mass concentration of free gas and adsorbed gas in shale matrix is respectively as follows: $${C}_{\text{fre}}=\frac{{m}_{\text{fre}}}{{V}_{r}}=\frac{{\rho}_{g}{V}_{g}}{{V}_{r}}=\frac{{\rho}_{g}{V}_{r}\phi \left(1{S}_{W}\right)}{{V}_{r}}={\rho}_{g}\phi \left(1{S}_{W}\right)\text{,}$$(4) $${C}_{\text{ads}}=\frac{{m}_{\text{ads}}}{{V}_{r}}=\frac{{\rho}_{r}{V}_{r}{V}_{\text{des}}{\rho}_{a}}{{V}_{r}}={\rho}_{r}{V}_{\text{des}}{\rho}_{a}\text{,}$$(5) $$C={C}_{\text{fre}}+{C}_{\text{ads}}\text{,}$$(6)
where C_{fre} and C_{ads} respectively refer to mass concentration of free gas and mass concentration of adsorption gas, kg/m^{3}; m_{fre} is mass of free gas kg; m_{ads} is mass of adsorption gas, kg; V_{r} is shale bulk volume and V_{g} is gas volume, m^{3}; ρ_{r} is the shale density, kg/m^{3}; ρ_{a} is gas density under standard state, kg/m^{3}; V_{des} is adsorption gas that involves in desorption, m^{3}/kg; S_{w} is the water saturation of shale, dimensionless; ϕis porosity, dimensionless. Many scholars generally believe that shale isothermal desorption curve and isothermal adsorption curve are reversible, and Langmuir model can be used to calculate isothermal adsorption and desorption curves (Lu et al., 1995; Chareonsuppanimit et al., 2012). But in our previous study, we found that shale adsorption isothermal curve and isothermal desorption curve were irreversible, and desorption curve should be calculated by desorption equation (Guo et al., 2013). Therefore V_{des} should be calculated by desorption model as: $${V}_{\text{des}}=\frac{{V}_{d}p}{{p}_{d}+p}+c\text{,}$$(7) where V_{d}, p_{d} and c are constant. V_{d} is the maximum adsorption capacity during desorption process, m^{3}/t; p_{d} is a comprehensive constant, Pa^{−1}; c is the residual adsorption capacity, m^{3}/t.
Gas density can be calculated by equation of state: $${\rho}_{g}=\frac{p{M}_{g}}{ZRT}\text{.}$$(8)
Substitute Equations (2) and (3) into (1), we obtain: $${v}^{T}=\left(\frac{k}{\mu}\nabla p+\frac{D}{{\rho}_{g}}\nabla C\right)\text{.}$$(9)
The continuity equation for gas flow is: $$\nabla \left({\rho}_{g}{v}^{T}\right)=\frac{\partial \left[{\rho}_{g}\phi \left(1{S}_{W}\right)\right]}{\partial t}\frac{\partial {q}_{d}}{\partial t}+Q\text{,}$$(10) where q_{d} is total adsorption gas that involves in desorption, Q is source or sink term, which means gas production or injection. $${q}_{d}=\frac{{\rho}_{r}{M}_{g}}{{V}_{\text{std}}}\left(\frac{{V}_{d}p}{p+{p}_{d}}+c\right)\text{,}$$(11)where V_{std} is gas mole volume in standard state, 22.4 × 10^{−3}m^{3}/mol. $$\frac{\partial {q}_{d}}{\partial t}=\frac{\partial {q}_{d}}{\partial p}\frac{\partial p}{\partial t}\text{,}$$(12) $$\frac{\partial {q}_{d}}{\partial p}=\frac{{\rho}_{r}{V}_{d}{M}_{g}{p}_{d}}{{V}_{\text{std}}{\left(p+{p}_{d}\right)}^{2}}\text{.}$$(13)
Substitute Equations (4)–(13) into Equation (10): $$\frac{k{M}_{g}}{RT}\nabla \left(\frac{p}{\mu Z}\nabla p\right)+\frac{D\varphi \left(1{S}_{W}\right){M}_{g}}{RT}\nabla \text{<ce:glyph name="e;rad"e;/>}\u2022\nabla \left(\frac{p}{Z}\right)+\frac{D{V}_{d}{p}_{d}{\rho}_{r}{\rho}_{a}}{{\left(p+{p}_{d}\right)}^{2}}\frac{\partial p}{\partial t}=\frac{\varphi \left(1{S}_{W}\right){M}_{g}}{RT}\frac{\partial}{\partial t}\left(\frac{p}{Z}\right)+\frac{{\rho}_{r}{V}_{d}{p}_{d}{M}_{g}}{{V}_{\text{std}}{\left(p+{p}_{d}\right)}^{2}}\frac{\partial p}{\partial t}+Q\text{.}$$(14)
Suppose one dimension flow, Equation (14) can be rewritten as: $$\frac{k{M}_{g}}{RT}\frac{\partial}{\partial x}\left(\frac{p}{\mu Z}\frac{\partial p}{\partial x}\right)+\frac{D\phi \left(1{S}_{W}\right){M}_{g}}{RT}\frac{{\partial}^{2}}{\partial {x}^{2}}\left(\frac{p}{Z}\right)=\frac{\phi \left(1{S}_{W}\right){M}_{g}}{RT}\frac{\partial}{\partial t}\left(\frac{p}{Z}\right)+\left(\frac{{\rho}_{r}{V}_{d}{p}_{d}{M}_{g}}{{V}_{\text{std}}{\left(p+{p}_{d}\right)}^{2}}\frac{D{V}_{d}{p}_{d}{\rho}_{r}{\rho}_{a}}{{\left(p+{p}_{d}\right)}^{2}}\right)\frac{\partial p}{\partial t}+Q\text{.}$$(15)
3.4 Mathematical model solution
Numerical solution method is used to solute Equation (15), which refers to Civan (Civan et al., 2011; Civan et al., 2012) With timeforward difference, the first term on the left side of Equation (15) can be discretized as: $$k\frac{\partial}{\partial x}\left[\frac{p}{\mu z}\frac{\partial p}{\partial x}\right]=k\left({\left(\frac{p}{\mu z}\frac{\partial p}{\partial x}\right)}_{i+1/2}^{n}{\left(\frac{p}{\mu z}\frac{\partial p}{\partial x}\right)}_{i1/2}^{n}\right)/\mathrm{\Delta}{x}_{i}\text{.}$$(16) Assume: $$\text{\hspace{0.17em}}{\lambda}_{i+1/2}={\left(\frac{p}{\mu z}\right)}_{i+1/2}^{n}=\left({\left(\frac{p}{\mu z}\right)}_{i+1}^{n}+{\left(\frac{p}{\mu z}\right)}_{i}^{n}\right)/2\text{,}$$(17) $${\lambda}_{i1/2}={\left(\frac{p}{\mu z}\right)}_{i1/2}^{n}=\left({\left(\frac{p}{\mu z}\right)}_{i}^{n}+{\left(\frac{p}{\mu z}\right)}_{i1}^{n}\right)/2\text{,}$$(18) $${\left(\frac{\partial p}{\partial x}\right)}_{i+1/2}^{n}=\frac{{p}_{i+1}^{n}{p}_{i}^{n}}{{x}_{i+1}{x}_{i}}=\frac{{p}_{i+1}^{n}{p}_{i}^{n}}{\mathrm{\Delta}{x}_{i+1/2}}\text{,}$$(19) $${\left(\frac{\partial p}{\partial x}\right)}_{i1/2}^{n}=\frac{{p}_{i}^{n}{p}_{i1}^{n}}{{x}_{i}{x}_{i1}}=\frac{{p}_{i}^{n}{p}_{i1}^{n}}{\mathrm{\Delta}{x}_{i1/2}}\text{,}$$(20) $$\mathrm{\Delta}{x}_{i1/2}=\mathrm{\Delta}{x}_{i+1/2}=\mathrm{\Delta}{x}_{i}=\mathrm{\Delta}x\text{.}$$(21)
The second term on the left side of the Equation (15) can be discretized as: $$D\phi \left(1{S}_{W}\right)\frac{{\partial}^{2}}{\partial {x}^{2}}\left(\frac{p}{z}\right)=D\phi \left(1{S}_{W}\right)\left[{\left(\frac{p}{z}\right)}_{i+1}^{n}2{\left(\frac{p}{z}\right)}_{i}^{n}+{\left(\frac{p}{z}\right)}_{i1}^{n}\right]/\mathrm{\Delta}{x}_{i}^{2}\text{.}$$(22)
The first term on the right side of Equation (15) can be discretized as: $$\phi \left(1{S}_{W}\right)\frac{\partial}{\partial t}\left(\frac{p}{Z}\right)=\frac{\phi \left(1{S}_{W}\right)}{\mathrm{\Delta}t}\left[{\left(\frac{p}{Z}\right)}_{i}^{n+1}{\left(\frac{p}{Z}\right)}_{i}^{n}\right]\text{.}$$(23)
The second term on the right side of Equation (15) can be discretized as: $$\left(\frac{BA}{{\left(p+{p}_{d}\right)}^{2}}\right)\frac{\partial p}{\partial t}=\left\{\left(\frac{BA}{{\left({p}_{i}^{n}+{p}_{d}\right)}^{2}}\right){p}_{i}^{n+1}\left(\frac{BA}{{\left({p}_{i}^{n}+{p}_{d}\right)}^{2}}\right){p}_{i}^{n}\right\}/\mathrm{\Delta}t\text{,}$$(24) where: $$A=\frac{{\rho}_{r}{\rho}_{a}D{V}_{d}{p}_{d}RT}{{M}_{g}}\text{,}$$(25) $$B=\frac{{\rho}_{r}{V}_{d}{p}_{d}RT}{{V}_{\text{std}}}\text{.}$$(26)
Substituting Equations (16)–(26) into (15), we can get the finite difference equations with timeforward difference. It can be seen from Figure 3 that pressure change at the inlet is very slow and pressure propagates to the inlet after 100 days, so we presume a constant pressure boundary of 29 MPa in our calculation. The diameter of sample is 10.5 cm, and the length of sample L is 14.2 cm. The core is divided equally in steps of Δx = 0.02 cm, the time step is Δt = 2 s, and the flow at the outlet of the first grid is gas production Q.
The initial condition of Equation (15) is: $$p={p}_{0}=29\text{\hspace{0.17em}}\text{Mpa},\text{}0\le x\le L,\text{}t=0\text{.}$$(27)
The boundary conditions are given by: $$p=29\text{\hspace{0.17em}}\text{Mpa},\text{}x=0,\text{}t>0\text{,}$$(28) $$p=2\text{\hspace{0.17em}}\text{Mpa},\text{}x=L,0<t\le 37\times 24\times 3600\text{,}$$(29) $$p=0.1\text{\hspace{0.17em}}\text{Mpa},\text{}x=L,t>37\times 24\times 3600\text{.}$$(30)
Using “Matlab” software to calculate Equation (15), we obtain daily gas production and pressure changes over time. Parameters used in the calculation are displayed in Table 2. The diffusion coefficient D is a fitted value. Other parameters are experimental results of parallel sample 2. Permeability k was measured by pulse decay method of core plug, and both porosity Φ and density ρ_{r} were measured by helium expansion method. Sw was measured by DeanStark method. The desorption parameters was determined by supercritical adsorption and desorption experiment results, as is shown in Figure 5. The experiment was conducted at 85 °C with methane as adsorbent. The desorption curve and adsorption curve were different, and the desorption parameters (V_{d}, p_{d} and c) were regressed by desorption Equation (7), which refers to Guo (Guo et al., 2013).
Figure 5 shows the comparison between numerical calculation gas production result and experimental result. It can be seen that the numerical result and experimental result are in a good agreement. The initial gas production is high and declines very fast, while the late gas production is low but declines very slowly. Figure 7 shows pressure propagation as time changes. As is shown, pressure propagation is very slow and pressure propagates to the inlet after 100 days, which corresponds to experimental result. Pressure drop near the outlet is steep and pressure drop near the inlet is smooth.
Parameters of numerical calculation.
Fig. 5 Supercritical adsorption and desorption experiment results of sample 2. 
Fig. 7 Diagram of pressure propagation. 
3.5 Analysis of influence factors
3.5.1 Influence of adsorption gas
The study object is still the sample used in the experiment. When we calculate the influence of different parameters on gas production, the outlet pressure is set to 0.1 MPa. Adsorption gas accounts for a large proportion of shale gas. How adsorption gas affects gas production and pressure propagation is very important for effective development of shale gas reservoirs. Figure 8 shows the calculated results of gas production with and without adsorption gas. As is seen from the figure, gas production decreases very quickly when there is no adsorbed gas. The presence or absence of adsorbed gas has little effect on gas production at the early stage. However, at the later stage, gas production with adsorbed gas is far greater than that without adsorbed gas, and adsorbed gas can maintain gas production at the later stage. Figure 9 shows the comparison of pressure propagation with and without adsorbed gas. According to the figure, pressure propagation process becomes faster when there is no adsorbed gas and the presence of adsorbed gas can delay pressure propagation. Adsorbed gas can replenish gas pressure in pores by desorption. Gas desorption supplies additional gas to pores, thereby reservoir pressure can be maintained. Gas production composition analysis is very important during shale gas exploration. Adsorption gas production and free gas production can be analyzed through our proposed mathematical model. Figure 10 shows the ratio of adsorption gas production to free gas production during the experiment. We can see that the ratio of adsorption gas production to free gas production increases with time, and adsorption gas provides a large proportion of gas production during the later stage. Here it should be mentioned that the porosity of the sample is 0.4% and the free gas in the shale matrix is less compared with adsorption gas. So in first several days of production, the ratio of adsorption gas to free gas reaches 0.58, and in the real gas field production, this value will be smaller, but the trend will be the same.
Desorption constants p_{d} and V_{d} also affect the gas production. Figure 11 shows gas production with different V_{d} (V_{d} = 1.2 × 10^{−3}m^{3}/kg, V_{d} = 2.2 × 10^{−3}m^{3}/kg, V_{d} = 3.2 × 10^{−3}m^{3}/kg) when other parameters in Table 2 remain the same. It can be seen that as V_{d} increases, gas production increases and production decline rate declines slower. Figure 12 shows gas production with different p_{d} (P_{d} = 2.52 MPa, P_{d} = 3.52 MPa, P_{d} = 4.52 MPa) when other parameters in Table 2 remain the same. It shows that gas production increases as p_{d} increases, and p_{d} has less impact on gas production than V_{d}.
Fig. 8 Influence of adsorption gas on production. 
Fig. 9 Influence of adsorption gas on pressure propagation. 
Fig. 10 Ratio of adsorption gas production to free gas production. 
Fig. 11 Influence of V_{d} on gas production. 
Fig. 12 Influence of p_{d} on gas production. 
3.5.2 Influence of permeability and diffusion coefficient
Gas flow in shale matrix is the combination of seepage and diffusion. Permeability and diffusion coefficient will affect gas production. Figure 13 shows the gas production with different permeability (k = 2 × 10^{−25}m^{2}, k = 4 × 10^{−25} m^{2}) when other parameters in Table 2 remain the same. It is shown that gas production increases as permeability increases. Figure 14 is gas production with different diffusion coefficient (D = 3 × 10^{−11} m^{2}/s, D = 3 × 10^{−12}m^{2}/s) when other parameters in Table 2 remain the same. It shows that gas production increases as diffusion coefficient increases.
Fig. 13 Influence of permeability on gas production. 
Fig. 14 Influence of diffusion coefficient on gas production. 
3.5.3 Influence of porosity
Free gas accounts for a large proportion of shale gas and free gas is closely related to porosity. Figure 15 shows the shale gas production with different porosity (Φ = 0.4%, Φ = 2%, Φ = 4%) when other parameters in Table 2 remain the same; we can see that gas production increases as porosity increases. So in the development of shale gas reservoirs, porosity is a very important parameter which is directly related to shale gas production.
Fig. 15 Influence of porosity on gas production. 
4 Conclusion

Physical modeling of gas production in shale matrix result shows that gas production displays high initial production with rapid decline rate and a low late production with a smooth decline rate; pressure propagation in shale matrix is slow. The shape of gas production curve is consistent with actual shale gas production data.

The proposed mathematical model can describe gas flow in shale matrix; the numerical result of mathematical model and experimental result are in a good agreement.

Adsorbed gas will affect gas production and pressure propagation speed. The presence of adsorbed gas can delay pressure propagation. Gas production decreases very quickly when there is no adsorbed gas. The presence or absence of adsorbed gas has little effect on gas production in the early stage. However, at the later stage, gas production with adsorbed gas is far greater than that without adsorbed gas, and adsorbed gas can maintain gas production at the later stage.

Desorption constant, permeability, diffusion coefficient and porosity will affect gas production. The physical parameters (gas desorption capacity, diffusion coefficient, permeability and porosity) need to be considered in shale gas development. The mathematical model provides a more accurate calculating method for shale gas production forecast.
Nomenclature
v^{T}: Total flow velocity in shale matrix, m^{3}/s
v^{D}: Darcy flow velocity, m^{3}/s
v^{K}: Diffusion flow velocity, m^{3}/s
D: Diffusion coefficient, m^{2}/s
C: Mass concentration, kg/m^{3}
ρ_{r}: Shale matrix density, kg/m^{3}
M_{g}: Gas molecular weight, mol/kg
C_{fre}: Free gas mass concentration, kg/m^{3}
C_{ads}: Adsorption gas mass concentration, kg/m^{3}
V_{r}: Shale bulk volume, m^{3}
ρ_{a}: Gas density under standard state, kg/m^{3}
V_{des}: Adsorption gas involves in desorption, m^{3}/kg
S_{w}: Water saturation, dimensionless
V_{d}: Maximum adsorption capacity during desorption process, m^{3}/t
p_{d}: Comprehensive constant, Pa^{−1}
c: Residual adsorption capacity, m^{3}/t
q_{d}: Total adsorption gas involves in desorption, m^{3}
Q: Source or sink term, m^{3}/s
V_{std}: Gas mole volume in standard state, 22.4 × 10^{−3}m^{3}/mol
Acknowledgments
This research was financially supported by The National Municipal Science and Technology Project (No. 2016ZX05037006005).
References
 Albinali A., Holy R., Sarak H., Ozkan E. (2016) Modeling of 1D anomalous diffusion in fractured nanoporous media, Oil Gas Sci.Technol. – Rev. IFP Energies nouvelles, 71, 4, 56. [Google Scholar]
 Alharthy N.S., Al Kobaisi M., Kazemi H., Graves R.M. (2012) Physics and modeling of gas flow in shale reservoirs, SPE Abu Dhabi International Petroleum Conference and Exhibition, Abu Dhabi, UAE, 1114 November. [Google Scholar]
 Beskok A., Karniadakis G.E. (1999) Report: A model for flows in channels, pipes, and ducts at micro and nano scales, Microscale Thermophys. Eng. 3, 1, 43–77. [Google Scholar]
 Cao C., Li T., Shi J., Zhang L., Fu S.X., Wang B.T. (2016) A new approach for measuring the permeability of shale featuring adsorption and ultralow permeability, J. Nat. Gas Sci. Eng. 30, 548–556. [CrossRef] [Google Scholar]
 Chareonsuppanimit P., Mohammad S.A., Robinson R.L., Khaled J., Gasem K.A.M. (2012) Highpressure adsorption of gases on shales: Measurements and modeling, Int. J. Coal Geol. 95, 2, 34–46. [CrossRef] [Google Scholar]
 Civan F. (2010) Effective correlation of apparent gas permeability in tight porous media, Transp. porous media 82, 2, 375–384. [CrossRef] [Google Scholar]
 Civan F., Rai C.S., Sondergeld C.H. (2011) Shalegas permeability and diffusivity inferred by improved formulation of relevant retention and transport mechanisms, Transp. porous media 86, 3, 925–944. [CrossRef] [Google Scholar]
 Civan F., Rai C.S., Sondergeld C.H. (2012) Determining shale permeability to gas by simultaneous analysis of various pressure tests, SPE J. 17, 3, 717–726. [CrossRef] [Google Scholar]
 Curtis J.B. (2002) Fractured shalegas systems, AAPG Bull. 86, 11, 1921–1938. [Google Scholar]
 Deng J., Zhu W., Ma Q. (2014) A new seepage model for shale gas reservoir and productivity analysis of fractured well, Fuel 124, 15, 232–240. [CrossRef] [Google Scholar]
 de Swaan A. (2014) Pressure transients in a fractalcluster model of porous media, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 71 1, 9. [Google Scholar]
 Fathi E., Tinni A., Akkutlu I.Y. (2012) Shale gas correction to Klinkenberg slip theory, SPE Americas Unconventional Resources Conference, Pittsburgh, Pennsylvania, USA, 5–7 June. [Google Scholar]
 Freeman C.M., Moridis G.J., Blasingame T.A. (2011) A numerical study of microscale flow behavior in tight gas and shale gas reservoir systems, Transp. porous med. 90, 1, 253–268. [CrossRef] [Google Scholar]
 Guo W., Xiong W., Gao S.S., Hu Z.M., Liu H.L., Yu R.Z. (2013) Impact of temperature on the isothermal adsorption/desorption characteristics of shale gas, Petrol. Explor. Dev. 40, 4, 481–485. [CrossRef] [Google Scholar]
 Javadpour F. (2009) Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone), J. Can. Petrol. Technol. 48, 8, 16–21. [Google Scholar]
 Javadpour F., Fisher D., Unsworth M. (2007) Nanoscale gas flow in shale gas sediments, J. Can. Petrol. Technol. 46, 10, 55–61. [Google Scholar]
 Liu J., Chen Z., Elsworth D., Qu H., Chen D. (2011) Interactions of multiple processes during CBM extraction: A critical review, Int. J. Coal Geol. 87, 3, 175–189. [CrossRef] [Google Scholar]
 Loucks R.G., Reed R.M., Ruppel S.C., Jarvie D.M. (2009) Morphology, genesis, and distribution of nanometerscale pores in siliceous mudstones of the Mississippian Barnett Shale, J. Sediment. Res. 79, 12, 848–861. [CrossRef] [Google Scholar]
 Lu X.C., Li F.C., Watson A.T. (1995) Adsorption measurements in Devonian shales, Fuel 74, 4, 599–603. [CrossRef] [Google Scholar]
 Ozkan E., Raghavan R.S., Apaydin O.G. (2010) Modeling of fluid transfer from shale matrix to fracture network, SPE Annual Technical Conference and Exhibition, Florence, Italy, 19–22 September. [Google Scholar]
 Pan Z., Ma Y., Connell L.D., Down D.I., Camilleri M. (2015) Measuring anisotropic permeability using a cubic shale sample in a triaxial cell, J. Nat. Gas Sci. Eng. 26, 336–344. [CrossRef] [Google Scholar]
 Passey Q.R., Bohacs K., Esch W.L., Klimentidis P., Sinha S. (2010) From oilprone source rock to gasproducing shale reservoirgeologic and petrophysical characterization of unconventional shale gas reservoirs, SPE International Oil and Gas Conference and Exhibition in China, Beijing, China, 8–10 June. [Google Scholar]
 RomeroSarmiento M.F., Pillot D., Letort G., LamoureuxVar V., Beaumont V., Huc, A.Y., Garcia B. (2015) New rockeval method for characterization of unconventional shale resource systems, Oil Gas Sci.Technol. – Rev. IFP Energies nouvelles 71, 37. [CrossRef] [Google Scholar]
 Shabro V., TorresVerdin C., Sepehrnoori K. (2012) Forecasting gas production in organic shale with the combined numerical simulation of gas diffusion in kerogen, Langmuir desorption from kerogen surfaces, and advection in nanopores, SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 8–10 October. [Google Scholar]
 Swami V., Clarkson C.R., Settari A. (2012) NonDarcy flow in shale nanopores: do we have a final answer? SPE Canadian Unconventional Resources Conference, Calgary, Alberta, Canada, 30 October–1 November. [Google Scholar]
 Villazon M., German G., Civan F., Devegowda D. (2011) Parametric investigation of shale gas production considering nanoscale pore size distribution, formation factor, and nonDarcy flow mechanisms, SPE Annual Technical Conference and Exhibition, Denver, Colorado, USA, 30 October–2 November. [Google Scholar]
All Tables
All Figures
Fig. 1 Experimental apparatus of gas production in shale matrix. 

In the text 
Fig. 2 Liquid nitrogen adsorption and desorption experimental results. 

In the text 
Fig. 3 Poresize distribution of shale samples. 

In the text 
Fig. 4 Inlet and outlet pressure curves. 

In the text 
Fig. 6 Comparison of experimental result and numerical calculation result. 

In the text 
Fig. 5 Supercritical adsorption and desorption experiment results of sample 2. 

In the text 
Fig. 7 Diagram of pressure propagation. 

In the text 
Fig. 8 Influence of adsorption gas on production. 

In the text 
Fig. 9 Influence of adsorption gas on pressure propagation. 

In the text 
Fig. 10 Ratio of adsorption gas production to free gas production. 

In the text 
Fig. 11 Influence of V_{d} on gas production. 

In the text 
Fig. 12 Influence of p_{d} on gas production. 

In the text 
Fig. 13 Influence of permeability on gas production. 

In the text 
Fig. 14 Influence of diffusion coefficient on gas production. 

In the text 
Fig. 15 Influence of porosity on gas production. 

In the text 