A new approach to model shale gas production behavior by considering coupled multiple flow mechanisms for multiple fractured horizontal well

Fuel 237 (2019) 283–297

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

A new approach to model shale gas production behavior by considering coupled multiple flow mechanisms for multiple fractured horizontal well Ting Lua,b,c,d, Shimin Liuc, , Zhiping Lia,b,d, ⁎

T

⁎

a

School of Energy Resource, China University of Geosciences (Beijing), Beijing 100083, China Beijing Key Laboratory of Unconventional Natural Gas Geological Evaluation and Development Engineering, Beijing 100083, China c Department of Energy and Mineral Engineering, G3 Center and Energy Institute, The Pennsylvania State University, University Park, PA 16802, USA d Key Laboratory of Strategy Evaluation for Shale Gas, Ministry of Land and Resources, Beijing 100083, China b

ARTICLE INFO

ABSTRACT

Keywords: Coupled flow mechanisms Production analysis Sorption Diffusion Stress-dependent permeability

Production decline analysis of multiple fractured horizontal wells (MFHW) is crucial for long-term shale gas development. Analytical solutions of production decline model accounting for sorption, diffusion (slip flow and Knudsen diffusion) and non-static (stress-dependent) permeability can precisely predict the long term production behavior, especially for the late time production period. However, little work has simultaneously incorporated all these mechanisms into production decline analysis for shale gas wells. In this work, a new production decline model for MFHW in shale gas reservoirs incorporating multiple flow mechanisms is established. To weaken the strong nonlinearities of seepage mathematical equation caused by combining multiple mechanisms, perturbation technology is employed to establish the point source solution considering stress-dependent permeability of MFHW and little effort has been done on this before. Besides, Laplace transformation, numerical discrete method, Stehfest numerical inversion algorithm and Gaussian elimination method are employed to solve the new model’s mathematical problem. Estimated inversion values of reservoir parameters are consistent with the reported data from Barnett shale. Further, comparisons between production behaviors with and without multimechanics flows were made in three flow periods and the production discrepancy increases with continuous depletion, which is attributed that desorption and diffusion is increasingly important as pressure depleting. Finally, the effects of major parameters on production decline curves are analyzed by using the proposed model and it was found that different parameters have their own influence period and sensitivity intensity.

1. Introduction The conventional natural gas productions are depleting at an alarming rate and the gas production from shale becomes one of the primary energy sources for United State. The shale gas revolution is only possible because of the technology advancements including multistage fracturing and horizontal drilling [1]. Gas production in transient flow regime of multiple fractured horizontal well (MFHW) is known to be extensively elongated because of the ultra-low porosity and permeability of gas shale formation. This extended transient flow period in shale reservoirs makes production decline analysis and prediction become a challenging task, but the detailed knowledge on the production decline will be required for the gas shale production planning. It is widely reported that the pore characterization of shales is

challenging because of the wide spectrum of pore sizes and types within shales [2]. It can be widely classified as two major types according to Warren-Root fractured reservoir model [3]: (1) micro-/mesco-pores within shale matrices serving as gas storage house by adsorption and the gas transport is defined by diffusive flow due to concentration gradient; (2) macropores and/or fractures at which the gas can be stored as compressed free gas and the gas transport is dominant by pressure gradient. As a dual/triple porosity reservoir, the shale fracture will provide the gas pathways with Darcian flow and the mass transport [4]. Because of the multiscale pore structure in shale, the overall gas deliverability of shale falls into a multimechanics flow regime including Darcy flow, Knudson diffusion, slippage flow and sorption [5,6]. In addition to the complex pore structure, the shale formation also involves the stimulated fractures which will provide main fluid flow pathway during both the flowback and production stages. Therefore,

⁎ Corresponding authors at: Department of Energy and Mineral Engineering, G3 Center and Energy Institute, The Pennsylvania State University, University Park, PA 16802, USA (S. Liu), School of Energy Resource, China University of Geosciences (Beijing), Beijing 100083, China (Z. Li). E-mail addresses: [email protected] (S. Liu), [email protected] (Z. Li).

https://doi.org/10.1016/j.fuel.2018.09.101 Received 27 June 2018; Received in revised form 18 September 2018; Accepted 22 September 2018 0016-2361/ © 2018 Elsevier Ltd. All rights reserved.

Fuel 237 (2019) 283–297

T. Lu et al.

Nomenclature Bg C CD Cg Cgi D h km kfi I0(x) K0(x) K1(x) Lref Lh M Mg n p pi psc q (t ) q* qij q¯ij (s ) qsf qD r rm rw rD R Rm s

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

volume factor, dimensionless wellbore storage coefficient, m3Pa−1 dimensionless wellbore storage coefficient, dimensionless gas compressibility, Pa−1 gas compressibility at initial condition, Pa−1 diffusion coefficient, m2/s reservoir thickness, m permeability of matrix system, m2 permeability of 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 number of hydraulic fractures apparent molecular weight of shale gas, kg/kmol molar quantity of shale gas, kmol pressure of 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 sandsurface flow rate, m3/s dimensionless production rate of the line sink, dimensionless radial radius, m radial radius in spherical matrix blocks, m well radius of multiple fractured horizontal well, m dimensionless radius, dimensionless gas constant, J/(mol·K) external radius of matrix block, m variable of Laplace transformation, dimensionless

LfUi LfLi ΔLfij ΔLfDij Z ρ ρsc ϕ μg μgi σ β ψ ψL ψi Δψ Δψs ψD ω λ γD

the gas production behavior characterization and prediction is expected to be a multimechanistic flow process. The multimechanistics flow model will be required to define the unique production behavior of shale reservoirs. The pore structures of natural rocks are known to be fractal and gas flow for real state gas(es) through rocks is very complex [7]. Especially in shale, the properties of single-component diffusion of hydrocarbons in zeolites were known to be a concentration gradient drive diffusive flow [8]. It was proposed by previous researchers that the gas flow in nanopores of shale can be modeled with a diffusive transport regime with a constant diffusion coefficient [9]. A formulation was presented [9] for gas flow in the nanopores of mudrocks based on Knudsen diffusion and slippage flow and pointed out that the contributions of Knudsen diffusion increase as pores become smaller and smaller. The permeability of tight/shale gas reservoir (k < 0.1 md) is commonly pressure dependent [10,11]. With consideration of diffusive flow in shale matrices and the stress-sensitivity of fracture system, a dual-mechanism dual-porosity model was established [12] for gas transport in shale reservoirs. However, desorption of shale gas was ignored in the proposed model. It’s [13] firstly proved that the desorption phenomenon of shale gas could be described by Langmuir isotherm theory based on experimental data. Other scholars [14–21] also arrived at the

skin factor, dimensionless time, s dimensionless time, dimensionless reservoir temperature, K temperature at standard condition, K flow velocity of shale gas in 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 yaxis, m length of upper wing of ith fracture, m length of lower 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, dimensionless gas viscosity at current condition, Pa·s gas viscosity at initial condition, Pa·s adsorption index, dimensionless a parameter related to permeability modulus, Pa−1·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 dimensionless permeability modulus, dimensionless

similar conclusion that the gas sorption in shale can be described by the Langmuir model. MFHW has been proved to be the most effective well type and completion technique for the economical shale gas development. Numerical simulator (dusty-gas model) with an adjusted Knudsen diffusion coefficient [22] was used to characterize the flow regimes for tight and shale gas reservoirs, but desorption in matrix and stress-sensitivity of fracture system were ignored. Incorporating gas adsorption, stress dependence, non-Darcy flow and surface diffusion, an integrative model for shale-gas-reservoir simulation was proposed [23]. Meanwhile various semi-analytical models have been proposed to diagnose the pressure and production for MFHW shale wells. Amongst these semianalytical methods, linear flow model is widely used. Many previous studies [23–32] focused on linear gas flow modeling. These proposed linear models are simple and they can be practically used to identify the linear flow regime. However, these models ignore the fracture interference effect which could be a major decline driving parameter at low reservoir pressures. To overcome the shortcomings of linear models, other researchers [33–36] proposed improved semi-analytical models combining point source function which can obtain transient pressure solution of MFHW easily. Even though these improved models are capable of matching a portion of gas production accurately, some of the

