Impact of Local Effects on the Evolution of Unconventional Rock Permeability

Abstract

When gas is extracted from unconventional rock, local equilibrium conditions between matrixes and fractures are destroyed and significant local effects are introduced. Although the interactions between the matrix and fracture have a strong influence on the permeability evolution, they are not understood well. This may be the reason why permeability models in commercial codes do not include the matrix-fracture interactions. In this study, we introduced the local force to define the interactions between the matrix and the fracture and derived a set of partial differential equations to define the full coupling of rock deformation and gas flow both in the matrix and in the fracture systems. The full set of cross-coupling formulations were solved to generate permeability evolution profiles during unconventional gas extraction. The results of this study demonstrate that the contrast between the matrix and fracture properties controls the processes and their evolutions. The primary reason is the gas diffusion from fractures to matrixes. The diffusion changes the force balance, mass exchange and deformation.

1. Introduction

The eastern Ordos basin of China, where shale and coal are rich in organic matter and favorable for gas accumulations, has become one of the most important gas development areas for PetroChina. Unconventional reservoirs within this area have an extremely low intrinsic permeability and low porosities. For most low permeability reservoirs, hydraulic fracturing and horizontal drilling are the key techniques to extract natural gas. Gas production in shale reservoirs is attributed to the conductivity of the matrix and fracture systems [1,2]. However, the long-term gas production from these reservoirs is known to be a function of fluid transport in the shale matrix and strongly influenced by the fluid transport process in the inorganic matrix, kerogen and fractures [3,4,5,6]. Therefore, a better understanding of effective permeability evolution in matrixes and fractures and interactions between them is important to guide the industrial production.

Gas transport in shale reservoirs is a combination of desorption and diffusion within the micropores, and Darcy flow (pressure driven volume flow) within the macro-pores, micro-fractures and fractured network system [7]. Because of the high contrasts of matrix properties and fracture ones, they have different mechanical behaviors that directly affect the permeability evolution. Many scholars established models to investigate the evolution of the permeability of different parts over the past few decades. The models were developed from a single porosity/permeability model [8] to dual porosity/permeability model [9,10,11,12]. For a single porosity model, the effects of effective stress variation and the matrix shrinkage/swelling were taken into consideration [13,14], which ignored the interactions between different range radius pores. Then researchers paid attention to the effects of adsorption-induced strain [15,16]. Based on the poroelasticity theory, Zhang et al. [15] developed a strain-based porosity model and a permeability model under variable stress conditions. These models include the coupling interaction between gas flow/diffusion and rock mechanic behavior.

In previous studies, many scholars ignored the dynamic behaviors of matrix and fracture properties, especially the interactions between them. When gas is extracted from the reservoir, the gas pressure in the inorganic matrix, kerogen and microfracture will decrease to a lower magnitude. The effective stress within them will change. The variation of effective stress will affect the pore radius of the matrix and kerogen/micro-fracture, which means that the intrinsic permeability is a variable [17,18]. Chen et al. [19] and Masoudian [20] studied the impact of effective stress and/or strain on permeability in shale fractures. However, the impact of shale matrix effective stress and/or strain on fracture permeability was not considered. Some models considered the variation of stress and just focused on gas flow principles. Cao et al. [2] and Peng et al. [21] developed a model considering the deformation induced by the changes in effective stress. In these studies, only the matrix mechanical deformation is taken into account and the mechanical interactions between the matrix and kerogen/micro-fracture induced by the differential pressure are ignored. However, matrix and micro-fracture properties are different, and their mechanical behaviors are also different under the same loading conditions. In this study, we defined mechanical equilibrium equations for the inorganic matrix and the kerogen/micro-fracture, respectively, to control the shale deformation based on our previous study [18]. Through the full coupling of two solid deformation systems and two gas flow systems in them, we studied the impact of local transient behaviors on the evolution of rock permeability.

2. Conceptual Model

In this section, an overlapping approach is introduced to analyze the full shale-gas interactions. Shale is multi-pore media including micropores, macro-pores, micro-fractures and fractured network systems. When gas is extracted from or injected into shale, local equilibrium conditions between matrixes and fractures are disturbed and significant local effects are introduced into the porous medium. Because of the high contrasts of matrix properties and fracture ones, it may take a much longer time to reach a new state of equilibrium. A schematic diagram of the conceptual overlapping approach is shown as Figure 1. For a specific domain, as shown in Figure 1a, it was meshed, generating many nodes, as shown in Figure 1b and each mathematical node can be overlapped by four physical nodes corresponding to four different physical fields representing the solid deformation in the fracture, solid deformation in the matrix, gas flow in the fracture and gas flow in the matrix, respectively, shown as Figure 1c. The physical behaviors of four physical points are not isolated but connected to each other by a set of coupling relations as shown in Figure 1d.

