Fluid flux throughout matrix-fracture interface: Discretizing hydraulic fractures for coupling matrix Darcy flow and fractures non-Darcy flow

Journal Pre-proof Fluid Flux Throughout Matrix-Fracture Interface: Discretizing Hydraulic Fractures for Coupling Matrix Darcy Flow and Fractures non-Darcy flow Salam Al-Rbeawi, Jalal F. Owayed PII:

S1875-5100(19)30313-0

DOI:

https://doi.org/10.1016/j.jngse.2019.103061

Reference:

JNGSE 103061

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 10 August 2019 Revised Date:

7 November 2019

Accepted Date: 7 November 2019

Please cite this article as: Al-Rbeawi, S., Owayed, J.F., Fluid Flux Throughout Matrix-Fracture Interface: Discretizing Hydraulic Fractures for Coupling Matrix Darcy Flow and Fractures non-Darcy flow, Journal of Natural Gas Science & Engineering, https://doi.org/10.1016/j.jngse.2019.103061. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier B.V. All rights reserved.

Fluid Flux Throughout Matrix-Fracture Interface: Discretizing Hydraulic Fractures for Coupling Matrix Darcy Flow and Fractures non-Darcy flow Salam Al-Rbeawi, METU-Northern Cyprus Campus, Turkey Jalal F. Owayed, Kuwait University, Kuwait

Abstract The objective of this paper is focusing deep insights on reservoir fluid flux from a structurally complicated matrix to non-uniformly propagated hydraulic fractures. A hypothetical unconventional reservoir is considered with stimulated and un-stimulated reservoir volume. A Semi-analytical multilinear flow regimes model is developed for pressure distribution by coupling matrix Darcy flow model and hydraulic fracture non-Darcy flow model. Hydraulic fractures are discretized to several segments with a specific fluid flux to these segments. The study has reached to several conclusions such as: fluid flux from SRV matrix to hydraulic fractures may have a significant impact on reservoir performance and discretizing hydraulic fractures to several segments may give different pressure behavior, flow rate, and productivity index than single segment fracture. The novel point presented in this study is presenting a semi-analytical model that couples matrix Darcy flow and hydraulic fracture non-Darcy flow using trilinear dual porosity model and discretizing hydraulic fractures.

1- Introduction Oil and gas production from unconventional resources has been given great attention during the last couple of decades. Part of this attention is focused on simulating and modeling fluid transportation in structurally complicated reservoirs that are typically characterized as dual porosity and permeability media. Knowing that these reservoirs are commonly depleted by initiating hydraulic fractures because of their ultralow permeability, developing mathematical models for pressure distribution, flow rate behavior, and productivity index becomes very challenge tasks. The reason for that is the differences in the petrophysical characteristics of hydraulic fractures and reservoir matrices as well as the differences in the flow patterns that dominate fluid transportation in the matrices and fractures. Categorically and unequivocally, the reservoir matrix between hydraulic fractures may undergo some changes in its characteristics because of the impact of the stimulation process by hydraulic fractures. Microfracture networks could be developed in this area in addition to naturally existed fracture that are originally embedded in the matrix. Accordingly, the porosity and permeability of this area may grow bigger than the original porosity and permeability of the reservoir. This part of the reservoir is called stimulated reservoir volume (SRV). However, part of the reservoir might not be affected by the fracturing process i.e. hydraulic fractures may not exist in this part and only natural fractures are embedded in the matrix. This part of the reservoir is called un-stimulated reservoir volume (USRV). Because these two reservoir volumes and hydraulic fractures have different characteristics, flow patterns that could be developed in the three parts are expected to be different also. Darcy flow is considered the dominant flow pattern where reservoir fluid flows from the matrix to natural fractures in the USRV and SRV. While non-Darcy flow is considered the dominant flow pattern in hydraulic fractures where cross-section area of flow, determined by hydraulic fracture width and height, is small enough for fluid velocity to grow up significantly. Increasing fluid velocity would lead to increase inertial forces that in turn lead to non-Darcy flow. Starting from early 1900s Forchheimer, 1901 had pointed out that Darcy’s law is not applicable for such conditions and proposed the following model:





(1)



= −  +    

where   is the inertial factor that has numerous empirical models tabulated in the literature (Economides M. et al.,

2013). Several analytical solutions for pressure behaviors (Camacho-V et al. 1996, Spivey et al. 2004, Zeng and Zhao 2008) demonstrated the significant changes that could be observed due to non-linear trends of non-Darcy flow. This non-linear

scheme is most represented by the Forchheimer approach for Darcy flow given by Eq. (1). However, Barree and Conway 2004, 2009 stated that Forchheimer approach might not be applied for predicting the flow rate if high potential gradients exist. They proposed a new model for the pressure drop at high potentials and flow rates and their model converges to the Forchheimer model when high potentials are no longer present. Hydraulic fractures may provide excellent free flow paths, however, the great impact of non-Darcy flow on pressure profile is expected to occur inside these fractures. Pressure drop inside these fractures is highly affected by a lot of parameters such as cross-section area of flow, fracture conductivity, and reservoir fluid properties. The narrowing in crosssection area of flow causes an increase in pressure drop while increasing fracture conductivity leads to small pressure drop. The flow rate inside these fractures increases as the fluid flux from the matrix to the fractures increases, therefore, non-Darcy flow develops as reservoir fluid inside hydraulic fractures approaches wellbores. This fact was confirmed by (Zhang and Yang 2014; Luo and Tang 2015) who stated that non-Darcy flow in the fractures reduces the effective conductivity and change pressure profile. Not only the pressure profile that could be affected by non-Darcy flow, productivity index is very sensitive to non-Darcy flow because of the extra pressure drop required to overcome the inertial forces. Mohan et al. 2009 indicated that the productivity improvement could be overestimated by a factor on the order of two or three if non-Darcy flow is neglected. Considering Darcy flow in the matrix may not reflect the real flow pattern of reservoir fluid. According to Javadpour et al. 2007, gas flow in the shale matrix, for example, can be described by Knudson flow and slip flow in nanopores and Darcy flow in microspores while desorption and diffusion flow mechanisms are the flow patterns in solid kerogen such as the flow in Coalbeds. Recent studies have revealed that the main contributor to transferring reservoir fluid from the matrix to the fracture is the Darcy flow in the matrix driven by the pressure difference between the matrix and fractures (Ozkan et al. 2010). While the main contributor of the fluid transfer from hydraulic fractures to the wellbore is the non-Darcy flow driven by the pressure difference between fracture tips and wellbore pressure (Holditch and Morse 1976; Guppy et al. 1982 a,b; Martins et al. 1990; Vincent et al. 1999; Settari et al. 2000). Several attempts have been conducted frequently to formulate reservoir fluid transfer from the matrix to the fractures. The early attempts were focused on the flow throughout the matrix/fracture interface in a naturally fractured reservoir where a shape-based function is used to incorporate this interface into flow equations (Hughes and Blunt 2001). In these equations, saturations, flow rates, and history dependent of transfer are all explained by relative permeability. While de-Swan 1976, Kazemi et al. 1992, and Gupta and Civan 1994 modified transfer functions that have an exponential-like behavior of immersion experiments. However, these functions were unable to capture rate-dependent effects (Beckner et al. 1988; Dutra and Aziz 1992). Several theoretically generated transfer functions are obtained by developing analytical models of reservoir fluid transfer throughout matrix/fracture interface considering simple shapes and prescribed sets of matrix/fracture interface boundary conditions (Civan 1994; Rasmussen and Civan 1998; Jamili et al. 2011; Civan and Rasmussen 2012).

For the flow from the matrix with imbedded natural fractures to hydraulic fractures in unconventional reservoirs, simulating and modeling transfer phenomena via matrix/fracture interface may need different approach than shape-based functions approach that typically used in naturally fractured reservoirs. The challenges here are represented by two points. The first is the change in the flow pattern i.e. moving from Darcy to non-Darcy flow and the second is the changes in the porous media i.e. moving from structurally complicated porous media (the matrix) characterized by the existence of unknown discontinuities and other heterogeneities of varying sizes at unknown location (Raghavan et al. 2017) to simple proppantpacked porous media (hydraulic fractures). Ozkan et al. 2010 presented for the first time mathematical formulations for dual flow mechanisms in dual porosity media where hydraulic fractures are used for depleting shale matrix characterized by embedded natural fractures. In these formulations, spherical matrix block was considered and stress-dependent permeability inside natural fractures was considered also for the closure of these fractures under pressure drop. The transfer function in these formulations is defined

by matrix storativity  and interporosity flow coefficient . The proposed models in this study were solved analytically

with no consideration given to the change in the flow pattern. While Zhang and Yang 2014 introduced a semi-analytical

model that couples Darcy flow in the matrix and non-Darcy flow in hydraulic fractures. They discretized hydraulic fractures depleting reservoirs with infinite boundaries and used dual-porosity models (Warren and Root 1963; Kazemi 1969) instead of trilinear flow models (Brown et al. 2011; Ozkan et al. 2011; Al-Rbeawi 2018) to calculate wellbore pressure behavior with production time. In this paper, the trilinear dual-porosity (TDP) model is adopted to estimate pressure behavior in stimulated and unstimulated reservoir volume. Darcy flow is assumed the dominant flow pattern in these two parts of the reservoirs. The matrix blocks in the USRV and SRV are assumed having cylindrical shape with height equals to the reservoir thickness and they are fully enveloped by fractures. Inside hydraulic fractures, non-Darcy flow is considered the dominant flow pattern. Fluid flux from the matrix to the fractures is presented to the flow models by discretizing hydraulic fractures along their halflengths to several segments. To predict wellbore pressure drop, the discrete hydraulic fracture network is used for taking into

account the details of all expected changes in fluid flux from the matrix to the fractures   and fluid flux inside these

fractures   as well as fracture characteristics such as fracture width, conductivity  , minimum relative permeability to Darcy flow  , and Reynolds number  . The novel points presented in this paper are represented by proposing

three new models for three parameters that have a significant impact on the fluid flux throughout the interference between the

matrix and hydraulic fractures. The first is the characteristic length   that is defined in the literature as a function of the

matrix pore size and pore size distribution while in this study it is defined as a function of the open space to flow inside hydraulic fractures or the reciprocal of the closed space from hydraulic fracture width by the proppant sand grains. The

second is the equivalent hydraulic diameter of hydraulic fractures !  that is an indication of the cross-section area of flow inside these fractures. The third one is the Reynold number of hydraulic fractures wherein the common definition of this

number as the ratio of the inertial to viscous force has been used in this study. The three parameters are introduced to the pressure distribution model, decline rate model, and cumulative production model to predict the reservoir performance at different production time.

2- Mathematical Modeling

Consider a rectangular shape closed reservoir as the one shown in Fig. (1). This reservoir has two side boundaries

" , $  and pay zone thickness ℎ. This reservoir is depleted by multiple hydraulic fractures that drain reservoir fluid to a

horizontal wellbore and consists of two porous media. The first is the outer drainage area, called un-stimulated reservoir volume (USRV), where no hydraulic fractures and only natural fractures are embedded in the matrix. This area is

characterized by matrix permeability and porosity &'( , ∅  respectively. While the second is the inner drainage area,

called stimulated reservoir volume (SRV), where hydraulic fractures propagate and matrix permeability and porosity are

'( , ∅  respectively. The following assumptions will be utilized for the derivation of the mathematical models of

pressure distribution:

1- Hydraulic fractures symmetrically propagate in the two sides of the horizontal wellbore. 2- Hydraulic fractures fully penetrate the reservoir in the vertical direction ℎ* = ℎ. 3- Fracture half-length "*  is equal for all hydraulic fractures. 4- Reservoir thickness is uniform everywhere.

