Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms

Journal of Hydrology 510 (2014) 299–312

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms Hai-Tao Wang ⇑ State Key Laboratory of Oil & Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 0086.610500, Sichuan Province, China

a r t i c l e

i n f o

Article history: Received 18 July 2013 Received in revised form 25 November 2013 Accepted 9 December 2013 Available online 21 December 2013 This manuscript was handled by Peter K. Kitanidis, Editor-in-Chief, with the assisstance of Jean-Raynald de Dreuzy, Associate Editor Keywords: Shale gas reservoirs Multiple fractured horizontal wells Diffusion Desorption Stress-sensitivity Well test

s u m m a r y Gas flow in shales is believed to result from a combination of several mechanisms, including desorption, diffusion, viscous flow and the effect of stress-sensitivity of reservoir permeability. However, little work has been done in literature to simultaneously incorporate all these mechanisms in well testing models for shale gas reservoirs. This paper presents a new well testing model for multiple fractured horizontal wells (MFHW) in shale gas reservoirs with consideration of desorption, diffusive flow, viscous flow and stresssensitivity of reservoir permeability. Comparing with current well testing models for MFHW, the model presented here takes into consideration more mechanisms controlling shale gas flow, which is more in line with the actual reservoir situation. Laplace transformation, point source function, perturbation method, numerical discrete method and Gaussian elimination method are employed to solve the well testing model. The pressure transient responses are then inverted into real time space with Stehfest numerical inversion algorithm. Type curves are plotted, and different flow regimes in shale gas reservoirs are identified. The effects of relevant parameters are analyzed as well. The presented model can be used to interpret pressure data more accurately for shale gas reservoirs and provide more accurate dynamic parameters which are important for efficient reservoir development. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction There are multiple types of pores in shale gas reservoirs. According to Wang et al. (2009), four different types of pores are present in shale gas reservoirs: pores in organic matrix, pores in nonorganic matrix, natural fractures and hydraulic fractures. Organic-matter pores, ranging from a few nanometers to a few micrometers, are especially important because they can adsorb shale gas, as well as store free gas. Compared to the pores in conventional gas reservoirs, for the same pore volume, the exposed surface area in organic-matter pores is larger, thus they can absorb more shale gas. Generally speaking, pores in shale gas reservoirs can be classified as two major types like Warren–Root model (Warren and Root, 1963): (1) pores in shale matrix whose diameter is very small. Shale gas in this kind of pore space is mainly stored by adsorption, and gas flow is believed to be diffusive flow driven by concentration difference and (2) fracture which is not only storage space for free shale gas, but also connection between different pores. Like conventional gas reservoirs, shale gas flow in fractures is seepage flow driven by pressure difference.

⇑ Tel.: +86 13408600476. E-mail address: [email protected] 0022-1694/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2013.12.019

Kucuk and Sawyer (1980) first studied the pressure transient behavior of shale gas reservoirs by using analytical method and numerical simulation method. However, the analytical model presented in their paper did not take into account the effects of desorption and diffusion; the numerical model took into consideration the effect of desorption, but the effect of diffusive flow was still ignored. Some researchers (Bumb and McKee, 1988; Lane et al., 1989; Gao et al., 1994; Spivey and Semmelbeck, 1995) represented by Bumb and McKee (1988) proved that the desorption behavior of shale gas could be described by Langmuir isotherm theory based on experimental data. Carlson and Mercer (1991) investigated the behavior of gas flow in shale gas reservoir by coupling conventional dual-porosity model and the effects of desorption and diffusion. However, in this paper the hydraulic fractured vertical well in shale gas reservoirs was treated as a non-fractured vertical well with magnified well radius, thus pressure responses calculated by this model could not reflect characteristic flow regime for fractured wells, such as linear flow regime. In addition, the proposed model did not take into account the stress-sensitivity of natural fracture system. Ozkan et al. (2010) established a dual-mechanism dual-porosity model for shale gas reservoirs, taking into account the diffusive flow in shale matrix and the stress-sensitivity of natural fracture system. However, desorption of shale gas, which is an important

300

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

Nomenclature Bg C CD Cg Cgi D h k ki I0(x) K0(x) K1(x) Lref Lh M Mg n p pi psc ^ðtÞ q q⁄ qij qij ðsÞ qsc qsf ^D q r rm rw rD R Rm s S t tD

volume factor, dimensionless wellbore storage coefficient, m3 Pa1 dimensionless wellbore storage coefficient, dimensionless gas compressibility, Pa1 gas compressibility at initial condition, Pa1 diffusion coefficient, m2/s reservoir thickness, m permeability of natural fracture system, m2 permeability of natural fracture system at initial condition, m2 modified Bessel function of first kind, zero order modified Bessel function of second kind, zero order modified Bessel function of second kind, first order reference length, m length of horizontal well, m number of hydraulic fractures apparent molecular weight of shale gas, kg/kmol molar quantity of shale gas, kmol pressure of natural fracture system, Pa initial pressure of shale gas reservoirs, Pa pressure at standard condition, Pa surface production rate of the line sink, m3/s mass flow rate per unit reservoir between shale matrix and fracture, kg/(m3 s) flux density of the jth segment in the ith fracture, m3/ (s m) Laplace transformation of qij constant surface production rate of the multiple fractured horizontal well, m3/s sandsurface flow rate, m3/s dimensionless production rate of the line sink, dimensionless pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi radial distance, r ¼ x2 þ y2 , m radial distance in spherical matrix blocks, m well radius of horizontal well, m qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dimensionless radial distance, r D ¼ x2D þ y2D , dimensionless gas constant, J/(mol K) external radius of matrix block, m variable of Laplace transformation, dimensionless skin factor, dimensionless time, s dimensionless time, dimensionless

source of production in shale gas reservoirs, was ignored in their model. Freeman (2010) and Cipolla et al. (2010) employed a numerical simulator to study the flow regimes for shale gas reservoirs incorporating the effect of gas desorption, but the diffusive flow in shale matrix and stress-sensitivity of natural fracture system were not taken into account. Guo et al. (2012) established a well testing model for multistage fractured horizontal wells in shale gas reservoirs. In their model, the diffusion and desorption effects were considered, but the stress-sensitivity effect was not considered. The permeability of shale is ultralow, thus the stress-sensitivity of shale gas reservoirs should not be ignored. In addition, in Guo’s paper, hydraulic fractures were assumed to be perpendicular to the horizontal well, which is not always true in actual reservoirs because the minimum principal stress may not be parallel to the horizontal well. Multiple fractured horizontal well is proved to be the most effective well type for the development of shale gas reservoirs,

T Tsc v V VD VE Vi VL x, y yi

Dyi XfLi XfRi DXfi,j DXfDi,j Z

q qsc /

l li r aj b

w wL wi Dw Dws wD x k K

cD

reservoir temperature, K temperature at standard condition, K flow velocity of shale gas in natural fracture system, m/s volumetric gas concentration, sm3/m3 dimensionless gas concentration, dimensionless equilibrium volumetric gas concentration, sm3/m3 volumetric gas concentration at initial condition, sm3/ m3 Langmuir volume (at standard condition), sm3/m3 x- and y-coordinates, m y-coordinate of the intersection of the ith fracture and y-axis, m difference between yi and yi1, Dyi = yi  yi1 length of left wing of ith fracture, m length of right wing of ith fracture, m length of discrete segment (i, j), m dimensionless length of discrete segment (i, j), dimensionless Z-factor of shale gas, dimensionless shale gas density, kg/m3 shale gas density at standard condition, kg/m3 porosity, fraction gas viscosity, Pa s gas viscosity at initial condition, Pa s adsorption index, dimensionless angle between jth fracture and y-axis, degree a parameter related to permeability modulus, Pa1 s pseudo-pressure, Pa/s Langmuir pseudo-pressure, Pa/s pseudo-pressure at initial condition, Pa/s pseudo-pressure difference, Pa/s additional pseudo-pressure drop, Pa/s dimensionless pseudo-pressure, dimensionless storativity ratio, dimensionless interporosity flow coefficient, dimensionless total storage capacity, Pa1 dimensionless permeability modulus, dimensionless