To explain the interactions between the matrix and fracture, we use gas injection as an example. When gas is injected into the shale as shown in Figure 1a, gas flows into micropores and fractures quickly while there is no gas in the matrix. Subsequently, the pressure in the matrix increases gradually as the gas diffuses from the fractures into the matrixes. The non-synchronization of the pressure change between the matrix and fracture generates a local force which can cause interactions. Specifically, because of the increase of the pressure in the fracture, the effective stress in the fractures firstly decrease, which results in the swelling of the fracture. At the same time, the matrix localized in the vicinity of the fracture compartment shrinks under an extra local force applied by the fracture. The swelling and shrinkage, especially induced by a local force, directly affects the permeability evolution of the matrix and fracture. When the pressure in the fracture reaches a specific magnitude, the gas diffuses into the matrix. The differential pressure starts to decline slowly to zero at a final equilibrium state. In the same way, the matrix swells because of the decrease in effective stress. Finally, the fracture and matrix both swell and the matrix changes from local swelling to macro swelling [22]. Additionally, the permeability of the matrix and fracture will have a net increase.

3. Governing Equations

In this section, a set of partial differential equations are defined. These equations govern the deformation of the matrix and fracture deformation, and control the transport of gas flow. Originally, we develop the governing equations of single pore media based on previous studies [15] which describe the interactions in the two kinds of solid media. These derivations are based on the following assumptions:

• Shale is a homogeneous, isotropic, dual poroelastic continuum.
• Strains are much smaller than the length scale.
• Gas contained within the pores is ideal, and its viscosity is constant under isothermal conditions.
• Gas flow through the shale fracture is defined by Darcy’s law and defined by Knudsen diffusion in the matrix.

Shale contains a matrix and fracture which have different mechanical properties and interact with each other. Therefore, we derived mechanical equations for the matrix and fracture, respectively. They are fully coupled with local force. Darcy’s law is used for both the flow in the matrix and the flow in the fracture.

3.1. Formulation of Solid Deformation

The linear constitutive relations can be obtained by extending the known poroelasticity [23]. In previous studies, the gas sorption-induced strain is assumed to result in volumetric strain only [2,15]. However, in our new model, the volumetric strain includes both the gas sorption strain and local strain. Additionally, their effects on all three normal components of strain are the same. By making an analogy between thermal contraction and matrix shrinkage, the constitutive relations for the deformed shale matrix and fracture can be defined as

$${\epsilon }_{mij}=\frac{1}{2G_m}{\sigma }_{ij}-\left(\frac{1}{6G_m}-\frac{1}{9K_m}\right){\sigma }_{kk}{\delta }_{ij}+\frac{\alpha }{3K_m}{P}_m{\delta }_{ij}-\frac{1}{3K_m}\Delta P{\delta }_{ij}+\frac{{\epsilon }_{ms}}{3}{\delta }_{ij}$$

(1)

$${\epsilon }_{fij}=\frac{1}{2G_f}{\sigma }_{ij}-\left(\frac{1}{6G_f}-\frac{1}{9K_f}\right){\sigma }_{kk}{\delta }_{ij}+\frac{\beta }{3K_f}{P}_f{\delta }_{ij}-\frac{1}{3K_f}\Delta P^{\prime }{\delta }_{ij}+\frac{{\epsilon }_{fs}}{3}{\delta }_{ij}$$

(2)

where $G_m=E_m/2\left(1+{\upsilon }_m\right)$ is the shear modulus of matrix, $E_m$ and ${\nu }_m$ are the Young’s modulus values of shale matrix and the Poisson’s ratio of the matrix, respectively; $G_f=E_f/2\left(1+{\upsilon }_f\right)$ is the shear modulus of the fracture, $E_f$ and ${\nu }_f$ are the Young’s modulus values of the shale fracture and the Poisson’s ratio of fracture, respectively; $K_m=E_m/3\left(1-2{\upsilon }_m\right)$ is the bulk modulus of matrix, $K_f=E_f/3\left(1-2{\upsilon }_f\right)$ is the bulk modulus of fracture; $\alpha$ and $\beta$ are the Biot coefficients; P_m is the fluid pressure in the matrix, P_f is the fluid pressure in the fracture; ${\epsilon }_{ms}$ is the gas sorption-induced strain in the matrix, ${\epsilon }_{fs}$ is the gas sorption-induced strain in the fracture; ${\sigma }_{kk}$ is the total stress; ${\delta }_{ij}$ is the Kronecker delta; $\Delta P$ is the differential pressure between the fracture and the matrix. In addition, $\Delta P=P_m-P_f$ for the constitutive relation of the matrix while $\Delta P^{\prime }=P_f-P_m$ for the constitutive relation of the fracture.

