An Analytical Flow Model for Heterogeneous Multi-Fractured Systems in Shale Gas Reservoirs
Honghua Tao1,2,3, Liehui Zhang1, Qiguo Liu1, Qi Deng1, Man Luo4 and Yulong Zhao1,
1State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China. 2John and Willie Leone Family Department of Energy and Mineral Engineering, The Pennsylvania State University, University Park, PA 16802, USA. 3Key Laboratory of Shale Gas Exploration, Ministry of Land and Resources Engineering, Chongqing 400042, China. 4Petro China West Pipeline Company, Urumqi 830013, China
Received: 14 October 2018 / Accepted: 4 December 2018 / Published: 6 December 2018
The development of shale gas in North America has achieved large-scale commercial success [1,2,3], which has set off a shale gas revolution worldwide. As a key technology in shale gas exploration and development, well testing plays an irreplaceable role. The characteristics of shale gas reservoirs can be obtained through the transient pressure analysis of multiple fractured horizontal wells (MFHWs) in shale gas reservoirs.
In order to describe the random and complex fractures, some works [4,5] have investigated discrete fracture networks through numerical simulation approaches.
Tang et al.  established a three-dimensional numerical model based on the construction of spatial discretization by the finite volume method. Wang  proposed a unified model for shale gas reservoirs based on discrete fracture networks to investigate shale gas production by rate transient analysis. However, this requires numerical simulation, and the process is time-consuming and occupies a large amount of computing resources.
Fortunately, the analytical approach is a convenient and rapid method for the evaluation of dynamic characteristics of the shale gas reservoir, which takes less time and needs less reservoir data compared with numerical simulation approaches. Thus, the analytical approach has been used in more applications in recent years.
Two types of analytical model are used to analyze transient pressure behaviors. One type is the detailed analytical model, which is based on the source function and superposition principle [6,7,8]. This characterizes the stimulated reservoir volume (SRV) region in a shale gas reservoir as a circular or rectangular zone and extends the one-region model to a dual-region composite model. Similarly, the shortcomings of the detailed analytical model also cause a large increase in the amount of calculation required. In order to describe the SRV region more concisely and conveniently, the other type, which is linear models, such as the tri-linear flow model  and the five-region flow model , was developed. The five-region flow model was established based on the tri-linear flow model and takes into account not only the stimulated region, but also the nearby unstimulated region. These two models represent a rapid way to capture key characteristics in shale gas reservoirs.
Based on these two analytical models (the detailed analytical model and linear model), other improved models were developed, e.g., models considering the effects of fractures in the SRV region , non-equal spacing fractures , fracture networks in the shale matrix [13,14], the non-Darcy high-speed flow inside the hydraulic fracture , the shale matrix diffusion and dual porosity model , a transient flow approach , and non-Darcy flow with a threshold pressure gradient in tight gas reservoirs . Recently, Zeng et al. , Zeng , and Zeng et al.  proposed a seven-region flow model, which takes into account the spatial heterogeneity and typical seepage features, such as ad-desorption and diffusion in shale gas reservoirs. Unfortunately, all of the models described above only consider the linear flow in all regions, and thereby neglect the fractal features and sub-diffusive flow in the SRV region.
In order to capture the features of fractal geometry and sub-diffusive flow in highly heterogeneous porous media, an analytical flow model that considers anomalous diffusion and other significant features to describe the flow characteristics in the SRV region has been proposed [22,23,24,25,26].
Chen and Raghavan  utilized fractional derivatives to characterize the process of anomalous diffusion in the complex fractures and took into account the continuous-time random walk in hydraulically fractured reservoirs with single porosity. Subsequently, Ren and Guo  presented a dual porosity and anomalous diffusion model for shale gas reservoirs. Unfortunately, they did not consider the heterogeneity of multi-fractured systems by applying a three-region or five-region model. Later, Albinali and Ozkan  proposed a tri-linear anomalous diffusion and dual-porosity model that uses fractional calculus to account for non-uniform velocity in porous media. However, the fractal geometry features of the induced fractures in the SRV region are not considered in the model. Wang et al.  considered the fractal characteristics in the complex system by coupling fractal relations to account for the heterogeneity in the SRV region. Fan and Ettehadtavakkol  applied micro-seismic data to verify the fractal flow model and proposed a semi-analytical model for rate transient analysis in shale gas reservoirs.
All the models described above do not fully consider the various diffusion of shale gas in the shale matrix, the dual porosity in the SRV region, or the stress sensitivity of permeability and fractal-anomalous diffusion in complex fractures. Table 1 demonstrates the differences by comparing previous analytical flow models with the present model. Previous models only considered homogeneous properties and simple transport mechanisms in shale gas reservoirs.
Table 1. Feature comparisons of analytical models for multiple fractured horizontal wells (MFHWs). SRV: stimulated reservoir volume.
Based on the above, this work proposes a new analytical model based on fractal-anomalous diffusion. Firstly, the present model is coupled with anomalous diffusion and other significant features, such as ad-desorption, slip flow, surface flow, pressure-dependent permeability, and fractal geology. Using the Laplace transformation method and Duhamel’s theorem , the analytical solution of the present model is obtained. Then, the flow regimes are identified, and the effects of relevant parameters are analyzed.
Therefore, the present model can effectively describe the complex fracture networks in the SRV region and more accurately account for the various transport mechanisms of MFHWs in shale gas reservoirs. Due to the lack of well-testing data in shale gas reservoirs, the present model has only been applied to one case, but more cases will be studied in the future.
Figure 1 is a schematic of the typical five-region flow model and the improved five-region flow model (new model) in a shale gas reservoir. Higher fractal permeability, dual-porosity, and anomalous diffusion in the SRV region are taken into account around each fracture. The other three regions occupy the remaining space between adjacent fractures. One-quarter of each hydraulic fracture is taken into account due to the assumption of symmetry in the reservoir.
Figure 1. Schematic of physical models for hydraulically fractured horizontal wells. (a) The typical five-region flow model proposed by Stalgorova and Mattar . (b) The improved five-region flow model (new model). Fracture half-length: xf; width of the hydraulic fracture: wD; distance from the hydraulic fracture to stimulated reservoir volume (SRV): y1; no flow bound: x2,y2.
As shown in Figure 1, the reservoir between two adjacent fractures is subdivided into five regions in the improved five-region flow model. There is vertical linear flow from region 4 to region 2 and from region 3 to region 1 (SRV). Similarly, horizontal linear flow exists from region 2 to region 1 and from region 1 to each hydraulic fracture. Compared with the typical five-region flow model, ad-desorption and various diffusion in the shale matrix, dual-porosity (shown as spherical matrix in Figure 1), the fractal geometry (shown as a power-law type in Figure 1) and anomalous diffusion (sub-diffusion) in the SRV region, and stress-sensitive permeability in each region are considered in this work. The main assumptions of this new model are as follows:
- A hydraulically fractured horizontal well is at the center of a closed shale gas reservoir;
- Each hydraulic fracture is perpendicular to the horizontal well, spaced uniformly along the horizontal wellbore, and has the same length;
- Fluid flow in each region is a one-dimensional single-phase flow;
- Desorption in shale matrix yields to the Langmuir isotherm adsorption law;
- The continuity of flux and pressure at interfaces is used to couple the adjacent regions.
Mechanisms and Properties
Adsorption/Desorption and Apparent Permeability
Shale gas adsorption in the shale matrix typically yields to the Langmuir isotherm adsorption law, and pseudo-pressure can be written as follows [28,29]:
where VE is defined as the adsorption equilibrium concentration (sm3/m3), the Langmuir adsorption concentration is represented by VL (sm3/m3), the Langmuir pressure is represented by PL (MPa), and P means the pressure in the reservoir (MPa).
where σm is the adsorption factor.
where m(p) is the pseudo-pressure (MPa2/(mPa·s)), the gas viscosity is represented by μ (mPa·s), and the real gas deviation factor is represented by z.
The main transport mechanisms in the shale matrix are surface diffusion, Knudsen diffusion, and slip flow. Based on the results of previous research, the expression of total equivalent permeability (apparent permeability) is as follows :
where kma is defined as an apparent permeability which is related to surface diffusion, Knudsen diffusion, and slip flow (m2); ke is the equivalent slip rate of slip flow (m2); the Knudsen diffusion equivalent permeability is represented by kd (m2); the surface diffusion equivalent permeability is represented by ks (m2); and βt is the matrix comprehensive diffusion factor that considers the slip flow, Knudsen, and surface diffusion.
Fractal Permeability and Porosity in Induced Fractures
The distribution of induced fractures is extremely complex and irregular, and therefore, it is not accurate enough to describe the porosity of induced fractures in Euclidean geometry. Fractal geometry has been verified as an effective method to describe the complex pore structure of fibrous porous media [31,32,33,34]. Based on fractal geometry, fractal permeability and fractal porosity in induced fractures comply with a power-law type as follows [35,36,37,38]:
where Kfr is the permeability at the reference length, Lref is the reference length; the mass fractal dimension of the inducec fractures is represented by df, the Euclidean dimension is represented by de, the radial coordinate value at any point is represented by r, and the tortuosity index is represented by θ.
where ∅fr is the porosity at the reference length.
Anomalous Diffusion in Induced Fractures
In induced fractures, the disorder, non-local, and memory features should be considered in the SRV region. This complex transport process is anomalous diffusion, which is described by fractional calculus. The modified Darcy flow velocity is given by the following form :
The fractional derivative is defined as follows :
where the Gamma function is represented by Γ(x). The Laplace transform of the fractional derivative is
when α=1, Equation (7) is reduced to the classical Darcy’s law as follows :
The permeability in hydraulically fracturing shale gas reservoirs is sensitive to pore pressure, according to previous experiments [3,40]. Given the relationship with pore pressure, fractal permeability is introduced by permeability modulus as follows:
where ki is the permeability under the initial pseudo-pressure (mi), the corresponding pseudo-pressure in the reservoir is represented by m, and γ is the stress sensitivity factor.
Governing Flow Equations and Solutions
In order to obtain the final solution, the governing diffusivity equations for each region are written with the relevant initial and boundary conditions. Definitions of all dimensionless terms are given in Appendix A.
Unstimulated Regions (Region 4 + Region 3 + Region 2)
Starting with the fourth region, the diffusivity equation that considers the ad-desorption and various diffusion can be written in a dimensionless form:
where γ∗D is the dimensionless stress-sensitive factor.
The perturbation inversion proposed by Pedrosa Jr.  is applied to pseudo-pressure, as presented in Equation (13).
Additionally, a zero-order approximation is performed to linearize the diffusivity equation. Then, the diffusion Equation (12) can be approximately written in a Laplace form, as follows:
The outer boundary condition (no-flow) is
The inner boundary condition (pressure continuity) is
Therefore, the general form of the solution in the fourth region can be given as follows:
and η4D is the dimensionless conductivity in region 4.
Region 3, which has low permeability, can only flow vertically to region 1. Similarly, a general form of the solution for the third region can be given as follows:
Also, the governing equation of region 2 becomes
The outer boundary condition (no-flow) is
The inner boundary condition (pressure continuity) is
Therefore, the solution for region 2 becomes
Region 1 (SRV)
Region 1 represents the SRV region in which the transient inter-porosity flow from the matrix to fracture subsystem is applied. Moreover, the anomalous diffusion, fractal permeability, and porosity in induced fractures are also considered.
- Matrix subsystem:
Similarly, the pressure solution in the matrix subsystem of region 1 can be obtained:
- Induced fractures subsystem:
The diffusivity equation of the complex fractures networks can be derived in the following dimensionless form. More detailed derivations are given in Appendix B.
The outer boundary condition (flow continuity) is
The inner boundary condition (pressure continuity) is
Therefore, the general form of the pressure solution in the SRV is
Hydraulic Fracture Region
Considering that the stress sensitivity of permeability and flow exchange is directly related to the quality dimension, the diffusivity equation in hydraulic fractures becomes
The perturbation inversion  and zero order approximation in the Laplace form are applied, and the diffusivity equation then becomes
Equation (35) can be written as follows:
Boundary condition 1 is
Boundary condition 2 is
The pressure solution for the hydraulic fracture region is
Thus, the pressure solution at the wellbore can be given as follows:
However, by applying the superposition principle and Duhamel’s principle , the final solution for wellbore pressure considering convergence and storage is written as follows:
Then, the perturbation inversion  and Stehfest numerical inversion  are applied. Finally, the pressure solution at the downhole can be written with the real-time data as