Subscript D dimensionless Superscript  Laplace transform

and some work has been done to study the pressure transient dynamics of this kind of well type. Based on linear flow assumption, many researchers (Aboaba and Cheng, 2010; Bello and Watenbargen, 2010; Al-Ahmadi et al., 2010; Brohi et al., 2011; Ozkan et al., 2011) proposed linear flow model, linear composite model or tri-linear flow model to study production from a multiple fractured horizontal well in a shale gas reservoir. These models are easy to solve, however, they did not take into account desorption or diffusive flow which is typical in shale gas reservoirs. In addition, these models can only calculate the pressure responses of certain regime, such as early-time linear-flow regime; they can not reflect the complete pressure dynamics and flow regimes through the production of multiple fractured horizontal wells, such as pseudo-radial flow regime and interference between fractures. In view of this, this paper presents a semi-analytical model for multiple fractured horizontal wells in shale gas reservoirs which takes into consideration multiple flow mechanisms of shale gas,

301

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

such as desorption of adsorbed shale gas, diffusive flow in shale matrix, viscous flow and stress-sensitivity in natural fracture system. To account for the stress-sensitivity of natural fracture system caused by the closure of natural fracture during production, the stress-dependent permeability of natural fracture system is adopted. The model proposed here is more comprehensive and can give the complete pressure responses during production. It can be used to accurately interpret well testing data and thus to obtain important dynamic parameters for shale gas reservoirs.

2. Physical model for multiple fractured horizontal wells in shale gas reservoirs The schematic illustration in Fig. 1 show a horizontal well intercepted by multiple fractures (the number of fractures is M) in shale gas reservoirs. Fig. 2 is the profile of Fig. 1. The y-coordinate of the intersection of the ith fracture and y-axis is yi. Other assumptions of the presented mathematical model are as follows: (1) Shale gas reservoirs are dual-porosity systems, containing natural fractures and micro- and nano-pores in shale matrix. (2) Gas flow in natural fracture system can be described by Darcy’s law, and permeability of natural fractures is stressdependent. Considering ultralow permeability of shale matrix, gas flow in shale matrix is assumed to be diffusive flow driven by concentration difference instead of pressure difference. (3) The multiple fractured horizontal well produces at a constant rate qsc; however, the flow rate of each fracture is distinct, and even for the same fracture the flow rate is also a function of coordinates. (4) The effects of gravity and capillary force are negligible. (5) The initial pressure throughout the reservoir is uniform which equals to pi. (6) Hydraulic fractures may not always be perpendicular to the horizontal well. There may be diagonal fractures under control of different ground stress. It is very difficult to directly establish a well testing model for multiple fracture horizontal wells in shale gas reservoirs due to multiple controlling mechanisms and complex inner boundary condition. Source function method is employed in this paper to solve the problem. Based on the basic line sink solution and by discretizing hydraulic fractures, the pressure responses of multiple fractured horizontal wells in shale gas reservoirs can be obtained.

3. Basic line sink solutions in shale gas reservoirs 3.1. Seepage model for shale gas in fracture system 3.1.1. Partial differential equation Shale gas flow in natural fracture system is described by the following partial differential equation which can be obtained by combining mass conservation equation, equation of motion and equation of state (derivation details are given in Appendix A):

 2   @ 2 w 1 @w @w 2psc T @V ðwi wÞb /li C gi @w þ þ b ¼ e þ @r 2 r @r @r @t ki T sc @t ki

where w is the pseudo-pressure which is defined by Eq. (A.10) in Appendix A, Pa/s; r is the radial distance, m; b is a parameter related to permeability modulus, a characteristic parameter of stress-sensitivity, which is defined by Eq. (A.12) in the Appendix A, Pa1 s; wi is the initial pseudo-pressure, Pa/s; / is the reservoir porosity, fraction; li is the gas viscosity at initial condition, Pa s; Cgi is the gas compressibility at initial condition, Pa1; ki is the permeability of natural fracture system at initial condition, m2; t is the production time, s; psc is the pressure at standard condition, Pa; T is the reservoir temperature, K; Tsc is the temperature at standard condition, K; V is the average gas concentration in shale matrix, s m3/m3. 3.1.2. Initial and boundary conditions The inner boundary condition for a line sink can be expressed as following:

limeðwi wÞb r e!0

 ^ðtÞT @w p q ¼ sc @r r¼e pkhT sc

ð2Þ

^ðtÞ is production rate of the line sink, s m3/s; k is the permewhere q ability of natural fracture system, m2; h is the reservoir thickness, m. The outer boundary condition for an infinitely large shale gas reservoir is:

wjr!1 ¼ wi

ð3Þ

The initial condition for natural fracture system is:

wjt¼0 ¼ wi

ð4Þ

Eqs. (1)–(4) make up the seepage model for natural fracture system in shale gas reservoirs. 3.1.3. Dimensionless mathematical model For convenience of calculation, some dimensionless variables are defined which are shown in Table 1. With the definition in Table 1, the dimensionless forms of Eqs. (1)–(4) can be obtained as follows:

z 1

XfLi

i

M Horizontal well

XfRi

y

O

x

h

ð1Þ

Hydraulic fracture

Fig. 1. Physical model for MFHW in shale gas reservoirs.

302

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

z

Impermeable top boundary

h

y1

yi

Horizontal well

yM

y

x Impermeable bottom boundary Fig. 2. Profile of impermeable top and bottom boundaries.

Table 1 Definitions of dimensionless variables. Dimensionless pseudo-pressure

wD ¼ ppkiqhT scT Dw, Dw ¼ wi  w

Dimensionless radial distance

r r D ¼ Lref

Dimensionless gas concentration difference Dimensionless time

VD = Vi  V

Storativity ration

x ¼ /CKgi cD ¼ b ppsckiqhTscscT

sc sc