284

Fuel 237 (2019) 283–297

T. Lu et al.

early or late production behaviors cannot be modeled and predicted because some of the critical reservoir parameters were ignored in their models, such as stress-sensitivity of permeability in fracture system. Compared with the works above, stress-sensitivity was added into a well testing model and the improved model [29] incorporated all microflow mechanisms (gas desorption, diffusion and stress-sensitivity Darcian flow) for shale gas reservoirs. This work is based on Wang’s mircoflow theory to present a semi-analytical production analysis model for shale gas reservoirs. Some improvements are also made: extending diffusion to Knudsen diffusion and slippage flow; and establishing a point sink instead of a line sink. The proposed model in this article is capable of comprehensively defining the gas production over the extended production history of the well with continuous depletion. It can also serve as a powerful tool to inversely estimate the reservoir parameters when the production data set is available through the history matching. More importantly, the validated model can be used to forecast the long-term gas production for the future shale wells with reservoir parameters as direct inputs.

Fig. 2. Micro unit of matrix pore (Modified from Liu et al. [43]).

2. Proposed gas flow modeling for gas shale wells pores. Gas is stored as compressed free gas in pores and void spaces and as adsorbed gas within the shale matrices. As shown in Fig. 1, gas flows through a network of pores with different length scales ranging from nanometers to fracture-scale. Gas desorption and adsorption happens on the internal surfaces of shale matrices, Knudsen diffusion and slip flows are the primary gas flux within the matrix and these non-Darcian gas fluxes flow towards natural fractures. As soon as the gas arrives at natural fractures, the gas flow will follow the Darcy’s law dominated by the pressure gradient. (2) Hydraulic stimulated fracture network is infinite-conductivity because the permeability of fracture network is known to be orders of

In the formation, the gases need to flow through both shale matrix system and fracture system. The matrix system is composed with microscale pores which have very low permeability. The fracture system includes both naturally occurred fracture and hydraulic stimulated fractures. It is known that the gas flow regime varies at different length scales, which makes overall gas deliverability description quite challenging. In order to mathematically define the gas flow behavior for the tight shale formation, the following assumptions are made: (1) As a dual-porosity anisotropic medium, a shale gas reservoir consists of fracture system including hydraulic stimulated fractures network and naturally occurred fracture system and micro-/meso-

Fig. 1. Gas storage and flow in shale from nanoscale to macroscale: (a) Macro-/fracture-scale: Darcy flow in fractures including induced and natural fracture; (b) Microscale: gas diffusion from matrix to fracture system; (c) Nanoscale: (I) dissolved gases transfer from organic kerogen to pore surface; (II) gas desorption serving as gas source for diffusive flows (Modified from Guo et al. [37]).

285

Fuel 237 (2019) 283–297

T. Lu et al.

Fig. 3. Schematic of a multiple fractured horizontal well: (a) physical model of MFHW with dual porosity reservoir; (b) MFHW’s discretion sketches and (c) boundary conditions.

magnitudes higher than that of natural fracture. And these stimulated fractures are assumed to be perpendicular to the lateral horizontal wellbore as illustrated in Fig. 3. (3) MFHW produces at a constant rate qsc when calculating transient pressure responses. However, each fracture’s flow rate can be different within the system with continuous pressure depletion. The initial reservoir pressure throughout the control volume is constant at pi. (4) MFHW is located in the center of cylindrical reservoir with a nonflow outer boundary whose top and bottom boundaries are also impermeable.

1 (r f vf ) = r r

(

f

f

t

)

+ (1

f

)

SC

V t

(1) 3

where ρf is gas density of fracture system, kg/m ; ρsc is gas density at standard condition, kg/m3; vf is gas velocity in fracture, m/s; ϕf is porosity of natural fracture system, fraction; and V is the volumetric gas concentration of adsorbed gas from shale matrix [13], sm3/m3. The state of natural gas in shale can be described by the equation of state (EOS) for real natural gas, thus, gas density can be defined as follows:

=

Shale gas reservoir has multiple flow mechanisms due to its multiscale pore/fracture structure from the molecular scale (nanometer) to macro scale as illustrated in Fig. 1. In order to accommodate all the flow mechanisms, each flow mechanism should be properly coupled into a holistic flow model through the governing equation. In this subsection, we will try to couple all these involved flow regimes into one integrated flow model.

Mg p

(2)

ZRT

where Mg is molar mass, kg/mol; p is gas pressure, Pa; Z is Z-factor, dimensionless; R is gas constant, 8.314 Pa⋅m3/(mol⋅K); T is absolute temperature, K. The high overburden stress can result in the closing of fractures with an exponential decrease of fracture permeability. Fractures’ closure which is caused by high overburden stress would decrease permeability in fracture system. Here, we introduce the theory of Klkanl and Pedrosa [38] to account for the counter effects. The radial velocity for gas flow in fracture system can be described by:

2.1. Seepage model for naturally occurred fracture system

vf =

Considering the mass flow rate between shale matrix and fracture in unit volume reservoir being equal to (1 − ϕf)ρsc V (kg/(m3·s)), t Continuity equation of naturally occurred fracture system in radial coordinate is as follows:

kfi e

(pi p)

pf

µg

r

(3) 2

where kfi is the initial permeability of fracture system, m ; pi is the initial reservoir pressure, Pa; γ is the permeability modulus describing

286

Fuel 237 (2019) 283–297

T. Lu et al.

stress-dependent permeability in the naturally occurred fracture system, dimensionless. μg is gas viscosity, Pa⋅s; pf is the pressure of fracture system, Pa. The pseudo pressure can be defined as: p

=

p0

2p dp µg Z

As shown in Fig. 2, an idealized spherical matrix pore will be used to define the gas transport in matrix shale pores. The mass conservation equation of this system considering diffusion and slippage can be described as Eq. (11) based on equation of continuity [44,45]:

(vrm A

(4)

f r2

+

1 r

f

r

f

+

2

r

= e(

fi

f)

f

µgi Cgi kfi

f

t

+

2pksc T (1 kfi TSC

f

)

V t

vrm A

where β is a parameter related to permeability modulus and its definition is given by Eq. (6), pa−1s; ψf is pseudo pressure of fracture system, Pa/s; ψfi is pseudo pressure of fracture system at initial condition, Pa/s; Cg is gas compressibility, expressed as Eq. (7), pa−1; psc is the pressure at standard condition, Pa; Tsc is the temperature at standard condition, K.

(

Cg =

2rm 8RT 3 M

pavg rm

m m4

r2 r)

t

=

pm M pm ZRT r

M

m

ZRT

(12)

Cgm pm 4 r 2 r

pm

(13)

t

p Fkm r2 m µg Cgm r

=

pm t

(14)

(15)

(16)

rm = Rm

where Rm is average equivalent pore radius for shale matrix, m. According to Fick’s second law, the partial differential equation of transient diffusion in shale matrix can be written as:

(8)

2 1 a

(

m

Fkm µg Cgm

V 3D V = t Rm rm

1 V V rm2 D = rm2 rm rm t

(17)

The inner boundary of matrix blocks is non-flow as gas diffusion in the spherical matrix blocks is symmetric, so

(9)

µg

=

Fkm 4 r2 µg

where D is the apparently diffusion coefficient, m /s. In Eq. (1), V represents desorption component from shale matrix. t Based on the assumption that the matrix blocks are spherical, transient diffusive flow in shale matrix is as follows according to the Fick’s second law.

0.5

0.5

(11)

2

V rm

And slippage coefficient can be mathematically defined as follows [41]:

8 RT F=1+ M

Dk +

D = Dk +

where vrm is radical velocity in matrix, m/s; Dk is the Knudsen diffusion coefficient, m2/s; pm is the matrix pressure, Pa; F is slip coefficient which is the theoretical dimensionless coefficient for a slippage effect, dimensionless; km is matrix permeability, md; μg is gas viscosity, mPa⋅s. Knudsen diffusion coefficient is defined by [5]

Dk =

VV )

t

So the apparent diffusivity can be introduced as follows:

The equivalent permeability of shale matrix is in nano-Darcy scale and Knudsen diffusion is the dominate flow in shale gas matrix [5,9]. Shale matrix is rich in nano-scaled pore whose diameter is compatible to the mean free path of gas molecules, so slippage can’t be ignored and should be properly incorporated [39]. Knudsen diffusion and slip flows contribute to the overall gas flow within the matrix pores [40], yielding gas flow velocity in matrix:

vrm

= Dk Cgm +

VV )

m m

1 r2 r

2.2. Gas flux model for shale matrix

pm r

m m

According to Eq. (7), Cg is assumed to be p in this study because the simulated reservoir pressure is less than 15 MPa at which the ideal gas law still holds in reasonable accuracy. Substituting Eqs. (12) and (13) into Eq. (11) and dividing by unit volume, 4πr2Δr yields

(7)

Dk k = +F m pm µg

(

=

-1

(6)

1 1 dp p Z dZ

m m

t

µg Z 2p

m m )r

where A is the external surface area of the sphere, m ; ρm is gas density of matrix system, kg/m3; ϕm is the matrix porosity, fraction. VV is the volume of micro unit of matrix pore, m3. Substituting gas flow velocity in matrix, i.e. Eq. (8) in to Eq. (11), we can arrive:

(5)

=

(vrm A

2

where ψ is pseudo pressure, Pa/s. Combining Eqs. (1)–(4), we can obtain the governing equation for fracture system by evaluating gas viscosity (μg) and gas compressibility (Cg) at initial condition (μg = μgi and Cg = Cgi) helps solve equations. 2

m m )r + r

=0

(18)

rm 0

The outer boundary conditions can be described as that gas concentrations in the fracture system and on the external surface of matrix blocks are approximately the same.

(10)

(19)

V |rm= Rm = VE

where rm is matrix pore radius, m; pavg is the average pressure, Pa; α is the tangential momentum accommodation coefficient [42], dimensionless. In gas shale reservoir, the major gas storage space is micro- and nano-pores. These pores have nano-scale sized diameters resulting in large specific surface areas which can be served as gas sorbing sites. Generally, gas in micro- and nano-pores is driven by concentration gradient which can be simplified using Fick’s law. In this current work, we try to model this micro-scale flow using transient diffusion method.

3

3

where VE is unit volumetric gas concentration [13], sm /m (volume of gas adsorbed per unit volume of the reservoir in equilibrium at pressure p). The initial condition is given as:

V |t = 0 = Vi

(20)

The desorption behavior of the matrix international surface can be described as:

287

Fuel 237 (2019) 283–297

T. Lu et al.

2.3. Fracture system and shale matrix coupling flow model

Table 1 Definitions of dimensionless variables. Dimensionless variables

Based on the definitions of dimensionless variables in Table 1, we can get the governing flow equations of fracture system:

Definitions

Dimensionless pseudo-pressure

D

Dimensionless radius

=

rD =

Dimensionless gas concentration difference Dimensionless unit volumetric gas concentration difference Dimensionless time

kfi hTsc psc qsc T r Lref

Dimensionless unit flow rate

=

2

i

D rD2

VD = Vi V VED = Vi VE

D

=

µg Z 2p

= Cgi +

Inter-porosity coefficient

µi

= Desorption index

=

2 kfi h (1 qsc µgi

)

psc qsc T kfi hTsc

D (rD, tD )

q¯ij Lref qsc

¯ = qDij

Storativity ratio

Cgi +

Cgi +

qsc µgi

)

(1

2

D

L

=

D

+

L

psc qsc T L

+

i )( L

drD2

+

i)

(25)

1

d 2 D¯0

VD rmD 2 rmD

rmD = 1

VD 1 VD = rmD tD

+ (1

)

VD tD

(32)

ln[1

D D (rD, tD )]

(33)

D

tD

D

+ (1

VD tD

)

(34)

+

D D1

+

2 D D2

D D (rD , tD )]

=1+

=

(35)

+ D (rD, tD)

D D (rD, tD )

+

+

1 2

2 D

(rD, tD ) +

2 D D (rD, tD )

+

(36) (37)

+

1 D0 = rD rD

D0

tD

+ (1

)

VD tD

(38)

+

1 d D¯0 = s D¯0 + (1 rD drD

) sV¯D

(39)

Substituting Eq. (31) into Eq. (39), we can get the partial differential equation of transient flow in fracture system and matrix coupling model as follows

According to the definitions of dimensionless variables in Table 1, Eqs. (16)–(20) can be rewritten as follows

VD =3 tD

D0

d 2 D¯0

L VL

