Composite linear flow model for multi-fractured horizontal wells in heterogeneous shale reservoir

Journal of Natural Gas Science and Engineering 38 (2017) 527e548

Contents lists available at ScienceDirect

Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse

Composite linear flow model for multi-fractured horizontal wells in heterogeneous shale reservoir Jie Zeng a, b, Xiangzeng Wang c, Jianchun Guo b, **, Fanhua Zeng a, * a

University of Regina, Petroleum Systems Engineering, Regina, SK S4S 0A2, Canada State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China c Yanchang Petroleum Group, Xi'an, Shaanxi, 710614, PR China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 31 August 2016 Received in revised form 1 January 2017 Accepted 2 January 2017 Available online 9 January 2017

Multi-stage fracturing is currently key technique to develop shale reservoirs. Different analytical models have been proposed to fast investigate post-fracturing pressure- and rate-transient behaviors, and hence, estimate key parameters that affect well performance. However, previous analytical models generally neglect reservoir heterogeneity or typical seepage characters of shale, such as adsorption/desorption, gas slippage, and diffusion effects. This paper presents an analytical model for pressure- and rate-transient analysis of multi-stage fractured shale reservoir, considering heterogeneity, typical seepage characters and, specifically, fluids flow from upper/lower reservoir when vertical fractures partially penetrate the formation. This model is similar to five-region-flow model, but subdivides the reservoir into seven parts, namely, two upper/lower-reservoir regions, two outer-reservoir regions, two inner-reservoir regions, and hydraulic fracture region, which are all transient dual porosity media except the hydraulic fracture. As reservoir heterogeneity along the horizontal wellbore is included, the fracture distribution can be various. Fracture interference is simulated by locating a no-flow boundary between two adjacent fractures. The actual locations of these no-flow boundaries of a specific heterogeneous reservoir are determined based on the pressure value which varies with time and space. Thus, the two sides of this boundary has minimum pressure difference, satisfying the no-flow assumption. Adsorption/desorption, gas slippage and diffusion effects are included for rigorous modeling of flow in shale. This model is validated by comparing with commercial well testing software, obtaining a good match in most flow regimes. Log-log dimensionless pressure, pressure derivative and production type curves are generated to conduct sensitivity analysis. Results suggest that larger desorption coefficient causes smaller pressure and its derivative value as a larger proportion of gas is desorbed in formation and contributes to productivity. Solutions from Azari's (1990, 1991) work, where the effect of fracture height is merely treated as a skin factor are investigated as well. Results show that our model is more accurate in partially penetrating cases, and errors of Azari's method become particularly noticeable in early-middle time response. The influence of other parameters, such as matrix permeability, matrix block size, secondary fracture permeability and hydraulic fracture conductivity, are also discussed. Optimal fracture pattern is selected based on cumulative production. Besides, field data are analyzed and compared graphically with modeling solutions, and reliable results are obtained. As numerical and semianalytical methods require extensive computing processing, this model is a practical alternative to predict well-testing results and select optimal well pattern of shale reservoirs. Reservoir heterogeneities in vertical direction can be further added to our model by vertically subdividing the reservoir into more parts. © 2017 Elsevier B.V. All rights reserved.

Keywords: Shale heterogeneous reservoir Multifractured horizontal well Transient rate and pressure behavior Partially penetrating fracture

* Corresponding author. Petroleum Systems Engineering, Faculty of Engineering and Applied Science, University of Regina, Regina, Saskatchewan, S4S 0A2, Canada. ** Corresponding author. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, 8 Xindu Avenue, Xindu, Chengdu, Sichuan, China. E-mail addresses: [email protected] (J. Guo), [email protected] (F. Zeng). http://dx.doi.org/10.1016/j.jngse.2017.01.005 1875-5100/© 2017 Elsevier B.V. All rights reserved.

528

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

1. Introduction

1.1. Fracture networks in shale reservoir

Shale reservoirs, which have recently gained extensive interests, bring many unique engineering challenges that may not commonly exist in conventional reservoirs. The natural gas within shale reservoir can be stored in the matrix nanopore space, natural fracture systems and adsorbed in the shale organic matter (Schein and Mack, 2007). Currently, multi-stage fractured horizontal wells (MFHWs), which significantly increase the contact area of fracture surface and payzone, are proved to be highly efficient in shale gas/ oil extraction. In published literatures, a number of analytically based solutions have been proposed to fast predict post-frac performance and evaluate critical factors that dominate well behavior since numerical models are time-consuming and need detailed description of fracture and reservoir properties which are usually not available or complete. Brown et al. (2011) proposed a new trilinear flow model with hydraulic fracture region, inner stimulated region (dual porosity) and outer reservoir region to facilitate investigating pressure behavior of MFHWs in shales. They also stated that transient dual-porosity model is more applicable for shale reservoir compared with pseudo dual-porosity idealization. Based on trilinear flow model, Ozkan et al. (2011) made comparison of MFHWs performance in tight sand and shale reservoirs. They concluded that increasing natural fracture density is an effective method to enhance productivity of shale. Later, Stalgorova and Matter (2013) extended Brown et al. (2011) model into a fiveregion model to simulate the fracture with a stimulated surrounding area. Deng et al. (2015) modified five-region model to deal with non-uniform distributed fractures cases, however, the location of impermeable boundary between two adjacent fractures is hard to determine in heterogeneous reservoir. To optimize fracturing treatment size and well spacing, Nobakht and Clarkson (2012) analyzed how the fractures in one well affect the adjacent well production in a zipper-shape fracture pattern. Clarkson et al. (2013) investigated the stress sensitivity of matrix permeability and fracture conductivity and found that they have significant impact on well production, especially for overpressured shale play. Recently, Ozcan et al. (2014), and Wang et al. (2015) combined linear flow with fractal theory and analytically studied the pressure behavior of multifractured unconventional wells, they provided alternative method to simulate well performance in fractured reservoir. In view of shale's multiple flow mechanisms, Zhao et al. (2013), Wang (2013), Liu et al. (2014), Zeng et al. (2015), Zhang et al. (2016), and Zhao et al. (2016) developed analytical and semi-analytical models accounting for desorption, diffusion, and viscous flow. In Wang's (2013) and Liu et al. (2014) models, the stress dependent natural fracture permeability and different angles between fracture and horizontal wellbore are also included. And by adding diffusion and desorption effects, Tian et al. (2014) improved Ezulike and Dehghanpour's (2013) quadrilinear flow model. Ren and Guo (2015) presented a more general semi-analytical model, incorporating various fracture length, diverse fracture intervals, different fracture-wellbore angles, fracture asymmetry and partially penetrating fracture. In summary, these models are versatile enough to model predominant flow regimes for shale but neglect one or several of the following typical characteristics and possible situations in shale reservoir: (1) dual porosity feature caused by natural fractures and artificially created fracture networks, (2) gas slippage in the nanopores, (3) diffusion due to concentration gradient in the matrix pore, (4) desorption from the surface of organic matter, (5) various fracture spacing, (6) partially penetrating fracture, and specifically (7) reservoir heterogeneity.

Natural fractures, including discrete macro fractures and continuous micro fractures (Tian, 2014), are generally observed in shale gas plays which consist of thin layers of easily broken and fine grained rock. The existence of these fractures is a critical factor that enhances well production in such low-permeability formations with nanometer to micrometer pores. Normally, after fracturing operation, some of the pre-existing natural fractures can be reopened, conductive with proppant or intersected by hydraulic fractures, thus, forming complex fracture network which can be detected through micro-seismic mapping (Barree et al., 2002; Mayerhofer et al., 2010). From analytical point of view, dualporosity models, triple-porosity idealizations and even fractional diffusion models are used to simulate the fluid transport in fractured unconventional reservoir. Although, discrete fracture network method is capable to account for detailed properties of individual fracture, it is computationally expensive owing to the requirement of thorough characterization data (Ozcan et al., 2014). Dual porosity models (Barenblatt et al., 1960; Warren and Root, 1963; Kazemi et al., 1969; de Swaan O, 1976) involve applying effective mean-characteristics to describe fractured formation are developed based on uniform fracture and matrix property assumption which is appropriate for continuous-fracture systems. Although, dual porosity models may not always be authentic in actual reservoir condition, they are practical for simulating fractured reservoirs on account of computer's run-time efficiency and minimum data requirement. As an alternative approach, tripleporosity model (Abdassah and Ershaghi, 1986) was proposed to consider two groups of different matrix blocks in addition to the fracture. Later, several investigators, such as Al-Ghamdi and Ershaghi (1996), Liu et al. (2003), Wu et al. (2004), Al-Ghamdi et al. (2010), and Dehghanpour and Shirdel (2011), extended the triple porosity model. By dividing the matrix volume into two subdomains, Ezulike and Dehghanpour (2013) developed quadrilinear flow model to conduct simultaneous matrix to microfracture and matrix to hydraulic fracture depletion. And fractal theory provides an optional way of sophisticated analytical study of fractured reservoirs (Cossio et al., 2012; Ozcan et al., 2014). In fact, dual porosity and triple porosity models are both special cases of multiporosity model (Bai et al., 1993), and the fundamental formulas of our model, similar with Ertekin et al. (1986) and Ozkan et al. (2010) works, are based on de Swaan O's (1976) transient dual porosity assumption combined with linear flow model. 1.2. Multi-mechanism flow of shale gas Gas transfer mechanism in shale reservoirs seems to be significantly different compared with conventional gas formations because of the nanopore structure and adsorption/desorption effect. Therefore, it is crucial to understand how the fluids are stored and transported in the nanopore network before developing the gas-bearing formation. Although, a thorough understanding of gas transfer in shale reservoir has not been obtained, the following commonly accepted mechanisms are summarized, as shown in Fig. 1. 1.2.1. Slip flow, Knudsen diffusion and Darcy flow Generally, the nanopore size of shale organic matrix (average size less than 4 nm-5nm (Kang et al., 2011)) is smaller than the gas molecular free path, thus, apart from gas molecules' collision, the collision between gas molecules and nanopore walls is essential as well and should not be ignored (Wu et al., 2016). In this situation, slip flow (random molecular flow) will occur in nanopores, where the assumption of continuum Darcy flow is no longer valid and

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

529

Fig. 1. Schematic of gas flow in matrix nano pores.

modified Darcy's equation is required. Moreover, as gas is very compressible fluid, a pressure gradient, initiated by well production, can immediately generate a gas concentration gradient (Sigal and Qin, 2008), which drives the diffusion within nanopores (Ertekin et al., 1986). Ordinarily, Knudsen diffusion theory is frequently utilized to model diffusion effects (Javadpour et al., 2007, Javadpour, 2009). And the corresponding diffusion velocity in nanopores of a simplified spherical matrix block can be expressed by the following form governed by Fick's law (Ertekin et al., 1986):

vdm ¼ 

Mg Dm vCm rmg vr

(1)

Where vdm is the velocity of Knudsen diffusion, Mg is the gas molecular weight, Cm is molar concentration, rmg is matrix gas density, and Dm is the diffusivity coefficient and the approximate value can be calculated from Equation (2) proposed by Ertekin et al. (1986) in ft2/D:

Dm

31:54 ¼ pffiffiffiffiffiffiffi k0:67 m Mg