ki t (where K ¼ /C gi þ 2qpkli hÞ tD ¼ l K 2 sc i i Lref

Dimensionless permeability modulus Dimensionless production rate of the line sink

^D ðt D Þ ¼ q ^ðtÞ=qsc ^D ¼ q q

Note: Lref is an arbitrary reference length which can be rw, Lh, etc.

 2   @ 2 wD 1 @wD @wD @w @V D þ  cD ¼ ecD wD x D þ ð1  xÞ 2 r D @r D @r D @tD @t D @r D   @w ^D ¼ q lim ewD bD r D D  eD !0 @r

ð5Þ ð6Þ

D r D ¼ eD

wD jrD !1 ¼ 0

ð7Þ

wD jtD ¼0 ¼ 0

ð8Þ

Eqs. (5)–(8) make up the dimensionless seepage model for natural fracture system in shale gas reservoirs. 3.1.4. Perturbation technique and Laplace transformation The dependence of permeability of natural fracture system on pore pressure makes the above seepage model strongly nonlinear. With Pedrosa’s substitution (Pedrosa, 1986),

wD ðr D ; t D Þ ¼ 

1

cD

ln½1  cD nD ðr D ; tD Þ

ð15Þ ð16Þ

D r D ¼ eD

nD0 jrD !1 ¼ 0

ð17Þ

nD0 jtD ¼0 ¼ 0

ð18Þ

Laplace transformation with respect to tD is then adopted to invert the above seepage model into Laplace domain.

nD0 ¼ VD ¼

Z

þ1

0 Z þ1

nD0 estD dt D

ð19Þ

V D estD dtD

ð20Þ

0

2

ð10Þ

Taking the Laplace transformation of Eqs. (15)–(18), one can obtain the following system of equations:

ð11Þ

d nD0 1 dnD0 þ ¼ xsnD0 þ ð1  xÞsV D r D drD dr 2D  dn  ^D ¼ q lim r D D0  eD !0 dr 

D r D ¼eD

nD jrD !1 ¼ 0

ð12Þ

nD jtD ¼0 ¼ 0

ð13Þ

The strong nonlinearities in Eqs. (5) and (6) are considerably weakened in Eqs. (10) and (11). Performing a parameter perturbation in cD by defining the following series

nD ¼ nD0 þ cD nD1 þ c2D nD2 þ    1 1 ln½1  cD nD ðr D ; t D Þ ¼ nD ðrD ; t D Þ þ cD n2D ðr D ; t D Þ þ     2 cD 1 ¼ 1 þ cD nD ðr D ; t D Þ þ c2D nD ðr D ; tD Þ þ    1  cD nD ðr D ; tD Þ

@ 2 nD0 1 @nD0 @n @V D þ ¼ x D0 þ ð1  xÞ r D @r D @tD @t D @r 2D  @nD0  ^D ¼ q lim r D eD !0 @r 

ð9Þ

After some algebraic manipulations, Eqs. (5)–(8) become

@ nD 1 @nD 1 @n @V þ ¼ x D þ ð1  xÞ D rD @r D 1  cD nD @t D @tD @r 2D  @nD  ^D ¼ q lim rD eD !0 @r 

and considering that the value of cD, dimensionless permeability modulus, is always small, the zero-order perturbation solution can satisfy accuracy requirement. Thus Eqs. (10)–(13) become

2

ð21Þ ð22Þ

D r D ¼ eD

nD0 jrD !1 ¼ 0

ð23Þ

Eqs. (21)–(23) represent the zero-order dimensionless seepage model for natural fracture system in shale gas reservoirs.

ð14aÞ 3.2. Transportation model for shale gas in matrix system

ð14bÞ ð14cÞ

Micro-scale and nano-scale pores account for a large proportion in shale gas reservoir. On the one hand, the large exposed surface area in these pores is the main storage space for adsorbed shale

303

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

gas due to their small pore diameters. On the other hand, gas flow in these pores is different from the Darcy flow because of their ultralow permeability. It is commonly recognized that shale gas flow in micro-scale and nano-scale pores in shale matrix is diffusive flow driven by concentration difference, which can be describe by Fick’s law. 3.2.1. Pseudo-steady diffusive flow model For convenience, in demonstrating the derivation of the well testing model, this paper considers spherical geometry for the shale matrix blocks. However, the spherical matrix geometry chosen here is not a necessity for the mathematical model; the well testing model and solution presented in this paper can be easily extended to other matrix shapes. According to the Fick’s first law, pseudo-steady diffusive flow in shale matrix is as follows:

@V 6Dp2 ¼ 2 ðV E  VÞ @t Rm

ð24Þ

where VE is the equilibrium volumetric gas concentration, s m3/m3 2 l KL2 Defining V ED ¼ V i  V E , k ¼ i k sref , s ¼ 6DRmp2 and with the definii tion of VD in Table 1, Eq. (24) becomes

@V D ¼ kðV ED  V D Þ @t D

ð25Þ

The desorption behavior of shale gas can be described by Langmuir isotherm:

w wL þ w wi Vi ¼ VL wL þ wi

VE ¼ VL

ð26Þ ð27Þ

where VL is the Langmuir volume, s m3/m3; wL is the Langmuir pseudo-pressure, Pa/s. With Eqs. (26) and (27) and the definition of VED, one can get

V ED ¼ V i  V E ¼ rwD

ð28Þ

In Eq. (28) r is adsorption index, representing the effect of desorption of shale gas, and it is defined as

r¼

psc qsc T

wL V L

pkhT sc ðwL þ wÞ½wL þ wi 

where D is the Diffusion coefficient, m2/s. The partial differential equation describing transient state diffusive flow in shale matrix can be obtained by using Fick’s second law:

  1 @ @V @V 2 ¼ r D m r2m @r m @r m @t

The diffusive flow of shale gas in the spherical matrix blocks is symmetric, thus the center of matrix blocks can be treated as a no-flow boundary, which gives the following inner boundary condition:

 @V  ¼0 @r m rm !0

ð35Þ

Gas concentration on the external surface of the matrix blocks is evaluated at the gas pressure in the natural fracture system. So the outer boundary condition can be written as

Vjrm ¼Rm ¼ V E

ð36Þ

Initially,

Vjt¼0 ¼ V i

ð37Þ

With the dimensionless parameters given in Table 1 and rmD ¼ Rrmm , Eqs. (33)–(37) can be written in the following dimensionless form:

 @V D @V D  ¼ 3k @t D @r mD rmD ¼1   1 @ @V D 1 @V D 2 ¼ r mD k @t D @r mD r 2mD @r mD  @V D  ¼0 @r mD rmD !0

ð38Þ ð39Þ ð40Þ

V D jrmD ¼1 ¼ V ED

ð41Þ

V D jtD ¼0 ¼ 0

ð42Þ

l KL2

where k ¼ i ki sref and s ¼ RD2 . m The solution of Eqs. (38)–(42) in the Laplace domain is (see Appendix B for details)

ð29aÞ VD ¼

rcan be considered as a constant within the common pressure range of interest and can be evaluated at w = wi. Eq. (29a) becomes p q T wL V L r ¼ sc sc pkhT sc ðwL þ wi Þ½wL þ wi 

ð34Þ

r 

sh

pffiffiffiffiffiffiffi s=k

pffiffiffiffiffiffiffi sh s=kr mD r mD

wD

Combination of Eqs. (43) and (38) gives

ð29bÞ sV D ¼ 3k

Taking Laplace transformation of Eqs. (25) and (28) with respect to tD, one can obtain the following equations

 dV D   dr mD 

¼ 3rk

pffiffiffiffiffiffiffi hpffiffiffiffiffiffiffi i s=k coth s=k  1 wD

ð30Þ

3.3. Natural fracture-shale matrix coupling model

V ED ¼ rwD

ð31Þ

3.3.1. Pseudo-steady state coupling model Substitution of Eq. (32) into Eq. (21) yields

VD ¼

rk sþk

wD

ð32Þ

ð44Þ

r mD ¼1

sV D ¼ kðV ED  V D Þ

Substituting Eq. (31) into Eq. (30) and after some algebraic manipulations, one can obtain

ð43Þ

2

d nD0 1 dnD0 rkð1  xÞ þ ¼ xsnD0 þ swD r D dr D sþk dr 2D

ð45Þ

By using perturbation technique, the zero-order form of Eq. (45) is as following:

3.2.2. Transient diffusive flow model While flow of shale gas in natural fracture system is still represented by Eq. (1), the internal source item @V representing the effect @t of desorption of shale gas is given by the following equation based on the spherical matrix blocks assumption:

d nD0 1 dnD0 rkð1  xÞ þ ¼ xsnD0 þ snD0 r D dr D sþk dr 2D

 @V 3D @V  ¼ @t Rm @rm rm ¼Rm

f ðsÞ ¼ xs þ

ð33Þ

2

ð46Þ

If we define the function f(s) by

rkð1  xÞ sþk

s

ð47Þ

304

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

Then Eq. (46) can be rewritten as

4. Pressure transient characteristics for multiple fractured horizontal wells in shale gas reservoirs

2

d nD0 1 dnD0 þ ¼ f ðsÞnD0 r D dr D dr 2D

ð48Þ 4.1. Pressure responses

Eq. (48) is the pseudo-steady state coupling differential equation of natural fracture system and shale matrix. The general solution of Eq. (48) is

nD0

qffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffi  f ðsÞr D þ BI0 f ðsÞr D ¼ AK 0

ð49Þ

With the outer boundary condition Eq. (23) and the theory of Bessel’s function, the coefficient in Eq. (49), B, should equal zero, thus one can get

nD0 ¼ AK 0

qffiffiffiffiffiffiffiffi  f ðsÞr D

^D K 0 nD0 ¼ q where ( f ðsÞ ¼

^D ¼ q

ð51Þ

Noting that xK1(x)jx?0 = 1, we can obtain the expression of the coefficient A.

^D A¼q

ð52Þ

ð60Þ

qffiffiffiffiffiffiffiffi  ^D K 0 f ðsÞr D ¼q

nD ¼

ð53Þ

Eq. (53) is the basic zero-order line sink solution of the pseudosteady state model for shale gas reservoirs. 3.3.2. Transient state coupling model Substituting Eq. (44) into Eq. (21), one can obtain

qffiffiffiffiffiffiffiffi  ^ q f ðsÞr D K0 qsc

ð62Þ

ð54Þ

Similarly, we can get the zero-order form of Eq. (54) by using perturbation technique: 2 pffiffiffiffiffiffiffi hpffiffiffiffiffiffiffi i d nD0 1 dnD0 s=k  1 nD0 þ ¼ xsnD0 þ 3rkð1  xÞ s=k coth 2 rD dr D drD

ð56Þ

Eq. (55) becomes

(

(

^xi;j ¼  2N2jþ1 X fLi sin ai 2N

^xi;j ¼ 2ðjNÞ1 X fRi sin ai 2N ^i;j ¼ yi  2ðjNÞ1 X fRi cos ai y 2N

ð64Þ

;

16j6N

ð65Þ

;

N þ 1 6 j 6 2N

ð66Þ

The coordinates of the jth end point in the ith fracture (xi,j, yi,j), are:

(

xi;j ¼  Njþ1 X fLi sin ai N yi;j ¼ yi þ Njþ1 X fLi cos ai N

ð57Þ

ð58Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxD  xwD Þ2 þ ðyD  ywD Þ2

^i;j ¼ yi þ 2N2jþ1 X fLi cos ai y 2N

2

Eq. (57) is the transient state coupling differential equation of natural fracture system and shale matrix. It can be observed that Eqs. (57) and (48) have identical expressions except that the expressions of f(s) are distinct in these two equations. That implies that one can easily get the solution for Eq. (57) just by replacing the expression of f(s) in the solution for Eq. (48), i.e., Eq. (53). Thus the basic zero-order line sink solution of the transient state model for shale gas reservoirs can be obtained as follows:

ð63Þ

The pressure responses of multiple fractured horizontal wells in shale gas reservoirs are obtained by discretizing fractures and applying superposition principle. As shown in Fig. 3, the top view of Fig. 1, both the left wing and right wing of each fracture are divided into N equal segments, and then we can get 2N discrete segments in each fracture. The coordi^i;j Þ, nates of midpoint for the jth segment in the ith fracture, ð^ xi;j ; y are:

ð55Þ The function f(s) in the transient state coupling model is defined as following

qffiffiffiffiffiffiffiffi  ^ q K0 f ðsÞRD qsc

where

RD ¼

2

pffiffiffiffiffiffiffi hpffiffiffiffiffiffiffi i f ðsÞ ¼ xs þ 3rkð1  xÞ s=k coth s=k  1

ð61Þ

nD represents the zero-order solution. where  ^D in Table 1, Eq. (61) becomes With the definition of q

nD ¼

d nD0 1 dnD0 þ ¼ xsnD0 þ 3rkð1 r D dr D dr 2D pffiffiffiffiffiffiffi hpffiffiffiffiffiffiffi i s=k  1 wD  xÞ s=k coth

qffiffiffiffiffiffiffiffi  f ðsÞr D

If the line sink is located in (xw, yw) instead of the center of the shale gas reservoir, Eq. (62) can be rewritten in a more general form:

Substitution of Eq. (52) into Eq. (50) yields

nD0

xÞ xs þ rkð1 s; pseudo-steady diffusive flow sþk pffiffi

pffiffis  f ðsÞ ¼ xs þ 3rkð1  xÞ ks coth  1 ; transient diffusive flow k

^D K 0 nD ¼ q

r D !0

qffiffiffiffiffiffiffiffi  ^D K 0 f ðsÞr D ¼q

ð59Þ

For convenience, Eq. (59) is rewritten as

 qffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffi   f ðsÞr D  A  f ðsÞr D K 1

d nD0 1 dnD0 þ ¼ f ðsÞnD0 r D dr D dr 2D

qffiffiffiffiffiffiffiffi  f ðsÞr D

ð50Þ

Substitution of Eq. (50) into Eq. (22) yields

nD0

For a line sink in a laterally infinite shale gas reservoir, the basic zero-order solution is obtained in the previous section, which is

(

xi;j ¼ jN1 X fRi sin ai N yi;j ¼ yi  jN1 X fRi cos ai N

;

16j6N

ð67Þ

;

N þ 1 6 j 6 2N þ 1

ð68Þ

Dyi ¼ yi  yi1

ð69Þ

The mathematical equation of the ith fracture is:

ywi ¼ ctgðai Þxwi þ yi

ð70Þ

Eq. (70) can be written in dimensionless form as following:

ywDi ¼ ctgðai ÞxwDi þ yDi where xwDi = xwi/Lref, ywDi = ywi/Lref and yDi = yi/Lref.

ð71Þ

305

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

1

( xi , 2 , yi , 2 ) ( xi , N , yi , N )

α1

( xi , N+1 , yi , N+1 ) (0, y1)

i

( xi , 1 , yi , 1 )

( xi , N+2 , yi , N+2 )

M

( xˆi , 1 , yˆ i , 1 )

αj

αM

( xˆi , N , yˆ i , N ) (0, yi)

(0, yM)

( xˆi , N+1 , yˆ i , N+1 )

y

( xi , 2 N , yi , 2 N ) ( xi , 2 N+1 , yi , 2 N+1 )

x

( xˆi , 2 N , yˆ i , 2 N )

Fig. 3. Schematic of discretion of a multiple fractured horizontal well.

Each discrete segment can be seen as a surface source with a uniform flux density qij, and the pressure responses in shale gas reservoir caused by the jth segment in the ith fracture can be obtained by integrating Eq. (63) over the surface source:

nDi;j ðxD ; yD Þ ¼

Z

i;j ðsÞ q qsc

K0

qffiffiffiffiffiffiffiffi  f ðsÞRD ðxD ; yD ; xwD ; ywD Þ dl

ð72Þ

Ci;j

where

RD ðxD ; yD ; xwD ; ywD Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dl ¼ dx þ dy

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxD  xwD Þ2 þ ðyD  ywD Þ2

By writing Eq. (82) for midpoints for all discrete segments (k = 1, 2, . . . , M, t = 1,2 . . . , 2N), we can get 2N  M equations with Di;j ðuÞ and nwD . In order to suc(2N  M + 1) unknowns which are q cessfully obtain the pressure responses, one more equation is necessary. The assumption of constant flow rate gives the following equation that must hold at any time.

ð76Þ

K0

pffiffi

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sRD ðxD ;yD ;xwDi Þ 1 þ ctg2 ðai ÞdxwDi ð77Þ

xDi;j

where

RD ðxD ; yD ; xwDi Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxD  xwDi Þ2 þ ½yD þ ctgðai ÞxwDi  yDi 2

ð78Þ

By applying the principle of superposition, the pressure responses at (xD, yD) caused by 2N  M segments can be obtained as follows:

nD ðxD ; yD Þ ¼

M X 2N X nDi;j ðxD ; yD Þ

ð79Þ

i¼1 j¼1

Thus, the pressure response at the discrete (k, v)(1 6 k 6 M,1 6 v 6 2N) can be obtained as follows:

^Dk;t Þ ¼ nD ð^xDk;t ; y

M X 2N X i¼1 j¼1

ð82Þ

i¼1 j¼1

M X 2N X q i;j DLfi;j  ¼ sc ½q s i¼1 j¼1

ð83Þ

where DLfi,j is length of the discrete segment (i, j). With some algebraic manipulations, Eq. (83) becoms

 M X 2N  X i;j Lref DLfi;j q 1 ¼ s q L ref sc i¼1 j¼1

ð84Þ

Or more simply,

Then Eq. (75) becomes xDi;jþ1

^Dk;t Þ nDi;j ð^xDk;t ; y

ð74Þ

The dimensionless flux per unit length of discrete segment (i, j) is defined as

qi;j ðsÞLref qsc

M X 2N X

nwD ¼

qffiffiffiffiffiffiffiffi  Z i;j ðsÞ xi;jþ1 q f ðsÞRD ðxD ; yD ; xwD ; ywD Þ K0 qsc xi;j Z qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i;j ðsÞLref xDi;jþ1 q 1 þ ctg2 ðai Þdxwi ¼ qsc xDi;j qffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f ðsÞRD ðxD ; yD ; xwD ; ywD Þ K0 1 þ ctg2 ðai ÞdxwDi ð75Þ

Z

ð81Þ

ð73Þ

nDi;j ðxD ; yD Þ ¼

Di;j ðsÞ nDi;j ðxD ;yD Þ ¼ q

^Dk;t Þ ¼ nwD nD ð^xDk;t ; y Thus, Eq. (80) becomes

With Eqs. (70) and (71), the integration in Eq. (72) becomes

qDi;j ðsÞ ¼

Permeability of shale gas reservoirs is much smaller than that of hydraulic fractures, which causes the flowing resistance in hydraulic fractures is much smaller than that in shales. Compared to low permeable shales, hydraulic fractures can be considered as infinite conductive fractures, thus we can have

^Dk;t Þ nDi;j ð^xDk;t ; y

segment

ð80Þ

M X 2N X

Di;j DLfDi;j  ¼ ½q

i¼1 j¼1

1 s

ð85Þ

Eqs. (82) and (85) compose a (2N  M + 1)-order system of linear algebraic equations which can be expressed in the following matrix form.

EX ¼ F

ð86Þ

Integrations in the coefficient matrix, E, can be calculated by numerical integration method. By applying Gaussian elimination method, we can get the pressure responses of multiple fractured horizontal wells in shale gas reservoirs. The nwD in Eq. (82) do not take into consideration the effects of wellbore storage and skin. Since the solution is given in Laplacetransform domain, wellbore storage and skin effects can be easily added as a final step. According to Van Everdingen (1953), the effects of wellbore storage and skin can be incorporated in the above pressure solutions by Kucuk and Ayestaran (1985) and Thomas and Gringarten (2009):

306

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

nwDH ¼

snwD þ S s þ C D s2 ðsnwD þ SÞ

ð87Þ

where S is the skin factor defined by Eq. (88), dimensionless; CD is the dimensionless wellbore storage coefficient defined by Eq. (89).

S¼

pki hT sc

ð88Þ

Dws psc qsf T C CD ¼ 2phKL2ref

ð89Þ

where qsf is the subsurface flow rate, m3/s. Dws is a extra pseudopressure drop, Pa/s; C is the wellbore storage coefficient, m3 Pa1; nwDH in Eq. (87) is zero-order dimensionless bottom hole pressure in the Laplace domain, and by Stefest numerical inversion algorithm we can get the zero-order pressure responses nwDH(tD) in real time domain. Then with Eq. (90), we can obtain the bottom hole pressure for MFHW in shale gas reservoirs which takes into account the stress-sensitivity of natural fracture system.

wwD ðtD Þ ¼ 

1

cD

ln½1  cD nwDH ðtD Þ

ð90Þ

4.2. Type curves and discussions In this section, the dimensionless pseudo-pressure and derivative responses for a multiple fractured horizontal well in shale gas reservoir is calculated with the well testing model proposed above. Type curves are plotted and the effects of relevant parameters on the type curves are analyzed. The upper sets of lines in Figs. 4 and 6–11 are dimensionless pressure curves, and the lower sets of lines are dimensionless pressure derivative curves. Fig. 4 shows the type curves for a multiple fractured horizontal well in shale gas reservoirs with different values of dimensionless permeability modulus, cD. In addition to dimensionless permeability modulus, the values of other relevant parameters in Fig. 4 are listed as follows: CD = 10, S = 0.3, M = 3, XfLi = XfRi = 30 m, Dyi = 400 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, x = 0.12, k = 1E9, r = 1.2. The upper sets of lines represent dimensionless pressure curves, while the lower sets of lines represent dimensionless pressure derivative curves. It can be seen that there are nine possible flow regimes for a multiple fractured horizontal well in shale gas reservoirs:

(1) Pure wellbore storage period. During this period, gas flow is mainly dominated by wellbore storage effect, and both the dimensionless pressure and pressure derivative curves align in a unit slope trend. (2) Transition flow period after the wellbore storage period, which is characterized by a hump in the pressure derivative curves. (3) Early-time linear flow period perpendicular to each fracture (Fig. 5a), which is marked by a half slope trend in pressure derivative curve. During this period, each fracture produces independently, and the interference between different fractures has not been felt yet. (4) Transition flow period between early-time linear flow period and intermediate-time pseudo-radial flow period. During this period, type curves with different values of cD begin to separate from each other. (5) Intermediate-time pseudo-radial flow around each fracture (Fig. 5b). A constant pressure derivative on the type curve can be observed, and the value of the pressure derivative is ‘‘1/2M’’ where M is the number of hydraulic fractures. The existence of this period depends primarily on the fracture length and fracture spacing. If fracture spacing is much smaller compared to fracture length, this period may not be observed in the type curves. (6) Intermediate-time linear flow period (Fig. 5c). During this period, pressure wave reaches adjacent fractures, and the interference between different fractures gradually becomes obvious. Two parallel straight lines can be observed in type curves. However, due to the effect of stress-sensitivity of natural fracture system, the value of the slope of the straight lines is a little larger then 0.5, and the slope increases with increasing value of cD. (7) Transition flow period after the Intermediate-time linear flow period. During this period, the departure between type curves with different values of cD becomes obvious. The bigger the dimensionless permeability modulus, the bigger the departure between the type curve for cD – 0 from the type curve for cD = 0. (8) Diffusive flow period, representing gas diffusion from shale matrix to natural fracture system, which is marked by a ‘‘dip’’ in the pressure derivative curve. Due to the pressure difference between natural fracture system and shale matrix, the adsorbed gas in shale matrix starts to desorb from the surface of solid material and then diffuse from

1E+2

1E+1

1E+0

1E-1

1E-2 1E-2

1E-1

1E+0

1E+1

1E+2

1E+3

1E+4

1E+5

1E+6

1E+7

1E+8

Fig. 4. Type curves with different dimensionless permeability modulus cD.

1E+9

1E+10

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

Fig. 5a. Early-time linear flow.

Fig. 5b. Intermediate-time pseudo-radial flow.

Fig. 5c. Intermediate-time linear flow.

Fig. 5d. Late-time pseudo-radial flow.

shale matrix to natural fracture system driven by concentration difference. Stress-sensitivity mainly affects the shape of the ‘‘dip’’ in type curves during this period. The smaller the value of cD, the deeper the ‘‘dip’’. (9) Late-time pseudo-radial flow period (Fig. 5d). For no stresssensitivity case, the derivative curve during this period is a horizontal line with a value of 0.5; however, considering the effect of stress-sensitivity, the derivative curve is no longer a horizontal line but bends upward. It can be observed that the bigger the value of cD, the bigger the slope of the derivative curves, representing more pressure loss. Fig. 6 shows the effect of storativity ratio, x, on pressure transient responses of a multiple fractured horizontal well in shale gas reservoirs. In addition to storativity ratio x, the values of other relevant parameters in Fig. 6 are listed as follows: CD = 10, S = 0.3, M = 3, XfLi = XfRi = 30 m, Dyi = 400 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, cD = 0.1, k = 1E9, r = 1.2. The storativity ratio has significant effect on the early-time linear flow period, intermediate-time pseudo-radial flow around each fracture, intermediate-time linear flow period and diffusive flow period. The smaller the storativity ratio, the deeper and wider the ‘‘dip’’ in derivative curve, and the higher the derivative curve during early-time linear flow period and intermediate-time linear flow period. Fig. 7 shows the effect of adsorption index, r, on pressure transient responses of a multiple fractured horizontal well in shale gas reservoirs. In addition to adsorption index r, the values of other relevant parameters in Fig. 7 are listed as follows: CD = 10, S = 0.3, M = 3, XfLi = XfRi = 30 m, Dyi = 400 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, x = 0.12, k = 1E9, cD = 0.1. It can be seen that adsorption index has obvious effect on diffusive flow period. The larger the value

307

of adsorption index, the deeper and wider the ‘‘dip’’ in the derivative curve, meaning less pressure loss in reservoir. This is because that adsorption index represents the amount of absorbed gas in shale matrix, and larger r means more adsorbed gas, so more desorbed gas would diffuse from shale matrix to natural fracture system to compensate for the pressure loss in natural fracture system. Fig. 8 shows the effect of interporosity flow coefficient, k, on pressure transient responses of a multiple fractured horizontal well in shale gas reservoirs. In addition to interporosity flow coefficient k, the values of other relevant parameters in Fig. 8 are listed as follows: CD = 10, S = 0.3, M = 3, XfLi = XfRi = 30 m, Dyi = 400 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, x = 0.05, cD = 0.1, r = 1.2. Fig. 8 indicates that the larger the value of interporosity flow coefficient, the earlier the occurrence time of diffusive flow period. This is because interporosity flow coefficient is directly proportional to the diffusion coefficient of shale gas, and larger value of k corresponds to larger value of diffusion coefficient, D, so the reflection time of diffusive flow is relatively earlier. Fig. 9 shows the effect of half length of fractures, Xfi, on pressure transient responses of a multiple fractured horizontal well in shale gas reservoirs. In addition to Xfi, the values of other relevant parameters in Fig. 9 are listed as follows: CD = 10, S = 0.3, M = 3, Dyi = 400 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, x = 0.05, k = 109, cD = 0.1, r = 1.2. It can be seen that Xf has primary effect on early-time and intermediatetime pressure responses. For type curves with smaller the half length of fractures, the position of the derivative curves during early-time linear flow period is higher, and the duration time of intermediate-time pseudo-radial flow period is longer. Fig. 10 shows the effect of fracture spacing, Dyi, on pressure transient responses of a multiple fractured horizontal well in shale gas reservoirs. In addition to Dyi, the values of other relevant parameters in Fig. 10 are listed as follows: CD = 10, S = 0.3, M = 3, XfLi = XfRi = 30 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, x = 0.12, k = 109, cD = 0.1, r = 1.2. We can observe that different values of Dyi mainly have effect on the existence and duration time of intermediate-time pseudo-radial flow period. When the fracture length is kept constant, larger fracture spacing corresponds to longer intermediate-time pseudo-radial flow period. It should be pointed out that when fractures are close to each other enough, the intermediate-time pseudo-radial flow period may not be observed in type curves, such as the Dyi = 100 m case in Fig. 10. Fig. 11 is comparison of type curves calculated by different models proposed in this paper, i.e., pseudo-steady model and transient model. The values of relevant parameters in Fig. 11 are listed as follows: CD = 10, S = 0.3, M = 3, XfLi = XfRi = 30 m, Dyi = 500 m, a1 = a2 = a3 = 75°, Lref = 0.1 m, x = 0.12, k = 109, cD = 0.1, r = 1.2. It can be seen that type curves calculated by different models yield almost same pressure responses at early and late times. The intermediate-time pressure responses calculated by these two models, however, are quite distinct. Type curves calculated by pseudo-steady model display a characteristic ‘‘dip’’ in derivative curves, while type curves calculated by transient model only shows a sign of flat concave. Another conclusion can be drawn from Fig. 11 is that the occurrence time of diffusive flow period is a little earlier when calculated by transient model. That is because the desorption of shale gas in shale matrix is more sensitive to pressure change in transient model.

308

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

Fig. 6. Effect of storativity ratio x on type curves.

Fig. 7. Effect of adsorption index r on type curves.

Fig. 8. Effect of interporosity flow coefficient k on type curves.

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

Fig. 9. Effect of half length of fractures Xf on type curves (Xfi = XfLi = XfRi, and Xf1 = Xf2 = Xf3).

Fig. 10. Effect of fracture spacing Dyi on type curves (Dy1 = Dy2 = Dy3).

Fig. 11. Comparison of type curves calculated by pseudo-steady model and transient model.

309

310

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

(4)

(3) (1)

(2)

Fig. 12. Type curves calculated by a simplified case of the complex model proposed in this paper.

4.3. Simplified model and comparison In the case x = 1, Eq. (47) is simplified into f(s) = s, and the model presented in the above sections is reduced into a simplified well testing model for MFHWs in shale gas reservoirs which takes into consideration stress-sensitivity, but the existence of natural fracture and the effects of desorption and diffusive flow are not taken into account. If we further set dimensionless permeability modulus cD equal to 0, the model presented in the above sections will be further simplified into a well testing model for MFHWs in conventional gas reservoirs. To ensure the validity of the model presented in this paper, type curves for MFHWs in conventional gas reservoir are calculated with the simplified case of our model. Fig. 12 shows the calculated type curves without consideration of wellbore storage effect and skin effect. The values of relevant parameters in Fig. 12 are listed as follows: CD = 0, S = 0, M = 2, XfLi = XfRi = 20 m, Dyi = 200 m, a1 = a2 = a3 = 90°, Lref = 20 m, x = 1, cD = 0. Four possible flow regimes can be observed in Fig. 12: (1) First linear flow period perpendicular to each fracture; (2) first pseudo-radial flow period around each fracture; (3) second linear flow period; and (4) second pseudo-radial flow period. Compared to type curves shown in Fig. 6 in Horne and Temeng’s paper (1995), the identification of flow regimes and corresponding characteristics of different flow regimes are in good agreement. The verification justified the validity of the well testing model presented in this paper. 5. Conclusions (1) Gas flow in shales is a result of multiple transportation mechanisms, including desorption, diffusion, viscous flow and stress-sensitivity of reservoir permeability. So far, little work has been done in literature to simultaneously incorporate all these mechanisms in well testing models for multiple fracture horizontal wells (MFHW) in shale gas reservoirs. This paper presents a new well testing model for MFHWs in shale gas reservoirs with consideration of desorption, diffusive flow, viscous flow and stress-sensitivity of reservoir permeability. Compared to the existing well testing models for shale gas reservoirs, the model presented in this paper takes into consideration more transportation mechanisms and thus is more in accordance with actual seepage situation in shales. (2) The mathematical model is solved successfully. Type curves are plotted, and different flow regimes in shale gas

reservoirs are identified. The effects of relevant parameters are analyzed as well. Comparison between type curves calculated by the simplified case of our model and corresponding well testing model for conventional gas reservoirs shows good agreement, which also indicates the validity of our model. (3) The presented model can be used to interpret pressure data more accurately for shale gas reservoirs and provide more accurate dynamic parameters which are important for efficient reservoir development. Note: Visual Basic programs for plotting type curves based on the proposed model in this paper can be provided upon request.

Acknowledgements Support for this work is provided by Open Fund (PLN1211) of State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (Southwest Petroleum University), the National Natural Science Foundation of China (Grant No. 51304165), the National Science Fund for Distinguished Young Scholars of China (Grant No. 51125019), and the Doctoral Fund of Ministry of Education of China (Grant No. 20125121120003). We thank the Editors, for their timely and exemplary handing of the review process. We thank the reviewers for their detailed and very constructive comments on the manuscript. Appendix A. Derivation of partial differential equation for natural fracture system (Eq. (1)) Partial differential equation for natural fracture system in shale gas reservoirs is established by combining mass conservation equation, equation of motion and equation of state. (1) Mass conservation equation Mass conservation law on fracture system in shale gas reservoir yields the following equation in radial coordinates:

1 @ðr qv Þ @ðq/Þ ¼  q r @r @t

ðA:1Þ

where q⁄ is mass flow rate between shale matrix and fracture in unit-volume reservoir, and it can be expressed as

q ¼ qsc

@V @t

ðA:2Þ

311

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

where V is the volumetric gas concentration in shale matrix, s m3/ m3 (2) Equation of motion In Eq. (A.1), the radial velocity for gas flow in natural fracture system is given by

v

k @p ¼ l @r

where b is a parameter related to permeability modulus and is defined by Eq. (A.12), Pa1 s; Cg is the gas compressibility and is defined by Eq. (A.13), Pa1; psc is the pressure at standard condition, Pa; Tsc is the temperature at standard condition, K.



ðA:3Þ

where v is the flow velocity in natural fracture system, m/s; k is the permeability of natural fracture system, m2; p is the pressure of natural fracture system, Pa The increase in the effective overburden pressure as a result of production in shale gas reservoirs often results in the closing of natural fractures. To account for this, a stress-dependent permeability is adopted in the natural fracture flow model. Following Klkanl and Pedrosa (1991), the permeability modulus c which is used to describe the stress-sensitivity of natural fracture system is defined by

lZ

b¼

2p



c

ðA:12Þ

lZ in Eq. (A.12) can be treated as a constant when pressure is rela2p tively high.

Cg ¼

1 1 dp  p Z dZ

ðA:13Þ

Both gas viscosity l and gas compressibility Cg in Eq. (A.11) are functions of pressure p, however, in order to solve the model successfully, both l and Cg are evaluated at initial condition, i.e., l = l(pi) = li and Cg = Cg(pi) = Cgi. Eq. (A.11) becomes

ðA:4Þ

 2   @ 2 w 1 @w @w 2psc T @V ðwi wÞb /li C gi @w ðA:14Þ þ þ b ¼ e þ @r2 r @r @r ki @t ki T sc @t

By using the method of separation of variables and integration, Eq. (A.4) becomes

where li is the gas viscosity at initial condition, Pa s; Cgi is the gas compressibility at initial condition, Pa1. Eq. (A.14) is Eq. (1) in the main body of this paper.

c¼

Z

ki

k

1 dk k dp

1 dk ¼ k

Z

pi

cdp

ðA:5Þ

p

where pi is the reservoir pressure at initial condition, Pa; ki is the permeability of natural fracture system at initial condition, m2. Eq. (A.5) can be further written as

k ¼ ki ecðpi pÞ

ðA:6Þ

Eq. (A.6) is one of the most common relationships which are used to describe stress-dependent permeability. This definition corresponds to Type I rocks discussed by Raghavan and Chin (2004). Other relationships describing the stress-dependent permeability could also be used in the derivation of the mathematical model. Substituting Eq. (A.6) into Eq. (A.3), we can obtain

v¼

ki eðpi pÞc @p @r l

ðA:7Þ

Eq. (A.7) is the final expression of gas flow velocity in natural fracture system which takes into account the stress-sensitivity of permeability. (3) Equation of state Real natural gas can be described by the following equation of state:

pV ¼ nZRT

ðA:8Þ

Eq. (A.8) can also be written as

pM q¼ g ZRT

ðA:9Þ

From Section 3.3.2 we can obtain the complete transient diffusive flow model in shale matrix as follows:

 @V D @V D  ¼ 3k @t D @r mD rmD ¼1   1 @ @V D 1 @V D 2 ¼ r k @t D r 2mD @r mD mD @r mD  @V D  ¼0 @r mD rmD !0

Z

p

p0

2p dp lZ

ðA:10Þ

we can obtain the partial differential equation for natural fracture system in terms of pseudo-pressure as follows:

 2   @ 2 w 1 @w @w 2psc T @V ðwi wÞb /lC g @w ðA:11Þ þ þ b ¼ e þ @r 2 r @r @r ki @t ki T sc @t

ðB:1Þ ðB:2Þ ðB:3Þ

V D jrmD ¼1 ¼ V ED

ðB:4Þ

V D jtD ¼0 ¼ 0

ðB:5Þ

Taking Laplace transformations of Eqs. (B.1)–(B.5) one can obtain

 dV D  sV D ¼ 3k  dr mD  r mD ¼1   1 d dV D 1 2  sV D ¼ 0 r k r 2mD @r mD mD dr mD   dV D  ¼0  dr mD 

ðB:6Þ ðB:7Þ ðB:8Þ

r mD !0

V D jrmD ¼1 ¼ V ED

where Z is the compressibility factor of shale gas, dimensionless; Mg is the apparent molecular weight, kg/kmol; R is the universal gas constant, J/(mol K); T is the reservoir temperature, K. (4) Partial differential equation for natural fracture system With Eqs. (A.1), (A.2), (A.7) and (A.9) and the definition of pseudo-pressure defined in Eq. (A.10),

w ¼ wðpÞ ¼

Appendix B. Solution of transient diffusive flow model in matrix system

ðB:9Þ

Let E ¼ r D V D , then Eq. (B.7) becomes 2

d E s  E¼0 dr 2D k

ðB:10Þ

The general solution of Eq. (B.10) is

E ¼ A1 e

pffiffiffiffiffi

s=kr mD

þ A2 e

pffiffiffiffiffi

s=kr mD

ðB:11Þ

or

V D ¼ A1

e

pffiffiffiffiffi

s=kr mD

r mD

pffiffiffiffiffi e s=krmD þ A2 r mD

Substituting Eq. (B.12) into Eq. (B.8), we can obtain

ðB:12Þ

312

H.-T. Wang / Journal of Hydrology 510 (2014) 299–312

 dV D   dr mD 

r mD !0

pffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffi  s=ke s=krmD r mD  e s=krmD ¼ A1 r 2mD pffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffiffi s=kr mD e s=krmD  e s=krmD þ A2 ¼0 r 2mD

ðB:13Þ

To satisfy the inner boundary condition expressed as Eq. (B.13), the following relationship can be obtained

A2 ¼ A1

ðB:14Þ

Thus Eq. (B.12) can be rewritten as

e

V D ¼ A1

pffiffiffiffiffi

s=kr mD

r mD



e

pffiffiffiffiffi

s=kr mD

! ðB:15Þ

r mD

Substituting Eq. (B.15) into Eq. (B.9), we can obtain

A1

e

pffiffiffiffiffi

s=kr mD



r mD

e

pffiffiffiffiffi

s=kr mD

! ¼ V ED

r mD

ðB:16Þ

Substituting Eq. (31) in the main body of this paper into Eq. (B.16), we can obtains

r

pffiffiffiffiffi wD A1 ¼  pffiffiffiffiffi e s=k  e s=k

ðB:17Þ

Substituting Eq. (B.17) into Eq. (B.15), we can obtain the dimensionless gas concentration in shale matrix in the Laplace domain as follows:

VD ¼

sh

r 

pffiffiffiffiffiffiffi s=k

pffiffiffiffiffiffiffi s=kr mD sh r mD

wD

ðB:18Þ

where sh(x) is the hyperbolic sine function. We can further obtain the derivation of Eq. (B.18) as following:

 @V D   @r mD 

¼r

pffiffiffiffiffiffiffi hpffiffiffiffiffiffiffi i s=k coth s=k  1 wD

ðB:19Þ

r mD ¼1

Substituting Eq. (B.19) into Eq. (B.6), we can get

sV D ¼ 3k

 dV D   dr mD 

¼ 3rk

hpffiffiffiffiffiffiffi i pffiffiffiffiffiffiffi s=k cothð s=kÞ  1 wD

ðB:20Þ

r mD ¼1

References Aboaba, A., Cheng, Y., 2010. Estimation of fracture properties for a horizontal well with multiple hydraulic fractures in gas shale. In: SPE Paper 138524 Presented at the SPE Eastern Regional Meeting. Morgantown, West Virginia, USA.

Al-Ahmadi, H.A., Almarzooq, A.M., Watenbargen, R.A., 2010. Application of linear flow analysis to shale gas wells-field cases. In: SPE Paper 130370 Presented at the SPE Unconventional Gas Conference. Pittsburgh, Pennsylvania, USA. Bello, R.O., Watenbargen, R.A., 2010. Multi-stage hydraulically fractured horizontal shale gas well rate transient analysis. In: SPE Paper 126754 Presented at the North Arica Technical Conference and Exhibition. Cairo, Egypt. Brohi, I., Pooladi-Darvish, M., Aguilera, R., 2011. Modeling fractured horizontal wells as dual porosity composite reservoirs – application to tight gas, shale gas and tight oil cases. In: SPE Paper 144057 Presented at the SPE Western North American Region Meeting. Anchorage, Alaska, USA. Bumb, A.C., McKee, C.R., 1988. Gas-well testing in the presence of desorption for coalbed methane and Devonian shale. SPE Form. Eval. 3 (1), 179–185. Carlson, E.S., Mercer, J.C., 1991. Devonian shale gas production: mechanisms and simple models. J. Petrol. Technol. 43 (4), 476–482. Cipolla, C.L., Lolon, E.P., Erdle, J.C., et al., 2010. Reservoir modeling in shale-gas reservoirs. SPE Reser. Eval. Eng. 13 (4), 638–653. Freeman, C.M., 2010. A numerical study of microscale flow behavior in tight gas and shale gas. In: SPE Paper 141125 Presented at the SPE Annual Technical Conference and Exhibition. Florence, Italy. Gao, C., Lee, W.J., Spivey, J.P., et al., 1994. Modeling multilayer gas reservoirs including sorption effects. In: SPE Paper 29173 Presented at the SPE Eastern Regional Meeting. Charleston, West Virginia. Guo, J.J., Zhang, L.H., Wang, H.T., et al., 2012. Pressure transient analysis for multistage fractured horizontal wells in shale gas reservoirs. Trans. Porous Media 93 (3), 635–653. Horne, R.N., Temeng, K.O., 1995. Relative productivities and pressure transient modeling of horizontal wells with multiple fractures. In: SPE Paper 29891 Presented at the SPE Middle East Oil Show. Bahrain. Klkanl, J., Pedrosa Jr., O.A., 1991. Perturbation analysis of stress-sensitive reservoirs. SPE Form. Eval. 6 (3), 379–386. Kucuk, F., Ayestaran, L., 1985. Analysis of simultaneously measured pressure and sandface flow rate in transient well testing (includes associated papers 13937 and 14693). J. Petrol. Technol. 37 (2), 323–334. Kucuk, F., Sawyer, W.K., 1980. Transient flow in naturally fractured reservoirs and its application to Devonian gas shales. In: SPE Paper 9397 Presented at the SPE Annual Technical Conference and Exhibition. Dallas, Texas. Lane, H.S., Watson, A.T., Lancaster, D.E., 1989. Identifying and estimating desorption from Devonian shale gas production data. In: SPE Paper 19794 Presented at the SPE Annual Technical Conference and Exhibition. San Antonio, Texas. Ozkan, E., Raghavan, R., Apaydin, O.G., 2010. Modeling of fluid transfer from shale matrix to fracture network. In: SPE Paper 134830 Presented at the SPE Annual Technical Conference and Exhibition. Florence, Italy. 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. Pedrosa, Jr., O.A., 1986. Pressure transient response in stress-sensitive formations. In: Paper SPE 15115 Presented at the 1986 SPE California Regional Meeting, Oakland, Aprial 2–4. Raghavan, R., Chin, L.Y., 2004. Productivity changes in reservoirs with stressdependent permeability. SPE Reser. Eval. Eng. 7 (4), 308–315. Spivey, J.P., Semmelbeck, M.E., 1995. Forecasting long-term gas production of dewatered coal seams and fractured gas shales. In: SPE Paper 29580 Presented at the Low Permeability Reservoirs Symposium. Denver, Colorado. Thomas, V.S., Gringarten, A.C., 2009. Superposition principle and reciprocity for pressure transient analysis of data from interfering wells. SPE J. 4 (3), 488–495. Van Everdingen, A.F., 1953. The skin effect and its influence on the productive capacity of a well. Trans. AIME 198, 171–176. Wang, F.P., Reed, R.M., John, A., et al., 2009. Pore networks and fluid flow in gas shales. In: SPE Paper 124253 Presented at the 2009 SPE Annual Technical Conference and Exhibition. New Orleans, Louisiana. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3 (3), 245–255.