Applying Langmuir isotherm, the sorption-induced volumetric strains of matrix and fracture can be defined as [15,18]

$${\epsilon }_{ms}={\epsilon }_L\frac{P_m}{P_L+P_m}$$

(3)

$${\epsilon }_{fs}={\epsilon }_L\frac{P_f}{P_L+P_f}$$

(4)

where ${\epsilon }_L$ is the Langmuir strain constant and P_L is the Langmuir pressure. When gas flows from the matrix into the fracture, the local deformations of the matrix and fracture are controlled by [18]

$$\sigma +\Delta P={\sigma }_{\mathit{me}}+{\alpha P}_m$$

(5)

$$\sigma +\Delta P^{\prime }={\sigma }_{fe}+\beta P_f$$

(6)

where $\Delta P=P_f-P_m$, $\Delta P^{\prime }=P_m-P_f$ are local forces induced by differential pressures between the fracture and matrix systems; $\sigma +\Delta P$ and $\sigma +\Delta P^{\prime }$ are the dynamic effective stress values. ${\sigma }_{me}$ and ${\sigma }_{fe}$ are the effective stress component of the matrix and the fracture, respectively; $\alpha$ and $\beta$ are the Biot coefficients of the matrix and the fracture. From Equations (1)–(6), the volumetric strains of matrix and fracture can be expressed as

$${\epsilon }_{mv}=-\frac{1}{K_m}\left(\overline{\sigma }-\alpha P_m+\Delta P\right)+{\epsilon }_{ms}$$

(7)

$${\epsilon }_{fv}=-\frac{1}{K_f}\left(\overline{\sigma }-\beta P_f+\Delta P^{\prime }\right)+{\epsilon }_{fs}$$

(8)

where $\overline{\sigma }=-{\sigma }_{kk}/3$ is the mean compressive stress. Combining Equations (1), (2), (7) and (8) yields the general Navier-type equations for the matrix and fracture, respectively

$$G_m{u}_{i,kk}+\frac{2G_m}{1-2\nu }{u}_{k,ki}-\alpha P_{m,i}+\Delta P_{,i}-K_m{\epsilon }_{ms,i}+F_i=0$$

(9)

$$G_f{u}_{i,kk}+\frac{2G_f}{1-2\nu }{u}_{k,ki}-\beta P_{f,i}+\Delta P_{,i}^{\prime }-K_f{\epsilon }_{fs,i}+F_i=0$$

(10)

Equation (9) and Equation (10) are the governing equations for the shale matrix and fracture deformation. They are cross-coupled by the local force, $\Delta P$, which reflects the mechanical interaction between the matrix and the fracture.

3.2. Formulation of Gas Flow in the Fracture

The gas flow within the natural fractures obeys Darcy’s law. The equation for the mass balance of the gas is defined as

$$\frac{\partial m_f}{\partial t}+\nabla ·\left(-\frac{k_f}{\mu }{\rho }_{gf}\nabla P_f\right)=-Q_{mf}$$

(11)

where $\mu$ is the gas dynamic viscosity, $m_f={\varphi }_f{\rho }_{gf}+{\rho }_g{\rho }_c\frac{V_L{P}_f}{P_f+P_L}$ is the gas content in the fracture including the free-phase gas and adsorbed gas, ${\varphi }_f$ is fracture porosity, ${\rho }_g$ is the gas density at standard conditions, ${\rho }_{gf}=\frac{M_g}{RT}{P}_f$ is the gas density, $k_f$ is the permeability of the fractures and $-Q_{mf}$ is mass the transfer from the matrix to the fractures.

3.3. Formulation of Gas Flow in the Matrix

Gas flow in the matrix follows Darcy’s law, so the equation for the mass transfer of the gas in the matrix is defined as

$$\frac{\partial m_m}{\partial t}+\nabla ·\left(-\frac{k_m}{\mu }{\rho }_{gm}\nabla P_m\right)=Q_{mf}$$

(12)

where $m_m={\varphi }_m{\rho }_{gm}+{\rho }_{ga}{\rho }_c\frac{V_L{P}_m}{P_m+P_L}$ is the gas content in the matrix including free-phase gas and adsorbed gas, $k_m$ is the permeability of the matrix, $\mu$ is the dynamic viscosity of the gas, ${\rho }_g$ is the gas density at standard conditions, ${\rho }_{gm}=\frac{M_g}{RT}{P}_m$ is the gas density in the matrix ($M_g$ is the molecular mass of the gas, R is the universal gas constant, T is the absolute gas temperature), and $Q_{mf}$ is the gas mass transfer from the fracture to the matrix.