(2)

Where km is the absolute permeability of matrix nanopore. Simultaneously, the pressure forces result in Darcy flow, combined with slip flow, the modified Darcy velocity is written as:

vsm ¼ 

kam vpm

(3)

mmg vr

Where vsm is modified Darcy velocity, kam is the measured matrix gas permeability with slippage effects, mmg is gas viscosity in the matrix and pm is matrix pressure. To take slippage effect into account, Klinkenberg's (1941) equation is applied to describe actual gas permeability:

kam ¼ km 1 þ

bk pm;ave

! (4)

Where bk is gas slippage factor, and pm;ave is average matrix pressure. Klinkenberg (1941) gave the expression of slippage factor. Then, the above equation can also be written in the following form:

kam ¼ km F And

(5)

4clave rp

F ¼1þ

(6)

Where c is a collision proportionality constant which is normally equal to 1 (Zhang et al., 2015), lave is the gas molecular mean free path, and rp is the nanopore radius.

rffiffiffi

p mmg;ave

lave ¼

2 pm;ave

sffiffiffiffiffiffiffi RT Mg

(7)

Where mmg;ave is matrix average gas viscosity, and T is temperature in K. Javadpour (2009) also used the following form of c based on Brown et al. (1946):

c¼

2

at

1

(8)

Javadpour (2009) pointed out that at here is the tangential momentum accommodation coefficient with a value ranging from 0 to 1, determined by the smoothness of nanopore wall. Normally, experimental measurement is needed to determine at for specific shale systems Javadpour (2009). As no experimental data are available to obtain at , here, at is assumed to equal 1, and then, c becomes 1. Therefore, Klinkenberg's (1941) and Brown et al. (1946) slippage factors have the same expression where c equals 1 and either of them can be used to calculate correction factor for slippage. Substituting Equation (7) into Equation (6):

 F ¼1þ

8pRT Mg

0:5

cmmg:ave pm;ave rp

(9)

Note the pressure unit in Equation (9) is in Pa. Next, the total flow velocity through the nanopore in the matrix is assumed to be the combination of parallel diffusion velocity and modified Darcy velocity and is given by:

vtm ¼ 

kam vpm

mmg vr



Mg Dm vCm rmg vr

(10)

The total velocityvtm, as shown in Equation (10), has quite analogous form to that proposed by Ertekin et al. (1986), however, they merely used Fick's law to characterize slip velocity. Similar to conventional formations, Darcy flow is suitable for gas flow in macro fractures (secondary fractures), which will be

530

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

discussed later. For gas transfer from shale matrix to fracture, we follow the derivation of Ertekin et al. (1986) and Ozkan et al. (2010) as shown in Appendix A.

partially penetrating fracture and reservoir heterogeneity.

2. Methodologies 1.2.2. Adsorption/desorption and surface diffusion Due to the large surface area of organic nanoporous materials and their affinity to methane (Yu et al., 2016), a certain part of gas within shale formation is adsorbed in the organic-rich rocks. Once the well produces, reservoir pressure starts decreasing, and the gas molecules on the surface of organic matter desorb to reduce pressure decline. As gas desorption carries on, a gas concentration gradient forms from organic bulk to its surface and controls the surface diffusion of kerogen (Javadpour, 2009). Traditionally, the classic Langmuir isotherm (Langmuir, 1916) has been applied to represent gas single layer desorption for shale formations. Here, the equivalent volume concentration concept is used to describe Langmuir theory, as shown below:

VEsc ¼

VL pm pL þ pm

(11)

Where VEsc is the equilibrium volume at standard condition that can be adsorbed by a rock of unit volume at pressurepm , VL is Langmuir volumetric concentration (constant), and pL is Langmuir pressure (constant). In our model, the effects of surface diffusion are assumed to be included in desorption term since the surface diffusion rate is slow (Javadpour, 2009).

In this paper, de Swaan O's (1976) transient dual-porosity idealization combined with linear flow assumption, as shown in Fig. 2, is applied to model fluids flow within fractured formation. Thus, hydrocarbons within the spherical matrix block would first transfer to the fracture networks (secondary fractures), and then, flow to hydraulic fracture through connected secondary fractures. The entire reservoir is considered to be transient dual-porosity because shale gas reservoir generally has extremely low permeability which conduces long-time transient flow regimes to produce gas (Brown et al., 2011; Kim and Lee. 2015). The spherical matrix element assumption has been made as Kucuk and Sawyer (1980) studied how the choice of matrix element shape (including spherical, cylindrical, Warren and Root and Kazemi models) and size impacts on dual-porosity reservoir behavior, and found that the shape of element makes a little difference during the early-time period and the size of element dominates transient-time behavior. They used cylindrical and spherical element model to analyze Devonian shale gas reservoir and pointed out the limitation of Warren and Root's cubic element model. Besides, many researchers, in the literatures, used spherical matrix assumption to simulate gas flow in tight formation (sands, Devonian shales. coal seams) (Ertekin et al., 1986; Ozkan et al., 2010; Apaydin et al., 2012). Other assumptions that used in this work are listed below:

1.3. Reservoir heterogeneity Reservoir heterogeneity provides uncertainty in both initial and long-term gas production from shale formation. The reservoir properties not only vary from basin to basin, but also they can be quite different even within the same field (Schein and Mack, 2007). Currently, numerical models with fine gridding are primarily chosen to calculate pressure- and rate-transient response of wells in heterogeneous formation. Other approaches, including analytical, semi-analytical and boundary element methods, only deal with simple heterogeneous cases (Medeiros et al., 2010). Thus, for more complicated cases, the limitations of existing analytical and semianalytical models are evident. This paper establishes a fast and practical analytical model for MFHWs in shale reservoir, with consideration of diffusion, slip flow, viscous flow, desorption,

(1) The MFHW is located in the center of a rectangular fractured shale reservoir with impermeable outer boundaries. Reservoir blocks containing different permeability and porosity along the horizontal wellbore are included to present heterogeneity, as shown in Fig. 3. Further investigations can be done to simulate layered reservoirs which are commonly observed due to the nature of deposit. (2) Hydraulic fractures along the wellbore are all vertical and can partially penetrate the formation in both vertical and horizontal directions. Fracture interference is simulated by way of assuming a no-flow boundary between two adjacent hydraulic fractures. For homogeneous reservoir, the no-flow is certainly located in the mid-distance as pressure wave propagate at the same velocity within the formation.

Fig. 2. Schematic of dual-porosity idealization and mass transfer mechanism.

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

hmi ¼

531

2:637  104 kmappi;ave ði fmi aave cmg;ave mmg;ave

¼ 1; 2; 3; 4; 5 or 6 represents flow regionsÞ

Fig. 3. Schematic of MFHWs in a heterogeneous box-shape reservoir.

However, pressure wave propagates faster and pressure drops quicker in high-permeability region, the no-flow boundary for heterogeneous reservoir may not be at the mid-distance. These no-flow boundaries are determined based on reservoir pressure, in other words, during the well life, the two sides of the boundary has minimum pressure difference, and hence, the fluids flow cross the boundary can be ignored. (3) Fluids flow from reservoir to hydraulic fracture is simulated by seven linear flow regions (Fig. 4a) which depict reservoir blocks with identical or different properties. Regions 3 and 4 are needed to simulate the non-penetrated reservoir (regions beyond the fracture tip) while regions 1 to 4 are required to simulate reservoir in the left side of no-flow boundary as shown in Fig. 4b. Because this model is developed to deal with multi-fracture horizontal wells, the geometry of these systems is not likely to lead to radial flow between two adjacent fractures (Stalgorova and Matter, 2013). Gas-phase viscous flow, combined with diffusion, slippage and desorption effects, occurs within the matrix, while in hydraulic fracture and secondary fracture, the gas flow is merely subject to Darcy's law.

2.1. Dimensionless variables For convenience, the solutions of this model are founded on dimensionless variables. The dimensionless pseudopressure and time are defined respectively as:

ppD ¼

tD ¼

i kref href h pp ðpini Þ  pp ðpÞ : 1422qFt T

href t d2ref

(12)

(13)

href ¼

2:637  104 kref fref ctref mref

hfi ¼

2:637  104 kfi ffi cfg;ave mfg;ave

(16)

hF ¼

2:637  104 kF fF cFg;ave mFg;ave

(17)

Where ctref , fref , andmref represent reference values of compressibility, porosity and viscosity respectively. kmappi;ave is apparent matrix permeability of ith region at average pressure (if average pressure is available), or following Ozkan et al. (2010) and Brown et al. (2011) who used initial condition to obtain linear flow equation, although this may result in calculation error especially when the pressure drop is significant. The subscripts m, f , and F denote matrix, secondary fracture and hydraulic fracture properties.

  mmg;ave cmg;ave Dm kmapp;ave ¼ kma;ave 158 þ1 kam;ave

(14)

(18)

kfi andkF are secondary fracture permeability of ith region and hydraulic fracture permeability. 158 is the value of unit conversion factor. The unit conversion process can be found in Appendix A. And cmg;ave , cfg;ave , and cFg;ave are gas compressibility in matrix, secondary fracture and hydraulic fracture at average pressure, therefore, their values are equal. Similarly, mfg;ave and mFg;ave are average gas viscosity in secondary fracture and hydraulic fracture. And fmi , ffi and fF are the corresponding porosity of matrix, secondary fracture in ith region and hydraulic fracture. Zhang et al. (2016) put forward a parameter called adsorption factor to describe the effects of adsorption/desorption. The bigger adsorption coefficient value is, the larger amount of gas are absorbed in the rock. Here, a dimensionless parameter, namely desorption coefficient, which is the reciprocal of the above mentioned adsorption factor, has been proposed to describe the effects of gas desorption.

aave ¼

mÞ cmg;ave þ ð1f f m

psc TZave Zsx Tsc pm;ave

cmg;ave

VL pL

ðpL þpm;ave Þ

2

(19)

The larger desorption coefficeint is, the more gases are desorbed from the organic matter during production. The dimensionless forms of above diffusivity are shown by:

.

hmDi ¼ hmi href .

Where dref , href and kref are the reference length, diffusivity and permeability respectively. By separately defining constant reference parameters rather than using reservoir properties as reference value, one can perform sensitivity analysis in dimensionless form, as changing reservoir properties will not affect analysis results. qFt is the total production rate of MFHW. And T is temperature in  R. The reference diffusivity, matrix diffusivity in ith region, secondary fracture diffusivity in ith region, and hydraulic fracture diffusivity are given by:

(15)

hfDi ¼ hfi href

(20)

.

hFD ¼ hF href According to Al-Hussainy and Ramey (1966), the pseudopressure term in Equation (12) is defined by:

Zp pp ðpÞ ¼ 2 pref

p

dp

mg ðpÞZðpÞ

(21)

In the dual porosity model, the flow from the matrix block is instantly and evenly distributed in one-half the fracture volume (de Swaan O, 1976) and the fracture heighthf is small compared with matrix radius rm , the half fracture volume is:

532

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

Fig. 4. Schematici of the seven-flow-region model: (a) fluids flow in seven regions; (b) possible flow region distribution in a heterogeneous multi-fractured reservoir.

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

Vf 2 hf ¼ 4prm 2 2

