An Analytical Flow Model for Heterogeneous Multi-Fractured Systems in Shale Gas Reservoirs

: The use of multiple hydraulically fractured horizontal wells has been proven to be an efﬁcient and effective way to enable shale gas production. Meanwhile, analytical models represent a rapid evaluation method that has been developed to investigate the pressure-transient behaviors in shale gas reservoirs. Furthermore, fractal-anomalous diffusion, which describes a sub-diffusion process by a non-linear relationship with time and cannot be represented by Darcy’s law, has been noticed in heterogeneous porous media. In order to describe the pressure-transient behaviors in shale gas reservoirs more accurately, an improved analytical model based on the fractal-anomalous diffusion is established. Various diffusions in the shale matrix, pressure-dependent permeability, fractal geometry features, and anomalous diffusion in the stimulated reservoir volume region are considered. Type curves of pressure and pressure derivatives are plotted, and the effects of anomalous diffusion and mass fractal dimension are investigated in a sensitivity analysis. The impact of anomalous diffusion is recognized as two opposite aspects in the early linear ﬂow regime and after that period, when it changes from 1 to 0.75. The smaller mass fractal dimension, which changes from 2 to 1.8, results in more pressure and a drop in the pressure derivative.


Introduction
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. [4] established a three-dimensional numerical model based on the construction of spatial discretization by the finite volume method. Wang [5] 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.  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.

Physical Model
Energies 2019, 12 3 of 18 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. ; width of the hydraulic fracture: ; distance from the hydraulic fracture to stimulated reservoir volume (SRV): 1 ; no flow bound: 2 , 2 .