3.4. Formulation of Cross-Couplings

The mechanisms of mass transfer for a dual porosity media are fluid expansion and viscous displacement. The final form of the transfer function for a single-phase flow from the matrix to the fracture is given as [24]

$$Q_{mf}=aV{\rho }_g\frac{k_m}{\mu }\left(P_m-P_f\right)$$

(13)

where ${\rho }_g$ is the density of the gas, $k_m$ is the permeability of the matrix, $\mu$ is the viscosity, a is called the matrix-fracture transfer shape factor and has dimensions of L⁻² that equal 1, $P_m$ is the matrix pressure, and $P_f$ is the fracture pressure.

We derived the general permeability model of the shale matrix and fracture. Shale rock contains a fracture system with well-connecting macropores and a matrix system with micropores. For each system, considering it contains a solid volume $V_s$ and pore volume $V_p$, the shale bulk volume can be defined as $V=V_p+V_s$ and the porosity can be defined as $\varphi =V_p/V$. According to Equation (5) and (6), the volumetric evolution of the porous medium loaded by $\overline{\sigma }$, $p_m$ or $p_f$ and $\Delta p=p_f-p_m$ can be described in terms of $\Delta V/V$ and $\Delta V_p/V_p$, the volumetric strain of the shale matrix/fracture and the volumetric strain of the pore space, respectively [23]. The relations are

$$\frac{\Delta V_i}{V_i}=-\frac{1}{K_i}\Delta \overline{\sigma }+\left(\frac{1}{K_i}-\frac{1}{K_{is}}\right)\Delta P_i+\frac{1}{K_i}\Delta P+\Delta {\epsilon }_s$$

(14)

$$\frac{\Delta V_{ip}}{V_{ip}}=-\frac{1}{K_{ip}}\Delta \overline{\sigma }+\left(\frac{1}{K_{ip}}-\frac{1}{K_{is}}\right)\Delta P_i+\frac{1}{K_{ip}}\Delta P+\Delta {\epsilon }_s$$

(15)

$$\Delta P=P_f-P_m$$

(16)

where subscript i = 1 and 2 represent the shale fracture and matrix, respectively. If we apply ${\alpha }_i=1-K_i/K_{is}$ and ${\beta }_i=1-K_{ip}/K_{is}$, then the above equations can be expressed as

$$\frac{\Delta V_i}{V_i}=-\frac{1}{K_i}\left(\Delta \overline{\sigma }-{\alpha }_i\Delta P_i-\Delta P\right)+\Delta {\epsilon }_s$$

(17)

$$\frac{\Delta V_{ip}}{V_{ip}}=-\frac{1}{K_{ip}}\left(\Delta \overline{\sigma }-{\beta }_i\Delta P_i-\Delta P\right)+\Delta {\epsilon }_s$$

(18)

where $K_p$ and $K_s$ are the bulk moduli of the pore and the bulk modulus of the solid. We assume that the sorption-induced strain for shale is the same as for the pore space. Applying the definition of porosity, the following expressions can be defined as [15]

$$\frac{\Delta V_i}{V_i}=\frac{\Delta V_{is}}{V_i}+\frac{\Delta {\varphi }_i}{1-{\varphi }_i}$$

(19)

$$\frac{\Delta V_{ip}}{V_{ip}}=\frac{\Delta V_{is}}{V_i}+\frac{\Delta {\varphi }_i}{{\varphi }_i\left(1-{\varphi }_i\right)}$$

(20)

By solving Equations (16)–(20), we can obtain the relationship as

$$\Delta {\varphi }_i={\varphi }_i\left(\frac{1}{K_i}-\frac{1}{K_{ip}}\right)\left(\Delta \overline{\sigma }-\Delta P_i-\Delta P\right)$$

(21)

Substituting $K_{ip}={\varphi }_i{K}_i/{\alpha }_i$ [25] into the above equation yields

$${\varphi }_i-{\varphi }_{i0}={\varphi }_i\left(1-\frac{{\alpha }_i}{{\varphi }_i}\right)\frac{\Delta \overline{\sigma }-\Delta P_i-\Delta P}{K_i}$$

(22)

Rearranging Equation (22) gives

$${\varphi }_i=\frac{{\varphi }_{i0}}{1-\frac{\Delta \overline{\sigma }-\Delta P_i-\Delta P}{K_i}}-\frac{{\alpha }_i}{1-\frac{\Delta \overline{\sigma }-\Delta P_i-\Delta P}{K_i}}\frac{\Delta \overline{\sigma }-\Delta P_i-\Delta P}{K_i}$$

(23)

Because generally $\left(\Delta \overline{\sigma }-\Delta P_i-\Delta P\right)/K_i\ll 1$, the above equation can be simplified into

