Physical and mathematical modeling of gas production in shale matrix

. 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 coef ﬁ cient, permeability and porosity) also need to be considered in shale gas development.


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 A., 2014;Romero-Sarmiento 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 dusty-gas 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 dual-porosity and a triple-porosity finite-difference model for single-phase 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 multi-scale 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.

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.

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

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.

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 adsorption-desorption 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 slit-like 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. 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 5 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.

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 : The seepage term can be described by Darcy's law: where k is permeability of shale matrix, m 2 ; m is viscosity of methane, Pa · s; p is gas pressure, Pa. The diffusion term can be described by Knudsen diffusion equation: where D is diffusion coefficient, m 2 /s; C is mass concentration of gas in shale matrix, kg/m 3 ; r g is gas density, kg/m 3 . Shale gas exists in adsorption gas and free gas state. Take one-dimensional case for an example, Assume shale volume is V r , m 3 , the shale rock density is r 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: 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 r is the shale density, kg/m 3 ; r 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; fis 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: 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: Substitute Equations (2) and (3) into (1), we obtain: The continuity equation for gas flow is: where q d is total adsorption gas that involves in desorption, Q is source or sink term, which means gas production or injection. where V std is gas mole volume in standard state, 22.4 Â 10 À3 m 3 /mol.

Mathematical model solution
Numerical solution method is used to solute Equation (15), which refers to Civan Civan et al., 2012) With time-forward difference, the first term on the left side of Equation (15) can be discretized as: Assume: The second term on the left side of the Equation (15) can be discretized as: The first term on the right side of Equation (15) can be discretized as: The second term on the right side of Equation (15) can be discretized as: where: Substituting Equations (16)-(26) into (15), we can get the finite difference equations with time-forward 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 Dx = 0.02 cm, the time step is Dt = 2 s, and the flow at the outlet of the first grid is gas production Q.
The initial condition of Equation (15) is: The boundary conditions are given by: 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 G. Wei et al.: Oil & Gas Science and Technology -Rev. IFP Energies nouvelles 73, 12 (2018) 5 parameters are experimental results of parallel sample 2. Permeability k was measured by pulse decay method of core plug, and both porosity F and density r r were measured by helium expansion method. Sw was measured by Dean-Stark method. The desorption parameters was determined by supercritical adsorption and desorption experiment results, as is shown in Figure 6. 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.

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 .

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.

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 (F = 0.4%, 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.   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.

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 k Permeability, m 2 m Viscosity, Pa · s p Pressure, Pa D Diffusion coefficient, m 2 /s C Mass concentration, kg/m 3 r g Gas density, kg/m 3 V r Shale volume, m 3 r 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 V g Gas volume, m 3 r a Gas density under standard state, kg/m 3 V des Adsorption gas involves in desorption, m 3 /kg S w Water saturation, dimensionless Ф Porosity, 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 t Time, s