Laplace-domain multiwell convolution for simulating pressure interference response of multiple fractured horizontal wells by use of modified Stehfest algorithm

Accepted Manuscript Laplace-domain multiwell convolution for simulating pressure interference response of multiple fractured horizontal wells by use of modified Stehfest algorithm Junlei Wang, Ainlin Jia, Yunsheng Wei, Yadong Qi, Yu Dai PII:

S0920-4105(17)30956-7

DOI:

10.1016/j.petrol.2017.11.074

Reference:

PETROL 4490

To appear in:

Journal of Petroleum Science and Engineering

Received Date: 8 September 2017 Revised Date:

31 October 2017

Accepted Date: 30 November 2017

Please cite this article as: Wang, J., Jia, A., Wei, Y., Qi, Y., Dai, Y., Laplace-domain multiwell convolution for simulating pressure interference response of multiple fractured horizontal wells by use of modified Stehfest algorithm, Journal of Petroleum Science and Engineering (2018), doi: 10.1016/ j.petrol.2017.11.074. 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.

ACCEPTED MANUSCRIPT

Laplace-domain multiwell convolution for simulating pressure interference response of multiple fractured horizontal wells by use of modified Stehfest algorithm Junlei Wang1, Ainlin Jia1, Yunsheng Wei1, Yadong Qi1, Yu Dai2

*corresponding author: [email protected], [email protected]

RI PT

1. PetroChina Research Institute of Petroleum Exploration and Development; 2. PetroChina Soutwest Oil & Gasfield Company

M AN U

SC

Abstract: In the multiwell-pad-production scheme, it is a critical task for petroleum engineers to provide an understanding of various interwell communications to quantify well-to-well interference. The focus of this study is put on transient pressure simulation of interwell interference communicated by connecting fracture (i.e. fracture hit) during interference testing of wells with variable operating conditions. On the basis of a framework generalized model of capturing discrete-fracture-network (DFN) flowing, this paper established a rigorous Laplace-transform solution for multiwell pad interfered through fracture-hit connection. In Laplace-transform domain, the variable-rate operating conditions are directly incorporated to integrate with constant-unit-rate solution by using convolution method. Based on the fundamental Stehfest algorithm modified by Chen and Raghavan (1994), which was originally applicable to single-well model, an inverse algorithm was devised to provide the results of transient pressure response of multiwell pad in the real-time domain.

TE D

The approach is validated by comparing with reservoir numerical simulators with regard to interference pressure response of two-well case producing in piecewise-continuous rate sequences. The results are all in an excellent agreement. By conducting sensitivity analysis where the effects of connecting fracture properties (i.e., conductivity, azimuth, number, and geometrical complexity etc.) and rate-sequence changing are considered, detailed discussions on interference diagnostic and identification are presented to capture the physics of interwell interference by utilizing various visualized approaches including pressure profile, flux distribution along fracture face and pressure field. Finally, a field example of multiwell-pad completion, communicated by fracture hit, is used to demonstrate the advantage and improvement of this model over single-isolated-well model and matrix-communication-multiwell model in achieving accurate history matching and providing reasonable fracturing-parameter explanation.

AC C

EP

The new approach combines the flexibility of presenting complex fracture-hit connection and the convenience of convolution of Laplace transformation. The corresponding inverse numerical algorithm has the ability of overcoming difficulties in handling sudden changes in flow rate in well-test-analysis applications. In summary, the study provides a well understanding and insight into the physics of interference mechanisms by using semi-analytical approach. Key words: multiwell pad; pressure interference; connecting fracture communication; Laplace transformation; Stehfest algorithm; discontinuous rate sequence

1. Introduction

With the large scale development of unconventional resources, it is of significant importance to quantify the timing and magnitude of interwell interference in determining optimal well-placement design and assessing the potential of infill drilling or re-fracturing (Lawal et al., 2013; Awada et al., 2016; Esquivel and Blasingame, 2017). In the widely-applied practice of high well density and extensive hydraulically fracturing, the created fractures of an offset well have a tendency to propagate into other existing well through complex fracture connections, commonly referred to as fracture hit, which is explained that hydraulic fractures extend and hit offset well. Officially beginning with the work of Lawal et al. (2013), the fracture hit has been receiving widespread attention in the well-to-well production interference due

ACCEPTED MANUSCRIPT to pressure depletion (Sardinha et al. 2014; Yu et al., 2016, 2017), and fracturing interference due to the changing of in-situ field and associated stress-shadow effect (Morales et al., 2016; Marongiu-Porcu et al., 2016; Wu et al. 2017). According to the classification presented by Sardinha et al. (2012), well interference can be observed during three time periods, including stimulation period (i.e., preflowback phase), productivity testing period (i.e., initial-flowback phase), and long-term established-production period. In this work, the focus is put on facilitating pressure transient analysis during productivitytesting period to detect the impact of pressure-depletion interference on responses of multiwell pad through fracture-hit communication.

SC

RI PT

Interference testing provides qualitative indications of communication between two or more wells in a pad. Such analysis techniques as pressure drawdown and build-up are some of the most frequently used methods in production engineering to estimate the degree of inter-well connectivity between different wells. The impact of interference was evaluated by assessing the change in the productivity index of wells due to offset wells being added, put on production or shut-in (Yaich et al., 2014; Yu et al., 2016). From the viewpoint of filed application, identifying characterization of well-to-well interference is the critical content to achieve an insight into understanding the various mechanisms of interference. Once the DFN flowing of well producing under variable operating condition was captured, we could obtain the best practice of understanding.

AC C

EP

TE D

M AN U

The numerical and semi-analytical models are the most common approaches to detect the interwell-pressure interference throughout the life of multiwell production. Although numerical models have remarkable advantage in obtaining accurate results of interference simulation, the huge time-consuming computation makes this approach limited in pressure-interference analysis (Gupta et al. 2012; Tang et al. 2017). As compared, the advanced semianalytical approaches are applicable compromises in capturing the physics of fracture geometry and pressure-interference mechanisms. As a pioneer in multiwell pressure interference analysis, Economides and Ogbe (1985) provided guidelines for the selection of the appropriate analytical model, and outlined criteria for interference well test. A serious of studies of interference testing analysis was advanced on hydraulically fractured wells and horizontal wells (Meehan et al. 1989; Fokker and Verga, 2002; Cumming et al. 2014; Tung et al. 2016). Noted that the prior studies were built on the assumption of matrix communication. Until to the presentation of fracture hits, Awada et al. (2016) presented a simple analytical model of multi-well pad with fracture hits to identify interwell interference through connected fractures and reservoir matrix. Whereas it was restricted of being applied into more complicated case because of the oversimple assumption of uniform fracture conductivity and planar fracture geometry. Based on the assumption of infinite conductivity, Akinluyi and Hazlett (2016) put the focus on calculating interference time for inter-fracture communication in recycled lean gas injection to displace reservoir fluids between zipper fracs in liquid-rich shales. By use of the time-of-flight-based streamline model, Yu et al. (2016, 2017) built a comprehensive model to describe the details of natural/connecting fracture properties for visualizing the interwell interference process. To accelerate the computation speed and avoid the convergence pitfall caused by time discretization, Jia et al. (2017a) developed a Laplace-domain hybrid model to simulate the transient flow behavior of multiwell communication through interconnected secondary fractures based the pre-existing model (Jia et al. 2015; 2016). The hybrid model was improved by taking non-uniform initial pressure distribution in the fracture, variable operating conditions and two-phase flow in the fracture into account (Jia et al. 2017b). According to the requirement of interference testing, all wells are repeatedly needed to shut-in and put back on production in sequence. Therefore, this case of piecewise-continuous rate sequences is consequently encountered. On the basis of Duhamel’s theorem, variable rate effect is incorporated into an integral with regard to time variable, known as the convolution. The total time interval needs to be partitioned into subinterval to obtain transient response. Fortunately, Laplace transform has the ability of eliminating the requirement for time discretization by convolution, and then is numerically or analytically computed by deconvolution. However, the piecewise-continuous nature of variable

ACCEPTED MANUSCRIPT operating condition would cause another problem for deconvolution in Laplace domain. As a result, the most common numerical inversion algorithms (e.g., Stehfest 1970) used in petroleum engineering cannot be directly implemented. As Roumboutsos and Stewart (1988) stated: the Stehfest inversion algorithm is strictly limited to continuous functions and it will fail drastically if, for example, steps are present in rate schedule.

RI PT

Stehfest numerical inversion algorithm would cause significant oscillations, scattering or completely erroneous results. On the other hand, although other feasible algorithms such as Crump (1976) and Iseger (2006) are able to overcome the issue of discontinuities in pressure or rate, they are limited in the application due to the requirement of the Laplace-transformed function inverted in complex plane, not convenient. Therefore, it is necessary to present a rapid and rigorous approach to simulate transient response during interference testing in multiwell pad with complex connection mechanism.

M AN U

SC

Different from the works of Jia et al. (2017a, 2017b) and Yu et al. (2016, 2017), taking advantage of benefits of Laplace transform and convolution, the focus of this study is to put on incorporating the operating condition of discontinuous rate data into a generalized framework model of multiwell connected through fracture hit. The main contribution of this work is 1) to establish a generalized Laplace-domain model using classical instantaneous source and Green’s functions to describe the multiwell flowing, 2) to incorporate the existence of discontinuities (i.e. changing in flow rate) in the construction of inverse transforms by Stehfest algorithm with the Laplace-domain model. Additionally, synthetic cases were constructed to validate the accuracy by comparing with the results from a numerical simulator, and a field example demonstrated the practical application of the new approach.

2. Multiwell convolution

TE D

Multiwell convolution refers to dual implications, namely the spatial superposition and temporal superposition. Convolution transforms pressure data measured at varying rates into an equivalent unit rate single drawdown. This work considers a multiwell problem involving Nw active wells in a reservoir. To distinguish the data between different wells, the pressure drop and production rate at well ξ are denoted by △p<ξ>(t) and q<ξ>(t). According to the Superposition Principle and Duhamel’s Principle, the pressure drop at any well is the summation of the convolution of flow rates and unit pressure response functions for all wells in the reservoirs, satisfying the following spatial and temporal superposition: ξ

EP

pi − pw (t ) = ∆pw (t ) = ∑ς =w1  ∫ qw (t − τ ) pg 0 ξ

N

t

ς

ξ ,ς

