A multi-scale flow model for production performance analysis in shale gas reservoirs with fractal geometry

Shale gas reservoirs can be divided into three regions, including hydraulic fracture regions, stimulating reservoir volume regions (SRV regions), and outer stimulating reservoir volume regions (OSRV regions). Due to the impact of hydraulic fracturing, induced fractures in SRV regions are often irregular. In addition, a precise description of secondary fractures in SRV regions is of critical importance for production analysis and prediction. In this work, the following work is achieved: (1) the complex fracture network in the SRV region is described with fractal theory; (2) a dual inter-porosity flow mechanism with sorption and diffusion behaviors is considered in both SRV and OSRV regions; and (3) both multi-rate and multi-pressure solutions are proposed for history matching based on fractal models and Duhamel convolution theory. Compared with previous numerical and analytic methods, the developed model can provide more accurate dynamic parameter estimates for production analysis in a computationally efficient manner. In this paper, type curves are also established to delineate flow characteristics of the system. It is found that the flow can be classified as six stages, including a bi-linear flow regime, a linear flow regime, a transition flow regime, an inter-porosity flow regime from the matrix to the fractures in the inner region, inter-porosity flow regime from matrix to fractures in the outer region, and a boundary dominant flow regime. The effects of the fracture and matrix properties, fractal parameters, inter-porosity flow coefficients, and sorption characteristics on type curves and production performance were studied in detail. Finally, production performance was analyzed for Marcellus and Fuling shale gas wells, in the U.S.A. and China, respectively.

Scientific RepoRts | (2018) 8:11464 | DOI: 10.1038/s41598-018-29710-1 as region 2 (inner region); and flow in the OSRV region defined as region 3 (outer region). To establish mathematical models, the following basic assumptions are made: (1) Darcy flow occurs in hydraulic fractures; (2) Random fractures in SRV regions are considered using fractal geometry; (3) The OSRV region could be homogenous or double porosity with sorption behaviors; (4) Sorption characteristics are considered in the rock matrix; (5) The compressibility factor and viscosity change with pressure; (6) The Z-factor changes with pressure; (7) Knudsen diffusion and slippage effects are considered in the matrix model. Model establishment. Different from conventional reservoirs, shale gas seepage often exhibits strong nonlinearity 8 . In addition, physical parameters, such as the viscosity, compressibility factor and deviation factor, are functions of pressure. However, these methods are effective for only specific problems. The viscosity, compressibility and deviation factor of shale gas change with pressure, which leads to a non-linear seepage equation. The pseudo-time and pseudo-pressure [26][27][28] are defined and a detailed deduction is provided in Appendix A of the Supplementary File. These definitions can be written in the following forms: Pseudo-pressure: where p is pressure, Pa; μ is viscosity, Pa▪s; c ti is initial total compressibility, Pa −1 ; c t is total compressibility, Pa −1 ; t is time, s; Z is deviation factor; and P i is initial pressure, Pa. The spherical element was assumed to flow in the matrix, and a comprehensive continuity equation with sorption can be formulated as: where V i is sorption volume, m 3 /m 3 ; ρ im is gas density in the matrix, kg/m 3 ; ρ sc is standard condition gas density, kg/m 3 ; v im is velocity in the matrix, m/s; t a is pseudo-time, s. i = o, i can represent the outer reservoir (region 3) or inner reservoir (region 2), respectively. Pseudo-time t a is defined in Eq. 2. The total compressibility in shale has been previously presented 29,30 and can be expressed as: where c g is gas compressibility, Pa −1 ; S gi is gas saturation; c ep is intermediate variables; c s is sorption compressibility, Pa −1 ; c f is formation compressibility, Pa −1 ; c w is water compressibility, Pa −1 ; S wi is initial water saturation; B g is gas volume coefficient; V L is Langmuir volume, m 3 /kg; p L is Langmuir pressure, Pa; ρ B is density of shale, kg/m 3 ; and ϕ is porosity. Accounting for gas slippage and gas diffusion, gas velocity in the matrix can be expressed as: where c im is matrix compressibility, Pa −1 . The initial condition and boundary conditions can be given as: where k of is the fracture permeability in the outer region, m 2 ; ϕ of is fracture porosity in the outer region; k om is matrix permeability in the outer region, m 2 ; x f is fracture half-length, m; x e is boundary length, m; c oft is total compressibility in the outer region, Pa −1 ; and m of is fracture pseudo-pressure in the outer region, Pa. Eqs 10 through 12 constitute the natural fracture model of the OSRV region. Natural fracture model of the SRV region (region 2). The change in porosity and permeability in the inner region are large and heterogeneous. Therefore, the fractal model can be utilized to describe SRV behavior. Furthermore, we work exclusively in Cartesian coordinates. Thus, the fractal relationships are given as 17,18 : Assuming one-dimensional flow in the y direction, the diffusivity equation is given by: Substituting Eq. 13 into Eq. 14 yields: where k If is the fracture permeability in the inner region, m 2 ; ϕ If is fracture porosity in the inner region; k im is matrix permeability in the inner region, m 2 ; x f is fracture half-length, m; c If is fracture compressibility in the inner region, Pa −1 ; c Ift is total compressibility of the inner region, Pa −1 ; m If is fracture pseudo-pressure of the inner region, Pa; d is fractal dimension representing the dimension of fractal fracture network embedded in the Euclidean matrix; and θ is connectivity index characterizing the diffusion process.  Hydraulic fracture model (region 1). The flow in hydraulic fractures can be expressed as: Inner boundary condition can be given as: Outer boundary condition can be expressed as: where Q is flow rate, m 3 /s; c F is fracture compressibility, Pa −1 ; w is fracture width, m; ϕ F is hydraulic fracture porosity; and k F is the fracture permeability, m 2 .