kfi hTsc (

tD

Taking the Laplace transformation of Eq. (38) with respect to tD, we can obtain the following seepage model of fracture system in Laplace domain:

where σ is desorption index [29], dimensionless, which represents desorption effect and its definition is given as:

=

D0

ln[1

rD2

(24)

D

1 D

D D (rD, tD )

2

(23)

Combining Eqs. (21)–(23), we can obtain

VED =

D

D D

The zero-order perturbation solution can satisfy accuracy requirement because dimensionless permeability modulus is always very small. So Eq. (34) can be rewritten as follows

where VL is Langmuir volume, sm3/m3; ψL is Langmuir pseudo-pressure, Pa/s. Here we identify a dimensionless variant VED

VED = Vi VE

=

1 1

(22)

i

=e

rD

1 D 1 = rD rD 1 D

D

(21)

i

Vi = VL

2

The degree of nonlinearity in Eq. (34) decreases apparently compared to Eq. (32). We can perform a parameter perturbation in γD by defining the following series:

ki D psc qsc T L VL kfi hTsc ( L + i )( L + i )

1

+

+

rD2

2L2 ) Rm ref

Note: Lref is an arbitrary length and here we consider it is to be equal to rw.

VE = VL

D

D

Substituting Eq. (33) into Eq. (32) yields

Cgi 2 kfi h (1 qsc µgi 2 kfi h

1 D rD rD

+

where γD is dimensionless permeability modulus, dimensionless. Permeability of naturally occurred fracture system is a pressure dependent parameter which leads to a nonlinear seepage model. In order to solve this nonlinear equation, we can use Pedrosa’s substitution to solve that problem [46].

kfi t

tD = µgi L 2 ref

Dimensionless permeability modulus

,

drD2

(26)

+

1 d D¯0 = f (s ) D¯0 rD drD

(40)

where (27)

f (s ) = s + 3

(28)

3. Reservoir pressure and gas production relationship

VD |rmD = 1 = VED

(29)

3.1. Shale gas well pressure responses

VD |tD= 0 = 0

(30)

There are m vertical fractures with vertical fracture plane respected to formation horizontal plane. Orthogonal coordinate system is shown in Fig. 3. As permeability in hydraulic stimulated fractures is much larger than formation permeability, gas flow in stimulated fractures can be assumed as infinitely conductive.

2 rmD

VD rmD

rmD

=0 rmD 0

The solution of Eqs. (26)–(30) in Laplace space is as follows (the detailed solving procedure can be seen in Appendix A).

sV¯D = 3

[ s / coth(s / ) 1] ¯D

(31)

288

(1

)[ s / coth( s / ) 1]

(41)

Fuel 237 (2019) 283–297

T. Lu et al.

3.1.1. Point source solution with boundary conditions Ozkan and Raghavan’s method [47] made an assumption that gas of finite volume is instantaneously flowing out from a tiny sphere’s surface at t = 0. Combining with mass conservation, the cumulative reduction through the sphere surrounding the source is equal to the instantaneous flowing out volume. The well produces at a constant production q (t ) and the inner boundary condition can be expressed by: t 0

lim+ 2 kfi e

( i

)

r

0

B T 2 g sc

f

psc T

r

dt =

q

lim+ 0

0

2qsc h

e

D

rD2

D D Lref

dt =

rD

q

lim+ 0

0

2qsc h

rD

dt =

(r ) ¯

D0

(

rD

3 lim+2Lref rD2 0

¯

2 D

d D¯ 0 drD

D0

rD

(rD

)

rD =

A1 Ak

Ai

q

(44)

wDH

(1

where xD =

x

) qsc

Lref

(xD x wD )2 + (yD ywD ) 2 + (zD zwD )2 , yD =

Lref

, zD =

z Lref

, x wD =

xw , ywD Lref

S=

=

yw , Lref

CD =

(46)

zwD =

1 1 0

Am 2n, m 2n LfDm2n

2n, k k

LfDij

I0 (rD (xD , yD ; x wD , ywD ) f (s ) ) K1 (reD f (s ) )

(xD x wD )2 + (yD ywD )2 , reD =

=

=

xDi j ) 2 + (yDk

K 0 ( (xDk

v

I0 ( (xDk v

xDi j )2 + (yDk v

qDm ¯

0 0 = 0 1

2n

wD

v

) 2 f (s ) )+

)2 f (s ) ) K1 (reD f (s ) )

d (51)

I1 (reD f (s ) )

¯

s+

s wD + S 2 ¯ s CD (s wD

(52)

+ S)

kfi hTsc s

psc qsf T

(53)

C 2 2 h Lref

wD (tD )

I1 (reD f (s ) )

qDij

(54)

where qsf is the subsurface flow rate, m /s. s is an extra pseudopressure drop, Pa/s; C is the wellbore storage coefficient, m3Pa−1. ¯ in Eq. (52) can be By Stehfest numerical inversion algorithm, wDH transformed to be in real space from Laplace space. Combining with perturbation technology, the bottom hole pressure for MFHW in shale gas reservoir can be obtained as follows

¯ ¯ = ¯ = q K 0 (rD (xD , y ; x wD , y ) f (s ) ) D wD D D0 qsc

where rD (xD , yD ; x wD , ywD ) =

Am

1, m 2n

qD11

3

zw . Lref

In order to estimate the gas shale well production decline curves, it is reasonable to assume an elliptical non-flow boundary as illustrated in Fig. 3. The boundaries at z = 0 and h are both impermeable, so applying the method of images to the point source solution Eq. (46) in such case yields:

+

Ak

1, k k

1 1 1

1, m 2n

where S is the skin factor of fracture defined by Eq. (53), dimensionless; CD is the dimensionless wellbore storage coefficient defined by Eq. (54).

(45)

f (s ) (xD x wD ) 2 + (yD ywD )2 + (zD zwD ) 2 ] y

A1

1, k k

where Ai j, k v is coefficient of pseudo-pressure-drop of vth discrete unit of kth fracture caused by continuous line source of unit strength at the jth discrete unit of the ith fracture. ¯ in Eq. (50) doesn’t consider the effects of wellbore storage The wD and skin effect. The effects of wellbore storage and skin effect can be added into the above pressure solutions [48,49]:

We can get the zero-order of dimensionless pseudo-pressure as follows:

q¯ hD exp[ ¯ D0 = 2qsc

j, k v

¯

q sc µ gi

1,1 1

Ak

yDi, j + 1

~ kfi q h Cgi +

(49)

(50)

reD) = 0

µgi

A1

Am 2n,1 1 LfD11

f (s ) D¯0 = 0

2 kfih

1,1 1

yDi, j

rD =

=

¯ LfDij ] = 1 [qDij s

Eqs. (48) and (49) consist a (2n × m + 1)-order system of linear algebraic equations which can be expressed in the following matrix form.

Combining Eq. (41) and applying perturbation technique, we can get: 1 2 rD rD

(48)

(43)

rD =

D0

Lref rD2

2n

i=1 j =1

By using perturbation technique and having zero-order solution, we can obtain: t

¯

Dij (xDk, v , yDk , v )

The constant flow rate equation can be written as follows: m

where q is production rate of the line sink, m3/s; kfi is the initial permeability of the fracture system, m2; h is the effective thickness, m. The dimensionless of Eq. (42) can be expressed as follows t

2n

i=1 j =1

(42)

r=

m

¯ =

wD

=

1 D

ln [1

D wDH

(tD )]

(55)

3.2. Gas production responses based on the well pressure responses

(47)

Based on Van Everdingen and Hurst’s study [50], there is a relationship between dimensionless production, q¯D , at a constant bottom hole pressure condition and dimensionless bottom pseudo-pressure, ¯ wD , at a constant production condition and it’s expressed as follows:

re . rw

3.1.2. Bottom pressure solution of MFHW The pressure response of MFHW in shale gas reservoir will be obtained by discretizing fractures and applying principle of superposition (detailed derivations can be seen in Appendix B). The pressure responses at (xD, yD) caused by the 2n × m segments can be estimated by applying the principle of superposition.

q¯D =

1 ¯ s 2 wD

(56)

With Stehfest numerical inversion algorithm, we can get well production rate, qsc, in real space as follow

289

Fuel 237 (2019) 283–297

T. Lu et al.

qsc =

kfi hTsc

4.2. Modeling of production behavior with and without multimechanics flows

qD

psc T

(57)

In order to understand the unique production behavior with multimechanistic flows, we tried to compare the shale MFHW’s production behavior with and without consideration of sorption effect, Knudsen diffusion, slippage effects and stress-dependent permeability. Different flow regime analysis was conducted to investigate the impacts of various flow mechanisms on the overall production behaviors. Fig. 6 shows the type curves and flow regime division for a MFHW in shale gas reservoirs. Gas flow is mainly dominated by well storage effect in wellbore storage flow period. Between wellbore storage flow and early-time linear flow regimes, it is transition flow period characterized by a hump in the pressure derivative curve. Then it reaches early-time linear flow in each hydraulic fracture, which is marked by slope = 0.5 in pressure derivative curve as shown in Fig. 6. Fracture interference has not achieved and each fracture produces independently during regime III. Followed this, it reaches pseudo-radial flow around each fracture flow period, whose existence depends on whether fracture spacing is long enough or not. In intermediate-time linear flow period, pressure drop wave reaches the outside of fracture system and the stimulated reservoir can be treated as an entirety surrounding by the linear flow line. During this period, gas flow is dominated by fracture interference in fracture system and diffusion in matrix system cannot be ignored as more and more gas flow from matrix to fracture resulting in an elevated concentration gradient. Diffusion flow represents gas diffusion from matrix to fracture system since regime Ⅴ, but it dominates gas flow in regime VI. In regime VI, the adsorbed gas starts to desorb from the shale matrix and then diffuse in matrix driven by gas concentration gradient. Boundary dominated flow period can be reached after a long production time since the low permeability property. Production in this regime is often very small and production behaviors of shale will not be discussed in this article. Based on type flow regimes division and gas flow characteristics in each period, three key times (in early-time linear flow, pseudo-radial flow around each fracture and intermediate-time linear flow, respectively) are selected to analysis shale’s production behaviors. Fig. 7 shows the production decline curves of a MFHW in shale gas reservoir at constant bottom-hole pressure at three key times. Overall, it can be found that production rate of shale gas reservoirs is higher with multimechanics flow involvements (adsorption/desorption effect, Knudsen diffusion and slippage effects) than that without considering these effects. It is apparent that the early production difference is much smaller at the early-time linear flow stage and production discrepancy becomes bigger and bigger as continuous reservoir depletion. This is attributed

4. Results and discussion of shale gas production modeling for a MFHW The proposed model is a semi-analytical model which can inverse reservoir’s parameters and also predict production decline curve. By analyzing the daily production versus time using the proposed models, it was found that the production behavior of a MFHW is most sensitive with the controlling equivalent radius. Then a workflow of interpreting production data is proposed here. During the modeling exercise, the value of controlling equivalent radius is initially assumed and estimating fracture’s half-length which is strongly associated with controlling equivalent radius. Then, the production data were fitted to estimate the best input modeling parameters. Finally, the validated model can be used either to inversely estimate the reservoir parameters or predict the future production behaviors of the existing well. 4.1. Proposed model application and validation According to parametric sensitivity level, we initially validate the model then provide a simple workflow to invert unknown parameters and/or predict production rate decline curves. To analyze the production decline curve, we took the published production data from Barnett Shale as the working data set [51]. As shown in Fig. 3, the physical model for the simulated example can be treated as corresponding to a multiple fractured horizontal well where the infinitely conductive fractures extend to the top and bottom boundary of the reservoir, also the circular outer boundary. The well was fracture-stimulated in 6 stages, one frac stage contains about 5 perforation clusters, so there is total 30 clusters i.e. 30 fracture number. The matrix permeability of shale gas reservoir is assumed as 0.1 md. And the values of other parameters of simulation are given in Table 2. After doing sensitivity analyses, it turns out that equivalent controlling radius, re, and half-length of fracture, xf, have a great impact on the production decline behavior. We assume that the imagined drainage area (shown in bottom of Fig. 4) is equal to the draining area of MHFW (shown in top of Fig. 4). To obtain an estimate of xf, an estimate of re must be made based on the assumption that the two areas inside dashed in Fig. 4 are equal, as show in Eq. (58). Then a series of xf and re can be made.

4x f × L + x f2 = re2

Table 2 The values of relevant parameters [51].

(58)

As shown in Fig. 5, we set a series of xf and re combination, i.e. re = 400 m, xf = 187 m; re = 350 m, xf = 153 m; re = 300 m, xf = 120 m and run the model to get the production rate decline curves with reservoir and well parameters listed in Table 2. Comparing simulation data and real production rate data from Barnett shale in Fig. 5, it was found that there was a good fitting with the actual production data when re = 400 m, xf = 187 m except simulated results are a litter bigger than real production rate between 20 days and 100 days. Besides, we can see the trend of production rate in the following 10,000 days and even longer time which helps to predict when to shut in. Based on this modeling exercise, it is confirmed that the proposed model can be used to describe the production behavior of gas shale MHFW with considering all the possible flow regimes (mainly including transient and boundary dominated flow regimes). Therefore, the model along with the estimated input parameters will be used as base model for the sensitivity analyses and the different well production prediction.

290

Parameters

Value

Formation thickness, h (m) Porosity, ϕ (decimal) Initial fracture permeability, kfi (m2) Gas gravity, γg (decimal) Gas compressibility, Cg (Pa−1) Initial pressure, pi (Pa) Langmuir volume, VL (m3/m3) Langmuir pressure, pL (Pa) Bottom-hole pressure, pwf (Pa) The heal-to-toe perforated length of the well, L (m) Reservoir temperature, T (K) Wellbore diameter, rw (m) Gas viscosity, μi (Pa·s) Gas compressibility factor, Z (m3/m3) Apparent diffusivity, D (m2/s) Matrix equivalent pore radius, Rm (m) Constant production, qsc (m3/s)

30 0.03 1 × 10−16 0.69 4 × 10−8 2.4 × 107 6.22 3,693,498 1,724,000 700 366.5 0.2 2 × 10−5 0.97 1 × 10−12 1 × 10−8 2.78

Fuel 237 (2019) 283–297

T. Lu et al.

key times, the productions of the multimechanistics-based modeled results show 9.8%, 45.1% and 60.9%. Production increases compared the well without considering the multimechanistics flow models. This suggests that multimechanistics-based approach is capable of predicting the production of the later production behavior with flat extended production tails. 4.3. Model sensitivity analyses and parametric studies In this section, the production sensitivity analyses for different key reservoir parameters were conducted and discussed. The key reservoir parameters include initial absolute permeability in fracture system, kfi, effective controlling radius, re, fracture number, M, fracture half-length, xf and fracture spacing, df. The influences of these parameters are apparent and properly estimating these parameters are very important for future production analysis and forecasting. From Fig. 8, it can be seen that half-length of fracture primarily impacts early stage of gas production. Besides, with the increase of halflength of fracture, the declining rate of production becomes smaller. This is because the greater the fracture half-length, the bigger of the stimulated reservoir volume, and more gas can serve as source. Figs. 9 and 10 suggest that fracture number and fracture spacing also impact early stage but have little effect. The decline rate becomes smaller with the increase of fracture number as Fig. 9. In Fig. 10, it was found that the production rate seems has no apparent monotonic increase or decrease with the increase of fracture spacing. It’s easy to understand that production rate is bigger with more intensive fracture configuration, that’s why when df = 100 m, the production rate is bigger than it when df is equal to 150 m. It doesn’t mean that the smaller value of fracture spacing, the bigger of production rate is. Because there will be fractures interference when fracture’s arrangement is much more intensive. For example the production rate with df = 50 m is smaller than df = 100 m. As discussed as above, fracture shape (half-length of fracture, fracture number and fracture spacing) impacts early stage of gas production (when t < 100 days) and producing gas for this period mainly comes from fracture depletion. When the fracture depletion is done, the microscale flow will dominate the overall gas production which will be independent of the fracture permeability as shown in Fig. 11. It can be seen that permeability influences the whole period of production, bigger initial permeability of fracture determines bigger initial production rate, meanwhile production rate curve declines quicker with the increase of initial permeability of fracture. Fig. 12 shows that how equivalent controlling radius affects the production decline curve. The difference of production occurs at later producing period and becomes larger and larger. It can be seen that with the increase of equivalent controlling radius, the production curves decline slower, which means the bigger controlling area is, the bigger productivity is. As shown in Fig. 13, dimensionless permeability modulus, γD, mainly affects the later producing period. With dimensionless permeability modulus larger, production rate decreases more rapidly. According to Figs. 8–13, production decline curve is affected by many parameters and the trends are different. By extracting production rate data at certain specific time (t = 10 days, 100 days and 1000 days), we can easily obtain production difference, Δq, and difference percentage, Δqper of adjacent parameters as shown in Table 3. Averaging the difference percentage at the same time could be helpful to compare sensitivity of parameters. Trends of average difference percentage of parameters are shown in Fig. 14. Apparently the bigger value of difference percentage is, the more sensitive of relevant parameters is. From Fig. 14, we can see that fracture half-length, xf, fracture number, M, fracture spacing, df are relatively sensitive when t < 100 days, meanwhile permeability, kfi,

Fig. 4. Workflow of calculation equivalent controlling radius with half-length of fracture.

Fig. 5. Comparing simulated results and production rates of Barnett Shale MFHW.

that desorption and diffusive gas influx is increasingly important as the reservoir production. Besides, production is too low for economic producing after t > 10,000 days, so we did not analyze production behavior in diffusion and boundary dominated flow periods. Based on the modeled results, it is understood that the early gas production is dominated by the free gas from fracture system before pseudo-radial flow around each fracture and when the reservoir pressure was drawdown to a critical pressure, the micro-scale flow influx will override that pressure driven gas flow. Because of the moderate gas influx from matrix, the production period of a gas shale reservoir can be extensively extended and the moderate production rate can be achieved without the additional treatment. This unique feature can be an incentive for the investors to explore shale gas resource by a long term investment. In Fig. 7, quantitative production analyses were given. At three different

291

Fuel 237 (2019) 283–297

T. Lu et al.

Fig. 6. Typical pressure curves for MFHW in shale and its flow regimes dividing.

and equivalent controlling radius, re, are more sensitive when t > 100 days. That is to say, fracture shape mainly effects early period of production. Dimensionless permeability modulus, γD, and equivalent controlling radius, re mainly effect later period of production. Besides, initial Permeability of fracture, kfi and equivalent controlling radius, re seem to be the most two sensitive parameters, dimensionless permeability modulus, γD is the third, and fracture half-length, xf, fracture spacing, df and fracture number, M, are the last three sensitive parameters.

1000000 Production of xf=50m

Production rate, m3/d

Production of xf=100m

5. Conclusions This work conducts production decline behavior of multiple fractured horizontal wells (MFHW) in shale gas reservoirs by incorporating multiple flow mechanisms. The main conclusions can be drawn as follows:

Production of xf=150m

100000

re=300m, M=30, df=100m,

kfi=0.1md, γD=0.015

10000

1000

1

10

100

1000

10000

Time, days Fig. 8. Effect of half-length of fracture on production decline curves.

(1) The proposed coupled production model can be used to predict the long-term gas production behavior. The inclusion of desorption,

Fig. 7. Comparison of production rate and production difference percentage of MFHW with and without multimechanics in shale at three key times (with constant down hole pressure, pwf = 1.742 Mpa): Δq1 at t = 10 days in earlytime linear flow; Δq2 at t = 100 days in pseudo radial flow around each fracture; Δq3 at t = 7000 days in intermediate-time linear flow.

292

Fuel 237 (2019) 283–297

T. Lu et al.

1000000

1000000

Production of M=30

Production of re=300m Production of re=350m Production of re=400m

Production of M=25

Production rate, m3/d

Production rate, m3/d

Production of M=20

100000

10000

re=300m, kfi=0.1md, df=100m, xf=100m, γD=0.015

1000

1

10

100

1000

100000

kfi=0.1md, M=30, df=100m,

1000

10000

Time, days

1

10

100

1000

10000

Time, days

Fig. 9. Effect of fracture number on production decline curves.

Fig. 12. Effect of controlling radius on production decline curves.

1000000

1000000 Production of df=50m

Production of γD=0.015

Production of df=100m

Production of γD=0.05

Production of df=150m

Production rate, m3/d

Production rate, m3/d

xf=100m, γD=0.015

10000

100000

10000

re=300m, M=30, kfi=0.1md,

Production of γD=0.085

100000

kfi=0.1md, M=30, df=100m, xf=100m, re=300m

10000

xf=100m, γD=0.015

1000

1

10

100

1000

1000

10000

Time, days

Production rate, m3/d

Production of kfi=0.01md

100000

re=300m, M=30, df=100m,

xf=100m, γD=0.015

1

10

100

Time, days

1000

10000

(2) Meanwhile, multiple mechanisms will bring strong nonlinearities to seepage mathematical equation. To conquer this problem, we use perturbation method and the strong nonlinearities are considerably weakened. Perturbation technology is used to establish the point source solution with consideration of stress-dependent permeability which little research has done before. (3) Based on equivalent area method, we set a series of xf and re to predict the production decline curves, and combining with the field production data. Then we can inversely estimate the xf and re which are consistent with the reported data from the field. (4) Type pressure curves for shale MFHW were divided into seven flow periods and comparisons between production behaviors with and without multimechanics flows were made in three flow periods. The production discrepancy becomes larger and larger with continuous depletion and the elevated productions were found to be 9.8%, 45.1% and 60.9% in early-time linear flow, pseudo-radial flow around each fracture and intermediate-time linear flow, respectively. This is attributed that desorption and diffusion is increasingly important as depleting. (5) Using the proposed model, the sensitivity analyses were conducted and it was found that fracture half-length, xf, fracture number, M, fracture spacing, df mainly affect early period of production. Initial permeability of fracture, kfi, determines the initial production rate

Production of kfi=0.05md

1000

100

Fig. 13. Effect of dimensionless permeability modulus on production decline curves.

Production of kfi=0.1md

10000

10

Time, days

Fig. 10. Effect of fracture spacing on production decline curves.

1000000

1

1000

10000

Fig. 11. Effect of initial permeability of fracture on production decline curves.

diffusive flow, slippage effect and stress-dependent permeability is the novel method to ensure the model can precisely predict the long term production behavior, especially for the late time production stage.

293

Fuel 237 (2019) 283–297

T. Lu et al.

Table 3 Difference value and difference percentage of production rate affected by parameters at t = 10 days, 100 days and 1000 days. parameters

10 days

100 days

3

3

1000 days

Δq, m /d

Δqper, %

Δq, m /d

Δqper, %

Δq, m3/d

Δqper, %

xf, m

100/50 150/100

17,559 10,254

11.9 6.2

1625 2105

3.3 4.1

134 136

2.0 2.0

M, count

25/20 30/25

8405 3624

5.0 2.0

872 1337

1.6 2.5

134 82

2.0 1.2

df, m

100/50 150/100

8204 7147

5.0 4.1

1519 1111

2.9 2.0

160 137

2.5 2.0

kfi, md

0.1/0.05 0.05/0.01

91633 82263

46.7 78.7

40734 53666

36.0 74.2

2196 11413

54.5 32.0

re, m

350/300 400/350

6211 3182

3.3 1.6

14,196 12,153

22.4 15.7

3056 3220

37.1 28.5

γD, dimensionless

0.015/0.0505 0.05/0.85

4255 4318

2.0 2.2

4309 4442

4.1 4.4

5132 8683

30.2 73.3

80

Difference percentage, %

70

radius, re, and dimensionless permeability modulus, γD, mainly effects the later period of producing, besides, initial permeability of fracture, kfi, seems to be the most sensitive parameter, equivalent controlling radius, re, is the second, dimensionless permeability modulus, γD, is the third, and fracture half-length, xf, fracture spacing, df and fracture number, M, are the last three sensitive parameters.

Controlling radius, re Fracture number, M

Dimensionless permeability modulus, γD Initial permeability of fracture, kfi Half length of fracture, xf

Fracture spacing, df

60 50 40 30

Acknowledgment

20

This article is supported by development scale prediction and development pattern study of shale gas reservoirs, China (Grant No. 2016ZX05037-006) which are gratefully acknowledged.

10 0

0

200

400

600

800

1000

1200

Time, days Fig. 14. Trends of average difference percentage of parameters by extracting production rate data at t = 10 days, 100 days and 1000 days.

and affects the whole producing period. Equivalent controlling Appendix A The dimensionless equations for transient diffusion model for matrix is expressed as follows based on Section 2.2:

VD =3 tD

VD rmD

(A-1)

rmD = 1

1 VD 1 VD 2 rmD = 2 rmD rmD rmD tD

VD rmD

(A-2)

=0 rmD

(A-3)

0

(A-4)

VD |rmD = 1 = VED

(A-5)

VD |tD= 0 = 0 µi

2 kfi h Cgi + (1 q sc µ gi

2 ) Lref

R2

where = , = Dm , rmD = Rm , VED = Vi VE ki m Taking Laplace transformations of Eqs. (A-1)–(A-5) yields r

294

Fuel 237 (2019) 283–297

T. Lu et al.

dV¯D sV¯D = 3 drmD 1 2 rmD

dV¯D drmD

V¯D 1 = sV¯D rmD

2 rmD

rmD

(A-6)

rmD = 1

(A-7)

=0 rmD

(A-8)

0

¯ V¯D |rmD = 1 = VED

(A-9)

Eq. (A-7) becomes by making

d 2E¯ E¯ = rmD V¯D 2 drmD

s ¯ E=0

(A-10)

The general solution of Eq. (A-10) is as follows s / rmD

E = A1 e

+ A2 e

(A-11)

s / rmD

Then s / rmD

e V¯D = A1

+ A2

rmD

s / rmD

e

(A-12)

rmD

Substituting Eq. (A-12) into Eq. (A-8) yields

dV¯D drmD

s/ e

= A1

s / rmD

rmD 2 rmD

0

rmD

e

s / rmD

+ A2

s/ e

s / rmD

rmD 2 rmD

e

s / rmD

(A-13)

The following relationship can be obtained combining the inner boundary condition Eq. (A-13):

A2 =

(A-14)

A1

Then Eq. (A-12) can be expressed as

e V¯D = A1

s / rmD

s / rmD

e

rmD

rmD

(A-15)

Substituting Eq. (A-15) into Eq. (A-9) yields

¯ = A1 (e VED

s/

e

s/

(A-16)

)

After taking Laplace transformation with respect to tD, Eq. (27) can be further rewritten as follows:

¯D

¯ = VED

(A-17)

Combing Eqs. (A-16) and (A-17) yields

A1 =

(e

s/

e

s/

)

¯D

(A-18)

The dimensionless gas concentration in shale matrix in the Laplace space can be obtained by substituting Eq. (A-18) into Eq. (A-15):

V¯D =

sh ( s / )

sh ( s/ rmD ) ¯ D rmD

(A-19)

where sh ( s / ) is the hyperbolic sine function. Eq. (A-19) can be further derived as follows:

V¯D rmD

=

[ s/ coth s /

1] ¯D

(A-20)

rmD = 1

Thus Eq. (A-6) can be rewritten as

dV¯D sV¯D = 3 drmD

=3

[ s / coth s /

1] ¯D

(A-21)

rmD = 1

Appendix B As shown in Fig. 3, the upper and lower wings are both discreted into n equal segments, and we can get 2n segments in each fracture. The coordinates of mid-point for the jth segment in the ith fracture, (xij , yij ) are

x ij =

2n

2j + 1 LfUi 2n

yij = yi

,1

j

n

(B-1) 295

Fuel 237 (2019) 283–297

T. Lu et al.

x ij =

2(j

n) 2n

1

LfLi

yij = yi

,n+1

j

2n

(B-2)

where LfUi is the length of upper ith fracture, m; and LfLi is the length of lower ith fracture, m. The coordinates of the jth end point in the ith fracture (xij, yij) are n

j+1

x ij = LfUi n ,1 yij = yi j

x ij =

n n

1

LfLi

yij = yi

j

n

,n+1

(B-3)

j

2n

(B-4)

As a MFHW produces at a constant production rate and discrete unit number, n, is big enough, each discrete segment can be assumed as producing at a uniform flux density qij. The pressure response caused by the jth segment in the ith fracture then can be obtained by integrating Eq. (47).

¯ =

Dij

q¯ij qsc

K 0 (rD (xD , yD ; x wD , ywD ) f (s ) ) +

ij

I0 (xD , yD ; x wD , ywD ) f (s ) ) K1 (reD f (s ) )

dl

I1 (reD f (s ) )

(B-5)

where, l is the length of the jth segment in the ith fracture, m. As the hydraulic fracture plane is vertical to y-axes, then dl = dx. Eq. (B-5) becomes

¯

Dij (xD , yD )

q¯ij

=

xi, j + 1

Lref

qsc

K 0 (rD (xD , yD ; x wD , ywD ) f (s ) ) +

I0 (rD (xD , yD ; x wD , ywD ) f (s ) ) K1 (reD f (s ) ) I1 (reD f (s ) )

x i, j

dx wDi

Here we define the dimensionless flux per unit length of the discrete segment (i, j) as follows q~ij Lref q~Dij = qsc

(B-6)

(B-7)

So Eq. (B-6) can be rewritten as xDi, j + 1

¯ ¯ Dij (xD , yD ) = qDij

K 0 (rD (xD , yD ; x wD , ywD ) f (s ) ) +

I0 (rD (xD , yD ; x wD , ywD ) f (s ) ) K1 (reD f (s ) )

xDi, j

I1 (reD f (s ) )

dx wDi

(B-8)

So m fractures can be discrete into 2n × m segments, and the pressure responses at (xD, yD) caused by the 2n × m segments can be estimated by applying the principle of superposition. m

¯ (xD , y ) = D D

2n

¯

Dij (xD , yD )

(B-9)

i=1 j =1

Then the pressure response at the discrete segment (xDk, v , yDk, v ), (1 m

¯ (xDk, v , y ) = Dk , v D

2n

k

m, 1

v

2n) can be obtained by the following equation:

¯

Dij (xDk, v , yDk , v )

(B-10)

i=1 j =1

As we consider that hydraulic fractures are infinitely conductive, so

¯ (xDk, v , y ) = ¯ Dk , v D wD

(B-11)

Hence, Eq. (B-9) effectively becomes:

¯ =

m

2n

¯

Dij (xDk, v , yDk , v )

wD

(B-12)

i=1 j =1

By writing Eq. (B-12) for all midpoints for discrete segments k, v (k = 1, 2, m ; k = 1, 2, 2n) , we can get 2n × m equations. Solving fracture ~ ¯ requires 2n × m + 1 equations and we need one more additional independent equation. This can be constant flow rate flux qDij and pressure wD equation as follows: m

2n

[q¯ij Lfij ] =

qsc

i=1 j =1

s

(B-13)

where ΔLfij is the length of the discrete segment (i, j). Eq. (B-13) can be rewritten as: m

2n

i=1 j =1

q¯ij Lref qsc

Lfij Lref

=

1 s

(B-14)

296

Fuel 237 (2019) 283–297

T. Lu et al.

References

[26] Al-ahmadi HA, Almarzooq AM, Wattenbarger RA. Application of linear flow analysis to shale gas wells-field cases. SPE Unconventional Gas Conference. Society of Petroleum Engineers; 2010. [27] Brohi I, Pooladi-darvish M, Aguilera R. Modeling fractured horizontal wells as dual porosity composite reservoirs-application to tight gas, shale gas and tight oil cases. SPE Western North American Regional Meeting. Society of Petroleum Engineers; 2011. [28] Ozkan E, Brown M, Raghavan R, et al. Comparison of fractured-horizontal-well performance in tight sand and shale reservoirs. SPE Reserv Eval Eng 2011;14(2):248–59. [29] Wang HT. Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms. J Hydrol 2014;510:299–312. [30] Sang Y, Chen H, Yang S, et al. A new mathematical model considering adsorption and desorption process for productivity prediction of volume fractured horizontal wells in shale gas reservoirs. J. Nat Gas Sci Eng 2014;19:228–36. [31] Guo J, Zhang L, Zhu Q. A quadruple-porosity model for transient production analysis of multiple-fractured horizontal wells in shale gas reservoirs. Environ Earth Sci 2015;73(10):5917–31. [32] Deng Q, Nie RS, Jia YL, et al. A new analytical model for non-uniformly distributed multi-fractured system in shale gas reservoirs. J Nat Gas Sci Eng 2015;27:719–37. [33] Zhao Y, Zhang L, Zhao J, et al. ‘Triple porosity’ modeling of transient well test and rate decline analysis for multi-fractured horizontal well in shale gas reservoirs. J Petrol Sci Eng 2013;110:253–62. [34] Zhao Y, Zhang L, Luo J, et al. Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir. J Hydrol 2014;512:447–56. [35] Zhang D, Zhang L, Zhao Y, et al. A composite model to analyze the decline performance of a multiple fractured horizontal well in shale reservoirs. J Nat Gas Sci Eng 2015;26:999–1010. [36] Su Y, Zhang Q, Wang W, et al. Performance analysis of a composite dual-porosity model in multi-scale fractured shale reservoir. J Nat Gas Sci Eng 2015;26:1107–18. [37] Guo C, Xu J, Wu K, et al. Study on gas flow through nano pores of shale gas reservoirs. Fuel 2015;143:107–17. [38] Klkanl J, Pedrosa QA. Perturbation analysis of stress-sensitive reservoirs. SPE Form Eval 1992;6(3):379–86. [39] Swami V, Settari A. A pore scale gas flow model for shale gas reservoir. SPE Americas Unconventional Resources Conference. Society of Petroleum Engineers; 2012. [40] Darabi H, Ettehad A, Javadpour F, et al. Gas flow in ultra-tight shale strata. J Fluid Mech 2012;710:641–58. [41] Brown GP, Dinardo A, Cheng GK, et al. The flow of gases in pipes at low pressures. J Appl Phys 1946;17(10):802–13. [42] Agrawal A, Prabhu SV. Survey on measurement of tangential momentum accommodation coefficient. J Vac Sci Technol A 2008;26(4):634–45. [43] Liu P, Qin Y, Liu S, et al. Non-linear gas desorption and transport behavior in coal matrix: Experiments and numerical modeling. Fuel 2018;214:1–13. [44] Huang T, Guo X, Chen F, et al. Modeling transient flow behavior of a multiscale triple porosity model for shale gas reservoirs. J Nat Gas Sci Eng 2015;23:33–46. [45] Chen F, Duan Y, Wang K, et al. A novel pressure transient response model considering multiple migration mechanisms in shale gas reservoir. J Nat Gas Sci Eng 2015;22:321–34. [46] Pedrosa OA. Pressure transient response in stress-sensitive formations. SPE California Regional Meeting. Society of Petroleum Engineers; 1986. p. 203–14. [47] Ozkan E, Raghavan R. New solutions for well-test-analysis problems: Part 2 computational considerations and applications. SPE Form Eval 1991;6(3):369–78. [48] Kucuk F, Ayestaran L. Analysis of simultaneously measured pressure and sanface flow rate in transient well testing. J Petr Technol 1985;37(2):323–34. [49] Schroeter TV, Gringarten AC. Superposition principle and reciprocity for pressure transient analysis of data from interfering wells. SPE J 2009;14(3):488–95. [50] Van Everdingen AF. The skin effect and its influence on the productive capacity of a well. J Pet Technol 1953;5(6):171–6. [51] Clarkson CR. Production data analysis of unconventional gas wells: workflow. Int J Coal Geol 2013;109–110:147–57.

[1] Dong Z, Holditch S, McVay D. Resource evaluation for shale gas reservoirs. SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers; 2012. [2] Ma L, Slater T, Dowey PJ, et al. Hierarchical integration of porosity in shales. Sci Rep 2018;8:1683. [3] Warren JE, Root PJ. The behavior of naturally fractured reservoirs. SPE J 1963;3(3):245–55. [4] Wang Y, Zhu Y, Liu S, et al. Pore characterization and its impact on methane adsorption capacity for organic-rich marine shales. Fuel 2016;181:227–37. [5] Javadpour F, Fisher D, Unsworth M. Nanoscale gas flow in shale gas sediments. J Can Pet Technol 2007;46(10):55–61. [6] Zhang R, Liu S, Wang Y. Fractal evolution under in situ pressure and sorption conditions for coal and shale. Sci Rep 2017;7(1):1–11. [7] Al-Hussainy R, Ramey HJ, Crawford PB. The flow of real gases through porous media. J Pet Technol 1966;18:624–36. [8] Xiao J, Wei J. Diffusion of hydrocarbons theory. Chem Eng Sci 1992;47:1123–41. [9] Javadpour F. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J Can Pet Technol 2009;48(8):16–21. [10] Franquet M, Ibrahim M, Wattenbarger RA, et al. Effect of pressure-dependent permeability in tight gas reservoirs, transient radial flow. Can Int Pet Conf 2004. [11] Mckernan R, Mecklenburgh J, Rutter E, et al. Microstructural controls on the pressure-dependent permeability of Whitby Mudstone. Geol Soc 2017;454:SP454-15. [12] Ozkan E, Raghavan R, Apaydin OG. Modeling of fluid transfer from shale matrix to fracture network. SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers; 2010. [13] Bumb AC, McKee CR. Gas-well testing in the presence of desorption for coalbed methane and Devonian shale. SPE Form Eval 1988;3(1):179–85. [14] Lane HS, Watson AT, Lancaster DE. Identifying and estimating desorption from Devonian shale gas production data. SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers; 1989. [15] Gao C, Lee J, Spivey JP, et al. Modeling multilayer gas reservoirs including sorption effects. SPE Eastern Regional Conference & Exhibition. Society of Petroleum Engineers; 1994. [16] Spivey JP, Semmelbeck ME. Forecasting long-term gas production of dewatered coal seams and fractured gas shales. SPE Rocky Mountain Regional/Low Permeability Symposium. Society of Petroleum Engineers; 1995. [17] Lu XC, Li FC, Watson AT. Adsorption measurements in Devonian shales. Fuel 1995;74(4):599–603. [18] Weniger P, Kalkreuth W, Busch A. High-pressure methane and carbon dioxide sorption on coal and shale samples from the Paraná Basin, Brazil. Int J Coal Geol 2010;84(3–4):190–205. [19] Hu H, Zhang T, Wiggins-Camacho JD, et al. Experimental investigation of changes in methane adsorption of bitumen-free Woodford Shale with thermal maturation induced by hydrous pyrolysis. Mar Pet Geol 2015;59:114–28. [20] Wang Y, Zhu Y, Liu S, et al. Methane adsorption measurements and modeling for organic-rich marine shale samples. Fuel 2016;172:301–9. [21] Guo H, Jia W, Peng P, et al. Evolution of organic matter and nanometer-scale pores in an artificially matured shale undergoing two distinct types of pyrolysis: a study of the Yanchang Shale with Type II kerogen. Org Geochem 2017;105:56–66. [22] Freeman CM. A numerical study of microscale flow behavior in tight gas and shale gas reservoir systems. SPE International Student Paper Contest at the SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers; 2010. [23] Wang J, Luo H, Liu H, et al. An integrative model to simulate gas transport and production coupled with gas adsorption, non-Darcy flow, surface diffusion, and stress dependence in organic-shale reservoirs. SPE J 2017;22:244–64. [24] Aboaba A, Cheng Y. Estimation of fracture properties for a horizontal well with multiple hydraulic fractures in gas shale. SPE Eastern Regional Meeting. Society of Petroleum Engineers; 2010. [25] Bello RO, Wattenbarger RA. Multi-stage hydraulically fractured shale gas rate transient analysis. SPE North Africa Technical Conference and Exhibition. Society of Petroleum Engineers; 2010.

297