(τ )dτ  

(1)

AC C

where, pi is the initial pressure, pg(t) is the unit response function in the reservoir, and superscript <ξ,ζ> is the response for the effect of well ζ on well and ξ. According to the principle of Laplace transformation, convolution of Eq.(1) is rewritten as the following Laplace-transform form: N ξ ξ ξ ,ζ ∆~ pw ( s ) = ∑ ς =w1 q~w ( s ) ~ pg (s)

(2)

Equation 2 is the well-known convolution equation in Laplace domain. Besides the limit case of constant-rate and –bottomhole-pressure (BHP) condition, the realistic operating condition of production well cannot remain constant in practice. For real-time-domain solution, the variable operating condition can be directly incorporated into the convolution, as shown in Eq.(1). For Laplace-transform solution, it is the critical procedure to construct feasible analytical formulation and reliable inverse algorithm to solve Laplace transformations. Here, the step-constant condition generally refers to the step-rate or pressure schedules, which is common due to the requiring of shut-in in productivity testing and the changing of choke sizes during production management. Furthermore, continuous profile of operating condition can be discretized into

ACCEPTED MANUSCRIPT a set of discontinuous constant-rate or –BHP cases using rectangular discretization. For example, the convoluted response of ξ-th well is written as the following Laplace expression:

~ξ

ϑw (s) =

Nτ

ξ

ϑ ∑ κ =1

ξ

ξ w ,κ

L[ H (t − t

ξ p ,κ −1

1 Nτ ξ ξ )] = ∑ ϑ w,κ exp( −t p ,κ −1 s ) s κ =1

(3)

where H( ) is the Heaviside unit function, tp is the time of response changing, and ϑ =p or q.

~ ξ exp( − st p ,κ −1 ) f ( s ) =

t t −t

ξ p ,κ −1

RI PT

Noted that the Laplace form of Eq.(3) is equivalent to the real-time in a mathematical context, but it would cause error when directly using Stehfest algorithm to Eq.(3) in the computational context. Chen and Raghavan (1994) presented a modified Stehfest algoritm to overcome the above difficulties, and the equivalent relationship was given as

~ f ( s' )

(4)

M AN U

3. Modeling interference response

SC

For t>tp, where s’ is based on t-tp, and detailed derivation is provided in Appendix A. As shown in Fig.1, based on Stehfest algorithm, the results directly using Eq.(3) are not convergent around the drastically changing of time domain, but using Eq.(4) can achieve better matching.

3.1 Physical model

In a multiwell pad, horizontal wells are drilled in parallel and completed with multi-staged transverse fractures. During fracturing procedure, it is of high probability for lateral neighboring wells to communicate each other through fracture connection. Here, we regard the description of physical model in Jia et al. (2017a, 2017b) and Yu et al. (2016, 2017) as the basic assumption in this work.

TE D

As shown in Fig.2, the horizontal well intercepted with multiple finite conductivity fractures is regarded as basic well construction. The neighboring wells are adjoined and communicated through connecting fracture or hydraulic fracture, and wells are designed in aligned or staggered pattern. The emphasis in this study is put on the following two parts: 1) the impact of the properties of connecting fracture on interference response; 2) the impact of operating condition of one well on the other well. In addition, other assumptions are also illustrated next:

EP

The reservoir is isotropic, homogeneous and infinite slab strata with thickness h, porosity φ, permeability km, and compressibility ct;

AC C

Horizontal well is intercepted by Nf hydraulic fractures with a full-length Lf,n, conductivity kfnwfn, angle θn with x-axis, and the located position (xofn, yofn), n=1….Nf; The fractures are assumed to fully penetrate the reservoir thickness; The flux mass from the reservoir to unfractured horizontal segment is not considered, and all production rate is provided by a set of fractures; No pressure loss inside horizontal wellbore is considered, and the radial-convergence effect of flow towards the wellbore within the transverse fractures is ignored; All fractures are grouped as two sets: fractures with mass proppant or completely propped are noted as hydraulic fractures, and fractures with less proppant or partly propped as connecting fractures. Wells are adjoined by several connecting fractures, connected from the fracture tip of one well to the other well. Besides, two coordinate systems are established with regard to reservoir flow of 2D coordinate x-y and fracture flow of 1D coordinate xn. It is noted that W1 is assigned to the pre-existing well, and W2 is assigned to the offset well in Section 4~5 for the convenience of discussion.

ACCEPTED MANUSCRIPT

3.2 Mathematical model The fundamental work of single-well-scale model with DFN is provided in the previous work (Wang et al., 2017). Therefore, the details of single-well derivation are not provided, and only the extension into multiwell case and the consideration of various operating condition are highlighted in this work. For simplicity, the variables are defined in the dimensionless form presented in Appendix A.

3.2.1 Coupled reservoir-fracture flow model

 ∂ 2 p mD ∂ 2 p mD  + 2 ∂y D2  ∂x D

RI PT

Reservoir flow model. In x-y coordinate system as shown in Fig.3, hydraulic or connecting fractures are distributed with arbitrary angles. We define the angle between the fracture panel and the x-axis as θ. Thus, the dimensionless pressure governing equation caused by Nf fractures is expressed as δ m (x D − xofDm − u cos θ m , τ )  ∂p mD N f t D xofDm + L fDm  − q fDm (u , t − τ ) =∑∫ ∫  dud τ 0 x ofDm m =1  ∂t D δ m ( y D − y ofDm − u sin θ m , τ )

(5)

SC

and is subject to initial and boundary condition. Noted that δ( ) is the instantaneous source function.

Nw  Nς  N m ς ξ ς ; n , j  ~ pmD ( n , j ) = ∑ ∑  ∑ q fD ( m, i ) ~ p gD ( m ,i )  ς =1    m =1  i =1

M AN U

After using Laplace transformation, Eq.(5) is written as the product form like Eq.(2). On the basis of the principle of superposition, the reservoir pressure solution close to j-th fracture segment of n-th fracture connected to ξ-th well can be expressed as the discretized form: (6)

ς ς ;n, j In Eq.(5), q fD ( m ,i ) is the dimensionless rate of fracture segment, ~ p gD ( m , i ) is the dimensionless

pressure drop of unit-rate source caused by fracture segment. According to Ozkan and Raghavan (1991) work, the Laplace-domain solution of unit-rate source is given as

TE D

ς ς ς ς  ς ς ( xofDn + xDn, j − xofDm − u D cosθ m ) 2  xDm ,i +0.5 ∆xDm  ς ;n , j ~ pgD ( m,i ) = ∫ ς K0  s ς duD ς ς ς ς xDm ,i −0.5 ∆xDm  + ( yofDn + xDn, j − xofDm − u D sin θ m ) 2   

(7)

EP

Based on the principle of superposition, the multiwell convolution with discretized fracture segments can be expressed as a line system of Nm×Nζ×Nw equations. We vectorize the rate and pressure variables, and Eq.(6) can be rewritten as the following matrix form:

~p + − A + q~ + = 0 mD R wD

(8)

AC C

The element in coefficient matrix and unknown vector refers to the single-well scale, next:

+

AR

 AR1   O  ς = AR  O  N AR w 

1  ~pmD     M   ~ + ~ς  , pmD =  pmD  M     ~p N w   mD

  q~wD+     M  ~ + ~ ς  , qwD =  qwD   M     q~ Nw   wD

        

(9)

ς ς where ~ pmD and q~wD respectively indicate the subvector of dimensionless pressure and influx caused

ς

by ζ-th well. For the detailed expressions of submatrix AR for individual well, please refer to Appendix B. Coefficient matrix is able to reflect the effect of fracture geometries (i.e., number, angle, length and location), independent of conductivity of fracture and operating condition of well.

ACCEPTED MANUSCRIPT Fracture flow model. As shown in Fig.4, in 1D coordinate of planar fracture panel, the pressure equation governing flow within n-fracture connected to the ζ-th well is expressed as ς

∂ 2 p fDn ∂x

2 ς Dn

−

2π ς

C fDn

∫

ς

L fDn

0

ς

ς

q fDn (u , t )δ ( xDn − u ) du +

2π ς

C fDn

ς

N n ,w

∑ [q ς

ς

wfDn ,v

v =1

ς

(t )δ ( xDn − xwfDn ,v )] = 0

(10)

Through using the integral treatment along fracture and Laplace transformation, the final solution for n-th fracture could be generated in the following closed form:

∫

ς

x Dn

0

dζ ∫

ζ

0

2π ς q~fDn (ξ )dξ − ς C fDn

ζ

N n ,w

∑[q~ v =1

ς

ς

wfDn , v

ς

ς

ς

( xDn − xwfDn , v ) H ( xDn − xwfDn , v )]

RI PT

2π ς ς ς ~ p fDn ( xDn ) = ~ p fDn (0) + ς C fDn

(11)

After dividing the n-th fracture into segments of equal length, the fracture-flowing model can be transformed into a discretized equation, as follows:

and

ζ x'  x Dnζ , j ~ ζ ( a ) dx ' ' − x wfDn ,w dx ' ζ q~ ζ ( x ' ' ) dx ' ', T = R dx ' q fDn  ∫0 ∫0 ∫0 ∫0 fDn  ζ =  N n ,w q~ ζ ( x ζ − x ζ ) H ( x ζ − x ζ )  wfDn , v Dn , j wfDn , v Dn wfDn , v ∑  , T = U ζ ζ ζ ζ ζ  v =1 − q~wfDn  ( x − x ) H ( x − x ) ,v wfDn , w wfDn , v wfDn , w wfDn , v   

M AN U

ζ

Tn , j , w

(12a)

SC

2π ς ς ς ς ~ pwfDn,w − ~ p fDn, j = ς (U n, j ,w − Rn, j ,w ) C fDn

(12b)

where xwDfn,w is the location of w-th wellbore on n-th fracture, and pwfDn,w is the pressure of the w-th interconnection on n-th fracture. Eq.(12a) can be further written in the following matrix form (Appendix B): ζ  q~wfDn  ,v   ζ N n ,w  q~ ζ  ζ  wfDn ,v  A ∑ Un , v  M  v =1  ζ   q~   wfDn ,v 

TE D

ζ ~   p wfDn , w  ζ ~  p 2π 2π ζ ζ ζ  wfDn , w  − ~p fDn = ζ AFn , w q~ fDn + ζ  M  C fDn C fDn  ζ  ~ p   wfDn , w 