$$\frac{{\varphi }_i}{{\varphi }_{i0}}=1-\frac{{\alpha }_i}{{\varphi }_{i0}}\frac{\Delta \overline{\sigma }-\Delta P_i-\Delta P}{K_i}=1+\frac{{\alpha }_i}{{\varphi }_{i0}}\left(\Delta {\epsilon }_{iet}-\Delta {\epsilon }_{il}\right)$$

(24)

where $\Delta {\epsilon }_{iet}=-\left(\Delta \overline{\sigma }-\Delta P\right)/K_i$ is defined as the total effective volumetric compressive strain, and $\Delta {\epsilon }_{il}=-\Delta P/K_i$ is defined as the local strain induced by differential pressures between the two systems.

The typical relationship between porosity and permeability follows the cubic law [26]

$$\frac{k}{k_0}={\left(\frac{\varphi }{{\varphi }_0}\right)}^3$$

(25)

By substituting Equation (24) into Equation (25) we obtain shale permeability model

$$\frac{k_i}{k_{i0}}={\left(1-\frac{{\alpha }_i}{{\varphi }_{i0}}\frac{\Delta \overline{\sigma }-\Delta P_i-\Delta P}{K_i}\right)}^3={\left[1+\frac{{\alpha }_i}{{\varphi }_{i0}}\left(\Delta {\epsilon }_{iet}-\Delta {\epsilon }_{il}\right)\right]}^3$$

(26)

The total effective volumetric strain can be written as

$$\Delta {\epsilon }_{iv}=\Delta {\epsilon }_{iet}-\Delta {\epsilon }_{is}+c_l\Delta {\epsilon }_{\left(i+1\right)s}$$

(27)

where $\Delta {\epsilon }_{iv}$ is the volumetric strain, $\Delta {\epsilon }_{is}$ is the volumetric strain induced by sorption, c_l is the local strain coefficient.

By substituting Equation (27) into Equation (26), we can obtain the permeability model for the shale matrix and fracture.

$$\frac{k_i}{k_{i0}}={\left[1+\frac{{\alpha }_i}{{\varphi }_{i0}}\left(\Delta {\epsilon }_{iv}+\Delta {\epsilon }_{is}-\Delta {\epsilon }_{il}-c_l\Delta {\epsilon }_{\left(i+1\right)s}\right)\right]}^3$$

(28-a)

$$\frac{k_m}{k_{m0}}={\left[1+\frac{{\alpha }_m}{{\varphi }_{m0}}\left(\Delta {\epsilon }_{mv}+\Delta {\epsilon }_{ms}-\frac{p_f-p_m}{K_m}-c_l\frac{K_f}{K_m}\Delta {\epsilon }_{fs}\right)\right]}^3$$

(28-b)

$$\frac{k_f}{k_{f0}}={\left[1+\frac{{\alpha }_f}{{\varphi }_{f0}}\left(\Delta {\epsilon }_{fv}+\Delta {\epsilon }_{fs}-\frac{p_m-p_f}{K_f}-c_l\frac{K_m}{K_f}\Delta {\epsilon }_{ms}\right)\right]}^3$$

(28-c)

$$c_l=\eta \frac{P_f-P_m}{\Delta P_{max}}$$

(28-d)

where subscript i = 1 and 2 represent the shale fracture and matrix, respectively. c_l is the local strain coefficient which is in proportion to differential pressures and $\left(P_f-P_m\right)/\Delta P_{max}$, $\eta$ is a constant.

Shale has many nanopores. The flow regimes of the gas also strongly affect the apparent permeability [26,27]. The relation between apparent permeability, $k_{mapp}$, and intrinsic permeability of matrix, $k_{m0}$, is

$$k_{mapp}=k_{m0}g\left(K_n\right)$$

(29)

Kn is the Knudsen number and can be expressed as

$$g\left(K_n\right)=\left(1+\zeta K_n\right)\left(1+\frac{4K_n}{1+K_n}\right)$$

(30)

$\zeta$ is a dimensionless rarefaction coefficient. Its value varies: $0<\zeta <{\zeta }_0$ for $0<Kn<\infty$. ${\zeta }_0$ is an empirical parameter and the dimensionless rarefaction correlation is presented by Civan et al. [28]

$$\zeta =\frac{{\zeta }_0}{1+\frac{A}{K_n^B}}$$

(31)

where A = 0.17, B = 0.4348, and ${\zeta }_0$ = 1.358.

$K_n$ is defined as the ratio of the molecular mean free path, $\lambda \left(nm\right)$ and pore radius $r\left(nm\right)$.

$$K_n=\frac{\lambda }{r}$$

(32)

The mean-free-path of molecules $\lambda$ is given by [24]

$$\lambda =\frac{K_{B}T}{\sqrt{2}\pi {\tilde{d}}^2{p}_m}$$

(33)

where $K_B$ is Boltzmann constant, T is the temperature of shale reservoir, $\tilde{d}$ is the collision diameter for molecules. Based on the research of Wei et al. [29], The nanopore radius can be obtained as

$$r=r_0-\Delta r$$

(34)

where r is the average nanopore radius, r₀ is the initial nanopore radius, $\Delta r$ is the thickness of the adsorbed layer. The average thickness of the adsorbed layer can be expressed as [29]

$$\Delta r=t_a\exp \left(-D{\left[ln\left(\frac{{\rho }_{van}}{{\rho }_g}\right)\right]}^2\right)$$

(35)

where D is a constant and equals to 0.07 [29,30], and t_a is the thickness of the adsorbed layer at extremely high pressures, ${\rho }_g$ is the density of the gas at the specific temperature and pressure, ${\rho }_{van}$ is the gas density of the adsorbed phase (generally assumed to be the van der Waals density of the gas) and is 370 kg/m³. Substituting Equation (37) into Equation (36) yields

$$r=r_0-t_a\exp \left(-D{\left[ln\left(\frac{{\rho }_g}{{\rho }_G}\right)\right]}^2\right)$$

(36)

Therefore, the Knudsen number becomes

$$K_n=\frac{K_{B}T}{\sqrt{2}\pi {\tilde{d}}^2{p}_m}\frac{1}{r_0-t_a\exp \left(-D{\left[ln\left(\frac{{\rho }_g}{{\rho }_G}\right)\right]}^2\right)}$$

(37)

Therefore, the final formulation of apparent permeability of inorganic matrix can be expressed as

$$k_{mapp}=k_{m0}{\left[1+\frac{\alpha }{{\varphi }_{m0}}\left(\Delta {\epsilon }_{mv}+\Delta {\epsilon }_{mI}-\Delta {\epsilon }_{ms}\right)\right]}^3\left(1+\zeta K_n\right)\left(1+\frac{4K_n}{1+K_n}\right)$$

(38)

4. Evolution of Shale Permeability under Stress-Controlled Conditions

The above complete set of formulations, four field equations (Equations (9)–(12)) and two permeability models (Equation (28-a)–(28-c), Equation (38)), are implemented into COMSOL MULTIPHYSICS, a commercial PDE solver. The model geometry of 5 cm × 10 cm is shown in Figure 2. The new model considers mechanical deformations, sorption-induced volumetric strain, local strain induced by local force and interactions between two systems. The simulations were conducted under the constant confining stress condition, 15 MPa. Methane gas (absorbing gas) was used in the simulation and the pore pressure increased from 4 MPa to 8 MPa. The extended material properties are shown in

Table 1. The property parameters of the simulation sample.

Parameter	Value	Physical Meaning	Units
Em	10	Young’s modulus of the matrix	GPa
Ef	2	Young’s modulus of the fracture	GPa
vm	0.35	Poisson’s ratio of the matrix	-
vf	0.2	Poisson’s ratio of the fracture	-
α	0.8	Biot coefficient of the fracture	-
β	0.4	Biot coefficient of the matrix	-
μ	1.11×10−5	Viscosity of Methane	Pa·s
φm0	0.08	Initial matrix porosity	-
φf0	0.04	Initial fracture porosity	-
km0	0.5×10−20	Initial matrix permeability	m²
kf0	1×10−19	Initial fracture permeability	m²
PL	6.109	Langmuir pressure constant	MPa
εL	0.02	Langmuir volumetric strain constant	-
$\rho$m	1250	Matrix density	kg/m³
$\rho$f	1000	Fracture density	kg/m³
Pa	0.1	Atmosphere pressure	MPa
$\rho$g	0.178	Density of gas at standard condition	kg/m³
M	0.016	Molar mass of methane	kg/mol
R	8.314	Gas constant	J/(mol·K)
T	298.15	Temperature of the reservoir	K

.

From Figure 3, it can be seen that the different mechanical properties of the shale matrix and fracture play an important role in controlling the gas flow process. Pore pressure in the fracture equals that in the matrix at the initial equilibrium station. When the gas is injected into the subject, the pore pressure in the fracture increases firstly because of the higher intrinsic permeability while the pore pressure in the matrix keeps stable because of extremely low intrinsic permeability. More interestingly, the asynchronization of gas flow in the matrix and fracture generates a differential pressure which will induce local strain that has an effect on permeability. When the gas diffuses into the matrix, the pore pressure in the matrix starts to increase gradually. Until the gas diffuses into the whole shale, the subject reaches the equilibrium station again and the differential pressure becomes zero again too.

A typical permeability evolution during gas injection is illustrated in Figure 4. It consists of two main features: the permeability ratio at point D is higher than it at the initial station; there is a wave crest B and a trough of wave C. During gas injection, the shale permeability ratio evolution can be divided into three stages. At the first stage (SI), the shale permeability ratio increases by around 20% in the short term. This is because the pore pressure in the fracture increases dramatically while that in the matrix is still the initial value. Therefore, the differential pressure will generate a local compress strain on the matrix located in the vicinity of the fracture which causes the swell of the fracture. In addition, the gas-sorption induced strain also results in the improvement of the permeability. Prior to the diffusion in the matrix, the strain induced by differential pressures and sorption both reach the maximum and the permeability reaches a crest as well. At the second stage (SII), the shale permeability ratio switches from an increase to a decrease to the initial value and then continues to decline to the minimum at point C. The pore pressure in the matrix starts to increase as the gas diffuses into the matrix. The differential pressure declines gradually. The sorption-induced swelling strain has a negative impact on the fracture and this effect will be accumulated with the expanding of the local strain region. Because the porosity of the matrix is higher than the fracture’s, the accumulated sorption-induced strain of the matrix has a significant influence on the fracture permeability and this process lasts for a relatively long time. At the third stage (SIII), the permeability ratio recovered and increases and after 30 days, it finally increased by 30%. This phenomenon is contributed by two reasons: (1) the effective strain of shale changed from a local strain to a global strain when the gas spreads uniformly within the shale; (2) the effective stress decreases due to the increase of the pore pressure under the constant confining stress condition. Therefore, the permeability ratio has a net growth at the end equilibrium station.

4.1. Impact of Local Strain on Permeability

In this study, the local strain includes two components: a local strain induced by differential pressure; the local strain induced by an interaction of the sorption strain. In order to investigate the impacts of local strain on permeability, we conducted four scenarios by considering the mechanism of local strains induced by the differential pressure and the interaction of the sorption-induced strain, both of them and none of them, respectively. The results of the resultant permeability ratio of shale are shown in Figure 5.

From Figure 5, it can be seen that the local strains play an extremely important role in the shale permeability evolution. The first scenario represents the permeability ratio profile without the local strain effects. In this case, the permeability increases monotonically and is controlled by the effective strain which follows the principle of effective stress. The second scenario represents the permeability profile with the impact of local strain induced by differential pressure only. Based on the above analysis in Section 2, the local strain increases the fracture aperture. The reason is that the differential pressure generated a compress strain in the matrix located in the vicinity of the fracture. At the end equilibrium station, the local strain transforms into a global strain with the disappears of differential pressure. Therefore, the permeability declines by a small portion. The third scenario represents the permeability profile with the impacts of local strain induced by an interaction of the sorption strain. The fourth scenario represents the permeability profile with the impact of local strain induced by differential pressure and the interaction of sorption strain. From the third and fourth scenarios, we can know that with the impact of the local strain, the resultant permeability ratio increases over 1 at the initial stage and then declines dramatically under 1. With the transformation of the local strain to the global strain, the local effect finally disappears. Therefore, the resultant permeability ratio rebounds to a value over 1.

4.2. Impact of Modulus Ratios on Permeability

In order to investigate the influence of shale mechanical properties on permeability, a simulation case was conducted with different bulk modulus ratios (K_m/K_f) and the same injection pressure and confining pressure. Results corresponding to three cases (K_m/K_f = 4, 6 and 10) are shown in Figure 6. Permeability ratios for all the cases follow the same pattern. However, when the difference of the matrix and the fracture modulus are small like K_m/K_f = 4, there is a net increase of the permeability. When the modulus ratio is large enough such as K_m/K_f = 6 and 10, there is a significant decrease in the permeability ratio. Additionally, the higher the bulk modulus ratio is, the more the permeability ratio generated decreases. These results show that the difference of the mechanical properties between the matrix and the fracture is positively related to the impact of the local strain on permeability. These phenomena may reveal that there are more fracture or macropores in the shale with large different bulk modulus ratios and the more fractures or macropores, the more significant the local strain effects.

4.3. Impact of Pore Pressure on Permeability

Figure 7 presents the evolution of shale permeability with different injection pressures. The crest and trough of the permeability ratio are higher and lower when the injection pressure is higher. Figure 7b is shale permeability ratio with pore pressure and time. It is a three dimensions graph. From this figure, we can know that the permeability ratio profile is not only related to the pore pressure but also to the time the value is obtained at.

4.4. Impact of Klinkenberg Effects on Permeability

The impact of the Klinkenberg effect (Slip effect) on matrix permeability and the resultant permeability of shale is shown in Figure 8. The parameters used in this simulation case are collected from the literature [31] and are listed in

