## Abstract

In contrast to conventional gas-bearing rocks, gas shale has extremely low permeability due to its nano-scale pore networks. Organic matter which is dispersed in the shale matrix makes gas flow characteristics more complex. The traditional Darcy’s law is unable to estimate matrix permeability due to the particular flow mechanisms of shale gas. Transport mechanisms and influence factors are studied to describe gas transport in extremely tight shale. Then Lattice Boltzmann simulation is used to establish a way to estimate the matrix permeability numerically.

#### Authors:

##### Xiao-Ling Zhang, Li-Zhi Xiao, Long Guo, Qing-Ming Xie

State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China. Key Laboratory of Shale Gas Exploration, Ministry of Land and Resources, Chongqing Institute of Geology and Mineral Resources, Chongqing 400042, China. Received: 22 January 2014, Published online: 13 January 2015

The results show that net desorption, diffusion, and slip flow are very sensitive to the pore scale. Pore pressure also plays an important role in mass fluxes of gas. Temperature variations only cause small changes in mass fluxes. The Lattice Boltzmann method can be used to study the flow field in the micropore spaces and then provides numerical solutions even in complex pore structure models. Understanding the transport characteristics and establishing a way to estimate potential gas flow is very important to guide shale gas reserve estimation and recovery schemes.

## 1 Introduction

Shale gas, which is hard to exploit, will play an increasingly critical role in China’s natural gas production in the future (Jia et al. 2012; Zou et al. 2010; Chen et al. 2012; Sun et al. 2013). The matrix of gas shale consists of fine clay minerals and organic matter, and pore space reduces largely due to the compaction and cementation of these very fine grains (Wang et al. 2013; Li et al. 2007).

Organic-rich shale is treated as source rock which has potential to produce natural gas. The pore scale of gas shale is some two orders of magnitude smaller than that of conventional sandstone (Bustin et al. 2008; Loucks et al. 2009). The pore throat diameter of shale is also very small, only dozens of times the size of methane molecules (Curtis et al. 2010). Therefore, permeability of gas shale is extremely low, normally only about tens to hundreds of nanodarcy. The micropore structure and flow characteristics are very complex and not well understood (Huang et al. 2012; Li et al. 2013).

The shale gas production mechanism after fracturing is as follows. First, free gas in the formation fracture networks is transported to the wellbore driven by the pressure difference between the borehole and formation pressures. Then free gas in the pore network flows to the fracture network due to the gas concentration difference between the pore network and fracture network, resulting in reduction of pore pressure. Finally, gas adsorbed on the surface of organic matter will desorb to pore spaces to increase the pore pressure. Overall, natural gas moves from the pore network to the fracture network and finally flows to the wellbore to be collected (Zou et al. 2014).

To estimate how much gas can be produced, it is necessary to evaluate the gas transport ability of the formation rock. It is very difficult to estimate matrix permeability with traditional experimental methods and to predict flow ability after fracturing due to the complex flow mechanisms of shale gas. In recent years, a lot of studies have been carried out to study the micro-scale flow characteristics of shale gas, which focus on two aspects as, (1) Laboratory measuring methods consist of the pulse pressure decay method for cores and the gas expansion method for crushed samples. They were developed to overcome long measuring times through applying pulse pressure and reducing the volume of the samples. (2) Numerical simulation methods consist of the molecular dynamics method, finite difference method (Shabro et al. 2009, 2011), finite element method (Roy et al. 2003; Yao et al. 2013), level set method (Prodanovic and Bryant 2007; Prodanovic et al. 2008), and the Lattice Boltzmann method (Tolke et al. 2010; Fathi and Akkutlu 2011; Maier and Bernard 2010).

These methodologies resulted in some research progress which established the theoretical foundation of shale gas flow mechanisms, but quantitative research into adsorption mechanisms and influencing factors for micro-scale gas transport characteristics is still needed. Based on micro-scale gas flow mechanisms, this paper focuses on various possible factors that affect shale gas flow characteristics, and quantitatively analyzes the influence of each factor. Finally, Lattice Boltzmann simulation is used with complex models to estimate the matrix permeability.

## 2 Micro flow mechanisms