(13)

AC C

EP

It is noted that 1) pwfD is determined by operating condition and qfD is determined by pressure interference caused by all fracture flowing, and conductivity caused by own fracture flowing; 2) the source/sink terms are contained in the fracture segment connected to wellbore or the interconnected fracture segment connected to another fracture.

Coupled condition. According to the continuity condition, the reservoir pressure pmD for Eq.(8) is equivalent to the fracture pressure pfD for Eq.(13) on the surface of fracture. The Continuity relation between reservoir and n-th fracture model is given as ς ς ς ~ p mDn , j = ~ p fD ( x Dn , j )

(14)

For the case of connection to wellbore, it is characterized by the constrain of flux for individual well, satisfying Nς N m

ς ς q~wD = ∑∑ q fD ( m,i )

(15)

m=1 i =1

Restrain of vector sum of sink/source on the intersection of fracture and fracture enables the flow redirection automatically rather than artificially. For each interconnection in fracture network, the inflow

ACCEPTED MANUSCRIPT rate must equal to outflow rate to ensure the mass balance (assuming steady-state fracture flow and no fluid accumulation or storage within the fracture volume), as

(inflow ) i + ( outflow ) i = q wf , i , i=1…Nυ

(16)

where the positive qwf,i represents sink term, while negative conversely. In addition, the pressures on the interconnection are equivalent each other, which is presented as ~ (17) p =~ p =~ p wfDn , j

wfDh , k

RI PT

wfDm , i

In summary, the analytical reservoir solutions given by Eq.(8) and the closed form of fracture panels given by Eq.(13) have the ability to describe the closed mathematical model.

3.2.3 Solution methodology

+

M

⋅X

+

=Y

SC

Basic solution. By combing Eq.(8) and Eq.(13) on the basis of Eq.(14)~(17), we can obtain the basic solution for the multiwell model. The associated matrix equation is written as +

(18)

M AN U

where vector of X represents unknown pressure and flux distribution, vector of Y represents the known operating condition of well, and the detailed M is presented in Appendix C. Theoretically, Eq.(18) can be solved by combining Newton iteration and Stehfest numerical inverse algorithm to provide the unknowns in the real-time domain.

Superposition solution. Due to the variable-rate operating condition with discontinuous rate sequence, Eq.(18) has the inability to be solved by directly using Stehfest inversion algorithm as demonstrated in Fig.1. To facilitate the application of Stehfest inversion algorithm, the Laplace-transformed expression for unknown variables is given as by using the equivalence of the relationship provided in Eq.(4), which is modified as the matrix equations +

(s) = C 1

+

~ X

+

 Nτ ς ~ + ( s1 ) + t D ∑  ∑ C κ ⋅ X  ς =1  κ = 2 Nw

TE D

~ X

+

 + ( sκ )   

(19)

where the coefficient matrix is a diagonal matrix with regard to discontinuous rate and times, as 1 1 N N ς ς  qwD qwDκ − qwD (κ −1) qwDw1 − qwDw(κ −1)  κ − q wD (κ −1) = ,L, ,L,  1 N ς t D − t pD (κ −1) t D − t pDw(κ −1)   t D − t pD (κ −1)

(20)

EP

Cκ

+

AC C

The final expression is divided into a series of inverse Laplace-domain solution. Laplace variable sk is based on tD-tpD(k) for tD-tpD(k) with tD0=0 and qwD(0)=0 in each expression, which is consistent with Duhamel’s principle in the mathematical context. Here, to solve Eq.(18) provides the basic solution to Eq.(19), and then the elementary solution in bracket of Eq.(19) is given as

~1 1 X ( sκ )   M ~ς ς  X ( sκ )  M  ~ Nw Nw X ( sκ 

−1 ς 1 N   M1,1 ( sκ ) L M1,ς ( sκ ) L M1, N w ( sκ w ) Y~ 1 ( sκ1 )       M M M   M    ~ς ς ς 1 N  =  M ς ,1 ( sκ ) L M ς ,ς ( sκ ) L M ς , N w ( sκ w )  ⋅ Y ( sκ )     M M M   M  ~ N N )   M N ,1 ( sκ1 )L M N ,ς ( sκς )L M N , N ( sκN w )  Y w ( sκ w   w w w w  

       ) 

(21)

Note that the discontinuous-rate operating condition is incorporated into the known vector of Y. The element in Y with regard to operating condition is rewritten as

ACCEPTED MANUSCRIPT ~ς ς ς ς ς ς Y ( sκ ) = ( qwD (κ +1) − qwD (κ ) ) exp( −t pD (κ ) sκ ) ς

and the Laplace variable sκ

(22)

of submatrix in M is correspondingly transformed.

4. Results and discussions 4.1 Model validation

RI PT

In this section, two cases have been provided to verify our approach by comparing with the simulations from numerical simulator. We defined basic well-placement pattern, where the fracture geometries and dimension are illustrated as shown in Fig.2a. The input parameters for two cases are summarized in Table 1, and well spacing is 459.8m. Both flow and shut-in periods have been taken into account, and both increasing and decreasing rate sequences are considered to verify the accuracy of our convolution solution.

M AN U

SC

It is noted that the Stehfest parameter used in this paper are calculated based on Eq.(A4) (the results are list in Table A1), where even integer N is set to be 8. To obtain accurate early-transient flow behavior, the local-grid-refinement (LGR) presented by Bennett et al. (1986) is used to refine the block dimensions in the along-fracture direction and normal-to-fracture direction, and the top schematic is shown in Fig.5. The LGR configuration for global cells is set 408×187×7 to generate refine cells with dimension of 0.1×0.1×0.1m to describe the hydraulic fracture. And all grid dimension are log distributed, ranging from 0.1 to 12.8m. Here, actual fracture width is 0.0127m, and fracture porosity is 0.35. In numerical simulator, equivalent fracture width is set to be smallest grid block, 0.1m. Therefore, the equivalent fracture porosity is calculated as 0.0445, according to the relationship of φe=φfwf/we, and equivalent fracture permeability is calculated to be 120CfD according to kf=CfD(kmxf/wf).

TE D

Two multi-fractured horizontal wells are placed in parallel. Interference through reservoir matrix would occur when the connecting fractures are not connected each other between adjacent wells. Fig.5a shows the comparison of dimensionless pressure between our method and numerical simulations in the case of matrix communication. The spot lines indicate the results of numerical simulation, while the unbroken line represents the response computed in this paper. The dimensionless pressure throughout the whole period was matched in an excellent agreement between two different approaches.

AC C

EP

Based on the above case, the adjacent wells are assumed to be communicated by fracture-hit connection, where one connecting fracture with the same conductivity as hydraulic fracture is set. Fig.5 b provides a comparison of the dimensionless pressure of our model and numerical simulations. The curves were also in excellent agreement during the whole time period. Noted that the method reproduces responses at every rate change; put another way, it handles each discontinuity without loss in accuracy. Therefore, the ability of our algorithms is demonstrated. In addition, the pressure field under two cases at t=200day are provided in Fig.5. It is shown that the characterization of pressure distribution is distinguished remarkably, and the following section would further investigate the mechanisms.

4.2 Simulation of interference response As we know, a change in rate at one well can cause a pressure response change at the second well during productivity testing. In detail, the pressure at the shut-in well is examined to determine whether there is a change in the buildup trend after other well is put on production. In this subsection, we build a basic two-well configuration in aligning pattern: 1) one well with transverse fractures are in extremely close proximity with the other. 2) Fracture tips are connected either with [connecting-fracture communication] or without fracture hits (matrix communication); 3) only one stage is simulated; and 4) fracture number of one stage is 3, dimensionless conductivity of hydraulic fracture is 50, connecting fracture is 10, and dimensionless lengths of hydraulic and connecting fractures is 2.

ACCEPTED MANUSCRIPT In addition, the transient interference response under different operating conditions is analyzed by devising three step-rate cases as shown in Fig.6, where Case 1 was used to detect the fundamental mechanisms and physics of interwell interference, Case 2 was to demonstrate the performance deterioration due to an offset well being put back on production, and Case 3 was to demonstrate the performance improvement due to an offset well being shut in.

4.2.1 Interference mechanisms

RI PT

In interference testing, as described that “all wells are first to build up pressure, and after a certain period, some wells are brought back on production in sequence bottomhole pressure (BHP) ore wellhead pressure of the shut-in well is measured (Lindner and Bello, 2015)”, pre-existing well is set to be on constant-rate production scheme throughout all-life period while the offset well is shut in to observe the pressure response, as shown in Case 1.

M AN U

SC

Connecting fracture & matrix communication. To detect the communicating mechanism, the length of connecting fracture is set to be zero (i.e., well spacing equals to two half-lengths of hydraulic fractures between adjacent wells), and only the inner fracture is connected with each other. Under this assumption, the transient response of fracture communication is equivalent to matrix communication when two wells are on concurrent production schemes with same rate. However, when two wells are controlled by different operating conditions, the response in fracture-communication case would be different from the matrix case in physics.

EP

TE D

Fig.7 illustrates the dimensionless pressure response of two wells with respect to time, where solid line represents pre-existing well and dash line represents offset well. It can be seen that the connection mechanism determines the magnitude of interference time, extremely dependent on conductivity of connecting fracture or matrix permeability. In the matrix communication, the interference would make a difference at tD=0.02. At this time point, the curve of pre-exiting well begins to deviate from single-well production situation, and the dimensionless pressure drop of offset well is about 0.0025 as shown in Fig.7a. If regarding △pD=0.0025 on offset well as the criteria for identifying interference time, interference caused by the onset of pre-exiting well takes effect on advance (tD=0.002) in fracture communication. If regarding the deviation on pre-existing well as the criteria, the interference would appear earlier than tD=10-5. To highlight the influence of fracture-hit communication, Fig.7b illustrates the changing of flux distribution along fracture at different dimensionless times. In fracture-hit communication, the flux along two adjacent fractures is continuously distributed, similar to the case of one asymmetrical fracture; in matrix communication, the flux distribution is piecewise continuous with an obvious gap, where offset well serves as a role similar to isolated natural fractures (Izadi and Yildiz, 2007).

AC C

Fig.8 demonstrates the evolution of pressure-distribution information, where Fig.8a~b indicate the pressure filed for different time steps, and Fig.8c~d indicate the corresponding pressure profile along adjacent connecting fracture. In the case of matrix communication, the reservoir depletion around the pre-existing well is remarkably stronger than the offset well as shown in Fig.8a, verified by a step-like pressure profile in Fig.8c (spot lines indicate the pressure profile of outermost fracture, and solid lines indicate the profile of inner fracture in Fig.8c~d). It needs to emphasize that two-well system with one well producing and one well being shut-in is advantageous over the single MFHW in performance improvement. Since the fractures of offset well plays a role of natural fractures in enhancing the conduction ability of surrounding reservoir. In the case of fracture-hit communication, pressure distribution of two wells behaves like a unified system, similar to a horizontal well with 3 asymmetrical fractures. It is consistent with the regulation of flux profile in Fig.7b. The adjacent connecting fractures make the fluid flow faster from the offset well to pre-existing well. For detailed illustration in Fig.8d, the pressure profile along inner adjacent fractures in fracture-hit communication is continuous, rather than piecewise continuous along outermost fracture in matrix communication.

Effect of properties of connecting fracture. According to the above demonstration, connecting

ACCEPTED MANUSCRIPT fracture plays a significant role in determining the intensity and characteristics of interwell interference response. Therefore, the effect of conductivity (CfDs), number (Nsf) and azimuth (θsf) of connecting fractures are discussed in detail.

SC

RI PT

The effect of connecting-fracture number is detected in Fig.9a, where CfDs=10. Before interference, the pressure in pre-existing well decreases as a result of reservoirs depletion, just as an individual well is on production alone. Once interference takes place, the pressure curves of pre-existing well would deviate from the individual case, and the pressure response of offset well could be also detected. It is noted that 1) the detecting time on offset well in fracture-communication case is obviously smaller than matrix communication (0 fracture), 2) the detecting times of the cases with 1~3 connecting fractures are be an order of magnitude, tD=0.001~0.01. In physics, once adjacent wells are communicated by connecting fracture, the pressure would propagate from pre-existing well to offset well along connecting fracture, nearly independent on number of connecting fracture. On the other hand, more connecting fractures generate a larger drainage area in multiwell system, so a smaller pressure drop is required to maintain same production rate. As demonstrated by pressure filed illustrated in Fig.10a, dimensionless pressure drop of the case with 3 connecting fractures is in a smaller range of 0~0.29, while the range of the case with 0 connecting fracture is 0~0.45.

TE D

M AN U

Fig.9b illustrates the pressure response of two wells under different conductivities of connecting fracture. Once one well is in conjunction with the adjacent well, sudden and significant pressure perturbations are induced through pressure interference. In detail, a higher value of connecting fracture conductivity indicates that pressure propagates from one well to the other well along a highly-effective pathway easily. Correspondingly, the pressure field with different connecting-fracture conductivities at tD=1 are given in Fig.10b, which supports the statement. Taking CfDs=10 for example, two adjacent wells are communicated by one connecting fracture. The reservoirs surrounded by offset well are exploited by the withdrawal of pre-existing well, which is conducted by fracture-hit communication. Put another way, a high degree of connectivity indicates smaller reservoir pressure drop or consumption. As shown in pressure color bar of Fig.10b, the largest dimensionless pressure drop is 0.43 for CfDs=1, 0.40 for CfDs=5, 0.37 for CfDs=10, and 0.30 for CfDs=50.

AC C

EP

Fig.9c displays a group of curves showing varying fracture azimuths which range from 0~150°. The flexibility of considering arbitrary azimuth of fractures is also our advantage over numerical simulation using orthogonal grid. Here, the length of connecting fracture is set to be 2 on the assumption of 3 connecting fractures with CfDs=10. When the angle is exceeding 90°, the two wells would be in the staggered pattern as described by Yu and Sepehrnoori (2014). The azimuth determines the relative position of two wells. When communicated through matrix, the interference time on offset well is dependent on relative position of adjacent wells. Interference time is almost independent on azimuth of connecting fracture when communicated through fracture. We further put the focus on the case of fracture-hit communication. It is emphasized that the matrix communication begins to play a more significant role after interference occurs. As shown in Fig.9c, small variations of pressure response are detected for 0°<θsf<90° due to weak matrix communication (larger well spacing); remarkable variation are detected for θsf>90° due to strong matrix communication (smaller well spacing). The dimensionless pressure drop in Fig.10c also demonstrates that the range for θsf=0° is 0~0.29, while the range for θsf=0° is 0~0.38. In summary, as to be expected, the observations mentioned previously are consistent with the work of Yu et al. (2017) assuming one constant-BHP well and one shut-in well.

4.2.2 Interference measurement In an event when a well’s production was interrupted because of operational events, the degree of pressure response in offset wells was recorded. The simplest way of quantifying interwell interference is to look at the impact on BHP. Four scenarios are presented to investigate the essence of interference by using various operating conditions.

ACCEPTED MANUSCRIPT

SC

RI PT

Using the operating condition as in Case 2, pre-existing well is always on constant-rate production, offset well is put back on production at tD=10. Fig.11a presents the pressure response of two wells by fracture-hit communication. On the beginning when offset well is shut-in, hydraulic fractures connected to offset well play a similar role as isolated fractures. When tD=0.03, the pressure wave arrives at the fractures of offset well, and pwD of pre-existing well noted by blue dashed line is larger than the case of individual well denoted by scatter line. After tD>10 when offset well was put back on production, the pre-existing and offset wells rapidly reacts to the changing of operating condition. After offset well is open to produce, the pre-existing and offset wells were always lower than individual cases in terms of 1/pwD. Put another way, the performance of the pre-existing well deteriorated, and the performance of offset well was lower than the expectations set by the individual production. As the comparison, Fig.11b reflects the response by matrix communication. It is noted that the response in pre-existing well did not rapidly react to the withdrawal of fluid from the offset well which was put back on production at tD=10. Due to the weak interference, the curve of pre-existing well noted by blue dashed line overlaps with individual case until tD=15, smaller than tD=0.03 in fracture-hit communication. It means that there exists delay-time effect of reacting to the changing of operating condition. The comparison between Fig.11a and Fig.11b illustrates the impact of interference intensity on responses of two wells.

TE D

M AN U

It is interesting to point out that pre-existing well would quickly react to the changing of operating condition of offset well if interwell interference has been detected obviously, vice versa. Fig.12 provides two scenarios to confirm the view. When offset well was put back on production at tD=1 as shown in Fig.12a, the pressure propagation caused by pre-existing well did not arrive at the offset well. Hence, interference effect was detected at tD=2.5, indicating strong delay-time effect. Due to weak pressure propagation, the offset well denoted by red dash line is overlapped with individual offset well denoted by red scatter line at the beginning of tD=1. On the contrary, before the offset well was put back on production at tD=100, the pressure propagation had arrived at the offset well at tD=15 when the curve of pre-existing well began to deviate from the individual at a small scale. Afterwards, the response of pre-existing and offset wells would rapidly make a reaction to the offset well startup, as shown in Fig.12b. In addition, the massive withdrawal of fluid from pre-existing well contributed to a lower performance of offset well than the expectation set by the individual.

AC C

EP

Using the operating condition as in Case 3, the pre-existing well is on constant-rate production until tD=1000, and the offset well is on limited-period production and then shut in at tD=10. On the onset of production, two wells were initially on the production with same rate. As shown in Fig.13a, the response curves were overlapped with each other. When tD=0.2, the interference between the communication wells would take effect, and the curves of pre-existing and offset well exhibit a downward trend compared with the individual case. Additionally, there are several findings needed to emphasize: 1) Different from the deviation in Fig.11a in physics, the pressure wave caused by pre-existing well arrived at the fracture network connected to offset well being shut in, and the performance of pre-existing well improved in Fig.11a; the pressure wave caused by pre-existing well encountered wave caused by offset well in the midpoint of well spacing, and then interference took effect, which made the performance of pre-existing and offset wells deteriorate in Fig.13a. 2) After offset well was shut in when tD=10, the performances were improved for a limited period, and then deteriorated due to the strong interference. Noted that the transformed occasions on two wells were same, all at tD=20. 3) After pre-existing well was also shut in when tD=103, the response curves of two well were overlapped again. Generally, one well would rapidly react to the changing of operating condition on the other well once interference occurred. Compared with the fracture-hit communication, there exist two distinguishing features: 1) The delay-time effect was obvious after offset well was shut-in in matrix communication as provided in Fig.13b. It is explained that the pre-existing well reacts to the operating-condition changing of offset well due to the weak intensity of interwell interference. 2) The improvement of fractures of offset well on performance of pre-existing well became weak, which was indicated by a smaller gap denoted in

ACCEPTED MANUSCRIPT Fig.13b. The larger the size is, the stronger the interference would be. On the whole, the interference is detected when the performance of a well improves due to an offset well being shut-in shown in Case 3 or deteriorate when an offset well that was shut-in is put back on production shown in Case 2.

4.3 Application to field case

RI PT

A multiwell pad containing six wells in SouthChina is used to facilitate inter-well interference analysis. All wells are drilled in same layers of the Longmaxi number. The plan view of the pad layout for these wells is shown in Fig.14a. Three wells, Well A, B and C, were used to perform productivity testing, and the rates and BHPs for three wells are presented in Fig.15.

M AN U

SC

As analyzed from Fig.15, Well B was put on production from 3.2 to 7.2 day, while both well A and C were shut in. Pressure drop of well A is decreased by 2.16 MPa during 4-day production of adjacent well B, while well C is decreased by 0.42MPa. Hence, it is of substantial possibility that well A and well B were communicated through connecting fracture. To capture the characterization of the effect of fracture-hit communication on history matching, we build a simplified multiwell model of three multiple fractured horizontal wells communicated by one connecting fracture per one stage, shown in Fig.14b. Considering the limited testing time, well C being shut in is not interference with other wells, so the multiwell model only contains well A and B. Table 2 provides the basic input parameters used in the simulation model, the well spacing between well A and B is about 400m, and the well spacing between well B and C is 300m. For gas well, pseudo-pressure variable is introduced to account for gas property variation with regard to pressure, given as:

p p ( p) =

µ gi Z gi pi

∫

p

p ref

p' dp' µ g ( p' ) Z g ( p' )

(23)

TE D

The approach in this study was then used for matching the interference response in multiwell model. The procedure of history matching is twofold.

EP