(22)

Then, the flow-capacity and storativity ratios of ith region are defined as (Serra et al., 1983; Ozkan et al., 2010):

li ¼ sd2ref

ui ¼

kmappi;ave Vm k rm 2 . ¼ sd2ref mappi;ave kfi kfi hf Vf 2 3

fmi aave cmg;ave Vm 2 fmi aave cmg;ave rm ¼ 3 ffi cfg;ave Vf ffi cfg;ave hf

(23)

Wherehf ands are secondary fracture height and shape factor. Chang (1993) derived the shape factor for spherical matrix flow as:

15 s¼ 2 rm

(29)

2.2.2. Fracture-reservoir system 2.2.2.1. Region 6. The diffusivity equation for 1D flow in z-direction within the upper reservoir region is:

vz2D

 cðsÞ6 sppfD6 ¼ 0

(30)

Where

cðsÞ6 ¼

1

hfD6



l6 5s

" 1

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15su6

hfD6 l6

coth

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 15su6

hfD6 l6

(31)

(25) The boundary condition at the top of reservoir is:

vppfD6 ¼0 vzD zD ¼z2D

The dimensionless distances are:

x1D ¼ xf

ppmD ðrD ¼ 1; xD ; sÞ ¼ ppfD ðxD ; sÞ

v2 ppfD6 (24)

533

.

Pressure continuity condition on the interfaces of regions 2 and 6, and regions 4 and 6 can be expressed by:

dref

. xeD ¼ xe dref . y1D ¼ y1 dref . y2D ¼ y2 dref . .  2dref z1D ¼ z1 dref ¼ hF .  . z2D ¼ z2 dref ¼ h 2dref . wD ¼ wF dref rD ¼ r=rm

(32)

ppfD6

zD ¼z1D

¼ ppfD2

zD ¼z1D

¼ ppfD4

zD ¼z1D

(33)

(26) 2.2.2.2. Region 5. Similarly, the diffusivity equation for region 5 is:

v2 ppfD5 vz2D

 cðsÞ5 sppfD5 ¼ 0

(34)

Where

cðsÞ5 ¼

2.2. Method and algorithm An eighth of the fracture-reservoir system is considered due to the symmetry, as shown in Fig. 4, to model flow mechanism for shale reservoir with a MFHW. Seven subsystems are presented including two upper/lower flow regions (regions 5 and 6), two outer flow regions (regions 3 and 4), two inner flow regions (regions 1 and 2) and hydraulic fracture region. In each flow region, only 1D flow occurs in secondary fracture, and the flow directions in these region are proved to be reasonable as we have compared our solution with commercial software's calculations. Flow equations with boundary conditions are provided below in Laplace domain. The detailed derivation of matrix spherical flow and fracture linear flow solutions are given respectively in Appendix A and Appendix B.

ppfD5

zD ¼z1D

And at the matrix spherical surface:

1

hfD5 l5

coth

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 15su5

hfD5 l5

(35)

(36)

¼ ppfD1

zD ¼z1D

¼ ppfD3

zD ¼z1D

(37)

2.2.2.3. Region 4. In this region, 1D flow in x-direction is applied and the diffusivity equation is

Where

vx2D

cðsÞ4 ¼ (28)

5s

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15su5

On the interfaces of regions 5 and 1, and regions 5 and 3, the boundary condition is:

  vppmD 1 v s 2 r ¼ p ði ¼ 1; 2; 3; /; 6Þ D 2 vr vrD hmDi pmD rD D

ppmD ðrD ¼ 0; xD ; sÞ ¼ ppmD ð0; xD ; sÞ



"

vppfD5 ¼0 vzD zD ¼z2D

v2 ppfD4

The inner boundary condition at rD ¼ 0 is given by:

hfD5

l5

At the reservoir top, the boundary condition is:

2.2.1. Matrix spherical flow Diffusivity equation of a spherical matrix block in ith region can be written as:

(27)

1

kf 6;v vppfD6 þ kf4 ;v z1D vzD

1

hfD4



l4 5s

" 1

 cðsÞ4 sppfD4 ¼ 0

(38)

zD ¼z1D

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15su4

hfD4 l4

coth

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 15su4

hfD4 l4

(39)

The subscript v represents properties in vertical direction. Closed boundary condition at outer reservoir is given by:

534

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

vppfD4 ¼0 vxD xD ¼xeD

(40)

vy2D

Pressure continuity condition on the interfaces of regions 4 and 2 is:

ppfD4

xD ¼x1D

v2 ppfD1

¼ ppfD2

(41)

xD ¼x1D

v2 ppfD3 vx2D

kf 5;v vppfD5 þ kf 3;v z1D vzD

 cðsÞ3 sppfD3 ¼ 0

(50)

Where

cðsÞ3 ¼

1

hfD3



l3 5s

" 1

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15su3

hfD3 l3

coth

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 15su3

(43)

hfD3 l3

1

hfD1



l1 5s

" 1

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15su1

hfD1 l1

coth

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 15su1

hfD1 l1

(51)

Based on the flux-continuity at the interface of regions 2 and 1 and pressure-continuity conditions between hydraulic fracture and region 1, the boundary conditions are:

kf 1

Where

zD ¼z1D

¼0

(42)

zD ¼z1D

xD ¼x1D

kf 5;v vppfD5 þ kf 1;v z1D vzD

 cðsÞ1 sppfD1

cðsÞ1 ¼ 2.2.2.4. Region 3. In a similar way, diffusivity equation for region 3 is:

kf 3 vppfD3 þ kf 1 x1D vxD

vppfD1 vppfD2 ¼ k f2 vyD yD ¼y1D vyD yD ¼y1D

ppfD1

yD ¼

wD 2

¼ ppFD

wD 2

yD ¼

(52)

(53)

No-flow condition at outer boundary:

vppfD3 ¼0 vxD xD ¼xeD

(44)

2.2.2.7. Fracture region. Diffusivity equation for hydraulic fracture is given by:

Pressure continuity on the interface of regions 3 and 1:

ppfD3

xD ¼x1D

¼ ppfD1

(45)

xD ¼x1D

2.2.2.5. Region 2. As this region is adjacent to regions 4 and 6, the diffusivity equation becomes:

v2 ppfD2 vy2D

þ

kf 4 vppfD4 kf 2 x1D vxD

þ xD ¼x1D

kf 6;v vppfD6 kf 2;v z1D vzD

zD ¼z1D

 cðsÞ2 sppfD2 ¼0

(46)

1

hfD2



l2 5s

" 1

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 15su2

hfD2 l2

coth

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi # 15su2

hfD2 l2

(47)

The no-flow boundary approximation between two adjacent fracture results in:

vppfD2 ¼0 vyD yD ¼y2D

(48)

And the pressure continuity on the interface of regions 1 and 2 is given by:

ppfD2

vx2D

þ

vppfD1 kf 1b s  p ¼0 kF ðwD =2Þ vyD yD ¼wD =2 hFD pFD

(54)

Where kf 1b is fracture bulk permeability in region 1. Strictly speaking, bulk permeability concept should be applied in every dual porosity section. However, because the fracture-to-reservoir volume ratio is the same in all dual porosity regions, it is equivalent to directly use intrinsic properties to describe boundary conditions of two dual porosity regions. Fluids flow from hydraulic fracture tip are neglected:

vppFD ¼0 vxD xD ¼x1D

(55)

Applying Darcy's law at fracture and wellbore interface, the following equation is obtained:

Where

cðsÞ2 ¼

v2 ppFD

yD ¼y1D

¼ ppfD1

yD ¼y1D

vppFD q pkref href dref ¼  Fi vxD xD ¼0 qFt kF wF hF

(56)

Where qFi is the flow rate of ith fracture in Laplace domain. Solving the above diffusivity equations with their boundary conditions, the solution of constant rate production problem is given by:

ppFD

xD ¼0

¼

pkref href dref qFi

pffiffiffiffiffiffi  pffiffiffiffiffiffi qFt kF wF hF aF tanh aF x1D

(57)

In terms of constant pressure case, the definition of dimensionless pressure is different:

(49) ppD

  ppini  pp  ¼ ppini  ppwf

(58)

ref

2.2.2.6. Region 1. Diffusivity equation for the region adjacent to hydraulic fracture is:

In Laplace domain, the dimensionless pressure under constant pressure condition at fracture and wellbore interface is:

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

ppFD

xD ¼0

¼

1 s

L½pF2 ðqF2 Þ  pF1 ðqF1 Þ ¼ 0 /  i h L pFn ðqFn Þ  pFðn1Þ qFðn1Þ ¼ 0

(59)

The relationship between dimensionless rate of ith fracture under constant pressure condition and dimensionless pressure under constant rate condition for the same well is then given by:

qDi ¼

qFi 1 qFt p pFDi

1 s

n X

(60)

qFi ¼

xD ¼0

qDi

n  X

(65)

   i h L pF1;1 qF1;1  pF1;2 qF1;2 ¼ 0 h  /  i L pFn;1 qFn;1  pFn;2 qFn;2 ¼ 0

(66)

(62) (3) Equations (63), (65) and (66) can be described by the following matrix:

By applying Stehfest numerical inversion method (Stehfest,

a11 6 0 6 6 6 6 0 6 6 a11 6 6 0 6 6 6 4 0 1

 qFi;1 þ qFi;2 ¼ 2qFt

Pressure drops caused by the production from two sides of hydraulic fracture are equal as well:

Cumulative production in Laplace form is consequently obtained:

2

(64)

i¼1

(61)

qD s

qFi;1 þ qFi;2 2

Where qFi;1 =2 and qFi;2 =2 represent flow rate from left-side drainage area and right-side drainage area of hydraulic fracture. Fig. 5 demonstrates rate distribution of a MFHW. Therefore:

i¼1

QD ¼

(63)

Due to the reservoir heterogeneity, the actual production contribution from two fracture sides can be different, therefore:

It is important to notice that the meaning of ppFD in Equations (59) and (60) is totally different. The term ppFD , in Equation (60), is obtained from Equation (57). Thus Equation (60) is used to repre sent the relationship between ppFDi (constant rate) and qDi xD ¼0 (constant flowing bottomhole pressure) from the same well. And Equation (59) shows the value of constant dimensionless flowing bottomhole pressure in Laplace domain. The total dimensionless rate for constant pressure condition is the sum of each fracture's dimensionless rate:

qD ¼

535

a12 0

0 a21

0 a22

0 a12 0

0 a21 a21

0 a22 a22

0 1

0 1

0 1

0 0 / 0 0 a31 / 0 1

0 0

/ /

0 0

0 0

0 0

0 0 a32

/ / /

0 0 0

0 0 0

an1 0 0

0 1

0 1

an1;1 1

an1;2 1

an;1 1

1970), dimensionless pressure, rate and cumulative production can be transformed into real-time domain.

2

3 3 2 qF1;1 0 7 6 qF1;2 7 6 0 7 7 6 7 7 6 7 6 qF2;1 7 6 / 7 7 6 7 7 6 6 7 7 6 an2 7 7 6 qF2;2 7 6 0 7 6 6 7 7 0 7  6 qF3;1 7 ¼ 6 0 7 7 6 7 7 6 0 7 7 6 qF3;2 7 6 0 7 7 6 / 7 6 / 7 7 6 7 7 6 an;2 5 4 qFn;1 5 4 0 5 qFn;2 L½2qFt  1 0 0

3

(67)

Where

!

pkref href dref

pffiffiffiffiffiffi pffiffiffiffiffiffi aij ¼ qFt kF wF hF tanh aF x1D aF 2.2.3. Coupling solutions Because pressure drop within the horizontal wellbore is negligible for gas flow, the assumption of infinite conductivity wellbore has been made. Thereupon, at every fracture-wellbore interfaces, the same pressure is maintained. Then, follow the frame work presented below, the problem can be solved. (1) In heterogeneous reservoir constant-rate condition, the flow rates of fractures in different reservoir blocks can be diverse. But initially, each fracture has equal flow rate qFi ¼ qFt =n. n is fracture number. (2) As each hydraulic fracture has an initial value of production rate, the corresponding pressure at each fracture-wellbore interface is obtainable. Under infinite conductivity wellbore condition, no pressure drop occurs:

(68) ij

(4) After solving the linear system presented in Equation (67), new flow rates for each fracture are obtainqF1 0 ; qF2 0 ; / qFi 0 ; /qFn 0 .Where

qFi ¼

qFi;1 þ qFi;2 2

(69)

Then Stehfest algorithm is utilized to obtain real-time solutions qF1 0 ; qF2 0 ; /; qFi 0 ; /qFn 0 . (5) Comparing qF1 0 ; qF2 0 ; /; qFi 0 ; /; qFn 0 with the initial flow rates qF1 ; qF2 ; /; qFi ; /; qFn .

536

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

(6) If the absolute value of minimum difference between qFi andqFi 0 is small enough (in this case 1011), then, qF1 0 ; qF2 0 ; /; qFi 0 ; /; qFn 0 are considered as exact fractures' flow rates. Otherwise, qF1 0 ; qF2 0 ; /; qFi 0 ; /; qFn 0 will be used as new initial input parameters and steps (1) to (5) are repeated until the absolute value of minimum difference is within the range of acceptable value.

2.2.4. No-flow boundary determination The no-flow boundary, as mentioned before, is not a real physical boundary but rather is an interference boundary between two adjacent fractures. To determine the exact location of no-flow boundary, the pressure value of region 2 at yD ¼ y2D should be calculated first from the following equations:

ppfD2

y2D

2ppfD1 y1D pffiffiffiffiffiffi

pffiffiffiffiffiffi

¼ exp a2 ðy2D  y1D Þ þ exp  a2 ðy2D  y1D Þ (70)

And

ppfD1

y1D

3. Model validation Previously, Ozkan et al. (2010) have verified the dual porosity equations for spherical matrix blocks by comparing with conventional slab-matrix transient dual porosity model. Even so, partially penetrating fracture case, where upper and lower reservoir regions’ contributions are included, has not been validated yet. In this section, to systematically verify our model, the commercial well testing software KAPPA Ecrin v4.30 is applied. The analytical test design model with fractured horizontal well, two porosity spherical matrix reservoir and rectangle boundary has been selected to conduct this comparison. As Olarewaju and Lee (1989) and Azari et al. (1990, 1991) also considered and examed the effects of different fracture height on pressure response for linear flow model before, using the dimensionless fracture height concept that resembles a skin factor, their methods are included in our verification process. The desorption, slippage and diffusion effects are not taken into account in this analytical test design model. For a fair comparison, results without diffusion, slippage and desorption effects are compared with those from commercial software and previous methods. It is worth noting that the definitions of storativity ratio and interporosity flow parameter in the software, as shown in

pffiffiffiffi pffiffiffiffi

pffiffiffiffiffiffi 



pffiffiffiffiffiffi  pffiffiffiffi pffiffiffiffiffiffi a k k a tanh½pa pffiffiffiffi1 f 1 f 2 pffiffiffiffi2 ffiffiffiffi2 ðy2D y1D Þ  w exp  a1 y1D þ exp  2 a1 y1D exp a1 y1D a1 kf 1 þkf 2 a2 tanh½ a2 ðy2D y1D Þ  ¼ ppFD D    pffiffiffiffiffiffi  p ffiffiffiffi p ffiffiffiffi p ffiffiffiffi

 pffiffiffiffiffiffi pffiffiffiffiffiffi a1 kf 1 kf 2 pa 2 ffiffiffiffi2 tanh½paffiffiffiffi2 ðy2D y1D Þ  exp  a1 w2D þ exp  2 a1 y1D exp a1 w2D pffiffiffiffi a k þk a tanh½ a ðy y Þ  1 f1

In our model, the value of ppFD at yD ¼ wD =2 in Equation (71) cannot be directly calculated because hydraulic fracture pressure is only function of xD and time, not yD . Average hydraulic fracture, defined in Equation (72), could be a kind of method to describe average pressure behavior within hydraulic fracture.

Z ppFD

ave

¼

2

2

2D

1D

Equations (74) and (75), are quite different from our model, thus, adjustment are required and the equivalent values of them have been used in commercial software.

u¼

Vf ff ctf Vm fm ctm þ Vf ff ctf

(74)

x1D 0

ppFD dxD x1D

(72)

Where

ppFD

f2

(71)

2 l ¼ arw

(75)

Where rw is wellbore radius, kmb and kfb are bulk matrix perme-

pffiffiffiffiffiffi

  pffiffiffiffiffiffi qFi pkref href dref exp aF ðx1D  xD Þ þ exp  aF ðx1D  xD Þ 

pffiffiffiffiffiffi 

pffiffiffiffiffiffi ¼ pffiffiffiffiffiffi qFt kF wF hF aF exp aF x1D  exp  aF x1D

However, the calculation of average pressure, involving integration, requires cumbersome computational efforts. Comparing pressure at wellbore and fracture interface with average pressure, it is easy to find out that conspicuous difference occurs in early time behavior (tD ¼ 1e10, t ¼ 0.883hre8.831hr) while in late time the difference becomes much smaller, as shown in Fig. 6. Therefore, in this work, directly using ppFD at xD ¼ 0 to approximately represent pressure within hydraulic fracture. Then, we obtain the pressure values at both sides of no-flow boundary, and compare the percentage of difference between them. Finally, this boundary is located in a place where there is minimum cumulative pressure difference during production.

kmb kfb

(73)

ability and fracture permeability. The input parameters are listed in Figs. 7 and 8. From the comparison results, the early-time behavior shows an infinite fracture pressure response with 1/2 slope as the hydraulic fracture conductivity is large. This regime is followed by a transient regime, and then, a bilinear flow regime with 1/4 slope and a very short period 1/2 formation linear flow regime. Next, after a transient region, the unit slope boundary dominated regime finally appears. It is also shown that the agreement of the three methods is excellent under fully penetrating condition. However, if the fracture penetrating ratio is only 70%, our prediction matches better with commercial software's outcome than that from Olarewaju and lee (1989) and Azari et al. (1990, 1991) methods during early-middle time response. This can be caused by the fact

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

537

Fig. 5. Schematic of fracture flow rate distribution.

Fig. 6. Comparison of average hydraulic fracture pressure and pressure at fracture and wellbore interface.

Fig. 8. Comparison of our model, commercial software, and previous methods (70% penetrating ratio).

that the geometry of linear flow changes on the interface of fracture and inner formation in partially penetrating cases. Pressure behaviors at the very beginning (t < 0.01 h) are not available in

commercial software. Overall, our model and its algorithm are trustable. However, the model also has limitations, errors may occur when the dimensionless hydraulic fracture conductivity is

Fig. 7. Comparison of our model, commercial software and previous methods (fully penetrating fracture).

Fig. 9. Effects of gas desorption on pressure behavior.

538

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

small. 4. Results and discussions 4.1. Effects of gas desorption Gas desorption plays a crucial role in shale gas production. As mentioned earlier, desorption coefficient, defined by Equation (19), is applied here to analyze the effects of gas desorption. As illustrated in Fig. 9, the larger desorption coefficient is, the lower pressure drop will be because the desorbed gas provides an extra source for well production. It also puts off the arrival of boundary dominated regime. However, fracture dominated linear flow period is not affected by gas desorption as only fluids within the hydraulic fracture contribute to well production at that time. Figs. 10 and 11 show that the contribution of desorption is essential. With the amount of desorbed gas increases, the dimensionless rate declines slower, resulting in longer production time and significant increment in cumulative production.

Fig. 12. Effects of matrix radius on pressure behavior.

4.2. Effects of matrix radius Variation of spherical matrix radius mainly influences reservoir

Fig. 13. Effects of matrix radius on rate behavior.

Fig. 10. Effects of gas desorption on rate behavior.

Fig. 14. Effects of matrix radius on cumulative production.

Fig. 11. Effects of gas desorption on cumulative production.

bulk properties. According to Apaydin et al. (2012), the bulk or effective properties can be inferred from well performance data and are also dependent on the interpretation model. They gave a certain range of matrix radius, from 3 to 6 ft. For secondary fracture, the

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

thickness or height was assumed to be 5  104ft to 5  103ft. Kucuk and Sawyer (1980) studied how the choice of matrix element shape impacts on dual-porosity behavior, and found that the size of element dominates transient-time behavior. In our model, the variation of matrix radius mainly influences the secondary fracture to reservoir bulk volume ratio, and hence, affects fracture bulk properties, flow-capacity and storativity ratios. Pressure and rate behaviors under different matrix radius are studies and compared. Results indicate that matrix radius mainly controls transient flow region. The greater the radius is, the later hydraulic fracture linear flow will vanish and the lager pressure drop will occur in transient regime (see Fig.12). Since the increase in matrix radius leads to smaller fracture-reservoir ratio and fracture bulk permeability, the production rate is smaller under larger matrix radius with a slower declining pace (see Fig. 13). It seems that the choice of matrix radius has a critical impact on cumulative production, as shown in Fig. 14. But in practical, the intrinsic and bulk properties are normally first determined, and then, the corresponding matrix radius is obtained. Under given intrinsic and bulk properties, the simulation results are reliable.

539

Fig. 16. Effects of matrix permeability on rate behavior.

4.3. Effects of matrix permeability In reality, shale reservoir typically has ultra-low permeability ranging from 10 to 100 nanodarcies (Cipolla et al., 2010). To highlight the influence of matrix permeability, a wider range of matrix permeability has been analyzed in this section. Fig. 15 demonstrates how matrix permeability affects pressure response. Smaller matrix permeability results in heavier pressure drop during transient time, longer fracture linear flow period and later arrival of boundary response. As for rate behavior, dimensionless rate increases with the growth of matrix permeability with a weakening trend, but drops quicker. For large matrix permeability cases, dimensionless rate declines at similar pace (see Fig. 16). This is accordance with pressure response where there are marginal differences on the appearance of boundary dominated flow under high matrix permeability. The shape of cumulative production curve is less subject to matrix permeability, when it reaches 104md (see Fig. 17). One can conclude that cumulative production is more sensitive to the variation of matrix permeability with ultra-low value.

Fig. 17. Effects of matrix permeability on cumulative production.

4.4. Effects of secondary fracture permeability Because virtually all gases are transferred to hydraulic fracture

