1 Introduction

With the growing shortage of domestic and foreign conventional reservoir, the production of gas from shale formation has played an increasingly important role over recent years, and is gradually becoming the key component in the energy supply all over the world. Since shale gas reservoir is tight with nanopores, it is improper to characterize the gas flow with conventional Darcy flow or Fick diffusion. Due to the ultra-low porosity and permeability, the multi-flow mechanisms, e.g., Darcy flow, slippage, diffusion, and desorption, coexist in shale gas reservoirs. These complexities bring out significant challenges in mathematical modeling and simulation approaches for the prediction of shale gas productivity.

Scholars have presented many models to describe the multi-scale structures and nonlinear percolation mechanism for multi-stage fracturing horizontal wells (MFHWs). Ozkan et al.^[1] used the tri-linear flow model to simulate the pressure transient and production behavior of MFHWs in unconventional reservoirs, e.g., the fracture network around the hydraulic fractures. Brown et al.^[2] presented an analytical trilinear-flow solution to simulate the MFHWs. In their model, the internal zone between the hydraulic fractures was filled with natural fractures, and the linear flow was assumed to be the only flow regime. Brohi et al.^[3] developed the composite model with the assumption of a linear dual porosity flow for the inner reservoir and a linear single porosity for the outer reservoir combined with the continuity of pressure and flux at their interface. Fan et al.^[4] suggested a composite model by dividing the shale reservoir into three main areas, i.e., hydraulic fractures, stimulated reservoir volume (SRV), and unstimulated reservoir, and found that, according to the arrival time of pressure wave to pseudo-steady flow of transition area, the pressure wave arrival to the SRV boundary could be decided, and the size of the SRV could be determined by the time. Zhang et al.^[5] proposed a dual-porosity model with the consideration of the Langmuir adsorption isotherm and Fick's law to analyze the production decline performance, and pointed out that the transition flow regime and production performance of an MFHW in shale are mainly dominated by the Langmuir volume, radius, and permeability of the SRV.

In fact, the gas flow in extremely low permeability shale formations is a complex process with the coexistence of multi-mechanisms, such as viscous flow, Knudsen diffusion, slip flow, and desorption. Zhao et al.^[6] used the Knudsen diffusion to describe the gas flow in nano-pores with the consideration of the slippage effect, and concluded that the Knudsen diffusion coefficient and slippage factor mainly affected the interporosity flow period. Beskok and Karniadakis^[7] developed a unified multi-scale permeability model for tight porous media, including all these flow regimes (0 ≤ K_n ≤ ∞). Michel et al.^[8] simplified the Beskok-Karniadakis model by assuming the slip coefficient as a theoretical value -1. However, the new model is applicable only for the case when K_n ≤ 0.1. Javadpour^[9] proposed the concept of apparent permeability with the consideration of the Knudsen diffusion, slippage flow, and advection flow, and found that the apparent permeability was higher than the Darcy permeability for smaller pores and lower pressures.

The numerical method and analytical solution for the shale gas nonlinear performance have also been studied. Kabir and Hasan^[10] used the perturbation theory to treat the nonlinearity of the product of viscosity and compressibility. Mireles^[11] considered the nonlinear compressibility by building convolution, where the product of viscosity and compressibility was simplified as a parameter changing with the average pressure and the nonlinear compressibility was transformed from real space to Laplace space. By dividing the fracture wings into several segments with different flow rates and using the boundary element method to solve the partial different equations in the Laplace domain, Zerzar et al.^[12] obtained a comprehensive solution for MFHWs. Zhang et al.^[13] proposed an effective mean-free path to address the Knudsen layer effect so that the lattice Boltzmann methods could be extended beyond the slip-flow regime. They showed that the proposed approach could significantly improve the near-wall accuracy of the lattice Boltzmann method, and provide a computationally economic solution technique over a wide range of Knudsen numbers. Barreto et al.^[14] treated the effect of the viscosity-compressibility product change as a nonlinear source term, and solved it by one iteration with the Green function.

This work aims to establish a unified mathematical model under the nonlinear multi-scale effects, such as slippage, diffusion, and desorption. The continuous inhomogeneous models of equivalent porosity and permeability are proposed for the whole shale gas reservoir. Then, the related semi-analytical approach is developed. By comparing the simulation results with the field data from MFHWs, the applicability of the present mathematical model and the implemented approach are verified. The nonlinear multi-scale effects, such as slippage, diffusion, and pressure dependent desorption, on the shale gas production are investigated.

2 Nonlinear mathematical modeling with inhomogeneity 2.1 Mathematical formulation for unsteady percolation

The shale gas is stored by either free or adsorption conditions. The free gas is stored in both the micro-pores in the matrix and the natural fractures, and the adsorption gas is adsorbed on the surface of the matrix particles. With the decrease in the producing pressure, the adsorption gas will be desorbed into the micro-pores, and then flow into the natural fractures^[15].