Step 1: Initial guess of the properties of hydraulic fractures are obtained based on the single-well model without any interference from other well. These guesses are in a feasible range of reasonable parameters. For history matching using single-well model, it was not possible to obtain an adequate history-matching throughout the whole testing period.

AC C

Step 2: Initial guess serves as the benchmark for the accurate matching. The multiwell models using the initial guesses of fracture properties are established, including both with and without fracture-hit communications. And then the BHPs of both well A and well B are history-matched. The parameters of fracture properties are list in Table 3. The values of fracture properties obtained from the without-frac-hit are larger the with-frac-hit. It is explained that one well could communicate with the other in the frac-hit case, enhancing the performance of multiwell system. If without frac-hit (i.e., matix communication), larger length and conductivity of fractures are required to obtain same matching results, i.e., offsetting the weak communication. Using the parameters from the frac-hit case to match BHP data again, the comparisons between frac-hit and matrix communication are presented in Fig.16. It shows that simulations of pressure drop from matrix communication is larger than the actual in well B matching, and pressure drop from matrix communication is smaller than the actual in well A matching, which are the result of weak communication. These findings are consistent with the insights we obtained by using alternative approaches, such as microseismic mapping and rate transient analysis, etc. The above example justifies the reliability of the new model in this work and corroborates its application potential in the practical field development.

ACCEPTED MANUSCRIPT

5. Conclusion In this work, a semi-analytical approach was developed to simulate pressure interference response in multiwell pad and verified with the results of numerical simulations. The corresponding algorithm was also presented to overcome difficulties in handling discontinuities in flow rate in step-rate testing. Here, some important observations are further emphasized below:

RI PT

(1) The approach has the flexibility to account for various interference mechanisms by utilizing a generalized model of discrete fracture network, and step-rate sequence of testing well by integrating Laplace transformation and modified Stehfest algorithm. (2) Interwell interference could be quantified by observing the pressure responses of neighboring two wells. Sudden and significant pressure perturbations are induced when fractures are connected. When fractures are not connected but are in close proximity, pressure perturbations are so weak that it is difficult to be detected, not existed in essence.

SC

(3) Higher conductivity, larger number and azimuth of connecting fracture could improve the well performance. However, the conductivity plays a more significant role in determining interference time than the number and azimuth of connecting fracture.

M AN U

(4) If the interference was detected, the response of one well would rapidly react to the rate-changing of other well; if not, delay-time effect was remarkable, depending on the time interval between interference time and rate-changing time. (5) When the offset well was shut-in, the performance of pre-existing well would be improved after interference. When the offset well was put back on production, the performance of the pre-existing well would deteriorate and the offset well would be lower than the expectations set by the individual well. The variation trend would be more remarkable in the fracture-hit communication.

EP

Acknowledgement

TE D

(6) Without taking the effect of fracture-hit communication into account, the explained parameters of hydraulic fracture would be larger than the actual, which is the associated result of offsetting the improvement of performance from connecting fracture under the constraint of obtaining same history matching.

AC C

This work was supported by National Science and Technology Major Project of China (Grant No. 2017ZX05063, 2017ZX05037).

Nomenclature Field variables B c h k Lfn pp(p) Np Nυ Nf p pg qwfn,v

Volume factor Compressibility, Pa-1 Formation height, m Permeability, m2 Length of n-th fracture, m Pseudo-pressure function, Pa Number of fracture segment Number of interconnection Number of hydraulic fracture Pressure, Pa Pressure of unit response, Pa Production rate of v-th source in n-th fracture, m3/s

ACCEPTED MANUSCRIPT Flux density along n-th fracture face, m2/s Time, s Initial location of n-th fracture in the x-y coordinate, m Initial location of n-th fracture in the x-y coordinate, m v-th sink/source location within n-th fracture, m Hydraulic fracture width, m Equivalent hydraulic fracture width, m Gas factor Fluid viscosity, Pa·s Fracture porosity Equivalent fracture porosity Fluid density, g/m3

RI PT

qfn t xofn yofn xwfn,v wf we Zg µ φf φe ρ

Dimensionless fracture conductivity Dimensionless time Dimensionless pressure Intersection angle between m-fracture and x-y coordinate

Subscripts f i m D w ref

Fracture property Initial condition Reservoir property Dimensionless Wellbore property Reference variable

M AN U

SC

Dimensionless variables CfD tD pD θm

Reference

TE D

Akinluyi, O. and Hazlett, R. 2016. EOR potential for lean gas reinjection in zipper fracs in liquid-rich basins. SPE-179577-MS presented at the SPE Improved Oil Recovery Conference, Tulsa, Oklahoma, USA, 11-13 April.

EP

Awada, A., Santo, M., Lougheed, D., et al. 2016. Is that interference? A work flow for identifying and analyzing communication through hydraulic fractures in a multiwell pad. SPE Journal. 21(5): 1554-1566. SPE-178509-PA. Bennett, C. O., Reynold, A. C., Raghavan, R., et al., 1986. Performance of finite conductivity, vertical fractured wells in single-layer reservoirs. SPE Formation Evaluation, 8, 399-412.

AC C

Chen, C.C. and Raghavan, R. 1994. An approach to handle discontinuities by the Stehfest algorithm. SPE 28419 presented at the SPE 59th Annual Technical Conference and Exhibition, Louisiana, USA, 25-26 September. Cumming, J.A., Wooff, D.A., Whittle, T., et al. 2014. Multiwell deconvolution. SPE Reservoir Evaluation & Engineering. 17(4): 457-465. Crump, K.S. 1976. Numerical inversion of Laplace transforms using a Fourier series approximation. Journal of ACM. 233: 89-96. Economides, M.J. and Ogbe, D.O. 1985. Single-well and multiwell response interference analysis. SPE 13665 presented at the SPE California Regional Meeting, California, USA, 27-29 March. Esquivel, R. and Blasingame, T.A. 2017. Optimizing the development of the Haynesville shale-lessons-learned from well-to-well hydraulic fracture interference. URTeC: 2670079 presented at the Unconventional Resources Technology Conference, Austin, Texas, USA, 24-26 July. Fokker, P.A. and Verga, F. 2008. A semianalytic model for the productivity testing of multiple wells. SPE Reservoir Evaluation & Engineering. 11(3): 466-477. Gupta, J.K., Zielonka, M.G., Albert, R.A., et al., 2012. Integrated methodology for optimizing development of

ACCEPTED MANUSCRIPT unconventional gas resource. SPE 152224 presented at the SPE Hydraulic Fracturing Technology Conference, Texas, USA, 6-8 Feburary. Iseger, P.D. 2006. Numerical transform inversion using Gaussian Quadrature. Probability in the Engineering and Informational Sciences. 20: 1-44. Izadi, M., and Yildiz, T., 2007. Transient flow in discretely fractured porous media. SPE 108190 presented at the 2007 SPE Rocky Mountain Oil & Gas Technology Symposium, Colorado, USA, 14-18 April.

RI PT

Jia, P., Cheng, L.S., Huang, S.J., et al. 2015. Transient behavior of complex fracture networks. Journal of Petroleum Science and Engineering, 132: 1-17. Jia, P., Cheng, L.S., Huang, S.J., et al. 2016. A semi-analytical model for the flow behavior of naturally fractured formations with multi-scale fracture networks. Journal of Hydrology, 537: 208-220. Jia, P., Cheng, L.S., Clarkson, C.R., et al., 2017a. A Laplace-domain hybrid model for representing flow behavior of multifractured horizontal wells communicating through secondary fractures in unconventional reservoirs. SPE Journal, 1-21

SC

Jia, P., Cheng, L.S., Clarkson, C.R., et al., 2017b. Flow behavior analysis of multi-well communication through secondary fractures in tight oil reservoirs using Laplace domain hybrid model: A field example from the Western Canadian sedimentary basin. URTeC: 2671483 presented at the Unconventional Resources Technology Conference, Texas, USA, 24-26 July.

M AN U

Lawal, H., Jackson, G., Flores, C. 2013. A novel approach to modeling and forecasting frac hits in shale gas wells. SPE 164898 presented at the EAGE Annual Conference & Exhibition incorporating SPE Europec, London, UK, 10-13 June. Lindner, P. and Bello, H. 2015. Eagle Ford well spacing: a methodology to integrate, analyze and visualize multisource data in solving a complex value-fouced problem. URTeC 2174709 presented at Unconventional Resources Technology Conference, San Antonio, Texas, USA, 20-22 July. Marongiu-Porcu, M., Lee, D., Shan, D., et al., 2016. Advanced modeling of interwell-fracturing interference: An Eagle Ford shale-oil study. SPE Journal, 21(5): 1-16.

TE D

Meehan, D.N., Horne, R.N., Ramey, H.J. 1989. Interference testing of finite conductivity hydraulically fractured wells. SPE 19784 presented at the 64th Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 8-11 October.

EP

Morales, A., Zhang, K., Gakhar, K., et al., 2016. Advanced modeling of interwell fracturing interference: An Eagle Ford shale oil study-refracturing. SPE-179177-MS presented at the SPE Hydraulic Fracturing Technology Conference, Texas, USA, 9-11 Feburary. Ozkan, E., Raghavan, R., 1991. New solutions for well-test-analysis problems: Part 1-analytical consideration. SPE Formation Evaluation, 6 (3): 359-368.

AC C

Roumboutsos, A. and Stewart, G. 1988. A direct deconvolution for convolution algorithm for well test analysis. SPE 18157 presented at the 1988 SPE Annual Technical Conference and Exhibition, Houston, Texas, 2-5 October. Sahai, V., Jackson, G., Rai, R., 2012. Optimal well spacing configurations for unconventional gas reservoirs. SPE 155751 presented at the Americas Unconventional Resources Conference, Pennsylvania, USA, 5-7 June. Sardinha, C., Lehmann, J., Pyecroft, J. 2014. Determining interwell connectivity and reservoir complexity through frac pressure hits and production interference analysis. SPE-171628-MS presented at the SPE/CSUR Unconventional Resources Conference, Alberta, Canada, 30 September-2 October. Stehfest, H. 1970. Algorithm 368: numerical inversion of Laplace transforms. Communication of the ACM, 13(1): 47-49. Tang, H.W., Chai, Z., Yan, B.C., et al. 2017. Application of multi-segment well modeling to simulate well interference. URTeC: 2668100 presented at the Unconventional Resource Technology Conference, Austin, Texas, USA, 24-26 July. Tung, Y., Virues, C., Cumming, J., et al., 2016. Multiwell deconvolution for shale gas. SPE-180158-MS presented at the SPE Europec Featured at 78th EAGE Conference and Exhibition, Vienna, Austria, 30 May-2 June.

