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

ABSTRACT 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. SIMILAR CONTENT BEING VIEWED BY OTHERS EFFECT OF PORE-THROAT STRUCTURE ON MOVABLE

FLUID AND GAS–WATER SEEPAGE IN TIGHT SANDSTONE FROM THE SOUTHEASTERN ORDOS BASIN, CHINA Article Open access 05 March 2025 DEVELOPMENT OF COUPLED FLUID-FLOW/GEOMECHANICS MODEL CONSIDERING

STORAGE AND TRANSPORT MECHANISM IN SHALE GAS RESERVOIRS WITH COMPLEX FRACTURE MORPHOLOGY Article Open access 20 August 2024 RESEARCH ON THE MICROSCOPIC PORE-THROAT STRUCTURE AND RESERVOIR

QUALITY OF TIGHT SANDSTONE USING FRACTAL DIMENSIONS Article Open access 01 October 2024 INTRODUCTION The successful exploitation of shale gas heavily depends on the combination of horizontal

drilling, completions, and fracturing technology1,2. It has been proven that multi-stage fracturing horizontal wells constitute a very effective way to exploit low-permeability shale gas

reservoirs. Due to the complex fracture network created in stimulated reservoir volume regions, it is a great challenge to accurately analyze production performance and evaluate

post-fracture performance in such complex fractured reservoirs3. Since 1972, many analytical and semi-analytical methods on production and flow behaviors in conventional gas and shale gas

reservoirs have been utilized to simulate pressure transient behavior for horizontal wells with hydraulic fractures. Soliman _et al_.4 analyzed fluid flow mechanisms for multi-fractured

horizontal wells and presented an approximate model to investigate the influence of fracture number, fracture orientation, fracture direction, fracture conductivity, and horizontal well

location on flow behaviors. Larsen and Hegre5 presented two kinds of models for fractured horizontal wells with multiple fractures with finite conductivity: circular fractures with radial

flow and vertical fractures with linear flow. For circular fractures, a cylindrical coordinate system is used with the perpendicular axis of the fracture coinciding with the z axis and the

center of the fracture located at z = 0. The fracture size is given by the width, the well radius rw and fracture outer radius rf. For a vertical fracture, a Cartesian coordinate system is

used with the x axis chosen as vertical axis, the y axis coinciding with the axis of the wellbore, and the center of the fracture located at z = 0. The fracture size is given by the width,

the fracture half-length xf away from the wellbore, and the fracture half-length yf along the wellbore. Guo and Evans6 presented a randomly distributed vertical fractures model to predict

production performance, and generated type curves to estimate reservoir properties and fracture characteristics. Bin _et al_.7 presented an improved analytical solution for fractured

horizontal wells with different fracture intensities to investigate the effect of fracture properties on flow characteristics. However, due to constant rate assumption, the model cannot be

applied to interpret production data. Wang _et al_.8 proposed a semi-analytical model in tight oil reservoirs with fractal geometries, which was used to describe complex pore sizes and

random fractures in porous media. Since sorption characteristics were not incorporated into their model, the model cannot consider the effect of sorption on production. A technique was

introduced by Rbeawi and Tiab9 to interpret the pressure transient behavior of fractured horizontal wells. The distribution of hydraulic fractures could be longitudinal or transverse,

vertical or inclined, symmetrical or asymmetrical. Brown _et al_.10 and Stalgorova and Mattar11 proposed a classical tri-linear flow model to simulate production performance of fractured

horizontal wells in unconventional reservoirs. However, complex geometries and sorption characteristics were not taken into account in these studies, so that they were not suitable for

production performance analysis in shale gas reservoirs. The above studies4,5,6,7,8,9,10,11 mainly focused on fluid flow in primary fractures, and ignored fluid flow inside secondary

fractures and matrix. Due to the influence of induced fractures, uniform fracture models were not applicable to describe heterogeneous SRV zones. Some scholars attempted to utilize numerical

models with a discrete fracture network to describe secondary fractures in SRV zones12,13,14,15. Although the production performance of a complex fracture network can be discerned through

the use of numerical methods, numerical simulation is less attractive and pragmatic due to the large amount of requisite knowledge and time-consumption during the history matching process3.

In contrast, the analytical model is an alternative for accurately forecasting flow behaviors and calculation efficiency. Fractal geometry has shown potential in the analysis of flow and

transport properties in porous media. Katz and Thompson16 were the pioneers to find experimental evidence that indicated that the pore spaces of a set of sandstone samples were fractals and