Considering the loss of gas mass by adsorption, we can write the continuity equation as follows^[16]:

(1)

where t is the time, and the shale porosity ϕ is considered as the inhomogeneous distribution in the whole shale gas reservoir.

According to the gas state equation, the density of real gas ρ_g under the isothermal condition is

(2)

where the pressure function p=p(r, z, t), and the constants ρ_gsc, Z_sc, T_sc, and p_sc are the known density, compressibility factor, temperature, and pressure of gas at the standard state, respectively. T is the formation temperature, and Z denotes the compressibility factor of real gas.

The volumetric flux vector of the flowing gas v and the pressure p satisfies the following non-Darcy seepage law:

(3)

where μ is the viscosity of the shale gas. The apparent permeability K is considered to be inhomogeneous and pressure dependent in the whole shale gas production area.

In terms of the Langmuir isotherm adsorption model, the gas mass by the adsorption per unit volume of the shale gas q_ad can be expressed as

(4)

where V_L is the Langmuir gas volume, and p_L is the Langmuir gas pressure.

Substituting Eqs. (2)–(4) into Eq. (1) yields the following equation:

(5)

where

For a fractured horizontal well in the shale gas reservoir, the boundary conditions of the pressure are

(6a)

(6b)

where the wellbore radius r_w is tiny, and it can be approximated as r_w ≈ 0. p_w and p_e are the given wellbore pressure and the supply boundary pressure, respectively. ε is the half-width of the gas outlet. L is the half-distance between two adjacent hydraulic fractures. v_r and v_z are the velocities perpendicular and parallel to the wellbore, respectively.

The initial condition of the pressure is

(7)

2.2 Inhomogeneity of equivalent porosity

After the hydraulic fracturing of the shale gas reservoir, a large number of artificial macro-fractures are created near the wellbore, and many micro-fractures correspondingly get extended around these macro-fractures (hydraulic fractures), which develop a huge fractured network. The shale gas flows from the matrix to the micro-fractures, then through the hydraulic macro-fracture, and finally into the horizontal wellbore. When the hydraulic fractures are periodically arranged, the symmetry plane between two adjacent hydraulic fractures is a normal non-flow boundary (see Fig. 1).

The most used partition model usually supposes that the porosity and intrinsic permeability have different values in the hydraulic fracture, the micro-fracture, and the matrix zones^[17], i.e., the porosity and intrinsic permeability are homogeneous and assigned as the average value in each partitioned zone. In fact, the distributions of the porosity and intrinsic permeability in the whole shale gas reservoir should be continuous and inhomogeneous. In view of the energy dissipation caused by fracture propagation, the values of porosity and intrinsic permeability would reduce with the increases in the distances parallel and perpendicular to the hydraulic fractures. Based on the three partition models (see Refs. [1], [15], [17], and [18]), the present model intends to characterize the continuous inhomogeneity of porosity in the whole gas production area, in which a uniform mathematical expression of the porosity distribution suited for the whole gas production area is proposed. Its advantage is that the difficulties from satisfying the continuous conditions of gas flux and pressure can be eliminated at the interfaces of the partitioned regions.

To characterize the continuous inhomogeneity of porosity in the whole shale gas reservoir, the porosity ϕ can be approximately assumed as follows:

(8)

The porosity model in the z=0 plane should satisfy the following constrained conditions:

(9)

This yields the constants A and B in Eq. (8) as follows:

(10)

where ϕ_f, ϕ_mf, and ϕ_m are the porosities of the hydraulic fracture zone, the micro-fracture zone, and the matrix zone, respectively. r_c is the radius of the fracture half-length, and n is the dimensionless constant to describe the porosity distribution aggregation along the r-direction.

Figure 2 reveals the continuous inhomogeneous distribution of porosity in the rz-plane. The parameters used for the computation are

It can be seen that the distribution of the porosity along the radial and axial directions is continuous, and decreases gradually with the increasing distance to the wellbore.

2.3 Equivalent permeability with inhomogeneity

For the fractures and matrix, the intrinsic permeability can be expressed as the cubic law of the corresponding porosity as follows^[19]:

(11)

where K₀ is the intrinsic permeability, and

K_f, K_mf, and K_m are the permeabilities of the hydraulic fracture zone, the micro-fracture zone, and the matrix zone, respectively.

The constant C is required to approach to C_f and C_m when the porosity ϕ is close to ϕ_f and ϕ_m, respectively. To characterize the continuous inhomogeneity of the intrinsic permeability in the whole reservoir, the weighted average values of C_f and C_m can be introduced as follows:

(12)

where the weight coefficient γ(0≤ γ ≤1) can be expressed as

(13)

Figure 3 illustrates visually the distribution curve of the intrinsic permeability obtained from Eqs. (11)-(13) in the rz-plane. It can be observed that the distribution of intrinsic permeability decreases along the radial and axial directions, respectively, and the symmetry plane between two adjacent hydraulic fractures is a non-flow boundary because of the assumption of identical hydraulic fractures.

On the basis of the Beskok-Karniadakis model, the apparent permeability of the gas-bearing shale is given as follows^[18]:

(14)

where K is the apparent permeability, f(Kn) is the permeability correction factor, the dimensionless rarefaction coefficient α can be obtained by , and b denotes the slip coefficient which is usually taken as the value -1. The Knudsen number Kn can be expressed as^[18]

(15)

where D_k is the Knudsen diffusivity, and r is the character pore size of the porous media. It should be pointed out that the apparent permeability model expressed by Eq. (14) can be reduced into the classic Darcy seepage model if the intrinsic permeability is homogeneous and the permeability correction factor f(Kn)=1.

3 Semi-analytical solution

By the Boltzmann transformation, the variable is introduced as follows:

(16)

Then, the pressure p can be thought as a function of the new variable u. We divide the domain Ω ={(r, z)| r_w ≤ r≤ r_e, 0≤ z≤ L} into N× M subdomains, where the outer boundary radius r_e is selected to be sufficiently large. In the subdomains

the dimensionless intrinsic permeability

and the porosity ϕ (r, z)≈ ϕ_ij =ϕ (r_i, z_j). For a given time t, the pseudopressure function and the masssource function can be, respectively, approximated in Ω_ij as follows:

(17)

(18)

Thus, the partial differential equation (5) in Ω_ij can be converted into the following ordinary differential equation:

(19)

The boundary conditions in Eq. (6) and the initial condition in Eq. (7) for the pressure p can be rewritten as

(20)

Integrating twice to the two sides of Eq. (19) with respect to the variable u results in

(21)

where the integral constant D₀ =0 because Eq. (17) implies that m(u_w)=0. D is the other integral constant to be determined.

By exchanging the integral order and then integrating by parts, we can simplify the expression (21) into the single integration form as follows:

(22)

By eliminating the pseudopressure function m from Eqs. (17) and (22), the pressure field p(u) in the subdomain Ω_ij satisfies the following integral equation:

(23)

Implementing the total differentiation to the two sides of Eq. (23) yields the relationship between the pressure increment Δ p and the variable increment Δ u as follows:

(24)

For a given time t and space joint (r_i, z_j) in the subdomains Ω_ij (i=0, 1, ..., N; j=0, 1, ..., M), the corresponding variable u and pressure p are denoted as u_{i, j} =(r_i²+z_j²)/(4t) and p_{i, j} =p(u_{i, j}), respectively. For a sufficient small step Δ u_{i, j}, the pressure has an increment Δ p_{i, j}. From Eq. (24), by applying the rectangular rule to the integration, the explicit iterative scheme of pressure can be expressed as follows:

(25)

where ξ_{i, j} =ξ (p_{i, j}), and F_{i, j} =F(p_{i, j}).

For a given time t, the procedure of the numerical calculation for the pressure field p(r, z, t) is given as follows.

First, set j = 0 in the z = 0 plane and take the ith-node along the r-direction, i.e.,

(26)

where the variable step size Δ u_{i, 0} =u_{i+1, 0} -u_{i, 0}, and enlarges with the increase in the sequence number i. Denote the pressure value at the node of u_{i, 0} as p_{i, 0}. Select a suitable value D, and run the iterative scheme (25) with the known initial value of pressure p_{1, 0} =p_w until the error criterion |(p_{N, 0}-p_{e, 0})/p_{e, 0}| < 10^-R, in which R is the error control index, is satisfied. Increase the number N gradually, and repeat the above step until the constant D is convergent. Then, the pressure value p_{i, 0} at each node u_{i, 0} can be determined. In the z=0 plane, with the relationship between u and r in Eq. (16), we can obtain the pressure value p_{i, 0} for the given time at each node (r_i, 0) (i=0, 1, ..., N).

Next, in the plane of the given radial coordinate r_i, take the jth-node along the z-direction as follows:

(27)

The step Δ u_{i, j} =u_{i, j+1} -u_{i, j}, and enlarges with the increasing sequence number j. Beginning with the known pressure values p_{i, 0} (i=0, 1, ... N) obtained from the first step, we select a suitable value D, and run the iterative scheme (25) until the error criterion |(p_{N, M}-p_{e, e})/p_{e, e}| < 10^-R, in which R is the error control index, is satisfied. Increase the number N gradually, and repeat the above step until the constant D is convergent. Then, we can generate the value of pressure at each node p_{i, j}.

Finally, renew the time t, and repeat the above two steps. Then, we can obtain the pressure value p_{i, j} =p(r_i, z_j, t).

4 Results and discussion

The proposed mathematical model and semi-analytical method are effectively verified by comparing the simulation results with the field data from the fractured horizontal well in Refs. [20] and [21]. The nonlinear multi-scale effects, such as slippage, diffusion, and desorption, on the gas production are discussed in detail. In practical applications, the product of μ and Z is approximately equal to its average value under the average pressure of the shale gas reservoir. The fundamental parameters used in the simulation are given as follows.

(ⅰ) The gas state parameters are

(ⅱ) The gas desorption parameters are

(ⅲ) The seepage parameters are^[18]

Figure 4 shows the convergence validation of the parameter D. It can be seen that when N>20, the parameter D is convergent. For this study, we take N=100.

Figure 5 exhibits the comparison between the simulation results and the field data from the Marcellus shale with six hydraulic fractures. The practical parameters are

In the present numerical simulation, the other parameters are taken as follows:

As shown in Fig. 5(a), the simulation results with the consideration of the nonlinear multi-scale effects of slippage, diffusion, and desorption agree well with the 210-day field test data^[15].

This verifies the applicability of the proposed mathematical model and the semi-analytical method. Figure 5(b) shows the curves of the gas production rate and the cumulative gas production for the 25-year period. It can be seen that the gas production rate is relatively high but drops quickly at the early period of production. At the later period, the production rate is relatively low with a small decline rate.

Figure 6 shows the comparison between the simulation results by the present model and the results from Ref. [21] by the numerical reservoir simulator of CMG-IMEX for the gas production in the Marcellus shale reservoir. The length of the well is 426.7 m with 28 hydraulic fractures. The practical parameters are

In the present numerical simulation, the other parameters are taken as follows:

The results indicate that the present model agrees well with the numerical simulation results in Ref. [21]. It can also be seen that the gas desorption plays a positive role on the gas production rate, but has a smaller effect during the early period of production.

Except for the aforementioned fundamental parameters, the parameters used in the following simulation examples in Figs. 7-9 are listed as follows:

To examine the effects of the slippage on the production performance, Fig. 7 shows the production curves for different slippage factors. For the 25-year production process, the relative difference of the cumulative gas yield between the cases b=1 and b=-2 approaches around 11%, and the cumulative gas production between the case with the consideration of the nonlinear multi-scale effect when b=-2 and the Darcy seepage model is 17%, which is in the evidence that the slippage factor is sensitive to the gas production.

Figure 8 shows that, the gas diffusion coefficient has a great effect on the gas production. This may be because that a smaller diffusion coefficient makes it harder for the shale gas adsorbed in the shale matrix to diffuse into the micro-fracture zone, and the adsorbed gas needs a very long time to compensate for the pressure drop. The cumulative gas yield between the cases D_k= 2×10⁻⁷ m²/s and and D_k= 8×10^-7 m²/s approaches around 13% in 25 years. The relative difference of the cumulative gas production between the proposed model with the consideration of the nonlinear multi-scale effect with D_k=8×10⁻⁷ m²/s and the Darcy model is 19%, which indicates that the classic Darcy seepage model underestimates the gas production to the present model.

The comparison between the curves of the cumulative gas production with and without the consideration of the desorption effect is shown in Fig. 9. The gas desorption has a positive effect on the gas production, a smaller effect during the early period of production, while a higher contribution during the late period of production because of the substantial pressure depletion and larger gas drainage area, which finally contributes to 18% of the total gas production in around 25 years.

5 Conclusions

By the continuization of the traditional homogeneous partition models of porosity and permeability, we propose the continuous inhomogeneous distribution models of equivalent porosity and permeability for the whole shale gas reservoir including the hydraulic fractures, the micro-fractures, and the matrix regions. On the basis of the Beskok-Karniadak model, a unified mathematical formulation is established with the consideration of the nonlinear multi-scale effects, such as slippage, diffusion, and desorption. Then, to simulate the nonlinear unsteady percolation of the shale gas reservoir, the related semi-analytical numerical method is developed. By means of the Boltzmann transformation, the initial-boundary value problem of the nonlinear partial differential equation for shale gas transport is converted and simplified into the equation in integration form. A concise iteration scheme to solve the pressure field is constructed by the numerical discretization of the integral equation.

The numerical simulation results match well with the field data from the fractured horizontal wells, which illustrates the availabilities of the present mathematical model and the implemented approach. The nonlinear multi-scale effects, such as slippage, diffusion, and desorption, on the gas production are investigated. The slippage factor is sensitive to the gas production, and the gas diffusion coefficient has a great effect on gas production, which indicates that the classic Darcy seepage model underestimates the gas production of the present model. The gas desorption has a smaller effect during the early period of production while a higher contribution during the late period of production.