Transient responses of multifractured systems with discrete secondary fractures in unconventional reservoirs

Accepted Manuscript Transient responses of multifractured systems with discrete secondary fractures in unconventional reservoirs Linsong Cheng, Pin Jia, Zhenhua Rui, Shijun Huang, Yongchao Xue PII:

S1875-5100(17)30092-6

DOI:

10.1016/j.jngse.2017.02.039

Reference:

JNGSE 2092

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 3 August 2016 Revised Date:

11 January 2017

Accepted Date: 22 February 2017

Please cite this article as: Cheng, L., Jia, P., Rui, Z., Huang, S., Xue, Y., Transient responses of multifractured systems with discrete secondary fractures in unconventional reservoirs, Journal of Natural Gas Science & Engineering (2017), doi: 10.1016/j.jngse.2017.02.039. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Transient responses of multifractured systems with discrete secondary fractures in unconventional

ACCEPTED MANUSCRIPT reservoirs Linsong Cheng a, Pin Jia a, , Zhenhua Rui b, , Shijun Huang a, Yongchao Xue a *

*

a

College of Petroleum Engineering, China University of Petroleum, Beijing, China

b

Independent Project Analysis, Inc.

Road, Changping, Beijing, China. E-mail address: [email protected] (P. Jia)

*Corresponding author at: Independent Project Analysis, Inc.

M AN U

E-mail address: [email protected] (Z. Rui)

SC

RI PT

*Corresponding author at: College of Petroleum Engineering, China University of Petroleum, 18 Fuxue

AC C

EP

TE D

.

1 / 33

Abstract

ACCEPTED This paper provides a semi-analytical model MANUSCRIPT to analyze the transient responses of multifractured systems with discrete secondary fractures in unconventional reservoirs. The model idealizes the stimulated reservoir volume as two regions: a primary hydraulic-fracture and discrete secondary fractures. Transient pressure and derivative responses are obtained by coupling an analytical model for matrix flow and a numerical model for fracture flow. This approach can consider the fracture characteristics in detail. For

RI PT

example, the discrete secondary fractures may interconnect the primary hydraulic-fracture with different inclination, nonuniform fracture spacing and varying conductivity. The semi-analytical model was verified in comparison with the results of a fully numerical simulation model.

SC

Transient flow behavior of the system was analyzed in detail. The study has shown that a multifractured system with discrete secondary fractures may exhibit four flow regimes: bilinear flow, fluid feed flow, matrix

M AN U

linear flow and pseudosteady-state flow. In the fluid feed flow period, the secondary fractures act as supply sources to feed the primary hydraulic-fracture to generate a dip on the type curves. Local solutions for the matrix linear flow are similar to that for a finite conductivity vertical fracture and can be correlated with the ratio of the total length of the secondary fractures to that of the primary hydraulic-fracture and an intercept term. With the increasing of number of secondary fracture, the ending time of the bilinear flow will be

TE D

advanced and the fluid feed flow period will come earlier. The dimensionless secondary fracture length determines the slope of the straight line of the matrix linear flow on square-root-dimensionless-time plot. As the secondary fracture inclination decreases, the duration of the matrix linear flow period is shortened and

EP

the beginning time of the pseudosteady-state flow period is delayed. In addition, the secondary fracture conductivity affects the depth of the dip and the beginning time of the matrix linear flow period. However,

AC C

the influence is not appreciable on the log/log graph.

Keywords: transient responses; multifractured system; secondary fracture; semi-analytical model; unconventional reservoirs

2 / 33

1. Introduction

MANUSCRIPT The technology of horizontal ACCEPTED wells with multi-stage hydraulically fracturing enables commercial development of hydrocarbons from unconventional reservoirs, such as shale gas reservoirs. Application of these technologies generates stimulated reservoir volumes (SRVs) around the well, which form a multifractured system. It has been postulated that, in the multifractured system, along with the primary hydraulic-fractures, the branched-fracture networks, may be activated during the stimulation process,

RI PT

creating secondary fractures (Fisher et al., 2004; He, et al., 2016). These secondary fractures can be the result of reactivation of healed natural fractures, for example (Weng et al., 2014). To provide practical evaluation of well performance and asymptotic approximations about potential flow regimes, many researchers simplified

SC

the system to obtain analytical/semi-analytical solutions to analyze the transient flow behavior.

The most common approaches to model the flow process are based on the assumption that each SRV

M AN U

can be idealized as two regions: a primary hydraulic-fracture which is essentially a planar fracture, and an enhanced region, which can be modeled by introducing a high permeability zone or a dual-porosity region (see diagrams in Figs. 1a and 1b). Using a representative element of the system (boxed in the red dashed rectangle in Fig. 1b), the trilinear flow model (Brown et al., 2011; Ozkan et al., 2011) and its variations (Brohi et al., 2011; Stalgorova and Mattar, 2012; Stalgorova and Mattar, 2013; Xu et al., 2013; Sureshjani

TE D

and Clarkson, 2015a) are developed by coupling several linear flows in continuous flow regions: the primary hydraulic-fracture, the rectangular enhanced region and the un-stimulated region, to derive analytical pressure/rate solutions. Semi-analytical models are set up to consider the simultaneous effect of multiple

EP

SRVs (Medeiros et al., 2008; Sun et al., 2016) and a circular shape enhanced region (Zhao et al., 2014; Zhao et al., 2015a). These models have provided a better understanding of the performance of the multifractured

AC C

system and offered practical approaches to forecast production and evaluate stimulation effectiveness. Most importantly, local solutions for specific flow regimes are derived to facilitate the interpretation of the production data. However, if there are some discrete secondary fractures with high conductivity, which will control the principle flow pathways in SRV, these fractures should be represented explicitly to capture the flow behavior more rigorously. There are some models taking account of discrete secondary fractures with the assumption that secondary fractures are orthogonal to primary hydraulic-fracture at uniform fracture spacing. Assuming that formation fluids flow from matrix to microfracture (refer to secondary fracture), then to macrofracture (refer 3 / 33

to primary hydraulic-fracture), Al-Ahmadi and Wattenbarger (2011) and Tivayanonda et al. (2012) coupled

MANUSCRIPT the three linear flows to obtain ACCEPTED analytical pressure/rate solutions. Chen et al. (2016) provided a semi-analytical model to analyze the effect of secondary fractures on the flow behavior of multifractured horizontal wells by discretizing the fractures to grids. However, the orthogonal fracture geometry limits above models. In practical, the secondary fractures may interconnect the primary hydraulic-fracture with different inclinations. In addition, although many numerical models can accurately simulate production of

RI PT

fracture system with complex geometry, most of them focus on the long-term production profile (Noorishad and Mehran, 1982; Granet et al., 2001; Karimi-Fard et al., 2003; Z Li and Lee, 2008; Mayerhofer et al. 2010; Cipolla et al. 2011; Mi et al., 2014; Wu et al. 2014; Zhao et al., 2015b; Ou et al., 2016), few studies have

SC

been done about transient responses and flow regimes for a primary hydraulic-fracture with discrete secondary fractures.

M AN U

Therefore, there is still lack of a rigorous model and detailed discussion about the transient responses of multifractured system with discrete secondary fractures. In this paper, we provide a semi-analytical model to the transient flow behavior of such system. An analytical model for matrix flow is dynamically coupled with a numerical model for fracture flow. In this approach, discrete secondary fractures can interconnect primary hydraulic-fracture with different inclinations, nonuniform fracture spacing and varying conductivity.

TE D

Transient responses characteristics of the system are identified and effects of the key parameters of the secondary fractures on the flow behavior are also investigated. 2. Physical model

EP

Fig. 1a shows a horizontal well treated by four-stage hydraulic fracturing. Fracture branching forms stimulated reservoir volumes, which consist of complex fracture systems. To model this system, we assume

AC C

that all the SRVs are located uniformly along the horizontal well with the same properties. This assumption is realistic because it is common field practice to design equally spaced hydraulic fracturing with same properties. As Fig. 1c illustrates, each SRV around the horizontal well is idealized as two regions: a primary hydraulic-fracture and several discrete secondary fractures. Therefore, there is no contrast between matrix permeability of un-stimulated region and that of SRV. Due to the symmetry of the system, we consider a representative element of the system, for example the region specified by the red dashed lines in Fig. 1c. For the unconventional reservoirs with microdarcy matrix permeability, it can be assumed that the contribution of the region beyond the primary hydraulic-fracture is insignificant and can be neglected. Therefore, the system 4 / 33

shown in Fig. 2a can be used for modeling. The validation of this assumption will be discussed later.

MANUSCRIPT Because the matrix permeabilityACCEPTED of un-stimulated region is same with that of SRV, the boundary of SRV (black dashed line in Fig.2a) can be removed. Fig. 2b shows the schematic and dimensions for the suggested physical model. The half-length of the primary hydraulic-fracture on the both sides of the well is equal to Lpf and the size of rectangular drainage area is Lpf×ye. The number of secondary fractures is Msf and these fractures can have different inclination, length and conductivity. However, to facilitate the development of

RI PT

the flow model, we assume that the length of the secondary fractures is same and the conductivity is uniformly distributed. In Fig. 2c, the primary hydraulic-fracture and each secondary fracture are divided into Npf and Nsf segments, respectively. Thus, the fracture system is totally discretized into Nf (Nf=Npf+Nsf·Msf)

SC

fracture segments, which are numbered sequentially. Based on the Star-Delta transformation, we eliminate the intersections caused by the interconnected fracture segments. The application of this transformation is

M AN U

discussed later.

In addition, fully penetrating fractures are assumed and 2D solution is used. And we ignore the flow convergence to the wellbore within the primary hydraulic-fracture. Reservoir fluids only flow through the primary hydraulic-fracture to the wellbore. The porosity and compressibility in the matrix and the fractures

3. Mathematical model

TE D

are dependent of pressure. The system considers a single phase compressible fluid.

As noted above, the flow process can be envisaged as two distinct flows governed by different physics: flow within the matrix, towards the fractures and flow within the fractures, towards the wellbore. We derive

EP

the solutions for the matrix and the fractures, and then couple the two solutions by imposing pressure and flux continuity on the fracture surfaces.

AC C

3.1. Analytical matrix flow model

The basic idea of the mathematical formulation follows that provided by Brown et al. (2011). The difference is that our model accounts for 2D flow in the matrix, instead of 1D flow perpendicular to the primary hydraulic-fracture in Brown et al. (2011). If we consider the primary hydraulic-fracture and secondary fracture segments as line sources, then the flow in the matrix can be regarded as flow process in a closed rectangle, along with Nf continue sources of varying strength. Then, the transient flow equation for the matrix flow can be given by

5 / 33

∂ 2 pD ∂ 2 pD ∂pD + = ∂xD2 ∂y D2 ∂t D

(1)

pD ( xD , y D , t D → 0 ) = 0

(2)

∂pD = 0, xD = 0, 0 ≤ yD ≤ yeD ∂xD

(3)

ACCEPTED MANUSCRIPT Initial condition

∂pD = 0, xD = 1,0 ≤ yD ≤ yeD ∂xD ∂pD = 0,0 ≤ xD ≤ 1, yD = 0 ∂y D

RI PT

The closed boundary conditions

SC

∂pD = 0,0 ≤ xD ≤ 1, yD = yeD ∂y D

(4)

(5)

(6)

M AN U

And Eq. 1 is also subject to the boundary condition at the surface of Nf fracture segments corresponding to the continue withdrawal of different amount of fluid. Definitions for dimensionless parameters are given in Appendix A.

To solve Eq. 1, we first take the Laplace transform to Eqs. 1 through 6, and then use line source

TE D

functions to obtain the analytical solution. Source functions are powerful tools for the flow model to obtain transient responses of fractured wells. Ozkan and Raghavan (1991) derived line source solutions for different boundary conditions in Laplace domain. Thus, the dimensionless pressure caused by Nf fracture segments in

EP

a closed rectangular drainage region can be expressed as pDM = ∑ j =f1 sqDj pDM , j N

(7)

AC C

A detailed derivation of Eq. 7 is presented in Appendix B. In Eq. 7, qDj is the flow rate entering the fracture segment j from matrix. s is the Laplace tranform variable. pDM , j is line source function of segment j. Then, the dimensionless pressure of the segment i can be obtained by pDi = ∑ j =f1 sqDj pDi , j N

for i =1…Nf

(8)

By use of Eq. 8, we can build a linear system of Nf equations which contains the pressure and flow rate for all the segments

A⋅ qD = pD 6 / 33

(9)

spD1,2 L spD1,Nf   spD11,  qD1   pD1    q  p  sp sp sp L  ACCEPTED D21 , D2,2 D2,N f  MANUSCRIPT  D2   D2    A =  spD31, spD3,2 L spD3,Nf  ,qD =  qD3  , pD =  pD3         M   M  O   q  p  spDNf ,1 spDNf ,2 L spDNf ,Nf   DNf   DNf   

We can find that Eq. 9 contains the information of the matrix and the geometry properties of the primary hydraulic-fracture and secondary fractures.

RI PT

3.2. Numerical model for fracture flow

The primary hydraulic-fracture and secondary fractures are represented by homogeneous, finite, slab, porous media. We first present the flow equation for the fractures; and then discuss the approach to deal with

flow can be described by ∂2 plfD

Initial condition

−

π

⋅ qlfD (εD,tD) =

clfD

1 ∂plfD ηlfD ∂tD

M AN U

∂ε

2 D

SC

the flow in the intersections. For the primary hydraulic-fracture and each secondary fracture, the transient

plfD (ε D , t D = 0) = 0

(10)

(11)

The closed boudanry of fracture tips can be described by

TE D

∂plfD ∂ε D

=0

(12)

ε D = Tips

If we consider the constant production rate condition of the horizontal well, thus

EP

∂plfD ∂εD

εD =Wellbore

=−

π

(13)

clfD

AC C

In Eqs. 10 through 13, we introduce the symbol l to derive a general formulation for the fracture flow. For the flow equation of the primary hydraulic-fracture, l=p , and for the secondary fractures，l=s. qlfD has the same physical meaning as qDj in real domain. The symbol ε represents the coordinate along the fractures. In the Laplace domain, Eqs. 10 through 13 have the following form

7 / 33

∂2 p  lfD − π ⋅ qlfD (εD, s) = s plfD ACCEPTED  ∂εD2 cMANUSCRIPT ηlfD lfD  ∂plfD =0   ∂εD εD =Tips  1π ∂plfD =−  ∂εD s clfD εD =Wellbore 

(14)

Eq. 14 is discretized to yield the following finite-difference formulation

RI PT

TDi−1,i plfDi−1(s) −(TDi−1,i +TDi+1,i + βDi ) plfDi (s) +TDi+1,i plfDi+1(s) − qlfDi (s) = 0

(15)

where, TD is the dimensionless transmissibility of the adjacent segments. TDi-1,i and βDi are defined as

(16)

M AN U

SC

  clfD  λ λ TDi −1.i = Di −1 Di , λDi =   ∆LlfD / 2  λDi −1 + λDi  i    clfD ∆ LlfD     β Di = s  η lfD  i 

The dimensionless definitions in Eqs. 10 through 16 are all defined in Appendix A. For the flow in the intersections, the situation has been comprehensively discussed in our form paper (Jia et al. 2015). Here, we give a brief introduction about the approach. In the fracture flow model, the pressure of the intersection cells created by fracture interconnecting must be obtained to determine flow

TE D

condition between the interconnected fractures. However, the intersections cannot be represented by line sources in the matrix flow model because of their significant small size compared with that of the fracture segments. The intersections in the numerical model for fracture flow should be removed. Here, we adopt the

EP

Star-Delta transformation to solve the problem. After the transformation, the intersection cells are eliminated and the interconnected fracture segments will be direct-connected. And the transmissibility between the

AC C

direct-connected segments is expressed as TDi∗ , j =

TDi ,0TDj ,0

∑

4

(17)

T k =1 Dk ,0

where, TDk,0 is the transmissibility between the intersection and segment k and TDk,0=γDk. By use of Eq. 17, we modify the Eq. 15 for the interconnected segments to eliminate the unknowns of the intersections. Based on Eqs. 15 through 17, we combine the NS equations of all fracture segments shown in Fig. 2c to obtain a linear system containing the pressure and flow rate of all the fracture segment and the wellbore pressure 8 / 33

T ⋅ plfD = qlfD − pwD ⋅ b

(18)

ACCEPTED MANUSCRIPT

T

where, b = TD1, w ,0,L,0  . TD1, w is the dimensionless transmissibility between the first fracture segment and the wellbore, and TD1,w=γD1. In addition, constant production rate condition can be expressed as

π

TD1,w ⋅  ppfD1(s) − pwD (s) =−

(19)

s

RI PT

3.3. Solution for transient responses

Transient responses are solved by coupling the two models for matrix flow and fracture flow. In our semi-analytical model, the fracture segments are considered as planar sources. Withdrawal of fluid from

SC

matrix into segments causes pressure transients in the system. Hence, the dimensionless pressure plfD ( εD ) and flux density qlfD ( εD ) in the fracture model must equal the dimensionless pressure drop pD ( xD, yD ) and

M AN U

flux density qD ( xD, yD ) on the plane source of the matrix model. Thus, by of the continuity of pressure and flow rate on the surface of the fracture, the following conditions must be held on the fracture surface.

plfD = pD

(20)

qlfD = qD

(21)

flow rate and wellbore pressure

TE D

Combining Eqs. 9, 18, and 19 generates a linear system to be solved for the fracture segment pressures,

EP

A  − I  oT 

T

−I o   qlfD   o  T b  ⋅  plfD  =  o b T TD1, w   pwD   −π

   s 

(22)

T

AC C

where, qlfD = qlfD1,L, qlfDNf  , plfD =  plfD1,L, plfDNf  , o is Nf×1 null vector and I is Nf×Nf identical matrix. Eq.     22 can be solved by Gauss elimination method. And then taking Stehfest numerical inversion algorithm (Stehfest 1970) to the unknowns will yield solutions in real domain. Following the idea of Stalgorova and Mattar (2012), we provide some possible model modifications to extend the model to gas flow. For gas flow, the dimensionless pressure should be defined in terms of gas pseudopressure. However, for the case of gas, the diffusivity term ηm = km / (φmµgcmt ) is varying with pressure, so the flow governing equation becomes nonlinear. The common way to deal with the problem is to use pseudo-time (Anderson and Mattar, 2005) to linearize the equation. 9 / 33

Finally, in the aforementioned sections, we adopted 1D solution to model the flow in the primary

ACCEPTED hydraulic-fracture and neglect the flow convergence MANUSCRIPT to the wellbore within the fracture. The effect has been taken as a skin factor by Mukherjee and Economides (1991)

sc =

km h  h π −   ln k pf w pf  2rw 2 

(23)

This skin factor, along with the wellbore storage coefficient (dimensionless CD = C / 2πφm cmt hL2pf ), can

(1949)

spwD + sc s + CD s 2 ( spwD + sc )

(24)

SC

pwD ,storage =

RI PT

be incorporated into the solutions by use of the following equation derived by Van Everdingen and Hurst

4. Results and Discussion

M AN U

Before proceeding to the validation of the model, we first take a sensitivity analysis to obtain accurate results. From Eq. 22, wen can see that for the solution in Laplace domain, the factor affecting the accuracy of resluts is only the fracture gridding issue. Hence, we primarily conduct sensitivity anlaysis of the number of fracture segments. The tested range of the number of segments for each fracture in the sensitivity analysis is from 5 to 30. For the caes with lowest conductivity (cpfD<10 and csfD<10), we also test up to 40 segments.

TE D

The results show that the solutions do not change appreciably when half of the primary hydraulic-fracture is discretized into 15 segments and each secondary fracture is divided into 10 semgents. Therefor, in all cases studied, we set the value of Npf and Nsf to 15 and 10, respectivly. Then, we verify the semi-analytical model,

EP

and investigate the characteristics of transient flow behavoir and the parameters effect analysis, including the inclintion, number, length, conductivity and geometry of the secondary fractures.

AC C

4.1. Model validation

A fully numerical simulation model was built using discrete fracture model presented by Mi et al. (2014) to validate the semi-analytical model. By use of finite element scheme, Mi et al. (2014) developed a numerical production simulation model for hydraulic fracture interconnected with discrete natural fractures in unconventional reservoirs. In the fully numerical model, the matrix domain is discretized by Delaunay triangulation algorithm. The fracture width, which is negligible compared to the matrix block size, is taken into account in the flux calculation but not in the geometrical representation. Therefore, the model adopts linear elements for the fractures. The cell dimensions of the matrix around the fractures are refined to capture 10 / 33

the characteristics of the fractures more rigorously. Fig. 3 shows a top view of the grid configuration for the

ACCEPTED system in Fig. 2. The fractures are indicated by the MANUSCRIPT black bold linear elements with a width of 0.01m. The well is completed in the primary hydraulic-fracture element, which is located in the left-bottom of the global grid system. The basic parameters are listed in Table 1. In the semi-analytical model, we adopt 55 (Nf=55) fracture segments and 10 points in each log cycle to simulate the system.

Fig. 4 shows that except at very early times, the transient pressures of the semi-analytical model are in

RI PT

good agreement with that of the fully numerical model. This may be caused by different methods to tackle with the radial flow within the primary-hydraulic fracture. In the fully numerical model, the radial flow is dealed with well flow model. However, the semi-analytical model uses the equivalent skin factor. The two

SC

different approaches generate different skin factors for the flow chocking in the primary hydraulic-fracture. The transient derivatives agree well because derivatives of pressures with respect to time are independent of

M AN U

the skin factor.

One of the assumptions made for the model is that the contribution of the un-stimulated region beyond the primary hydraulic-fracture is negligible. Contribution from this region depends on many parameters, mainly including the matrix permeability (km), primary hydraulic-fracture length (Lpf), secondary fracture length (Lsf), number of secondary fracture (Msf) and the representative element size in y-direction (ye). To

TE D

validate this assumption, we compare the results of the model in Fig. 2 (refered to no flow boundary case) with that of the element specified by red dashed line in Fig. 1c (refered to normal flow boundary case). The tested parameters are dimensionless secondary fracture length (LsfD), dimensionless representative element

EP

size in y-direction (yeD) and number of secondary fracture (Msf). And we set the value of matrix permeability to 4×10-5×10-3µm2. Some parameters are refered to Table 1.

AC C

Fig. 5a shows the effect of LsfD. We can find that for large values of LsfD, the plot of no flow boundary case has an excellent fit with that of normal flow boundary case. But with decreasing of LsfD, the curve of normal flow boundary case shifts to the right; that is, the un-stimulated region beyond primary hydraulic-fracture starts to feed into the SRV and inter-fracture region. This is mainly because as LsfD decreases, the distance from the secondary fracture tips to the boundary (yD=yeD) becomes larger and the ability of the secondary fracture to convey formation fluid from the matrix near this boundary will be weak. Therefore when the pressure pulse reaches the boundary (yD=yeD), the transients from the un-stimulated region beyond the primary hydraulic-fracture will be significant enough to prevent the system of a 11 / 33

pseudosteady-state flow period. The influence of yeD is illustrated in Fig. 5b. when yeD equals to 0.9, the accuracy of the assumption of the noACCEPTED flow boundary MANUSCRIPT decreases. The reasons may be that as yeD increases, the contribution of the un-stimulated region beyond the primary hydraulic-fracture is increased. The above two situations are similar to the cases discussed by Sureshjani and Clarkson (2015a). In Fig .5c, the differences between the plots of no flow boundary case and normal flow boundary case will appear when Msf=1. With decreasing of number of secondary fracture, the fracture system in Fig.2 will behave as a single fracture

RI PT

perpendicular to wellbore. Because there is no contrast of matrix permeability throughout the region, the inter-fracture region and un-simulated region beyond the primary hydraulic-fracture will supply production simultaneously.

SC

Although the accuracy of assumption of neglecting the contribution of the un-stimulated region beyond the primary hydraulic-fracture decreases in some tests, the largest starting time of the derivation of the plots

M AN U

is no more than 0.3 (tD=0.3). Based on the parameters listed in Table 1, this dimensionless time in real case is 15.5 years, which is enough long for production analysis of a typical well in unconventional reservoirs. Hence, the assumption of neglecting the contribution of the un-stimulated region beyond the primary hydraulic-fracture is practically valid.

4.2. Transient responses characteristics

TE D

The matrix permeability, km, is set to 7×10-5×10-3µm2. And the fractures, including primary-hydraulic fracture and secondary fractures are totally discretized into 55 segments. Some of the parameters are referred to Table 1. For the system in Fig. 2, the dimensionless transient derivatives of the system are shown in Fig. 6.

EP

And Fig. 7 illustrates the pressure distribution pictures for the region in each defined flow period. Analysis of Fig. 6 indicates that the system may exhibit the following four flow periods.

AC C

Bilinear Flow. It is caused by simultaneous linear flow in the primary hydraulic-fracture and the secondary fractures, as shown in Fig. 8a. At this stage, the production is mainly supported by the fracture system. The dimensionless derivatives show a straight line with a one-quarter-slope on the type curves.

Fluid Feed Flow. In this period, a dip on the curve of the transient derivative responses will be identified. This can be inspired by the origin of the dual-porosity signature. As shown in Fig. 8b, the ability of the primary hydraulic-fracture support to production diminishes gradually and the secondary fractures will act as supply sources to feed the primary hydraulic-fracture. Thus, the fluid feed takes effect and becomes stronger. Essentially, most of the production fluid comes from the matrix near the secondary fractures. This 12 / 33

situation is similar to the cases discussed in the literatures (Al-Ahmadi and Wattenbarger 2011; Johns et al.

ACCEPTED MANUSCRIPT

2013; Chen et al. 2016).

Matrix Linear Flow. It is called linear flow because fluid from the matrix flows into the fracture system, including the primary hydraulic-fracture and secondary fractures, in a direction perpendicular to the fracture surface, as shown in Fig. 8c. Each fracture segment haves its own drainage area. In the log/log graph, the dimensionless transient responses are indicated by half-slope straight lines. This flow behavior is

RI PT

essentially similar with the linear flow of a finite conductivity vertical fracture with constant rate condition provided by Sureshjani and Clarkson (2015b). At this stage, the flow resistance due to the finite conductivity fracture can be expressed as a constant skin in linear flow solution.

Sureshjani and Clarkson (2015b) is given by

(25)

M AN U

pwD = π t D + b

SC

The pressure solution for a finite conductivity vertical fracture during linear flow period presented by

where, parameter b is the intercept of the square-root-dimensionless-time plot, which has a positive value. This value is a function of dimensionless fracture conductivity, FCD. The definition of FCD can be refered to the work of Sureshjani and Clarkson (2015b).

In this paper, the system behaves as a finite conductivity vertical fracture with the total length of the

TE D

fractures, including the primary hydraulic-fracture and the secondary fractures. Following the idea of Sureshjani and Clarkson (2015b), we can derive linear flow solution for the system in Fig.2

EP

pwD = π

kmt

φm cmt µ ( L pf + M sf Lsf )

2

+ b'

(26)

AC C

where, parameter b' is the skin factor for flow resistance in the fracture system. For a certain system, b' is a constant and is a function of the dimensionless primary hydraulic-fracture conductivity, dimensionless secondary fracture conductivity, the ratio of the total length of the secondary fractures to that of the primary

hydraulic-fracture and also the geometry of the fracture system. It is noted that the resistance due to flow convergence to the wellbore within the primary hydraulic-fracture is also contributes to b' . However, the expression of b' is out of the discussion scope in this paper. Then, Eq. 26 can be rewritten as

pwD =

L pf Lpf + M sf Lsf 13 / 33

π

km t

φm cmt µ L2pf

+ b'

(27)

If we define the ratio of the total length of the secondary fractures to that of the primary

MANUSCRIPT hydraulic-fracture as Rf and use the ACCEPTED definition of dimensionless time; then, the solution for the matrix linear flow can be given by

pwD =

1 1 + Rf

π tD + b'

(28)

or

qw B µ q Bµ ' 1 t1 2 + w b 12 2π km h 2 π (1 + R f ) hL pf (φm cmt µ km ) 1

RI PT

pw = pi −

(29)

A comparison of Eq. 25 and Eq. 28 indicates that on the square-root-dimensionless-time plot, the differences between these two curves are clearly the slop and the intercept. Eq. 29 also suggests that the total

SC

length (Msf · Lsf) of the secondary fractures can be interpreted, if the primary hydraulic-fracture length (Lpf)is known.

M AN U

Pseudosteady-state Flow. As shown in Fig. 8d, once the closed boundaries begin to influence the well responses, a bend becomes evident and the transient response curves follow a unit-slope straight line reflecting boundary-dominated flow in the log/log plot. Pseudosteady-state flow has also been discussed by many researchers (Brown et al. 2011; Stalgorova and Mattar 2013; Sureshjani and Clarkson 2015a).

TE D

As shown in Fig. 9, wellbore storage may mask the characteristics of bilinear flow, whereas other flow periods are identical to those shown in Fig. 7. If the wellbore storage coefficient increases, the afterflow period will even mask fluid feed flow. We also compare the results in Fig. 7 with that of Chen et al. (2016).

EP

All the input parameters are the same, excluding the secondary fracture location. In the model of Chen et al. (2016), the secondary fractures are orthogonal with the primary hydraulic-fracture at uniform fracture

AC C

spacing. As shown in Fig. 10, the deviances of the two curves are mainly at middle times, which is between the linear flow period and pseudosteady-state flow period. The reason may be that in the first three flow periods, each fracture segment behaves independently of others. Hence, the fracture location has little influence on the flow behavior. After the end of the linear flow period, the pressure interference of the fractures takes effect. Different fracture geometries generate different interference situations, which lead to the contrasts of the transient responses of the two curves. Finally, the two curves merge because the characteristics of the pseudosteady-state flow are mainly determined by the properties of matrix and the representative element size. Also, it can be found that compared with the parallel secondary fractures in Chen et al. (2016), the interconnected secondary fractures make the interference come earlier and thus, result in a 14 / 33

shorter duration of the linear flow period. In addition, we validate linear ACCEPTED flow solution Eq.MANUSCRIPT (28) by plotting it against the semi-analytical model. We consider three cases with different dimensionless secondary fracture conductivity and secondary fracture length, which are referred to as Cases 1 to 3 listed in Table 2. The rest of the input data for these three cases are refered to Table 1. By use of the semi-analytical model, we first generate a synthetic production profile of

pwD =

1 1 + Rf

π tD

RI PT

Case 1. Because the value of b' is not determined in the paper, we use the following equation (29)

which is the linear solution for the fracture system in Fig.2 with infinite conductivity, to validate the slope of

SC

Eq. (28). Fig. 11 shows the comparison results in the square-root-dimensionless-time plot for Case 1 and Eq. (29). It can be found that the synthetic data are parallel to the line corresponding to the fracture system with

M AN U

infinite conductivity. The slope of the straight line of the synthetic data is 0.59, which can be obtained from the figure or Eq. (28). Fig. 12 illustrates the synthetic data for all the three cases. By comparing the data of Case 1 and Case 2, we can find that for a given fracture system, as the dimensionless conductivity of secondary fracture decreases, the intercept of plotted data in the square-root-dimensionless-time plot increases, while the slope of the straight line remains unchanged. In the plot of Case 3, a smaller

TE D

dimensionless secondary fracture length generates a larger slope of the straight line and a larger intercept, compared with that of Case 1. These observations can be reasonably interpreted by Eq. (28). The reasons are that the slope of the straight line is only inverse to Rf and the decreasing of dimensionless fracture

EP

conductivity will cause enlargement of the flow resistance in the fracture system, which results in a larger value of parameter b' . However, as the dimensionless length of secondary fractures decreases, the flow

AC C

resistance will also decrease.

4.3. Effect of secondary fracture number To facilitate the comparison, we use a limited case that the secondary fractures are orthogonal to the primary hydraulic-fracture at uniform fracture spacing. Transient derivatives of the system on the log/log plot for different secondary fracture number are shown in Fig. 13. It can be seen that the secondary fracture number has great influence on the early to middle times, including the bilinear flow, fluid feed flow and linear flow periods. In the linear flow period, a larger secondary fracture number will generate a lower dimensionless derivative; that is, a smaller value of the y-intercept of the straight line of the local solution in 15 / 33

the log/log plots. This is also indicated by Eq. 28 for the reason that on the log/log plot, the y-intercept of the equation is inversely proportional toACCEPTED Rf (the ratio of MANUSCRIPT the total length of the secondary fracture to that of the primary hydraulic-fracture). In addition, as the number of secondary fracture number increases, the duration of the bilinear flow will be shortened and the starting time of the fluid feed flow period will be advanced. It is mainly because the secondary fracture number determines the quantity of the sources to feed the primary hydraulic-fracture. The increased sources will timely maintain the production.

RI PT

4.4. Effect of secondary fracture length To investigate the influence of the secondary fracture length on the transient flow behavior, we use the system in Fig. 2 and data in Table 1, but with different LsfD. The results in Fig. 14 indicate that the effect of

SC

the secondary fracture length is predominantly on the linear flow and the subsequent transient flow. Similar to the effect of the fracture number, with the increasing of the fracture length, Rf increases and the y-intercept

M AN U

of the asymptotic approximations for the linear flow becomes small in the log/log plot. A samller LsfD will result in a larger distance from the secondary fracture tips to the boundary that is parallel to the primary-hydraulic fracture. Hence, if LsfD is small enough, the pseudosteady-state flow will also be influenced and even be delayed.

4.5. Effect of secondary fracture inclination

TE D

Here, we design four cases to investigate the effect of the secondary fracture inclination on the flow behavior of the fracture system in Fig.2. The value of dimensionless secondary fracture length, LsfD, is set to 0.2 to guarantee that the secondary fracture with smallest inclination (10°) is limited in the closed rectangular

EP

area. The dimensionless derivatives are shown in Fig. 15. We can find that with the decreasing of the secondary fracture inclination, the duration of the linear flow period is shortened and the beginning time of

AC C

the pseudosteady-state flow periods is delayed. The reason may be that a smaller fracture inclination will generate smaller fracture spacing, which results in an earlier pressure interference. Thus, the linear flow will end earlier. As the tips of the fracture with smaller fracture inclination will have a larger distance from the boundary (yD =yeD), the time to reach pseudosteady-state flow for the closed zone will be prolonged.

4.6. Effect of secondary fracture conductivity

Fig. 16 shows the transient derivatives on the log/log graph for the system in Fig. 2 with different secondary fracture conductivities. It can be seen that the secondary fracture conductivity mainly impacts the depth of the dip and the starting time of the linear flow period. However, the influence is not appreciable on 16 / 33

the log/log plot. The larger the secondary fracture conductivity is, the deeper the dip gets and the earlier the

MANUSCRIPT linear flow comes. The reason mayACCEPTED be that with the increasing of the csfD, the ability of the secondary fractures to convey the fluid in the matrix around these fractures to maintain the production will be stronger. As a result, the value of derivative responses becomes smaller. These findings are very similar to the work of John et al. (2013) and Chen et al. (2016). In their studies, increasing the secondary fracture conductivity results that the beginning time of the linear flow is advanced and the dip becomes profound.

RI PT

5. Conclusions In this paper, we investigate the transient responses of multifractured system with discrete secondary fractures in unconventional reservoirs by use of a semi-analytical model. The following conclusions can be

SC

made:

1. Existence of discrete secondary fractures introduces a fluid feed flow that the secondary fractures behave

M AN U

as supply sources to feed the primary hydraulic-fracture to maintain the production. This flow period often develops after the bilinear flow and a dip can be characterized on the log/log plot. Thus, a multifractured system with discrete secondary fractures may exhibit four flow regimes: bilinear flow, fluid feed flow, matrix linear flow and pseudosteady-state flow.

2. Local solutions for the matrix linear flow are similar to that for a finite conductivity vertical fracture and

TE D

can be correlated with the ratio of the total length of the secondary fractures to that of the primary hydraulic-fracture and an intercept term. Compared with the orthogonal fracture geometry, the secondary fractures interconnected the primary hydraulic-fracture with different inclinations make the pressure

EP

interference come earlier and thus, the ending time of the matrix linear flow will be advanced. 3. With the increasing of number of secondary fracture, the ending time of the bilinear flow will be

AC C

advanced and the fluid feed flow period will come earlier. The dimensionless secondary fracture length determines the slope of the straight line of the matrix linear flow on square-root-dimensionless-time plot. As the secondary fracture inclination decreases, the duration of the matrix linear flow period is shortened and the beginning time of the pseudosteady-state flow periods is delayed. In addition, the secondary fracture conductivity affects the depth of the dip and the beginning time of the matrix linear flow period. However, the influence is not appreciable on the log/log graph.

Acknowledgments 17 / 33

This study was partially supproted by the Major State Basic Research Development Program (973

ACCEPTED MANUSCRIPT Program) (No. 2013CB228000 and No. 2015CB250902). And we also would like to thank financial support of the PetroChina Innovation Foundation (No. 2014D-5006-0206) and the Beijing Natural Science Foundation (No. 3144033).

Nomenclature Liquid formation volume factor, m3/m3

y

y-coordinate, m

C

wellbore storage coefficient, m3/Pa

ye

representative element size in y-direction

yw

y-coordinate of fracture segment midpoint, m

∆Lpf

length of primary hydraulic-fracture segment, m

∆Lsf

length of secondary fracture segment, m

cpft

compressibility of matrix, Pa

-1

compressibility of primary hydraulic-fracture, Pa -1

csft

compressibility of secondary fracture , Pa

h

formation thickness, m

km

matrix permeability, m2

kpf

primary hydraulic-fracture permeability, m2

βD

ksf

secondary fracture permeability, m2

γD

Lpf

primary hydraulic-fracture half-length, m

Lsf

secondary fracture length, m

Msf

number of secondary fracture

Npf

Greek symbols

SC

cmt

-1

RI PT

B

parameter for fracture flow dimensionless transmissibility between fracture

M AN U

segment and its interface

fracture direction, m

θi

inclination of fracture segment to the x-axis

µ

fluid viscosity, Pa·s

number of primary hydraulic-fracture segment

φm

matrix porosity, fraction

Nsf

number of secondary fracture segment

φ pf

primary hydraulic-fracture porosity, fraction

Nf

number of total fracture segment

φsf

secondary fracture porosity, fraction

p

pressure, Pa

pi

initial reservoir pressure, Pa

q

flow rate per length entering the fracture, m2/s 3

production rate, m /s

rw

wellbore radius, m

sc

choking skin factor

TDi,j * TDij

time, s dimensionless

transmissibility

between

fracture

between

fracture

segment i and j

AC C

t

EP

qw

TE D

ε

dimensionless

transmissibility

segment i and j after Star-Delta transformation

Superscripts

−

laplace transform

Subscripts D

Dimensionless

f

Fracture

i

fracture segment

m

matrix

wpf

primary hydraulic-fracture width, m

p

primary

wsf

secondary fracture width, m

s

secondary

x

x-coordinate, m

w

wellbore

xw

x-coordinate of fracture segment midpoint, m

Appendix A – Dimensionless Parameters For the sake of simplicity, we define the dimensionless pressures in the matrix, the primary 18 / 33

hydraulic-fracture and the secondary fracture as pD =

2π km h ACCEPTED 2π kMANUSCRIPT h 2π k h ( pi − p ) , p pfD = m ( pi − p pf ) , psfD = m ( pi − psf qw B µ qw B µ qw B µ

)

(A.1)

For the dimensionless time tD =

km t

(A.2)

φm cmt µ L2pf

the secondary fracture are given by

qD =

2qL pf

qw

, q pfD =

2q pf Lpf

qw

, qsfD =

RI PT

The dimensionless flow rate entering the fracture from the matrix, the primary hydraulic-fracture and

2qsf L pf

qw

(A.3)

η pfD =

km Lpf

, csfD =

k pf (φm cmt )

km (φ pf c pft )

Other dimensionless definitions are

k sf wsf km L pf

, η sfD =

ksf (φm cmt ) km (φsf csft )

∆L ∆L x y y x y ε , yD = ,εD = , xwDi = wi , ywDi = wi , yeD = e , ∆L pfD = pf , ∆LsfD = sf L pf L pf L pf L pf L pf L pf Lpf L pf

TE D

xD =

k pf wpf

M AN U

c pfD =

SC

The dimensionless fracture conductivity and diffusivity are defined by

(A.4)

(A.5)

(A.6)

Appendix B – Matrix Flow Equations

EP

Taking the Laplace transform with respect to time to Eqs. 1 through 6, we have the following form ∂ 2 pD ∂ 2 pD + = spD ∂xD2 ∂yD2

(B.1)

AC C

And, the closed boundary conditions are transformed to  ∂pD  ∂x  D  ∂pD   ∂xD   ∂pD  ∂y D   ∂pD  ∂y D

= 0, xD = 0,0 ≤ y D ≤ yeD = 0, xD = 1, 0 ≤ yD ≤ yeD

(B.2) = 0, 0 ≤ xD ≤ 1, y D = 0 = 0, 0 ≤ xD ≤ 1, y D = yeD

Eqs. B.1 and B.2 describe a 2D flow in the closed rectangular zone. For the secondary fractures, the segments may not align with the x or y-axis. For example, the fracture segment j with the midpoint of (xwDj,

ywDj) has an inclination of θj to the x-axis. A schematic is shown in Fig. B.1. 19 / 33

Thus, by use of the point source function in a rectangular reservoir given by Ozkan and Raghavan

ACCEPTED MANUSCRIPT (1991), we can obtain an inclined line source function in a closed rectangular reservoir by integrating the point source solution with respect to zw over the interval o to h and then with respective to α over the internal –∆LlfDj/2 to ∆LlfDj/2. Then, the appropriate expression of pDM , j in Eq. 7 can be given by

pDM , j =

(

π

2s ∫

∆LlfDj / 2

−∆LlfDj /2

)

(

)

(

)

RI PT

  cosh s y% + cosh s y%   D1 D2      s sinh s yeD   dα  cosh ( ε k y% D1 ) + cosh ( ε k y% D 2 )   ∞ +2∑ cos ( kπ xD ) cos ( kπ x%D )   k =1 ε k sinh ( ε k yeD )  

, y% D 2 = yeD − ( yD + ywDj + α sin θ j ) and

SC

Here, x%D = xwDj + α cos θ j , y% D1 = yeD − yD − ywDj − α sin θ j

(B.3)

M AN U

ε k = s + k 2π 2 hD2 . We evaluate the integral in Eq. B.3 by two parts. The first is the integral on the left hand of plus sign

FA = ∫

∆LlfDj /2

cosh

(

)

s y% D1 + cosh

−∆LlfDj / 2

Eq. B.4 can be derived as

s sinh

(

(

s yeD

)

s y% D 2

)dα

(B.4)

TE D

FA = YM1 , j + YM2 , j

(

)

 sinh s y% D1 1  1 Y =  M,j s sin θ j s sinh s yeD    sinh s y% D 2 1 YM2 , j =  s sin θ j s sinh s yeD 

AC C

EP

where

(B.5)

(

( (

)

∆LlfDj / 2

)

−∆LlfDj / 2

(B.6)

∆LlfDj / 2

)

−∆LlfDj / 2

The second is the integral on the right hand of plus sign FB = ∫

∆LlfDj / 2

−∆LlfDj / 2

cosh ( ε k y% D1 ) + cosh ( ε k y% D 2 )  ∞ 2∑ k =1 cos ( kπ xD ) cos ( kπ x%D )  dα ε k sinh ( ε k yeD )

(B.7)

We integrate the first term in the summation of series XYMk ,1,j = 2cos ( kπ xD ) ∫

∆LlfDj / 2

−∆LlfDj /2

Using the integration by parts, we can derive Eq. B.8 as 20 / 33

cos ( kπ x% D )

cosh ( ε k y% D1 )

ε k sinh ( ε k yeD )

dα

(B.8)

XYMk ,1,j

∆LlfDj /2    cosh ( ε k y% D1 )  1   sin ( kπ x%D ) +   ACCEPTED  kπ cos θ j  MANUSCRIPT  ε k sinh ( ε k yeD )  −∆LlfDj / 2   = 2cos ( kπ xD )  ∆LlfDj / 2  2      sinh ( ε k y% D1 )  1   cos ( kπ x% D ) ε k sin θ j     kπ cosθ j   ε k sinh ( ε k yeD )    −∆LlfDj / 2   

 2 1  / 1 + ( ε k sin θ j )   kπ cosθ j  

  

2

(B.9)

   

XYMk ,2,j = 2cos ( kπ xD ) ∫

∆LlfDj / 2

−∆LlfDj /2

cos ( kπ x% D )

Eq. B.10 can be derived as

cosh ( ε k y% D 2 )

ε k sinh ( ε k yeD )

dα

(B.10)

SC

∆LlfDj /2    cosh ( ε k y% D 2 )  1   % sin k x + π  (  D)  kπ cos θ j   ε k sinh ( ε k yeD )  −∆LlfDj / 2   = 2cos ( kπ xD )  ∆LlfDj / 2  2      sinh ( ε k y% D 2 )  1   cos ( kπ x% D ) ε k sin θ j    kπ cosθ j   ε k sinh ( ε k yeD )     −∆ L / 2 lfDj  

M AN U

XYMk ,2,j

RI PT

In a similar way, we can get

 2 1  / 1 + ( ε k sin θ j )   kπ cosθ j   Hence,

  

2

   

TE D

FB = ∑ k =1 ( XYMk ,1, j + XYMk ,2, j ) ∞

EP

Then, Eq. B.3 can be given by

pDM , j =

(B.11)

π

2s

( FA + FB )

(B.12)

(B.13)

AC C

For pDi , j , we can change (xD,yD) to (xwDi,ywDi) in Eq. B.3.

References

Al-Ahmadi, H.A., and Wattenbarger, R.A., 2011. Triple-porosity models: one further step towards capturing fractured reservoirs heterogeneity. Paper SPE 149054-MS presented at SPE/DGS Saudi Arabia Section Technical Symposium and Exhibition. Al-Khobar, Saudi Arabia. http://dx.doi.org/10.2118/149054-MS. Anderson, D.M., and Mattar, L., 2005. An improved pseudo-time for gas reservoirs with significant transient flow. Paper PETSOC 2005-114 presented at the Canadian International Petroleum Conference. Calgary, Alberta. http://dx.doi.org/10.2118/2005-114. 21 / 33

Brohi, I., Pooladi-Darvish, M., Aguilera, R., 2011. Modeling fractured horizontal wells as dual porosity

ACCEPTED MANUSCRIPT composite reservoirs-application to tight gas, shale gas and tight oil cases. Paper SPE 144057-MS Presented at the SPE Western North American Region Meeting. Anchorage, Alaska, USA. http://dx.doi.org/10.15530/144057-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 Eng. Form. Eval. 14(6), 235–245.

RI PT

http://dx.doi.org/10.2118/125043-PA. Chen, C.C., Raghavan, R., 1997. A multiply-fractured horizontal well in a rectangular drainage region. SPE J. 2(4), 455-465. http://dx.doi.org/10.2118/37072-PA.

SC

Chen, Z., Liao, X., Zhao, X., Lv, S., Zhu, L., 2016. A semianalytical approach for obtaining type curves of multiple-fractured horizontal wells with secondary-fracture networks. SPE J. 21(4), 538–549.

M AN U

http://dx.doi.org/10.2118/178913-PA.

Cipolla, C.L., Fitzpatrick, T., Williams ,M.J., Ganguly, U.K., 2011. Seismic-to-simulation for unconventional reservoir development. Paper SPE 146876-MS presented at the Reservoir Characterization and Simulation Conference and Exhibition. Abu Dhabi, UAE. http://dx.doi.org/10.2118/146876-MS. Granet, S., Fabrie, P., Lemonnier, P., Quintard, M., 2001. A two-phase flow simulation of a fractured using

a

new

fissure

element

method.

J.

Pet.

Sci.

Eng.

32,

35–52.

TE D

reservoir

http://dx.doi.org/10.1016/S0920-4105(01)00146-2. Fisher, M.K., Heinze, J.R., Harris, C.D., Davidson, B.M., Wright, C.A., Dunn, K.P., 2004. Optimizing

EP

horizontal completion techniques in the Barnett Shale using microseismic fracture mapping. Paper SPE 90051-MS presented at the SPE Annual Technical Conference and Exhibition. Houston, TX, USA.

AC C

http://dx.doi.org/10.2118/90051-MS.

He, J., Rui, Z., Ling, K., 2016. A new method to determine Biot’s coefficients of Bakken samples. J. Nat. Gas Sci. Eng. 35, 259-264. http://dx.doi.org/10.1016/j.jngse.2016.08.061. Jia, P., Cheng, L., Huang, S., Liu, H., 2015. Transient behavior of complex fracture networks. J Petrol. Sci. Eng. 132, 1-17. http://dx.doi.org/10.1016/j.petrol.2015.04.041. Jones, J.R., Volz, R., Diasmari, W., 2013. Fracture complexity impacts on pressure transient responses from horizontal wells completed with multiple hydraulic fracture stages. Paper 167120-MS presented at SPE Unconventional Resources Conference. Calgary, Alberta, Canada. http://dx.doi.org/10.2118/167120-PA. 22 / 33

Karimi-Fard, M., Durlofsky, L.J., Aziz, K., 2003. An efficient discrete fracture model applicable for general

ACCEPTED MANUSCRIPT purpose reservoir simulators. Paper SPE 79699-MS presented at the SPE Eastern Regional Meeting, Canton, Ohio, USA. http://dx.doi.org/10.2118/79699-MS. Li, L., Lee, S.H., 2008. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete

fracture

networks

and

homogenized

media.

SPE.

J.

11

(4),

750–758.

http://dx.doi.org/10.2118/103901-PA.

stimulated

reservoir

volume(SRV)?

SPE

Prod.

http://dx.doi.org/10.2118/119890-PA.

RI PT

Mayerhofer, M.J., Lolon, E.P., Warpinski, N.R., Cipolla, C.L., Walser, D., Rightmire, C. M., 2010. What is Oper.

25

(1),

89–98.

SC

Medeiros, F., Ozkan, E., Kazemi, H., 2008. Productivity and drainage area of fractured horizontal wells in tight gas reservoirs. SPE Reserv. Eval. Eng 11 (5), 902–911. http://dx.doi.org/10.2118/108110-PA.

transport

using

discrete

M AN U

Mi, L., Jiang, H., Li, J., Li, T., Tian, Y., 2014. The investigation of fracture aperture effect on shale gas fracture.

J.

Nat.

Gas

Sci.

Eng.

31,

631–635.

http://dx.doi.org/10.1016/j.jngse.2014.09.029.

Mukherjee, H., and Economides, M.J., 1991. A Parametric Comparison of Horizontal and Vertical Well Performance. SPE Form. Eval. 6(2), 209-216. http://dx.doi.org/10.2118/18303-PA.

equation

in

fractured

TE D

Noorishad, J., Mehran, M., 1982. An upstream finite element method for solution of transient transport porous

media.

Water

Resour.

Res.

18

(3),

588–596.

http://dx.doi.org/10.1029/WR018i003p00588.

EP

Ozkan, E., Brown, M., Raghavan, R., Kazemi, H., 2011. Comparison of fractured-horizontal-well performance in tight sand and shale reservoirs. SPE Eng. Form. Eval. 14(2), 235–245.

AC C

http://dx.doi.org/10.2118/121290-PA. Ozkan, E., Raghavan, R., 1991. New Solutions for Well-Test-Analysis Problems: Part I-Analytical Consideration. SPE Form. Eval. 6(3), 359-368. http://dx.doi.org/10.2118/18615-PA. Ou, C., Rui, R., Li, C., Yong, H., 2016. Multi-index and two-level evaluation of shale gas reserve quality. J. Nat. Gas Sci. Eng. 35, 1139-1145. http://dx.doi.org/10.1016/j.jngse.2016.09.056. Stalgorova, E. and Mattar, L., 2012. Practical analytical model to simulate production of horizontal wells with branch fractures. Paper SPE 162515-MS presented at the SPE Canadian Unconventional Resources Conference. Calgary, Alberta, Canada. .http://dx.doi.org/10.2118/162515-PA. 23 / 33

Stalgorova, E. and Mattar, L., 2013. Analytical Model for Unconventional Multifractured Composite

ACCEPTED MANUSCRIPT Systems. SPE Reserv. Eval. Eng. 16(3), 246–256. http://dx.doi.org/10.2118/162516-PA. Stehfest, H., 1970. Numerical inversion of Laplace transforms. ACM Commun. 13(1), 47–49. Sureshjani, M.H. and Clarkson, C.R., 2015a. An Analytical Model for Analyzing and Forecasting Production From Multifractured Horizontal Wells With Complex Branched-Fracture Geometry. SPE Reserv. Eval. Eng. 18(3), 356–374. http://dx.doi.org/10.2118/176025-PA.

RI PT

Sureshjani, M.H. and Clarkson, C.R., 2015b. Transient linear flow analysis of constant-pressure wells with finite conductivity hydraulic fractures in tight/shale reservoirs. J Petrol. Sci. Eng. 133, 455-466. http://dx.doi.org/10.1016/j.petrol.2015.06.036.

SC

Sun, J., Gamboa, E., Schechter, D., Rui, Z., 2016. An intergrated workflow for characterization and simulation of complex fracture networks utilizing microseismic and horizontal core data. J. Nat. Gas Sci.

M AN U

Eng. 34, 1347-1360. http://dx.doi.org/10.1016/j.jngse.2016.08.024.

Tivayanonda, V., Apiwathanasorn, S., Ehlig-Economides, C., Wattenbarger, R.A., 2012. Alternative interpretations of Shale Gas/Oil rate behavior using a triple porosity model. Paper SPE 159703-MS presented at the SPE Annual Technical Conference and Exhibition. San Antonio, TX, USA. http://dx.doi.org/10.2118/ 159703-MS.

TE D

Van Everdingen, A.F., Hust, W., 1949. The application of the Laplace transformation flow problem in reservoirs. J. Pet. Technol. 1(12), 305-324. http://dx.doi.org/10.2118/949305-G. Weng, X., Kresse, O., Chuprakov, D., Cohen, C.E., Prioul, R., Ganguly, U., 2014. Applying complex fracture

EP

model and integrated workflow in unconventional reservoirs. J Petrol. Sci. Eng. 124, 468-483. http://dx.doi.org/10.1016/j.petrol.2015.04.041.

AC C

Wu, Y.S., L, J., Ding, D., Wang, C., Di, Y., 2014. A generalized framework model for the simulation of gas production

in

unconventional

gas

reservoirs.

SPE

J.

19(5),

845-857.

http://dx.doi.org/10.2118/163609-PA. Xu, B., Haghighi, M., Li, X., Cooke, D., 2013. Development of new type curves for production analysis in naturally

fractured

shale

gas/tight

gas

reservoirs.

J

Petrol.

Sci.

Eng.

105,

107-115.

http://dx.doi.org/10.1016/j.petrol.2013.03.011. Zhao, Y.L., Zhang, L.H., Luo, J.X., Zhang, B.N., 2014. Performance of fractured horizontal well with simulated

reservoir

volume

in

unconventional 24 / 33

gas

reservoir.

J.

Hydrol.

512,

447–456.

http://dx.doi.org/10.1016/j.jhydrol.2014.03.026.

ACCEPTED Zhao, X., Rui, Z., Liao, X., Zhang, R. 2015a. MANUSCRIPT The qualitative and quantitative fracture evaluation methodology

in

shale

gas

reservoir.

J.

Nat.

Gas

Sci.

Eng.

27,

486-495.

http://dx.doi.org/10.1016/j.jngse.2015.08.009. Zhao, X., Rui, Z., Liao, X., Zhang, R. 2015b. A simulation method for modified isochronal well testing to determine

shale

gas

well

productivity.

J.

Nat.

Gas

Sci.

Eng.

27,

479-485.

EP

TE D

M AN U

SC

RI PT

http://dx.doi.org/10.1016/j.jngse.2015.08.035.

AC C

Fig. 1. Basic representation of a multifractured system. (a) Real multifractured system. (b) A variation of the trilinear flow model (Stalgorova and Mattar, 2013). (c) Idealized multifractured system in this study.

25 / 33

RI PT

ACCEPTED MANUSCRIPT

Fig. 2. The physical model. (a) Schematic. (b) Dimensions. (b) Discrete fracture segments with application

EP

TE D

M AN U

SC

of Star-Delta transformation.

AC C

Fig. 3. A top view of the global grid configuration of the numerical model

26 / 33

RI PT

ACCEPTED MANUSCRIPT

TE D

M AN U

SC

Fig. 4. Comparison of the resluts of the semi-analytical model and numerical model

(b) Effect of yeD.

AC C

EP

(a) Effect of LsfD.

(c) Effect of Msf.

Fig. 5. Effect of secondary fracture parameters on the assumption of neglecting the un-stimulated region beyond the primary hydraulic-fracture.

27 / 33

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Fig. 6. Transient derivative responses of the system

Fig. 7. Pressure distribution of the region for different flow periods. (a) Bilinear flow (tD=1.E-08). (b) Fluid feed flow (tD=1.E-06). (c) Matrix linear flow (tD=1.E-03). (e) Pseudosteady-state flow (tD=1.E-01). 28 / 33

RI PT

ACCEPTED MANUSCRIPT

Fig. 8. Schematic of flow periods of the system shown in Fig. 2. (a) Bilinear flow. (b) Fluid feed flow. (c)

TE D

M AN U

SC

Matrix linear flow. (e) Pseudosteady-state flow.

AC C

EP

Fig. 9.. Transient derivative responses of the system with consideration of wellbore effect

29 / 33

Fig. 10. Comparison of the results of our model with that of Chen et al. (2016).

SC

RI PT

ACCEPTED MANUSCRIPT

Fig. 11. Square-root-dimensionless-time plot for Eq. (28) with the linear flow solution for the system with

AC C

EP

TE D

M AN U

infinite conductivity

Fig. 12. Square-root-dimensionless-time plot for Eq. (28) with different LsfD and csfD

30 / 33

RI PT

ACCEPTED MANUSCRIPT

TE D

M AN U

SC

Fig. 13. Effect of the secondary fracture number on the transient derivative responses

AC C

EP

Fig. 14. Effect of the secondary fracture length on the transient derivative responses

31 / 33

Fig. 15. Effect of the secondary fracture inclination on the transient derivative responses

SC

RI PT

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

Fig. 16. Effect of the secondary fracture conductivity on the transient derivative responses

AC C

Fig. B.1. Schematic of fracture segment with inclination θj to the x-axis

Table 1 Data used for comparison of the semi-analytical model and numerical model Well, formation and fluid data

Liquid formation volume factor, B (m3/m3) Compressibility of matrix, cmt (MPa-1) Compressibility of primary hydraulic-fracture, cpft (MPa-1) Compressibility of secondary fracture, cpft (MPa-1) Formation thickness, h (m)

1.1

Primary hydraulic-fracture width, wpf (m)

0.01

1.3×10-4

Secondary fracture width, wsf (m)

6.7×10-4

Reservoir size in y-direction, ye (m)

2.3×10-4

Fluid viscosity, µ (mPa·s)

0.87

Matrix porosity, φm (fraction)

0.09

43

32 / 33

0.01 56

Matrix permeability, km (10-3µm2)

2×10-5

Primary hydraulic-fracture porosity, φ pf

(fraction) ACCEPTED MANUSCRIPT Secondary

Primary hydraulic-fracture permeability, kpf (10-3µm2)

fracture

porosity,

80

φsf

0.45

0.21

(fraction) Secondary fracture permeability, ksf (10-3µm2)

20

Primary hydraulic-fracture half-length, Lpf (m)

80

Dimensionless

primary

hydraulic-fracture conductivity, cpfD Dimensionless

secondary

fracture

500

125

conductivity, csfD

4

Initial formation pressure, pi (MPa)

30

Well production rate, qw (m3/s)

0.3

Wellbore radius, rw (m)

0.1

hydraulic-fracture diffusivity, ηpfD Dimensionless

1 2

Dimensionless

Dimensionless

size

secondary

Skin factor for flow chocking, sc

0.5

36

0.5

2

0.1

36

AC C

reservoir

length, LsfD

csfD

33 / 33

fracture

in

y-directon, yeD

LsfD

EP

TE D

3

secondary

diffusivity, ηsfD

Table 2 Cases for the matrix linear flow anlysis Case No.

primary

RI PT

Number of secondary fracture, Msf

Dimensionless

SC

40

M AN U

Secondary fracture length, Lsf (m)

fracture

2.15×104 3.35×104

0.7

0.5 8.19×10-4

ACCEPTED MANUSCRIPT Highlights: New models are developed based on line source solutions in the Laplace domain. Existence of discrete secondary fractures introduces a fluid feed flow. Local solutions for the matrix linear flow are similar to that for a vertical fracture. Secondary fracture inclinations shorten the duration of the matrix linear flow period.

AC C

EP

TE D

M AN U

SC

RI PT

(1) (2) (3) (4)