Fig. 18. Effects of secondary fracture permeability on pressure behavior.

Fig. 15. Effects of matrix permeability on pressure behavior.

through secondary fracture (fracture systems) with relative high conductivity, the permeability of fracture system could be a parameter of more significance than the matrix permeability. Though, the intrinsic secondary fracture permeability seems to be

540

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

large and ranges from 2md to 2000md, the bulk value of fracture permeability is in the range of 6  104md to 0.6md. Experimental results show that the secondary fracture values used here are reasonable (Ning et al., 1993; Cho et al., 2013; Zhang et al., 2014). Besides, several researchers used comparative secondary/natural fracture permeability value (Ozkan et al., 2010; Brown et al., 2011; Apaydin et al., 2012; Ozcan et al., 2014). The results presented in Fig. 18 reveal that secondary fracture permeability has stronger influence on pressure transient behavior than matrix permeability. When the intrinsic secondary fracture permeability is small, four flow regimes including fracture linear flow, transient flow, formation linear flow and boundary-dominate flow (not shown in Fig. 18). And the hump in pressure derivative after 1/2-slpoe fracture linear flow becomes more prominent. With secondary fracture permeability increase, the formation linear flow regime gradually disappears. The corresponding flow regimes can also be observed in dimensionless rate curve. Compared with matrix permeability, the secondary fracture permeability affects rate behavior in a similar manner. The most noticeable difference is that the increase in fracture permeability accelerates rate decline even the fracture permeability is relatively high (see Fig.19). Results in Fig. 20 indicate that cumulative production decreases dramatically when fracture permeability changes from 2000md to 2md. The secondary fracture permeability has no effects on ultimate cumulative production but

Fig. 19. Effects of secondary fracture permeability on rate behavior.

Fig. 20. Effects of secondary fracture permeability on cumulative production.

ultra-low permeability cases take an unrealistic long period to reach ultimate cumulative production. Therefore, effectively creating secondary fractures or enhancing existing natural fracture permeability could be one of the key factors to ensure a successful fracturing treatment in shale. 4.5. Effects of hydraulic fracture conductivity Figs. 21e23 present pressure, rate and cumulative production data under different hydraulic fracture conductivities. Different from single porosity tight reservoir where fracture conductivity influence early time pressure behavior (Yao et al., 2013), fracture conductivity in our dual porosity model dominates the whole production period except boundary flow regime. However, with fracture conductivity increase, the dimensionless rate and cumulative production growth tend to slow down. Therefore, for a specific reservoir, there exists an optimal fracture conductivity. Further incremental in fracture conductivity is unnecessary. 4.6. Fracture pattern optimization example for heterogeneous reservoir The objective of this section is to provide an application example

Fig. 21. Effects of hydraulic fracture conductivity on pressure behavior.

Fig. 22. Effects of hydraulic fracture conductivity on rate behavior.

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

541

Fig. 23. Effects of hydraulic fracture conductivity on cumulative production.

Fig. 25. Cumulative production of case 1 to case 5 in around 10 years.

of this model to demonstrate the effectiveness of this model in dealing with reservoir heterogeneity. In this part, the exact locations of no-flow boundaries between adjacent factures for each case, calculated from this model, are shown for the first time in this paper, which helps the readers have a better understanding of our model. Meanwhile, some fracture pattern optimization strategies are obtained. Since multi-stage hydraulic fracturing is a costly reservoir stimulation technique, the fracture pattern optimization is an indispensable part in treatment design to achieve economical rates and reduce financial risk. Reservoir heterogeneity is one of the most significant factors that affect the design of fracture patterns. The example utilized as a template framework in optimization of heterogeneous reservoir and its properties are illustrated in Fig. 24. It is a reservoir composed of five blocks with equal size but different permeability. A 10-year cumulative production is applied to select optimal fracture pattern. Compared with matrix permeability, it is shown that for heterogeneous reservoir, the overall performance

are mainly controlled by secondary fracture permeability. Thus, in this section, only heterogeneity caused by different secondary fracture permeability is discussed. Although numerous cases can be run using our model, the following five cases are enough to provide general ideas in fracture pattern design. From Figs. 24 and 25, it is easy to find out that in heterogeneous reservoir, the location of fracture plays a more significant role in effectively stimulating the formation than fracture number. In other words, inappropriately placed fracture is not only ineffective but cost additional initial investment. Comparing cases 2, 3 and 4, it is clear that case 4 with 3 fracture performs as well as the case with 4 fractures. It is particularly noticeable that case 3 with 3 fractures even obtains better cumulative production than the 4-fractrue case. Under the same fracture number, case 3 has more than twice the cumulative production compared with case 5. Therefore, cases 3 and 1 could be the optimal designs here. Further analysis is required to select the most economical one. It is recommended that

Fig. 24. Schematic of reservoir model and the calculated locations of no-flow boundary.

542

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

Fig. 26. Rate and pressure data of a MFHW in Barnett shale (Apiwathanasorn, 2012).

diffusion, slippage and desorption on interpretation results are anslyzed through production history matching, using the same basic input data as Apiwathanasorn (2012) who used pressure and derivative in history matching. Four cases are considered, namely, model without diffusion, slippage and desorption (case 1), model with slippage and diffusion (case 2), model with desorption (case 3), and model with all three typical mechanisms (case 4). The interpreted values of fracture bulk permeability for best matching, as shown in Fig. 27, range from 1.45nd to 3.25nd, which are reasonable compared with Apiwathanasorn's (2012) results. Of course, transient rate analysis normally has more than one solutions, results from our model and Apiwathanasorn's model are in the same order of magnitude. In fact, the above bulk fracture permeabilities are calculated from the originally estimated intrinsic fracture permeabilities. Since the matrix radius and secondary fracture thickness in our model are assumed to be 5 ft and 0.001 ft respectively, the fracture to reservoir bulk volume ratio is fixed and can be used to obtain bulk properties. The intrinsic properties for four cases are listed below: C Case 1: secondary fracture permeability: 0.013md, matrix permeability: 1.5  107md, fracture bulk permeability: 3.25nd. C Case 2: secondary fracture permeability: 0.012md, matrix permeability: 9  108md, fracture bulk permeability: 3nd. C Case 3: secondary fracture permeability: 0.0062md, matrix permeability: 4  107md, fracture bulk permeability: 1.55nd. C Case 4: secondary fracture permeability: 0.0058md, matrix permeability: 1.8  107md, fracture bulk permeability: 1.45nd.

Fig. 27. Type curve matching of field data.

priority should be put on reservoir blocks with high permeability when selecting candidate blocks for fracturing treatment in a heterogeneous reservoir. For vertical heterogeneity within upper/lower reservoir region, the reservoir can be subdivided into more parts. And in outer/inner reservoir region, layer thickness weighted average properties can be used to approximately simulate the entire system.

5. Field example A MFHW, with 39 stages, completed in Barnett shale is applied here to perform rate transient analysis (Apiwathanasorn, 2012). The actual rate and corresponding pressure data covered a period of about 600 days are given in Fig. 26. All input parameters used in our model are list in Fig. 27. For convenience, each stage is treated as a hydraulic fracture. From Fig. 26, it is shown that the gas pressure overall tends to be stable (approaching constant pressure condition), while at early time, there is relatively larger pressure fluctuation. Therefore, as the corresponding pressure for a given rate is difficult to obtain, the production data are directly utilized in history matching instead of utilizing pressure normalized rate. Previously, Apiwathanasorn (2012) has investigated the well performance using different models to conduct history matching with good agreement. Specifically, for dual-porosity model matching, three possible flow regime were used in that study. Results suggest that the bulk fracture permeability (working permeability) is estimated at 6.5nd to 8.4nd. In this work, the effects of

From field case analysis, although the exact reservoir properties remain unknown, the following conclusions can be obtained from our data: neglecting diffusion and slippage effects can overestimate both matrix and secondary fracture permeability. If ignoring desorption effect, secondary fracture permeability would be overestimated while the matrix permeability tends to be underestimated. However, compared with diffusion and slippage, desorption seems to have heavier influence on reservoir property evaluation. 6. Conclusions (1) This paper presents an analytical for model multi-fractured horizontal wells in heterogeneous shale reservoir. It is validated by comparing with commercial software. By comparing with Azari et al. (1990, 1991) work, where the effect of fracture height is merely treated as a skin factor, it is shown that our model provides more accurate results in partially penetrating fracture cases. (2) Sensitivity analysis results suggest that larger desorption coefficient results in smaller pressure and its derivative value because a larger proportion of gas is desorbed in formation and contributes to productivity. Neglecting gas desorption or inaccurately estimate desorption related parameters can cause significant error in investigating the resource potential. Matrix permeability affects cumulative production only within a certain range. If the matrix permeability value is larger than this range, no obvious difference appears in cumulative production with matrix permeability change. And matrix element size controls fracture-reservoir volume ratio which determines reservoir bulk properties. Selecting appropriate matrix and fracture dimensions plays a key role

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

in obtaining reliable simulation results. Compared with matrix permeability, cumulative production is more sensitive to secondary fracture permeability scince it is hard or impossible to change the matrix permeability. Therefore, reopening natural fractures or artificially creating fracture networks, which will not only increase secondary fracture permeability but also increase fracture-reservoir volume ratio, is highly recommended. There is an optimal fracture conductivity for each specific case. Insufficient hydraulic fracture conductivity cannot make fully used of well production potential. However, producing under excessive fracture conductivity results in additional investment and is not economical. (3) The location of no-flow boundaries in heterogeneous reservoir may not be located in the mid-distance of two adjacent fractures. In fracture pattern optimization, fractures are better placed in high-permeability blocks to receive shorttime return in heterogeneous reservoir. (4) Through field data matching demonstrate the reliability of our model. It is also shown that the existence of diffusion, gas slippage, and desorption brings more uncertainties and complexity in field data analysis. Knowing how these mechanisms influence evaluation results can reduce the uncertainty and obtain data with higher accuracy. Acknowledgements The authors would like to acknowledge the financial support from National Science Fund for Distinguished Young Scholars (Grant No. 51525404) and National Natural Science Foundation of China (Grant No. 51374178).

xf xe x1D xeD y y1 y2 yD y1D

543

fracture half-length, ft reservoir size in x direction, ft dimensionless fracture half-length dimensionless reservoir size in x direction y-coordinate, ft distance from fracture to permeability boundary, ft distance to no-flow boundary between two fractures, ft dimensionless y-coordinate dimensionless distance from fracture to permeability boundary dimensionless distance to no-flow boundary between two fractures z-coordinate, ft half-fracture height, ft half reservoir thickness, ft dimensionless z-coordinate dimensionless half-fracture height dimensionless half reservoir thickness gas compressibility factor fluid density, lbm/ft3 fluid viscosity, cp diffusivity, ft/hour dimensionless fracture diffusivity porosity, fraction shape factor, fraction flow capacity ratio, fraction gas molecular mean free path, ft intermediate variable for this model tangential momentum accommodation coefficient, fraction

y2D z z1 z2 zD z1D z2D Z

r m h hD

f

s l lave a at

Appendix A. Derivation process and solutions of matrix block Nomenclature aij aave Bg bk c cg Cm Dm dref h k Lw Mg N p ppD pp q qD QD r s t tD T v V wF wD x xD