Darcy’s law assumes that fluid properties are stable without any physical or chemical reaction with the rock. Only the viscous interaction forces between molecules control the flow without change in the fluid properties (viscosity) and pore geometry (cross-sectional area, length) in the regime of laminar flow. But physical adsorption of methane occurs and gas changes states on the surface of kerogen in organic-rich shale. Shale gas flow is also controlled by molecular diffusion. Due to these inherent characteristics of gas shale such as, (1) widely dispersed organic matter and (2) nano-scale pore system (Zou et al. 2013), shale gas transport mechanisms do not fully follow Darcy’s law. We cannot simply use Darcy’s law to calculate matrix permeability of gas shale because of these special phenomena and problems (Civan 2010; Civan et al. 2011).

Shale gas flow is controlled by flow mechanisms at different scales such as, (1) adsorption and desorption, (2) diffusion, (3) slip flow, and (4) nonlinear Darcy flow. All these flow mechanisms make the prediction of shale gas transport more and more complicated. Three flow mechanisms in nano-scale pores will be discussed in the following sections. First a parameter which distinguishes different flow regimes will be introduced.

### 2.1 Knudsen number

In 1934, Knudsen defined a dimensionless parameter (Knudsen parameter, *K*_{n}) which reflects the degree of collision between gas molecules when rarefied gas flows in narrow channels. Therefore, gas flow regimes at different scales can be determined by the value of *K*_{n} expressed as (Ziarani and Aguilera 2012)

where *l*_{char} is the characteristic value of the geometry through which fluids flow (such as channel height, tube radius, etc.), m; *λ* is the mean free path of gas molecules [namely, the average distance traveled by a moving particle (such as an atom or a molecule) between successive impacts (collisions)], m.

Different flow regimes classified by the value of *K*_{n} are shown in Table 1. The rarefied effect of gas molecules becomes more and more significant with an increase of *K*_{n}. For normal shale gas reservoir temperatures (350–450 K) and pressures (25–28 MPa), the mean free path of methane molecules would be between 0.9 and 4.4 nm. If the pore size is 50 nm diameter, *K*_{n} will be 10^{−3}–10^{−2} then the flow regime is slip flow. So in nano-scale pore spaces the continuum assumption of fluid flow breaks down and the Navier–Stokes (for short, N–S) equation with no-slip boundaries is no longer valid.

**Table 1** Flow regimes of Knudsen number.

### 2.2 Adsorption and desorption

As shown in Fig. 1a, before well drilling there is dynamic equilibrium in pores of organic matter which means the adsorption and desorption rates are equal. But when the pore pressure reduces, this equilibrium breaks down and net desorption item occurs which will play a complementary role in the mass flux of the pores as shown in Fig. 1b (Parker et al. 2009).

**Fig. 1** **a** Dynamic equilibrium, **b** Net desorption.

Adsorption mechanisms are very important for shale gas flow. Physical adsorption of supercritical gases at high pressures is known to have some unusual features compared to subcritical fluid adsorption. This is important in the understanding of many adsorption processes and presents a challenge for fundamental theoretical research. All methods used for obtaining isotherms (gravimetric, volumetric, piezometric, and total desorption methods, column breakthrough methods, closed-loop recycle methods, and isotope exchange methods) measure the Gibbsian excess rather than the absolute amount adsorbed. In order to evaluate the amount of adsorbed gas in shale, it is necessary to know the absolute adsorption isotherm.

where *n*_{ex} is the excess adsorbed amount; *ρ*(*z*) is the density at coordinates *z*.

The main idea in this study is the application of the non-ideal EOS for the corresponding bulk phase to an adsorbed phase, but the parameters of the EOS are modified using the local density and the mean field approximations. The behavior of the adsorption phase obeys the EOS

where *f*^{(a)} and *f* are the fugacities of gas in a pore and in the bulk phase, respectively; *u* is the overall potential energy of the gas–solid interaction; *T* is the temperature, K; *R* is the universal gas constant. This equation is valid for subcritical fluids as well.

The force model consists of surface force and long-range Van der Waals forces, which are simulated by a Shan–Chen model (Shan and Chen 1993). On the basis of this operation, 0.5 nm (approximately the CH_{4} molecule diameter) lattice unit length is selected. When more adsorbed molecules remain on the surface, other molecules are less likely to be attracted. Therefore, surface force is proportional to *θ* which can be written as

where *P* is the pressure of gas, Pa; *P*_{L} is the Langmuir pressure constant, Pa. By fitting a Langmuir isotherm to experimental data, constants *P*_{L} and *A* can be obtained. The bond force between the molecules and the surface is confined to the first layer, and the long-range Van der Waals forces hold other layers. The interactions between molecules are governed by

where *ψ* is a function of apparent density; *ρ* is the density of each lattice unit. For high Knudsen numbers, a combination of bounce back and specular reflection was used to generate the slip effect at the wall. The temperature of the whole gas flow is assumed to be constant. Thus, the pressure can be expressed as

The drop of pore pressure leads to a higher desorption velocity compared to adsorption velocity which causes an added flux in the pores as shown in Fig. 2.

**Fig. 2** Net desorption adds gas in the pore.

### 2.3 Diffusion

From the perspective of molecular dynamics, gas diffusion is actually the result of irregular thermal motion of gas molecules. According to the value of *K*_{n}, diffusion can be divided into Fickian diffusion, Knudsen diffusion, and transitional diffusion. In nano-scale pores with high Knudsen number, Knudsen diffusion is the dominant flow mechanism.

Knudsen diffusion in nanometer pores can be written in the form of a pressure gradient (Roy et al. 2003). In nanometer pores ignoring the viscosity effect, mass flux due to diffusion caused by molecular concentration difference (Javadpour et al. 2007) can be expressed as

The Knudsen diffusion coefficient *D*_{k} can be written as

where *R* is the universal gas constant, J/kmol K; *M* is the molar mass of gas, kg/kmol; *T* is the temperature, K; Δ*P* is the pressure difference, Pa; *R*_{0} is the pore radius, m; *L* is the length of the capillary, m; *D*_{k} is the diffusion coefficient, m^{2}/s. It indicates that the Knudsen diffusion coefficient *D*_{k} is proportional to the pore radius *R*_{0} and the square root of the temperature *T.*

### 2.4 Slip flow

Free gas in the formation fracture network flows to the borehole due to the difference between the borehole and formation pressures.

The N–S equation is the momentum conservation equation which is used to describe the flow of incompressible viscous fluid (continuous flow). For long and straight capillary tubes with a circular cross section, the N–S equation can be simplified as

where *ν* is the flow velocity per unit area, m/s; Δ*p* is pressure difference between the ends of the capillary, Pa; *μ* is the gas viscosity under atmosphere pressure, Pa s; *L* is the length of the pore. The solution of the above nonlinear difference equations is given by

where *C*_{0} is an unknown constant which can be calculated under different boundary conditions. Velocity equations are solved under two different boundary conditions; continuum flow (Fig. 3a) and slip flow (Fig. 3b) (Zhang et al. 2012; Javadpour 2009).

**Fig. 3** Velocity profile. **a** Continuum flow, **b** Slip flow.

For continuum flow, fluid can be treated as a unit. The velocities on the pore wall are zero. Mass flux *J*_{a_no_slip} caused by pressure difference can be expressed as

and absolute permeability is given by

In the case of continuum flow, the mass flux of the fluid depends on the proportional constant *K*, the so-called absolute permeability of the rock, mD. When the parameters such as cross-sectional area, length of the capillary, fluid viscosity, and the pressure difference are fixed, *K* is an inherent characteristic of the porous medium independent of the fluid used in the measurement.

For slip flow in nano-scale pores, the velocities on the pore wall are no longer zero. So the mass flux *J*_{a_slip} caused by pressure difference can be expressed as

And absolute permeability is given by

where *ρ*_{avg} is the average density of the fluid; *R*_{0} is the radius of the pore; *F* is the correction factor; *R* is the universal gas constant (Boltzmann’s constant) = 8.314, J/mol K; *M* is the molar mass of the gas, g/mol; *T* is the temperature, K; *p*_{avg} [=(*P*_{1} + *P*_{2})/2] is the average pressure in the capillary, Pa; *α* is the dimensionless adjustment coefficient of angular momentum (affected by the surface smoothness of the pore walls, gas type, temperature, and pressure, the range is from 0.65 to 0.85).

When the pore scale reduces to nanometers, the slip effect can no longer be ignored and the mass flux due to pressure difference needs to be corrected. The smaller the pore radius is, the larger the value of the correction factor *F* will be. The flow velocity per unit area considering the slip effect will increase compared to the situation of continuum flow alone.