self-similar over three to four orders of magnitude in length. Fractal theory can also describe complex fracture shape in shale gas and is easy to apply to obtain analytical solutions. Chang

and Yortsos17 introduced fractal theory into petroleum engineering, and reported that the permeable fractures embedded within the matrix would be arranged in a fractal dimension. Cossio _et

al_.18 presented a model with fractal geometry to develop a rapid and accurate semi-analytical solution for flow in a single vertical fracture that fully penetrates a homogeneous

infinite-acting reservoir. Yao _et al_.19 proposed a fractal double-porosity model for the transient flow analysis of fluids. Kong _et al_.20 developed the basic formulae of seepage

velocity, permeability, and porosity in both porous and fractured fractal media. Wang _et al_.3 presented semi-analytical modelling of flow behavior in fractured media with fractal geometry.

Zhang _et al_.21 proposed an approach for estimating the size and equivalent permeability of fractal SRV zones for vertically fractured wells in tight reservoirs. The proposed fluid flow

models with fractal geometry3,16,17,18,19,20,21 are convenient for conducting pressure transient analysis for understanding flow in complex reservoirs22,23,24,25. Random fractures will be

produced during the process of hydraulic fracturing. Fractal theory can exhaustively describe the randomness of pore sizes and gas reservoir heterogeneity24. However, most of the previous

models are only involved in the aspects of pressure transient analysis in tight reservoirs, while shale multi-scale flow characteristics are not considered4,5,6,7,8,9,10,11,12. Second, gas

PVT is assumed to be constant, so that the final results will lead to an assessment of errors. However, gas PVT changes with pressure in the development of shale gas reservoirs. Third, the

characteristics of fractal geometry in SRV regions are not reflected in their models4,5,6,7,8,9,10,11,12. In this paper, an improved tri-linear flow model in shale gas reservoirs with

fractal geometry is presented to analyze and predict production performance of multi-fractured wells. First, multi-scale flow mechanism in shale gas reservoirs was considered in the

analytical model (Fig. 1). Second, the main fractures were considered into longitudinal rectangular fractures and secondary fractures in the SRV zone were described using fractal theory.

Third, gas PVT changed with reservoir pressure in this paper, and a multi-rate solution with variable pressure was also proposed to interpret real well production data measured under various

production systems. Finally, the downhill Simplex method was employed to form the history matching procedure. METHODOLOGY CONCEPTUAL MODEL Figure 1 presents the conceptual model of

stimulated reservoir volume (SRV) with fractal geometry. Many authors have presented analytical models to investigate regular fracture networks, such as the single planar hydraulic fractures

network and orthogonal hydraulic fractures network4,5,6,7,8,9,10,11,12,13,14,15. A schematic diagram of the presented model is shown in Fig. 2. Shale gas flow can be divided into three flow

regions: flow in hydraulic fractures defined as region 1 (hydraulic fracture region); flow in the SRV region defined 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 nonlinearity8. 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-pressure26,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: $$m(p)=\frac{{\mu }_{i}{Z}_{i}}{{p}_{i}}{\int }_{0}^{p}\frac{p}{\mu Z}\,dp$$ (1) Pseudo-time:

$${t}_{a}={\int }_{0}^{{\rm{t}}}\frac{{\mu }_{gi}{c}_{ti}}{{\mu }_{g}(p){c}_{t}(p)}dp$$ (2) where _p_ is pressure, Pa; _μ_ is viscosity, Pa▪s; cti 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. MATRIX MODEL Matrix inter-porosity flow occurs in both SVR regions and OSRV

regions, as shown in Fig. 2. The spherical element was assumed to flow in the matrix, and a comprehensive continuity equation with sorption can be formulated as:

$$-\frac{1}{{r}_{im}^{2}}\frac{\partial }{\partial {r}_{im}}({r}_{im}^{2}{\rho }_{im}{v}_{im})=\frac{\partial ({\rho }_{im}{\varphi }_{im})}{\partial {t}_{a}}+{\rho }_{sc}\frac{\partial

{V}_{i}}{\partial {t}_{a}}$$ (3) where _V_i is sorption volume, m3/m3; _ρ_im is gas density in the matrix, kg/m3; _ρ_sc is standard condition gas density, kg/m3; _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 presented29,30 and can be expressed as: $${c}_{t}=({c}_{g}[{S}_{gi}-{c}_{ep}]-\frac{\partial {c}_{ep}}{\partial p}+{c}_{d})$$ (4&x02010;1)