Physical Model
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,  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: (1) A hydraulically fractured horizontal well is at the center of a closed shale gas reservoir; (2) Each hydraulic fracture is perpendicular to the horizontal well, spaced uniformly along the horizontal wellbore, and has the same length; (3) Fluid flow in each region is a one-dimensional single-phase flow; (4) Desorption in shale matrix yields to the Langmuir isotherm adsorption law; (5) The continuity of flux and pressure at interfaces is used to couple the adjacent regions.

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 V E is defined as the adsorption equilibrium concentration (sm 3 /m 3 ), the Langmuir adsorption concentration is represented by V L (sm 3 /m 3 ), the Langmuir pressure is represented by P L (MPa), and P means the pressure in the reservoir (MPa).
where σ m is the adsorption factor.
where m(p) is the pseudo-pressure (MPa 2 /(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 [30]: where k ma is defined as an apparent permeability which is related to surface diffusion, Knudsen diffusion, and slip flow (m 2 ); k e is the equivalent slip rate of slip flow (m 2 ); the Knudsen diffusion equivalent permeability is represented by k d (m 2 ); the surface diffusion equivalent permeability is represented by k s (m 2 ); 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 K fr is the permeability at the reference length, L re f is the reference length; the mass fractal dimension of the inducec fractures is represented by d f , the Euclidean dimension is represented by d e , 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 [22]: The fractional derivative ∂t ∝ is defined as follows [39]: where the Gamma function is represented by Γ(x). The Laplace transform of the fractional derivative when α = 1, Equation (7) is reduced to the classical Darcy's law as follows [23]:

Pressure-Dependent Permeability
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 k i is the permeability under the initial pseudo-pressure (m i ), 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. [41] 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: where 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: where 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 where

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: where -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. where 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 where

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 where The perturbation inversion [41] 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 [27], the final solution for wellbore pressure considering convergence and storage is written as follows: Then, the perturbation inversion [41] and Stehfest numerical inversion [42] are applied. Finally, the pressure solution at the downhole can be written with the real-time data as

Flow Regimes
In order to obtain the main flow regimes of the improved five-region flow model, the type curves of the pressure-transient response were plotted by employing pseudo-steady inter-porosity flow in the SRV region.
The related parameters are listed in Table 1. Figure 2 shows

Sensitivity Analysis
In the corresponding sensitivity analysis, firstly, one relevant parameter was changed while keeping the other parameters at their original values. Then, all the relevant parameters were changed at the same time. Model parameters were given values in the simulation by referring to relevant literature [6,12,16,17,23,25,26], and they are listed in Table 2.

Parameter Name
Parameter Value

Sensitivity Analysis
In the corresponding sensitivity analysis, firstly, one relevant parameter was changed while keeping the other parameters at their original values. Then, all the relevant parameters were changed at the same time. Model parameters were given values in the simulation by referring to relevant literature [6,12,16,17,23,25,26], and they are listed in Table 2.

Parameter Name Parameter Value
Dimensionless half fracture length, x fD 1 Dimensionless fracture conductivity, F CD 2 Inter-porosity flow coefficient, λ λ = 0.2 Storage capacity coefficient, ω ω = 0.2 Dimensionless distance in x direction, x nfD /x eD x nfD = 1, x eD = 50 Dimensionless distance in y direction, y nfD /y eD y nfD = 1, y eD = 50 Ratio of permeability, k i /k j k 3a /k nf = 0.0005, k 2a /k nf = 0.1, k 4a /k 2a = 0.02 Absorption factor, σ m 5 Diffusion factor (apparent permeability coefficient), β t 1.  Figure 3 shows that the fracture conductivity mainly affects the early flow stages. The greater the fracture conductivity is, the smaller the gas flow resistance is, and the smaller pressure consumption is with the same production. It is not difficult to see that the fracture conductivity mainly influences the pressure and pressure derivative curves in the bilinear flow and first linear flow stages. With an increase in the fracture conductivity, the duration of the bilinear flow stage decreases and the duration of the first linear flow stage increases. As seen in Figure 3, when F CD = 25, only the first linear flow regime can be observed.
Energies 2019, 12 10 of 18 Figure 3 shows that the fracture conductivity mainly affects the early flow stages. The greater the fracture conductivity is, the smaller the gas flow resistance is, and the smaller pressure consumption is with the same production. It is not difficult to see that the fracture conductivity mainly influences the pressure and pressure derivative curves in the bilinear flow and first linear flow stages. With an increase in the fracture conductivity, the duration of the bilinear flow stage decreases and the duration of the first linear flow stage increases. As seen in Figure 3, when FCD = 25, only the first linear flow regime can be observed.  Figure 4 demonstrates the type curves of the pressure and pressure derivative for MFHWs in a shale gas reservoir with various anomalous diffusion exponent (α) and tortuosity index (θ) values. As can be seen, one intersection point exists between the anomalous diffusion and classical diffusion pressure derivative curves. At the early bilinear and linear flow stages, the pressure and pressure derivative for α < 1 or θ > 0 (anomalous diffusion) are smaller than those for α = 1 or θ = 0 (classical diffusion). When the value of α increases (θ decreases), the pressure and its derivative will also increase. The reason for this is that anomalous diffusion delays the performance of pressure derivative behaviors. However, after the inter-porosity flow stage, with different α values, the difference will be more obvious, and the trend is the opposite. In other words, a decrease in α (θ increasing) causes the pressure and its derivative to increase over time. This accounts for the characteristic of sub-diffusion (slower flow) when α < 1 or θ > 0 (anomalous diffusion).  Figure 4 demonstrates the type curves of the pressure and pressure derivative for MFHWs in a shale gas reservoir with various anomalous diffusion exponent (α) and tortuosity index (θ) values. As can be seen, one intersection point exists between the anomalous diffusion and classical diffusion pressure derivative curves. At the early bilinear and linear flow stages, the pressure and pressure derivative for α < 1 or θ > 0 (anomalous diffusion) are smaller than those for α = 1 or θ = 0 (classical diffusion). When the value of α increases (θ decreases), the pressure and its derivative will also increase. The reason for this is that anomalous diffusion delays the performance of pressure derivative behaviors. However, after the inter-porosity flow stage, with different α values, the difference will be more obvious, and the trend is the opposite. In other words, a decrease in α (θ increasing) causes the pressure and its derivative to increase over time. This accounts for the characteristic of sub-diffusion (slower flow) when α < 1 or θ > 0 (anomalous diffusion).
Energies 2019, 12 11 of 18 Figure 4. Effect of the anomalous diffusion exponent on type curves. Figure 5 shows that the mass fractal dimension of induced fractures (Hausdorff index) has a significant effect on the pressure behavior at almost all the stages, except for the wellbore storage stage. Overall, the smaller the mass fractal dimension is, the larger the gas flow resistance is and the greater the pressure consumption is with the same production. As can be seen, the locations of the type curves are higher with a smaller df. The reason for this is that a smaller df value represents more resistance in the complex induced fractures.  Figure 6, the adsorption factor mainly influences the position of the type curves at the inter-porosity flow stage. A larger adsorption factor represents a stronger adsorption and production capacity and therefore makes the "concave" appear wider and deeper on the type curves. Figure 7 shows the effect of the  Figure 5 shows that the mass fractal dimension of induced fractures (Hausdorff index) has a significant effect on the pressure behavior at almost all the stages, except for the wellbore storage stage. Overall, the smaller the mass fractal dimension is, the larger the gas flow resistance is and the greater the pressure consumption is with the same production. As can be seen, the locations of the type curves are higher with a smaller d f . The reason for this is that a smaller d f value represents more resistance in the complex induced fractures.  Figure 5 shows that the mass fractal dimension of induced fractures (Hausdorff index) has a significant effect on the pressure behavior at almost all the stages, except for the wellbore storage stage. Overall, the smaller the mass fractal dimension is, the larger the gas flow resistance is and the greater the pressure consumption is with the same production. As can be seen, the locations of the type curves are higher with a smaller df. The reason for this is that a smaller df value represents more resistance in the complex induced fractures.   Figure 6, the adsorption factor mainly influences the position of the type curves at the inter-porosity flow stage. A larger adsorption factor represents a stronger adsorption and production capacity and therefore makes the "concave" appear wider and deeper on the type curves. Figure 7 shows the effect of the  6-8 demonstrate the influences of the adsorption factor, apparent permeability coefficient, and inter-porosity flow coefficient on the type curves of MFHWs. As shown in Figure 6, the adsorption factor mainly influences the position of the type curves at the inter-porosity flow stage. A larger adsorption factor represents a stronger adsorption and production capacity and therefore makes the "concave" appear wider and deeper on the type curves. Figure 7 shows the effect of the apparent permeability coefficient on the transient pressure response. The apparent permeability has a similar effect to that of the inter-porosity coefficient in Figure 8. The total seepage and diffusion ability of the shale matrix is represented by the apparent permeability coefficient. The smaller the apparent permeability coefficient or inter-porosity coefficient is, the later the "depression" appears on the type curves.
Energies 2019, 12 12 of 18 apparent permeability coefficient on the transient pressure response. The apparent permeability has a similar effect to that of the inter-porosity coefficient in Figure 8. The total seepage and diffusion ability of the shale matrix is represented by the apparent permeability coefficient. The smaller the apparent permeability coefficient or inter-porosity coefficient is, the later the "depression" appears on the type curves.   apparent permeability coefficient on the transient pressure response. The apparent permeability has a similar effect to that of the inter-porosity coefficient in Figure 8. The total seepage and diffusion ability of the shale matrix is represented by the apparent permeability coefficient. The smaller the apparent permeability coefficient or inter-porosity coefficient is, the later the "depression" appears on the type curves.    Figure 9 shows the impact of the stress sensitivity factor on the pressure-transient response of MFHWs. It can be seen that stress sensitivity affects the whole flow stage, and it has a greater impact in the late time period. The reason for this is that the pressure drop becomes greater in the late time period. The greater the stress sensitivity is, the higher the positions of the pressure and pressure derivative curves are. This depicts the weaker seepage capacity. As shown in Figure 10, when all the factors are changed at the same time from a smaller parameter group ① to a larger parameter group ②, the positions of type curves for parameter group ② are obviously lower than the positions of type curves for parameter group ① . This indicates that when all the factors become larger, the final pressure drop becomes smaller. The reason for this is that most factors with greater values, such as , , , , and , can have positive  Figure 9 shows the impact of the stress sensitivity factor on the pressure-transient response of MFHWs. It can be seen that stress sensitivity affects the whole flow stage, and it has a greater impact in the late time period. The reason for this is that the pressure drop becomes greater in the late time period. The greater the stress sensitivity is, the higher the positions of the pressure and pressure derivative curves are. This depicts the weaker seepage capacity.  Figure 9 shows the impact of the stress sensitivity factor on the pressure-transient response of MFHWs. It can be seen that stress sensitivity affects the whole flow stage, and it has a greater impact in the late time period. The reason for this is that the pressure drop becomes greater in the late time period. The greater the stress sensitivity is, the higher the positions of the pressure and pressure derivative curves are. This depicts the weaker seepage capacity. As shown in Figure 10, when all the factors are changed at the same time from a smaller parameter group ① to a larger parameter group ②, the positions of type curves for parameter group ② are obviously lower than the positions of type curves for parameter group ① . This indicates that when all the factors become larger, the final pressure drop becomes smaller. The reason for this is that most factors with greater values, such as , , , , and , can have positive Figure 9. Effect of the stress sensitivity factor on type curves.
As shown in Figure 10, when all the factors are changed at the same time from a smaller parameter group 1 to a larger parameter group 2 , the positions of type curves for parameter group 2 are obviously lower than the positions of type curves for parameter group 1 . This indicates that when all the factors become larger, the final pressure drop becomes smaller. The reason for this is that most factors with greater values, such as F CD , α, d f σ m , β t , and λ, can have positive effects by making the pressure consumption smaller, and only γ * D has the opposite influence on pressure and pressure derivatives. effects by making the pressure consumption smaller, and only * has the opposite influence on pressure and pressure derivatives.

Case Study
This section shows an application of the presented model in a fractured horizontal well (A1) of an actual shale gas field in the Sichuan basin, which has 12 fractures evenly distributed along its horizontal wellbore. The depth of well A1 is 880 m and the thickness of the shale layer is 76 m. The production was 2400 cubic meters per day for 16 h, and then it was shut down for 73 h during the pressure build-up test. For more details, refer to the related literature [16]. After transferring the build-up testing data to dimensionless forms, the actual log-log curves were plotted.
As shown in Figure 11, the improved five-region flow model proposed in this work was applied to match the build-up testing data and was able to perfectly match the real testing data by adjusting the relevant parameters. The results of the interpretation are listed in Table 3. The results reveal that hydraulic fracturing greatly increases the permeability of the fractured zone and produces complex induced fractures with fractal features.

Case Study
This section shows an application of the presented model in a fractured horizontal well (A1) of an actual shale gas field in the Sichuan basin, which has 12 fractures evenly distributed along its horizontal wellbore. The depth of well A1 is 880 m and the thickness of the shale layer is 76 m. The production was 2400 cubic meters per day for 16 h, and then it was shut down for 73 h during the pressure build-up test. For more details, refer to the related literature [16]. After transferring the build-up testing data to dimensionless forms, the actual log-log curves were plotted.
As shown in Figure 11, the improved five-region flow model proposed in this work was applied to match the build-up testing data and was able to perfectly match the real testing data by adjusting the relevant parameters. The results of the interpretation are listed in Table 3. The results reveal that hydraulic fracturing greatly increases the permeability of the fractured zone and produces complex induced fractures with fractal features.

Conclusions
In order to describe the flow retardation in complex fractures in a way that considers the SRV region with anomalous diffusion and fractal features, an improved five-region model was established in this work by introducing the time-fractional flux law. Based on the present model, type curves of pressure and pressure derivative without wellbore storage were plotted and five flow stages were identified: bilinear flow, first linear flow, inter-porosity and fractal-anomalous flow, second linear flow, and boundary control flow. The sensitivity analysis revealed that fractal-anomalous diffusion has a significant impact on pressure-transient behaviors. When the anomalous diffusion exponent decreased from 1 to 0.75, which indicates Darcy flow changing to anomalous diffusion, the pseudo-pressure had less depletion at the early linear flow stages, but this subsequently became greater. When the Hausdorff index changed from 2 to 1.8, greater pressure consumption was needed to achieve the same production. Additionally, stress sensitivity, absorption, and Knudsen diffusion showed non-negligible influences on the pressure-transient response. These effects cannot be ignored. Therefore, the typical five-region flow model which does not take the fractal-anomalous diffusion into account cannot be applied for heterogeneous multi-fractured systems. The present model can be used to provide a more accurate and appropriate interpretation of well-testing data to guide exploration and development.

Conflicts of Interest:
The authors declare no conflict of interest.