Solution of Mathematical Model
Bottom-hole pressure solution at constant rate conditions. A constant-rate solution can be obtained by using the Laplace transform method (Appendix B of the supplementary file), and we provide the solution of the non-linear difference equation using pseudo-pressure and pseudo-time definitions. The dimensionless definitions of the model and dimensionless solution at a constant rate condition are given as: In the above results, one-dimensional linear flow in the hydraulic fracture is considered. Flow choking within the fracture must be considered by using the following equation that was presented by Mukherjee and Economides 31 : S C is flow skin, and expressed as: Bottom-hole rate solution at constant pressure conditions. Note that the constant-pressure solution and constant-rate solution conform to the relation suggested by Van-Everdingen and Hurst 32 : Therefore, the rate solution at a constant pressure condition can be obtained as: The solutions at variable rate or variable pressure conditions. During the process of production, both the rate and pressure change with the time. Constant rate and constant pressure solutions cannot be applied to a real reservoir. Through using convolution theory 29,33 , the variable rate and variable pressure solutions can be given by Calculation process. The material balance equation must be established to calculate pseudo-time and pseudo-pressure. The modified gas compressibility factor and material balance equation that were proposed by Moghadam 30 for shale gas reservoirs were adopted in this paper. The modified gas compressibility factor Z** is defined as: where the parameters c ep and c d are defined in Eq. 4. The modified material balance equation can be written as: The average pressure can be determined by taking advantage of Eqs 32 and 33, and pseudo-pressure and pseudo-time can be calculated. The entire calculation process is illustrated in Fig. 3, which shows the procedures of variable rate solution and variable pressure solution. The procedures above could be applied to interpret production data from shale gas reservoirs.