$${c}_{ep}=({c}_{f}+{c}_{w}{S}_{wi})({p}_{i}-p)$$ (4&x02010;2) $${c}_{d}=\frac{{\rho }_{B}{B}_{g}{V}_{L}}{\varphi }(\frac{{p}_{L}}{{({p}_{L}+p)}^{2}})$$ (4&x02010;3) 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, m3/kg; _p_L is Langmuir pressure, Pa; _ρ_B is density of shale, kg/m3; and

_φ_ is porosity. Accounting for gas slippage and gas diffusion, gas velocity in the matrix can be expressed as: $${v}_{im}=-\,{k}_{im}{F}_{iD}\frac{\partial {p}_{im}}{\partial

{r}_{im}}=-\,{k}_{im}[(1+8\alpha {k}_{n})+{\mu }_{g}{D}_{g}{c}_{g}/{k}_{im}]\frac{\partial {p}_{im}}{\partial {r}_{im}}$$ (5) The equation of gas state is given as: $${\rho

}_{im}=\frac{{p}_{im}M}{nZRT}$$ (6) Substituting Eq. 1, and Eqs 4 through 6 into Eq. 3 produces the governing equation of matrix: $$\frac{1}{{r}_{im}^{2}}\frac{\partial }{\partial

{r}_{im}}({r}_{im}^{2}{k}_{im}[(1+8\alpha kn)+\frac{{\mu }_{g}{D}_{g}{c}_{g}}{kim}]\frac{\partial {m}_{im}}{\partial {r}_{im}})={\varphi }_{im}{c}_{imt}\frac{\partial {m}_{im}}{\partial

{t}_{a}}$$ (7) and $${c}_{imt}={c}_{g}+{c}_{im}+{c}_{id}$$ (8) where _c_im is matrix compressibility, Pa−1. The initial condition and boundary conditions can be given as:

$${{m}_{im}|}_{{t}_{a}=0}={m}_{i}({p}_{i})$$ (9&x02010;1) and $${{m}_{im}|}_{{r}_{im}=R}={m}_{if}$$ (9&x02010;2) $${\frac{\partial {m}_{im}}{\partial {r}_{im}}|}_{{r}_{im}=0}=0$$

(9&x02010;3) Equation 7 through 9 constitute the matrix seepage governing equation in both the SRV and OSRV regions. The dimensionless equations and solution method are provided in

Appendix B. NATURAL FRACTURE MODEL OF THE OSRV REGION (REGION 3) According to the above definitions of pseudo-pressure and pseudo-time, assuming one-dimensional flow in the _x_ direction, as

shown in Fig. 3, the diffusivity equation and the associated boundary conditions for the outer reservoir are given by: $$\frac{{\partial }^{2}{m}_{of}(p)}{\partial {x}^{2}}=\frac{{\varphi

}_{of}{c}_{oft}{\mu }_{gi}}{{k}_{of}}\frac{\partial {m}_{of}(p)}{\partial {t}_{a}}+\frac{3}{R}\frac{{k}_{om}{F}_{oD}}{{k}_{of}}{\frac{\partial {m}_{om}(p)}{\partial

{r}_{om}}|}_{{r}_{om}=R}$$ (10) $${(\frac{\partial {m}_{of}(p)}{\partial x})}_{x={x}_{e}}=0$$ (11) $${m}_{of}{(p)}_{x={x}_{f}}={m}_{If}{(p)}_{x={x}_{f}}$$ (12) where _k_of is the fracture

permeability in the outer region, m2; _φ_of is fracture porosity in the outer region; _k_om is matrix permeability in the outer region, m2; _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 as17,18:

$${k}_{If}^{Fractal}(y)={k}_{If}{(\frac{y}{{x}_{f}})}^{d-\theta -2}={k}_{If}{y}_{D}^{d-\theta -2}$$ (13&x02010;1) $${\varphi }_{If}^{Fractal}(y)={\varphi

}_{If}{(\frac{y}{{x}_{f}})}^{d-2}={\varphi }_{If}{y}_{D}^{d-2}$$ (13&x02010;2) Assuming one-dimensional flow in the _y_ direction, the diffusivity equation is given by:

$$\begin{array}{rcl}\frac{\partial }{\partial y}({k}_{If}^{Fractal}\frac{\partial {m}_{If}(p)}{\partial y})+{(\frac{y}{{x}_{f}})}^{d-2}\frac{{k}_{If}}{{x}_{f}}({\frac{\partial

{m}_{If}(p)}{\partial x}}_{x={x}_{f}}) & = & {(\frac{y}{{x}_{f}})}^{d-2}{\varphi }_{If}{c}_{Ift}\mu \frac{\partial {m}_{If}(p)}{\partial {t}_{a}}\\ & &

+\,{(\frac{y}{{x}_{f}})}^{d-2}\frac{3}{R}\frac{{k}_{Im}{F}_{ID}}{1}{\frac{\partial {m}_{\text{Im}}(p)}{\partial {r}_{Im}}|}_{{r}_{Im}=R}\end{array}$$ (14) Substituting Eq. 13 into Eq. 14

yields: $$\begin{array}{l}{(\frac{y}{{x}_{f}})}^{d-\theta -2}{k}_{If}[\frac{{\partial }^{2}{m}_{If}}{\partial {y}^{2}}+\frac{d-\theta -2}{y}\frac{\partial {m}_{If}}{\partial

y}]+{(\frac{y}{{x}_{f}})}^{d-2}\frac{{k}_{If}}{{x}_{f}}({\frac{\partial {m}_{If}(p)}{\partial x}}_{x={x}_{f}})\\ \,={(\frac{y}{{x}_{f}})}^{d-2}{\varphi }_{If}{c}_{Ift}\mu \frac{\partial

{m}_{If}(p)}{\partial t}+\,{(\frac{y}{{x}_{f}})}^{d-2}\frac{3}{R}\frac{{k}_{Im}{F}_{ID}}{1}{\frac{\partial {m}_{\text{Im}}(p)}{\partial {r}_{Im}}|}_{{r}_{Im}=R}\end{array}$$ (15) The

boundary conditions for the inner reservoir are given by: $${(\frac{\partial {m}_{If}}{\partial y})}_{y={y}_{e}}=0$$ (16&x02010;1) and $${{m}_{If}|}_{y=\frac{\omega

}{2}}={{m}_{F}|}_{y=\frac{\omega }{2}}$$ (16&x02010;2) where _k_If is the fracture permeability in the inner region, m2; _φ_If is fracture porosity in the inner region; _k_im is matrix

permeability in the inner region, m2; _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. Eqs 15 and 16 are the natural fracture model of the OSRV region, and the solution method is presented in Appendix B. HYDRAULIC

FRACTURE MODEL (REGION 1) The flow in hydraulic fractures can be expressed as: $$\frac{\partial ({k}_{F}\frac{\partial {m}_{F}(p)}{\partial x})}{\partial x}+{\frac{2{k}_{If}}{\omega

}\frac{\partial {m}_{If}(p)}{\partial y}|}_{y=\frac{\omega }{2}}={\varphi }_{F}{c}_{F}{\mu }_{gi}\frac{\partial {m}_{F}(p)}{\partial {t}_{a}}$$ (17) Inner boundary condition can be given as:

$${\frac{\partial {m}_{F}}{\partial x}|}_{x=0}=\frac{Q{\mu }_{gi}}{2wh{k}_{F}}$$ (18&x02010;1) Outer boundary condition can be expressed as: $${\frac{\partial {m}_{F}}{\partial

x}|}_{x={x}_{f}}=0$$ (18&x02010;2) where _Q_ is flow rate, m3/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, m2. 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: $${\bar{m}}_{FD}(p)=\frac{\pi }{{C}_{FD}s\sqrt{{\alpha }_{F}}}\frac{\cosh [\sqrt{{\alpha

}_{F}}(1-{x}_{D})]}{\sinh (\sqrt{{\alpha }_{F}})}$$ (19) where $${\alpha }_{F}=2{\beta }_{F}/{C}_{FD}+s/{\eta }_{FD}$$ (20) $${\beta }_{F}=-\,{(\frac{{w}_{D}}{2})}^{\gamma -1}\sqrt{{\alpha

}_{{\rm{o}}}}\frac{{K}_{n-1}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{y}_{eD}^{\gamma }]{I}_{n-1}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{(\frac{{w}_{D}}{2})}^{\gamma

}]-{I}_{n-1}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{y}_{eD}^{\gamma }]{K}_{n-1}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{(\frac{{w}_{D}}{2})}^{\gamma }]}{{I}_{n}[\frac{\sqrt{{\alpha }_{o}}}{\gamma

}{(\frac{{w}_{D}}{2})}^{\gamma }]{K}_{n-1}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{y}_{eD}^{\gamma }]+{I}_{n-1}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{y}_{eD}^{\gamma

}]{K}_{n}[\frac{\sqrt{{\alpha }_{o}}}{\gamma }{(\frac{{w}_{D}}{2})}^{\gamma }]}$$ (21) $${\alpha }_{o}=\frac{{\beta }_{o}}{{C}_{RD}{y}_{eD}}+{u}_{I}$$ (22) $${\beta

}_{o}=\sqrt{{u}_{o}}\,\tanh [\sqrt{{u}_{o}}({x}_{eD}-1)]$$ (23) $${u}_{I}=s[f(s)]=s[\frac{{\omega }_{I}}{{\eta }_{ID}}+\frac{{\lambda }_{I}{\eta }_{ID}}{5s}({\tau }_{I}\,\coth \,{\tau

}_{I}-1)]$$ (24) $${u}_{o}=s[g(s)]=s[\frac{{\omega }_{o}}{{\eta }_{oD}}+\frac{{\lambda }_{o}{\eta }_{oD}}{5s}({\tau }_{o}\,\coth \,{\tau }_{o}-1)]$$ (25) 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