while reservoir inner and outer boundary conditions are:

+,-

/

+., . 0. , 1,

= 0.0

(2)

4& |.,06.7 = 4' |.,06.7

(3)

+,8

/

+9, 9 09 , 1,

= 0.0

(4)

4:'( |9,0;, ⁄ = 4:* |9,0;,⁄ (5)

+,=

/

= 0.0

/

=

+., . 06.7 ,

(6)

+,=

+., . 07.7 ,

>

'*?,

(7)

Figure 1: Stimulated and un-stimulated reservoir drainage areas.

Considering that the matrix block in the un-stimulated and stimulated reservoir volumes, shown in Fig. (2), have

cylindrical shape of uniform and equal radius A . The height of the matrix block in the two volumes is equal to formation

thickness ℎ. Assuming that the pressure at the surface of the matrix block is uniform also such that radial fluid flux could

be developed from the matrix bock to the fractures d. Furthermore, the flux from the matrix and imbedded natural fractures in the stimulated reservoir volume to hydraulic fractures is uniform also. Reservoir fluid in the two volumes is subjected to two velocity components (Ertekin et al. 1986). The first is the Darcy flow velocity component given by:  =

(8)

C

B +B 

+B

D

while the second is the slip velocity component given by: ' = (9)

C

EF  +B GH

+B

D

where   is the gas density, I  is the concentration in the matrix, and ! is the diffusivity coefficient of the matrix

given by (Ertekin et al. 1986): !=

(10)

J

KEF

L 

Figure 2: Matrix blocks in stimulated and un-stimulated reservoir volumes.

The two constants in Eq. (10) depend on gas type and given by Ertekin et al. 1986 as M = 31.57 and R = 0.67 if !

is expressed in TU  ⁄!.

Accordingly, reservoir fluid may flow in the matrix with a velocity determined by the resultant of the two components of Darcy and slip velocity:  =



(11)

+B

(12)

=

and: +

B

C + D

EF  +B GH

B

Considering that:

+B

+B

C + D +

B +B

C

B

(13)

X

then: +B +B

6

=

+

VW +B

D=

+B +B

(15)

=

where:

B X

YH +B X +B

\]H , _ ^B