element of coefficient matrix desorption coefficient at average pressure gas volumetric factor, fraction slippage factor, psi collision proportionality constant gas compressibility, 1/psi molar concentration, lbm-mol/ft3 diffusivity coefficient, ft2/D reference length, ft height or thickness, ft permeability, md well length, ft gas molecular weight, lbm/lbm-mol ftacture number or stages pressure, psi dimensionless pressure pseudopressure, psi2/cp production rate, Mscf/D dimensionless rate dimensionless cumulative production radius, ft Laplace etransform parameter time, hours dimensionless time reservoir temperature,  R fluid velocity, ft/D volume, ft3 hydraulic fracture width, ft dimensionless hydraulic fracture width x-coordinate, ft dimensionless x-coordinate

Here, the derivation process of multi-mechanism dual-porosity model follows the similar line as Ertekin et al. (1986) and Ozkan et al. (2010). Starting with a mass balance equation for spherical matrix block where the uniform distribution assumption of flux and pressure is made on matrix surface:

 i  h  v fm rmg v ð1  fm Þrgsc VEsc 1 v 2 r rmg vtm ¼ þ  2 vt vt r vr

(A-1)

The second term on the right side of Equation A-1 represents the mass rate of desorption. And the total matrix transport velocity vtm consists of two velocity terms, one for modified Darcy flow with gas slippage, the other for diffusion effect:

vtm ¼ vsm þ vdm

(A-2)

The first term on the right side of Equation A-2 is modified Darcy velocity which is given by:

vsm ¼ 

kam vpm

mmg vr

(A-3)

The diffusion velocity formula is defined as:

Mg Dm vCm vdm ¼  rmg vr

(A-4)

Then, total velocity formula in matrix is the combination of Equations A-3 and A-4, and can be rewritten as:

Mg Dm vCm kam vpm  vtm ¼  rmg vr mmg vr In real gas case, the gas density can be expressed as:

(A-5)

544

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

rmg ¼

side:

pm Mg ZRT

(A-6)

For single phase flow of real gas, the molar concentration is:

Cm

pm ¼ ZRT

(A-7)



2



psi psi

kam ðmdÞ vpm ðpsiÞ      ft mmg ðcpÞ psi*D vrðftÞ cp ft

vtm

vðpm =ZÞ pm 1 1 vZ ¼  vr Z pm Z vpm



vpm  vr

1 kam vpm ¼ 158 mmg vr (A-18)

vpm pm cmg vpm ¼ vr Z vr

Therefore, the total matrix velocity is given by Equation A-11:

(A-9)

vpm kam vpm  vr mmg vr

vtm ¼ 

kam

mmg

Where cmg is gas compressibility in matrix. Because the fracture and matrix porosities are assumed to be constant, the value of gas compressibility also equals total compressibility. Finally, the following form of total velocity can be obtained:

vtm ¼ cmg Dm

1

6894:76*103 mPa*86400s 1mPa*s

(A-8)

Further extending the term vðpm =ZÞ:



¼



304800*304800mm2 0:987*103 mm2

mmg 

Substituting Equations A-6 and A-7 into Equation A-4 yields:

Mg Dm vðpm =ZÞ kam vpm  ¼ vr rmg RT mmg vr

1

kam 





ft md

(A-10)

 158

mmg cmg Dm kam

 kmapp vpm vpm ¼ þ1 vr mmg vr

(A-19)

And hence, Equation A-15 is obtained. Since both kmapp and ctmapp are pressure dependent variables, to linearize our matrix flow equations, the relative parameters are assumed to be calculated from average pressure (if available) and are defined by:

  mmg;ave cmg;ave Dm kmapp zkmapp;ave ¼ kma;ave 158 þ1 kam;ave

(A-20)

Equation A-10 can also be written as:

vtm ¼ 

kam

mmg

 158

mmg cmg Dm kam

 þ1

vpm vr

(A-11)

ctmapp zctmapp;ave ¼ cmg;ave þ

Making the substitution of Equations A-11 in Equation A-1:

" #   mmg cmg Dm 1 v 2 kam vpm 158 r þ 1 r mg mmg kam vr r 2 vr ¼ fm

vrmg vV þ rsc ð1  fm Þ Esc vt vt

pL þ pm;ave

2

Converting Equation A-14 into dimensionless form:

(A-12)

  1 v 1 vppmD 2 vppmD r ¼ D 2 vr vrD hmD vtD rD D

(A-22)

Linearized dimensionless matrix diffusivity and secondary fracture diffusivity are given below:

vVEsc VL pL vpm ¼ vt ðpL þ pm Þ2 vt

(A-13)

Applying the pseudopressure defined in Equation (21), Equation A-12 becomes:

  vppm vppm 1 v 2 r ¼ fm mmg ctmapp k mapp vr vt r 2 vr

  mmg cmg Dm kmapp ¼ kma 158 þ1 kam ð1  fm Þ psc TZ VL pL fm Zsx Tsc pm ðpL þ pm Þ2

fref mref ctref kmapp;ave d2ref 2 fm mmg;ave ctmapp;ave kref rm

¼

fref mref ctref kmapp;ave d2ref 2 fm mmg;ave acmg;ave kref rm

(A-23)

hfD ¼

ff mfg;ave ctf ;ave kref fref mref ctref kf

(A-24)

The initial condition is given by:

(A-15) ppmD ðrD ; xD ; tD ¼ 0Þ ¼ 0 (A-16)

(A-25)

The dimensionless pseudopressure at inner boundary (rD ¼ 0) is a finite value, that is:

ppmD ðrD ¼ 0; xD ; tD Þ ¼ ppmD ð0; xD ; tD Þ

Equation A-15 is obtained as follows: Recalling Equation A-10:



hmD ¼

(A-14)

Where

vtm ðft=DÞ ¼ cmg ð1=psiÞDm

VL pL

(A-21)

Where

ctmapp ¼ cmg þ

ð1  fm Þ psc TZave fm Zsx Tsc pm;ave 

(A-26)

And at the outer boundary:

.  vp ðpsiÞ m ft D vrðftÞ 2

kam ðmdÞ vpm ðpsiÞ  mmg ðcpÞ vrðftÞ

ppmD ðrD ¼ 1; xD ; tD Þ ¼ ppfD ð1; xD ; tD Þ (A-17)

It is clear that the unit of the first term in the right side of Equation (16) is ft/D. And for the second term of the right-hand

(A-27)

Where ppfD ðxD ; tD Þ is the average pressure at the matrix and secondary fracture interference. To solve Equations A-22 to A-27, substituting the following expression into Equations A-22, A-25 to A-27:

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

FmD ðrD ; xD ; tD Þ ¼ rD ppmD ðrD ; xD ; tD Þ

(A-28)

And applying Laplace transformation:

v F mD s  F ¼0 2 hmD mD vrD

(A-29)

F mD ðrD ¼ 0; xD ; sÞ ¼ 0

(A-30)

F mD ðrD ¼ 1; xD ; sÞ ¼ ppfD ðxD ; sÞ

(A-31)

v2 ppfD6

Solving Equations A-29 to A-30, yields:

pffiffiffiffiffiffiffiffiffiffiffiffiffi  sinh s=hmD rD ppmD ðrD ; xD ; sÞ ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffi ppfD ðxD ; sÞ s=hmD rD sinh

(A-32)

¼

15u

(A-33)

lhfD

Therefore, Equation A-32 can be rewritten as:

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi .   sinh 15su lhfD rD ppmD ðrD ; xD ; sÞ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi .  ppfD ðxD ; sÞ rD sinh 15su lhfD

k vp rg mapp m mmg vr

!

(A-34)

ppfD6

(A-35)

vx

k vp rmg mapp m mmg vr

2 þ 4 hf

!

3 5¼

  v rfg ff

ðrm ;x;tÞ

v ppfD vx2D



2kmapp;ave d2ref hf kf rm

vppmD v rD

¼ ð1;xD ;tD Þ

1 vppfD

hfD vtD

vt

(A-36)

(A-37)

vx2D





l vppmD

5

vrD ð1;xD ;sÞ

¼

s

hfD

ppfD

(A-38)

2

vx2D

 cðsÞsppfD ¼ 0

(A-39)

Where

cðsÞ ¼

1

hfD



l 5s

" 1

sffiffiffiffiffiffiffiffiffiffiffi 15su

hfD l

coth

zD ¼z1D

(B-3)

qffiffiffiffiffiffiffiffiffiffiffiffiffi hqffiffiffiffiffiffiffiffiffiffiffiffiffi i vppfD6 ¼ p cðsÞ6 stanh cðsÞ6 s ðz2D  z1D Þ pfD2 vzD z1D z1D

(B-4)

qffiffiffiffiffiffiffiffiffiffiffiffiffi hqffiffiffiffiffiffiffiffiffiffiffiffiffi i vppfD6 ¼ p cðsÞ6 stanh cðsÞ6 s ðz2D  z1D Þ pfD4 vzD z1D z1D

(B-5)

v2 ppfD5

 cðsÞ5 sppfD5 ¼ 0

(B-6)

vppfD5 ¼0 vzD zD ¼z2D ppfD5

zD ¼z1D

¼ ppfD1

(B-7)

zD ¼z1D

¼ ppfD3

zD ¼z1D

(B-8)

qffiffiffiffiffiffiffiffiffiffiffiffiffi hqffiffiffiffiffiffiffiffiffiffiffiffiffi i vppfD5 ¼ p cðsÞ5 stanh cðsÞ5 s ðz2D  z1D Þ pfD3 vzD z1D z1D

(B-9)

qffiffiffiffiffiffiffiffiffiffiffiffiffi hqffiffiffiffiffiffiffiffiffiffiffiffiffi i vppfD5 ¼ p cðsÞ5 stanh cðsÞ5 s ðz2D  z1D Þ pfD1 vzD z1D z1D In region 4, the diffusivity equation and boundary conditions are

Submit the derivative of Equation A-34 at fracture and matrix interface into Equation A-38:

v ppfD

zD ¼z1D

¼ ppfD4

(B-10)

Transformed into Laplace form:

v2 ppfD

¼ ppfD2

(B-2)

Solving Equations B-6 to B-8:

The dimensionless form of Equation A-36 is: 2

zD ¼z1D

vz2D

ðr¼rm ;x;tÞ

2

(B-1)

vppfD6 ¼0 vzD zD ¼z2D

Then, mass balance equation for fracture system is:

  v rfg vf

 cðsÞ6 sppfD6 ¼ 0

Similarly, for region 5, the diffusivity equation and boundary conditions are:

The matrix inflow mass rate per unit fracture volume per unit time can be expressed as:

2 f ðx; tÞ ¼  hf

vz2D

Solving Equations B-1 to B-3:

The relationship among dimensionless matrix and fracture diffusivities, storativity ratio and flow-capacity ratio can be summarized as:

1

Appendix B. Derivation process and solutions of seven-flowregion model (constant rate condition) The diffusivity equation and boundary conditions for region 6 are:

2

hmD

545

sffiffiffiffiffiffiffiffiffiffiffi # 15su

hfD l

v2 ppfD4 vx2D

þ