Validation of solution.
The accuracy of the coupled model was validated by comparing our results with those of an analytical solution using IHS Harmony software. The settings of the reservoir and fracture properties are listed in Table 1. We selected special cases of fractal parameters d = 2 and θ = 0, which represents constant permeability and constant porosity, as the fractal geometry has been not considered in commercial software. The total analysis time is 255 d.  Fig. 4. It is apparent that the production rate and cumulative gas production obtained from the proposed model agree very well with those from the analytical solution. A validation of the pressure solution under a variable rate condition is shown in Fig. 5. Figure 5 presents five cycles of open and shut well conditions. The pressure solution associated with a variable rate and the rate solution associated with variable pressure can be calculated using Eqs 30 and 31, respectively. The same parameters are listed in Table 1, and the production rate value for open well is set to 2 MMscf. Fractal parameters remain as d = 2 and θ = 0. The solution of the pressure under variable rate conditions is consistent with the solution from the commercial software IHS Harmony.

Results and Discussion
Flow characteristics analysis of type curves. Flow characteristics analysis of pressure and pressure derivative curves is shown in Fig. 6 and the basic data used are shown in Table 2. It is apparent that the fluid flow in shale gas reservoirs can be classified into six stages:   In this regime, the process of inter-porosity flow is terminated, and the pressure between the matrix and fractures have increased to a state of dynamic balance.
Effect of sensitive factors on flow characteristics. Figure 7 shows the influence of parameters, including the fractal dimension d; connectivity index θ; transfer coefficients of the inner region and outer region, λ I and λ O , respectively; storage coefficients of the inner region and outer region, ω I and ω O , respectively; and conductivity of fracture and region on pressure and pressure derivative curves. The same data are also listed in Table 2.
Influence of fractal geometry parameters on type curves. Figure 7a presents   significant impact on flow characteristics within all periods. Dimensionless pressure p D reflects pressure depletion during the production periods. From the pressure curve in Fig. 7a, it is indicated that a large fractal dimension will result in a small pressure depletion, which implies that the production rate will be enhanced under the same pressure drop. In addition, fracture morphology becomes more complex when fractal dimension d becomes large. Figure 7b presents the effect of the connectivity index of natural fractures in SRV region θ on plots of type curves.
The θ values are set to 0, 0.4, and 0.8. Different from the fractal dimension, a large connectivity index will result in a large pressure depletion, which indicates that the production rate will be reduced under the same pressure drop. Moreover, a larger θ value leads to worse connectivity in the formation. In other words, pressure depletion is proportional to the connectivity index of natural fractures in the SRV region and is inversely proportional to fractal dimension. Figure 7c through d show the effect of the inter-porosity coefficient between the matrix and natural fractures on type curves. It reflects gas flow ability from the matrix to natural fractures in the inner region. Similarly, the inter-porosity flow coefficient of outer region λ o is proportional to matrix permeability k om of the outer region, and is inversely proportional to fracture permeability k of of the outer region. Thus, it reflects the gas flow ability from the matrix to natural fractures in the outer region. The mechanism of the theory reveals that a large inter-porosity flow coefficient reflects relatively strong matrix flow ability and weak fracture flow ability. Figure 7d presents the effect of the inter-porosity flow coefficient λ o on type curves. The values are set to 0.03, 0.3, and 3. It is seen from Fig. 8d that a large λ o value corresponding to strong matrix flow ability will lead to a small pressure depletion as the matrix compensates for pressure loss by supplying natural fractures. From the derivative curves of Fig. 7c and d, it is observed that flow in intermediate time is mainly affected by the inter-porosity flow coefficient. Figure 7e and f show the effect of the storage coefficient of natural fractures on the type curves. The ω I values reflect natural fracture energy, which was set to 0.000001, 0.1, and 0.5 in the inner region. The ω I values reflect natural fracture energy and are set to 0.000001, 0.00001, and 0.0001in the inner region. It is known from the pressure curves of both figures that a large storage coefficient will lead to a small pressure depletion due to the storage of rich gas resources in the natural fractures. This results in an extension of the linear flow duration, which is indicated from the derivative curves of both figures. From the derivative curves of Fig. 7e and f, it could also be concluded that flow in early and medium time is mainly influenced by the storage coefficient.  Influence of fracture conductivity and region conductivity on type curves. Figure 7g shows the effect of the hydraulic fracture conductivity on the type curves. The C fD values reflect the flow ability of shale gas in hydraulic fractures and are set to 1, 5, and 20. It is observed from the pressure curves that a small fracture conductivity will lead to a large pressure depletion, which implies that the production rate will be enhanced by improving fracture conductivity. However, the production rate is improved only at an early time. It is apparent from the pressure derivative curves that there is no effect on flow in middle and late time. Figure 7h presents the effect of the inner region conductivity C RD on type curves. The C RD values indicate flow ability of shale gas in the inner region and are set to 1, 2, and 5. C RD has an important impact on fluid flow during all production periods. According to the definition of Eq. B-29, C RD is proportional to the fracture permeability of the inner region and is inversely proportional to the fracture permeability of the outer region. The pressure curves in Fig. 7h show that a large region conductivity will lead to a large pressure depletion. It is apparent from the derivative curves that the shapes of the curves do not change as the inner region conductivity changes in the bi-linear flow period.

Influence of matrix transfer and storage parameters on type curves.
Effect of sensitive factors on production performance. Figure 8 shows the effects of different parameters, including fractal dimension d, connectivity index θ, Langmuir volume V L , Langmuir pressure P L , fracture half-length x f , fracture conductivity C fD , inter-porosity flow coefficient of inner region and outer region λ I and λ O , respectively, and storage coefficients of the inner region and outer region ω I and ω O , respectively, on production rate and cumulative production. The data used are listed in Table 3. The bottom-hole pressure was set to 2,000 psi in all cases.  Figure 8a presents the effect of fractal dimension d on plots of production rate and cumulative production versus time. A larger d represents a more complex structure of the fractal natural fractures. It is observed that the fractal dimension d has a significant impact on the solutions throughout the entire production period. A larger d value leads to a high production rate and a high cumulative production. In addition, as time increases, the rate first declines rapidly and then declines gradually. Figure 8b presents the effect of the connectivity index of natural fractures in the SRV region θ on plots of production rate and cumulative production versus time. The θ values were set to 0, 0.01, and 0.2. The connectivity index reflects connectivity and tortuosity between natural fractures. It is apparent from Fig. 8b that the connectivity index of natural fractures in the SRV region has a significant impact on the solutions throughout the entire production period. As shown in Fig. 8b, a larger θ value leads to a low production rate and low cumulative production. This point implies that a larger θ reflects more complicated tortuosity, which will lead to worse connectivity in the formation. Figure 8c presents the effect of the Langmuir volume on production behavior. The Langmuir volume V L values are set to 100 scf/ton, 200 scf/ton, and 300 scf/ton. It is apparent from the figure that a larger V L value will result in a high production rate and cumulative production. However, the increase in cumulative production is not obvious. The cumulative productions associated with V L = 100 scf/ton, 200 scf/ton, and 300 scf/ton are G p = 622.35 MMscf, 653.76 MMscf, and 682.26 MMscf, respectively. The fact that shale gas adsorption volume is directly proportional to the Langmuir volume will cause the shale gas adsorption volume to decrease as the Langmuir volume increases. Figure 8d presents the effect of the Langmuir pressure p L on plots of production rate and cumulative production versus time. The p L values are set to 100 psi, 500 psi, and1200 psi. Simulation results show that production rate and cumulative production will increase as the Langmuir pressure increases. The fact that shale gas adsorption volume is inversely proportional to the Langmuir pressure will cause an increasing shale gas adsorption volume as the Langmuir pressure increases.