ACCEPTED MANUSCRIPT Wang, J.L., Jia, A.L., Wei, Y.S., 2017. A generalized framework model for simulating transient response of a well with complex fracture geometry using source and Green’s function. Submitted to Journal of Petroleum Science and Engineering for Peer Review. Wu, K., Wu, B.Y., Yu, W. 2017. Mechanism analysis of well interference in unconventional reservoirs: insights from fracture-geometry simulation between two horizontal wells. SPE Production & Operation. Yaich, E., Foster, R.A., Abou-Sayed, I., 2014. A methodology to quantify the impact of well interference and optimize well spacing in the Marcellus shale. SPE-171578-MS presented at the SPE/CSUR Unconventional Resources Conference-Canda, Alberta, 30 September-2 October.

RI PT

Yu, W., Sepehrnoori, K., 2014. Optimization of well spacing for Bakken tight oil reservoirs. URTeC: 1922108 presented at the Unconventional Resources Technology Conference, Colorado, USA, 25-27 August. Yu, W., Wu, K., Zuo, L.H., et al., 2016. Physical model for inter-well interference in shale reservoirs: relative impacts of fracture hits and matrix permeability. URTeC: 2457663 presented at the Unconventional Resources Technology Conference, San Antonio, Texas, USA, 1-3 August.

SC

Yu, W., Xu, Y.F., Weijermars, R., et al., 2017. Impact of well interference on shale oil production performance: A numerical model for analysis pressure response of fracture hits with complex geometries. SPE-184825-MS presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition, Texas, USA, 24-26 January.

A.1 Fundamental ξ

M AN U

Appendix A: Laplace transformation and Stehfest inversion algorithm

The real-domain solution ϑw (t ) can be transformed into Laplace-domain solution by using Laplace transformation, which is given by ∞

~

ϑwξ ( s) = L[ϑwξ (t )] = ∫ exp(− st )ϑwξ (t )dt , 0 ≤ t ≤ ∞

TE D

(A1)

0

where s is plural in mathematical context satisfying that s=σ+iω. The inversed equation of Eq.(A1) is presented by using analytical expression as following:

~

1 γ + i∞ ~ ξ ϑw ( s ) exp(st )ds 2πi ∫γ − i∞

EP

ϑwξ (t ) = L−1[ϑwξ ( s )] =

(A2) ξ

where Re(s)=γ. It is noted that it is difficult to obtain the analytical expression of ϑw (t ) by directly

AC C

using integral of Eq.(A2). Therefore, the numerical inverse transforms are chosen as an alternative. Stehfest algorithm is very popular in petroleum engineering because of its ease of implementation (Chen and Raghavan, 1994; Fokker and Verga, 2008; Wang et al., 2017). According to the Stehfest algorithm (Stehfest, 1970), the real-domain solution is calculated numerically

ϑwξ (t ) =

ln 2 N ~ξ Vi [ϑw ( s )] i ln 2 ∑ s= t i =1 t

(A3)

Here, N is an even integer, and the Vi’s are constants, which are calculated as

Vi = (−1) N / 2 + i

min( i , N / 2 )

k N / 2 (2k )! ∑ k = ( i +1) / 2 ( N / 2 − k )!k!( k − 1)!(i − k )!( 2 k − i )!

(A4)

Generally, N=8 and the results are list in Table A1.

A.2 Derivation Here, the buildup response is taken as example to derive the equivalence of the relationship presented in

ACCEPTED MANUSCRIPT Eq.(4). For the sake of simplicity, equations are rewritten in the dimensionless form. According to Duhamel’s Principle, Eq.(3) can be rewritten as

ϑsD (t D ) = ϑwD [1 − H (t D − t pD )] = ϑwD (t D ) − ϑwD (t D − t pD )

(A5)

Where ϑsD is the buildup pressure or rate, ϑwD is the dimensionless pressure or rate, tpD is the

RI PT

dimensionless producing time. The corresponding Laplace solution is obtained by using Laplace transformation, which is presented as ~ ~ (A6) ϑsD ( s ) = ϑwD [1 − exp( − st pD )] Application of Eq.(A3) to Eq.(A5) yields

t      ln 2   ln 2  N  1 − exp − i ln 2 pD    ∑ ViϑwD  i t D  i =1 t D      tD 

also, application of Eq.(A3) to Eq.(A6) yields

 ln 2  ln 2  N tD  − ∑ ViϑwD  i t D  i =1  t D  t D − t pD

 ln 2    V ϑ ∑ i wD  i  t − t i =1  D pD  N

M AN U

ϑsD (t D ) =

SC

ϑwD (t D ) =

(A7)

(A8)

It is observed that the right hand sides of Eq.(A7) are equivalent to Eq.(A8). The relation suggests the equivalence as following

~ exp( − st pD )ϑwD ( s ) =

~ tD ϑwD ( s ' ) t D − t pD

(A9)

TE D

For tD>tpD, where s’ is based on tD-tpD. This relation appears to be applicable to all variable-rate conditions by using Duhamel’s principle (Chen and Raghavan, 1994).

Appendix B: Dimensionless-Variable Definitions For the sake of simplicity, the dimensionless definitions are used in this paper.

2πk m h( pi − p f ) 2πkm h( pi − pm ) k mt , p fD = , tD = qref µB φm µctm L2ref qref µB

EP

pmD =

(B1)

AC C

For the fracture flow model, the n-th dimensionless fracture conductivity CfD, dimensionless n-th fracture coordinate in the 1D coordinate xDn, and the dimensionless flow rate qfDn are defined as

C fDn =

k fn w fn km Lref

, xDn =

q L xn ∈[0, LfDn ] , q fDn = fn ref Lref qref

(B2)

Appendix C: Coefficient matrix in coupled model In the reservoir flowing model, the coefficient matrix of ζ-th well in Eq.(9) would be finally transformed into the form of fracture-segment sets, as follows:

ACCEPTED MANUSCRIPT  ARς1,1 L ARς1,m L ARς1, N ς   M M M  ς ς ς =  ARn,1 L ARn, m L ARn, Nς  M M  M  ς  ARN ,1 L ARNς ,m L ARNς , N s s s ς 

ς

AR

ζ ; n ,1 ~ ζ ; n ,1 ~ ζ ; n ,1  ~    p gD ( m ,1) L p gD ( m , i ) L p gD ( m , M m )    M  M M    ζ ;n, j ζ ;n, j ζ ;n, j  , ARnζ , m =  ~ p gD ( m ,1) L ~ p gD ( m , i ) L ~ p gD ( m , M m )     M M   M    ~ ζ ; n , M n  ζ ;n, M ζ ;n, M  p gD ( m ,1) L ~ p gD ( m , i )n L ~ p gD ( m , Mn m )    

RI PT

ζ

(C1)

In the fracture flowing model, matrix AFn, w of w-th source term within n-th fracture of ζ-th well accounts for the effect of influx entering fracture on pressure drop, which is expressed as

α nζ, w;1, M

   M  α nζ, w; j , M n   M   ζ α n, w;M n , M m 

SC

n

M AN U

ζ

AFn , w

 α nζ, w;1,1 L α nζ, w;1,i L   M M  ζ ζ =  α n , w; j ,1 L α n , w; j ,i L  M  M  ζ ζ α L α n , w ; M n ,i L  n , w;M n ,1

ζ

ζ

(C2)

ζ

and the element in sub matrix is satisfied as α n, w; j ,i = RSn, j ;n,i + RTn, w;n,i , where ζ

RSn, j ;n,i

0, i > j  ζ ζ ζ ζ = 0.5(∆xDn ) 2 + ( xDn, j − i∆xDn )∆xDn , i < j  ζ 2 0.125(∆xDn ) , j = 1 or i = j

ζ

0, i > 1 + M wζ, n  ζ ζ ζ ζ = ( xwfDn , w − x Dn ,i ) ∆x Dn , i < 1 + M w,n  ζ ζ ζ ζ ζ ζ w ( xwfDn , w − M w,n ∆x Dn )(0.5 xwfDn , w − 0.25M w,n ∆x Dn ), i = 1 + M n

EP

RTn , w;n ,i

TE D

and

(C3)

(C4)

ζ

In addition, the coefficient matrix AUn,v represents the effect of source/sink term on pressure

AC C

distribution within n-th fracture, which is diagonal matrix next: ζ

ζ

ζ

ζ

AUn,v = [ β n,w,1 Lβ n, w, j Lβ n, w,M n ] ζ

where β n , w, j

ζ ζ ζ ( xDnζ , j − xwfDn  , v ) H ( x Dn , j − x wfDn ,v ) −  = ζ ζ ζ ζ ( xwfDn , w − xwfDn  ) H ( x − x ) ,v wfDn , w wfDn , v  

(C5)

(C6)

Appendix D: Coefficient matrix in Eq.(18) The structure of coefficient matrix M, representing the coupled model, is the sum of matrixes, which is

ACCEPTED MANUSCRIPT +

M

A + B  =D I  0 I 

+

I  0  0 

(D1)

where A represents the flux distribution of fracture segment as +

+

ζ

= AR + [ AF L AF L AF 1

Nw

]

(D2)

B represents the rate of source within fracture as

B

+

= [B L B 1

ζ

Nw

LB

] and B

ζ

ζ

ζ

ζ

= ( β n , w ,1 L β n , w, j L β n , w , M n ) T

D represents the rate relation as +

= [D LD 1

ζ

LD

Nw

] and D

ζ

ζ

ζ

ζ

= ( ∆ x Dn L ∆ x Dn L ∆ x Dn ) T

SC

D

RI PT

A

AC C

EP

TE D

M AN U

In addition, I is the unit matrix.

(D3)

(D4)

ACCEPTED MANUSCRIPT

Table Table 1 Basic parameters used for model validation

value

Basic model parameter

value

Formation thickness, h (m) Formation permeability, km (m2) Formation porosity, φ (%) Fracture width, wHF (m) Fracture half-length, xHF (m) Fracture number

2.9 10-16 10 0.0127 114.95 3

Fracture porosity, φf (%) Fracture distance, Ls(m) Production rate, qsc (m3/s) Initial pressure, pi (Pa) Viscosity of fluid, µ (Pa·s) Total compressibility, ct (Pa-1)

35 1149.5 6.944×10-5 1×108 0.001 4.35×10-10

RI PT

Basic model parameter

Well A

Fracture height (m) Lateral length (m) Fracture conductivity (mD·m) Fracture number Fracture spacing (m) Initial pressure (Pa) Initial gas viscosity (mPa·s) Initial gas compressibility (MPa-1) Initial gas factor Formation temperature (K) Formation porosity (%) Initial water saturation (%)

15 1610 20 45 33.3 38.78 2.65×10-2 1.86×10-2 1.006 389.44 6.75 8.5

Well B

15 1580 20 48 31.25 38.78 2.65×10-2 1.86×10-2 1.006 389.44 6.75 8.5

TE D

M AN U

Basic model parameter

SC

Table 2 Basic parameters used for history model

Well C

15 1600 20 48 31.25 38.78 2.65×10-2 1.86×10-2 1.006 389.44 6.75 8.5

EP

Table 3 Summary of fracture properties derived from history match with and without frac-hit communication

Explained parameter

Initial guess

Without frac-hit With frac-hit

102

133

125

Well B HF length (m)

119

141

135

Connecting SF length (m)

N/A

N/A

44

Well A HF conductivity (mD-m)

15.6

19.8

27.3

Well B HF conductivity (mD-m)

19.8

25.6

35.8

Connecting SF conductivity (mD-m)

N/A

N/A

8.7

AC C

Well A HF length (m)

* HF-hydraulic fracture, SF-secondary fracture Table A1 Stehfest parameters used in this paper when N=8

The constant Vi is calculated based on Eq.(A4) i

Vi

i

Vi

1

-1/3

5

-(14376+2/3)

ACCEPTED MANUSCRIPT 48+1/3

6

18730

3

-906

7

-(11949+2/3)

4

5464+2/3

8

2986+2/3

AC C

EP

TE D

M AN U

SC

RI PT

2

ACCEPTED MANUSCRIPT

6

4

5 3

4 3

2

2 1

Dimensionless pressure, pwD

7

1 0

0 10-2

10-1

100

101

102

103

104

7 6 5

4 4 3

3

2

2 1

1

0

0

-1 106

105

5

8 Numerical simulation Numerical simulation of step-rate Duhamel convolution Directly Stehfest Modified Stehfest

Dimensionless rate, qwD

6

8

5

10-3

(b) 7

9

Dimensionless rate, qwD

Dimensionless pressure, pwD

6

10 Numerical simulation Duhamel convolution Directly Stehfest Modified Stehfest

10-3

10-2

Dimensionless time, tD

10-1

RI PT

(a) 7

100

101

102

-1 103

Dimensionless time, tD

(b)

(a)

1 2 3

M AN U

Stage 2

W1

W2

(c)

W1

+

Time

W1

=0

TE D

Hydraulic fracture Connecting fracture

Rate

EP

Rate

SC

Fig.1 Illustration of the drawback of directly using Stehfest algorithm and the improvement of using modified Stehfest algorithm to the condition of shut-in and flow in the case of a) step-rate sequence and b) hybrid continuous and step-rate sequence (modified from the dual-porosity flow model with ω=10-1, λ=10-8, CD=0, S=0)

wellbore

Time

AC C

Fig.2 Physical model with fracture-hit communication in different well-placement patterns : (a) aligning pattern, (b) staggered pattern, (c) variable-rate operating condition

ACCEPTED MANUSCRIPT xn

y

(xofn+xDn,jcosθn, yfDn+xDn,jsinθn) θn

RI PT

(xofn, yfDn)

(xofm+xDm,icosθm, yfDm+xDm,isinθm)

θm (xofm, yfDm)

x

Fig.3 Location of multiple fractures in a reservoir

SC

(a)

(b)

M AN U

Interface between fracture and reservoir

L fDn

( qwfDn ,1 , pwfDn ,1 )

(qwfDn ,w , pwfDn , w )

q fDi

… xwfDn ,1

…

∆xDn

xDn ,i

x Dn,i− 0.5

xDn

…

i

xwfDn , w

xDn ,i + 0.5

4

EP

3 Well 1 Well 2

2

1

0 10-3

10-2

10-1

100

101

102

Dimensionless time, tD

103

104

(b) 6 rate sequences of adjacent wells

Dimensionless pressure, pwD

5

rate sequences of adjacent wells Well 1 Well 2 tD qwD1 tD qwD1 0-20 2 0-10 1 20-100 0.5 10-50 0.25 100-200 0 50-150 0 200-1000 0.7 150-500 0.35 1000-10000 1 500-5000 0.5 100000 5000 0

AC C

Dimensionless pressure, pwD

(a) 6

TE D

Fig.4 Discretization of fracture into N segments in the coupled model

5

4

3

2

Well 1 tD qwD1 0-20 2 20-100 0.5 100-200 0 200-1000 0.7 1000-10000 1 100000

Well 2 tD 0-10 10-50 50-150 150-500 500-5000 5000

qwD1 1 0.25 0 0.35 0.5 0

Well 1 Well 2

1

0 10-2

10-1

100

101

102

Dimensionless time , tD

103

104

ACCEPTED MANUSCRIPT Fig.5 Model verification against numerical simulations with piecewise-continuous rate sequence in the case of a) matrix communication and b) fracture-hit communication.

Pre-existing well

Pre-existing well

Rate

Rate

Rate

Pre-existing well

Time

Time

RI PT

Time

Offset well

Offset well

Rate

Time

Time

Case 1

Case 2

Time

SC

Rate

Rate

offset well

Case 3

M AN U

Fig.6 Design of step-rate sequence used in Section 4.2 101

By fracture

10

10-1

10

rve Cu

gin be

d s to

at evi

rom ef

gle sin

tD=10-3

Matrix communication Fracture-hit communication Single well

-2

criteria for identify interference time 10-3

10-4 10-5

10-4

10-3

10-2

10-1

100

tD=10-1

tD=101

tD=103

tD=105

0.3

Dimensionless flux distribution, qwD

0

TE D

Dimensionless pressure, pwD

0.4

se ca ell -w

101

102

103

104

0.2 0.1 0.0

1.2

By matrix

0.8 0.4 0.0

-0.4 -0.8

105

-1.2 0.0

0.5

1.0

1.5

Dimensionless time, tD

(a) Pressure perturbation

EP

2.0

2.5

3.0

3.5

4.0

Dimensionless distance, xD

(b) Flux distribution

Fig.7 Interference simulation under various connected mechanisms 6

-5.00E-04

1.45E-02

6

5

AC C

5

0.158

6

0.258

5

0.358

4

2.08E+00

7.45E-02

2

0.458

3

4

1.60E-01

-1

0

1

2

i) tD=10-1

3

4

5

6

2

2

2.68E+00

0.0595

2

0.658

0

1

2

ii) tD=101

3

4

5

0.4800 0.5400

2

1

2.515

0.6600

2.98E+00

0.0895

0

0.7200

0

2.715

1

-1

3.08E+00

0.958

-1

0.110

-1

0.8400

3.18E+00

-2

3.23E+00

-2

-1

0

1

2

3

iii) tD=103

(a) Pressure filed in matrix communication

4

5

6

0.120 0.124

-2 -2

-1

0

1

2

3

-1

i) tD=10

4

5

6

2.815

0

0.7800

0.0995

1.04

2.615

0.6000

1 0.0795

0.858

6

2

2.88E+00

0

-2

-1

2.415

3

0.0695

2.78E+00

1 0.758

0

-2

2.215

4

0.4200

0.0495

2.58E+00

-1

-2

4

3

1.50E-01

-2

0.0295

0.3600

0.0395

1.35E-01

-1

2.115

5 0.3000

0.558

1

1.20E-01

6

0.2400

0.0195

3

1.05E-01

0

0.1800

5

2.48E+00

3

8.95E-02

1

6

0.00950

2.315 2.38E+00

4.45E-02

5.95E-02

-5.00E-04

5 2.28E+00

4

3

6

2.18E+00

2.95E-02

4

2.915

-1 3.015

0.9000

-2

0.9380

-2

-1

0

1

2

3

1

ii) tD=10

4

5

6

-2

3.115 3.125

-2

-1

0

1

2

3

iii) tD=103

(b) Pressure filed in fracture communication

4

5

6

tD=10-1

1.2 1.0 0.8 0.6 0.4 0.2

tD=10

1

3.4 3.2 3.0 2.8 2.6 2.4 2.2

tD=103

5.6 5.4 5.2 5.0 4.8 4.6

tD=105

-2

-1

0

1

2

3

4

5

0.15

tD=10-1

0.10 0.05 0.00 1.0

tD=101

0.8 0.6 0.4 0.2 3.2 3.0 2.8 2.6 2.4

tD=103

5.4 5.2 5.0 4.8 4.6

tD=105

6

-2

-1

0

RI PT

0.20 0.15 0.10 0.05 0.00

Dimensionless pressure distribution, pD

Dimensionless pressure distribution, pD

ACCEPTED MANUSCRIPT

1

2

3

4

5

6

Dimensionless distance, xD

Dimensionless distance, xD

(c) Pressure profile in matrix communication

(d) Pressure profile in fracture communication

10

-1

10

-2

10

-3

10

-4

10-5

Dimensionless pressure, pwD

0

The effect of number of connecting fracture is relatively smaller in the fracture-hit communication

10-4

10-3

10-2

10-1

100

M AN U

10

(b) 101

connected by 0 frac connected by 1 frac connected by 2 frac connected by 3 frac

100

101

102

103

104

10-1

By matrix By frac (CfDs=1) By frac (CfDs=5)

10-2

By frac (CfDs=10) By frac (CfDs=50)

10-3

the larger connecting fracture conductivity is, the earlier interference occurs

-4

105

10

10-4

TE D

Dimensionless pressure, pwD

(a) 101

SC

Fig.8 Evolution of pressure information under different communication mechanisms

Dimensionless time, tD

10-3

10-2

10-1

100

101

102

103

104

105

Dimensionless time, tD

(c) 101

pre-existing well

10

By fracture: θsf=0° θsf=30°

AC C

EP

Dimensionless pressure, pwD

0

10

-1

10

-2

By matrix: θsf=0° θsf=150°

θsf=45° θsf=120° θsf=150°

101 100

offset well

10-1 10-2

Interference time on shut-in well is independent on azimuth of fracture in the fracture-hit communication

10-3 10-4 10-4

Interference time on shut-in well is dependent on azimuth of fracture in the matrix communication

10-3

10-2

10-1

100

101

102

103

104

105

Dimensionless time, tD

Fig.9 Effect of properties of a) number, b) conductivity and c) azimuth of connecting fracture on interference response

ACCEPTED MANUSCRIPT 7

0.0E+00

7

0.0E+00

7

6

3.5E-02

6

5

7.0E-02

0.0E+00

9.0E-02

6.0E-02

5

2.5E-02

5.0E-02

5

8.0E-02

1.2E-01

7.5E-02

1.1E-01

4

4

4

1.0E-01

1.5E-01 1.4E-01

1.8E-01

3

0.0E+00

6 4.0E-02

6.0E-02

5

7

2.0E-02

3.0E-02

6

3

4 1.0E-01

1.2E-01

3

3 1.4E-01

2.1E-01

2

2

2.4E-01

1.6E-01

2 2.1E-01

2.7E-01

1

1.3E-01

1.8E-01

1

2

1.5E-01

1.8E-01

1

2.0E-01

1

1.8E-01

0

2.0E-01

2.5E-01

3.0E-01

2.2E-01

0

0

3.3E-01

0 2.8E-01

2.4E-01

3.6E-01

-1

-1

3.2E-01

3.9E-01

-2

-2

3.5E-01

-3

-3 0

1

2

3

4

5

6

7

-2

3.0E-01

3.7E-01

4.5E-01

-1

2.3E-01

-1

2.8E-01

4.2E-01

-2

2.6E-01

-1

8

-3 -2

-1

0

1

2

3

4

5

6

7

8

2.8E-01 2.9E-01

-3 -2

-1

0

1

ii) Nsf=1

i) Nsf=0

2.5E-01

-2

3.2E-01

2

3

4

5

6

7

8

-2

-1

0

1

2

iii) Nsf=2

3

4

5

6

7

8

iv) Nsf=3

(a) Effect of number of connecting fracture on pressure filed -1.0E-02 2.0E-02

6

7

0.0E+00

5

8.0E-02

0.0E+00

7.5E-02

5

4.0E-02

1.3E-01

5

8.0E-02

4

1.2E-01

1.5E-01 1.7E-01

3

1.8E-01

2.0E-01

2

2.0E-01

2 2.3E-01

1

2.6E-01

2.5E-01

4.0E-02 6.0E-02 8.0E-02

4

1.0E-01

1.4E-01

3

1.6E-01 1.8E-01

2

1.2E-01

3

1.4E-01

2

1.6E-01

2.0E-01

2.3E-01

1

2.0E-02

5

1.0E-01

1.4E-01

3

0.0E+00

6

6.0E-02

1.0E-01

4

7

2.0E-02

6 5.0E-02

1.1E-01

4

7

2.5E-02

6

5.0E-02

RI PT

7

2.2E-01

1

1.8E-01

1

2.4E-01

2.9E-01

0

2.8E-01

0

3.0E-01

3.2E-01 3.5E-01

-1

3.3E-01

-1

4.1E-01 4.3E-01

-3 -2

-1

0

1

2

3

4

5

6

7

0

-2

3.8E-01

3.0E-01

-1

8

3.4E-01

-1

0

1

i) CfDs=1

2

3

4

5

6

7

2.4E-01 2.6E-01 2.8E-01

-2

3.6E-01 3.7E-01

-3 -2

2.2E-01

-1

3.2E-01

-2

4.0E-01

-3

2.0E-01

0

2.8E-01

3.5E-01

3.8E-01

-2

2.6E-01

8

3.0E-01 3.0E-01

-3

-2

-1

0

1

ii) CfDs=5

2

3

4

5

6

7

8

-2

-1

0

1

2

3

4

5

6

7

8

vi) CfDs=50

iii) CfDs=10

7

0.0E+00

6

2.5E-02

6

0.0E+00

SC

(b) Effect of conductivity of connecting fracture on pressure filed 6

0.0E+00

4.0E-02

5

5

4.0E-02

4

8.0E-02

4

8.0E-02

4 3

3

3

1.2E-01

1.2E-01

1.3E-01

1.5E-01

2

2

1.8E-01

1.8E-01

1

1.6E-01

2.2E-01

0

1

2.0E-01

2.5E-01

2.6E-01

-1

-2

2.8E-01

-2 -2

-1

0

1

2

3

4

5

6

7

8

-1

0

1

2

3

ii) θsf=30°

i) θsf=0

4

5

6

-1

0

1

2

3

4

5

2.1E-01 2.4E-01

1

0

3.0E-01 3.3E-01

2.6E-01

-1

-2

1.8E-01

2.4E-01

-1 2.8E-01

3.6E-01

2.9E-01

-2

-2

1.5E-01

2.7E-01

0

2.8E-01 2.9E-01

2.9E-01

-3

1.2E-01

2

2.2E-01

2.4E-01

-1

9.2E-02

3

1.8E-01

2.0E-01

2.0E-01

0 2.3E-01

M AN U

1.6E-01

1

6.2E-02

4

1.4E-01

1.4E-01

2

3.2E-02

1.0E-01

1.0E-01

1.0E-01

2.0E-03

5

6.0E-02

6.0E-02 7.5E-02

6

2.0E-02

2.0E-02

5 5.0E-02

3.8E-01

-2 6

-2

-1

0

1

2

3

4

iv) θsf=150°

iii) θsf=45°

(c) Effect of azimuth of connecting fracture on pressure filed

Fig.10 Dimensionless pressure field of different fracture properties at tD=1

TE D l el

EP

100

w nt ce ja d a ct of ffe e es ur tive t i c s fra po e a Th ve a h Pressure propagation arrives at the offset well at tD=0.03

10-1

-2

10

10-4

AC C

Dimensionless pressure, pwD

pre-existing well offset well pre-existing well (alone) offset well (alone)

10-3

10-2

10-1

100

Dimensionless pressure, pwD

1 (b) 10

(a) 101

pre-existing well offset well pre-existing well (alone) offset well (alone)

100

Pressure propogation arrives at the offset well at tD=15

10-1 large size of gap

Offset well is put back on production at tD=10

Offset well is put back on production at tD=100

10 101

Dimensionless time, tD

102

103

104

-2

10-4

10-3

10-2

10-1

100

101

102

Dimensionless time, tD

Fig.11 Interference response of two adjacent wells by a) fracture-hit communication and b) matrix communication when offset well that was shut-in is put back on production

103

104

ACCEPTED MANUSCRIPT 1 (b) 10

Dimensionless pressure, pwD

Dimensionless pressure, pwD

obvious interference is detected at tD=2.5

100

pre-existing well offset well pre-existing well (alone) offset well (alone)

10-1

-2

10

pre-existing well offset well pre-existing well (alone) offset well (alone)

100

Pressure propogation arrives at the offset well at tD=15

10-1

RI PT

1

(a) 10

large size of gap

Offset well is put back on production at tD=1

small size of gap

Offset well is put back on production at tD=100

-2

10

10-3 10-4

10-3

10-2

10-1

100

101

102

103

10-4

104

10-3

10-2

10-1

100

101

102

103

104

Dimensionless time, tD

Dimensionless time, tD

1.5

Obvious interference is detected at tD=0.2

TE D

0.5 0.0 10-1

100

101

102

103

2.5 2.0

's

th e fra im p ct ur rov es e is men we t ak of o

M AN U

to

2.0

Dimensionless pressure, pwD

2.5

pre-existing well offset well pre-existing well (alone) offset well (alone)

3.0

Th e f ha r ac ve tur e a po s of sit o ive ffse ef t w fe ct ell b e

3.0

delay-time effect

1.5 1.0 0.5 0.0

10-1

104

100

101

102

103

Dimensionless time, tD

Dimensionless time, tD

(a)

C

EP

Fig.13 Interference response of two adjacent wells by a) fracture-hit communication and b) matrix communication when offset well that was on production is shut in

(b)

B

C

A

B

A

AC C

Dimensionless pressure, pwD

gi ns

pre-existing well offset well pre-existing well (alone) offset well (alone)

ffs et

(b) 3.5

(a) 3.5

1.0

SC

Fig.12 Effect of startup time of offset well at a) tpD=1 and b) tpD=100 on interference response in the matrix communication

400m

…

…

…

…

…

…

300m

Fig.14 Information of multiwell pad a) plan view of layout b) conceptual model

104

ACCEPTED MANUSCRIPT 22 20 18 16 14 12 10 8 6 4 2 0

37 36 35 Well A Well B Well C

34 33 32 31 30 3.0

3.5

4.0

4.5

5.0

5.5

6.0

Testing time (day)

RI PT

Bottomhole pressure (MPa)

38

flowing rate (104m3/d)

.

6.5

7.0

38

18 ct o

f fra

36

c-hit

c om m un

35

icati

on

14 12

34

10

history BHP data matching with frac-hit matching without frac-hit

33

8 6

32

4

31

37.5

16

Bottomhole pressure (MPa)

effe

flow rate (104m3/d)

itive

M AN U

pos

37

37.0

weak interference through matrix communication

36.5

36.0

history BHP data matching with frac-hit matching without frac-hit

35.5

2

30 10

20

30

40

50

60

70

Testing time (hour)

80

90

0

35.0 0

10

20

30

40

50

60

Testing time (hour)

Fig.16 History pressure matching of a) well B and b) well A using multiwell model

EP

0

AC C

bottomhole pressure (MPa)

(b) 38.0

20

TE D

(a)

SC

Fig.15 History pressure and rate data of well A, B and C

70

80

90

ACCEPTED MANUSCRIPT

1. A semi-analytical approach using classical Green’s function and instantaneous source function is established to capture the characterization of complex fracture-hit connection in multiwell pad. 2. An integrated solution of Laplace transformation and modified Stehfest algorithm is extended into calculating the real-time-domain results in the application of multiple multifractured horizontal wells.

RI PT

3. Detailed discussions on interference diagnostic and identification are presented to detect the physics of interwell interference by utilizing various visualized approaches.

AC C

EP

TE D

M AN U

SC

4. A field example is provided to confirm the advantage of our model in achieving accurate history matching over single-well model and multiwell model without fracture hit.