B Z6[

(B  `B

J =  C1 + (16)

D

(B 

(14) or:

C

YH  B

D

J  is the apparent matrix permeability and depends on gas properties that in turn depend on reservoir condition

(pressure and temperature). It should be always greater than the original matrix permeability  , however, the difference

between the two permeabilities is barely observed for different reservoir conditions as it can be seen in Fig. (3). Fig. (3) demonstrates that the contribution of slip flow to the total pressure drop in the matrix may not be significant and Darcy flow is the main contributor to the pressure drop. This fact is true as long as reservoir pressure is not small ≅ 100 bcd and the

original matrix permeability is not significantly small < 1 ∗ 10gh ij.

1.0E+00 P=1000 psi P=3000 psi

1.0E-01

P=5000 psi 1.0E-02

|} , }~

1.0E-03

1.0E-04

1.0E-05

1.0E-06

1.0E-07

1.0E-08 1.0E-08

€ = . ‚ ƒ „ = …. † ‡ = ˆ‰ °‹ 1.0E-07

1.0E-06

1.0E-05

|} , }~

1.0E-04

1.0E-03

1.0E-02

1.0E-01

1.0E+00

Figure 3: Comparison of apparent and original matrix permeability.

Reservoir fluid flux from the matrix block to the fractures can be assumed instantaneously and uniformly distributed in the fracture volume. Furthermore, fractures are considered fully envelope the matrix block from four sides. Therefore, the flux continuity at the matrix-fracture interface is: ik&. = (17)



l

C

B +B  +B

D

where ik&.  is the mass flux per unit volume of fracture per unit time, ℎ is the formation thickness, and ℎ  is the fracture height. Accordingly, the volumetric flow rate given by cylindrical matrix to the fractures is: 'Y = 2nA ℎ Z (18)

6 B +B

oH  +B

_

It is worth saying that Eq. (14) is similar to the model proposed by Barree and Conway 2004, 2009 for the flow inside hydraulic fractures wherein: +l

=

+.

(= 

p qBr [

st^Br v H w= y \x

su

(19)

where  *  is the linear velocity of fluid inside hydraulic fractures, 
 is the Darcy permeability of hydraulic fractures, and   is the characteristic length. Eq. (19) can be written in a similar form to Eq. (15):

+= +.

=

(=  p z

where { is the hydraulic fracture permeability coefficient and given by: { =  +

(21)

6gBr 

v H w= D \x

C6[

(20)

or:

{ =  +

6gBr 

Z6[

(22) where:

 =

Œ1l

,1x

_

GH 1 (=

(23)



and: ! = Ž (24)

;= = >

where *  and ℎ*  are hydraulic fracture width and height respectively. Eq. (22) indicates that hydraulic fracture

permeability coefficient { depends on minimum relative permeability of fracture  , equivalent hydraulic diameter of hydraulic fractures !  that represents the diameter corresponds to the cross-section area of flow inside hydraulic fractures

i.e * ℎ* , Reynolds number of hydraulic fractures  , and characteristic length  . This coefficient can be used as an

indication for the flow pattern inside hydraulic fractures whether it is Darcy or non-Darcy flow. When the value of this

coefficient approaches { ≈ 1.0, the flow pattern is controlled by the Darcy flow equation, otherwise, it follows the nonDarcy flow equation such as the one given in Eq. (1).

The characteristic length in the abovementioned models represents the part of the hydraulic fracture width that is not filled by proppant i.e. the length from the fracture width that is open for free flow. It depends on proppant porosity and proppant grain size. It has no connection with pore size and pore size distribution of the porous media as it was defined by Barree and Conway 2004. Accordingly, high porosity proppant and uniformly distributed small proppant grain size may lead to increase the open part of hydraulic fracture to flow and increase the characteristic length and vice versa. Fig. (4) illustrates the concept of the characteristic length proposed in this study. The proposed mathematical model for the characteristic length is:

=∑ (25)

6

Figure 4: The concept of the characteristic length.

”
“”

where j•  is the diameter of the sand grains of the proppant in TU.

As much as hydraulic fracture width is occupied by proppant, the characteristic length will be small and hydraulic

fracture permeability coefficient { will be less than one as it is shown in Fig. (5). In this case, non-Darcy flow is the

dominant flow pattern inside hydraulic fractures. While increasing the free spacing of fluid flow in these fractures means

large characteristic length and hydraulic fracture permeability coefficient may be equal to one { ≈ 1.0 as it is shown in Fig. (6). In this case, Darcy flow could be the dominant flow pattern inside hydraulic fractures. Figs (5) and (6) also point out to the following comments:

1- Increasing the equivalent hydraulic diameter of hydraulic fractures !  may help developing Darcy flow as it

represents increasing in the cross-section area of flow inside hydraulic fractures and increasing in hydraulic fracture permeability coefficient {.

2- large Reynolds number of hydraulic fractures   leads to decrease hydraulic fracture permeability coefficient { and severe non-Darcy flow conditions are expected.

Series1 ‹˜™š = › ∗ › … œ˜™ = . … Series2 ‹˜™š = › ∗ › † œ˜™ = . … Series3 ‹˜™š = › ∗ ›  œ˜™ = . … Series4 ‹˜™š = › ∗ › ˆ œ˜™ = . … Series5 ‹˜™š = › ∗ › ž œ˜™ = . … Series6 ‹˜™š = › ∗ › … œ˜™ = . › ‹˜™š = › ∗ › † œ˜™ = . › Series7 ‹˜™š = › ∗ ›  œ˜™ = . › Series8 ‹˜™š = › ∗ › ˆ œ˜™ = . › Series9 ‹˜™š = › ∗ › ž œ˜™ = . › Series10

1.2

1

0.8

—

0.6

0.4

0.2 Ÿ = › ∗ ››

0 0

0.1

0.2

0.3

0.4

0.5

|}–

0.6

0.7

0.8

0.9

1

Figure 5: Hydraulic fracture permeability coefficient for small characteristic length. 1.2

1

0.8 Ÿ = › ∗ ›

Series1 ‹˜™š = › ∗ ›… œ˜™ = . … Series2 ‹˜™š = › ∗ ›† œ˜™ = . … Series3 ‹˜™š = › ∗ › œ˜™ = . … Series4 ‹˜™š = › ∗ ›ˆ œ˜™ = . … Series5 ‹˜™š = › ∗ ›ž œ˜™ = . … Series6 ‹˜™š = › ∗ ›… œ˜™ = . › ‹˜™š = › ∗ ›† œ˜™ = . › Series7 ‹˜™š = › ∗ › œ˜™ = . › Series8 ‹˜™š = › ∗ ›ˆ œ˜™ = . › Series9 ‹˜™š = › ∗ ›ž œ˜™ = . › Series10

—

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5

|}–

0.6

0.7

0.8

0.9

1

Figure 6: Hydraulic fracture permeability coefficient for large characteristic length.

2-1-Fluid flow in un-stimulated reservoir volume (USRV) Two types of flow occur in the un-stimulated reservoir volume. The first is the flow from the matrix to natural fractures embedded in the matrix and the second is the flow from these fractures to the fractures embedded in the matrix of stimulated reservoir volume. Fig. (7) depicts fluid flow in the un-stimulated reservoir volume. Assuming a cylindrical matrix block in

the un-stimulated reservoir volume with a radius of A . Furthermore, natural fractures are fully envelope the matrix block from four sides and the matrix block height is equal to the formation height. The flow equation from the matrix to the natural fractures is (Al-Rbeawi 2017b):

6

+

, +,

(26)

CA

where:  & =

(27)

A =

+B, +,

D =  &

+B, +¡,

∅Y¢ ⁄ B∅Y¢ ⁄ =



B

(28)

Using the following boundary conditions:

4 |,07.7 = 0.0

(29)

4 |,06.7 = 4& |,06.7 (30)

the solution of Eq. (26) is: £ K'¥   4: = ¤£ K'¥B- , 4& ¤

(31)

B-

While the flow equation from natural fractures in the USRV to the natural fractures in the SRV is:

+¦ ,+.,

+

+B, +,

(32)

/

, 06.7

=  

+,+¡,

Using the boundary conditions given in Eqs. (2) and (3), Eq. (32) can be solved to:

Y§'¨K©- .1, g., ª 4:&'( = 4:' . Y§'¨K©-

(33)

where:

1, g6ª

«& =   − Kc¬ &

(34) =

(35)

£s K'¥B- 

£¤ K'¥B- 

¦  .=

¦ l B

∅Y¢ ⁄ l

  = ∅Y

(36)

¢ ⁄ =

where ­6  and ­7 are modified Bessel function of the first kind.

Figure 7: Fluid flow in the un-stimulated reservoir volume.

2-2-Fluid flow in stimulated reservoir volume (SRV) Three types of flow are expected to be developed in the stimulated reservoir volume as it is shown in Fig. (8). The first is flow from the matrix to natural fractures in the SRV, the second is the flow from natural fractures of USRV to natural fractures of SRV, and fluid flow from natural fractures of SRV to hydraulic fractures. Assuming the matrix block has a cylindrical shape of similar radius and size to the one in USRV, the solution of flow equation from the matrix to the natural fractures in this volume is similar to Eq. (31) except that:  ' =

∅Y¢ ⁄ B8 ∅Y¢ ⁄ =

(37)

While the flow equation from natural fractures in the SRV to hydraulic fractures is:

+¦ ,8 +9,

+

(38)

+B,

/

+,  06.7 ,

+

6

+,-

/

91, V?, +., . 06.7 ,

=  

+,8 +¡,

Using the boundary conditions given in Eqs. (4) and (5), Eq. (38) can be solved to:

4:' =

(39)

Y§'¨K©8 91, g., ª 4: Y§'¨K©8 91, g;, ⁄ª *

where:

«' = Tc  + (40)

K©-

91, V?,

UM®ℎ¨K«& " − 1ª − Kc &

£s K'¥B- 

£¤ K'¥B- 

Figure 8: Fluid flow in the stimulated reservoir volume.

2-3-Fluid flow in hydraulic fractures (HF) The flow inside hydraulic fractures is controlled by fluid flux from the matrix of stimulated reservoir volume and its imbedded natural fractures and hydraulic fracture characteristics such as cross-section area of flow and fracture conductivity.

Fluid flux is the quantity of reservoir fluid in volume base RR¯/! that is transferred from the matrix to hydraulic fracture in

1.0 TU of hydraulic fracture half-length "* . Non-Darcy flow is recognized as the flow pattern inside hydraulic fractures

(Barree and Conway 2004, 2009; Zhang and Yang 2014; Al-Rbeawi 2017a; Al-Rbeawi 2018) because of the narrow crosssection area of flow. At the matrix-hydraulic fracture interface:

+= +.

=

(41)

+B +B

Then, the flux given by Eq. (18) can be used to calculate dimensionless flux:

 =

±B8] .= ²8]

(42) = +

The flow equation inside hydraulic fractures is written as:

& +.

(43)

C{

+= +.

D+

±B8] ;= 

= ∅³¡ *

+= +¡

whose dimensionless form is: +

+.,

C{

(44)

+,= +.,

D+

±B, *?,

=

+,= +¡,

Discretizing hydraulic fractures along fracture half-length to several segments ´ as it is shown in Fig. (9), Eq. (44) can

be written for each segment d as follows:

+

+.,”

(45)

C

+,=” +.,”

D+

±B,”

z” *?,”

=

6 +,=”

z” +¡,

The boundary conditions for each segment are:

+,=”

/

=

/

=−

+.,” . ,”ts

(46)

+,=”

+.,” . ,”

±l,”

/

±l,”

/

z” *?,” . ,”ts z” *?,” . ,”

(47)

where µ  is the fluid flux inside hydraulic fractures: ±l”

µ = ²

8] oH

Accordingly, the solution of Eq. (44) is:

∗ ∗ 4:; = ∑¶ :µg6 ∗ 4 "µg6  − :µ ∗ 4 "µ  − µ06.7 

where:

±:l,

∗ = :µg6

∗ :µ =

∗ :µ =

±:l, z

z

/

/

(48) 6

*?,” '”

∗ :µ

.,”ts

(50)

.,”

(51)

±:B,”

and:

z”

(49)

(52)

4: "  = '*

>

?, K©= ¡J·K©= 

(53)

Figure 9: Discretizing hydraulic fractures to several segments.

Eq. (53) is obtained by applying the boundary conditions given in Eqs. (6) and (7) and explained very well in the literature (Brown et al. 2011; Ozkan et al. 2011; Al-Rbeawi 2016). It gives wellbore pressure drop of horizontal well intersected by multiple hydraulic fractures assuming uniformly distributed fractures with uniform conductivity and

dimensions and uniform fluid flux along fracture half-length. The parameter «*  in Eq. (53) refers to the configurations of

the stimulated and un-stimulated reservoir volumes in addition to the petrophysical properties of the two volumes. Mathematically, it is written: «* =

o=

*?,

where:

+

'

¥l

(54)

¸* = K« UM®ℎ¨K« $ −  ⁄2ª

(55)

« = V

or

?, 91,

(56)

+ cTc

¸ = Ž UM®ℎ ¹Ž " − "µ º '

'

¥l

¥l

(57) Because of the ultralow permeability of unconventional reservoirs, pressure pulse resulting from the depletion process may take a very long time to reach pseudo-steady state. Therefore, the transient flow approach (dual porosity model) is used to describe the unsteady state flow of reservoir fluid wherein an increasing pressure drop starts at the matrix/fracture interface and moves with time further into the matrix. In this approach the storativity function is: Tc = 1 + Ž ½' UM®ℎŽ »¼

½¼' »



where  is the storativity of hydraulic fractures, given by (Kazemi 1969; de Swan O 1976; Serra et al. 1983): =

∅Y¢ B8 ∅Y¢ =

(59)

(58)



and  is the fluid interporosity flow coefficient, defined by:  = 12 (60)

B8 =

"*

The pressure drops in the abovementioned models are used in the pseudo-pressure form:  ¿

4 = 2 ¾7 , , j4Á À

(61)

Eq. (49) represents the solution for wellbore pressure drop assuming constant sandface flow rate. The solution for

sandface flow rate assuming constant wellbore pressure, in dimensionless form, can be obtained by Duhamel principle (van Everdingen 1949) : =

6

' ¦ :F,

(62)

while cumulative production rate can be calculated using the following model (Helmy and Wattenbarger 1998): Â = ´

6

' Ã :F,

(63)

The mathematical formulations of the abovementioned models are explained in Appendix-A while the definitions of

dimensionless parameters, used in these models, are given in Appendix-B.

3- Pressure behavior Pressure behavior in fractured reservoirs is affected by a lot of parameters. Some of them belong to reservoir characteristics such as reservoir configurations and petrophysical properties while others belong to hydraulic fracture characteristics such as fracture dimension, geometry, and conductivity. It is also controlled by the dominant flow pattern in the porous media and fractures whether it is Darcy or non-Darcy flow. Non-Darcy flow is expected to be developed inside hydraulic fractures because of the small cross-section area of flow offered by these fractures. The narrowing in the area of

flow leads to increase fluid velocity inside hydraulic fractures that in turn causes increasing the inertial forces and increasing pressure drop inside fractures. However, Darcy flow is considered the dominant flow pattern in porous media. Transferring reservoir fluid throughout the interface between SRV matrix and hydraulic fractures may have also a significant impact on pressure behavior. The reason for this impact is the continuous changes in the characteristics of both sides of the interface along hydraulic fracture half-length and the changes in the flow patterns between the two sides themselves.

Applying Eq. (49) for reservoir with rectangular drainage area of " = 2.0 , $ = 1.0 consists of equal size SRV

and USRV requires discretizing hydraulic fractures to several segments and thereafter different matrix-fracture fluid fluxes

  are used. The results of pressure and pressure derivative is shown in Fig. (10). In this figure, hydraulic fractures are assumed to have constant fracture width  = 5 ∗ 10g  and constant fracture conductivity  = 10.0 along hydraulic

fracture half-length while minimum relative fracture permeability is  = 0.1. Moderate non-Darcy flow inside hydraulic

fractures is applied and represented by Reynolds number  = 1 ∗ 10 . Hydraulic fractures are discretized to a different number of segments with equal length "µ . For example, ´ = 1.0,  = 1 ∗ 107  means that Eq. (53) can be applied

directly to calculate wellbore pressure drop. Pressure behavior and pressure derivative for a single segment are shown in Fig.

(10) by the yellow color curves. Increasing the number of the segments means decreeing fluid flux   through the matrixhydraulic fracture interface and decreasing fluid flux inside hydraulic fractures  . This would give a different pressure drop than the case of the single segment. Discretizing hydraulic fractures to more than ´ = 1000.0 segments may not

change the results in a similar trend if the number of segments is ´ = 100.0. It is important to emphasize that discretizing

hydraulic fractures does not change the time for reaching pseudo-steady state flow, but it decreases pseudo-steady state pressure. 1.0E+05

1.0E+04

ǘœ ∗ É ˜œ = ….  ∗ ›.  „ œ = ˆ. ∗ ›g ʘ = ›.  |}– = . › ‹˜™š = ›

1.0E+03

1.0E+02

Ąœ & Ɯ ÇÄÁ„œ

‡–ÍÎϘÍÆ šÐф

1.0E+01

Ąœ

1.0E+00

Ɯ ÇÄ′„œ

1.0E-02

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

Ɯ

1.0E+00

Series2 È}œ =

›.  ∗ ›g

È}œ = Series5

›.  ∗ ›g›

Series3 È}œ =

1.0E-01

Series4 È}œ =

È}œ = Series6 1.0E+02

1.0E+04

›.  ∗ ›g† ›.  ∗ ›g… ›.  ∗ ›

1.0E+06

1.0E+08

Figure 10: Pressure behavior of fractured reservoirs for different fluid flux.

Fig. (11) shows the impact of different non-Darcy flow conditions on the pressure behavior of the reservoirs mentioned

in Fig. (10). In this figure, hydraulic fractures are discretized to ´ = 100.0,  = 1 ∗ 10 . Applying Eq. (53) without

considering non-Darcy flow and Reynolds number gives pressure behavior and pressure derivative that are overlaid by

pressure behavior and pressure derivative of small Reynolds number  ≤ 1 ∗ 106  obtained by applying Eq. (49).

Wellbore pressure drop exhibits significant increases with increasing Reynolds number because of increasing non-Darcy flow impact. Reynolds number in hydraulic fractures is an indication of the growing inertial forces caused by the continuously increasing in fluid velocity as the cross-section area of flow offered by fractures is very narrow compared to the fluid flux inside these fractures. Practically, Reynolds numbers may not have a direct connection with the fluid flux through the matrixhydraulic fracture interface but it has a connection with the fluid flux inside these fractures. As reservoir fluid flows towards

wellbore, the flux   increases, hence, Reynolds number increases also as the fluid approaches wellbore and sever non-

Darcy flow conditions might be developed. Similar to the impact of fluid flux from the matrix, Reynolds number may not change the starting time of pseudo-steady state, but it increases pseudo-steady state pressure.

The impact of Reynolds number on pressure behavior is not affected by fluid flux from the matrix to the fractures.

Similar behaviors are obtained for two different fluid fluxes  = 0.1 and  = 0.01 for small and large Reynolds

number as it is depicted in Fig. (12). This fact would be explained physically be the dependency of Reynolds number on fluid flux inside hydraulic fractures but not fluid flux throughout the matrix-fracture interface. 1.0E+05

1.0E+04

ǘœ ∗ ɘœ = ….  ∗ ›.  „ œ = ˆ. ∗ ›g ʘ = ›.  |}– = . › È}œ = . ›

1.0E+03

Ąœ & Ɯ ÇÄÁ„œ

1.0E+02

1.0E+01

Ąœ

1.0E+00

‡–ÍÎϘÍÆ šÐф

1.0E-01

Series3= ›.  ∗ › › ‹ ˜™š

Ɯ ÇÄ′„œ

1.0E-02

1.0E-03 1.0E-08

Series2  ‹ ˜™š = ›.  ∗ ›

1.0E-06

1.0E-04

1.0E-02

Ɯ

1.0E+00

… Series4 ‹ ˜™š = ›.  ∗ › † Series5 ‹ ˜™š = ›.  ∗ ›  ‹ Series6 ˜™š = ›.  ∗ ›

1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 11: Pressure behavior of fractured reservoirs for different hydraulic fracture Reynold number.

1.0E+05

1.0E+04

ǘœ ∗ ɘœ = …. ∗ ›. „ œ = ˆ. ∗ ›g ʘ = ›.  |}– = .›

1.0E+03

‹˜™š = › ∗ ›

Ąœ & Ɯ ÇÄÁ„œ

1.0E+02

1.0E+01

Ąœ

1.0E+00

‡–ÍÎϘÍÆ šÐф ‹˜™š = › ∗ ›

1.0E-01

Ɯ ÇÄ′„œ

1.0E-02

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

Ɯ

1.0E+00

Series2 È}œ = Series21 È}œ =

1.0E+02

1.0E+04

›. ∗ ›g… ›. ∗ ›g›

1.0E+06

1.0E+08

Figure 12: Pressure behavior for different hydraulic fracture Reynolds number and fluid flux.

Hydraulic fracture minimum relative permeability to Darcy flow   may have some effects on pressure behavior

during early and intermediate production time but not at late pseudo-steady state production as it can be seen in Fig. (13). It causes a significant decrease in the value of hydraulic fracture permeability coefficient { as it is explained in Figs. (5) and (6). Decreasing minimum relative permeability indicates significant growth in non-Darcy flow conditions that causes an

increase in the pressure drop. When this parameter reaches  = 1.0, hydraulic fracture permeability coefficient becomes

{ = 1.0 also i.e. hydraulic fracture permeability is considered in its maximum value that allows for Darcy flow to be developed inside these fractures. However, when  < 1.0, hydraulic fracture permeability coefficient becomes also

{ < 1.0 and thereafter non-Darcy flow is developed because fracture permeability is in its minimum value. Similar to

Reynolds number, it has been found that the effects of hydraulic fracture minimum relative permeability on pressure response are not different for different matrix-fracture fluid flux   regardless of the value of   or {. This is true since there

is no connection between the amount of fluid fluxing to the fracture and the permeability inside hydraulic fractures. However, both   and { are affected by changing the flow pattern inside hydraulic fractures from Darcy to non-Darcy flow. This change occurs as a result of increased fluid flux inside hydraulic fractures  .

1.0E+04

1.0E+03

ǘœ ∗ ɘœ = ….  ∗ ›.  „œ = ˆ.  ∗ ›g ʘ = ›.  ‹˜™š = › È}œ = . ›

1.0E+02

Ąœ & Ɯ ÇÄÁ„œ

1.0E+01

1.0E+00

Ąœ

1.0E-01

‡–ÍÎϘÍÆ šÐф

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

Ɯ

1.0E+00

›. 

|}– = Series5

. ˆ

Series3 |}– =

Ɯ ÇÄ′„œ

1.0E-02

Series2 |}– =

|}– = Series4

. ˆ . ›

|}– = . › Series6 1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 13: Pressure behavior of fractured reservoirs for different minimum relative permeability of hydraulic fractures.

Hydraulic fracture conductivity   is the index of essay flow inside hydraulic fractures. It is a function of fracture

permeability and dimensions (width and length). In this study, fracture conductivity in all segments is assumed constant because equal length segments are used in Eq. (49). Fracture conductivity has significant influences on pressure behavior as it is shown in Fig. (14). Increasing this conductivity means increasing the easiness of fluid flow inside these fractures and

decreasing pressure drop. Similar to fracture conductivity, hydraulic fracture width   has significant influences on

pressure behavior as it increases the cross-section area of flow. Even though increasing fracture width leads to increase hydraulic fracture Reynolds number and causes increasing in non-Darcy flow impact, it also increases fracture conductivity. The resultant effect of these parameters positively decreases pressure drop inside these fractures as it is exhibited in Fig. (15). 1.0E+04

1.0E+03

ǘœ ∗ É ˜œ = ….  ∗ ›.  „ œ = ˆ. ∗ ›g È }œ = . › |}– = . › ‹˜™š = ›

1.0E+02

Ąœ & Ɯ ÇÄÁ„œ

1.0E+01

1.0E+00

Ąœ

‡–ÍÎϘÍÆ šÐф

1.0E-01

1.0E-02

1.0E-05 1.0E-08

1.0E-04

1.0E-02

ÊӜ = Series5

›.  ∗ › …

ÊӜ = Series4

Ɯ ÇÄ′„œ 1.0E-06

›.  ∗ › g›

Series3 ÊӜ =

1.0E-03

1.0E-04

Series2 ÊӜ =

ÊӜ = Series6

Ɯ

1.0E+00

1.0E+02

1.0E+04

›.  ∗ ›  ›.  ∗ › › ›.  ∗ › †

1.0E+06

1.0E+08

Figure 14: Pressure behavior of fractured reservoirs for different hydraulic fracture conductivity. ǘœ ∗ ɘœ = ….  ∗ ›.  |}– = . › È}œ = . › ‹˜™š = ›

1.0E+04

1.0E+03

1.0E+02 ‡–ÍÎϘÍÆ šÐф

Ąœ & Ɯ ÇÄÁ„œ

1.0E+01

Ąœ

1.0E+00

Ɯ ÇÄ′„œ

1.0E-01

Series2 „œ =

›.  ∗ ›gˆ

Series5 „œ =

ˆ.  ∗ ›g

Series3 „œ = Series4 „œ =

1.0E-02

1.0E-03 1.0E-08

„œ = Series6 1.0E-06

1.0E-04

1.0E-02

Ɯ

1.0E+00

1.0E+02

1.0E+04

ˆ.  ∗ ›gˆ ›.  ∗ ›g ›.  ∗ ›g†

1.0E+06

1.0E+08

Figure 15: Pressure behavior of fractured reservoirs for different hydraulic fracture width.

Hydraulic fracture width is not different than Reynolds number and minimum relative hydraulic fracture permeability regarding the responses of wellbore pressure for different values of matrix-fracture fluid flux. It has been found that these responses are similar for all fracture regardless of the value of fluid flux since there is no relationship between them. The only exception, in this case, is hydraulic fracture conductivity wherein different pressure behaviors are observed for different conductivities and different fluid flux as it is illustrated in Fig. (16). The reason for that is the dependency of the two parameters on the length of the segment "  that is a function of the number of segments and hydraulic fracture half-length

"* . Increasing the number leads to decrease "  as it is calculated by "* /´, therefore, fluid flux throughout matrix-

fracture interface is decreased and fracture conductivity is decreased also. Accordingly, the pressure drop of bad fracture

conductivity  = 1 ∗ 10g6  and large value of matrix-fracture fluid flux  = 1 ∗ 10g6  is higher than the pressure drop of good fracture conductivity  = 1 ∗ 10½  and small value of fluid flux  = 1 ∗ 10g .

1.0E+04

1.0E+03

ǘœ ∗ ɘœ = ….  ∗ ›.  „ œ = ˆ. ∗ › g È}œ = . › |}– = . ›

ʘ = › ∗ ›g›

1.0E+02

Ąœ & Ɯ ÇÄÁ„œ

1.0E+01

1.0E+00

‡–ÍÎϘÍÆ šÐф

Ąœ

1.0E-01 ʘ = › ∗ ›†

1.0E-02

1.0E-03

Ɯ ÇÄ′„œ

1.0E-04

1.0E-05 1.0E-08

1.0E-06

1.0E-04

1.0E-02

Series2 È}œ Series5 È}œ

Ɯ

1.0E+00

1.0E+02

1.0E+04

= ›.  ∗ ›g… = ›.  ∗ ›g›

1.0E+06

1.0E+08

Figure 16: Pressure behavior for different hydraulic fracture conductivities and fluid fluxes.

4- Flow rate and cumulative production Assuming constant wellbore pressure, flow rate and cumulative production behaviors are calculated by applying Eqs. (62) and (63) respectively. Considering fluid flux throughout matrix-fracture interface could change these two behaviors from those obtained by applying mathematical models that do not consider fluid flux such as the one given in Eq. (53). Using different fluid fluxes gives different behaviors of flow rate and cumulative production rate as it can be seen in Fig. (17). Fig.

(17) represents reservoir configuration of " = 2.0 and $ = 1.0 while hydraulic fracture width and conductivity are

 = 5.0 ∗ 10g  and  = 10.0 respectively. The flow pattern inside hydraulic fractures is assumed non-Darcy flow

with Reynolds number and minimum relative permeability of  = 100.0 and  = 0.1 respectively. It is clearly understood that discretizing hydraulic fracture half-length to several segments and smaller fluid flux  = 1 ∗ 10g  may

give more flow rate and cumulative production than the flow rate and cumulative production obtained by single segment

hydraulic fractures and large fluid flux = 1 ∗ 107  . This would be explained by the minimum pressure drop that can be

obtained for the maximum number of segments and minimum fluid flux from the matrix to the fractures.

The two trends of flow rate and cumulative production are not affected by changing the number of the segments and changing the amount of fluid flux. The flow rate, for all fluid fluxes, shows a sharp decline at early production time followed by gradual decline during intermediate time while very sharp decline is observed at late pseudo-steady state production. At late production time when pseudo-steady state flow is reached, there is no significant impact for the number of the segments or the amount of fluid flux on flow rate behavior as it is shown in Fig. (17). Cumulative production, for all fluid flux, exhibits gradual increase at early production time followed by sharp increase during intermediate production time while pseud-steady state at late production time shows constant cumulative production with time but different for different fluid flux. Flow rate and cumulative production behaviors for different number of segments and different fluid fluxes are not changed with changing other parameters such as hydraulic fracture width and

conductivity or changing the flow pattern inside these fractures from Darcy to non-Darcy flow as it is demonstrated in Fig. (18). 1.0E+4 „œ = ˆ.  ∗ ›

1.0E+3

ǘœ g

∗ ɘœ = ….  ∗ ›.  ʘ = ›.  |}– = . › ‹˜™š = ›

Ȝ

1.0E+2

1.0E+5 1.0E+4 1.0E+3 1.0E+2

1.0E+1

Ȝ 1.0E-1

‡–ÍÎϘÍÆ šÐф

1.0E+1 1.0E+0

×Ĝ

ÄΘÖ~Ñ − ÎƘ~É ÎÆƘ šÐф

1.0E+0

1.0E-1 1.0E-2

Series2 È }œ = Series3 È }œ =

1.0E-3

×Ĝ

1.0E-4

1.0E-5 1.0E-08

Series4 È }œ =

1.0E-06

1.0E-04

Series5 È }œ =

È }œ = Series6 1.0E-02

›.  ∗ ›g

1.0E-2

›.  ∗ ›g† ›.  ∗ ›g…

1.0E-3

›.  ∗ ›g›

1.0E-4

›.  ∗ ›

Ɯ

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E-5 1.0E+08

Figure 17: Flow rate and cumulative production behaviors for different fluid fluxes.

1.0E+5 „œ = ˆ.  ∗ › g

1.0E+4

Ȝ

1.0E+3 1.0E+2

ǘœ ∗ ɘœ = ….  ∗ ›.  ʘ = ›.  |}– = . › ‹˜™š = ›

1.0E+5 1.0E+4 1.0E+3 1.0E+2

Ȝ

1.0E+0 1.0E-1

1.0E+1 ÄΘÖ~Ñ − ÎƘ~É ÎÆƘ šÐф

‡–ÍÎϘÍÆ šÐф

1.0E-1

1.0E-2

1.0E-2 g ÈSeries2 }œ = ›.  ∗ ›

1.0E-3

1.0E-5

1.0E-06

1.0E-04

1.0E-3

g† ÈSeries3 }œ = ›.  ∗ ›

×Ĝ

1.0E-4

1.0E-6 1.0E-08

g… ÈSeries4 }œ = ›.  ∗ ›

ÈSeries5 }œ = ›.  ∗ ›

1.0E-4

g›

 ÈSeries6 }œ = ›.  ∗ ›

1.0E-02

Ɯ

1.0E+00

1.0E+02

1.0E-5

1.0E+04

1.0E+06

Figure 18: Flow rate and cumulative production behaviors for different fluid fluxes.

5- Transient and pseudo-steady state productivity index Productivity index of constant sandface flow rate is calculated by (Ozkan et al. 2011):

Ô± = (64)

6

F, gÕ,

1.0E+0

×Ĝ

1.0E+1

1.0E-6 1.0E+08

where: 4;  is the pressure drop obtained by converting Eq. (49) to real time domain and 4  is calculated by the

following model presented by Raghavan et al. 2017: 4 =

(65)

>¡,

¼.1, 91,

while productivity index of constant wellbore pressure is calculated by (Helmy and Wattenbarger 1998): Ô = (66)

6

s Ù g Ý Ø, ¦ÚÛ1, Ü1, Z ØÕ, _ ,

  is the flow rate calculated by converting Eq. (62) to real time domain and ´  is the cumulative production

calculated by converting Eq. (63) to real time domain.

Figs. (19) and (20) show transient and pseudo-steady state productivity index for constant sandface flow rate approach

Ô±  and constant wellbore pressure approach Ô . The first one for severe non-Darcy flow conditions represented by  = 0.1 and  = 100 while the second for weak non-Darcy flow conditions wherein  = 0.5 and  = 10.

The following comments are pointed out from these two figures:

1- Discretizing hydraulic fractures to several segments and using small values for fluid flux gives productivity index higher than applying single segment models and a large amount of fluid flux. This fact is confirmed by the maximum flow rate and minimum pressure drop that are obtained by large numbers of segments and small fluid flux compared with those obtained for a small number of segments and large fluid flux. 2- The impact of matrix-fracture fluid flux on the two types of productivity index is more significant when non-Darcy flow is the dominant flow in hydraulic fractures. 3- Transient and pseudo-steady state productivity indices increase as non-Darcy flow conditions are weakened especially for small and moderate number of segments, however, for a large number of segments and small fluid flux, the indices may not have significant differences. 4- Sharp declines are seen in the two productivity index behaviors at early production time and gradual declines are observed during intermediate production while constant behavior characterizes the two productivity indices during late production time when pseudo-steady state flow is reached.

5- For all cases and conditions, the productivity index obtained by the constant sandface flow rate approach Ô±  is greater than the index obtained by applying constant wellbore pressure approach Ô . The differences between the two productivity indices are more significant for strong non-Darcy flow conditions.

6- The starting time of pseudo-steady state flow is affected by fluid flux from the matrix to hydraulic fractures and by the dominant flow pattern inside these fractures. This time is reached faster when hydraulic fractures are discretized to a large number of segments and small fluid flux is used. However, this time could be very long for sever non-Darcy flow conditions and small fluid flux.

ßÑÐ~ ÐÏ͘ΠޜÈ

= ›. ∗ ›g

1.0E+03

Series5 È}œ

= ›. ∗ ›g† = ›. ∗ ›g›

1.0E+02

= ›. ∗ ›

1.0E+01

Series3 È}œ Series4 È}œ

1.0E+02

È}œ Series6

1.0E+01 œɘ~ ÐÏ͘ΠޜÄ

ÄΘÖ~Ñ − ÎƘ~É ÎÆƘ šÐф

‡–ÍÎϘÍÆ šÐф

1.0E-01

= ›. ∗ ›g…

1.0E+00

ޜÈ

1.0E+00

Series2 È}œ

1.0E-02

ޜÄ

1.0E+03

1.0E-01

1.0E-02

ǘœ ∗ ɘœ = ….  ∗ ›.  „œ = ˆ.  ∗ ›g ʘ = ›.  |}– = .› ‹˜™š = › 1.0E-03 1.0E-08 1.0E-06 1.0E-04 1.0E-02 1.0E+00

Ɯ

1.0E+02

1.0E+04

1.0E+06

1.0E-03 1.0E+08

Figure 19: Productivity index for different fluid flux and strong non-Darcy flow conditions. 1.0E+04

Series5 È}œ

= ›. ∗ ›g† = ›. ∗ ›g›

1.0E+03

= ›. ∗ ›

1.0E+02

È}œ Series6

1.0E+02

= ›. ∗ ›g…

œɘ~ ÐÏ͘ΠޜÄ

1.0E+01

ޜÄ

ޜÈ

1.0E+04

Series4 È}œ

ßÑÐ~ ÐÏ͘ΠޜÈ

1.0E+00

= ›. ∗ ›g

Series3 È}œ

1.0E+03

1.0E+01

Series2 È}œ

‡–ÍÎϘÍÆ šÐф

1.0E+00

ÄΘÖ~Ñ − ÎƘ~É ÎÆƘ šÐф

1.0E-01

1.0E-01

ǘœ ∗ ɘœ = ….  ∗ ›.  „ œ = ˆ. ∗ ›g ʘ = ›.  |}– = . ˆ ‹˜™š = › 1.0E-02 1.0E-08 1.0E-06 1.0E-04 1.0E-02

Ɯ

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E-02 1.0E+08

Figure 20: Productivity index for different fluid flux and weak non-Darcy flow conditions.

6- Analytical and numerical verifications The approach presented in this study for utilizing fluid flux throughout matrix-hydraulic fractures interface in the flow equations of fractured reservoirs is verified analytically and numerically. In this approach, fluid flux is assembled to the flow equations by discretizing the half-length of hydraulic fractures to several segments and applying trilinear dual porosity (TDP) model to generate pressure, flow rate, and cumulative production behaviors. Unfortunately, the topic of fluid flux from the matrix to hydraulic fractures has not been given significant attention until the moment. Ozkan et al. 2010 presented dualmechanism dual-porosity model and new transfer function for reservoir fluid flowing in ultralow permeability porous media. This model was solved analytically for different conditions considering spherical shape matrix blocks that are fully enveloped

by fractures. In this model, hydraulic fractures are considered similar to the natural fractures embedded in the matrix. For the analytical validation of the proposed approach, the data given by Ozkan et al. 2010 and shown in Table-1 are used. Pressure and pressure derivative for the abovementioned data in Table-1are calculated analytically by the proposed

models in this study for " = 2.0, $ = 1.0,  = 0.5,  = 250, and  = 4 ∗ 10gá . The results are

plotted and compared with the pressure and pressure derivative calculated analytically by Ozkan et al. 2010 as shown in Fig.

(21). It is very clear that good matching is obtained between the two analytical solutions with very minor differences, especially at late production time. I think the results obtained by the proposed models are more accurate than those obtained by Ozkan et al. 2010 at late production time as both pressure and pressure derivative curves tend to have identical straight

line of slope 1.0 indicating the starting of pseudo-steady state flow. This behavior is not developed very well in the two curves of Ozkan et al. 2010.

Numerical validation is conducted by using the numerical simulator (IMEX 2015.10, CMG). The 3D grid system is (20x40x1) with cell dimension of (100 ft x100 ft x 25 ft). The reservoir and hydraulic fractures configuration is shown in Fig. (22) while reservoir and fluid data are shown in Table-2. Hydraulic fractures are discretized to (1000) segments with (0.5 ft)

length for each segment. Fluid flux through matrix-hydraulic fractures interface is assumed constant  = 0.001. NonDarcy flow is introduced to the flow equation by the minimum relative hydraulic fracture permeability  = 0.5 and hydraulic fracture Reynolds number  = 1000. The input data for the proposed models are shown in Table-3.

The pressure behavior of the reservoir of interest is calculated by the proposed models in this study and plotted first.

Then, a comparison for the plotted pressure behavior with the results of the numerical simulator is shown in Fig. (23). The agreement between the two techniques is considered reasonable even though very minor differences are observed especially at very early production time. These differences might be resulted from the assumption of constant hydraulic fracture conductivity along these fractures in the proposed models. The impact of fracture conductivity is seen always at early production time when hydraulic fracture linear flow regime is developed. Therefore, considering varied conductivity along these fractures could help improving the result during early production time.

Figure 21: Comparison for the analytical solutions using the proposed models and Ozkan et al. 2010.

Figure 22: Reservoir and hydraulic fractures configuration.

1.0E+07

1.0E+06

Series1 ‡™ÏÎ ÎÆÖ~É

Series2 ×Ö}˜–ÏÓÐ ÎÏ}ÖÐÆі

∆ÄĄš , ãÎυ /Óã

1.0E+05

1.0E+04

1.0E+03

1.0E+02

1.0E+01 1.0E-06

1.0E-04

1.0E-02 1.0E+00 1.0E+02 ĖÑ~ÖÓÆÏÑÍ ÆÏ}˜,™–Î

1.0E+04

1.0E+06

Figure 23: Pressure behavior calculated by the proposed models and CMG.

7- Limitations Even though the proposed models in this study consider fluid flux through matrix-hydraulic fracture interface as an attempt to eliminate the errors for neglecting the impact on pressure behaviors, more considerations should be given also to the fluid flux from un-stimulated reservoir volume (USRV) to the stimulated volume (SRV). The reason for that is the

variability of stimulated matrix permeability and stimulated matrix block size and shape in the SRV because of the fracturing treatment. The characteristics of the SRV continuously change with the distance from hydraulic fracture face to the so-called no-flow boundary between these fractures (Fuentes-Cruz et al. 2014; Al-Rbeawi 2018). Accordingly, discretizing this distance may give more accurate results for pressure behavior, flow rate, and productivity index. The matrix in this study is assumed consisting of a single layer matrix block of cylindrical shape with a height equal to the formation thickness. More accuracy might be obtained by assuming multilayers of matrix blocks with smaller height than the single layer. Furthermore, hydraulic fracture dimensions (width and height) are assumed uniform along fracture halflength. In reality, both may change and this would cause changing non-Darcy flow conditions inside these fractures by changing the cross-section area of flow.

8- Conclusions 1- Fluid flux throughout the matrix-hydraulic fracture interface may have a significant impact on pressure behavior, flow rate, and productivity index of unconventional reservoirs. It is highly recommended considering this flux in the flow equations otherwise the accuracy of the results is questionable. Discretizing hydraulic fractures into several segments helps introducing fluid flux to the flow equations. Pressure behavior, flow rate, and productivity index are affected by the number of these segments because fluid flux changes with changing the number of the segments. 2- For all reservoir pressures, the contribution of slip flow to the total pressure drop in the matrix is not significant. Darcy flow is the main contributor to the pressure drop in the matrix.

3- The characteristic length   in the permeability coefficient model of hydraulic fractures { depends on the proppant porosity and proppant size. It does not depend on the matrix pore size and pore size distribution as it is defined before. In

this study, the characteristic length represents the open space for reservoir fluid to flow inside hydraulic fractures. It is a function of the proppant grain size and defined as the reciprocal of the closed space by the proppant inside hydraulic fractures. The smaller the proppant grain size is the bigger the characteristic length is. Accordingly, increasing this length leads to an increase in the permeability coefficient that in turn leads to reduce the wellbore pressure drop. 4- Avery slight increase in the hydraulic fracture width cause a significant reduction in the wellbore pressure drop. This would be explained by increasing the equivalent hydraulic diameter of hydraulic fracture that in turn causes an increase in the permeability coefficient. However, increasing the equivalent hydraulic diameter causes also an increase in the Reynold number of hydraulic fractures and the impact of non-Darcy flow. 5- Flow pattern in the matrix and inside hydraulic fractures as well as petrophysical characteristics of these two porous media may affect fluid flux throughout the interface between them. If they have the same flow pattern, fluid flux is in its minimum impact. 6- Reynolds number and minimum relative hydraulic fracture permeability increase the non-Darcy flow conditions and thereafter increase the pressure drop regardless of the number of the segments and the amount of fluid flux. The smaller value of Reynolds number gives the minimum pressure drop while increasing minimum relative permeability and approaching to (1.0) gives minimum wellbore pressure drop. 7- The calculated productivity index with discretizing hydraulic fractures and considering fluid flux is higher than the index obtained by assuming single segment hydraulic fractures. The stabilized pseudo-steady state productivity index exhibits a

higher value for small fluid flux throughout the interference between the matrix and fractures regardless of the impact of non-Darcy flow. ¸ ³ ³¡ ! j•  Ô Ô± ℎ ℎ*    ' & ê; ´ 4 4 4* 4& 4' 4 4* 4' 4& 4 ∆4 4; Á U "4; ì'Y    µ  A A c U U * " $ " "* $ ñ ò ∅

 * '

äMc TåAiMUdå® å¯¬iæ TM³UåA äMc ³åibAæccdRd¯dU$, bcd g6 çåUM¯ ³åibAæccdRd¯dU$, bcd g6 è$jAM¬¯d³ jdMiæUæA åT ℎ$jAM¬¯d³ TAM³U¬Aæ, TU 4AåbbM®U jdMiæUæA, TU è$jAM¬¯d³ TAM³U¬Aæ ³å®j¬³Ud dU$ 4Aåj¬³Ud dU$ d®jæ" åT ³å®cUM®U æ¯¯RåAæ bAæcc¬Aæ, jdiæ®cd宯æcc 4Aåj¬³Ud dU$ d®jæ" åT ³å®cUM®U T¯å AMUæ, jdiæ®cd宯æcc åAiMUdå® Uℎd³®æcc, TU è$jAM¬¯d³ TAM³U¬Aæ ℎædéℎU, TU 4æAiæMRd¯dU$, ij êMUAd" bæAiæMRd¯dU$, ij êd®di¬i Aæ¯MUd æ bæAiæMRd¯dU$ åT ℎ$jAM¬¯d³ TAM³U¬Aæ Uå !MA³$ T¯å 4æAiæMRd¯dU$ åT cUdi¬¯MUæj AæcæA ådA 寬iæ, ij 4æAiæMRd¯dU$ åT ¬® − cUdi¬¯MUæj AæcæA ådA 寬iæ, ij êå¯æ³¬¯MA ædéℎU !diæ®cd宯æcc ³¬i¬¯MUd æ bAåj¬³Udå® 4Aæcc¬Aæ, bcd êMUAd" bAæcc¬Aæ, bcd è$jAM¬¯d³ TAM³U¬Aæ bAæcc¬Aæ, bcd 4Aæcc¬Aæ d® ¬®cUdi¬¯MUæj AæcæA ådA 寬iæ, bcd 4Aæcc¬Aæ d® cUdi¬¯MUæj AæcæA ådA 寬iæ, bcd !diæ®cd宯æcc bAæcc¬Aæ !diæ®cd宯æcc bAæcc¬Aæ d® ℎ$jAM¬¯d³ TAM³U¬Aæc !diæ®cd宯æcc bAæcc¬Aæ d® ®MU¬AM¯ TAM³U¬Aæc åT cUdi¬¯MUæj AæcæA ådA 寬iæ !diæ®cd宯æcc bAæcc¬Aæ d® ®MU¬AM¯ TAM³U¬Aæc åT ¬® − cUdi¬¯MUæj AæcæA ådA 寬iæ !diæ®cd宯æcc iMUAd" bAæcc¬Aæ 4Aæcc¬Aæ jAåb, bcd ë毯RåAæ bAæcc¬Aæ jAåb, jdiæ®cd宯æcc 4Aæcc¬Aæ jæAd MUd æ, jdiæ®cd宯æcc äMc T¯å AMUæ, êí³T/! !diæ®cd宯æcc T¯å AMUæ ¯¬dj T¯¬" UℎAå¬éℎå¬U iMUAd" − ℎ$jAM¬¯d³ TAM³U¬Aæ d®UæATM³æ, jdiæ®cd宯æcc ¯¬dj T¯¬" d®cdjæ ℎ$jAM¬¯d³ TAM³U¬Aæ d®UæATM³æ, jdiæ®cd宯æcc ¯¬dj T¯¬" d®cdjæ ℎ$jAM¬¯d³ TAM³U¬Aæ d®UæATM³æ, êTU ½ /! è$jAM¬¯d³ TAM³U¬Aæ æ$®å¯jc ®¬iRæA Mjd¬c, TU êMUAd" R¯å³ AMjd¬c, TU

Nomenclatures:

îMb¯M³æ åbæAMUåA çdiæ, ℎAc çdiæ, jdiæ®cd宯æcc è$jAM¬¯d³ TAM³U¬Aæ djUℎ, TU ï − ³ååAjd®MUæ TåA M båd®U d® Uℎæ båAå¬c iæjdM ð − ³ååAjd®MUæ TåA M båd®U d® Uℎæ båAå¬c iæjdM æcæA ådA R嬮jMA$, TU è$jAM¬¯d³ TAM³U¬Aæ ℎM¯T − ¯æ®éUℎ, TU æcæA ådA R嬮jMA$, TU IåibAæccdRd¯dU$ TM³UåA ódc³åcdU$, ³b êMUAd" båAåcdU$ ­®æAUdM¯ TM³UåA ¯¬dj æ¯å³dU$, TU/íæ³ !MA³$ T¯å æ¯å³dU$ ³åibå®æ®U îd®æMA æ¯å³dU$ d®cdjæ ℎ$jAM¬¯d³ TAM³U¬Aæc í¯db T¯å æ¯å³dU$ ³åibå®æ®U

  {

êMUAd" T¯å æ¯å³dU$ äMc jæ®cdU$, ­R/TU ½ IℎMAM³UæAdcUd³ ¯æ®éUℎ, 1/TU è$jAM¬¯d³ TAM³U¬Aæ bæAiæMRd¯dU$ ³åæTTd³æ®U

Subscripts  T i¬ ic

è$jAM¬¯d³ TAM³U¬Aæc ´MU¬AM¯ TAM³U¬Aæc êMUAd" åT ¬® − cUdi¬¯MUæj AæcæA ådA 寬iæ êMUAd" åT cUdi¬¯MUæj AæcæA ådA 寬iæ

References: Al-Rbeawi, S., 2018. Bivariate Log-Normal Distribution of Stimulated Matrix Permeability and Block Size in Fractured Reservoirs: Proposing New Multilinear-Flow Regime for Transient-State Production. SPEJ, preprint. doi:10.2118/189993-PA. Al-Rbeawi, S., 2017. How much stimulated reservoir volume and induced matrix permeability could enhance unconventional reservoir performance. JNGSE, 46 (2017), pp.764-781. doi.org/10.1016/j.jngse.2017.08.017 Al-Rbeawi, S., 2017. Analysis of pressure behaviors and flow regimes of naturally and hydraulically fractured unconventional gas reservoirs using multilinear flow regimes approach. JNGSE, 45 (2017), pp.637-658. doi.org/10.1016/j.jngse.2017.06.026 Barree, R. D., & Conway, M. W., 2004. Beyond Beta Factors: A Complete Model for Darcy, Forchheimer, and Trans-Forchheimer Flow in Porous Media. SPE paper presented at the SPE annual technical conference and exhibition held in Houston, Texas, USA, Sep. 26-29. doi:10.2118/89325-MS. Barree, R. D., & Conway, M., 2009. Multiphase Non-Darcy Flow in Proppant Packs. SPE Production & Operation, 24 (02), pp. 257-268. doi:10.2118/109561-PA. Beckner, B. L., Firoozabadi, A., & Aziz, K., 1988. Modeling Transverse Imbibition in Double-Porosity Simulators. Paper presented at the SPE California regional meeting held in Long Beach, California, USA March 23-25. doi:10.2118/17414-MS Brown, M., Ozkan, E., Raghavan, R., & Kazemi, H., 2011. Practical Solutions for Pressure-Transient Responses of Fractured Horizontal Wells in Unconventional Shale Reservoirs. SPE Reservoir evaluation & engineering, 14(6), pp. 663-676. doi:10.2118/125043-PA. Camacho-V., R., Vasquez-C., M., Roldan-C., J., Samaniego-V., F., & Macias-C., L., 1996. New Results on Transient Well Tests Analysis Considering Non-Laminar Flow in the Reservoir. SPE Formation Evaluation, 11 (04), pp. 237-243.doi: 10.2118/26180-PA. Civan, F., 1994. Numerical Simulation by the Quadrature and Cubature Methods. Paper presented at the SPE International conference and exhibition of Mexico held in Veracruz, Mexico, Oct. 10-13. doi:10.2118/28703-MS Civan, F., & Rasmussen, M. L., 2012. Parameters of Matrix/Fracture Immiscible-Fluids Transfer Obtained by Modeling of Core Tests. SPE Journal, 17 (02), pp. 540-554. doi:10.2118/104028-PA De Swaan O., A., 1976. Analytic Solutions for Determining Naturally Fractured Reservoir Properties by Well Testing. SPE Journal, 16 (03), pp. 117-122. Trans., AIME 261.doi:10.2118/5346-PA Dutra, T. V., & Aziz, K., 1992. A New Double-Porosity Reservoir Model for Oil/Water Flow Problems. SPE Reservoir Engineering, 7 (024), pp. 419-425. doi:10.2118/21248-PA Economides, M. J. et al. 2013. Petroleum Production System. Second Edition. Prentice-Hall, Westford, MA, USA. Ertekin, T., King, G. A., & Schwerer, F. C., 1986. Dynamic Gas Slippage: A Unique Dual-Mechanism Approach to the Flow of Gas in Tight Formations. SPE Formation Evaluation, 1 (01), pp. 43-52. doi:10.2118/12045-PA Forchheimer, P. 1901. Wasserbewguang Durch Boden. ZVDI 45: 1781. Fuentes-Cruz, G., Gildin, E., & Valko, P. P., 2014. Analyzing Production Data from Hydraulically Fractured Wells: The Concept of Induced Permeability Field. SPE Reservoir Evaluation & Engineering, 17 (02), pp220-232. doi:10.2118/163843-PA. Guppy, K. H., Cinco-Ley, H., & Ramey, H. J., 1982. Pressure Buildup Analysis of Fractured Wells Producing at High Flow Rates. JPT, 34 (11), pp. 2656-2666. doi:10.2118/10178-PA Guppy, K. H., Cinco-Ley, H., Ramey, H. J., & Samaniego-V., F., 1982. Non-Darcy Flow in Wells with Finite-Conductivity Vertical Fractures. Society of Petroleum Engineers Journal, 22 (05), pp. 681-698. doi:10.2118/8281-PA Gupta, A., & Civan, F., 1994. An Improved Model for Laboratory Measurement of Matrix to Fracture Transfer Function Parameters in Immiscible Displacement. Paper presented at the SPE 69th Annual technical conference and exhibition held in New Orleans, Louisiana, Sep. 25-28. doi:10.2118/28929-MS Helmy, M. W., & Wattenbarger, R. A.,1998. New Shape Factors for Wells Produced at Constant Pressure. Paper presented at the SPE gas technology symposium held in Calgary, Canada, March 15-18. doi:10.2118/39970-MS.

Holditch, S. A., & Morse, R. A.,1976. The Effects of Non-Darcy Flow on the Behavior of Hydraulically Fractured Gas Wells (includes associated paper 6417). JPT, 28 (10), pp.1169-1179. doi:10.2118/5586-PA. Hughes, R. G., & Blunt, M. J., 2001. Pore-Scale Modeling of Multiphase Flow in Fractures and Matrix/Fracture Transfer. SPE Journal, 6 (02), pp. 126-136. doi:10.2118/71297-PA Jamili, A., Willhite, G. P., & Green, D., 2011. Modeling Gas-Phase Mass Transfer Between Fracture and Matrix in Naturally Fractured Reservoirs. SPE Journal, 16 (04), pp. 795-811. doi:10.2118/132622-PA Javadpour, F., Fisher, D., & Unsworth, M., 2007. Nanoscale Gas Flow in Shale Gas Sediments. JPCT, 46 (10), pp.55-61. doi:10.2118/0710-06. Kazemi, H., Gilman, J. R., & Elsharkawy, A. M., 1992. Analytical and Numerical Solution of Oil Recovery from Fractured Reservoirs With Empirical Transfer Functions (includes associated papers 25528 and 25818). SPE Reservoir Engineering, 7 (02), pp. 219-227. doi:10.2118/19849-PA Kazemi, H. 1969. Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution. SPE Journal, 9 (04), pp.451-462.doi:10.2118/2156-PA. Luo, W., & Tang, C., 2015. A Semianalytical Solution of a Vertical Fractured Well with Varying Conductivity Under Non-Darcy-Flow Condition. SPE Journal, 20 (05), pp. 1028-1040. doi:10.2118/178423-PA Martins, J. P., Milton-Tayler, D., & Leung, H. K., 1990. The Effects of Non-Darcy Flow in Propped Hydraulic Fractures. paper presented at the SPE 65th Annual technical conference and exhibition held in New Orleans, Louisiana, USA, Sep. 23-26. doi:10.2118/20709-MS Mohan, J., Pope, G. A., & Sharma, M. M., 2009. Effect of Non-Darcy Flow on Well Productivity of a Hydraulically Fractured GasCondensate Well. SPE Reservoir Evaluation & Engineering, 12 (04), pp. 576-585. doi:10.2118/103025-PA. Ozkan, E., Raghavan, R. S., & Apaydin, O. G., 2010: Modeling of Fluid Transfer from Shale Matrix to Fracture Network. Paper presented at the SPE Annual technical conference and exhibition held in Florence, Italy, Sep. 19-22. doi:10.2118/134830-MS. Ozkan, E., Brown, M. L., Raghavan, R., & Kazemi, H.,2011. Comparison of Fractured-Horizontal-Well Performance in Tight Sand and Shale Reservoirs. SPE Reservoir engineering & evaluation 14 (02), pp.248-256. doi:10.2118/121290-PA. Raghavan, R., Chen, C., & DaCunha, J. J., 2017. Nonlocal Diffusion in Fractured Rocks. SPE Reservoir evaluation & engineering, 20 (02), pp. 383-393. doi:10.2118/184404-PA. Rasmussen, M. L., & Civan, F., 1998. Analytical Solutions for Waterfloods in Fractured Reservoirs Obtained by an Asymptotic Approximation. SPE Journal, 3 (03), pp. 249-252. doi:10.2118/50969-PA Settari, A., Jones, J. R., & Stark, A. J., 2000. Analysis of Hydraulic Fracturing of High Permeability Gas Wells to Reduce Non-Darcy Skin Effects. JCPT, 39 (05), pp. 56-63. doi:10.2118/00-05-04 Spivey, J. P., Brown, K. G., Sawyer, W. K., & Frantz, J. H., 2004. Estimating Non-Darcy Flow Coefficient from Buildup-Test Data with Wellbore Storage. SPE Reservoir Evaluation & Engineering, 7(04), pp.256-269.doi: 10.2118/88939-PA. Van Everdingen, A.F., 1949. The application of the Laplace transformation to flow problems in reservoirs. Petroleum Transaction, AIME, 186, pp.305-324. doi:10.2118/949305-G.

Vincent, M. C., Pearson, C. M., & Kullman, J., 1999. Non-Darcy and Multiphase Flow in Propped Fractures: Case Studies Illustrate the Dramatic Effect on Well Productivity. Paper presented at the SPE Western regional meeting held in Alaska, USA, May 26-28. doi:10.2118/54630-MS Warren, J.E., and Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE Journal, 3(03), pp.245-255.doi:10.2118/426-PA. Zeng, F., & Zhao, G., 2008. Semi-analytical Model for Reservoirs with Forchheimer's Non-Darcy Flow. SPE Reservoir Evaluation & Engineering, 11 (02), pp. 280-291. doi: 10.2118/100540-PA. Zhang, F., & Yang, D.,2014. Determination of Fracture Conductivity in Tight Formations with Non-Darcy Flow Behavior. SPE Journal, 19 (01), pp. 34-44. doi:10.2118/162548-PA

Conversion factors bbl x 1.589873

E-01=m

3

cp x 1.0

E-03=Pa.s

ft

x 3.048

E-01=m

in

x 2.54

E-00=cm

Psi x 6.894757

E+00=kPa

Ibm x 0.453592

E+00=Kg

md x 0.9869

2

E+11=cm

Appendix-A 1-Fluid flow in un-stimulated reservoir volume (USRV) Two types of flow occur in the un-stimulated reservoir volume. The first is the flow from the matrix to natural fractures imbedded in the matrix and the second is the flow from these fractures to the fractures imbedded in the matrix of stimulated reservoir volume. Assuming cylindrical matrix block in the un-stimulated reservoir volume with a radius of A . The flow equation from the matrix to the natural fractures is (Brown et al., 2011; Al-Rbeawi 2017b):

6 +

 +

1)

CA

+B +

D = ∅ò³¡ ⁄ &

+B +¡

(A-

Knowing that:

4 = 2)

U =

= ∆

6²8] W 7.777ô½õ= ¡ ¦ ∅Y¢ = .=

3)

 & = 4)

A = 5)

(A-

∅Y¢ ⁄ B∅Y¢ ⁄ =

(AZ .B _ 

l



(A-



B

(A-

eq. (A-1) can be converted to dimensionless form as: 6

+

, +,

CA

+B, +,

D =  &

+B,

and in Laplace domain: 6

+

, +,

CA

+:B, +,

+¡,

(A-6)

D −  & c4: = 0.0

(A-7)

Using the following boundary conditions:

4: |,07.7 = 0.0

(A-

8)

4: |,06.7 = 4:& |,06.7

(A-

9)

the solution of Eq. (A-7) is: £ K'¥B- ,  4: = ¤ 4:& £¤ K'¥B- 

10)

+¦ +. ¦

11)

(A-

The flow equation from natural fractures in the USRV to the natural fractures in the SRV is: +

 6 +B

/

l  + 0B

= ∅ò³¡ ⁄ 

++¡

(A-

Knowing that: ∅Y¢ ⁄ l

  = ∅Y

¢ ⁄ =

12)

and: = 13)

(A-

¦  .=

¦ l B

(A-

then, Eq. (A-11) can be written in dimensionless form as:

+¦ ,+.,

+

+B, +,

/

14) +¦ :

=  

+,+¡,

(A-

In Laplace domain, Eq. (A-14) is: ,-

+.,

+

+:B, +,

/

15) The term ¹

+:B, +,

/

16) +¦ :

, 06.7

, 06.7

, 06.7

+B,

/

−   c4:& = 0.0

+,  06.7 ,

= Kc &

(A-

º in Eq. (A-15) can be determined from Eq. (A-10):

£s K'¥B-  4: £¤ K'¥B-  &

(A-

Substitute Eq. (A-16) in Eq. (A-15) and re-arrange: ,-

+.,

17)

− öc  − Kc &

£s K'¥B- 

£¤ K'¥B- 

÷ 4:& = 0.0

(A-

Using the boundary conditions given in Eqs. (2) and (3), Eq. (A-17) can be solved to:

Y§'¨K©- .1, g., ª 4:' 4:& = . Y§'¨K©-

18)

1, g6ª

where:

«& = c  − Kc & 19)

£s K'¥B- 

£¤ K'¥B- 

(A-

(A-

where ­6  and ­7 are modified Bessel function of the first kind. 2-Fluid flow in stimulated reservoir volume (SRV)

Assuming the matrix block has a cylindrical shape of similar radius and size to the one in USRV, the solution of flow equation from the matrix to the natural fractures in this volume is similar to Eq. (A-1) except that:  ' = 20)

Z _

∅Y¢ ⁄ B8 B ∅Y¢ ⁄ =

.l



while the flow equation from natural fractures in the SRV to hydraulic fractures is:

(A-

+¦ 8

+

+9 ¦

 6 +B

/

l  + 0B

+

6 6 +-

/

91 . +. .0.l

= ∅ò³¡ ⁄ 

+8 +¡

(A-

21) Eq. (A-21) can be converted to dimensionless form:

+¦ ,8 +9,

+

+B,

/

+,  06.7 ,

22)

+

and in Laplace domain: +¦ :,8 +9,

+

+:B,

23)

/

+,  06.7 ,

The term ¹

+B,

+

6

+,-

/

=  

6

+:,-

/

=   Tc4:'

91, V?, +., . 06.7 ,

91, V?, +., . 06.7 ,

/

+,  06.7 ,

+,8 +¡,

(A-

(A-

º in Eq. (A-23) can be determined from Eq. (A-10) with a consideration given to the stimulated

reservoir volume instead of the unstimulated reservoir volume.

+:B, +,

/

24)

, 06.7

= Kc '

while the term ¹

+,-

+:,-

/

+., . 06.7 ,

£s K'¥B8  4: £¤ K'¥B8  '

/

+., . 06.7 ,

(A-

º is determined from Eq. (A-18):

= −K«& UM®ℎ¨K«& " − 1ª4:'

(A-

25) +¦ :

Substitute Eqs. (A-24), and (A-25) in Eq. (A-23) and re-arrange: ,8

+9,

26)

− ö  c +

6

91, V?,

K«& UM®ℎ¨K«& " − 1ª − Kc ' £

£s K'¥B8  ¤ K'¥B8 

÷ 4:' = 0.0

(A-

Using the boundary conditions given in Eqs. (4) and (5), Eq. (A-26) can be solved to:

4:' =

27)

Y§'¨K©8 91, g., ª 4: Y§'¨K©8 91, g;, ⁄ª *

where:

«' = c  + 28)

K©-

91, V?,

UM®ℎ¨K«& " − 1ª − Kc '

(A-

£s K'¥B8 

£¤ K'¥B8 

(A-

3-Fluid flow in hydraulic fractures (HF) The flow inside hydraulic fractures is controlled by fluid flux from the matrix of stimulated reservoir volume and its imbedded natural fractures and hydraulic fracture characteristics such as cross section area of flow and fracture conductivity. At the matrix-hydraulic fracture interface:

+= +.

29)

=

+B +B

The flux given by Eq. (18) can be used to calculate dimensionless flux:

(A-

 =

±B8] .= ²8]

30)

= +

(A-

The flow equation inside hydraulic fractures is written as:

& +.

31)

C{

+= +.

D+

= ∅³¡ *

±B8] ;= 

+= +¡

(A-

whose dimensionless form is: +

+.,

32)

C{

+,= +.,

D+

=

±B, *?,

+,= +¡,

Discretizing hydraulic fractures along fracture half-length to several segments ´, Eq. (A-32) can be written for each

segment d as follows: +

+.,”

33)

(A-

C

D+

+,=” +.,”

±B,”

z” *?,”

=

6 +,=”

z” +¡,

(A-

The boundary conditions for each segments are:

+,=”

/

=

+.,” . ,”ts

34)

+,=”

/

±l,”

/

(A-

±l,”

/

(A-

z” *?,” . ,”ts

= −z *

+.,” . ,”

” ?,” . ,”

35)

where µ  is the fluid flux inside hydraulic fractures: ±l”

µ =

²8] oH

(A-36)

Accordingly, the solution of Eq. (A-33) is (Zhang and Yang 2014; Luo and Tang, 2015):

∗ ∗ 4:; = ∑¶ :µg6 ∗ 4 "µg6  − :µ ∗ 4 "µ  − µ06.7 

where:

±:l,

∗ :µg6 =

∗ :µ =

±:l,

∗ :µ =

and:

z

z

/

/

.,”ts

.,”

4: "  =

∗ :µ

(A-37)

(A-38) (A-39)

±:B,” z”

6

*?,” '”

(A-40) >

'*?, K©= ¡J·K©= 

(A-41)

Appendix-B  = 1)

= ;=

B8 .=

(B-

¢

´ =

6W ¾ù ²8]
ø

(B-

4 =

= ∆

(B-

2)

= ∆

6²8] W

3)

 =

6²8] W = ∆

(B-4)

 = B8 9= .

(B-

U =

7.777ô½õ= ¡

(B-

5)



B- 1

¦ ∅Y¢ = .=

6)

 =

8)

9)

10)

(B-

.1

.= 9

.=

$ = 11)

.

.=

" = $ =

(B-

.=

7)

" =

.=

91

.=

(B-

(B-

(B-

Table-1: data for analytical validation, Ozkan et al. 2010. Gas specific gravity, 𝜸𝒈

𝟎. 𝟓𝟓

Molecular weight, 𝑴𝒘

𝟏𝟔

Initial reservoir pressure, 𝑷𝒊

𝟐𝟑𝟎𝟎 𝒑𝒔𝒊

Reservoir temperature, 𝑻

𝟏𝟎𝟗 °𝑭

Formation thickness, 𝒉

𝟐𝟓𝟎 𝒇𝒕

Wellbore radius, 𝒓𝒘

𝟎. 𝟐𝟓 𝒇𝒕

Reservoir boundary, 𝒙𝒆

𝟓𝟎𝟎 𝒇𝒕

Half-distance between hydraulic fractures, 𝒚𝒆

𝟐𝟓𝟎 𝒇𝒕

Viscosity, 𝝁

𝟎. 𝟎𝟏𝟒𝟖 𝒄𝒑

Matrix porosity, ∅𝒖𝒔𝒓𝒗

𝟎. 𝟎𝟓

Matrix compressibility, 𝒄𝒕

𝟎. 𝟎𝟎𝟎𝟗 𝒑𝒔𝒊−𝟏

Matrix permeability, 𝒌𝒖𝒔𝒓𝒗

𝟏 ∗ 𝟏𝟎−𝟖 𝒎𝒅

Natural fracture permeability, 𝒌𝒇

𝟐𝟎𝟎𝟎 𝒎𝒅

Hydraulic fracture permeability, 𝒌𝑭

𝟏 ∗ 𝟏𝟎𝟓 𝒎𝒅

Natural fracture porosity, ∅𝒇

𝟎. 𝟒𝟓

Hydraulic fracture, ∅𝑭

𝟎. 𝟑𝟖 𝒎𝒅

Natural fracture compressibility, 𝒄𝒕

𝟎. 𝟎𝟎𝟎𝟗 𝒑𝒔𝒊−𝟏

Hydraulics fracture compressibility, 𝒄𝒕

𝟎. 𝟎𝟎𝟎𝟗 𝒑𝒔𝒊−𝟏

Hydraulic fracture half-length, 𝒙𝑭

𝟐𝟓𝟎 𝒇𝒕

Hydraulic fracture width, 𝒘𝑭

𝟎. 𝟏𝟓 𝒇𝒕

Number of hydraulic fractures

𝟐𝟎

Table-2: data for numerical validation. Gas specific gravity, 𝜸𝒈

𝟎. 𝟕

Molecular weight, 𝑴𝒘

𝟐𝟐

Initial reservoir pressure, 𝑷𝒊

𝟒𝟎𝟎𝟎 𝒑𝒔𝒊

Reservoir temperature, 𝑻

𝟏𝟎𝟎 °𝑭

Formation thickness, 𝒉

𝟐𝟓 𝒇𝒕

Wellbore radius, 𝒓𝒘

𝟎. 𝟐𝟓 𝒇𝒕

Reservoir boundary, 𝒙𝒆

𝟏𝟎𝟎𝟎 𝒇𝒕

Half-distance between hydraulic fractures, 𝒚𝒆

𝟐𝟓𝟎 𝒇𝒕

Viscosity, 𝝁

𝟎. 𝟎𝟏 𝒄𝒑

Matrix porosity, ∅𝒖𝒔𝒓𝒗

𝟎. 𝟎𝟓

Matrix compressibility, 𝒄𝒕

𝟎. 𝟎𝟎𝟎𝟗 𝒑𝒔𝒊−𝟏

Matrix permeability, 𝒌𝒖𝒔𝒓𝒗

𝟏 ∗ 𝟏𝟎−𝟖 𝒎𝒅

Natural fracture permeability, 𝒌𝒇

𝟐𝟎𝟎𝟎 𝒎𝒅

Hydraulic fracture permeability, 𝒌𝑭

𝟏 ∗ 𝟏𝟎𝟓 𝒎𝒅

Natural fracture porosity, ∅𝒇

𝟎. 𝟒𝟓

Hydraulic fracture, ∅𝑭

𝟎. 𝟑𝟖 𝒎𝒅

Natural fracture compressibility, 𝒄𝒕

𝟎. 𝟎𝟎𝟎𝟗 𝒑𝒔𝒊−𝟏

Hydraulics fracture compressibility, 𝒄𝒕

𝟎. 𝟎𝟎𝟎𝟗 𝒑𝒔𝒊−𝟏

Hydraulic fracture half-length, 𝒙𝑭

𝟓𝟎𝟎 𝒇𝒕

Hydraulic fracture width, 𝒘𝑭

𝟎. 𝟏 𝒇𝒕

Number of hydraulic fractures

𝟕

Table-3: Input data to the proposed models. Number of segments, 𝑵

𝟏𝟎𝟎𝟎

Matrix-fractures fluid flux, 𝒒𝒎𝑫

𝟎. 𝟎𝟎𝟏

Hydraulic fracture Reynolds number, 𝑹𝒆𝒉𝒇

𝟏𝟎𝟎𝟎

Hydraulic fracture conductivity, 𝑭𝑪𝑫

𝟒. 𝟎

Minimum relative permeability, 𝒌𝒎𝒓

𝟎. 𝟓

Hydraulic fracture width, 𝒘𝑫

𝟎. 𝟎𝟎𝟎𝟐

Reservoir boundary, 𝒙𝒆𝑫

𝟐. 𝟎

Half-distance between hydraulic fractures, 𝒚𝒆𝑫

𝟎. 𝟓

Figure 1: Stimulated and un-stimulated reservoir drainage areas.

Figure 2: Matrix blocks in stimulated and un-stimulated reservoir volumes.

1.0E+00 P=1000 psi 1.0E-01

P=3000 psi P=5000 psi

1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07 1.0E-08 1.0E-08

1.0E-07

1.0E-06

1.0E-05

1.0E-04

1.0E-03

1.0E-02

Figure 3: Comparison of apparent and original matrix permeability.

Figure 4: The concept of the characteristic length.

1.0E-01

1.0E+00

1.2

Series1 Series2 Series3 Series4 Series5 Series6 Series7 Series8 Series9 Series10

1

0.8

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

1

Figure 5: Hydraulic fracture permeability coefficient for small characteristic length.

1.2

1

0.8

Series1 Series2 Series3 Series4 Series5 Series6 Series7 Series8 Series9 Series10

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Figure 6: Hydraulic fracture permeability coefficient for large characteristic length.

Figure 7: Fluid flow in the un-stimulated reservoir volume.

Figure 8: Fluid flow in the stimulated reservoir volume.

Figure 9: Discretizing hydraulic fractures to several segments.

1.0E+05

1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00 Series2

1.0E-01

Series3 Series4

1.0E-02

Series5

Series6

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 10: Pressure behavior of fractured reservoirs for different fluid flux.

1.0E+05

1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00 Series2

1.0E-01

Series3 Series4

1.0E-02

Series5

Series6

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

Figure 11: Pressure behavior of fractured reservoirs for different hydraulic fracture Reynold number.

1.0E+08

1.0E+05

1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00

1.0E-01

1.0E-02

Series2 Series21

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 12: Pressure behavior for different hydraulic fracture Reynolds number and fluid flux.

1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00

Series2

1.0E-01

Series3 Series4

1.0E-02 Series5

Series6

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 13: Pressure behavior of fractured reservoirs for different minimum relative permeability of hydraulic fractures.

1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00

1.0E-01

1.0E-02

Series2 Series3

1.0E-03 Series4 Series5

1.0E-04

Series6

1.0E-05 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 14: Pressure behavior of fractured reservoirs for different hydraulic fracture conductivity. 1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00

1.0E-01

Series2 Series3 Series4

1.0E-02

Series5

Series6

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

Figure 15: Pressure behavior of fractured reservoirs for different hydraulic fracture width.

1.0E+08

1.0E+04

1.0E+03

1.0E+02

1.0E+01

1.0E+00

1.0E-01

1.0E-02

1.0E-03

Series2

1.0E-04

Series5

1.0E-05 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E+08

Figure 16: Pressure behavior for different hydraulic fracture conductivities and fluid fluxes. 1.0E+4

1.0E+5

1.0E+3

1.0E+4 1.0E+3

1.0E+2

1.0E+2 1.0E+1 1.0E+1

1.0E+0 1.0E+0 1.0E-1 1.0E-1 1.0E-2 1.0E-2

Series2

1.0E-3

Series3

1.0E-3

Series4

1.0E-4

1.0E-4

Series5 Series6

1.0E-5 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

Figure 17: Flow rate and cumulative production behaviors for different fluid fluxes.

1.0E-5 1.0E+08

1.0E+5

1.0E+5

1.0E+4

1.0E+4

1.0E+3

1.0E+3

1.0E+2

1.0E+2

1.0E+1

1.0E+1

1.0E+0

1.0E+0

1.0E-1

1.0E-1

1.0E-2

1.0E-2 Series2

1.0E-3

1.0E-3

Series3

1.0E-4

1.0E-4

Series4

Series5

1.0E-5

1.0E-5

Series6

1.0E-6 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E-6 1.0E+08

Figure 18: Flow rate and cumulative production behaviors for different fluid fluxes.

1.0E+03

1.0E+03 Series2 Series3

1.0E+02

Series4

1.0E+02

Series5 Series6

1.0E+01

1.0E+01

1.0E+00

1.0E+00

1.0E-01

1.0E-01

1.0E-02

1.0E-02

1.0E-03 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

Figure 19: Productivity index for different fluid flux and strong non-Darcy flow conditions.

1.0E-03 1.0E+08

1.0E+04

1.0E+04 Series2 Series3

1.0E+03

1.0E+03

Series4 Series5 Series6

1.0E+02

1.0E+02

1.0E+01

1.0E+01

1.0E+00

1.0E+00

1.0E-01

1.0E-01

1.0E-02 1.0E-08

1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

1.0E+04

1.0E+06

1.0E-02 1.0E+08

Figure 20: Productivity index for different fluid flux and weak non-Darcy flow conditions.

Figure 21: Comparison for the analytical solutions using the proposed models and Ozkan et al. 2010.

Figure 22: Reservoir and hydraulic fractures configuration.

1.0E+07 Series1

Series2

1.0E+06

1.0E+05

1.0E+04

1.0E+03

1.0E+02

1.0E+01 1.0E-06

1.0E-04

1.0E-02

1.0E+00

1.0E+02

Figure 23: Pressure behavior calculated by the proposed models and CMG.

1.0E+04

1.0E+06

Highlights: 1- Reservoir fluid flux from a structurally complicated matrix to non-uniformly propagated hydraulic fractures is studied I this paper. 2- Semi-analytical multilinear flow regimes model is developed by coupling matrix Darcy flow and hydraulic fracture non-Darcy flow. 3- New Reynolds number in hydraulic fractures and minimum hydraulic fracture relative permeability to Darcy flow are presented. 4- Fluid flux from SRV matrix to hydraulic fractures may have significant impact on reservoir performance during the entire production life.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: