Multiphase flowback rate-transient analysis of shale gas reservoirs

Multiphase flowback rate-transient analysis of shale gas reservoirs

International Journal of Coal Geology 217 (2020) 103315 Contents lists available at ScienceDirect International Journal of Coal Geology journal home...

2MB Sizes 5 Downloads 112 Views

International Journal of Coal Geology 217 (2020) 103315

Contents lists available at ScienceDirect

International Journal of Coal Geology journal homepage: www.elsevier.com/locate/coal

Multiphase flowback rate-transient analysis of shale gas reservoirs Fengyuan Zhang, Hamid Emami-Meybodi



T

Department of Energy and Mineral Engineering and EMS Energy Institute, The Pennsylvania State University, University Park, PA 16802, USA

A R T I C LE I N FO

A B S T R A C T

Keywords: Flowback Multiphase flow Rate transient analysis Shale reservoirs Multi-fractured horizontal well

Multi-fractured horizontal well (MFHW) is commonly used in the development of unconventional reservoirs. Flowback data after hydraulic fracturing is of critical importance in the characterization of hydraulic fracture, stimulation evaluation, and reservoir simulation. This study presents a two-phase diagnostic plot and a semianalytical flowback model for the early-time flowback period when fluid influx from matrix remains insignificant and production is mainly from the fracture network. The developed model considers the changes of water saturation in hydraulic fracture (HF) under fracture depletion mechanism and pressure-dependent permeability and porosity, and is able to predict HF attributes, such as fracture half-length, initial fracture permeability, and initial pore volume of fracture network. A new method of calculating average pressure in the fracture is developed based on the solution of water-phase diffusivity equation and is compared with the traditional material balance approach based on the tank model. Three iterative workflows are proposed using the successive substitution method to calculate fracture attributes under different production conditions. Validation of the developed flowback model, average pressure calculation method, and the analysis workflow, as well as their application in an MFHW drilled in Horn River Shale show the applicability of this study to interpret flowback data quickly and estimate HF properties accurately.

1. Introduction Multi-fractured horizontal wells (MFHWs) have been widely used to economically produce from unconventional resources, such as tight and shale gas. Rate transient analysis (RTA) is a commonly used technique to conduct quantitative analysis on the well performance of MFHWs. Many mathematical tools have been investigated to characterize reservoir and HF properties using RTA (Clarkson, 2013). However, traditional RTA methods usually focus on long-term production data that require the well producing several months or even several years in ultratight reservoirs. In recent years, more and more studies focus on flowback data right after stimulation to obtain an early estimation of fracture characteristics and stimulation effects, which guides efficient exploitation in the future (Clarkson and Williams-kovacs, 2013). Early studies of flowback analysis had focused on the qualitative insights and flow regime analysis. Crafton (1998) used the reciprocal productivity index method to analyze post-stimulation flowback data after hydraulic fracturing to evaluate productive capacity, and later studied the effect of flowback rate on well performance and recovery (Crafton, 2008). Ilk et al. (2010) gathered detailed time-rate-pressure flowback data from five wells in the same shale gas reservoir and proposed a procedure to perform qualitative analysis using numerous



plots. Although this approach is not a model-based workflow that leads to the characterization of reservoir properties, it provides an important empirical diagnostic plot: gas water ratio (GWR) vs. cumulative gas production (Gp) plot in log-log scale, which is widely used by many authors. Clarkson and Williams-Kovacs (2013) applied the half-slope straight line of this plot to identify multiphase depletion flow in the fracture to conduct flowback analysis. Abbasi et al. (2014) pointed out that GWR vs. Gp plot for shale gas wells shows two distinct regions of first decrease and then increase with Gp. Adefidipe et al. (2014) introduced early gas production and late gas production concepts to divide flowback period based on the “V” shape of GWR vs. Gp plot. In recent works, many authors investigated the mathematical development of flowback models to conduct quantitative analysis. During the early flowback period, the fluids (water or water and gas) inside the fracture network flow along the fractures to the wellbore. Meanwhile, formation gas may enter the fracture network creating a coupled fracture-matrix flow. Early attempts of flowback modeling focused on fluid flow in the fracture. Abbasi et al. (2014) developed a single-phase water model to analyze flowback data of MFHWs in shale gas or oil reservoir. But this model is only valid for the single-phase water dominated flow regime, which usually lasts a very short time. Alkouh et al. (2014) considered gas matrix influx into the fracture and proposed a governing

Corresponding author. E-mail address: [email protected] (H. Emami-Meybodi).

https://doi.org/10.1016/j.coal.2019.103315 Received 30 April 2019; Received in revised form 16 October 2019; Accepted 17 October 2019 Available online 05 November 2019 0166-5162/ © 2019 Elsevier B.V. All rights reserved.

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

equation that combines water and gas flow in the fracture. To obtain an analytical solution, simplification was applied to the definition of total compressibility and material balance equation. Uzun et al. (2016) applied a similar simplification on total compressibility in their work. However, this method has a limitation that it requires the diffusivity equations to be linear. Xu et al. (2016) developed a flowing material balance equation for flowback analysis based on gas-phase model, which might cause errors when is applied to early-time flowback data dominated by water-phase flow. Yang (2017) proposed two semi-analytical models, one for early gas production ignoring fluid influx into the fracture, and the other for late gas production assuming the same pressure drop in fracture and matrix. More recently, Jia et al. (2018) developed a one-dimensional (1D) two-phase flowback model based on constant rate solution using water flowing material balance equation in an iterative procedure. Their model is applicable to cases with small gradients of saturation and pressure, such as with small fracture halflength and high fracture permeability. However, the model may not be adequate for systems experiencing significant saturation and pressure gradients. To consider fluid influx in flowback modeling, advanced analysis methods continue to evolve by coupling fracture flow with single-phase or two-phase influx from the matrix. Ezulike and Dehghanpour (2014) developed a two-phase bilinear flow model to analyze flowback data with gas influx. They introduced Dynamic Relative Permeability (DRP) concept and applied Laplace transform to the coupled partial differential equations. Although this model is robust, the main drawback is that the proposed DRP function was obtained by curve fitting, which requires known matching parameters. Clarkson and Williams-kovacs, (2013) applied multiphase RTA model in CBM to determine fracture properties in multiphase fracture depletion flow during the boundarydominated flow period. In an extended study, Williams-Kovacs and Williams-Kovacs and Clarkson, (2016) divided the flowback period into three different regimes and modeled after-breakthrough flow period using a modified material balance equation. Later, Clarkson and Qanbari (2016a) adopted the Dynamic Drainage Area (DDA) concept from well testing to conduct flowback analysis based on boundary dominated flow (BDF) solution with material balance equation within the distance of investigation in the matrix. The DDA method for flowback modeling was later improved and applied to MFHWs in CBM, tight oil and shale reservoirs for history matching, production forecasting and RTA (Clarkson and Qanbari, 2016b; Williams-Kovacs and Clarkson, (2016); Qanbari and Clarkson, 2016). More recent studies pointed out the importance of incorporating complex fracture geometries to estimate fracture attributes using two-phase production data (Altwaijri et al., 2018, Chen et al., 2018, Jia et al., 2017, Yang et al., 2017). However, such complex models demand expensive computational time due to the existence of local grid refinement or unstructured grids making their field application impractical. Apart from the flow-regime analysis, Jones and Blasingame (2019) applied a decline curve analysis to empirically estimate initial fracture volume and forecast total liquid productivity index. The present study focuses on analyzing the flowback data during primary fracture depletion. We developed a 1D model to estimate HF attributes and forecast the flowback behavior during early gas production from the fracture network, and proposed a two-phase diagnostic plot to identify flow regimes. The developed model provides a quick estimation of fracture properties and can be used for further development of more complex scenarios. A new average pressure calculation approach is proposed based on the solution of diffusivity equation, which releases the need for gas-phase flowback data as input. The proposed flowback model, average pressure method, and the analysis workflow are validated against several numerical simulations and applied to an MFHW drilled in Horn River Shale in British Columbia, Canada.

1D Model (Top view) fracture

Hydraulic fracture h Water Gas

Wellbore

Shale matrix

Fig. 1. A schematic of the conceptual model for flowback depicting flow through hydraulic fracture only under fracture depletion mechanism.

2. Theory and methods We first describe the conceptual model and the assumptions made in this study. RTA flow regime analysis on flowback data is introduced as a basis for further mathematical development and model application. Then, a 1D two-phase flowback model for boundary dominated flow (BDF) is developed. We also proposed a new analytical method to calculate average pressure in the fracture which does not rely on gas production data as the input, and compared the new method with the material balance approach. At the end of this section, the method of successive substitution is applied to iteratively obtain key fracture properties, such as fracture half-length, initial fracture permeability, and initial pore volume of the fracture. 2.1. Conceptual model As shown in Fig. 1, we considered an element of symmetry from a planar fracture network for the modeling of 1D flow during the early gas production. Here, HF is treated as a porous medium initially saturated with fracturing liquid (water) and hydrocarbon (gas). Even though fracturing fluid will leak off during stimulation job, most of the water in the matrix cannot be recovered and will be trapped due to imbibition and capillary pressure (Alkouh et al., 2014). Gas production only sources from initial free gas stored in fracture and gas influx from the matrix is negligible. Based on these explanations, the main assumptions include:

• Water only sources from the hydraulic fracture and matrix con• • • • •

tribution to the water production is negligible due to low water mobility in the matrix. Gas only sources from the initial gas in fracture, and matrix contribution is negligible. Water is assumed to be slightly compressible, and water compressibility is constant. Water formation volume factor, water viscosity, fracture permeability, and fracture porosity are pressure-dependent variables, which are evaluated at average fracture pressure. While linear, elliptical, and radial flow may be observed for different fracture geometries, which are controlled by xf/h (Zhang and Emami-Meybodi, 2019), we only considered linear flow in the fracture. The gravity segregated relative permeability curves (Clarkson, 2012) are valid for the fracture.

2.2. Flow regime analysis Flow regime analysis is a useful tool to identify what production stages the well has experienced and to extract reservoir and fracture properties from specific flow regimes. The analysis is usually based on the long-term production data by generating the diagnostic plot, which is rate normalized pressure (RNP) and its derivative (DRNP) vs. superposition (material balance) time in log-log scale (Song and Ehlig2

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

104 (b) FR3

FR1 + FR2

DRNP (psia×day/STB)

103 Single-phase and Two-phase Water Transient Linear Flow

102

1 1

101

2 1

100

Two-phase Water BDF

10-1 100

101

102

103

104

105

106

Superposition Pseudotime (psia×day/cp)

Fig. 2. (a) Schematic of the three flow regimes and (b) Two-phase diagnostic plot modified from Zhang and Emami-Meybodi (2018) during early flowback period.

2.2.3. FR3 Two-phase boundary dominated flow (BDF) in the fracture is the most dominant regime for water flow during the early flowback period exhibiting a unit-slope line in the diagnostic plot. Alkouh et al. (2014) proposed a model for FR3 using a simplified mobility and compressibility definition. Williams-Kovacs and Clarkson (2016) developed an “after-breakthrough” model for two-phase fracture depletion in a radial system. Jia et al. (2018) used flowing material balance equation (constant rate solution of water-phase diffusivity equation) to analyze flowback data during FR3. In the present paper, we focus on this flow regime and present a new RTA model for both constant bottomehole pressure (BHP) and constant/variable-rate condition. The behavior of gas production during flowback period is different with water production. As pointed out by Zhang and Emami-Meybodi (2018), the hydraulic fracture is initially saturated with water and very less trapped gas after stimulation job. During the early time of flowback, the produced gas mainly comes from the fracture depletion, and gas influx from matrix can be overlooked especially for shale gas reservoirs. Because the permeability of shale matrix is much smaller than fracture permeability, which makes the duration of fracture depletion flow longer than other types of reservoirs. But in the late flowback period, the fracture pressure drops significantly and most of the gas production comes from matrix influx. In this work, we only focus on the multiphase water and gas flow during the early flowback period. In the next work (Zhang and Emami-Meybodi, 2019), we elaborate the gasphase flow regimes in the fracture and matrix during late flowback period and discuss the coupling of fracture flow with matrix flow in details.

Economides, 2011). By considering the nonlinearity caused by pressure-dependent permeability, Zhang and Emami-Meybodi, 2018 replaced RNP and superposition time with rate normalized pseudopressure and superposition pseudotime, and generated diagnostic plots for long-term production to analyze flowback data. They proposed three possible water-phase flow regimes that may occur for the conceptual model depicted in Fig. 1: (1) single-phase water fracture transient linear flow (FR1), (2) two-phase fracture transient linear flow (FR2), and (3) two-phase boundary dominated flow in fracture (FR3). In this work, we proposed a two-phase diagnostic plot (see Fig. 2b) based on multiphase pseudopressure and pseudotime. The definitions of pseudo-variables used in the diagnostic plot are introduced in Section 2.3. Compared to the traditional single-phase diagnostic plots, the proposed diagnostic plot is applicable to both single- and two-phase flow systems. This understanding is consistent with Clarkson et al. (2019) that using single-phase diagnostic plots to identify flow regimes in multiphase flow systems may lead to erroneous results. Fig. 2a and b show the schematic of these three flow regimes and the modified twophase diagnostic plot of water, respectively.

2.2.1. FR1 Single-phase water transient flow (infinite acting linear flow) is the first flow regime that may be observed during the flowback period. After well stimulation, the hydraulic fracture is mostly saturated with fracturing liquid. Therefore, single-phase water production dominates the fracture system at the beginning, which shows the half-slope straight line in the diagnostic plot. Clarkson and Williams-kovacs (2013) modeled a similar single-phase water transient flow in a radial system. Although the analysis on FR1 could calculate fracture permeability and has been validated by synthetic data, it may not be applicable to field data because the duration of FR1 is often too short.

2.3. Two-phase BDF model Due to the high conductivity of hydraulic fracture, we assumed 1D BDF for the fracture. This work divides basic production conditions during flowback period into constant BHP and constant rate, which can be extended to more realistic variable BHP/rate condition using the superposition principle.

2.2.2. FR2 During the two-phase transient flow, the diagnostic plot shows a half-slope straight line as well. The produced gas is the initial trapped gas inside the fracture network. Although half-slope straight line can also be observed, the single-phase transient flow model cannot be used for FR2 because gas production has a significant impact on total compressibility so that assuming only water expansion dominates flow will cause an enormous error (Zhang and Emami-Meybodi, 2018). Furthermore, FR2 is usually short making it difficult to conduct any meaningful analysis. Nonetheless, the end of the half-slope straight line could be used to identify the beginning of FR3.

2.3.1. Constant BHP We first developed a 1D two-phase flowback model under constant BHP conditions. The water-phase material balance equation for an infinitesimal fracture element is given by:

qρw |x + dx − qρw |x =

d (dx wf hϕf ρw s w ) dt

(1)

where x is the location in the fracture, dx wfh is the control volume, q is 3

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

the water-phase downhole flowrate at a specific location x, sw is water saturation in the fracture, ρw is the water-phase density, wf is the fracture width, h is the fracture height, ϕf is the fracture porosity, and t is the time. Dividing both sides of Eq. (1) by dx, the following partial differential equation (PDE) can be obtained:

∂ (qρw ) ∂ = (wf hϕf ρw s w ) ∂x ∂t

variables method (Yang, 2017): ∞

pp (x , tp) = (ppi − ppwf )

n=1

(2)

μ wi B wi kf krw ∂p kfi μ w B w ∂x

w

( ) qρ wsc



∂x

= (ppi − ppwf )

(

kf krw ∂p w B w ∂x

qw = 2wf h μ

∂x

=0

∑ n=1

Eq. (3) can be

RNP =

where pi is the initial fracture pressure and pwf is the bottomhole pressure. The exponential relationship of pressure-dependent fracture permeability and porosity (Jr Pedrosa, 1986) were considered to include permeability and porosity changes with respect to pressure drop in the fracture:

ϕf = ϕfi g (p) = ϕfi exp[−cf (pi − p)]

(7)

4wf hkfi x f B wi μ wi

tp =

μ wi B wi kfi

ϕfi kfi

∫0

p

kf (p)⋅krw (s w )

b

μ w (p) B w (p)

∫p

dp

mp =

x f B wi μ wi 4wf hkfi

s w μ w ⎛c w + cf + ⎝

1 ∞

kfi ∑ exp ⎛−βn2 ϕ x 2 tp ⎞ fi f n=1 ⎝ ⎠

(pf )

∂x 2

=

(9)

ϕfi ∂pp kfi ∂tp

pp (x , 0) = ppi , pp (0, tp) = ppwf ,

∂x

=0

(16)

ϕfi x f2

(17)

221.73x f B wi μ wi ⎞ wf kfi h

pp (x , 0) = ppi ,

(10)

∂pp (x f , tp )

(15)



(18)



2.3.2. Constant rate Constant rate case in flowback modeling is as important as constant BHP case, not only because it is the basic condition for well testing, but also it can be easily extended to variable rate/BHP case with superposition principle. With the same definition of pseudotime and pseudopressure in Eqs. (8) and (9), we used Eq. (10) in 1D fracture subject to the following initial and boundary conditions:

The governing Eq. (4), initial condition, and boundary condition in Eq. (5) can be linearized to be:

∂2pp



where mp and bp are the slope and intercept of ln(RNP) vs. tp specialty plot, respectively, and all the parameters have the field units.



1 dsw ⎞⋅ϕf sw dpf ⎠

(14)

0.0156kfi

bp = ln ⎛⎜ ⎝

(8)

krw (s w )⋅kf (pf )

t

kfi ⎞ ⎛ exp ⎜−βn2 tp ϕfi x f2 ⎟ ⎠ ⎝

ln(RNP ) = mp tp + bp

where γ is the permeability modulus, kfi is the initial fracture permeability, and ϕfi is initial fracture porosity. We used the pseudopressure (pp) and pseudotime (tp) concepts to linearize the partial differential equation by evaluating pressure-dependent properties at average pressure ( pf ) and average water saturation (s w ) in the fracture:

pp =

(13)

where qw is the water-phase flow rate from the whole fracture at surface condition, RNP = (ppi-ppwf)/qw is the rate normalized pseudopressure, and Bwi is the water formation volume factor evaluated at initial pressure. Applying the long-term approximation described by Wattenbarger et al. (1998) to Eq. (15), a linear relationship between ln(RNP) and tp (i.e., constant BHP specialty plot) can be obtained:

(5)

(6)

x=0

Eq. (14) can be rearranged by defining rate normalized pressure (RNP):

(4)

kf = kfi f (p) = kfi exp[−γ (pi − p)]

)

to (13):



∂p (x f , t )

xf

x=0 ∞

qw = (ppi − ppwf )

where cw is the water compressibility, cf is the fracture compressibility, p is the pressure in the fracture, and μw is the water viscosity. Eq. (4) is subject to the following initial and boundary conditions:

p (x , 0) = pi , p (0, t ) = pwf ,

β x

kfi ⎞ ∂ sin n 2 ⎛ exp ⎜−βn2 tp βn ∂x ϕfi x f2 ⎟ ⎠ ⎝

Applying Darcy's water surface flowrate for the whole fracture

where ρwsc is the water density at the standard condition and Bw is the water formation factor. Using the definitions of water and fracture compressibility (Jia

ϕf s w ∂ ⎛ krw kf ∂p ⎞ ⎛⎜c w + cf + 1 ds w ⎞⎟ ∂p ⎜ ⎟ = Bw ⎝ s w dp ⎠ ∂t ∂x ⎝ B w μ w ∂x ⎠

∑ n=1

ϕf s w ρwsc ⎛ 1 dϕf ⎞ ϕf s w ρwsc 1 ds w ⎤ ⎡ ϕf s w ρwsc ⎛ 1 dρw ⎞ ⎛ ⎞ = wf h ⎢ + ⎜ ⎟ ⎜ ⎟ + B w ⎝ ρw dp ⎠ B w ⎜ ϕf dp ⎟ B w ⎝ s w dp ⎠ ⎥ ⎝ ⎠ ⎣ ⎦ ∂p (3) ∂t

et al., 2018) and Darcy's water rate q = written as:

(12)

x=0

Bw

krw kf ∂p wf h μ ∂x , w

kfi ⎞ x 2 ⎛ sin ⎛⎜βn ⎞⎟ exp ⎜−βn2 t 2 p⎟ βn x ϕ f⎠ fi x f ⎠ ⎝ ⎝

where βn = (2n – 1)π/2. Taking the derivative of Eq. (12) with respect to x and evaluating it at x = 0 gives:

Applying chain rule to the accumulative term (right-hand side) and ρ substituting ρwf = Bwsc , Eq. (2) can be expanded as:





∂pp (0, tp ) ∂x

=−

qw B wi μ wi ∂pp (x f , tp ) , =0 2kfi wf h ∂x

(19)

The analytical solution of Eq. (10) subject to Eq. (19) is given by Miller (1962):

(11)

Eq. (10) subject to Eq. (11) can be solved using the separation of 4

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

relationship but between RNP and tsp (i.e., variable-rate specialty plot), where the slope and intercept of specialty plot are same as constant rate case (Eqs. (22)–(24)).

pp (x , tp) ⎡ qw x f μ wi B wi ⎢ 1 kfi tp x2 x − + = ppi − − ⎢ + 2kfi wf h ⎢ 3 xf 2x f2 ϕfi x f2 ⎢ ⎣ ∞

2 cos



( ) exp ⎛⎝ nπx xf

−(nπ )2kfi tp



⎞⎤ ⎠⎥ ⎥ ⎥ ⎥ ⎦

(20) 2.4. Calculation of average pressure

We evaluate Eq. (20) at x = 0 and substitute pp(0,tp) = ppwf:

ppi − ppwf

⎡ qw x f μ wi B wi ⎢ 1 kfi tp − = ⎢ + 2kfi wf h ⎢ 3 ϕfi x f2 ⎢ ⎣

−(nπ )2kfi tp



∑ n=1

2 exp ⎛ ϕ x 2 ⎝ fi f (nπ )2 ⎜

⎞⎤ ⎠⎥ ⎥ ⎥ ⎥ ⎦

Average pressure in the fracture is a critical parameter used to evaluate pressure dependent properties and calculate pseudopressure and pseudotime in the flowback model. Material balance equation is a simple and popular way to obtain average pressure, which requires both water and gas production data as the input (Moghadam et al., 2011). Behmanesh et al. (2018a) derived an explicit expression to calculate average pseudopressure within the distance of investigation, but they didn't provide the solution of average pressure. This section proposes a diffusivity equation approach to calculate average pressure in the fracture under constant BHP and constant/variable rate conditions and also derived a two-phase material balance equation in a tank system.



(21)

Taking the long-term approximation (Wattenbarger et al., 1998) and introducing RNP definition from the constant BHP model, a linear relationship between RNP and tp (i.e., constant rate specialty plot) can be obtained from Eq. (21) in field unit:

RNP = mq tp + bcq

(22)

2.807B wi μ wi wf hx f ϕfi

(23)

mq =

bq =

2.4.1. Diffusivity equation approach The new method is based on the solution of pseudopressure from the diffusivity equation. The volumetric average pseudopressure under constant BHP, constant rate, and variable rate conditions are calculated by taking the integral of Eqs. (12), (20), and (26) with respect to x:

147.82x f B wi μ wi wf hkfi

(24)

where mq and bq are the slope and intercept of RNP vs. tp specialty plot, respectively.

ppf (tp)|constant

2.3.3. Variable rate Here we extend the constant rate solution of pseudopressure to variable-rate condition by applying superposition principle (Dake, 1983). By defining superposition pseudotime tsp in Eq. (25), the solution of pseudopressure under variable rate condition is given in Eq. (26). N

tsp =

∑ j=1

pp (x , tp) = ppi −

2kfi wf h

⎡1 kfi tsp x2 x ⎢ + − + − xf ϕfi x f2 2x f2 ⎢3 ⎣

BHP

=

x ∫0 f pp (x , tp)|constant BHP dx

xf ∞

= (ppi − ppwf )

∑ n=1

ppf (tp)|constant

⎡ (qw, j − qw, j − 1 )(tp, N − tp, j − 1 ) ⎤ ⎢ ⎥ qw, N ⎣ ⎦ qw, N x f μ wi B wi

(28)

According to Eqs. (16), (22) and (28), fracture properties, such as fracture half-length (xf) and initial fracture permeability (kfi), are part of the slope and intercept of specialty plot. The procedure to obtain these attributes is discussed in Section 2.5.



ϕfi x 2f

(nπ )2

n=1

RNP = mq tsp + bcq

q

kfi ⎞ 2 ⎛ exp ⎜−βn2 tp βn2 ϕfi x f2 ⎟ ⎝ ⎠

x ∫0 f pp (x , tp)|constant q dx

B wi qw tp ⎞ ⎛ = ppi ⎜1 − ppi Vfi ⎟ ⎝ ⎠

(30)

⎛ B wi ∑ (qw, j − qw, j − 1)⋅(tp, N − tp, j − 1 ) ⎞ ⎜ ⎟ j=1 = ppi ⎜1 − ⎟ ppi Vfi ⎜ ⎟ ⎝ ⎠

(31)

=

xf

(25)

ppf (tp )|variable

×

q

(29)

=

x ∫0 f pp (x , tp )|variable q dx

xf N





(qw, j − qw, j − 1 ) 2 cos qw, N

n=1

(

nπx

xf

)

(nπ )2

2 ⎛ (nπ ) kfi (tp, N − tp, j − 1 ) ⎞ ⎤ ⎥ exp ⎜− ⎟⎥ ϕfi x f2 ⎝ ⎠⎦

where ppf is the average pseudopressure in fractures. According to the definition of pseudopressure in Eq. (8), real pressure can be converted from pseudopressure using the following relation:

(26)

where qw,j is the water-phase flow rate at the surface condition in the jth flow rate interval, tp,j is the pseudotime at the jth flow rate interval. Similarly, we evaluate Eq. (20) at x = 0 and substitute pp(0,tp) = ppwf:

dpp

ppi − ppwf qw, N x f μ wi B wi ⎡ 1 kfi tsp = + − 2kfi wf h ⎢ 3 ϕfi x f2 ⎣



∑ n=1

2 ⎛ (nπ ) kfi (tp, N − tp, j − 1 ) ⎞ ⎤ exp ⎜− ⎟⎥ ϕfi x f2 ⎝ ⎠⎦

μ wi B wi kf (p) krw (s w ) dp kfi μ w (p) B w (p) dt

(32)

kfi μ w (p) B w (p) dpp dp = dt μ wi B wi kf (p) krw (s w ) dt

(33)

dt

2(qw, j − qw, j − 1 )

=

qw, N (nπ )2

p = pi + (27)

kfi μ wi B wi

∫0

t

μ w (p) B w (p) dpp dτ kf (p) krw (s w ) dt

(34)

Assuming that pseudo average pressure can be approximated with average pseudopressure, Eq. (34) can be rewritten to calculate average pressure based on average pseudopressure:

Introducing the definition of RNP and applying the same long-term approximation to variable-rate solution, we can obtain the linear 5

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

Gp = Gi − Gr

pf = pi + ≈ pi +

kfi μ wi B wi kfi μ wi B wi

∫0

t

∫0

t

μ w (pf ) B w (pf ) dppf kf (pf ) krw (s w ) dt μ w (pf ) B w (pf ) dppf kf (pf ) krw (s w ) dt

where Gi is the initial gas storage in the fracture and Gr is the remaining gas in the fracture, all at the surface condition. The definition of Gi and Gr are given by:

dτ dτ

Gi =

(35)

where ppf is the pseudo average pressure in the fracture. With the relation between pressure and pseudopressure in Eq. (35), we take the derivative of Eqs. (29), (30), and (31) with respect to t:

dppf dt

2ζ (ppi − ppwf ) krw kf (pf )

=−

s w kfi μ w ⎛c w + cf + ⎝

constant BHP

⎛ 2 −β ∑ exp ⎜ 2n x ⎜ n=1 ⎜ f ⎝ ∞

∫0

dppf (tp)|constant

t

dt dppf (tp )|variable

q



(pf )

pf = pi − Bg =

(36)

krw (s w ) kf (pf ) s w μ w ⎛c w + cf + ⎝

B wi qw, N ϕfi ⋅ Vfi kfi

1 dsw ⎞ ϕf sw dpf ⎠



⎡ ⎢

BHP

∑ exp ⎢− ∫0

n=1

τ2

⎢ ⎣

pf (tp)|constant

q

∫0

= pi −

t

(pf )

0.0127ζ (ppi − ppwf ) B w (pf ) dsw ⎞ ϕf dpf



⎤ ⎥ dτ1 ⎥ dτ2 1 dsw ⎛ ⎞ kfi μ w c w + cf + s dp ϕf (pf ) ⎥ w f ⎠ ⎝ ⎦

pf (tp )|variable q = pi −

∫0 ∫0

t



t

dsw ⎞ ϕf dpf



dτ (pf )

5.615qw, N B w (pf ) ϕfi Vfi ⎛s w c w + s w cf + ⎝

dsw ⎞ ϕf dpf



(40)

dτ (pf )

(41)

(42)

B w = B wi exp[c w (pi − pf )] ≈ B wi [1 + c w (pi − pf )]

(46)

Vfi − ((Gi − Gp ) Bg + (Wfi − 5.615Wp ) B w ) Vfi cf

(47)

psc ZT Tsc p

(48)

w

(Wfi − 5.615Wp ) B w Vfi [1 − cf (pi − pf )]

B wi

f

evaluated at the average pressure and water saturation in the fracture because pressure changes in distance can be ignored due to the high conductivity of the fracture (Zhang et al., 2013). According to the definition of pseudopressure in Eq. (8), the initial pseudopressure drawdown (ppi-ppwf) changes with the lower limit of the integral pwf in constant/variable-rate flowback models. To calculate (ppi-ppwf), we have to obtain pf first from Eqs. (41) or (47). Then s w can be estimated corresponding to each pf from Eq. (42). By interpolating the relation between s w and pf , we can evaluate (ppi-ppwf) at each pwf point. Unlike constant/variable-rate case, (ppi-ppwf) is a constant for constant BHP condition, which is independent of average fracture pressure. To obtain this constant, we have to know how s w changes from pwf to pi. This paper uses a linear function to denote the relation between s w and pf as a simplification:

We use the simplified material balance equation proposed by King (1993) to calculate average fracture water saturation in Eqs. (39)–(41).

sw =

Vfi s wi

2.5.1. Calculation of pseudotime and pseudopressure The specialty plots described in Sections 2.3 require the calculations of pseudopressure and pseudotime. Pseudotime tp is a function of average pressure in the fracture and can be calculated by substituting either of the Eqs. (39), (41), and (47) into Eq. (9). Choosing the appropriate average pressure formula depends on the well production condition and the availability of two-phase flowback data, which will 1 ds be discussed in Section 2.5.3. The s dpw term in tp formulation can be

(pf )

(39)

5.615qw B w (pf ) ϕfi Vfi ⎛s w c w + s w cf + ⎝

(45)

Based on the proposed two-phase flowback model and the average fracture pressure calculation methods, we developed an analysis technique to interpret early-time flowback data and estimate fracture properties such as xf, kfi, and Vfi. The analysis technique consists of three steps: (1) calculation of pseudo-variables, (2) identification of two-phase BDF regime, and (3) calculation of fracture properties using successive substitution method.

(38)

0.0156ζ (2n − 1)2kf (pf )

= pi −

Bg

2.5. Analysis technique

1 dsw ⎞ ϕf sw dpf ⎠

μ wi B wi ⎛s w c w + s w cf + ⎝

Vf − Vwr

(37)

where ζ = kfi/xf2 is the fracture coefficient. Substituting Eqs. (36), (37), and (38) into Eq. (35), we arrive at the expression of average pressure in the fracture under constant BHP, constant rate, and variable rate conditions in field unit:

pf (tp)|constant

, Gr =

where psc is the pressure at standard condition and Tsc is the temperature at standard condition.

(pf )

krw (s w ) kf (pf ) s w μ w ⎛c w + cf + ⎝

Bgi

where Vwr is the remaining water volume in the fracture. Substituting Eqs. (45) and (46) into Eq. (44), we obtain an implicit expression for average fracture pressure which can be solved by secant method:



⎞ dτ ⎟ ⎟ 1 ds s w μ w ⎛c w + cf + s dpw ⎞⋅ϕf (pf ) ⎟ w f ⎠ ⎝ ⎠

=−

dt

1 dsw ⎞ ϕf sw dpf

sgi Vfi

Vf = Vfi (1 − cf (pi − pf )), Vwr = (Wfi − Wp ) B w , Wfi =

krw (s w )⋅kf (pf )

B wi qw ϕfi =− Vfi kfi

q

(44)

(43)

where Wfi is the initial water storage in the fracture at the surface condition, Wp is the cumulative water production at the surface condition, and Vfi = 2ϕfixfwfh is the initial pore volume of the fracture. With the average water saturation, dsw term can be approximated

s w = Apf + B

dpf

A=

using backward difference so that we can calculate average pressure in the fracture iteratively.

s wi − s wf pi − pwf

(49)

, B = s wi − Api

(50)

where, A and B are the coefficients of linear correlation, swf is the average water saturation inside fracture when pf reaches pwf and it depends on the estimated ultimate water recovery in the end of flowback. Here we used cumulative water production at the end of flowback Wp,max to approximate the ultimate water cumulative production so that swf can be approximated as:

2.4.2. Material balance approach Here, the fracture network is treated as a tank system in the early flowback without gas matrix influx acting on fracture flow. The material balance equation based on gas fracture storage can be written as (King, 1993): 6

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

(Wfi − 5.615Wp |wf ) B wf

s wf =

Vfi [1 − cf (pi − pwf )]

kfi, p =

(51)

By substituting Eqs. (49), (43), and (6) into Eq. (8), we arrive at the initial pseudopressure drawdown (ppi-ppwf) for constant BHP condition:

ppi − ppwf

μ B wi = wi kfi

pi

kfi exp[−γ (pi − p)]⋅(Ap + B )

wf

μ w (p) B wi exp[c w (pi − p)]

∫p

dp

xf ,q = (53)

kfi, q =

2.5.2. Identification of two-phase BDF regime To perform analysis with the proposed flowback model, the data interval in BDF should be identified first. According to flow regime analysis presented in Section 2.2, two-phase BDF flow (FR3) starts right after the end of two-phase transient flow in the fracture (FR2). Therefore, a diagnostic plot should be used to identify the end of FR2, which is the beginning of FR3. The new initial and boundary conditions for FR2 under constant BHP condition and constant rate condition are given by:

∂pp (0, tp ) ∂x

=−

qw B wi μ wi , pp (x → ∞, t ) = ppi 2kfi wf h

Vfi, q

RNP|constant

BHP

q

=

=

62.536B wi μ wi wf h ϕfi kfi

39.825B wi μ wi wf h ϕfi kfi

39.825B wi μ wi wf h ϕfi kfi

(55)

tp (56)

tp (57)

(63)

147.82x f B wi μ wi bq wf h

5.614B wi μ wi = 2ϕfi x f wf h = mq

(64)

(65)

The two workflows cover both constant BHP and constant/variablerate flowback model which are capable of analyzing complex flowback data in reality. Workflow I uses material balance approach to calculate average pressure in the fracture, which requires both water and gasphase production data as the inputs. Whereas, Workflow II only requires water-phase production data and releases the need for gas data for flowback analysis by applying the new average pressure calculation method. For the field application, it depends on the availability and quality of gas-phase production data to choose the most appropriate workflow. If the well is producing under constant BHP condition, Workflow III is an alternative way to analyze flowback data which uses the diffusivity equation approach to calculate average pressure in the fracture. It follows nine steps:

tsp (58)

2.5.3. Successive substitution method for xf, kfi, and Vfi This section presents three workflows using successive substitution method to extract xf, kfi, and Vfi from flowback data. For the constant BHP flowback model in Eq. (16), one can use the slope mp and intercept bp of ln(RNP) vs. tp plot to calculate xf,p, kfi,p, Vfi,p, and ζp using the following expressions:

3.459B wi μ wi mp ϕfi exp(bp ) wf h

2.807B wi μ wi mq wf hϕfi

1. Guess an initial value for Vfi, 2. Obtain average fracture pressure pf from Eqs. (47) for Workflow I or from Eq. (41) for workflow II. Calculate (ppi-ppwf) and evaluate tp at pf , 3. Generate two-phase diagnostic plot, DRNP vs. tsp in log-log scale, to pick the time when half-slope straight line and unit-slope straight line intersects, which is the start of BDF regime (FR3), 4. Generate specialty plot ln(RNP) vs. tp if using constant BHP flowback model or RNP vs. tsp for constant/variable-rate flowback model. Then, extract slope and intercept of the straight line after the end of half-slope interval in step 3, 5. Calculate a new value of Vfi using Eq. (61) or (65), 6. Repeat step 2 to 4 until the relative error of Vfi meets the tolerance we set, 7. Determine xf and kfi using Eqs. (59)–(60), or (63)–(64).

We defined the derivative of RNP (DRNP = dRNP/dlntsp) with the new definitions of pseudopressure and pseudotime in this work. The two-phase transient flow solution under general production condition in Eq. (58) reveals that log-log DRNP vs. tsp plot shows a half-slope straight line denoting FR1 + FR2. The two-phase boundary dominated flow solution under variable rate condition in Eq. (28) indicates a unitslope straight line in the same plot denoting FR3. Hence, we used loglog DRNP vs. tsp plot as the new two-phase diagnostic plot to identify the flow regimes.

xf ,p =

(62)

(54)

Eq. (57) can be easily extended to variable rate condition with the superposition principle by replacing tp with tsp:

RNP|variable q =

(61)

However, the calculation of average pressure in Eqs. (39), (41), and (47) requires Vfi as a given input, which is usually unknown in field application. Fig. 3 shows three workflows using successive substitution method to estimate Vfi iteratively, and then calculate xf and kfi. Workflow I and II have the same calculation procedures. The only difference is the method used to calculate average pressure in the fracture. They both follow seven steps:

Wattenbarger et al. (1998) gave the solutions of Eq. (10) together with the new outer boundary conditions. Rearranging the solutions in terms of RNP in the field unit gives:

RNP|constant

(60)

For the constant rate or variable rate flowback model in Eqs. (22) and (28), the slope mq and intercept bq of RNP vs. tp and RNP vs. tsp plot can be used to calculate xf,q, kfi,q, and Vfi,q using the following expressions:

Apwf + B 1 ⎧ A [1 − e−(γ + cw )(pi − pwf ) ] ⎫ (Api + B ) − (γ + c )(p − p ) − w i wf ⎬ γ + cw ⎨ γ + cw e ⎩ ⎭

pp (x , 0) = pi ,

6.918B wi μ wi mp exp(bp )

ζ p = 64.155mp ϕfi

(52)

ppi − ppwf

pp (x , 0) = ppi , pp (0, t ) = ppwf , pp (x → ∞, t ) = ppi

mp ϕfi (exp(bp ) wf h)2

Vfi, p = 2ϕfi x f wf h =

If we assume constant water viscosity and apply integration by parts, we can obtain the analytical expression for Eq. (52):

=

2 2 767.59Bwi μ wi

1. Guess an initial value for Vfi and ζ, 2. Obtain average fracture pressure pf from Eq. (39), calculate (ppippwf) and evaluate tp at each pf , 3. Generate two-phase diagnostic plot, DRNP vs. tsp in log-log scale, to pick the time when half-slope straight line and unit-slope straight line intersects, which is the start of BDF regime (FR3),

(59) 7

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

Fig. 3. Three workflows of the successive substitution method to evaluate Vfi, xf, and kfi. Workflow I applies the material balance approach to calculate average pressure in the fracture, whereas Workflow II and III use the diffusivity equation approach. Workflow III is an alternative method for constant BHP condition.

4. Generate specialty plot ln(RNP) vs. tp and extract slope and intercept of the straight line after the end of half-slope interval in step 3, 5. Calculate a new value of ζ using Eq. (62), 6. Repeat step 2 to 5 until the relative error of ζ meets the tolerance we set, 7. Calculate a new value of Vfi using Eq. (61), 8. Repeat step 2 to 7 until the relative error of Vfi meets the tolerance we set, 9. Determine xf and kfi using Eq. (59)–(60).

Table 1 Input parameters used in numerical simulation.

Workflow III only requires water-phase production data as the input, like Workflow II, but it is only valid for constant BHP case. Therefore for flowback analysis under constant BHP condition, Workflow I, II and III all work.

Parameter

Value

Initial pressure, pi(psi) Reservoir temperature, T (°F) Specific gravity, SGg Water compressibility, cw(psi−1) Initial fracture porosity, ϕfi (%) Fracture compressibility, cf (psi−1) Reservoir thickness, h (ft) Fracture width, wf (ft) Initial water formation volume factor, Bwi(bbl/STB) Water viscosity @pi, μw(cp) Water density @pi, ρw(lb/ft3) Fracture half-length, xf(ft) Initial fracture permeability, kfi(mD)

5000 160 0.65 8 × 10−6 50 5 × 10−6 100 0.02 1.03 0.6 62 1000 800

3. Model validation and results calculation methods against the numerical simulation for Case 1, Case 2, and Case 3. These cases are all under constant BHP condition so that the diffusivity equation approach for constant BHP and variable rate case, and material balance approach can be compared against each other and the numerical simulation results by considering Vfi as a known parameter. Fig. 5 shows the comparison of results for the three cases. The results show that there is a good agreement between numerical simulation and analytical methods. Although the nonlinearity caused by two-phase flow and pressure dependent fracture permeability is increasing from Case 1 to Case 3, the good match of average pressure for all the cases indicates that the nonlinearity does not affect the accuracy of average pressure calculation. The average pressure values obtained from the variable-rate method are sensitive to the quality and trend of production data. Further investigation is required to improve the calculation of average pressure using variable-rate method.

To validate the developed flowback model, average pressure calculation methods, and the proposed workflows, we conducted several numerical simulations for one fracture stage as shown in Fig. 1 using IMEX—CMG to generate production data. A 1D flow system with 1000 uniform grids was considered as the physical model. To model fracture flow with high conductivity, we set the grid block width to 0.02 ft. and permeability to 800 md. Considering that the early gas production period only lasts for a short time (Adefidipe et al., 2014), we conducted the simulation for 10 days, which are long enough to show boundary dominated flow during the flowback period. The permeability modulus values were obtained from the laboratory measurements from Ozkan et al. (2010). The input parameters used in the numerical simulations are shown in Table 1. A list of simulation runs with different degrees of nonlinearity and production conditions are given in Table 2. We used the gravity segregated relative permeability curve (Clarkson, 2012) for the fracture, as shown in Fig. 4a. Gas viscosity values and formation volume factor vs. pressure is presented in Fig. 4b.

3.2. Validation of flowback model After validation of average pressure calculation methods, we tested the accuracy of the two-phase flowback model against the production data obtained from the numerical simulation, by considering Vfi as a known parameter. Material balance approach was used to calculate

3.1. Validation of average pressure calculation First, we tested the accuracy of the average fracture pressure 8

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

Table 2 The numerical cases used for validation. Case #

Initial water saturation, swi

Fracture permeability modulus, γ(psi−1)

Production condition

Bottom-hole pressure, pwf(psi)

Water-phase surface flowrate, qw(bbl/ d)

1 2 3 4 5 6

0.95 0.9 0.7 0.95 0.9 0.7

1 × 10−5 1 × 10−5 5 × 10−4 1 × 10−5 1 × 10−5 5 × 10−4

Constant Constant Constant Constant Constant Constant

700 700 700 – – –

– – – 2 2 2

BHP BHP BHP rate rate rate

constant BHP model has much more error than constant rate model, which is expected. The error of both constant BHP model and constant rate model increases with the extent of nonlinearity. Jia et al. (2018) observed a similar behavior where fracture properties calculation errors increased with the nonlinearity caused by the strong water saturation and pressure gradients. Thus, the constant BHP flowback model is only valid when the nonlinearity is not significant, whereas the constant rate model always works. For Case 3 with the greatest extent of nonlinearity under constant BHP condition, the variable rate model provides more accurate results compared to constant BHP model. To further validate our model, we compared the analysis results with the two-phase flowing material balance model proposed by Xu et al. (2016) and Behmanesh et al. (2018b) in Table 3 and listed the differences of the three models in Table 4. We extended the original model of Xu et al. (2016) from constant rate conditions to include variable-rate conditions by introducing superposition pseudotime. The estimated xf from Xu et al. (2016) shows to be accurate for all the numerical cases, but the calculated kfi has much more error than our flowback model in this work, especially for Case 1–3. It is because their model assumes constant kf and ϕf, whereas in the numerical cases kf and ϕf are pressure-dependent properties. What's more, the pressure and saturation gradient in Case 1–3 are greater than Case 4–6 as shown in Fig. 9, which causes more error using the approach of Xu et al. (2016). We also extended the two-phase RTA model of Behmanesh et al. (2018b) from radial system to linear system by adopting the bpss term in Zhang and Emami-Meybodi (2019), and to consider pressure-dependent permeability. The analysis results from Behmanesh et al. (2018b) in Table 3 show good accuracy in calculating both kf and xf for all the numerical cases. Compared to these two flowback models, our proposed flowback model only requires water production data and does not require gas production data during flowback analysis.

average pressure in this section. We analyzed Case 1, 2, and 3 with constant BHP flowback model and Case 4, 5, and 6 with constant rate model. The two-phase diagnostic plot and specialty plot for all the six cases are shown in Figs. 6 and 7. We also analyzed flowback data in Case 3 with the variable rate flowback model. The diagnostic plot and specialty plot are shown in Fig. 8. The two-phase diagnostic plot for each case shows a clear half-slope straight line indicating the transient flow, which is expected based on Eqs. (56)–(58). We picked the intersection point between the half-slope and unit-slope straight line as the start of two-phase boundary dominated flow. After identifying the flow regimes, the corresponding data interval is selected to conduct flowback RTA in specialty plot. Based on the analytical two-phase BDF solution given by Eq. (16) for constant BHP condition and Eq. (22) for constant rate condition, the specialty plot is expected to show a straight line during the BDF. But the results show that only the specialty plots in Case 4–6 have the perfect straight line, whereas Case 1–3 are curved. That is because, applying pseudotime to linearize the governing equation has such an assumption that all the pressure dependent properties, such as sw, kf, and ϕf, are evaluated at average fracture pressure ( pf ) and average water saturation (s w ). Although these properties vary along the fracture, an average value for each of these properties is selected to represent all the values in the entire domain. The pressure profile under constant BHP condition is very sharp at the beginning and drops unevenly in Fig. 9a. Whereas the pressure profile under constant rate condition is gentler and drops with the same rate everywhere during BDF in Fig. 9b. Therefore, the less error of interpretation results and the less curved the specialty plot are expected in constant rate flowback model. With the slope and intercept of the specialty plot, we used Eqs. (59)–(65) to calculate xf and kfi for the six cases. For constant BHP cases, (ppi-ppwf) is estimated with Eq. (53). The calculated results and relative errors are compared with the inputs in the simulator in Table 3. The results show that the obtained fracture properties using constant rate model are very close to the set values in the numerical simulation and the relative errors are all < 10%, which proves the accuracy. Although the (ppi-ppwf) calculation is proved to be accurate enough, the

3.3. Validation of analysis workflow In this section, we validated the proposed workflows of the successive substitution method in Fig. 3. First, we validated Workflow I 1.4

1.0

0.035

(b)

(a)

Relative Permeability

0.6 Water Gas 0.4

0.2

0.030

1.0 0.025

0.8 0.6

0.020

0.4

0.0 0.0

0.015

Gas formation volume factor Gas viscosity

0.2 0.0 0.2

0.4 0.6 Water Saturation

0.8

0.010

0

1.0

1000

2000

3000 4000 Pressure (psia)

5000

6000

Fig. 4. (a) Relative permeability curves for the fracture and (b) gas viscosity and formation volume factor vs. pressure. 9

Gas Viscosity (cp)

Gas Formation Volume Factor

1.2 0.8

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

5000

5000

5000

(a)

(b) 4000

3000

2000

1000

0 0.001

Numerical simulation Constant BHP Variable rate Material balance 0.01

0.1

1

4000 Average Pressure (psia)

Average Pressure (psia)

Average Pressure (psia)

4000

(c)

3000

2000

1000

10

0 0.001

Time (day)

0.1 Time (day)

1

2000

1000

Numerical simulation Constant BHP Variable rate Material balance 0.01

3000

10

0 0.001

Numerical simulation Constant BHP Variable rate Material balance 0.01

0.1

1

10

Time (day)

Fig. 5. Comparison of the diffusivity equation approach and material balance approach for (a) Case 1, (b) Case 2, and (c) Case 3.

15 hydraulic fractures. The input parameters are consistent with those used by Xu et al. (2016) and Anderson et al. (2013) for the same field case. The fracture width used in this study adopts the primary fracture network width in the field case analyzed by Jia et al. (2018). The fracture permeability modulus refers to the empirical range of 1 × 10−4–1 × 10−3 psi−1 from laboratory measurements (Ozkan et al., 2010) and is set to be 5 × 10−4 psi−1 for this case. Fig. 12 shows production history of water flowrate, gas flowrate and bottomhole pressure for this MFHW during flowback and long-term production period. After hydraulic fracturing, this well was flowed back for > 24 days and then produced for 78 days. In the flowback period, water is dominantly produced with a high flowrate and then quickly decreases due to the fracture depletion. When it comes to long-term production, water flowrate is negligible. Whereas gas is mainly produced and matrix contribution of gas influx becomes more and more significant Table 7. We applied the proposed two-phase water flowback model to analyze flowback data and used the diffusivity approach under variable rate condition to calculate the average pressure in the fracture. We started the Workflow II with an initial guess of fracture volume Vfi = 8 × 105 ft3 and generated the average pressure in the fracture, two-phase diagnostic plot, and specialty plot for water during each iteration. The tolerance of the relative error between Vfi and Vfi,new is set to be 0.5%. After convergence, the essential plots of flowback analysis is shown in Fig. 13. Fig. 13a shows the average pressure in the fracture determined from the diffusivity approach and the average fracture permeability calculated using Eq. (6) during flowback. Fig. 13b shows the saturation vs. pressure path obtained from Eq. (42), which is later extrapolated using exponential function to calculate pseudopressure. The two-phase diagnostic plot of water is plotted in Fig. 13c to identify FR2 and FR3, which shows a half-slope and unit-slope straight line respectively. The flow regimes and sequence for this field case are consistent with those in Fig. 2. With identified FR3, the specialty plot is generated in Fig. 13d to calculate xf and kfi based on the slope and intercept of straight line. The flowback analysis results are listed in Table 8. We also analyzed the long-term production data to extract the fracture attributes in the late time using the RTA technique for transient linear flow and history matching approach. Considering that gas is dominantly produced in the long-term period with water production negligible, we first use the single-phase gas formation linear flow model (Clarkson, 2013) to analyze the long-term production data. The gasphase diagnostic plot in Fig. 14a shows a clear half-slope straight line which illustrates the transient linear flow in the shale matrix. The matrix permeability for this MFHW is in a range of 5 × 10−6–1 × 10−4

with Case 6. According to Table 1, Vfi used in the numerical simulation is 2000 ft3. However, here Vfi is an unknown parameter and an initial guess (i.e., 1500 ft3) is needed. Average fracture pressure is obtained first to generate the diagnostic plot and specialty plot for constant rate flowback model. The slope and intercept of specialty plot were used in Eq. (65) to update Vfi = 1635 ft3. The relative error between the initial guess Vfi = 1500 ft3 and output Vfi = 1635 ft3 is 8.3%, which is greater than the acceptable 5% tolerance. Then, Vfi is updated with the new calculated value 1635 ft3. After five iterations, the relative error of Vfi between two consecutive iterations was less than the acceptable tolerance with the final Vfi = 2125 ft3. Using the estimated Vfi, xf and kfi were calculated to be 1063 ft. and 728 md, respectively, with the relative error of < 10%. These calculations are tabulated in Table 5. The calculated average pressure and specialty plot gradually converge to the numerical simulation during the five iterations as shown in Fig. 10. Workflow II follows a similar procedure as Workflow I, with the difference in average pressure calculation. Therefore, the above-mentioned validation is valid for Workflow II. We validated Workflow III with the flowback data generated from Case 1. With the set parameters Vfi = 2000 ft3 and ζ = 0.0008 md/ft2 in the numerical simulation, we started the workflow with the initial value of 1000 ft3 for Vfi and 4 × 10−4 md/ft2 for ζ. The average pressure in the fracture is obtained at first to generate the diagnostic plot and specialty plot for constant BHP flowback model. With the slope and intercept of specialty plot, we used Eqs. (61) and (62) to calculate Vfi to be 1127 ft3 and ζ to be 4.27 × 10−4 md/ft2. The relative error of Vfi and ζ between input and output are 11.3% and 6.4%, respectively. A smaller tolerance of 0.3% is required in Workflow III to obtain a good match of average pressure and accurate interpretation results. Then, Vfi and ζ are updated. After 69 iterations, both relative errors meet the acceptable tolerance. Using the estimated Vfi, xf and kfi were calculated to be 1137 ft. and 875 md, respectively, with the relative error of < 10%. Five of the 49 calculations are tabulated in Table 6 and shown in Fig. 11. 4. Field application We applied the flowback model to an MFHW drilled and completed in the Muskwa and Otter Park members of the Horn River Shale in part of British Columbia, Canada. This field example was analyzed by Xu et al. (2016) to estimate the initial effective fracture volume. The purpose of analyzing this field case is to verify the practical use of the new two-phase diagnostic plot and flowback model, which could respectively identify the water flow regimes in Fig. 2 and quantitatively evaluate the fracture half-length and permeability. Table 7 shows the real field data from the Horn River MFHW with 10

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

104

107 Case 1-(b)

Case 1-(a) 0.19 day

106 103

BDF

RNP (psia×day/STB)

DRNP (psia×day/STB)

Transient Flow (TF)

0.19 day 1 1

2

102

1

TF

BDF

105 104 103 102

Numerical simulation Straight line 101 103

Numerical simulation Straight line 101

104

105

106

105

0

Superposition Pseudotime (psia×day/cp)

2x105

3x105

104

107 Case 2-(a)

Case 2-(b) 0.22 day

106 BDF

RNP (psia×day/STB)

DRNP (psia×day/STB)

Transient Flow (TF)

103

0.22 day 1 2

102

4x105

Pseudotime (psia×day/cp)

1

1

TF

BDF

105 104 103 102

Numerical simulation Straight line

Numerical simulation Straight line 101 103

101 104

105

106

105

0

Superposition Pseudotime (psia×day/cp)

2x105

3x105

4x105

Pseudotime (psia×day/cp)

104

105 Case 3-(a)

Case 3-(b)

Transient Flow (TF)

103

BDF

104 RNP (psia×day/STB)

DRNP (psia×day/STB)

0.14 day

0.14 day 1 2

1

102 1

TF

BDF

103

102 Numerical simulation Straight line

Numerical simulation Straight line 101 103

101 104

105

106

0

Superposition Pseudotime (psia×day/cp)

5.0x104

105

1.5x105

Pseudotime (psia×day/cp)

Fig. 6. Two-phase (a) diagnostic plot and (b) specialty plot for Case 1–3 using constant BHP flowback model.

24 days flowback. The average kf = 212 md from long-term production data analysis is slightly below the estimated ultimate kf = 394 md from flowback analysis, which is reasonable according to (6). Because the ultimate average pressure in the fracture during flowback period is lower than that in long-term production. The decrease in xf and kf is mainly caused by the partial closure of hydraulic fractures.

md (Anderson et al., 2013), which leads to the interpreted xf in the range of 103–462 ft. based on the slope of specialty plot in Fig. 14b. Then, we use history matching approach to analyze the long-term production data with the single-phase gas Horizontal-MultifracturedSRV model in IHS Harmony software. Fig. 14c and d show a good match of gas flowrate and bottomhole pressure history during long-term production. The estimated xf and kf are 430 ft. and 212 md, respectively. Table 8 summarized the flowback analysis results using the new flowback model and long-term production data analysis results using RTA technique and history matching approach for this field case. A significant decrease in xf from 750 ft. to 430 ft. and in kf from 883 md to 212 md is observed during production. It is consistent with the observation in Fig. 13a that the kf drops from 883 md to 394 md after

5. Discussion In this paper, we proposed a 1D two-phase model to analyze waterphase flowback data of MFHWs in shale reservoirs and a two-phase diagnostic plot to identify flow regimes during flowback period. Pressure-dependent fracture permeability and porosity are accounted for in this model. Compared to other flowback models, the proposed 11

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

104

1400 Case 4-(a)

Case 4-(b) 1200 BDF

0.26 day

103

1

2 2

10

TR

RNP (psia×day/STB)

DRNP (psia×day/STB)

Transient Flow (TF)

1

1

Numerical simulation Straight line 101 103

1000

BDF 0.26 day

800 600 400 Numerical simulation Straight line

200 0

104

105

106

105

0

Superposition Pseudotime (psia×day/cp)

2x105

3x105

4x105

Pseudotime (psia×day/cp)

104

1400 Case 5-(b)

Case 5-(a) 1200 BDF

TF

0.43 day

103

1

2 2

10

RNP (psia×day/STB)

DRNP (psia×day/STB)

Transient Flow (TF)

1 1

Numerical simulation Straight line 101 103

1000

BDF 0.43 day

800 600 400 Numerical simulation Straight line

200 0

104

105

106

105

0

Superposition Pseudotime (psia×day/cp)

2x105

3x105

4x105

Pseudotime (psia×day/cp)

104

1400 Case 6-(a)

Case 6-(b) 1200 BDF

0.24 day

103

1

2 2

10

TF

RNP (psia×day/STB)

DRNP (psia×day/STB)

Transient Flow (TF)

1

1

Numerical simulation Straight line 101 103

1000

BDF 0.24 day

800 600 400 Numerical simulation Straight line

200 0

104

105

106

0

Superposition Pseudotime (psia×day/cp)

105

2x105

3x105

4x105

Pseudotime (psia×day/cp)

Fig. 7. Two-phase (a) diagnostic plot and (b) specialty plot for Case 4–6 using constant rate flowback model.

is assumed that pseudo average pressure can be approximated by average pseudopressure. Therefore, the proposed method is not applicable to the systems with significant nonlinearity in two-phase flow and pressure dependent permeability. As the nonlinearity increases, the approximation is expected to be worse. Also, the variable rate method of calculating average pressure is very sensitive to the data quality, which needs careful attention to smoothen the production data. More work needs to improve this method to make it robust. The flowback model in this paper is developed for early gas production period, which assumes water and gas sourced from the matrix are negligible. This could be invalid for reservoirs with relatively high matrix permeability or secondary fractures. When the fluid influx from the matrix cannot be overlooked, our model could be inapplicable. The

models consider water-phase relative permeability into pseudopressure to linearize the governing equation. Furthermore, this paper first developed flowback model under constant BHP condition although it is only valid for cases with low extent of nonlinearity, whereas other researchers derived flowback models under constant rate condition. With the flowback model, we proposed a comprehensive analysis workflow to iteratively calculate Vfi, rather than considering it as a known parameter. The validation compares the analysis result with the twophase flowing material balance models proposed by Xu et al. (2016) and Behmanesh et al. (2018b). The results show the superiority and practical applicability of the proposed model to estimate fracture properties. In the derivation of the new average pressure calculation method, it 12

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

104

1400 Case 3-(b)

Case 3-(a) 1200 BDF

TF

RNP (psia×day/STB)

DRNP (psia×day/STB)

Transient Flow (TF) 0.14 day

103

1 2

1

2

10

1

Numerical simulation Straight line 101 103

1000

BDF 0.26 day

800 600 400 200

Numerical simulation Straight line

0 104

105

106

105

0

Superposition Pseudotime (psia×day/cp)

2x105

3x105

4x105

Pseudotime (psia×day/cp)

Fig. 8. Two-phase (a) diagnostic plot and (b) specialty plot for Case 3 using variable rate flowback model.

modeling of 2D two-phase flow during the late flowback period needs further investigations.

input. Workflow III provides an alternative analysis method to calculate fracture properties under constant BHP condition. 5. A field example in Horn River Shale is analyzed using the proposed flowback model to demonstrate the practical application of the mathematical developments presented in this study. The field application shows that the fracture half-length and permeability are reduced by twice and four times, respectively, after three months of production.

6. Conclusions This work presents a semi-analytical two-phase model to analyze early-time flowback data from MFHWs and evaluate fracture properties, such as fracture half-length, initial fracture permeability, and initial pore volume of fracture, under constant BHP, constant rate, and variable rate/BHP condition. The following important conclusions can be drawn:

Nomenclatures bp bq Bg Bw Bwf Bwi cf ct cw DDA DRNP Gi Gr Gp h HF kf kfi krw mp

1. The two-phase semi-analytical model is based on water diffusivity equation and considers pressure-dependent fracture permeability and porosity. By introducing new definitions of pseudopressure and pseudotime, the proposed two-phase diagnostic plot is able to identify flow regimes more accurately during flowback period. 2. A diffusivity equation approach is developed to calculate average pressure under constant BHP, constant rate, and variable rate/BHP condition. The method releases the demand for gas-phsae production data as input and is not affected by the extent of nonlinearity. 3. Comparing against input fracture parameters in numerical simulation, the developed constant BHP flowback model has more error and is only applicable to cases with slight nonlinearity. The constant/variable-rate model is more accurate with relative error < 10%. 4. The proposed three analysis workflows can closely predict fracture properties. The iterative procedure is robust and practical for the field applications. Workflow I is valid for cases with gas and waterphase flowback data both available, whereas Workflow II and III are more applicable and only need water-phase production data as

The intercept of ln(RNP) vs. tp plot. The intercept of RNP vs. tp plot. Gas formation factor, ft3/scf. Water formation factor, bbl/STB Bottomhole water formation factor, bbl/STB Initial water formation factor, bbl/STB Fracture compressibility, psi−1 Total compressibility, psi−1 Water compressibility, psi−1 Dynamic drainage area Derivative of rate normalized pseudopressure, psia × d/STB Initial gas storage in the fracture, scf Remaining gas in fracture, scf Cumulative gas production, scf Fracture height, ft Hydraulic fracture Fracture permeability, md Initial fracture permeability, md Water relative permeability The slope of ln(RNP) vs. tp plot

Fig. 9. Pressure profile in 1D fracture under (a) constant BHP condition and (b) constant rate condition. 13

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

Table 3 Set and calculated values of fracture properties and relative error using the new flowback model for the six cases and the comparison with the results analyzed by the two-phase flowing material balance method of Xu et al. (2016) and Behmanesh et al. (2018b). Case #

Inputs

1 2 3 4 5 6 3 (variable rate)

This work

Xu et al. (2016)

Behmanesh et al. (2018b)

xf (ft)

kfi (md)

ppi-ppwf (psia)

xf (ft)

kfi (md)

ppi-ppwf (psia)

xf (ft)

kfi (md)

xf (ft)

kfi (md)

1000 1000 1000 1000 1000 1000 1000

800 800 800 800 800 800 800

6542 6216 2025 -

941 845 702 1008 1010 1052 1019

752 564 515 796 792 739 904

6499 6194 2052 -

981 1035 1019 1040 1068 933 1019

3470 3263 173 846 655 681 173

999 999 1003 1000 1004 1009 1003

646 657 938 680 661 821 938

Table 4 Comparison of the proposed model with two RTA models.

Pseudopressure, pp Pseudotime, tp

RNP =

b

kfi ∙ f (p) ∙ krw (s w ) μw (p) B w (p) krw (s w ) ∙ f (pf )

dp

∫0t



f

ρgsc kf krg ∫pp ⎡⎢ ρ μ B + b



ref

g g

ρwsc kf krw ⎤ ρref μw B w ⎥ ⎦

dp a



μti Cti t qt ∫0 μ C dτ d qe t t

qg Bgi + qw B w krg

RNP = f

1 kfi

(qg∗, j − qg∗, j − 1 )(t p, N − t p, j − 1 ) c ∗ qg, N

qg∗ = ppi − ppwf qw

Behmanesh et al. (2018b)

krg b ∼ dτ μg Ct

ΣNj = 1

pressure dependent k and ϕ • Considers • Requires only water production data

Comments

d

∫pp

qw

RNP/PNR

c

μwi B wi kfi

2p ∫pp μ Z dp b g

1 ds w ⎞ ⎛ s w μw ⎜c w + cf + g (pf ) s w dpf ⎟ ⎝ ⎠ (qw, j − qw, j − 1 )(t p, N − t p, j − 1 ) ΣNj = 1 qw, N

Flowrate

b

Xu et al. (2016)

∫0t

Superposition pseudotime, tsp

a

This work

ρ

ρ gsc ⎞ qe = ⎛ qg + ⎛ wsc ⎞ qw ρ ρ ⎝ ref ⎠ ⎝ ref ⎠ ⎜

ppi − ppwf qg∗

PNR =







qe ppi − ppwf

constant k and ϕ • Assumes • Assumes constant k and ϕ • Requires both water and gas production data • Requires both water and gas production data f

f

f

Exteneded from original formulation to include pressure-dependant permeability. Gp Bg Wp ∼ 1 ∂V fi Ct = ⎛1 − G ⎞ B Sgi Cg + ⎛1 − W ⎞ Swi Cw + V ∂p . fi ⎠ fi ⎠ gi fi f ⎝ ⎝ Exteneded from original formulation to include variable-rate conditions. 1 μt = ρ k krg ρ k krw . f gsc f + wsc ρref μg Bg ρref μw B w

mq MFHWs p pb pf pi pp ppf ppf ppi ppwf

Table 5 Validation results of Workflow I. Iteration #

Input Vfi, ft3

Output xf, ft

Output kfi, md

Output Vfi, ft3

Error of Vfi, %

1 2 3 4 5

1500 1635 1774 1906 2025

817 887 953 1012 1063

1850 1149 905 790 728

1635 1774 1906 2025 2125

8.3 7.8 6.9 5.9 4.7

The slope of RNP vs. tp plot Multi-fractured horizontal wells Pressure, psia Base pressure, psia Average pressure in the fracture, psia Initial fracture pressure, psia Pseudopressure, psia Pseudo average pressure in the fracture, psia Average pseudopressure in the fracture, psia Initial pseudopressure, psia Bottomhole pseudopressure, psia

Fig. 10. (a) Calculated average pressure and (b) generated specialty plot during iterations of Workflow I. 14

f

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

Table 6 Validation results of Workflow III. Itera-tion #

Input Vfi, ft3

Input ζ, 10−4 md/ft2

Output xf, ft

Output kfi, md

Output Vfi, ft3

Output ζ, 10−4 md/ft2

Error of Vfi, %

Error of ζ, %

1 6 18 49 69

1000 1434 1795 2135 2270

4.00 4.85 5.50 6.32 6.74

563 742 906 1071 1137

135 272 454 727 875

1127 1484 1813 2143 2276

4.27 4.94 5.53 6.34 6.76

11.2 3.4 1.0 0.3 0.2

6.4 1.8 0.6 0.3 0.3

Fig. 11. (a) Calculated average pressure and (b) generated specialty plot during iterations of Workflow III. 7000

18

4000

6000

1400

(a)

2

4000 3000

2000 Water flowrate

2000

8 6

100

200

300

400

500

3000

2000

1000

800

Gas flowrate

600

4

0

0 0

4000

Water flowrate

2

Pressure 0

10

1000

1000

1000

12

Pressure

Flowback Period

4

3000 Gas flowrate

1200

14 Water Flowrate (bbl/d)

6

5000

Gas Flowrate (MMSCF/d)

8

5000 Bottomhole Pressure (psia)

10

6000 Water Flowrate (bbl/d)

Gas Flowrate (MMSCF/d)

12

(b)

16

600

0 20

400

200 40

60

80

100

Time (day)

Time (hour)

Fig. 12. Production history of gas flowrate, water flowrate, and bottomhole pressure for the field case during (a) flowback and (b) long-term production.

psc pwf qDw q qw RNP RTA sg sgi sw sw swf swi Vfi Vwr wf Wfi Wp

x xf t tp tsp Tsc ϕf ϕfi ρgsc ρref ρw ρwsc μw μg γ ζ

Pressure at standard condition, psia Bottom-hole pressure, psia Dimensionless water flowrate The water-phase downhole flowrate, bbl/d The water-phase surface flowrate from the whole fracture, STB/d Rate normalized pseudopressure, psia × d/STB Rate transient analysis Gas saturation Initial gas saturation Water saturation Average water saturation Average water saturation in fracture at the end of flowback Initial water saturation Initial pore volume of the fracture, ft3 Remaining water volume in the fracture, ft3 Hydraulic fracture width, ft Initial water storage in fracture, ft3 Cumulative water production, bbl 15

Location in the fracture, ft Fracture half-length, ft Time, day Pseudotime, psia × day/cp Superposition pseudotime, psia × day/cp Temperature at standard condition, °F Fracture porosity Initial fracture porosity Gas density at standard condition, lb/ft3 Reference density, lb/ft3 Water density, lb/ft3 Water density at standard condition, lb/ft3 Water viscosity, cp Gas viscosity, cp Permeability modulus, psi−1 Fracture coefficient, md/ft2

Bottomhole Pressure (psia)

14

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

Table 7 Real field data from an MFHW in horn river shale.

Table 8 Comparison of flowback analysis and long-term production data analysis results.

Parameter

Value

Initial pressure, pi(psi) Initial reservoir pressure, pir(psi) Reservoir temperature, T (°F) Specific gravity, SGg Water compressibility, cw(psi−1) Initial fracture porosity, ϕfi (%) Matrix porosity, ϕm (%) Matrix permeability, km (md) Fracture compressibility, cf (psi−1) Matrix compressibility, cm (psi−1) Reservoir thickness, h (ft) Horizontal wellbore length, Lw (ft) Number of fractures Fracture spacing, Lf (ft) Fracture width, wf (ft) Initial water formation volume factor, Bwi(bbl/STB) Initial fracture water saturation, swi (%) Matrix water saturation, swm (%) Initial water viscosity, μw(cp) Initial water density, ρw(lb/ft3) Fracture permeability modulus, γ(psi−1)

5900 4900 278 0.66 3.33 × 10−6 60 5 1 × 10−4 8.14 × 10−5 6.48 × 10−6 394 4593 15 306 0.08 1.03 85 10 0.33 62 5 × 10−4

Fracture properties

Flowback analysis

xf, ft kfi, md Average kf, md

750 883 –

Long-term production data analysis RTAa

History matching

103–462 – –

430 – 212

Matrix permeability for this MFHW is in a range of 5 × 10−6–1 × 10−4

a

md.

D p q

Dimensionless Pseudo or constant BHP model Pseudo or constant rate mode

Declaration of Competing Interest The authors declare no conflict of interest. Acknowledgments The authors acknowledge the financial support from the Penn State Institutes of Energy and the Environment and College of Earth and Mineral Sciences for the PSIEE Seed Grant and Wilson Initiation Grant. The first author gratefully acknowledges the support of SPE Foundation through the Nico van Wingen Memorial Graduate Fellowship in petroleum engineering and SME (Society for Mining, Metallurgy & Exploration) through the WAAIME Scholarship. This research was enabled with the use of the software packages provided by Computer Modelling Group, Inc. and IHS Markit.

Subscripts i j fi f g sc w

Initial The jth flow rate interval Initial fracture Fracture Gas Standard condition Water 1000

6000

100

800 Pressure 5000

600 4500

Permeability 400 4000

200

Average Pressure in Fracture (psia)

(b) 5500

Average Fracture Water Saturation

Average Fracture Permeability (md)

(a) 8x10

6x10-1

4x10-1

2x10-1

3500

0

5

10

15 Time (day)

20

-1

0 4000

25

102

4500 5000 5500 Average Pressure in Fracture (psia)

6000

14 (d)

(c)

+ unit slope

RNP (psia×day/STB)

DRNP (psia×day/STB)

12

+ 1/2 slope

101

DRNP Straight line 100 104

10 8 6 4 RNP Straight line

2 0

105 Superposition Pseudotime (psia×day/cp)

106

0

4x104 8x104 1x105 2x105 Superposition Psuedotime (psia×day/cp)

2x105

Fig. 13. Plots during flowback analysis after the convergence of Vfi: (a) Average pressure and permeability in the fracture during flowback period. (b) The path of average fracture water saturation vs. average fracture pressure. (c) Two-phase diagnostic plot of water where the pseudo-variables are calculated based on the average pressure from plot (a). (d) Specialty plot of two-phase water flowback model. 16

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi

102

35 (b) 30 Gas RNP (psia2×day/Mscf)

+ 1/2 slope

101

DRNP Straight line

20 15 10 RNP Straight line

5

100

0 1

10 100 Gas Superposition Pseudotime (day)

1000

0

20

8 6

(d) Bottomehole Pressure (psia)

10

Field data History match

Flowback Period

12

16

2000

16 14

4 8 12 Square Root of Gas Superposition Psuedotime (day2)

(c)

18 Gas Flowrate (MMSCF/d)

25

4

1500

1000

Flowback Period

2

Gas DRNP (psia ×day/Mscf/ln(day))

(a)

Field data History match

500

2 0 20

40

60 Time (day)

80

0 20

100

40

60 Time (day)

80

100

Fig. 14. Plots during long-term production data analysis: (a) Single-phase diagnostic plot of gas where the pseudo-variables are calculated based on the average pressure in the distance of inverstigation in the matrix. (b) Specialty plot of single-phase gas formation transient linear flow. (c) History match of gas flowrate and (d) history match of bottomhole pressure during long-term production using the Horizontal-Multifractured-SRV model in IHS harmony software. Clarkson, C.R., Williams-kovacs, J., 2013. Modeling two-phase flowback of multifractured horizontal wells completed in shale. SPE J. 18, 795–812. Clarkson, C.R., Yuan, B., Zhang, Z., Tabasinejad, F., Behmanesh, H., Hamdi, H., Anderson, D., Thompson, J., Lougheed, D., 2019. Anomalous diffusion or classical diffusion in an anomalous reservoir? evaluation of the impact of multi-phase flow on reservoir signatures in unconventional reservoirs. In: SPE/AAPG/SEG Unconventional Resources Technology Conference. Unconventional Resources Technology Conference. Crafton, J.W., 1998. Well evaluation using early time post-stimulation flowback data. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Crafton, J.W., 2008. Modeling flowback behavior or flowback equals" slowback". In: SPE Shale Gas Production Conference. Society of Petroleum Engineers. Dake, L.P., 1983. Fundamentals of Reservoir Engineering. Elsevier. Ezulike, O.D., Dehghanpour, H., 2014. Modelling flowback as a transient two-phase depletion process. J. Nat. Gas Sci. Eng. 19, 258–278. Ilk, D., Currie, S.M., Symmons, D., Rushing, J.A., Broussard, N.J., Blasingame, T.A., 2010. A comprehensive workflow for early analysis and interpretation of flowback data from wells in tight gas/shale reservoir systems. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Jia, P., Cheng, L., Clarkson, C.R., Williams-Kovacs, J.D., 2017. Flow behavior analysis of two-phase (gas/water) flowback and early-time production from hydraulically-fractured shale gas wells using a hybrid numerical/analytical model. Int. J. Coal Geol. 182, 14–31. Jia, P., Cheng, L., Clarkson, C.R., Huang, S., Wu, Y., Williams-kovacs, J.D., 2018. A novel method for interpreting water data during flowback and early-time production of multi-fractured horizontal wells in shale reservoirs. Int. J. Coal Geol. 200, 186–198. Jones, M.R., Blasingame, T.A., 2019. A direct method for short-term forecasting of multiphase production rates using flowback data. In: SPE/AAPG/SEG Unconventional Resources Technology Conference. Unconventional Resources Technology Conference. Jr Pedrosa, O.A., 1986. Pressure transient response in stress-sensitive formations. In: SPE California Regional Meeting. Society of Petroleum Engineers. King, G., 1993. Material-balance techniques for coal-seam and devonian shale gas reservoirs with limited water influx. SPE Reserv. Eng. 8, 67–72. Miller, F.G., 1962. Theory of unsteady-state influx of water in linear reservoirs. J. Inst. Pet. 48, 365–379. Moghadam, S., Jeje, O., Mattar, L., 2011. Advanced gas material balance in simplified format. J. Can. Pet. Technol. 50, 90–98.

References Abbasi, M.A., Ezulike, D.O., Dehghanpour, H., Hawkes, R.V., 2014. A comparative study of flowback rate and pressure transient behavior in multifractured horizontal wells completed in tight gas and oil reservoirs. J. Nat. Gas Sci. Eng. 17, 82–93. Adefidipe, O., Dehghanpour, H., Virues, C., 2014. Immediate gas production from shale gas wells: a two-phase flowback model. In: SPE Unconventional Resources Conference. Society of Petroleum Engineers. Alkouh, A., Mcketta, S., Wattenbarger, R.A., 2014. Estimation of effective-fracture volume using water-flowback and production data for shale-gas wells. J. Can. Pet. Technol. 53, 290–303. Altwaijri, M., Xia, Z., Yu, W., Qu, L., Hu, Y., Xu, Y., Sepehrnoori, K., 2018. Numerical study of complex fracture geometry effect on two-phase performance of shale-gas wells using the fast EDFM method. J. Pet. Sci. Eng. 164, 603–622. Anderson, D.M., Turco, F., Virues, C.J.J., Chin, A., 2013. Application of rate transient analysis workflow in unconventional reservoirs: horn river shale gas case study. In: SPE Unconventional Resources Conference and Exhibition-Asia Pacific. Society of Petroleum Engineers. Behmanesh, H., Hamdi, H., Clarkson, C.R., Thompson, J.M., Anderson, D.M., 2018a. Analytical modeling of linear flow in single-phase tight oil and tight gas reservoirs. J. Pet. Sci. Eng. 171, 1084–1098. Behmanesh, H., Mattar, L., Thompson, J.M., Anderson, D.M., Nakaska, D.W., Clarkson, C.R., 2018b. treatment of rate-transient analysis during boundary-dominated flow. SPE J. 23, 1145–1165. Chen, Z., Liao, X., Yu, W., Zhao, X., 2018. Transient flow analysis in flowback period for shale reservoirs with complex fracture networks. J. Pet. Sci. Eng. 170, 721–737. Clarkson, C., 2012. Modeling 2-phase flowback of multi-fractured horizontal wells completed in shale. In: Paper SPE 162593 presented at SPE Canadian Unconventional Resources Conference, https://doi.org/10.2118/162593-MS. Calgary, Alberta, Canada, 30 October-1 November. dx. Clarkson, C.R., 2013. Production data analysis of unconventional gas wells: review of theory and best practices. Int. J. Coal Geol. 109-110, 101–146. Clarkson, C.R., Qanbari, F., 2016a. History matching and forecasting tight gas condensate and oil wells by use of an approximate semianalytical model derived from the dynamic-drainage-area concept. SPE Reserv. Eval. Eng. 19, 540–552. Clarkson, C.R., Qanbari, F., 2016b. A semi-analytical method for forecasting wells completed in low permeability, undersaturated CBM reservoirs. J. Nat. Gas Sci. Eng. 30, 19–27.

17

International Journal of Coal Geology 217 (2020) 103315

F. Zhang and H. Emami-Meybodi Ozkan, E., Raghavan, R.S., Apaydin, O.G., 2010. Modeling of fluid transfer from shale matrix to fracture network. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Qanbari, F., Clarkson, C.R., 2016. Rate-transient analysis of liquid-rich tight/shale reservoirs using the dynamic drainage area concept: examples from north American reservoirs. J. Nat. Gas Sci. Eng. 35, 224–236. Song, B., Ehlig-Economides, C.A., 2011. Rate-normalized pressure analysis for determination of shale gas well performance. In: North American Unconventional Gas Conference and Exhibition. Society of Petroleum Engineers. Uzun, I., Kurtoglu, B., Kazemi, H., 2016. Multiphase rate-transient analysis in unconventional reservoirs: theory and application. SPE Reserv. Eval. Eng. 19, 553–566. Wattenbarger, R.A., El-Banbi, A.H., Villegas, M.E., Maggard, J.B., 1998. Production analysis of linear flow into fractured tight gas wells. In: SPE Rocky Mountain Regional/Low-permeability Reservoirs Symposium. Society of Petroleum Engineers. Williams-Kovacs, J.D., Clarkson, C.R., 2016. A modified approach for modeling two-phase flowback from multi-fractured horizontal shale gas wells. J. Nat. Gas Sci. Eng. 30,

127–147. Xu, Y., Adefidipe, O., Dehghanpour, H., 2016. A flowing material balance equation for two-phase flowback analysis. J. Pet. Sci. Eng. 142, 170–185. Yang, Y., 2017. Mathematical Development for Flowback Rate Transient Analysis. Pennsylvania State University. Yang, R., Huang, Z., Li, G., Yu, W., Sepehrnoori, K., Lashgari, H.R., Tian, S., Song, X., Sheng, M., 2017. A semianalytical approach to model two-phase flowback of shalegas wells with complex-fracture-network geometries. SPE J. 22 (6), 1–808. Zhang, F., Emami-Meybodi, H., 2018. Flowback fracture closure of multifractured horizontal wells in shale gas reservoirs. In: SPE/AAPG Eastern Regional Meeting. Society of Petroleum Engineers. Zhang, F., Emami-Meybodi, H., 2019. Flowback fracture closure of multi-fractured horizontal wells in shale gas reservoirs. J. Pet. Sci. Eng (In press). Zhang, J., Kamenov, A., Zhu, D., Hill, A., 2013. Laboratory measurement of hydraulic fracture conductivities in the Barnett shale. In: IPTC 2013: International Petroleum Technology Conference.

18