Sorption.
Fracture parameters. Figure 8e shows the effect of fracture half-length on production behavior. The fracture half-length x f was set to 141 ft, 171 ft, and 201 ft. It is apparent that a large x f will lead to a high production rate and cumulative production. A larger x f value has an important impact on enhancing production rate and  Fracture conductivity also has an important impact on enhancing production rate and cumulative production. Figure 8f presents the effect of fracture conductivity on production behavior. The fracture conductivity C fD values were set to 0.1, 1, and 10. It is apparent that a large C fD value will lead to a high production rate and cumulative production. The cumulative production of C fD = 0.1, 1, and 10 are G p = 88.48 MMscf, 272.70 MMscf, and 583.31 MMscf, respectively. In the early stage, effectively increasing fracture conductivity effectively can remarkably improve the production rate and cumulative production. Figure 8g presents the effect of the inter-porosity flow coefficient of the inner region on production behavior. The fracture inter-porosity flow coefficient λ I values are set to 0.1, 1, and 10. It is apparent from the figure that a large λ I value will lead to a high production rate and cumulative production. Figure 8h presents the effect of the storage coefficient of the inner region on production behavior. The fracture storage coefficient ω I values are set to 0.4, 0.7, and 1. It is apparent from the figure that a large ω I value will lead to a high production rate and cumulative production.

History Matching of Production Data
Automatic history matching. Automatic parameter estimation (APE) is a mathematical process, known as multi-variable optimization, which automatically adjusts a specified set of function parameters to minimize error between the function and measured data. Minimizing the mathematical function is called the objective function. With analytical modelling, the objective function is calculated as follows.
If the Calculate Pressure mode is used: If the Calculate Rate mode is used: where p calc is calculated pressure for the i-th point; p hist is historical pressure for the i-th point; q calc is calculated rate for the i-th point; q hist is historical rate for the i-th point; n hist is historical point number; and E avg reflects the difference between the calculated and historical data. Therefore, a smaller E avg indicates a better history match. The Simplex routine 34 is a non-linear regression algorithm used for APE for reservoir and well parameters (k If , k of , x f , C fD , S C , etc.) when modelling pressure and rate transient data. It requires only function evaluations of the objective function, and not the derivatives. Modification of the downhill Simplex method to achieve greater convergence is accomplished by imposing constraints on the parameters during the search. Estimates of the parameters are always checked against preset maximum and minimum values for each parameter. Once the routine has converged on some parameters, it is restarted with a slight perturbation away from the final values and is allowed to converge again. This ensures that the parameter estimates found are not the result of some local minimum in the residual, but rather a more global minimum.
Field examples. Case 1: Marcellus shale, U.S.A. Using our model, history matching of production data was performed based on field data. The well is a multi-stage fractured horizontal well with 12 fractures in a typical shale gas reservoir in the Marcellus shale gas field, U.S.A. The pressure data and rate data were acquired from production tests after almost one-year of production. Figure 9a shows the pressure data and rate data from Nov. 24 th 2009 to Aug 6 th , 2010. The basic input parameters are summarized in Table 4. We selected the fracture permeability of the inner region, the fracture permeability of the outer region, choking skin, fracture half-length, fracture conductivity, fractal dimension, inter-porosity flow coefficient of inner region, and storage coefficient of inner region as uncertain parameters of history matching, which are listed in Table 4 and marked in bold. Figure 9b presents the historical matching results of production data from Marcellus shale gas wells using our proposed model. Our model fits the pressure and production rate very well. The presented model can also be applied to production performance analysis in complex conditions. Table 4 gives the initial estimation of uncertain parameters in shale gas, and the final matching results are shown in Fig. 9b. The k If value is 0.00588 mD, the k Of value is 0.000463 mD, the S C value is 0.089, the x f value is 279.086 ft, the C fD value is 15.1856, the d value is 1.96849, the λ I value is 1.85974, and the ω I is 0.128578. The relative error of our model for all of the data is E avg = 13.58%.
With the reservoir and fracture parameters determined by history matching, the production performance of the well can be predicted. In this case, the BHP was set to 631 psi from the final production data point. Figure 9c presents the predicted gas production and cumulative production. The cumulative production of the shale gas multi-fractured well with fractal geometry is estimated to be 1381.87 MMscf and the rate of well will become 0.01219 MMscf/d after 20 years. In the first five years, the rate of the well decreases rapidly, and cumulative production of the well increases dramatically. However, in the next 15 years, the rate of the well gradually decreased, and the cumulative production of the well increased slowly. From the fifth year to the twentieth year, the rate of well decreases only from 0.081876 MMscf/d to 0.01219 MMscf/d and cumulative production of well increases Scientific RepoRts | (2018) 8:11464 | DOI:10.1038/s41598-018-29710-1 only from 1241.7 MMscf to 1381.87 MMscf, which does not constitute a remarkable improvement. In other words, it is not economically feasible to continue developing this shale gas well for up to 20 years.
Case 2: Fuling shale, China. The second well is also a multi-stage fractured horizontal well with 21 fractures in a typical shale gas reservoir in the Fuling high pressure shale gas field, China. The pressure and production rate data were acquired from production test performed over 474 d. Figure 10a presents the pressure and production rate profile from Jan. 7 th 2016 to Apr. 24 th , 2017. The basic input parameters are summarized in Table 5. The uncertain parameters of history matching are shown in Table 5 and are marked in bold.
As shown in Fig. 10b, our model matched the BHP pressure and production rate very well from the Fuling shale gas field. Table 5 provides an initial estimation of uncertain parameters in the Fuling shale gas field, and the final matching results are shown in Fig. 10b. The k If value is 0.00273 mD, the k Of value is 2.77e-5 mD, the S C value is 0.107, the x f value is 127.232 ft, the C fD value is 12.6512, the d value is 2.01746, the λ I value is 0.28484, and the ω I is 0.0287237. The relative error of our model for all of the data is E avg = 4.16%.
According to the results, the production performance from Fuling shale gas wells can be predicted. In this case, if the BHP is set at 6948.28 psi, Fig. 10c shows the predicted gas production and cumulative production. The cumulative production of the shale gas multi-fractured well with fractal geometry is estimated to be 1001.43 MMscf after 20 years. In the first four years, the rate of the well decreases rapidly. However, in the next 16 years, the rate of well becomes approximately stable.