Table 2. The parameters of the apparent permeability model of the shale matrix.

Symbol	Value	Physical Meanings	Units
KB	1.38 × 10−23	Boltzmann constant	J/K
T	298.15	Temperature	K
τh	1	Tortuosity of the matrix	-
A	0.178	First constant for ζ	-
B	0.4348	Second constant for ζ	-
ζ0	0.25	Asymptotic upper limit of ζ	-

. The Klinkenberg effect (slip effect) is closely related to the magnitude of the pore pressure and it has a significant influence on the matrix permeability especially under low pore pressures. From the figure, we can know that the higher the pore pressure, the lower the matrix permeability. However, compared to the impact of local strain, the Klinkenberg effect has a weak influence on the resultant permeability. In addition, the Klinkenberg effect cannot explain the net increase of the resultant permeability under the constant confining stress condition.

4.5. Model Evaluation and Discussions

In this section, simulations are conducted using our new model to illustrate the impact of local strains on the permeability evolution. We collect experimental data from Xiang Li at al. [32]. A series of experiments were conducted on the Green River Shale under the condition of a constant total confining stress of 20 MPa. Different gas, He, CH₄ and CO₂, were injected into the specimens with artificial fractures. The values of the mechanical parameters such as Poisson’s ratio and Young’s modulus are assumed based on the literature [33,34] and are listed in

Table 3. The mechanical parameters of the shale samples.

Sample	Mechanical Properties
	Em (GPa)	Ef (GPa)	${\mathit{\nu }}_{\mathit{m}}$	${\mathit{\nu }}_{\mathit{f}}$
Green River Shale	8	6	0.25	0.25

. The Langmuir constants of shale and dynamic viscosity of gases are listed in

Table 4. The Langmuir constants of shale and dynamic viscosity of gases at 300 K.

Gas	${\mathit{\epsilon }}_{\mathit{L}}$	PL (MPa)	$\mathit{\mu } (\mathit{\mu }\mathit{P}\mathit{a}\cdot \mathit{s})$	$\mathit{\rho }$ (kg/m3)
He	-	-	18.9	1.293
CO2	0.0353	3.82	14.932	1.784
CH4	0.0093	6.1	11.067	0.648

.

To investigate the impact of non-sorbing gas on permeability, researchers obtained the permeability data at pore pressures of 2 MPa, 4.2 MPa, 6.16 MPa, 8.16 MPa and 10.1 MPa. We conducted simulations with the same pore pressures under the same conditions using our new model. Figure 9 shows the permeability evolution against time and pore pressure and the comparisons between the experimental data and simulation results using non-sorbing gas (He). From the experiment, we can know that the shale permeability increases with pore pressure. From the simulation results, we can obtain a three-dimensional permeability evolution not only with pore pressure but also with time. All the experimental data are located in a special zone between 2 h and 7 h.

To investigate the impact of the local effect of sorption-induced strain on permeability evolution. Sorbing gases (CH₄ and CO₂) were used to conduct experiments. The permeability for sorbing gases CH₄ and CO₂ all show the typical U-shaped curve, as shown in Figure 10 and Figure 11. The permeability decreases firstly and then rebounds to an increase. Compared to the simulation results, the experimental data are in a certain zone with a diffusion time from 4.5 h to 7 h and from 9 h to 17.5 h for CH₄ and CO₂, respectively. All these experimental data are obtained in the process from the local strain to global strain.

It is obvious that the result calculated by the model with the impact of the local effects matches the experimental data well. And the model without the impact of the local effects induces significant errors that causes the apparent permeability to decrease slowly. The reason is that the local effects applied press stress on the matrix that caused an increase in the effective stress. The change of the effective stress decreased the permeability of the matrix. This phenomenon clearly illustrates the importance of local strain due to local force for the apparent permeability evolution of shale.

5. Conclusions

The full shale matrix-fracture interactions such as mass exchange and deformation compatibility are included into a fully coupled shale deformation and gas flow model. Based on the results of this study, the following conclusions can be drawn:

For shale, gas injection induced effects can last very long because of the huge contrast between the matrix and fracture properties. The primary reason for those added effects is the gas diffusion from the fractures to the matrixes. The diffusion changes the force balance, mass exchange and deformation. Therefore, the non-equilibrium processes are much more important than the equilibrium ones.

The shale matrix and fracture’s permeability experience three stages during gas injection: the initial stage of the fracture permeability increase and the matrix permeability decrease as the gas pressure in the fracture increases, the intermediate stage of the fracture permeability decrease and the matrix permeability increase as the gas diffuses from the fractures into the matrixes, and the later stage of the fracture permeability recovery as the gas desorption expands and the matrix permeability increase due to the slip effects.