kf 6;v vppfD6 kf4 ;v z1D vzD

 cðsÞ4 sppfD4 ¼ 0

vppfD4 ¼0 vxD xD ¼xeD ppfD4

xD ¼x1D

¼ ppfD2

(B-11)

zD ¼z1D

(B-12)

xD ¼x1D

(B-13)

Solving Equations B-11 to B-8:

(A-40)

Equation A-39 is the diffusivity equation for secondary fracture within each flow region.

pffiffiffiffiffiffi vppfD4 pffiffiffiffiffiffi ¼ p a4 tanh½ a4 ðxeD  x1D Þ pfD2 vxD x1D x1D Where

(B-14)

546

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

a4 ¼

hqffiffiffiffiffiffiffiffiffiffiffiffiffi i kf 6;v qffiffiffiffiffiffiffiffiffiffiffiffiffi cðsÞ6 stanh cðsÞ6 s ðz2D  z1D Þ þ cðsÞ4 s kf 4;v z1D

a2 ¼ (B-15)

Similarly, for region 3, the diffusivity equation and boundary conditions are:

v2 ppfD3 vx2D

þ

kf 5;v vppfD5 kf 3;v z1D vzD

 cðsÞ3 sppfD3 ¼ 0 zD ¼z1D

xD ¼x1D

¼ ppfD1

(B-25) In region 1, the diffusivity equation and boundary conditions are:

v2 ppfD1

vppfD3 ¼0 vxD xD ¼xeD ppfD3

(B-16)

kf 4 pffiffiffiffiffiffi pffiffiffiffiffiffi a4 tanh½ a4 ðxeD  x1D Þ kf 2 x1D hqffiffiffiffiffiffiffiffiffiffiffiffiffi i kf 6;v qffiffiffiffiffiffiffiffiffiffiffiffiffi cðsÞ6 stanh þ cðsÞ6 s ðz2D  z1D Þ þ cðsÞ2 s kf 2;v z1D

(B-17)

vy2D

þ

kf 3 vppfD3 kf 1 x1D vxD

þ xD ¼x1D

kf 5;v vppfD5 kf 1;v z1D vzD

zD ¼z1D

 cðsÞ1 sppfD1 ¼0

(B-18)

xD ¼x1D

kf 1

Solving Equations B-16 to B-18:

pffiffiffiffiffiffi vppfD3 pffiffiffiffiffiffi ¼ p a3 tanh½ a3 ðxeD  x1D Þ pfD1 vxD x1D x1D

(B-19)

(B-26)

vppfD1 vppfD2 ¼ kf 2 vyD yD ¼y1D vyD yD ¼y1D

ppfD1

w yD ¼ 2D

¼ ppFD

wD 2

yD ¼

(B-27)

(B-28)

Solving Equations B-26 to B-28:

Where

vppfD1 ¼ ppFD b1 vyD wD =2 wD =2

qffiffiffiffiffiffiffiffiffiffiffiffiffi hqffiffiffiffiffiffiffiffiffiffiffiffiffi i k cðsÞ5 stanh a3 ¼ f 5;v cðsÞ5 s ðz2D  z1D Þ þ cðsÞ3 s kf 3;v z1D (B-20)

(B-29)

Where

In region 2, the diffusivity equation and boundary conditions are

 pffiffiffiffiffiffi pffiffiffiffiffiffi exp  a1 b1 ¼ a1  pffiffiffiffiffiffi exp  a1

v2 ppfD2 vy2D

þ

kf 4 vppfD4 kf 2 x1D vxD



pffiffiffiffiffiffi

 pffiffiffiffiffiffi  exp  2 a1 y1D exp a1  pffiffiffiffiffiffi

 pffiffiffiffiffiffi wD þ exp  2 a1 y1D exp a1 2 wD 2

þ xD ¼x1D

kf 6;v vppfD6 kf 2;v z1D vzD

wD 2 wD 2

 

pffiffiffiffi pffiffiffiffi pffiffiffiffi a k k a tanh½pa pffiffiffiffi1 f 1 f 2 pffiffiffiffi2 ffiffiffiffi2 ðy2D y1D Þ a1 kf 1 þkf 2 a2 tanh½ a2 ðy2D y1D Þ pffiffiffiffi pffiffiffiffi pffiffiffiffi a k k a tanh½pa pffiffiffiffi1 f 1 f 2 pffiffiffiffi2 ffiffiffiffi2 ðy2D y1D Þ a1 kf 1 þkf 2 a2 tanh½ a2 ðy2D y1D Þ

a1 ¼ zD ¼z1D

 cðsÞ2 sppfD2 ¼0

(B-31) For fracture region, the diffusivity equation and boundary conditions are

(B-22) v2 ppFD vx2D

yD ¼y1D

¼ ppfD1

yD ¼y1D

(B-23)

Solving Equations B-21 to B-23:

pffiffiffiffiffiffi vppfD2 pffiffiffiffiffiffi ¼ p a2 tanh½ a2 ðy2D  y1D Þ pfD1 vyD y1D y1D Where

kf 3 pffiffiffiffiffiffi pffiffiffiffiffiffi a3 tanh½ a3 ðxeD  x1D Þ kf 1 x1D hqffiffiffiffiffiffiffiffiffiffiffiffiffi i kf 5;v qffiffiffiffiffiffiffiffiffiffiffiffiffi cðsÞ5 stanh þ cðsÞ5 s ðz2D  z1D Þ þ cðsÞ1 s kf 1;v z1D

(B-21)

vppfD2 ¼0 vyD yD ¼y2D ppfD2

(B-30)

(B-24)

þ

vppfD1 kf 1b s  p ¼0 kF ðwD =2Þ vyD yD ¼wD =2 hFD pFD

(B-32)

vppFD ¼0 vxD xD ¼x1D

(B-33)

vppFD q pkref href dref ¼  Fi vxD xD ¼0 qFt kF wF hF

(B-34)

Finally, solving Equations B-32 to B-34:

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

ppFD

xD ¼0

pkref href dref q

pffiffiffiffiffiffi  ¼ Fi pffiffiffiffiffiffi qFt kF wF hF aF tanh aF x1D

(B-35)

Where

aF ¼

kf 1b s b þ kF wD =2 1 hFD

(B-36)

References Abdassah, D., Ershaghi, I., 1986. Triple-porosity systems for representing naturally fractured reservoirs. SPE Form. Eval. 1 (02), 113e127. SPE-13409-PA. http://dx. doi.org/10.2118/13409-PA. Al-Ghamdi, A., Ershaghi, I., 1996. Pressure transient analysis of dually fractured reservoirs. SPE J. 1 (01), 93e100. SPE-26959-PA. http://dx.doi.org/10.2118/ 26959-PA. Al-Ghamdi, A., Chen, B., Qanbari, F., et al., 2010. An improved triple porosity model for evaluation of naturally fractured reservoirs. Present. at Trinidad Tobago Energy Resour. Conf. Port Spain, Trinidad 27e30. June. SPE-132879-MS. http:// dx.doi.org/10.2118/132879-MS. Al-Hussainy, R., Ramey Jr., J.H., 1966. Application of real gas flow theory to well testing and deliverability forecasting. J. Petroleum Technol. 18 (05), 637e642. SPE-1243-B-PA. http://dx.doi.org/10.2118/1243-B-PA. Apaydin, O.G., Ozkan, E., Raghavan, R., 2012. Effect of discontinuous microfractures on ultratight matrix permeability of a dual-porosity medium. SPE Reserv. Eval. Eng. 15 (04), 473e485. SPE-147391-PA. http://dx.doi.org/10.2118/147391-PA. Apiwathanasorn, S., 2012. Evidence of Reopen Microfractures in Production Data of Hydraulically Fractured Shale Gas Wells. MS thesis. Texas A&M University, College Station, Texas (August 2012). Azari, M., Wooden, W.O., Coble, L.E., 1990. A complete set of Laplace transforms for finite-conductivity vertical fractures under bilinear and trilinear flows. Present. at SPE Annu. Tech. Conf. Exhib. New Orleans, La. 23e26. September. SPE-20556MS. http://dx.doi.org/10.2118/20556-MS. Azari, M., Wooden, W.O., Coble, L.E., 1991. Further investigation on the analytic solutions for finite-conductivity vertical fractures. Present. at Middle East Oil Show. Bahrain 16e19. November. SPE-21402-MS. http://dx.doi.org/10.2118/ 21402-MS. Bai, M., Elsworth, D., Roegiers, J.C., 1993. Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs. Water Resour. Res. 29 (6), 1621e1634. Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, PMM. Sov. Appl. Math. Mech. 24 (5), 852e864. Barree, R.D., Fisher, M.K., Woodroof, R.A., 2002. A Practical Guide to Hydraulic Fracture Diagnostic Technologies. Presented at SPE Annual Technical Conference and Exhibition. Texas, San Antonio, 29 September-2 October. SPE-77442MS. http://dx.doi.org/10.2118/77442-MS. Brown, G.P., Dinardo, A., Cheng, G.K., et al., 1946. The flow of gases in pipes at low pressures. J. Appl. Phys. 17, 803e813. Brown, M., Ozkan, E., Raghavan, R., et al., 2011. Practical solutions for pressuretransient responses of fractured horizontal wells in unconventional shale reservoirs. SPE Reserv. Eval. Eng. 14 (06), 663e676. SPE-125043-PA. http://dx.doi. org/10.2118/125043-PA. Chang, M.M., 1993. Deriving the Shape Factor of a Fractured Rock Matrix. Technique Report NIPER-696, (DE93000170). National Institute for Petroleum and Energy Research, Bartlesville, Oklahoma, USA. Cho, Y., Apaydin, O.G., Ozkan, E., 2013. Pressure-dependent natural-fracture permeability in shale and its effect on shale-gas well production. SPE Reserv. Eval. Eng. 16 (02), 216e228. SPE-159801-PA. http://dx.doi.org/10.2118/159801PA. Cipolla, C.L., Lolon, E.P., Erdle, J.C., et al., 2010. Reservoir modeling in shale-gas reservoirs. SPE Reserv. Eval. Eng. 13 (04), 638e653. SPE-125530-PA. http://dx. doi.org/10.2118/125530-PA. Clarkson, C.R., Qanbari, F., Nobakht, M., et al., 2013. Incorporating geomechanical and dynamic hydraulic-fracture-property changes into rate-transient analysis: example from the haynesville shale. SPE Reserv. Eval. Eng. 16 (03), 303e316. SPE-162526-PA. http://dx.doi.org/10.2118/162526-PA. Cossio, M., Moridis, G.J., Blasingame, T.A., 2012. A Semi-analytic Solution for Flow in Finite-conductivity Vertical Fractures Using Fractal Theory. Presented at SPE Latin America and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, pp. 16e18. April. SPE-153715-MS. http://dx.doi.org/10.2118/153715-MS. de Swaan O, A., 1976. Analytic solutions for determining naturally fractured reservoir properties by well testing. SPE J. 16 (03), 117e122. SPE-5346-PA. http://dx. doi.org/10.2118/5346-PA. Dehghanpour, H., Shirdel, M., 2011. A Triple Porosity Model for Shale Gas Reservoirs. Presented at Canadian Unconventional Resources Conference, Calgary, Alberta, Canada, pp. 15e17. November. SPE-149501-MS. http://dx.doi.org/10.2118/ 149501-MS.