Economides31: $${\bar{m}}_{FD}(p)=\frac{\pi }{{C}_{FD}s\sqrt{{\alpha }_{F}}}\frac{\cosh [\sqrt{{\alpha }_{F}}(1-{x}_{D})]}{\sinh (\sqrt{{\alpha }_{F}})}+\frac{{S}_{C}}{s}$$ (26) SC is flow

skin, and expressed as: $${S}_{C}=\frac{{k}_{If}h}{{k}_{F}{w}_{F}}[\mathrm{ln}(\frac{h}{2{r}_{w}})-\frac{\pi }{2}]$$ (27) 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 Hurst32: $${\bar{m}}_{FD}\times {\bar{q}}_{FD}=\frac{1}{{s}^{2}}$$ (28)

Therefore, the rate solution at a constant pressure condition can be obtained as: $${\bar{q}}_{FD}=\frac{{C}_{FD}\sqrt{{\alpha }_{F}}}{s\pi }\frac{\sinh (\sqrt{{\alpha }_{F}})}{\cosh

[\sqrt{{\alpha }_{F}}(1-{x}_{D})]}$$ (29) 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 theory29,33, the variable rate and variable pressure solutions can be given by

$${{\rm{m}}}_{wD}({t}_{Dn})=\sum _{i=1}^{n-1}\frac{q({t}_{i})}{q({t}_{n})}[{{\rm{m}}}_{FD}({t}_{Dn}-{t}_{Di-1})-{{\rm{m}}}_{FD}({t}_{Dn}-{t}_{Di})]+{{\rm{m}}}_{FD}({t}_{Dn}-{t}_{Dn-1})$$

(30) $${q}_{wD}({t}_{Dn})=\sum _{\begin{array}{c}i=1\\ n > 1\end{array}}^{n-1}{{\rm{m}}}_{FDi}[{q}_{FD}({t}_{Dn}-{t}_{Di-1})-{q}_{FD}({t}_{Dn}-{t}_{Di})]+{q}_{FD}({t}_{Dn}-{t}_{Dn-1})$$

(31) 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 Moghadam30 for shale gas reservoirs were adopted in this paper. The modified gas compressibility factor Z** is defined as: $${Z}^{\ast \ast

}=\frac{{\rm{p}}}{[\frac{1}{{S}_{gi}}\frac{p}{Z}({S}_{gi}-{c}_{ep}-{c}_{d})+\frac{{p}_{i}}{{Z}_{i}}(\frac{G}{{G}_{f}}-1)]\frac{{G}_{f}}{G}}$$ (32) where the parameters _c__ep_ and _c__d_ are

defined in Eq. 4. The modified material balance equation can be written as: $$\frac{p}{{Z}^{\ast \ast }}=\frac{{p}_{i}}{{Z}_{i}^{\ast \ast }}(1-\frac{{G}_{p}}{G})$$ (33) 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. Figures 4 and 5 present the validation results by comparing the analytical results based on the proposed

fractal model with the analytical solution from commercial software10. The algorithm from IHS Harmony software was presented by Brown _et al_.10 in 2009. A validation of the rate solution

under a constant pressure condition of 2,000 psi is shown in 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: * (1) The bi-linear flow regime is a straight line with a slope of 1/4 on the pressure derivative curve. The flow of gas occurs simultaneously both within hydraulic fractures and

near them. * (2) The linear flow regime is a straight line with a slope of 1/2 on the pressure derivative curve. In this period, linear flow around the fractures occupies the dominant

position. * (3) The transition flow regime (a) shows an indefinite slope line on the pressure derivative curve. * (4) The first inter-porosity flow regime is a regime of supplementation from

the shale matrix to the fracture system in the inner region. The derivative curve exhibits the first V-shaped segment. * (5) The second inter-porosity flow regime represents the

inter-porosity flow from the matrix system to the fracture system in the outer region. The pressure derivative curve shows the second V-shaped segment. * (6) The boundary dominant flow

regime is a straight line with a slope of 1 on the pressure derivative curve. 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 the effect of fractal dimension _d_ on plots of type curves, wherein fractal dimension d values are set to 1.6, 2, and 2.4. A larger _d_ represents a more

complex structure of fractal natural fractures. It is concluded that the fractal dimension _d_ has a 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. INFLUENCE OF MATRIX TRANSFER AND STORAGE PARAMETERS ON TYPE CURVES 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 kom of the outer region, and is inversely proportional to fracture permeability kof 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 CfD 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. 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. FRACTAL GEOMETRY The fractal dimension values are set to 1.99, 2, and 2.01. 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. SORPTION 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 Gp = 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. 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 cumulative production. The cumulative production of _x__f_ = 141 ft, 171 ft, and 201 ft

are Gp = 570.35 MMscf, 653.76 MMscf, and 734.77 MMscf. 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 Gp = 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. MATRIX TRANSFER AND STORAGE COEFFICIENTS 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: $${E}_{{\rm{avg}}}=\frac{{\sum

}_{i}(\frac{|{{\boldsymbol{(}}{p}_{calc}{\boldsymbol{)}}}_{i}-{{\boldsymbol{(}}{p}_{hist}{\boldsymbol{)}}}_{i}|}{{{\boldsymbol{(}}{p}_{hist}{\boldsymbol{)}}}_{i}})}{{n}_{hist}}$$ (34) If the

Calculate Rate mode is used: $${E}_{{\rm{avg}}}=\frac{{\sum

}_{i}(\frac{|{{\boldsymbol{(}}{{\rm{q}}}_{calc}{\boldsymbol{)}}}_{i}-{{\boldsymbol{(}}{{\rm{q}}}_{hist}{\boldsymbol{)}}}_{i}|}{{{\boldsymbol{(}}{{\rm{q}}}_{hist}{\boldsymbol{)}}}_{i}})}{{n}_{hist}}$$

(35) 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 routine34 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. 24th 2009 to Aug 6th, 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 Eavg = 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 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. 7th 2016

to Apr. 24th, 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 Eavg = 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. REFERENCES * Cipolla, C. L. Modeling production and evaluating fracture performance in unconventional gas reservoirs. _J. Pet. Technol._ 9, 84–90 (2009).

Article  Google Scholar  * Mayerhofer, M. J., Lolon, E. P. & Warpinski, N. R. What is stimulated reservoir volume. _SPE Prod. Oper._ 2, 89–98 (2011). Google Scholar  * Wang, J. L., Wei,

Y. S. & Qi, Y. D. Semi-analytical modeling of flow behavior in fractured media with fractal geometry. _Transp Porous Med._ 112, 707–736 (2016). Article  MathSciNet  CAS  Google Scholar 

* Soliman, M. Y., Hunt, J. L. & El Rabaa, W. Fracturing aspects of horizontal wells. _J. Pet. Technol._ 42(8), 966–973 (1990). Article  Google Scholar  * Larsen, L. & Hegre, T. M.

Pressure transient analysis of multi-fractured horizontal wells. SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana. Society of Petroleum Engineers.

https://doi.org/10.2118/28389-MS (1994, September 25–28). * Guo G. L. & Evans R. D. Pressure-transient behavior and inflow performance of horizontal wells intersecting discrete

fractures. SPE Annual Technical Conference and Exhibition, Houston, Texas. Society of Petroleum Engineers. https://doi.org/10.2118/26446-MS (1993, October 3–6). * Bin, Y., Su, Y. L.,

Rouzbeh, G. M. & Zhen, H. R. A new analytical multi-linear solution for gas flow toward fractured horizontal wells with different fracture intensity. _Journal of Natural Science and

Engineering._ 23, 227–238 (2015). Article  Google Scholar  * Wang, W. D., Shahvali, M. & Su, Y. L. A semi-analytical fractal model for production from tight oil reservoirs with

hydraulically fractured horizontal wells. _Fuel._ 158, 612–618 (2015). Article  CAS  Google Scholar  * AI Rbeawi, S. & Tiab, D. Transient pressure analysis of horizontal well with

multiple inclined hydraulic fractures using type-curve matching. SPE International Symposium and Exhibition on Formation Damage Control, Lafayette, Louisiana, USA. Society of Petroleum

Engineers. https://doi.org/10.2118/149902-MS (2012, February 15–17). * Brown, M. L., Ozkan, E., Raghavan, R. S. & Kazemi, H. Practical solutions for pressure transient responses of

fractured horizontal wells in unconventional reservoirs. SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana. Society of Petroleum Engineers.

https://doi.org/10.2118/125043-MS (2009, October 4–7). * Stalgorova, E. & Mattar, L. Practical analytical model to simulate production of horizontal wells with branch fractures. SPE

Canadian Unconventional Resources Conference, Calgary, Alberta, Canada. Society of Petroleum Engineers. https://doi.org/10.2118/162515-MS (2012, 30 October–1 November). * Wu, Y. S.,

Ehlig-Economides, C. & Qin, G. A triple-continuum pressure-transient model for a naturally fractured vuggy reservoir. SPE Annual Technical Conference and Exhibition, Anaheim, California,

U.S.A. Society of Petroleum Engineers. https://doi.org/10.2118/110044-MS (2007, November 11–14). * Monifar, A. M., Varavei, A. & Johns, R. T. Development of a coupled dual continuum and

discrete fracture model for the simulation of unconventional reservoirs. SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA. Society of Petroleum Engineers.

https://doi.org/10.2118/163647-MS (2013, February 18–20). * Cipolla, C. L., Fitzpatrick, T. & Williams, M. J. Seismic-to-simulation for unconventional reservoir development. SPE

Reservoir Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE. Society of Petroleum Engineers, https://doi.org/10.2118/146876-MS (2011, October 9–11). * Weng, X.,

Kresse, O. & Cohen, C. Modeling of hydraulic-fracture-network propagation in a naturally fractured formation. _SPE Prod. Oper._ 26(4), 368–380 (2011). CAS  Google Scholar  * Katz, A. J.

& Thompson, A. H. Fractal sandstone pores - implications for conductivity and pore formation. _Physical Review Letters._ 54(12), 1325–1328 (1985). Article  ADS  PubMed  CAS  Google

Scholar  * Chang, J. & Yortsos, Y. Pressure transient analysis of fractal reservoirs. _SPE Formation Evaluation._ 5(1), 31–38 (1990). Article  Google Scholar  * Cossio, M., Moridis, G.

J. & Blasingame, T. A. A semi-analytic solution for flow in finite-conductivity vertical fractures by use of fractal theory. _SPE J._ 2, 83–96 (2013). Article  Google Scholar  * Yao, Y.

D., Wu, Y. S. & Zhang, R. L. The transient flow analysis of fluid in a fractal, double-porosity reservoir. _Transp. Porous Med._ 94, 175–187 (2012). Article  MathSciNet  CAS  Google

Scholar  * Kong, X. Y., Li, D. L. & Lu, D. T. Transient pressure analysis in porous and fractured fractal reservoirs. _Sci. China Ser. E-Tech. Sci._ 52(9), 2700–2708 (2009). Article 

MATH  Google Scholar  * Zhang, Q., Su, Y. L., Wang, W. D. & Sheng, G. L. A new semi-analytical model for simulating the effectively stimulated volume of fractured wells in tight

reservoirs. _Journal of Natural Science and Engineering._ 27, 1834–1845 (2015). Article  Google Scholar  * Beier, R. A. Pressure transient field date showing fractal reservoir structure.

CIM/SPE International Technical Meeting, Calgary, Alberta, Canada. Society of Petroleum Engineers. https://doi.org/10.2118/21553-MS (1990, June 10–13). * Aprilian, S., Abdassah, D. &

Mucharan, L. Application of fractal reservoir model for interference test analysis in Kamojang geothermal field (Indonesia). _SPE Form. Eval._ 27(4), 7–14 (1993). Google Scholar  * Cai, J.

C. & Yu, B. M. A discussion of the effect of tortuosity on the capillary imbibition in porous media. _Trans. Porous Media._ 89(2), 251–263 (2011). Article  CAS  Google Scholar  *

Moinfar, A., Varavei, A., Sepehrnoori, K. & Johns, R. T. Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs.

_SPE J._ 19(2), 289–303 (2014). Article  Google Scholar  * Kikani, J. & Pedrosa, O. A. Perturbation analysis of stress-sensitive reservoirs. _SPE Formation Evaluation._ 6(3), 379–386

(1991). Article  CAS  Google Scholar  * Pedrosa, O. A. Pressure transient response in stress-sensitive formations. SPE California Regional Meeting, Oakland, California. Society of Petroleum

Engineers. https://doi.org/10.2118/15115-MS (1986, April 2–4). * Ansah, J., Knowles, R. S. & Blasingame, T. A. A Semi-analytic (p/z) rate-time relation for the analysis and prediction of

gas well performance. _SPE Reservoir Evaluation & Engineering._ 3(6), 525–533 (2000). Article  CAS  Google Scholar  * Wang, L., Dai, C. & Xue, L. A semi-analytical model for pumping

tests in finite heterogeneous confined aquifers with arbitrarily shaped boundary. _Water Resour Res_. https://doi.org/10.1002/2017WR022217 (2018). * Moghadam, S., JEJE, O. & Mattar, L.

Advanced gas material balance, in simplified format. _Journal of Canadian Petroleum Technology._ 50(1), 1–10 (2009). Google Scholar  * Mukherjee, H. & Economides, M. J. A parametric

comparison of horizontal and vertical well performance. _SPE Formation Evaluation._ 6(2), 209–216 (1991). Article  Google Scholar  * Van Everdingen, A. F. & Hurst, W. The application of

the Laplace transformation to flow problems in reservoirs. _J Petrol Technol._ 1(12), 305–324 (1949). Article  Google Scholar  * Stehfest, H. Numerical inversion of Laplace transforms.

_Commun. ACM._ 13(1), 47–49 (1970). Article  Google Scholar  * Nelder, J. A. & Mead, R. A simplex method for function minimization. _Computer Journal._ 7, 308–313 (1965). Article 

MathSciNet  MATH  Google Scholar  Download references ACKNOWLEDGEMENTS This work is partially funded by the National Natural Science Foundation of China (Grant Nos U1663208 and 51520105005),

the State Major Science and Technology Special Project of China during the 13th Five-Year Plan (Grant Nos 2016ZX05014-004, 2016ZX05025-003-007, 2016ZX05034001-007), and China Postdoctoral

Science Foundation (Grant no. 2017M620527). AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * ERE & BIC-ESAT, College of Engineering, Peking University, Beijing, 100871, China Lei Wang, 

Xiang Li & Zunyi Xia * Energy Innovation Software Co. Ltd., Beijing, 100094, China Zhenzhen Dong & Xiang Li Authors * Lei Wang View author publications You can also search for this

author inPubMed Google Scholar * Zhenzhen Dong View author publications You can also search for this author inPubMed Google Scholar * Xiang Li View author publications You can also search

for this author inPubMed Google Scholar * Zunyi Xia View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS L. Wang established the novel

mathematical model. L. Wang and Z. Dong wrote the main manuscript and prepared all the figures. X.L. provided production data in oil fields. Z.X. designed the computer program. Both X.L. and

Z.X. supervised the work and revised the manuscript. All of authors discussed the results and critically reviewed the manuscript. CORRESPONDING AUTHORS Correspondence to Lei Wang or Xiang

Li. ETHICS DECLARATIONS COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION PUBLISHER'S NOTE: Springer Nature remains neutral with regard to

jurisdictional claims in published maps and institutional affiliations. ELECTRONIC SUPPLEMENTARY MATERIAL SUPPLEMENTARY INFORMATION RIGHTS AND PERMISSIONS OPEN ACCESS This article is

licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give

appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in

this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative

Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a

copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Wang, L., Dong, Z., Li, X. _et al._ A multi-scale flow

model for production performance analysis in shale gas reservoirs with fractal geometry. _Sci Rep_ 8, 11464 (2018). https://doi.org/10.1038/s41598-018-29710-1 Download citation * Received:

06 April 2018 * Accepted: 17 July 2018 * Published: 30 July 2018 * DOI: https://doi.org/10.1038/s41598-018-29710-1 SHARE THIS ARTICLE Anyone you share the following link with will be able to

read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing

initiative