- Email: [email protected]
Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir
Journal of Hydrology 512 (2014) 447–456
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir Yu-Long Zhao ⇑, Lie-Hui Zhang, Jian-Xin Luo, Bo-Ning Zhang State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, PR China
a r t i c l e
i n f o
Article history: Received 30 October 2013 Received in revised form 25 February 2014 Accepted 10 March 2014 Available online 19 March 2014 This manuscript was handled by Peter K. Kitanidis, Editor-in-Chief, with the assistance of Markus Tuller, Associate Editor Keywords: Fractured horizontal well Unconventional gas reservoir Stimulated reservoir volume Well testing model Performance analysis
s u m m a r y This paper extended the conventional multiple hydraulic fractured horizontal (MFH) well into a composite model to describe the stimulated reservoir volume (SRV) caused by hydraulic fracturing. Employing the Laplace transform, Source function, and Dirac delta function methods, the continuous linear source function for general composite dual-porosity is derived, and the solution of the MFH well in a composite gas reservoir is obtained with the numerical discrete method. Through the Stehfest numerical algorithm and Gauss elimination method, the transient pressure responses for well producing at a constant production rate and the production rate vs. time for constant bottomhole pressure are analyzed. The effects of related parameters such as natural permeability and radial of the SRV region, formation permeability and interporosity coefficient on transient pressure and production performance are analyzed as well. The presented model and obtained results in this paper not only enrich the well testing models of such unconventional reservoir, but also can use to interpret on-site data which have significance on efficient reservoir development. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction Horizontal well with massive hydraulic fractures has been proved to be a common way to produce hydrocarbons from ultra-low permeability reservoirs, such as shale gas reservoir, and the objective of it has not only been to create several highconductivity flow paths, but also activate and connect existing discrete or micro-seismic natural fractures so as to generate large fracture networks (Clarkson, 2013). The volume of a reservoir containing the fracture network and the main hydraulic fractures with which effectively stimulated to increase the well performance is defined as SRV (Ozkan et al., 2009, 2011; Stalgorova and Mattar, 2012b; Clarkson, 2013). As transient pressure analysis been an important method for petro-engineers to study the formation characteristics, many models have been proposed to describe the theory of fluid seepage flow from the formation to the well and many detailed analytical solutions are developed. Guo et al. (1994) developed methods for the performance prediction of multiple fractured horizontal wells without considering the effect of interface between fractures. Horne and Temeng (1995) accounted for the interference between multiple fractures by the superposition of influence function. ⇑ Corresponding author. Tel.: +86 13568841671. E-mail address: [email protected] (Y.-L. Zhao). http://dx.doi.org/10.1016/j.jhydrol.2014.03.026 0022-1694/Ó 2014 Elsevier B.V. All rights reserved.
Raghavan et al. (1997) developed a mathematical model to discern the characteristic response of multiple-fractured horizontal wells. Wan and Aziz (1999) developed general solution for horizontal wells with multiple fractures. Crosby et al. (2002) and Wan and Aziz (2002) derived semi-analytical solutions to respectively describe the pressure response for transverse fractures and fractures rotated at any angle to the horizontal well. Zerzar et al. (2003) combined the boundary element method and Laplace transformation to present a comprehensive solution for multiple vertical fractures horizontal wells. Al-Kobaisi and Ozkan (2004) presented a hybrid numerical-analytical model for the pressure transient response of horizontal wells intercepted by a vertical fracture. Rbeawi and Tiab (2013) introduced new analytical models to investigate the pressure behavior and flow regimes of a horizontal well with multiple hydraulic fractures regarding them as being either longitudinal or transverse, either vertical or inclined, and either symmetrical or asymmetrical. Zhao et al. (2012, 2013) used source function to analyze the transient pressure response of a MFH well in shale gas reservoir. Although the above researchers have made much works on transient pressure analysis of MFH well in both conventional and unconventional reservoir, most of the models they proposed are not take in account the effect of the induced fracture around the main hydraulic fractures (called the fractured network). If we want to describe the fracture network in well testing models similar
448
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
Nomenclature cg ct1m,ct1f ct2m,ct2f
ct2
C h k2 k1m, k1f k2m, k2f k2 DLfi L LH M M12 p1m, p1f p2m, p2f psc _ q qscL qfi qsc rm r R s
gas compressibility, MPa1 total compressibility in the matrix and natural fracture system of the inner region respectively, MPa1 total compressibility in the matrix and natural fracture system of the outer region respectively (for composite dual porosity model), MPa1 total compressibility of the outer reservoir(for inner natural fracture and outer homogeneous composite reservoir model), MPa1 wellbore storage coefficient, m3/MPa formation thickness, m reservoir permeability, D matrix and fracture network permeability in SRV region, D matrix and micro-fracture permeability in outer region, D formation permeability of the outer reservoir, D the length of the ith discrete fracture, m half horizontal well length, m horizontal well length, m molar mass of nature gas, mol/kg mobility ratio, dimensionless, defined in Table 1 matrix/fracture network system pressure in SRV, MPa matrix/micro-fracture system pressure in outer region, MPa pressure at standard condition, psc = 0.10325 MPa instantaneous line source strength, sm3 continuous line source flow rate, m3/h the production rate of the ith discrete fracture, m3/d well production rate, m3/d the radial distance of the SRV, m radius, m the universal gas constant, R = 0.008314 (MPa m3)/ (K kmol) Laplace variable
with the real geological situation, we should implement lots of CT scanning and core experience and then establish the mathematical models of the induced fractures one by one without any simplification. Although this idea is very good, the implementations of such huge works are impossible. In order to concisely describe the fracture network (natural or induced) around the hydraulic fractures, assuming uniform distribution of identical hydraulic fractures along the length of the horizontal well, Ozkan et al. (2009, 2011) utilizes the concept of Lee and Brockenbrough (1986) of tri-linear model with inner reservoir of naturally fractured to represent the MFH well performance in unconventional reservoirs (see diagram in Fig. 1A). Brown et al. (2011) presented an analytical tri-linear flow model that incorporates transient fluid transfer from matrix to fracture to simulate the pressure transient and production behavior of fractured horizontal wells in unconventional reservoirs. Using the tri-linear flow model, Apaydin et al. (2012) presented an analytical model with composite matrix blocks to describe the effect of matrix micro-fractures on the effective matrix permeability. With the assumption that the enhanced region occupies only part of the space between the fractures and neglects the region beyond the fractures, Stalgorova and Mattar (2012a) improved the tri-linear flow model to simulate fractures that is surrounded by a simulated region of limited extendibility. Later,
Skin t T Tsc v W12
skin factor, dimensionless production time, h temperature, K temperature at standard condition, Tsc = 293 K flow velocity, m/h storability ratio between inner and outer region, dimensionless, defined in Table zg gas deviation factor, sm3/m3 w pseudo-pressure, MPa2/cp / formation porosity, fraction £1m, £1f matrix and fracture network porosity in SRV, fraction £2m, £2f matrix and micro-fracture porosity, fraction q fluid density in reservoir condition, g/cm3 q1m fluid density in the matrix system, g/cm3 q1 fluid density in the fracture system, g/cm3 l gas viscosity, cp lgi gas viscosity at the initial reservoir condition, cp a1 sharp factor of the inner region, 1/m2 a2 sharp factor of the outer region, 1/m2 x1f storability coefficient of natural fracture network in SRV, defined in Table 1 x2f storability coefficient of natural fractures in outer region, defined in Table 1 k1mf inter-porosity coefficient from matrix into natural fracture network in inner region, defined in Table 1 k2mf inter from matrix into natural fractures in outer region, defined in Table 1 c1, c2, f1(s), f2(s) variables groups Subscript D
dimensionless
Superscript – Laplace transform
Stalgorova and Mattar (2012b) extended the tri-linear flow model to a five regions flow model (see diagram in Fig. 1B), which can not only simulate the partial simulated reservoir volume but also take into account the fluid supplement beyond the fractures. Although the merits of the ‘‘tri-linear’’ and ‘‘five-flow’’ models are obvious, there are still some demerits of them. Both assume that all fractures have the same length and conductivity and are spaced uniformly along the horizontals, and the early-radial flow, later pseudo-radial flow period and interference phenomenon between the fractures flow cannot be observed, which makes the simplified linear flow models too simple. Therefore, we must take different geological conditions into consideration and propose different models to characterize the transient fluid flow into the MFH well with SRV. In this paper, we extend the conventional multi-fractured horizontal well (see Fig. 2b) into a novel composite model (see Fig. 2c), and derive a semi-analytical model to represent flow in an unconventional, ultra-tight formation, as described previously. The model is a dual-porosity composite system, consisting of a macro-fracture network (inner region, which compose of the natural and induced hydraulic fractures and represent the SRV) and a homogeneous matrix system (outer region). The MFH well locates in the center of the inner region with an arbitrary length and spacing (see Fig. 2a).
449
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
Region 3
Region 4
Hydraulic
Outer reservoir
Region 5
fracture
Inner reservoir naturally
Region 2
Region 1
Horizontal well
B—five-region model
A—tri-linear flow model
Fig. 1. Schematics of tri-linear and five-region flow model for MFH well.
Micro fractures
Micro fractures of
of formations
formations
Stimulation reservoir
Main hydraulically
volume (SRV)
fractures
A—(Complex geometry patterns)
SRV
Reservoir formation
B--(Conventional MFH model)
C—(Composite model in this paper)
Fig. 2. Hydraulic multiple fractured horizontal well with SRV and its presentation models.
2. Continuous line sources in composite dual-porosity gas reservoir To simulate the performance of stimulated reservoir volume in unconventional gas reservoir, we allow the SRV to be naturally fractured, and use the dual porosity idealization to represent the naturally fractured medium. Although the fluid transfer from the matrix to the fractures can be idealized by the method proposed
by Kazemi (1969) and de Swaan-O (1976), our choice is the pseudo-steady fluid transfer equation of Warren and Root (1963). Fig. 3 is a schematic of a continuous source located in the center of a composite dual-porosity gas reservoir. Although in some tight gas reservoirs, only the stimulated reservoir volume can be treated as a dual-porosity idealization, in order to derive the general source solution of this kind of reservoir, we treat the outer zone as the dual-porosity medium as well. This model can be applied
450
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
Fig. 3. Schematic of a continuous line source in a radial composite reservoir.
for these reservoirs with developed micro-fractures in shale and coaled methane gas reservoirs. If we further want to obtain the continuous source function with the homogeneous outer reservoir, the results can be easily derived from the former formula after simple transformation. The following assumptions are made during the derivation of the linear source function: (1) the reservoir is composite with the inner and outer region of dual-porosity, and the reservoir has a constant and uniform thickness with the upper and lower boundaries closed; (2) porosity and permeability is homogeneous in each formation medium system and region; and (3) fluid flow in the media follows the Darcy flow theory and neglects the gravitational and frictional effects (Zeng and Zhao, 2010).
2.2. Continuous line source function Utilizing the above governing equation for both inner and outer regions and combing the corresponding boundary conditions and initial conditions, the fully mathematical models of an instantaneous line source withdraw from the center of inner region can be established (see Appendix). The instantaneous line source function of composite reservoir with dual-porosity inner and outer reservoirs is obtained and the detailed derivations are presented in Appendix, and the dimensionless variables are defined in column 2 of Table 1. The instantaneous line source function can be obtained as _
Dw1f ¼ 2.1. Governing equations for this model
ð1Þ
1 @/ / @p
ð3Þ
1 @ q z @ p ¼ q @p p @p z
ð4Þ
Ignoring gravity effects, the Darcy’s Law in radial coordinate is given by 3
v ¼ 10
k @p l @r
Using the Duhamel theory (Everdingen and Hust, 1949), the continuous line source function can be obtained by
Dw1f ¼
The isothermal gas compressibility is
cg ¼
k1mf þ sx1f ð1 x1f ÞgrD grD ; k1mf þ sð1 x1f ÞgrD k2mf þ sx2f ð1 x2f Þ W12 ; grD ¼ f 2 ðsÞ ¼ k2mf þ sð1 x2f Þ M12
ð2Þ
The porosity compressibility is
cf ¼
M12 c1 K1 ðc1 r mD ÞK0 ðc2 r mD Þ c2 K0 ðc1 r mD ÞK1 ðc2 r mD Þ ; M12 c1 I1 ðc1 r mD ÞK0 ðc2 r mD Þ þ c2 I0 ðc1 r mD ÞK1 ðc2 r mD Þ qffiffiffiffiffiffiffiffiffiffiffi ffi qffiffiffiffiffiffiffiffiffiffiffiffi c1 ¼ sf1 ðsÞ; c2 ¼ sf2 ðsÞ AC ¼
f1 ðsÞ ¼
The real-gas law is simply the pressure/volume relationship predicted by a correction factor that accounts for the nonideal behavior of the gas, and the real-gas law is as follows
pM ¼ Z g qRT
ð6Þ
where
The governing equations of this composite reservoir can be derived by combining the continuity equation, the flux (Darcy) Law, and an equation of state (Badakhshan and Sobbi, 1995). The continuity equation in a homogeneous, cylindrical reservoir can be obtained by the mass conservation, which is
1 @ðr qv Þ 1 @ðq/Þ þ ¼0 r @r 3600 @t
psc T q k2f ½K ðc r Þ þ AC I0 ðc1 r D Þ pT sc lgi hð/2 ct2 Þfþm k1f L2 0 1 D
ð5Þ
For gas reservoir, both the gas viscosity l and compressibility cg are the function of reservoir pressure p, in order to make the problems more tractable, we make the assumptions that both l and cg are evaluated at the initial reservoir pressure pi.
psc T 1 k2f pT sc lgi hð/2 ct2 Þfþm k1f L2
Z
t
qscL ðsÞSc ðtD sÞds
ð7Þ
0
where Sc ¼ L1 ½Sc , Sc ¼ ½K0 ðc1 rD Þ þ AC I0 ðc1 rD Þ Taking the Laplace transform with tD in Eq. (7), the continuous line source function in inner region of composite dual-porosity reservoir with upper and lower boundaries closed is
~scL Dw1f ¼ q
psc T 1 ½K0 ðc1 r D Þ þ AC I0 ðc1 rD Þ T sc 86:4pk1f h
ð8Þ
Assuming that the line source with a constant rate qscL, we can define the dimensionless pseudo-pressure as follows
wD ¼
k2f hT sc 3:684 103 qscL psc T
Dw
ð9Þ
Combining Eq. (8) with Eq. (9), the continuous line source function is
w1fD ¼
1 ½K0 ðc1 r D Þ þ AC I0 ðc1 r D Þ M12 s
ð10Þ
451
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 Table 1 Definition of the dimensionless variables. Dimensionless variables Pseudo-pressure Dimensionless time
Both inner and outer region are dual-porosity Rp wl ¼ 2 0 l lpZg dp tD ¼ ð/
3:6k2f t
2 c t2 Þfþm
Dimensionless radial distance Storability ratio of fracture system in region 1
rD ¼
Inner region with dual-porosity and outer region is homogeneous Rp wl ¼ 2 0 l lpZ g dp 2t tD ¼ / 3:6k c l L2
lgi L2
2 t2
r L
gi
r D ¼ Lr
Þf x1f ¼ ð/ð/1 c1t1ct1Þfþm
Þf x1f ¼ ð/ð/1 c1t1ct1Þfþm
Storability ratio of fracture system in region 2
x2f ¼
Interporosity coefficient of region 1
k1mf ¼ a1 kk1m L2
k1mf ¼ a1 kk1m L2
Interporosity coefficient of region 2
k2mf ¼ 2 kk2m L2 2f M12 ¼ kk1f ==ll 2f ð/ c Þ W12 ¼ ð/12 ct1t2 Þfþm fþm
M12 ¼ kk1f2 ==ll
ð/2 ct2 Þf ð/2 ct2 Þfþm 1f
Mobility ratio between region 1 and region 2 Storability ratio between regions 1 and region 2 Dimensionless wellbore storage
1f
a
CD ¼
W12 ¼
C 2pð/2 ct2 Þfþm hL2
CD ¼
ð/1 ct1 Þfþm /2 ct2
C 2p/2 ct2 hL2
Note: when the inner region (SRV region) of reservoir is dual porosity and outer region is homogeneous, the continuous line source function is obtained by substituting the k2mf and x2f equal to 0 and 1 respectively. The naturally fracture’s permeability k2f will be replaced by k2.
When the outer reservoir does not have micro-fractures, the model may transform into an inner dual-porosity and outer homogeneous matrix system. With the same solving process described in Appendix and the new dimensionless variables defined in column 3 of Table 1, the source function can be obtained as
w1fD ¼
1 ½K0 ðc1 r D Þ þ AC I0 ðc1 r D Þ M12 s
ð11Þ
Where
AC ¼
pffiffi pffiffi pffiffi M12 c1 K1 ðc1 rmD ÞK0 ð sr mD Þ sK0 ðc1 r mD ÞK1 ð sr mD Þ pffiffi pffiffi M12 c1 I1 ðc1 r mD ÞK0 ð sr mD Þ þ sI0 ðc1 rmD ÞK1 ðr mD Þ
3.1. Analytical solutions The mathematical model assumes the wellbore to be intersected by N fractures, each fracture is assumed to have distinct property and can be at any position along the wellbore, and all the fractures are transverse to the well and fully-penetrating to the formation. We make the assumption that the flow from the reservoir to the wellbore sections is negligible comparing to the flow to the hydraulic fracture plane. The fractures are assumed to produce at a constant wellbore pressure (the horizontal well can always be treated as infinite in tight gas reservoir for low production rate). Giving that each half fracture is dispersed into MM sections of length DLfi, the total flow rate will be the summation of the fluid flow from each discrete fracture, which is:
qDi ¼ 1
wwDi ¼
MM2N X j¼1
qDfj wwDi;j s
Comparing to the mainly hydraulic fractures and wellbore, the formation permeability is too small, and the well production rate is usually very small for the unconventional gas reservoir. Therefore, we can treat them as infinite-conductivity, which will indicate zero pressure drops in them. So
ð12Þ
where qDi is the dimensionless discrete fracture rate, which is defined as
qfi DLfi qsc
2
wwD1;1
6 6 6 6 wwDk;1 6 6 6 6 4 wwD2N M;1 1
wwD1;k
wwD1;2N M
wwDk;k
wwDk;2N M
wwD2N M;k 1
wwD2N M;2N M
1
1
32
qD1
MM2N X Z tD j¼1
0
qDfj ðsÞw0wDi;j ðt D sÞds
ð14Þ
3
2
0
3
7 6 7 6 1 7 76 qD2 7 6 0 7 76 7 6 7 7 7 67 6 1 76 7¼6 7 7 7 67 6 1 76 7 6 7 76 7 6 7 1 54 qDM 2N 5 4 0 5 0
wwD
1 ð17Þ
In the above matrix, the symbol wwDi;j (both i and j vary in the range 1, 2N * MM) represents the effect of the production from element j at the location of element i. method to calculate the matrix elements for different configurations with the combination of continuous line source function can be referred to the paper of Zhao et al.(2012, 2013), which is 2
wwDi;j ¼
1 M12 sDLfD
Z j
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 2 2 xDi xD j a þ yDi ywD j þ7 6 K0 c1 6 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7da 6 j 4 2 2 5 AC I 0 c 1 xDi xD j a þ yDi ywD j
DLfD j 2 DLfD
2
ð18Þ
ð13Þ
The dimensionless pressure drop, pDi at element i in terms of pressure drops, pDi,j, at element i due to the production from element j = 1, . . . , MM * 2N, can be obtained by the following convolution:
wwDi ðtD Þ ¼
ð16Þ
There are 2N * MM + 1 equations when applying Eq. (15) at each element together with Eq. (12), which can solve the 2N * MM + 1 unknowns of wwD , qD1, qD2, . . . , qDM⁄2N1, qDM*2N. The matrix expression is
i¼1
qDi ¼
ð15Þ
wwD ¼ wwDi¼1;...;MM2N
3. Pressure behaviors for MFH well in composite gas reservoir
MM2N X
Assuming that the production rate in each discrete fracture is constant, in terms of Laplace transform, Eq. (14) can be written as
3.2. Computational consideration of the matrix coefficient According to the matrix described in Eq. (18), if the matrix coefficients are obtained, the dimensionless bottomhole pressure as well as the dimensionless production rate distribution in real time space along the fractures will be solved by the Gauss elimination method and Stehfest numerical inversion algorithm (Stehfest,
452
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
1970). To compute the coefficients efficiently, they can be calculated with different methods according to the observation point location i and the source element location j. If the observation and source element point locate at the same fracture, which means they have the same Y-axis, the calculation of Eq. (18) can be simplified as:
wwDi;j ¼
1 M12 sDLfDi
Z
DLfDj
2
DLfDj
ð19Þ
Using H0(x) to denote I0(x) and K0(x), the definite integration described in Eq. (19) can be simplified as follows, 8Rb Rb 1 > H0 ðaÞda 0 2 H0 ðaÞda if i > j > Z DLfDj 0 < 2 R b1 1 H ð c jx x a jÞd a ¼ 2 H ð a Þd a if i ¼ j 0 Di Dj 0 1 DLfDj 0 c1 > 2 > R b1 : R b2 H ð a Þd a H ð a Þd a if i < j 0 0 0 0 ð20Þ DLfDj where b1 ¼ c1 xDi xDj þ , b2 ¼ c1 xDi xDj 2 . Rx The definite integration of 0 H0 ðaÞda can be calculated with the numerical algorithm described by Milton and Irene (1964). When the observation point and the source element locate at different fractures, the result of Eq. (18) can be calculated directly by the numerical integration algorithm, such as adaptive algorithms, Gauss–Legendre and variable step size Runge–Kutta methods (Milton and Irene,1964; Escobar et al., 2013; Cai and Yu, 2011).
DLfDj 2
3.3. Transient pressure responses when well producing at a constant rate with considering wellbore storage and skin effects When we want to take into account the effects of wellbore storage as well as skin factor, using the Duhamel’s theorem (Everdingen and Hust, 1949), the dimensionless wellbore pressure in Laplace space becomes
swwDN þ Skin s þ C D s2 ðswwDN þ Skin Þ
ð21Þ
Combing Eq. (21) with Stehfest Numerical algorithm (Stehfest, 1970), the dimensionless pseudo-pressure and its derivative vs. time on log–log curves can be obtained. 3.4. Solution of well producing at a constant wellbore pressure When well producing at a constant bottomhole pressure, the dimensionless rate can be defined as follows:
qD ¼
qD ¼
1 s2 wwD
ð23Þ
The above equation is the solution for a multi-fractured horizontal well producing at a constant wellbore pressure.
½K0 ðc1 jxDi xDj ajÞ
2
þ AC I0 ðc1 jxDi xDj ajÞda
wwD ¼
Combing Eq. (22) with the expression of dimensionless bottom hole pressure in Table 1, the following equation can be obtained
3:684 103 psc qsc T kf hT sc Dw
ð22Þ
3.5. Discussion of the result We use the synthetic data in Table 2 to analyze the pressure response and the rate decline of the mathematical model derived in this paper, the results are shown in Figs. 4–13. Figs. 4–7 compare the differences in the transient pressure responses and the performances of a fractured horizontal well in homogeneous and fractured composite reservoirs with different parameters. The synthesis data are listed in Table 2. The solid and dash lines in these figures are obtained by the composite model with inner naturally fractured (SRV) and the outer homogeneous, and the dot lines with circle marker are obtained by the conventional model of MFH well in homogeneous reservoir (Raghavan et al., 1997; Wan and Aziz, 2002; Rbeawi and Tiab, 2013). Figs. 4 and 6 are the transient pressure responses of the well at a constant production rate. The effect of the fracture network permeability (k1f) in SRV region on type curves is shown in Fig. 4. It can be natural seen from the lines that the bigger k1f is, the smaller value of wwD and wwD * tD/CD will be in the early linear and early radial flow period in SRV region. Compared to the conventional MFH well model in homogeneous reservoir (Wan and Aziz, 2002; Rbeawi and Tiab, 2013), due to the existence of the SRV region, the pressure drop is smaller than that of the model without SRV when the well at a constant production rate. Fig. 6 shows the effect of interporosity coefficient (k1mf), which presents the capacity of the fluid transfer from the matrix to the natural fractures under same pressure difference. It can be seen that the bigger k1mf is, the early the transition flow between the matrix and natural fractures system will be, and also the corresponding pressure drop is much lower. Figs. 5 and 7 are the rate decline performance of the MFH well in composite reservoir under a constant bottomhole pressure. The comparison of the results of the homogeneous reservoir and the composite system with naturally fractured SRV displays significant contribution of the SRV to the production performance of the MFH. This contribution helps the well sustain a higher production rate during the all production flow period. In practical, the bigger the permeability of the fracture system is, the bigger production rate of the well is (as shown in Fig. 5). Fig. 7 shows the effect of interporosity coefficient (k1mf) between matrix and natural fracture
Table 2 Synthetic data used for the discussion of the results. Formation thickness, h, m Wellbore radius, rw, m Hydraulic fracture half-length, xF, m The radius of stimulated region, rm Bottomhole pressure, pwf, MPa Gas viscosity, l, cp Skin factor, skin
50 0.1 75 500 4 0.021 0.01
Inner region Matrix permeability, k1m, md Matrix porosity, £1m Total compressibility in matrix system, Ct1m, MPa1 Fracture permeability, k1f, md Fracture porosity, £1f Total compressibility in natural fracture system, Ct1f, MPa1 Interporosity coefficient, k1mf
Initial reservoir pressure, pi, MPa Horizontal well length, LH, m Number of fractures, N Reservoir temperature, T, K Molecular weight, Mg, kg/mol Wellbore storage coefficient, C, m3/MPa
35 1000 4 353 16 100
Outer region 0.0001 0.1 0.01 0.01 0.1 0.02 0.2
Permeability, k2, md Porosity, £2 Total compressibility, Ct2, MPa1
0.001 0.1 0.02
453
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 1E+1
1E+1
kkf=0.005mD f1 =0.005mD kkf=0.05mD f1 =0.05mD
Conventional model
1E-1
kk2=0.001mD 2=0.001mD
kk2=0.005mD 2=0.005mD
1E+0 ψwD, ψwD*tD/CD
ψwD, ψwD*tD/CD
1E+0
kk2=0.0001mD 2=0.0001mD
kkf=0.01mD f1 =0.01mD
1E-1
1E-2
1E-2 1E-3
pseudo-pressure, ψwD
pseudo-pressure, ψwD pseudo-pressure derivative, ψwD*tD/CD
pseudo-pressure derivative, ψwD*tD/CD 1E-3 1E-3
1E-2
1E-1
1E+0
1E+1 tD/CD
1E+2
1E+3
1E+4
1E-4 1E-3
1E+5
Fig. 4. Effect of induced fracture network permeability in SRV on pseudo-pressure and its derivative curves.
1E-2
1E-1
1E+0
1E+1 tD/CD
1E+2
1E+3
1E+4
1E+5
Fig. 8. Effect of outer reservoir permeability (k2) on pseudo-pressure and its derivative curves.
1000
100
kk2=0.001mD 2=0.001mD
kk2=0.0001mD 2=0.0001mD 100 k2=0.005mD k2=0.005mD q (104 m3/d)
q (104 m3/d)
10
10
1
1
0.1 1E-1
kkf=0.005mD f1 =0.005mD
kkf=0.01mD f1 =0.01mD
kkf=0.05mD f1 =0.05mD
Conventional model
1E+0
1E+1
0.1
1E+2
1E+3
0.01 1E-1
1E+4
1E+0
1E+1
1E+2
t (day)
λλ1m-f=1.00 1m-f =1.00
λλ1m-f=0.10 1m-f =0.10
Conventional model
1E-2
rrm=500m m=500m
rrm=750m m=750m
rrm=1000m m=1000m
rrm=1250m m=1250m
1E+8
1E-1
1E-2
pseudo-pressure, ψwD
pseudo-pressure, ψwD
pseudo-pressure derivative, ψwD*tD/CD 1E-2
1E-1
1E+0
1E+1
1E+2
1E+3
1E+4
pseudo-pressure derivative, ψwD*tD/CD
1E+5
1E-3 1E-3
1E-2
1E-1
1E+0
tD/CD
1E+1 1E+2 tD/CD
1E+3
1E+4
1E+5
1E+6
Fig. 10. Effect of SRV radius (rm) on pseudo-pressure and its derivative curves.
Fig. 6. Effect of k1mf on pseudo-pressure and its derivative curves.
100
100
rrm=500m m=500m
rrm=750m m=750m
rm=1000m rm=1000m
rrm=1250m m=1250m
10
q (104 m3/d)
q (104 m3/d)
10
0.1 1E-1
1E+7
1E+0 ψwD, ψwD*tD/CD
ψ wD, ψ wD*tD/CD
λλ1m-f=10.0 1m-f =10
1E-1
1
1E+6
1E+1
1E+1
1E-3 1E-3
1E+5
Fig. 9. Effect of outer reservoir permeability (k2) on well production performance.
Fig. 5. Effect of kf on well production performance.
1E+0
1E+3 1E+4 t (day)
1 λλ1m-f=0.10 1m-f =0.10
λλ1m-f=1.00 1m-f =1.00
λλ1m-f=10.0 1m-f =10
Conventional model
1E+0
1E+1
1E+2 t (day)
1E+3
1E+4
Fig. 7. Effect of k1mf on well production performance.
1E+5
0.1 1E-1
1E+1
1E+3
1E+5
1E+7
1E+9
t (day)
Fig. 11. Effect of SRV radius (rm) on well production performance.
454
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 1E+1 φ 0.05 1f=0.05
φ 0.1 1f=0.1
φ 0.5 1f=0.5
ψwD, ψwD*tD/CD
1E+0
1E-1
1E-2 pseudo-pressure, ψwD pseudo-pressure derivative, ψwD*tD/CD 1E-3 1E-3
1E-2
1E-1
1E+0
1E+1
1E+2
1E+3
1E+4
1E+5
1E+6
1E+7
tD/CD
Fig. 12. Effect of natural fracture porosity (£1f) on pseudo-pressure and its derivative curves.
4. Conclusions
100
q (104 m3/d)
0.05 φ 1f=0.05
0.1 φ 1f=0.1
0.5 φ1f=0.5
10
1 1E-1
and the transfer from the matrix to the natural fractures. As discussed in the preceding section, the interporosity coefficient (k1mf) has a limited effect on the well pressure and the production performance. So, we assume the k1mf is constant. According to the definition of storability ratio of natural fracture system in inner region (x1f) listed in Table 1, the corresponding values of x1f are 0.333, 0.5 and 0.833. The bigger the £1f is, the lower pressure drop will be in the intermediate flow period, which is welcome by the petroleum engineers. As shown in Fig. 13, increasing the natural fracture porosity significantly enhances the well production rate during the whole production period. So, combing with the previous analysis, it is a desirable to deduce a much more complex fracture network in the development of unconventional gas reservoir, which can increase both the fracture permeability and its porosity.
1E+0
1E+1
1E+2 t (day)
1E+3
1E+4
1E+5
Fig. 13. Effect of natural fracture porosity (£1f) on well production performance.
systems in SRV region, the k1mf mainly affects the interporosity flow period with limited time due to the finite SRV. At later time, all responses merge and the production rate for them are much larger than the conventional model (Wan and Aziz, 2002; Rbeawi and Tiab, 2013). Figs. 8 and 9 show the effects of outer reservoir permeability (k2) on well pressure and production performance in the range 104 mD 6 k2 6 5 103 mD. Because all the dimensionless variables defined in this paper are based on the outer reservoir characteristics, the pseudo-pressure change and its derivative values are becoming bigger with the increasing of k2. Although a higher k2 can increase the flow capacity of fluid flow from the formation to the SRV region in theory, the contribution of the outer reservoir will become evident after about 30 years for these case. If the formation permeability is much lower, it will take even longer. In Fig. 9, it can be found that a bigger k2 will lead to a more significant increase in productivity when the effect comes. Figs. 10 and 11 show the effect of the radius of SRV (rm) on the pseudo-pressure type curves as well as on the rate performance. The effect of the rm on type curves in mainly concentrated on the duration of radial flow in SRV region, which will been longer for a bigger rm (as shown in Fig. 10). Comparing the results with different rm in Fig. 11, bigger rm can cause a noticeable increment of production rate during the later flow period. Therefore, it is desirable to obtain a bigger rm to obtain a higher production rate in the development of unconventional gas reservoir. Figs. 12 and 13 show the effects of natural fracture porosity (£1f) on well transient pressure and production performance in the range of 0.05, 0.1 and 0.5. As the £1f mainly affect the capacity of the fluid stored in the natural fracture system, so a bigger £1f presents a higher density of natural fractures. In theory, increasing the £1f will increase the flow capacities both the natural fractures
In this paper, we investigated the transient pressure response and performance characteristics of fractured horizontal wells in a composite gas reservoir to describe the stimulated region volume around the well and discussed the effects of SRV and reservoir properties on them. From the above analysis, the following conclusions can be drawn: (1) The SRV has a significant importance on the well pressure and transient rate performance compared to that of the conventional model, especially the natural fracture permeability and porosity. (2) Although the contribution of SRV radius on well production rate need nearly 3 years in this synthesis case, the effect is significant when the effect comes. So, we should try our best to create a large SRV volume because of the life of unconventional gas well in much longer than 3 years. (3) Unlike conventional gas reservoir, relatively high outer formation permeability has limited effect on the well production performance for unconventional gas reservoir because of the effect is coming too later. (4) Permeability of the natural- or induced fracture system in the SRV region affects not only the early production rate, but the production rate in later flow period is also much bigger than that of the conventional model. So it is desirable to obtain larger induced-fracture densities in the SRV region.
Acknowledgments This work was supported by National Science Fund for Distinguished Young Scholars of China (Grant No. 51125019) and the Natural Science Foundation of China (Grant No. 51374181). The authors would also like to appreciate the reviewers and editors whose critical comments were very helpful in preparing this article. Appendix A. The instantaneous line source function derivation For a homogeneous porous medium, the differential equation that described the flow of a compressible fluid, the diffusivity equation, is obtained by combing: (1) the law of conservation of mass; (2) the equation of state and real-gas law; and (3) Darcy’s law (Badakhshan and Sobbi, 1995; Lee and Robert, 1996). Inner region The flow equations in the fractured system of inner region is
455
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
1 @ k1f p1f M @p1f / ct1f p1f M @p1f ¼ 1f r q1m r @r l Z g RT @r 3:6 Z g RT @t
ðA:1Þ
/1m ct1m p1m M @p1m þ q1m ¼ 0 Z g RT @t 3:6
ðA:2Þ
q1m is the pseudo-steady-state interporosity flow rate of the gas from the matrix to the fractures, which is (Warren and Root, 1963):
k1m
l
1
wðtÞestD dt D
ðA:12Þ
0
According to the Warren and Root’s approach (1963), it is assumed that a pseudo-steady state flow exists in the matrix system, so the seepage flow in the matrix can be negligible, yields
q1m ¼ a1
wðsÞ ¼
Z
ðq1m p1m q1f p1f Þ
ðA:3Þ
So the dimensionless mathematical models in Laplace space can be simplified as follows Inner region For fracture system
1 @ @ Dw1f ¼ f1 ðsÞsDw1f rD rD @r D @r D þsx1f ð1x1f ÞgrD where f1 ðsÞ ¼ k1mf grD ; k þsð1x Þg 1mf
1f
rD
ðA:13Þ 12 grD ¼ W M12
Outer region For fracture system
Outer region With the same method as inner region, the governing equations of the outer region can also be obtained. The flow equation in fractured system of outer regions is
/2m ct2m p2m M @p2m þ q2m ¼ 0 Z g RT @t 3:6
ðA:5Þ
The Interporosity flow equation of the fluid flow from the matrix to the natural fractures in outer region is (Warren and Root, 1963)
ðq2m p2m q2f p2f Þ
p1f jt¼0 ¼ p1m jt¼0 ¼ p2f jt¼0 ¼ p2m jt¼0 ¼ pi
ðA:7Þ _
Assume that at t = 0, a line source with a constant rate q is withdrawn from the center of the inner region instantaneously (as shown in Fig. 2). By forcing the cumulative flux through the side surface of a small cylinder surrounding the line source whose rate equals that of the fluid extracted from the source, the inner boundary condition is
0
t
limþ
e!0
k1f h
l
_ psc TZ @p 1 r 1f dt ¼ q 2p 3:6 p1f T sc @r r¼e
2f
Dw2f ðrD ; t D ÞjrD !1 ¼ 0
ðA:15Þ
@w1f 1_p T k2f limþ k1f h r D ¼ q sc 2 T @r p e!0 D sc ð/2 c t2 Þfþm lgi L r D ¼e Interface boundary conditions
ðA:8Þ
Dw1f ðrD ; t D ÞjrD ¼rmD ¼ Dw2f ðrD ; t D ÞjrD ¼rmD
ðA:17Þ
@ Dw1f 1 @ Dw2f ¼ @r D rD ¼rmD M12 @r D rD ¼rmD
ðA:18Þ
As Eqs. (A.13) and (A.14) have the form of modified Bessel’s equation, the general solution of them can be expressed as follows (Badakhshan and Sobbi, 1995; Lee and Robert, 1996)
Dw1f ¼ A1 I0
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞr D þ B1 K0 sf1 ðsÞr D
ðA:19Þ
Dw2f ¼ A2 I0
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf2 ðsÞr D þ B2 K0 sf2 ðsÞr D
ðA:20Þ
According to the outer boundary condition Eq. (A.15) and the pffiffiffiffiffiffiffiffiffiffiffiffi sf2 ðsÞrD ! 1), the result
The lateral infinite outer boundary condition is used for application in this paper, so
theory of Bessel’s function (limrD !1 I0
p2f ðr; tÞjr!1 ¼ pi
of coefficient A2 in Eq. (A.15) can be obtained as follows
ðA:9Þ
These two boundary conditions impose pressure and rate continuity at the interface between the two zones.
p1f ðr; tÞjr¼rm ¼ p2f ðr; tÞjr¼rm
qjp1f
qjp2f @p2f @p1f ¼ @r r¼rm M12 @r r¼rm
ðA:16Þ
ðA:6Þ
Providing the initial pressure in all over the place is uniform and equal to pi, the initial conditions are
Z
2mf
According to characteristics of the Dirac delta function (Zhao et al., 2013), the inner boundary condition becomes
is
l
þsx2f ð1x2f Þ where f2 ðsÞ ¼ k2mf . k þsð1x Þ
ðA:4Þ
The flow equation in matrix system with the steady-state flow
q2m ¼ a2
ðA:14Þ
Outer boundary condition
1 @ k1f p1f M @p1f / ct2f p2f M @p2f ¼ 2f r q2m r @r l Z g RT @r 3:6 Z g RT @t
k2m
1 @ @ Dw2f ¼ sf2 ðsÞDw2f rD rD @r D @r D
A2 ¼ 0
ðA:21Þ
ðA:10Þ
Substituting Eq. (A.19) into the inner boundary condition Eq. (A.16) yields
ðA:11Þ
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞrD B1 K 1 sf1 ðsÞr D limþ r D A1 I1
Substituting the pseudo-pressure and dimensionless variable defined in Table 1 into Eqs. (A.1)–(A.11), the differential equations and the boundary conditions above are all transformed into dimensionless form, and then taking the following Laplace transform with them
e!0
_ p T 1 k2f pffiffiffiffiffiffiffiffiffiffiffiffi q sc ¼ pk1f h sf1 ðsÞ T sc ð/2 ct2 Þfþm lgi L2
r D ¼e
ðA:22Þ
According to the characteristic of Bessel function (limx!0 xK1 ðxÞ ! 1, limx!0 xI1 ðxÞ ! 0), the coefficient B1 in Eq. (A.22) can be solved directly, which is
456
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 _
B1 ¼ q
psc T
k2f
1
ðA:23Þ
pT sc ð/2 ct2 Þfþm lgi L2 k1f h
Substituting Eqs. (A.19) and (A.20) into Eqs. (A.17) and (A.18) results in the following equations:
A1 I0
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞr mD þ B1 K0 sf1 ðsÞr mD ¼ B2 K0 sf2 ðsÞrmD ðA:24Þ
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞr mD B1 K1 sf1 ðsÞr mD A1 I1 pffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf ðsÞ p2ffiffiffiffiffiffiffiffiffiffiffiffi K1 sf2 ðsÞr mD ¼ B2 M12 sf1 ðsÞ
ðA:25Þ
Combining Eq. (A.24) with Eq. (A.25) and eliminating coefficient B2, the following relationship between A1 and B1 is
A1 ¼ B1
M12 c1 K1 ðc1 r mD ÞK0 ðc2 r mD Þ c2 K0 ðc1 rmD ÞK1 ðc2 rmD Þ M12 c1 I1 ðc1 r mD ÞK0 ðc2 r mD Þ þ c2 I0 ðc1 rmD ÞK1 ðc2 rmD Þ
ðA:26Þ
pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi where c1 ¼ sf1 ðsÞ, c2 ¼ sf2 ðsÞ. Taking Eqs. (A.23) and (A.26) into Eq. (A.19) yields _
1f ¼ q Dw
psc T
k2f
1
pT sc ð/2 ct2 Þfþm lgi L2 pk1f h
½K 0 ðc1 r D Þ þ AC I0 ðc1 r D Þ ðA:27Þ
where AC ¼
M12 c1 K1 ðc1 r mD ÞK0 ðc2 rmD Þ c2 K0 ðc1 r mD ÞK1 ðc2 rmD Þ M12 c1 I1 ðc1 r mD ÞK0 ðc2 rmD Þ þ c2 I0 ðc1 r mD ÞK1 ðc2 r mD Þ
References Al-Kobaisi, M., Ozkan, E., 2004. A hybrid numerical analytical model of finite conductivity vertical fractures intercepted by a horizontal well. SPE Paper 92040 presented at SPE International Petroleum Conference. Puebla, Mexico. Apaydin, O.G., Ozkan, E., Raghavan, R., 2012. Effect of discontinuous microfractures on ultratight matrix permeability of a dual-porosity medium. SPE Reser. Eval. Eng. 15 (4), 473–485. Badakhshan, A., Sobbi, F.A., 1995. Analytical composite models for estimation of real parameters of naturally fractured reservoirs. SPE Paper PETSOC-95-12 is Presented at the Annual Technical Meeting, Calgary, Canada. Brown, M., Ozkan, E., Raghavan, R., et al., 2011. Practical solutions for pressuretransient responses of fractured horizontal wells in unconventional reservoirs. SPE Reser. Eval. Eng. 14 (6), 663–676. Cai, J.C., Yu, B.M., 2011. A discussion of the effect of tortuosity on the capillary imbibition in porous media. Trans. Porous Media. 89 (2), 251–263. Clarkson, C.R., 2013. Production data analysis of unconventional gas wells: Review of theory and best practices. Int. J. Coal Geol. 109–110, 101–146. Crosby, D.G., Rahman, M.M., Rahman, M.K., et al., 2002. Single and multiple transverse fracture initiation from horizontal wells. J. Petrol. Sci. Eng. 35 (3–4), 191–204. de Swaan-O, A., 1976. Analytical solutions for determining naturally fractured reservoir properties by well testing. SPE J. 16 (3), 117–122.
Escobar, F.H., Martinez, J.A., Matilde, M.M., 2013. Pressure transient analysis for a reservoir with a finite-conductivity fault. CT&F – Ciencia, Tecnología y Futuro 5 (2), 5–18. Everdingen, V., Hust, A.F., 1949. The application of the Laplace transformation to flow problems in reservoirs. J. Petrol. Technol. 1 (12), 305–324. Guo, G.L., Evans, R.D., Chang, M.M., 1994. Pressure-transient behavior for a horizontal well intersecting multiple random discrete fractures. SPE Paper 28390 Presented at SPE Annual Technical Conference and Exhibition, New Orleans, USA. Horne, R.N., Temeng, K.O., 1995. Relative productivities and pressure transient modeling of horizontal wells with multiple fractures. SPE Paper 29891 Presented at SPE Middle East Oil Show, Bahrain, Bahrain. Kazemi, H., 1969. Pressure transient analysis of naturally fractured reservoirs with uniform fracture distributions. SPE J. 9 (4), 451–462. Lee, S.T., Brockenbrough, J.R., 1986. A new approximate analytic solution for finite conductivity vertical fractures. SPE Form. Eval. 1 (1), 75–88. Lee, W.J., Robert, A.W., 1996. Gas reservoir engineering. Society of Petroleum Engineers, SPE Textbook Series, vol. 5, United State, pp. 80–82. ISBN: 978-155563-073-7. Milton, A., Irene, A.S., 1964. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Courier Dover Publications, United States, pp. 355–435. Ozkan, E., Brown, M., Raghavan, R., et al. 2009. Comparison of fractured horizontalwell performance in conventional and unconventional reservoirs. SPE Paper 121290-MS Presented at SPE Western Regional Meeting, San Jose, California, USA. Ozkan, E., Brown, M., Raghavan, R., et al., 2011. Comparison of fractured-horizontalwell performance in tight sand and shale reservoirs. SPE Reser. Eval. Eng. 14 (2), 248–259. Raghavan, R., Chen, Chin.-Cheng., Agarwal, Bijin., 1997. An analysis of horizontal wells intercepted by multiple fractures. SPE J. 2 (3), 235–245. Rbeawi, S.A., Tiab, D., 2013. Pressure behaviours and flow regimes of a horizontal well with multiple inclined hydraulic fractures. Int. J. Oil, Gas Coal Technol. 6 (1/ 2), 207–241. Stalgorova, E., Mattar, L., 2012a. Practical Analytical model to simulate production of horizontal wells with branch fractures. SPE Paper 162515 Presented at SPE Canadian Unconventional Resource Conference, Calgary, Alberta, Canada. Stalgorova, E., Mattar, L., 2012b. Analytical model for history matching and forecasting production in multifrac composite systems. SPE Paper 162516-MS Presented at SPE Canadian Unconventional Resource Conference, Calgary, Alberta, Canada. Stehfest, H., 1970. Numerical inversion of Laplace transforms-algorithm 368. Commun. ACM. 13 (1), 47–49. Wan, J., Aziz, K., 1999. Multiple hydraulic fractures in horizontal wells. SPE Paper 54627-MS presented at SPE Western Regional Meeting, Anchorage, Alaska, USA, 26–27. Wan, J., Aziz, K., 2002. Semi-analytical well model of horizontal wells with multiple hydraulic fractures. SPE J. 7 (4), 437–445. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3 (3), 245–255. Zeng, F.H., Zhao, G., 2010. The optimal hydraulic fracture geometry under nonDarcy flow effects. J. Petrol. Sci. Eng. 72, 143–157. Zerzar, A., Bettam, Y., Tiab, D., 2003. Interpretation of multiple hydraulically fractured horizontal wells in closed systems. SPE Paper 84888 Presented at SPE International Improved Oil Recovery Conference in Asia Pacific, Kuala Lumpur, Malaysia. Zhao, Y.L., Zhang, L.H., Wu, F., 2012. Pressure transient analysis for multi-fractured horizontal well in shale gas reservoirs. J. Petrol. Sci. Eng. 90–91, 31–38. Zhao, Y.L., Zhang, L.H., Zhao, J.Z., et al., 2013. ‘‘Triple porosity’’ modeling of transient well test and rate decline analysis for multi-fractured horizontal well in shale gas reservoirs. J. Petrol. Sci. Eng. 110, 253–262.
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir Yu-Long Zhao ⇑, Lie-Hui Zhang, Jian-Xin Luo, Bo-Ning Zhang State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, PR China
a r t i c l e
i n f o
Article history: Received 30 October 2013 Received in revised form 25 February 2014 Accepted 10 March 2014 Available online 19 March 2014 This manuscript was handled by Peter K. Kitanidis, Editor-in-Chief, with the assistance of Markus Tuller, Associate Editor Keywords: Fractured horizontal well Unconventional gas reservoir Stimulated reservoir volume Well testing model Performance analysis
s u m m a r y This paper extended the conventional multiple hydraulic fractured horizontal (MFH) well into a composite model to describe the stimulated reservoir volume (SRV) caused by hydraulic fracturing. Employing the Laplace transform, Source function, and Dirac delta function methods, the continuous linear source function for general composite dual-porosity is derived, and the solution of the MFH well in a composite gas reservoir is obtained with the numerical discrete method. Through the Stehfest numerical algorithm and Gauss elimination method, the transient pressure responses for well producing at a constant production rate and the production rate vs. time for constant bottomhole pressure are analyzed. The effects of related parameters such as natural permeability and radial of the SRV region, formation permeability and interporosity coefficient on transient pressure and production performance are analyzed as well. The presented model and obtained results in this paper not only enrich the well testing models of such unconventional reservoir, but also can use to interpret on-site data which have significance on efficient reservoir development. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction Horizontal well with massive hydraulic fractures has been proved to be a common way to produce hydrocarbons from ultra-low permeability reservoirs, such as shale gas reservoir, and the objective of it has not only been to create several highconductivity flow paths, but also activate and connect existing discrete or micro-seismic natural fractures so as to generate large fracture networks (Clarkson, 2013). The volume of a reservoir containing the fracture network and the main hydraulic fractures with which effectively stimulated to increase the well performance is defined as SRV (Ozkan et al., 2009, 2011; Stalgorova and Mattar, 2012b; Clarkson, 2013). As transient pressure analysis been an important method for petro-engineers to study the formation characteristics, many models have been proposed to describe the theory of fluid seepage flow from the formation to the well and many detailed analytical solutions are developed. Guo et al. (1994) developed methods for the performance prediction of multiple fractured horizontal wells without considering the effect of interface between fractures. Horne and Temeng (1995) accounted for the interference between multiple fractures by the superposition of influence function. ⇑ Corresponding author. Tel.: +86 13568841671. E-mail address: [email protected] (Y.-L. Zhao). http://dx.doi.org/10.1016/j.jhydrol.2014.03.026 0022-1694/Ó 2014 Elsevier B.V. All rights reserved.
Raghavan et al. (1997) developed a mathematical model to discern the characteristic response of multiple-fractured horizontal wells. Wan and Aziz (1999) developed general solution for horizontal wells with multiple fractures. Crosby et al. (2002) and Wan and Aziz (2002) derived semi-analytical solutions to respectively describe the pressure response for transverse fractures and fractures rotated at any angle to the horizontal well. Zerzar et al. (2003) combined the boundary element method and Laplace transformation to present a comprehensive solution for multiple vertical fractures horizontal wells. Al-Kobaisi and Ozkan (2004) presented a hybrid numerical-analytical model for the pressure transient response of horizontal wells intercepted by a vertical fracture. Rbeawi and Tiab (2013) introduced new analytical models to investigate the pressure behavior and flow regimes of a horizontal well with multiple hydraulic fractures regarding them as being either longitudinal or transverse, either vertical or inclined, and either symmetrical or asymmetrical. Zhao et al. (2012, 2013) used source function to analyze the transient pressure response of a MFH well in shale gas reservoir. Although the above researchers have made much works on transient pressure analysis of MFH well in both conventional and unconventional reservoir, most of the models they proposed are not take in account the effect of the induced fracture around the main hydraulic fractures (called the fractured network). If we want to describe the fracture network in well testing models similar
448
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
Nomenclature cg ct1m,ct1f ct2m,ct2f
ct2
C h k2 k1m, k1f k2m, k2f k2 DLfi L LH M M12 p1m, p1f p2m, p2f psc _ q qscL qfi qsc rm r R s
gas compressibility, MPa1 total compressibility in the matrix and natural fracture system of the inner region respectively, MPa1 total compressibility in the matrix and natural fracture system of the outer region respectively (for composite dual porosity model), MPa1 total compressibility of the outer reservoir(for inner natural fracture and outer homogeneous composite reservoir model), MPa1 wellbore storage coefficient, m3/MPa formation thickness, m reservoir permeability, D matrix and fracture network permeability in SRV region, D matrix and micro-fracture permeability in outer region, D formation permeability of the outer reservoir, D the length of the ith discrete fracture, m half horizontal well length, m horizontal well length, m molar mass of nature gas, mol/kg mobility ratio, dimensionless, defined in Table 1 matrix/fracture network system pressure in SRV, MPa matrix/micro-fracture system pressure in outer region, MPa pressure at standard condition, psc = 0.10325 MPa instantaneous line source strength, sm3 continuous line source flow rate, m3/h the production rate of the ith discrete fracture, m3/d well production rate, m3/d the radial distance of the SRV, m radius, m the universal gas constant, R = 0.008314 (MPa m3)/ (K kmol) Laplace variable
with the real geological situation, we should implement lots of CT scanning and core experience and then establish the mathematical models of the induced fractures one by one without any simplification. Although this idea is very good, the implementations of such huge works are impossible. In order to concisely describe the fracture network (natural or induced) around the hydraulic fractures, assuming uniform distribution of identical hydraulic fractures along the length of the horizontal well, Ozkan et al. (2009, 2011) utilizes the concept of Lee and Brockenbrough (1986) of tri-linear model with inner reservoir of naturally fractured to represent the MFH well performance in unconventional reservoirs (see diagram in Fig. 1A). Brown et al. (2011) presented an analytical tri-linear flow model that incorporates transient fluid transfer from matrix to fracture to simulate the pressure transient and production behavior of fractured horizontal wells in unconventional reservoirs. Using the tri-linear flow model, Apaydin et al. (2012) presented an analytical model with composite matrix blocks to describe the effect of matrix micro-fractures on the effective matrix permeability. With the assumption that the enhanced region occupies only part of the space between the fractures and neglects the region beyond the fractures, Stalgorova and Mattar (2012a) improved the tri-linear flow model to simulate fractures that is surrounded by a simulated region of limited extendibility. Later,
Skin t T Tsc v W12
skin factor, dimensionless production time, h temperature, K temperature at standard condition, Tsc = 293 K flow velocity, m/h storability ratio between inner and outer region, dimensionless, defined in Table zg gas deviation factor, sm3/m3 w pseudo-pressure, MPa2/cp / formation porosity, fraction £1m, £1f matrix and fracture network porosity in SRV, fraction £2m, £2f matrix and micro-fracture porosity, fraction q fluid density in reservoir condition, g/cm3 q1m fluid density in the matrix system, g/cm3 q1 fluid density in the fracture system, g/cm3 l gas viscosity, cp lgi gas viscosity at the initial reservoir condition, cp a1 sharp factor of the inner region, 1/m2 a2 sharp factor of the outer region, 1/m2 x1f storability coefficient of natural fracture network in SRV, defined in Table 1 x2f storability coefficient of natural fractures in outer region, defined in Table 1 k1mf inter-porosity coefficient from matrix into natural fracture network in inner region, defined in Table 1 k2mf inter from matrix into natural fractures in outer region, defined in Table 1 c1, c2, f1(s), f2(s) variables groups Subscript D
dimensionless
Superscript – Laplace transform
Stalgorova and Mattar (2012b) extended the tri-linear flow model to a five regions flow model (see diagram in Fig. 1B), which can not only simulate the partial simulated reservoir volume but also take into account the fluid supplement beyond the fractures. Although the merits of the ‘‘tri-linear’’ and ‘‘five-flow’’ models are obvious, there are still some demerits of them. Both assume that all fractures have the same length and conductivity and are spaced uniformly along the horizontals, and the early-radial flow, later pseudo-radial flow period and interference phenomenon between the fractures flow cannot be observed, which makes the simplified linear flow models too simple. Therefore, we must take different geological conditions into consideration and propose different models to characterize the transient fluid flow into the MFH well with SRV. In this paper, we extend the conventional multi-fractured horizontal well (see Fig. 2b) into a novel composite model (see Fig. 2c), and derive a semi-analytical model to represent flow in an unconventional, ultra-tight formation, as described previously. The model is a dual-porosity composite system, consisting of a macro-fracture network (inner region, which compose of the natural and induced hydraulic fractures and represent the SRV) and a homogeneous matrix system (outer region). The MFH well locates in the center of the inner region with an arbitrary length and spacing (see Fig. 2a).
449
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
Region 3
Region 4
Hydraulic
Outer reservoir
Region 5
fracture
Inner reservoir naturally
Region 2
Region 1
Horizontal well
B—five-region model
A—tri-linear flow model
Fig. 1. Schematics of tri-linear and five-region flow model for MFH well.
Micro fractures
Micro fractures of
of formations
formations
Stimulation reservoir
Main hydraulically
volume (SRV)
fractures
A—(Complex geometry patterns)
SRV
Reservoir formation
B--(Conventional MFH model)
C—(Composite model in this paper)
Fig. 2. Hydraulic multiple fractured horizontal well with SRV and its presentation models.
2. Continuous line sources in composite dual-porosity gas reservoir To simulate the performance of stimulated reservoir volume in unconventional gas reservoir, we allow the SRV to be naturally fractured, and use the dual porosity idealization to represent the naturally fractured medium. Although the fluid transfer from the matrix to the fractures can be idealized by the method proposed
by Kazemi (1969) and de Swaan-O (1976), our choice is the pseudo-steady fluid transfer equation of Warren and Root (1963). Fig. 3 is a schematic of a continuous source located in the center of a composite dual-porosity gas reservoir. Although in some tight gas reservoirs, only the stimulated reservoir volume can be treated as a dual-porosity idealization, in order to derive the general source solution of this kind of reservoir, we treat the outer zone as the dual-porosity medium as well. This model can be applied
450
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
Fig. 3. Schematic of a continuous line source in a radial composite reservoir.
for these reservoirs with developed micro-fractures in shale and coaled methane gas reservoirs. If we further want to obtain the continuous source function with the homogeneous outer reservoir, the results can be easily derived from the former formula after simple transformation. The following assumptions are made during the derivation of the linear source function: (1) the reservoir is composite with the inner and outer region of dual-porosity, and the reservoir has a constant and uniform thickness with the upper and lower boundaries closed; (2) porosity and permeability is homogeneous in each formation medium system and region; and (3) fluid flow in the media follows the Darcy flow theory and neglects the gravitational and frictional effects (Zeng and Zhao, 2010).
2.2. Continuous line source function Utilizing the above governing equation for both inner and outer regions and combing the corresponding boundary conditions and initial conditions, the fully mathematical models of an instantaneous line source withdraw from the center of inner region can be established (see Appendix). The instantaneous line source function of composite reservoir with dual-porosity inner and outer reservoirs is obtained and the detailed derivations are presented in Appendix, and the dimensionless variables are defined in column 2 of Table 1. The instantaneous line source function can be obtained as _
Dw1f ¼ 2.1. Governing equations for this model
ð1Þ
1 @/ / @p
ð3Þ
1 @ q z @ p ¼ q @p p @p z
ð4Þ
Ignoring gravity effects, the Darcy’s Law in radial coordinate is given by 3
v ¼ 10
k @p l @r
Using the Duhamel theory (Everdingen and Hust, 1949), the continuous line source function can be obtained by
Dw1f ¼
The isothermal gas compressibility is
cg ¼
k1mf þ sx1f ð1 x1f ÞgrD grD ; k1mf þ sð1 x1f ÞgrD k2mf þ sx2f ð1 x2f Þ W12 ; grD ¼ f 2 ðsÞ ¼ k2mf þ sð1 x2f Þ M12
ð2Þ
The porosity compressibility is
cf ¼
M12 c1 K1 ðc1 r mD ÞK0 ðc2 r mD Þ c2 K0 ðc1 r mD ÞK1 ðc2 r mD Þ ; M12 c1 I1 ðc1 r mD ÞK0 ðc2 r mD Þ þ c2 I0 ðc1 r mD ÞK1 ðc2 r mD Þ qffiffiffiffiffiffiffiffiffiffiffi ffi qffiffiffiffiffiffiffiffiffiffiffiffi c1 ¼ sf1 ðsÞ; c2 ¼ sf2 ðsÞ AC ¼
f1 ðsÞ ¼
The real-gas law is simply the pressure/volume relationship predicted by a correction factor that accounts for the nonideal behavior of the gas, and the real-gas law is as follows
pM ¼ Z g qRT
ð6Þ
where
The governing equations of this composite reservoir can be derived by combining the continuity equation, the flux (Darcy) Law, and an equation of state (Badakhshan and Sobbi, 1995). The continuity equation in a homogeneous, cylindrical reservoir can be obtained by the mass conservation, which is
1 @ðr qv Þ 1 @ðq/Þ þ ¼0 r @r 3600 @t
psc T q k2f ½K ðc r Þ þ AC I0 ðc1 r D Þ pT sc lgi hð/2 ct2 Þfþm k1f L2 0 1 D
ð5Þ
For gas reservoir, both the gas viscosity l and compressibility cg are the function of reservoir pressure p, in order to make the problems more tractable, we make the assumptions that both l and cg are evaluated at the initial reservoir pressure pi.
psc T 1 k2f pT sc lgi hð/2 ct2 Þfþm k1f L2
Z
t
qscL ðsÞSc ðtD sÞds
ð7Þ
0
where Sc ¼ L1 ½Sc , Sc ¼ ½K0 ðc1 rD Þ þ AC I0 ðc1 rD Þ Taking the Laplace transform with tD in Eq. (7), the continuous line source function in inner region of composite dual-porosity reservoir with upper and lower boundaries closed is
~scL Dw1f ¼ q
psc T 1 ½K0 ðc1 r D Þ þ AC I0 ðc1 rD Þ T sc 86:4pk1f h
ð8Þ
Assuming that the line source with a constant rate qscL, we can define the dimensionless pseudo-pressure as follows
wD ¼
k2f hT sc 3:684 103 qscL psc T
Dw
ð9Þ
Combining Eq. (8) with Eq. (9), the continuous line source function is
w1fD ¼
1 ½K0 ðc1 r D Þ þ AC I0 ðc1 r D Þ M12 s
ð10Þ
451
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 Table 1 Definition of the dimensionless variables. Dimensionless variables Pseudo-pressure Dimensionless time
Both inner and outer region are dual-porosity Rp wl ¼ 2 0 l lpZg dp tD ¼ ð/
3:6k2f t
2 c t2 Þfþm
Dimensionless radial distance Storability ratio of fracture system in region 1
rD ¼
Inner region with dual-porosity and outer region is homogeneous Rp wl ¼ 2 0 l lpZ g dp 2t tD ¼ / 3:6k c l L2
lgi L2
2 t2
r L
gi
r D ¼ Lr
Þf x1f ¼ ð/ð/1 c1t1ct1Þfþm
Þf x1f ¼ ð/ð/1 c1t1ct1Þfþm
Storability ratio of fracture system in region 2
x2f ¼
Interporosity coefficient of region 1
k1mf ¼ a1 kk1m L2
k1mf ¼ a1 kk1m L2
Interporosity coefficient of region 2
k2mf ¼ 2 kk2m L2 2f M12 ¼ kk1f ==ll 2f ð/ c Þ W12 ¼ ð/12 ct1t2 Þfþm fþm
M12 ¼ kk1f2 ==ll
ð/2 ct2 Þf ð/2 ct2 Þfþm 1f
Mobility ratio between region 1 and region 2 Storability ratio between regions 1 and region 2 Dimensionless wellbore storage
1f
a
CD ¼
W12 ¼
C 2pð/2 ct2 Þfþm hL2
CD ¼
ð/1 ct1 Þfþm /2 ct2
C 2p/2 ct2 hL2
Note: when the inner region (SRV region) of reservoir is dual porosity and outer region is homogeneous, the continuous line source function is obtained by substituting the k2mf and x2f equal to 0 and 1 respectively. The naturally fracture’s permeability k2f will be replaced by k2.
When the outer reservoir does not have micro-fractures, the model may transform into an inner dual-porosity and outer homogeneous matrix system. With the same solving process described in Appendix and the new dimensionless variables defined in column 3 of Table 1, the source function can be obtained as
w1fD ¼
1 ½K0 ðc1 r D Þ þ AC I0 ðc1 r D Þ M12 s
ð11Þ
Where
AC ¼
pffiffi pffiffi pffiffi M12 c1 K1 ðc1 rmD ÞK0 ð sr mD Þ sK0 ðc1 r mD ÞK1 ð sr mD Þ pffiffi pffiffi M12 c1 I1 ðc1 r mD ÞK0 ð sr mD Þ þ sI0 ðc1 rmD ÞK1 ðr mD Þ
3.1. Analytical solutions The mathematical model assumes the wellbore to be intersected by N fractures, each fracture is assumed to have distinct property and can be at any position along the wellbore, and all the fractures are transverse to the well and fully-penetrating to the formation. We make the assumption that the flow from the reservoir to the wellbore sections is negligible comparing to the flow to the hydraulic fracture plane. The fractures are assumed to produce at a constant wellbore pressure (the horizontal well can always be treated as infinite in tight gas reservoir for low production rate). Giving that each half fracture is dispersed into MM sections of length DLfi, the total flow rate will be the summation of the fluid flow from each discrete fracture, which is:
qDi ¼ 1
wwDi ¼
MM2N X j¼1
qDfj wwDi;j s
Comparing to the mainly hydraulic fractures and wellbore, the formation permeability is too small, and the well production rate is usually very small for the unconventional gas reservoir. Therefore, we can treat them as infinite-conductivity, which will indicate zero pressure drops in them. So
ð12Þ
where qDi is the dimensionless discrete fracture rate, which is defined as
qfi DLfi qsc
2
wwD1;1
6 6 6 6 wwDk;1 6 6 6 6 4 wwD2N M;1 1
wwD1;k
wwD1;2N M
wwDk;k
wwDk;2N M
wwD2N M;k 1
wwD2N M;2N M
1
1
32
qD1
MM2N X Z tD j¼1
0
qDfj ðsÞw0wDi;j ðt D sÞds
ð14Þ
3
2
0
3
7 6 7 6 1 7 76 qD2 7 6 0 7 76 7 6 7 7 7 67 6 1 76 7¼6 7 7 7 67 6 1 76 7 6 7 76 7 6 7 1 54 qDM 2N 5 4 0 5 0
wwD
1 ð17Þ
In the above matrix, the symbol wwDi;j (both i and j vary in the range 1, 2N * MM) represents the effect of the production from element j at the location of element i. method to calculate the matrix elements for different configurations with the combination of continuous line source function can be referred to the paper of Zhao et al.(2012, 2013), which is 2
wwDi;j ¼
1 M12 sDLfD
Z j
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 2 2 xDi xD j a þ yDi ywD j þ7 6 K0 c1 6 7 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 7da 6 j 4 2 2 5 AC I 0 c 1 xDi xD j a þ yDi ywD j
DLfD j 2 DLfD
2
ð18Þ
ð13Þ
The dimensionless pressure drop, pDi at element i in terms of pressure drops, pDi,j, at element i due to the production from element j = 1, . . . , MM * 2N, can be obtained by the following convolution:
wwDi ðtD Þ ¼
ð16Þ
There are 2N * MM + 1 equations when applying Eq. (15) at each element together with Eq. (12), which can solve the 2N * MM + 1 unknowns of wwD , qD1, qD2, . . . , qDM⁄2N1, qDM*2N. The matrix expression is
i¼1
qDi ¼
ð15Þ
wwD ¼ wwDi¼1;...;MM2N
3. Pressure behaviors for MFH well in composite gas reservoir
MM2N X
Assuming that the production rate in each discrete fracture is constant, in terms of Laplace transform, Eq. (14) can be written as
3.2. Computational consideration of the matrix coefficient According to the matrix described in Eq. (18), if the matrix coefficients are obtained, the dimensionless bottomhole pressure as well as the dimensionless production rate distribution in real time space along the fractures will be solved by the Gauss elimination method and Stehfest numerical inversion algorithm (Stehfest,
452
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
1970). To compute the coefficients efficiently, they can be calculated with different methods according to the observation point location i and the source element location j. If the observation and source element point locate at the same fracture, which means they have the same Y-axis, the calculation of Eq. (18) can be simplified as:
wwDi;j ¼
1 M12 sDLfDi
Z
DLfDj
2
DLfDj
ð19Þ
Using H0(x) to denote I0(x) and K0(x), the definite integration described in Eq. (19) can be simplified as follows, 8Rb Rb 1 > H0 ðaÞda 0 2 H0 ðaÞda if i > j > Z DLfDj 0 < 2 R b1 1 H ð c jx x a jÞd a ¼ 2 H ð a Þd a if i ¼ j 0 Di Dj 0 1 DLfDj 0 c1 > 2 > R b1 : R b2 H ð a Þd a H ð a Þd a if i < j 0 0 0 0 ð20Þ DLfDj where b1 ¼ c1 xDi xDj þ , b2 ¼ c1 xDi xDj 2 . Rx The definite integration of 0 H0 ðaÞda can be calculated with the numerical algorithm described by Milton and Irene (1964). When the observation point and the source element locate at different fractures, the result of Eq. (18) can be calculated directly by the numerical integration algorithm, such as adaptive algorithms, Gauss–Legendre and variable step size Runge–Kutta methods (Milton and Irene,1964; Escobar et al., 2013; Cai and Yu, 2011).
DLfDj 2
3.3. Transient pressure responses when well producing at a constant rate with considering wellbore storage and skin effects When we want to take into account the effects of wellbore storage as well as skin factor, using the Duhamel’s theorem (Everdingen and Hust, 1949), the dimensionless wellbore pressure in Laplace space becomes
swwDN þ Skin s þ C D s2 ðswwDN þ Skin Þ
ð21Þ
Combing Eq. (21) with Stehfest Numerical algorithm (Stehfest, 1970), the dimensionless pseudo-pressure and its derivative vs. time on log–log curves can be obtained. 3.4. Solution of well producing at a constant wellbore pressure When well producing at a constant bottomhole pressure, the dimensionless rate can be defined as follows:
qD ¼
qD ¼
1 s2 wwD
ð23Þ
The above equation is the solution for a multi-fractured horizontal well producing at a constant wellbore pressure.
½K0 ðc1 jxDi xDj ajÞ
2
þ AC I0 ðc1 jxDi xDj ajÞda
wwD ¼
Combing Eq. (22) with the expression of dimensionless bottom hole pressure in Table 1, the following equation can be obtained
3:684 103 psc qsc T kf hT sc Dw
ð22Þ
3.5. Discussion of the result We use the synthetic data in Table 2 to analyze the pressure response and the rate decline of the mathematical model derived in this paper, the results are shown in Figs. 4–13. Figs. 4–7 compare the differences in the transient pressure responses and the performances of a fractured horizontal well in homogeneous and fractured composite reservoirs with different parameters. The synthesis data are listed in Table 2. The solid and dash lines in these figures are obtained by the composite model with inner naturally fractured (SRV) and the outer homogeneous, and the dot lines with circle marker are obtained by the conventional model of MFH well in homogeneous reservoir (Raghavan et al., 1997; Wan and Aziz, 2002; Rbeawi and Tiab, 2013). Figs. 4 and 6 are the transient pressure responses of the well at a constant production rate. The effect of the fracture network permeability (k1f) in SRV region on type curves is shown in Fig. 4. It can be natural seen from the lines that the bigger k1f is, the smaller value of wwD and wwD * tD/CD will be in the early linear and early radial flow period in SRV region. Compared to the conventional MFH well model in homogeneous reservoir (Wan and Aziz, 2002; Rbeawi and Tiab, 2013), due to the existence of the SRV region, the pressure drop is smaller than that of the model without SRV when the well at a constant production rate. Fig. 6 shows the effect of interporosity coefficient (k1mf), which presents the capacity of the fluid transfer from the matrix to the natural fractures under same pressure difference. It can be seen that the bigger k1mf is, the early the transition flow between the matrix and natural fractures system will be, and also the corresponding pressure drop is much lower. Figs. 5 and 7 are the rate decline performance of the MFH well in composite reservoir under a constant bottomhole pressure. The comparison of the results of the homogeneous reservoir and the composite system with naturally fractured SRV displays significant contribution of the SRV to the production performance of the MFH. This contribution helps the well sustain a higher production rate during the all production flow period. In practical, the bigger the permeability of the fracture system is, the bigger production rate of the well is (as shown in Fig. 5). Fig. 7 shows the effect of interporosity coefficient (k1mf) between matrix and natural fracture
Table 2 Synthetic data used for the discussion of the results. Formation thickness, h, m Wellbore radius, rw, m Hydraulic fracture half-length, xF, m The radius of stimulated region, rm Bottomhole pressure, pwf, MPa Gas viscosity, l, cp Skin factor, skin
50 0.1 75 500 4 0.021 0.01
Inner region Matrix permeability, k1m, md Matrix porosity, £1m Total compressibility in matrix system, Ct1m, MPa1 Fracture permeability, k1f, md Fracture porosity, £1f Total compressibility in natural fracture system, Ct1f, MPa1 Interporosity coefficient, k1mf
Initial reservoir pressure, pi, MPa Horizontal well length, LH, m Number of fractures, N Reservoir temperature, T, K Molecular weight, Mg, kg/mol Wellbore storage coefficient, C, m3/MPa
35 1000 4 353 16 100
Outer region 0.0001 0.1 0.01 0.01 0.1 0.02 0.2
Permeability, k2, md Porosity, £2 Total compressibility, Ct2, MPa1
0.001 0.1 0.02
453
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 1E+1
1E+1
kkf=0.005mD f1 =0.005mD kkf=0.05mD f1 =0.05mD
Conventional model
1E-1
kk2=0.001mD 2=0.001mD
kk2=0.005mD 2=0.005mD
1E+0 ψwD, ψwD*tD/CD
ψwD, ψwD*tD/CD
1E+0
kk2=0.0001mD 2=0.0001mD
kkf=0.01mD f1 =0.01mD
1E-1
1E-2
1E-2 1E-3
pseudo-pressure, ψwD
pseudo-pressure, ψwD pseudo-pressure derivative, ψwD*tD/CD
pseudo-pressure derivative, ψwD*tD/CD 1E-3 1E-3
1E-2
1E-1
1E+0
1E+1 tD/CD
1E+2
1E+3
1E+4
1E-4 1E-3
1E+5
Fig. 4. Effect of induced fracture network permeability in SRV on pseudo-pressure and its derivative curves.
1E-2
1E-1
1E+0
1E+1 tD/CD
1E+2
1E+3
1E+4
1E+5
Fig. 8. Effect of outer reservoir permeability (k2) on pseudo-pressure and its derivative curves.
1000
100
kk2=0.001mD 2=0.001mD
kk2=0.0001mD 2=0.0001mD 100 k2=0.005mD k2=0.005mD q (104 m3/d)
q (104 m3/d)
10
10
1
1
0.1 1E-1
kkf=0.005mD f1 =0.005mD
kkf=0.01mD f1 =0.01mD
kkf=0.05mD f1 =0.05mD
Conventional model
1E+0
1E+1
0.1
1E+2
1E+3
0.01 1E-1
1E+4
1E+0
1E+1
1E+2
t (day)
λλ1m-f=1.00 1m-f =1.00
λλ1m-f=0.10 1m-f =0.10
Conventional model
1E-2
rrm=500m m=500m
rrm=750m m=750m
rrm=1000m m=1000m
rrm=1250m m=1250m
1E+8
1E-1
1E-2
pseudo-pressure, ψwD
pseudo-pressure, ψwD
pseudo-pressure derivative, ψwD*tD/CD 1E-2
1E-1
1E+0
1E+1
1E+2
1E+3
1E+4
pseudo-pressure derivative, ψwD*tD/CD
1E+5
1E-3 1E-3
1E-2
1E-1
1E+0
tD/CD
1E+1 1E+2 tD/CD
1E+3
1E+4
1E+5
1E+6
Fig. 10. Effect of SRV radius (rm) on pseudo-pressure and its derivative curves.
Fig. 6. Effect of k1mf on pseudo-pressure and its derivative curves.
100
100
rrm=500m m=500m
rrm=750m m=750m
rm=1000m rm=1000m
rrm=1250m m=1250m
10
q (104 m3/d)
q (104 m3/d)
10
0.1 1E-1
1E+7
1E+0 ψwD, ψwD*tD/CD
ψ wD, ψ wD*tD/CD
λλ1m-f=10.0 1m-f =10
1E-1
1
1E+6
1E+1
1E+1
1E-3 1E-3
1E+5
Fig. 9. Effect of outer reservoir permeability (k2) on well production performance.
Fig. 5. Effect of kf on well production performance.
1E+0
1E+3 1E+4 t (day)
1 λλ1m-f=0.10 1m-f =0.10
λλ1m-f=1.00 1m-f =1.00
λλ1m-f=10.0 1m-f =10
Conventional model
1E+0
1E+1
1E+2 t (day)
1E+3
1E+4
Fig. 7. Effect of k1mf on well production performance.
1E+5
0.1 1E-1
1E+1
1E+3
1E+5
1E+7
1E+9
t (day)
Fig. 11. Effect of SRV radius (rm) on well production performance.
454
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 1E+1 φ 0.05 1f=0.05
φ 0.1 1f=0.1
φ 0.5 1f=0.5
ψwD, ψwD*tD/CD
1E+0
1E-1
1E-2 pseudo-pressure, ψwD pseudo-pressure derivative, ψwD*tD/CD 1E-3 1E-3
1E-2
1E-1
1E+0
1E+1
1E+2
1E+3
1E+4
1E+5
1E+6
1E+7
tD/CD
Fig. 12. Effect of natural fracture porosity (£1f) on pseudo-pressure and its derivative curves.
4. Conclusions
100
q (104 m3/d)
0.05 φ 1f=0.05
0.1 φ 1f=0.1
0.5 φ1f=0.5
10
1 1E-1
and the transfer from the matrix to the natural fractures. As discussed in the preceding section, the interporosity coefficient (k1mf) has a limited effect on the well pressure and the production performance. So, we assume the k1mf is constant. According to the definition of storability ratio of natural fracture system in inner region (x1f) listed in Table 1, the corresponding values of x1f are 0.333, 0.5 and 0.833. The bigger the £1f is, the lower pressure drop will be in the intermediate flow period, which is welcome by the petroleum engineers. As shown in Fig. 13, increasing the natural fracture porosity significantly enhances the well production rate during the whole production period. So, combing with the previous analysis, it is a desirable to deduce a much more complex fracture network in the development of unconventional gas reservoir, which can increase both the fracture permeability and its porosity.
1E+0
1E+1
1E+2 t (day)
1E+3
1E+4
1E+5
Fig. 13. Effect of natural fracture porosity (£1f) on well production performance.
systems in SRV region, the k1mf mainly affects the interporosity flow period with limited time due to the finite SRV. At later time, all responses merge and the production rate for them are much larger than the conventional model (Wan and Aziz, 2002; Rbeawi and Tiab, 2013). Figs. 8 and 9 show the effects of outer reservoir permeability (k2) on well pressure and production performance in the range 104 mD 6 k2 6 5 103 mD. Because all the dimensionless variables defined in this paper are based on the outer reservoir characteristics, the pseudo-pressure change and its derivative values are becoming bigger with the increasing of k2. Although a higher k2 can increase the flow capacity of fluid flow from the formation to the SRV region in theory, the contribution of the outer reservoir will become evident after about 30 years for these case. If the formation permeability is much lower, it will take even longer. In Fig. 9, it can be found that a bigger k2 will lead to a more significant increase in productivity when the effect comes. Figs. 10 and 11 show the effect of the radius of SRV (rm) on the pseudo-pressure type curves as well as on the rate performance. The effect of the rm on type curves in mainly concentrated on the duration of radial flow in SRV region, which will been longer for a bigger rm (as shown in Fig. 10). Comparing the results with different rm in Fig. 11, bigger rm can cause a noticeable increment of production rate during the later flow period. Therefore, it is desirable to obtain a bigger rm to obtain a higher production rate in the development of unconventional gas reservoir. Figs. 12 and 13 show the effects of natural fracture porosity (£1f) on well transient pressure and production performance in the range of 0.05, 0.1 and 0.5. As the £1f mainly affect the capacity of the fluid stored in the natural fracture system, so a bigger £1f presents a higher density of natural fractures. In theory, increasing the £1f will increase the flow capacities both the natural fractures
In this paper, we investigated the transient pressure response and performance characteristics of fractured horizontal wells in a composite gas reservoir to describe the stimulated region volume around the well and discussed the effects of SRV and reservoir properties on them. From the above analysis, the following conclusions can be drawn: (1) The SRV has a significant importance on the well pressure and transient rate performance compared to that of the conventional model, especially the natural fracture permeability and porosity. (2) Although the contribution of SRV radius on well production rate need nearly 3 years in this synthesis case, the effect is significant when the effect comes. So, we should try our best to create a large SRV volume because of the life of unconventional gas well in much longer than 3 years. (3) Unlike conventional gas reservoir, relatively high outer formation permeability has limited effect on the well production performance for unconventional gas reservoir because of the effect is coming too later. (4) Permeability of the natural- or induced fracture system in the SRV region affects not only the early production rate, but the production rate in later flow period is also much bigger than that of the conventional model. So it is desirable to obtain larger induced-fracture densities in the SRV region.
Acknowledgments This work was supported by National Science Fund for Distinguished Young Scholars of China (Grant No. 51125019) and the Natural Science Foundation of China (Grant No. 51374181). The authors would also like to appreciate the reviewers and editors whose critical comments were very helpful in preparing this article. Appendix A. The instantaneous line source function derivation For a homogeneous porous medium, the differential equation that described the flow of a compressible fluid, the diffusivity equation, is obtained by combing: (1) the law of conservation of mass; (2) the equation of state and real-gas law; and (3) Darcy’s law (Badakhshan and Sobbi, 1995; Lee and Robert, 1996). Inner region The flow equations in the fractured system of inner region is
455
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456
1 @ k1f p1f M @p1f / ct1f p1f M @p1f ¼ 1f r q1m r @r l Z g RT @r 3:6 Z g RT @t
ðA:1Þ
/1m ct1m p1m M @p1m þ q1m ¼ 0 Z g RT @t 3:6
ðA:2Þ
q1m is the pseudo-steady-state interporosity flow rate of the gas from the matrix to the fractures, which is (Warren and Root, 1963):
k1m
l
1
wðtÞestD dt D
ðA:12Þ
0
According to the Warren and Root’s approach (1963), it is assumed that a pseudo-steady state flow exists in the matrix system, so the seepage flow in the matrix can be negligible, yields
q1m ¼ a1
wðsÞ ¼
Z
ðq1m p1m q1f p1f Þ
ðA:3Þ
So the dimensionless mathematical models in Laplace space can be simplified as follows Inner region For fracture system
1 @ @ Dw1f ¼ f1 ðsÞsDw1f rD rD @r D @r D þsx1f ð1x1f ÞgrD where f1 ðsÞ ¼ k1mf grD ; k þsð1x Þg 1mf
1f
rD
ðA:13Þ 12 grD ¼ W M12
Outer region For fracture system
Outer region With the same method as inner region, the governing equations of the outer region can also be obtained. The flow equation in fractured system of outer regions is
/2m ct2m p2m M @p2m þ q2m ¼ 0 Z g RT @t 3:6
ðA:5Þ
The Interporosity flow equation of the fluid flow from the matrix to the natural fractures in outer region is (Warren and Root, 1963)
ðq2m p2m q2f p2f Þ
p1f jt¼0 ¼ p1m jt¼0 ¼ p2f jt¼0 ¼ p2m jt¼0 ¼ pi
ðA:7Þ _
Assume that at t = 0, a line source with a constant rate q is withdrawn from the center of the inner region instantaneously (as shown in Fig. 2). By forcing the cumulative flux through the side surface of a small cylinder surrounding the line source whose rate equals that of the fluid extracted from the source, the inner boundary condition is
0
t
limþ
e!0
k1f h
l
_ psc TZ @p 1 r 1f dt ¼ q 2p 3:6 p1f T sc @r r¼e
2f
Dw2f ðrD ; t D ÞjrD !1 ¼ 0
ðA:15Þ
@w1f 1_p T k2f limþ k1f h r D ¼ q sc 2 T @r p e!0 D sc ð/2 c t2 Þfþm lgi L r D ¼e Interface boundary conditions
ðA:8Þ
Dw1f ðrD ; t D ÞjrD ¼rmD ¼ Dw2f ðrD ; t D ÞjrD ¼rmD
ðA:17Þ
@ Dw1f 1 @ Dw2f ¼ @r D rD ¼rmD M12 @r D rD ¼rmD
ðA:18Þ
As Eqs. (A.13) and (A.14) have the form of modified Bessel’s equation, the general solution of them can be expressed as follows (Badakhshan and Sobbi, 1995; Lee and Robert, 1996)
Dw1f ¼ A1 I0
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞr D þ B1 K0 sf1 ðsÞr D
ðA:19Þ
Dw2f ¼ A2 I0
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf2 ðsÞr D þ B2 K0 sf2 ðsÞr D
ðA:20Þ
According to the outer boundary condition Eq. (A.15) and the pffiffiffiffiffiffiffiffiffiffiffiffi sf2 ðsÞrD ! 1), the result
The lateral infinite outer boundary condition is used for application in this paper, so
theory of Bessel’s function (limrD !1 I0
p2f ðr; tÞjr!1 ¼ pi
of coefficient A2 in Eq. (A.15) can be obtained as follows
ðA:9Þ
These two boundary conditions impose pressure and rate continuity at the interface between the two zones.
p1f ðr; tÞjr¼rm ¼ p2f ðr; tÞjr¼rm
qjp1f
qjp2f @p2f @p1f ¼ @r r¼rm M12 @r r¼rm
ðA:16Þ
ðA:6Þ
Providing the initial pressure in all over the place is uniform and equal to pi, the initial conditions are
Z
2mf
According to characteristics of the Dirac delta function (Zhao et al., 2013), the inner boundary condition becomes
is
l
þsx2f ð1x2f Þ where f2 ðsÞ ¼ k2mf . k þsð1x Þ
ðA:4Þ
The flow equation in matrix system with the steady-state flow
q2m ¼ a2
ðA:14Þ
Outer boundary condition
1 @ k1f p1f M @p1f / ct2f p2f M @p2f ¼ 2f r q2m r @r l Z g RT @r 3:6 Z g RT @t
k2m
1 @ @ Dw2f ¼ sf2 ðsÞDw2f rD rD @r D @r D
A2 ¼ 0
ðA:21Þ
ðA:10Þ
Substituting Eq. (A.19) into the inner boundary condition Eq. (A.16) yields
ðA:11Þ
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞrD B1 K 1 sf1 ðsÞr D limþ r D A1 I1
Substituting the pseudo-pressure and dimensionless variable defined in Table 1 into Eqs. (A.1)–(A.11), the differential equations and the boundary conditions above are all transformed into dimensionless form, and then taking the following Laplace transform with them
e!0
_ p T 1 k2f pffiffiffiffiffiffiffiffiffiffiffiffi q sc ¼ pk1f h sf1 ðsÞ T sc ð/2 ct2 Þfþm lgi L2
r D ¼e
ðA:22Þ
According to the characteristic of Bessel function (limx!0 xK1 ðxÞ ! 1, limx!0 xI1 ðxÞ ! 0), the coefficient B1 in Eq. (A.22) can be solved directly, which is
456
Y.-L. Zhao et al. / Journal of Hydrology 512 (2014) 447–456 _
B1 ¼ q
psc T
k2f
1
ðA:23Þ
pT sc ð/2 ct2 Þfþm lgi L2 k1f h
Substituting Eqs. (A.19) and (A.20) into Eqs. (A.17) and (A.18) results in the following equations:
A1 I0
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞr mD þ B1 K0 sf1 ðsÞr mD ¼ B2 K0 sf2 ðsÞrmD ðA:24Þ
qffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf1 ðsÞr mD B1 K1 sf1 ðsÞr mD A1 I1 pffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi sf ðsÞ p2ffiffiffiffiffiffiffiffiffiffiffiffi K1 sf2 ðsÞr mD ¼ B2 M12 sf1 ðsÞ
ðA:25Þ
Combining Eq. (A.24) with Eq. (A.25) and eliminating coefficient B2, the following relationship between A1 and B1 is
A1 ¼ B1
M12 c1 K1 ðc1 r mD ÞK0 ðc2 r mD Þ c2 K0 ðc1 rmD ÞK1 ðc2 rmD Þ M12 c1 I1 ðc1 r mD ÞK0 ðc2 r mD Þ þ c2 I0 ðc1 rmD ÞK1 ðc2 rmD Þ
ðA:26Þ
pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi where c1 ¼ sf1 ðsÞ, c2 ¼ sf2 ðsÞ. Taking Eqs. (A.23) and (A.26) into Eq. (A.19) yields _
1f ¼ q Dw
psc T
k2f
1
pT sc ð/2 ct2 Þfþm lgi L2 pk1f h
½K 0 ðc1 r D Þ þ AC I0 ðc1 r D Þ ðA:27Þ
where AC ¼
M12 c1 K1 ðc1 r mD ÞK0 ðc2 rmD Þ c2 K0 ðc1 r mD ÞK1 ðc2 rmD Þ M12 c1 I1 ðc1 r mD ÞK0 ðc2 rmD Þ þ c2 I0 ðc1 r mD ÞK1 ðc2 r mD Þ
References Al-Kobaisi, M., Ozkan, E., 2004. A hybrid numerical analytical model of finite conductivity vertical fractures intercepted by a horizontal well. SPE Paper 92040 presented at SPE International Petroleum Conference. Puebla, Mexico. Apaydin, O.G., Ozkan, E., Raghavan, R., 2012. Effect of discontinuous microfractures on ultratight matrix permeability of a dual-porosity medium. SPE Reser. Eval. Eng. 15 (4), 473–485. Badakhshan, A., Sobbi, F.A., 1995. Analytical composite models for estimation of real parameters of naturally fractured reservoirs. SPE Paper PETSOC-95-12 is Presented at the Annual Technical Meeting, Calgary, Canada. Brown, M., Ozkan, E., Raghavan, R., et al., 2011. Practical solutions for pressuretransient responses of fractured horizontal wells in unconventional reservoirs. SPE Reser. Eval. Eng. 14 (6), 663–676. Cai, J.C., Yu, B.M., 2011. A discussion of the effect of tortuosity on the capillary imbibition in porous media. Trans. Porous Media. 89 (2), 251–263. Clarkson, C.R., 2013. Production data analysis of unconventional gas wells: Review of theory and best practices. Int. J. Coal Geol. 109–110, 101–146. Crosby, D.G., Rahman, M.M., Rahman, M.K., et al., 2002. Single and multiple transverse fracture initiation from horizontal wells. J. Petrol. Sci. Eng. 35 (3–4), 191–204. de Swaan-O, A., 1976. Analytical solutions for determining naturally fractured reservoir properties by well testing. SPE J. 16 (3), 117–122.
Escobar, F.H., Martinez, J.A., Matilde, M.M., 2013. Pressure transient analysis for a reservoir with a finite-conductivity fault. CT&F – Ciencia, Tecnología y Futuro 5 (2), 5–18. Everdingen, V., Hust, A.F., 1949. The application of the Laplace transformation to flow problems in reservoirs. J. Petrol. Technol. 1 (12), 305–324. Guo, G.L., Evans, R.D., Chang, M.M., 1994. Pressure-transient behavior for a horizontal well intersecting multiple random discrete fractures. SPE Paper 28390 Presented at SPE Annual Technical Conference and Exhibition, New Orleans, USA. Horne, R.N., Temeng, K.O., 1995. Relative productivities and pressure transient modeling of horizontal wells with multiple fractures. SPE Paper 29891 Presented at SPE Middle East Oil Show, Bahrain, Bahrain. Kazemi, H., 1969. Pressure transient analysis of naturally fractured reservoirs with uniform fracture distributions. SPE J. 9 (4), 451–462. Lee, S.T., Brockenbrough, J.R., 1986. A new approximate analytic solution for finite conductivity vertical fractures. SPE Form. Eval. 1 (1), 75–88. Lee, W.J., Robert, A.W., 1996. Gas reservoir engineering. Society of Petroleum Engineers, SPE Textbook Series, vol. 5, United State, pp. 80–82. ISBN: 978-155563-073-7. Milton, A., Irene, A.S., 1964. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Courier Dover Publications, United States, pp. 355–435. Ozkan, E., Brown, M., Raghavan, R., et al. 2009. Comparison of fractured horizontalwell performance in conventional and unconventional reservoirs. SPE Paper 121290-MS Presented at SPE Western Regional Meeting, San Jose, California, USA. Ozkan, E., Brown, M., Raghavan, R., et al., 2011. Comparison of fractured-horizontalwell performance in tight sand and shale reservoirs. SPE Reser. Eval. Eng. 14 (2), 248–259. Raghavan, R., Chen, Chin.-Cheng., Agarwal, Bijin., 1997. An analysis of horizontal wells intercepted by multiple fractures. SPE J. 2 (3), 235–245. Rbeawi, S.A., Tiab, D., 2013. Pressure behaviours and flow regimes of a horizontal well with multiple inclined hydraulic fractures. Int. J. Oil, Gas Coal Technol. 6 (1/ 2), 207–241. Stalgorova, E., Mattar, L., 2012a. Practical Analytical model to simulate production of horizontal wells with branch fractures. SPE Paper 162515 Presented at SPE Canadian Unconventional Resource Conference, Calgary, Alberta, Canada. Stalgorova, E., Mattar, L., 2012b. Analytical model for history matching and forecasting production in multifrac composite systems. SPE Paper 162516-MS Presented at SPE Canadian Unconventional Resource Conference, Calgary, Alberta, Canada. Stehfest, H., 1970. Numerical inversion of Laplace transforms-algorithm 368. Commun. ACM. 13 (1), 47–49. Wan, J., Aziz, K., 1999. Multiple hydraulic fractures in horizontal wells. SPE Paper 54627-MS presented at SPE Western Regional Meeting, Anchorage, Alaska, USA, 26–27. Wan, J., Aziz, K., 2002. Semi-analytical well model of horizontal wells with multiple hydraulic fractures. SPE J. 7 (4), 437–445. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3 (3), 245–255. Zeng, F.H., Zhao, G., 2010. The optimal hydraulic fracture geometry under nonDarcy flow effects. J. Petrol. Sci. Eng. 72, 143–157. Zerzar, A., Bettam, Y., Tiab, D., 2003. Interpretation of multiple hydraulically fractured horizontal wells in closed systems. SPE Paper 84888 Presented at SPE International Improved Oil Recovery Conference in Asia Pacific, Kuala Lumpur, Malaysia. Zhao, Y.L., Zhang, L.H., Wu, F., 2012. Pressure transient analysis for multi-fractured horizontal well in shale gas reservoirs. J. Petrol. Sci. Eng. 90–91, 31–38. Zhao, Y.L., Zhang, L.H., Zhao, J.Z., et al., 2013. ‘‘Triple porosity’’ modeling of transient well test and rate decline analysis for multi-fractured horizontal well in shale gas reservoirs. J. Petrol. Sci. Eng. 110, 253–262.