547

Deng, Q., Nie, R., Jia, Y., et al., 2015. 2015. A new analytical model for non-uniformly distributed multi-fractured system in shale gas reservoirs. J. Nat. Gas Sci. Eng. 27, 719e737. http://dx.doi.org/10.1016/j.jngse.2015.09.015. Ertekin, T., King, G.A., Schwerer, F.C., 1986. Dynamic gas slippage: a unique dualmechanism approach to the flow of gas in tight formations. SPE Form. Eval. 1 (01), 43e52. SPE-12045-PA. http://dx.doi.org/10.2118/12045-PA. Ezulike, D.O., Dehghanpour, H., 2013. A model for simultaneous matrix depletion into natural and hydraulic fracture networks. J. Nat. Gas Sci. Eng. 16 (2014), 57e69. http://dx.doi.org/10.1016/j.jngse.2013.11.004. Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J. Can. Petroleum Technol. 48 (08), 16e21. PETSOC-09-0816-DA. http://dx.doi.org/10.2118/09-08-16-DA. Javadpour, F., Fisher, D., Unsworth, M., 2007. Nanoscale gas flow in shale gas sediments. J. Can. Petroleum Technol. 46 (10), 55e61. PETSOC-07-10-06. http://dx. doi.org/10.2118/07-10-06. Kang, S.M., Fathi, E., Ambrose, R.J., et al., 2011. Carbon dioxide storage capacity of organic-rich shales. SPE J. 16 (04), 842e855. SPE-134583-PA. http://dx.doi.org/ 10.2118/134583-PA. Kazemi, H., Seth, M.S., Thomas, G.W., 1969. The interpretation of interference tests in naturally fractured reservoirs with uniform fracture distribution. SPE J. 9 (04), 463e472. SPE 2156-B. http://dx.doi.org/10.2118/2156-B. Kim, T.H., Lee, K.S., 2015. Pressure-transient characteristics of hydraulically fractured horizontal wells in shale-gas reservoirs with natural- and rejuvenatedfracture networks. J. Can. Petroleum Technol. 54 (04), 245e258. SPE-176027PA. http://dx.doi.org/10.2118/176027-PA. Klinkenberg, L.J., 1941. The Permeability of Porous Media to Liquids and Gases. APIDrilling and Production Practice, New York, New York. Kucuk, F., Sawyer, W.K., 1980. Transient flow. In: Naturally Fractured Reservoirs and its Application to Devonian Gas Shales. Presented at SPE Annual Technical Conference and Exhibition, Dallas, Texas, pp. 21e24. September. SPE-9397-MS. http://dx.doi.org/10.2118/9397-MS. Langmuir, L., 1916. The constitution of fundamental properties of solids and liquids. J. Am. Chem. Soc. 38, 2221e2295. Liu, J., Bodvarsson, G.S., Wu, Y.-S., 2003. Analysis of flow behavior in fractured lithophysal reservoirs. J. Contam. Hydrology 62e63 (2003), 189e211. http:// dx.doi.org/10.1016/S0169-7722(02) 00169e9. Liu, M., Xiao, C., Wang, Y., et al., 2014. Sensitivity analysis of geometry for multistage fractured horizontal wells with consideration offinite-conductivity fractures in shale gas reservoirs. J. Nat. Gas Sci. Eng. 22 (2015), 182e195. http://dx. doi.org/10.1016/j.jngse.2014.11.027. Mayerhofer, M.J., Lolon, E., Warpinski, R., et al., 2010. What is stimulated reservoir volume? SPE Prod. Operations 25 (01), 89e98. SPE-119890-PA. http://dx.doi. org/10.2118/119890-PA. Medeiros, F., Ozkan, E., Kazemi, H., 2010. A semianalytical approach to model pressure transients in heterogeneous reservoirs. SPE Reserv. Eval. Eng. 13 (02), 341e358. SPE-102834-PA. http://dx.doi.org/10.2118/102834-PA. Ning, X., Fan, J., Holditch, S.A., et al., 1993. Property Measurement in Naturally Fractured Devonian Shale Cores Using a New Pressure Pulse Method. SCA Conference Paper. 9301. Nobakht, M., Clarkson, C.R., 2012. Multiwell Analysis of Multifractured Horizontal Wells in Tight/Shale Gas Reservoirs. Presented at SPE Canadian Unconventional Resources Conference. Calgary, Alberta, Canada, 30 October-11 November. SPE162831-MS. http://dx.doi.org/10.2118/162831-MS. Olarewaju, J.S., Lee, W.J., 1989. A New Analytical Model of Finite-conductivity Hydraulic Fracture in a Finite Reservoir. Presented at SPE Gas Technology Symposium. Texas, Dallas, pp. 7e9. June. SPE-19093-MS. http://dx.doi.org/10.2118/ 19093-MS. Ozcan, O., Sarak, H., Ozkan, E., et al., 2014. A Trilinear FlowModel for a Fractured Horizontal Well in a Fractal Unconventional Reservoir. Presented at SPE Annual Technical Conference and Exhibition. Amsterdam, The Netherlands, 27e29 October. SPE-170971-MS. http://dx.doi.org/10.2118/170971-MS. Ozkan, E., Raghavan, R., Apaydin, O.G., 2010. Modeling of Fluid Transfer from Shale Matrix to Fracture Network. Presented at SPE Annual Technical Conference and Exhibition. Italy, Florence, pp. 19e22. September. SPE-134830-MS. http://dx.doi. org/10.2118/134830-MS. Ozkan, E., Brown, M.L., Raghavan, R., et al., 2011. Comparison of fracturedhorizontal-well performance in tight sand and shale reservoirs. SPE Reserv. Eval. Eng. 14 (02), 248e259. SPE-121290-PA. http://dx.doi.org/10.2118/121290PA. Ren, J., Guo, P., 2015. A novel semi-analytical model forfinite-conductivity multiple fractured horizontal wells in shale gas reservoirs. J. Nat. Gas Sci. Eng. 24 (2015), 35e51. http://dx.doi.org/10.1016/j.jngse.2015.03.015. Schein, G.W., Mack, D.J., 2007. In: Economides, M.J., Martin, T., Huston, T.X. (Eds.), Chapter 11 Unconventional Gas. Modern Fracturing. Serra, K., Reynolds, A.C., Raghavan, R., 1983. New pressure transient analysis methods for naturally fractured reservoirs (includes associated papers 12940 and 13014 ). J. Petroleum Technol. 35 (12), 2271e2283. SPE-10780-PA. http://dx. doi.org/10.2118/10780-PA. Sigal, R.F., Qin, B., 2008. Examination of the importance of self diffusion in the transportation of gas in shale gas reservoirs. Petrophysics 49 (3), 301e305. Stalgorova, K., Matter, L., 2013. Analytical model for unconventional multifractured composite systems. SPE Reserv. Eval. Eng. 16 (03), 246e256. SPE-162516-PA. http://dx.doi.org/10.2118/162516-PA. Stehfest, H., 1970. Remark on algorithm 368: numerical inversion of Laplace transforms. Commun. ACM 13 (10), 624. http://dx.doi.org/10.1145/

548

J. Zeng et al. / Journal of Natural Gas Science and Engineering 38 (2017) 527e548

355598.362787. Tian, Y., 2014. Experimental Study on Stress Sensitivity of Naturally Fractured Reservoir. Presented at SPE Annual Technical Conference and Exhibition, Amsterdam, The Netherlands, pp. 27e29. October. SPE-173463-STU. http://dx. doi.org/10.2118/173463-STU. Tian, L., Guo, X., Liu, M., et al., 2014. Well testing model for multi-fractured horizontal well for shale gas reservoirs with consideration of dual diffusion in matrix. J. Nat. Gas Sci. Eng. 21 (2014), 283e295. http://dx.doi.org/10.1016/j. jngse.2014.08.001. Wang, H., 2013. Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms. J. Hydrology 510 (2014), 299e312. http://dx.doi.org/10.1016/j.jhydrol.2013.12.019. Wang, W., Su, Y., Sheng, G., et al., 2015. A mathematical model considering complex fractures and fractal flow for pressure transient analysis of fractured horizontal wells in unconventional reservoirs. J. Nat. Gas Sci. Eng. 23, 139e147. http://dx. doi.org/10.1016/j.jngse.2014.12.011. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3 (03), 245e255. SPE-426-PA. http://dx.doi.org/10.2118/426-PA. Wu, Y.-S., Liu, H.H., Bodvarsson, G.S., 2004. A triple-continuum approach for modeling flow. J. Contam. Hydrology 73 (2004), 145e179. http://dx.doi.org/ 10.1016/j.jconhyd.2004.01.002. Wu, K., Li, X., Guo, C., et al., 2016. A unified model for gas transfer in nanopores of shale-gas reservoirs: coupling pore diffusion and surface diffusion. SPE J. Prepr. 21 (05), 1583e1611. SPE-2014-1921039-PA. http://dx.doi.org/10.2118/20141921039-PA.

Yao, S., Zeng, F., Liu, H., et al., 2013. A semi-analytical model for multi-stage fractured horizontal wells. J. Hydrology 507 (2013), 201e212. http://dx.doi.org/10. 1016/j.jhydrol.2013.10.033. Yu, W., Sepehrnoori, K., Patzek, T.W., 2016. Modeling gas adsorption in marcellus shale with Langmuir and BET isotherms. SPE J. 21 (02), 589e600. SPE-170801PA. http://dx.doi.org/10.2118/170801-PA. Zeng, H., Fan, D., Yao, J., et al., 2015. Pressure and rate transient analysis of composite shale gas reservoir considering multiple mechanisms. J. Nat. Gas Sci. Eng. 27 (2015), 914e925. http://dx.doi.org/10.1016/j.jngse.2015.09.039. Zhang, J., Kamenov, A., Zhu, D., et al., 2014. Laboratory measurement of hydraulicfracture conductivityies in the Barnett shale. SPE Prod. Operations 29 (03), 216e227. SPE-163839-PA. http://dx.doi.org/10.2118/163839-PA. Zhang, P., Hu, L., Meegoda, J.N., et al., 2015. Micro/Nano-pore network analysis of gas flow in shale matrix. Scientific reports 5. Artic. number 13501 (2015). http:// dx.doi.org/10.1038/srep13501. Zhang, L., Gao, J., Hu, S., et al., 2016. Five-regionflow model for MFHWs in dual porous shale gas reservoirs. J. Nat. Gas Sci. Eng. 33 (2016), 1316e1323. http://dx. doi.org/10.1016/j.jngse.2016.05.063. Zhao, Y., Zhang, L., Zhao, J., et al., 2013. Triple Porosity”Modeling of transient well test and rate decline. J. Petroleum Sci. Eng. 110 (2013), 253e262. http://dx.doi. org/10.1016/j.petrol.2013.09.006. Zhao, Y., Zhang, L., Xiong, Y., et al., 2016. Pressure response and production performance for multi-fractured horizontal wells with complex seepage mechanism in box-shaped shale gas reservoir. J. Nat. Gas Sci. Eng. 32 (2016), 66e80. http://dx.doi.org/10.1016/j.jngse.2016.04.037.