Conclusions
In this study, a multi-scale flow model in shale gas reservoirs with fractal geometry was presented to investigate pressure transient response and analyze production performance, through which we examined the effects of shale rock properties and fractal characteristics. From the above analysis, the following conclusions are made: (1) The flow characteristics of type curves can be divided into six regimes: bi-linear flow regime, linear flow regime, transition flow regime, inter-porosity flow regime from the matrix to fractures in the inner region, inter-porosity flow regime from matrix to fractures in the outer region, and boundary dominant flow regime, respectively. (2) A large fractal dimension represents a complex structure of fractal natural fractures and leads to a small pressure depletion. A large connectivity index, which reflects connectivity and tortuosity among natural fractures, results in a large pressure depletion. (3) A large inter-porosity flow coefficient of the outer region and inner region corresponding to strong matrix flow ability will lead to a small pressure depletion because the matrix compensates for pressure loss by supplying natural fractures. A large storage coefficient results in a small pressure depletion due to rich gas resources that are stored by natural fractures. (4) A small fracture conductivity and a large region conductivity leads to a large pressure depletion. (5) A large fractal dimension, small connectivity index, large Langmuir volume, large Langmuir pressure, large fracture conductivity, large inter-porosity flow coefficient, and large storage coefficient can enhance the production rate and cumulative production of shale gas wells. (6) Based on a downhill Simplex algorithm, multiple uncertain parameters can be interpreted well by using the presented multi-rate solutions and multi-pressure solutions. The presented fractal model is well validated by matching real field examples.