Flowback fracture closure of multi-fractured horizontal wells in shale gas reservoirs

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: http://www.elsevier.com/locate/petrol

Flowback fracture closure of multi-fractured horizontal wells in shale gas reservoirs Fengyuan Zhang, Hamid Emami-Meybodi * Department of Energy and Mineral Engineering and EMS Energy Institute, The Pennsylvania State University, University Park, PA, 16802, USA

A R T I C L E I N F O

A B S T R A C T

Keywords: Flowback Multiphase flow Production data analysis Fracture closure Shale reservoirs Multi-fractured horizontal wells

In multi-fractured horizontal wells (MFHW), fracture properties such as permeability and fracture half-length significantly deteriorate during production, which negatively affects gas production from shale reservoirs. Therefore, it is crucial to evaluate the temporal changes in fracture properties using production data. This paper presents two multiphase flowback models for gas and water phase under boundary dominated flow (BDF) condition by considering gas influx from matrix. In addition, a workflow is proposed to quantitatively evaluate hydraulic fracture closure and fracture properties using both water-phase and gas-phase flowback data. A ma­ terial balance equation (MBE) is derived to obtain average pressure in the fracture by calculating gas influx from matrix under variable bottomhole pressure (BHP) condition based on Duhamel’s principle. The flowback models, average pressure calculation method, and workflow are validated against numerical simulations and applied to an MFHW drilled in Marcellus Shale. The results show that the developed models are capable of predicting fracture properties, and MBE could accurately calculate average pressure in the fracture. Furthermore, the proposed workflow provides quantitative insights into the performance of stimulation jobs by accurately pre­ dicting fracture half-length, fracture permeability, and permeability modulus using flowback data. By reconciling flowback analysis with long-term production data analysis, a significant decrease in fracture properties due to fracture closure is observed in the Marcellus MFHW.

1. Introduction Multi-fractured horizontal wells (MFHWs) are widely used to acquire economic production from unconventional resources such as shale res­ ervoirs. Evaluating hydraulic fracturing performance is very critical for MFHWs. One way to assess stimulation performance is rate transient analysis (RTA), which provides important matrix and fracture proper­ ties. In unconventional reservoirs such as shales, RTA usually focuses on analyzing long-term production data obtained over several months to even several years where various flow regimes can be identified (Clarkson, 2013). Apart from long-term RTA, fracture and reservoir properties could also be estimated from short-term data obtained during the flowback period. Flowback is the process of allowing injected fluids to flow from well to the surface after hydraulic fracturing. Although fracturing fluid leaks off to the formation, the flow in the matrix is still dominated by single-phase gas flow due to the low water relative permeability (Clarkson et al., 2017). Thus, two-phase flow mainly oc­ curs in hydraulic fractures. Compared to long-term RTA, flowback RTA provides insights on fracture information at the early times and hence

could be used to evaluate fracture closure during production. Over the past decade, numerous studies have qualitatively analyzed flowback data. In 2008, Crafton studied the effects of flowback rate, fracture complexity, and shut-ins on well performance and recovery (Crafton, 2008). Ilk et al. (2010) presented a diagnostic log-log plot showing that gas water ratio (GWR) vs. cumulative gas production (Gp) plot has a half slope, which was later used by others (Ghanbari et al., 2013; Abbasi et al., 2014; Zhang and Ehlig-Economides, 2014). How­ ever, there is no mathematical proof for this empirical diagnostic plot. Clarkson and colleagues used GWR vs. Gp to identify multiphase fracture depletion flow regime and the linear flow regime (Williams-Kovacs and Clarkson, 2016; Clarkson and Williams-Kovacs, 2013b). Later, Clarkson and Williams-Kovacs, (2013)a divided the flowback process into before-breakthrough (BBT) and after-breakthrough (ABT) according to the time when gas enters the fracture from matrix. During BBT, only single-phase water flows in fracture. When it comes to ABT, formation gas starts entering fracture and the two-phase flow occurs. Adefidipe et al. (2014) extended BBT period to two-phase flow with flowing fluids in-situ water and gas in the hydraulic fracture. They divided flowback

* Corresponding author. E-mail address: [email protected] (H. Emami-Meybodi). https://doi.org/10.1016/j.petrol.2019.106711 Received 5 June 2019; Received in revised form 21 October 2019; Accepted 17 November 2019 Available online 21 November 2019 0920-4105/© 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Fengyuan https://doi.org/10.1016/j.petrol.2019.106711

Zhang,

Hamid

Emami-Meybodi,

Journal

of

Petroleum

Science

and

Engineering,

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

period to early gas production (EGP) and late gas production (LGP) using GWR vs. time plot. Numerous studies have addressed the modeling of single-phase flowback (Abbasi et al., 2012; Wattenbarger and Alkouh, 2013; Kurto­ glu et al., 2015). However, single-phase model is only valid for the very early period because the flowback is multiphase flow in most of the time. To model two-phase water and gas flow in the fracture, Ezulike et al. (2013) developed a two-phase flowback model for MFHWs which con­ siders gas influx. To solve coupled equations, they proposed the Dy­ namic Relative Permeability (DRP) concept using Mellin transform. In an extended study, they improved their approach by simplifying DRP with an empirical function and used Laplace transforms to get the semi-analytical solution (Ezulike and Dehghanpour, 2014). Yang (2017) developed two-phase flowback models for both EGP and LGP, where the LGP model assumes the equal average pressure dropping rate inside fracture and matrix, which is far from the real case. Clarkson and Williams-Kovacs, (2013)a extended RTA method designed for multi-phase flow in CBM to analyze early two-phase flowback data in cylindrical-source MFHWs. They coupled pseudo-steady state boundary dominated equation with MBE to develop a two-phase model for frac­ ture flow (Williams-Kovacs et al., 2015). In order to model two-phase formation linear flow towards fracture, Clarkson and Qanbari (2016a) developed a semi-analytical model for the multi-phase flow of MFHWs in tight reservoirs using Dynamic Drainage Area (DDA) concept adopted from welltest analysis (Shahamat et al., 2015). This concept was first applied to coal-bed methane (CBM) wells (Clarkson and Qanbari, 2016b) and extended to analyze multi-phase flowback period for tight oil reservoirs (Clarkson et al., 2016). Later, the DDA flowback model was modified by including fracture propagation before flowback and frac­ turing liquid leak-off during stimulation (Clarkson et al., 2017). Recent studies have included complex fracture networks and fracture closure in flowback modeling. Jia et al. (2017) modeled flow through complex fracture networks in shale gas reservoirs for two-phase flow­ back using a hybrid numerical and analytical model. Yang et al. (2017) developed a semi-analytical LGP model for MFHWs in shale gas reser­ voir, which considered both gas macroscopic transport mechanism and complex fracture network geometries. However, this model can only be applied to forward simulation not fracture properties inversion. Fu et al., (2019) applied decline curve analysis to determine initial fracture vol­ ume and used fracture compressibility to evaluate fracture volume loss during flowback. Hossain et al. (2019) developed a post-flowback water production model to estimate the average fracture volume and observed a significant decrease in fracture volume during flowback by analyzing the field case. Jones and Blasingame (2019) proposed hyperbolic and modified-hyperbolic models to empirically history match and forecast total liquid productivity index and bottomhole pressure during flowback. Flowing material balance method (FMB) have been widely used in flowback modeling (Abbasi et al., 2014; Clarkson and Williams-Kovacs, 2013a, 2013b; Jia et al., 2018). The FMB method could be traced back to the work by Palacio and Blasingame (1993), in which they developed a single-phase gas BDF model in a radial system and proposed material balance pseudotime to model variable rate production based on a con­ stant rate (pseudo-steady state) solution. Agarwal et al., (1999) rear­ ranged the BDF solution to the pressure normalized rate form. Clarkson et al. (2008) extended the single-phase FMB method to two-phase gas flow in CBM reservoirs by modifying the definition of pseudopressure and total compressibility in pseudotime. However, their approach has many assumptions; for example, the production of CBM well is domi­ nated by gas and total compressibility is approximated by gas compressibility to apply the Blasingme’s material balance pseudotime. Later, Clarkson, (2013) generated a two-phase FMB model and a workflow for horizontal CBM wells. Sun et al. (2018) considered pressure-dependent permeability into the two-phase FMB model. Recently, Shahamat and Clarkson (2018) developed a general FMB model for multiphase (oil, gas, and water) flow by redefining multiphase

Fig. 1. Schematic of the conceptual model for flowback depicting water flow through hydraulic fracture and gas flow through fracture and matrix.

pseudopressures. However, FMB still approximates total compressibility with gas compressibility. Behmanesh et al. (2018b) relaxed the total compressibility assumption and derived a new generalized FMB equa­ tion based on the governing equation of total flowrate. This paper presents semi-analytical models for gas and water flow through fractures by incorporating pressure-dependent permeability, gas influx from matrix, and an MBE to obtain average pressure in frac­ tures. A workflow is proposed to calculate fracture properties and evaluate fracture closure by combining water flowback model and gas flowback model. The flowback models, MBE, and the workflow are validated against numerical simulations and applied to an MFHW drilled in Marcellus Shale in Pennsylvania, USA. Accordingly, the contributions of this work are in: (1) development of a two-phase gas FMB equation for the fracture network, in which the assumption of approximating total compressibility with gas compressibility is relaxed, (2) development of an MBE by considering the variable inflow pressure in the fracture, and (3) proposing a workflow for calculation of fracture half-length, fracture permeability, and permeability modulus, as well as evaluation of frac­ ture closure using flowback data. 2. Theory and methods This section describes two-phase semi-analytical models to analyze water-phase and gas-phase flowback data. First, the conceptual model and flow regime analysis are introduced. Then, mathematical formula­ tions are presented. Last, a workflow is proposed that reconciles waterphase and gas-phase flowback data to evaluate hydraulic fracture properties and its closure. 2.1. Conceptual model Recently, we developed a 1-D two-phase model for the flow of gas and water in fracture toward a wellbore to analyze early-time flowback data (Zhang and Emami-Meybodi, 2020). We extend our previous work by considering two-phase water and gas flow in the fracture and single-phase gas influx from the matrix. The conceptual model is shown in Fig. 1 and is comparable to that assumed by Jia et al. (2018). The main assumptions are as follows: � Gas sources from the initial gas in fracture and matrix contribution due to the influx. Water only sources from the fracture (Clarkson and Williams-Kovacs, 2013a) and water production from matrix is negligible due to the low mobility. Water saturation in the matrix is

2

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 2. (a) Different flow regimes during flowback and (b) two-phase diagnostic plot of water phase (Zhang and Emami-Meybodi, 2020).

� �

�

�

�

�

�

considered to be constant, which is a reasonable assumption during the flowback period (Williams-Kovacs and Clarkson, 2016). Water and gas flow in fractures is assumed to follow Darcy’s law. Fracture is an elastic porous medium and is assumed more compressible than matrix. Fracture permeability and porosity are pressure-dependent variables, whereas matrix permeability and porosity are constants. All the hydraulic fractures have the same attributes, i.e. length, width, and permeability. While linear, elliptical, and radial flow may be observed for different fracture geometries, which are controlled by xf/h (Al-Kobaisi et al., 2006), we only considered linear flow in the fracture. Fracture width and half-length are assumed to be constant during flowback period. However, the effect of fracture closure on produc­ tion is considered through variation of fracture permeability quan­ tified by permeability modulus. Water formation volume factor, water viscosity, gas compressibility, and gas viscosity are pressure-dependent variables. Water compressibility is constant. The variables used in the fracture model are evaluated at average fracture pressure. The variables used in the matrix model are eval­ uated at average pressure within the distance of investigation in the matrix. The gravity segregated relative permeability curves (Clarkson, 2013; Kanfar et al., 2014) are valid for the fracture.

phase water IALF, and FR3_W two-phase water boundary dominated flow. FR1_W appears in the very early time of flowback when the frac­ ture is saturated with mostly water, i.e., high water saturation. In water diagnostic plot, FR1_W can be detected by a half-slope. Theoretically, the corresponding data for FR1_W could be used to estimate fracture permeability, but the short duration of this regime sometimes restricts the feasibility of using a single-phase model. If there is any gas left inside fractures after the stimulation, FR2_W would be the first flow regime to occur. Although a half-slope trend can be observed in the two-phase water diagnostic plot (Fig. 2b), analyzing of FR2_W is impractical due to few numbers of data. When the effect of pressure drop reaches the fractures boundaries, FR3_W is observed with a unit-slope straight line in the two-phase water diagnostic plot (Fig. 2b). Therefore, water flow in our conceptual model is always 1-D flow in spite of the complex gas flow, which is the coupling of matrix influx with the flow in the fracture. Here we used the two-phase diagnostic plot proposed by Zhang and Emami-Meybodi, (2020), which is based on multiphase pseudopressure and pseudotime definitions. As pointed out by Clarkson et al. (2019), using single-phase diagnostic plots to identify flow regimes in multi­ phase flow systems may lead to erroneous results. The behavior of gas production over time is different as compared with that of water. We generated the gas-phase diagnostic plot in Fig. 3b, where the definitions of pseudo variables are listed in Zhang and Emami-Meybodi (2018) but are evaluated at average pressure in the distance of investigation in the matrix. The three flow regimes are ex­ pected to occur: gas bilinear flow (FR1_G), gas formation IALF flow (FR2_G), and gas formation boundary dominated flow (FR3_G). FR1_G corresponds to FR2_W and occurs at the beginning of production during which both water (in fractures) and gas (in fractures and influx from the matrix) flow toward wellbore. As water is produced from fracture, gas flow becomes dominant and can be identified by one-fourth slope in the

2.2. Flow regime analysis As mentioned above, water is the dominant flowing phase during the early flowback period and only flows inside fracture. As shown in Fig. 2a, three flow regimes of water are expected to occur in the fracture: FR1_W single-phase water infinite acting linear flow (IALF), FR2_W two-

Fig. 3. (a) Different flow regimes during flowback and (b) gas phase diagnostic plot (Zhang and Emami-Meybodi, 2018). 3

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

gas-phase diagnostic plot. However, single-phase gas bilinear flow model is not accurate enough to extract fracture permeability from FR1_G data due to the influence of water flow in the fracture. When the effect of pressure reaches the boundary of fractures, FR2_G will appear with a half-slope (see Fig. 3b), which usually takes a long time in shales due to the extremely low matrix permeability. As pressure further propagates in the matrix and reaches the symmetric line between two fractures, FR3_G appears and will show a unit-slope straight line in the diagnostic plot. Numerous studies have investigated gas flow in the matrix during FR2_G and FR3_G to develop formation linear flow and boundary dominated flow models (Nobakht and Clarkson, 2012a, 2012b, 2012c). In recent years, Williams-Kovacs and Clarkson (2016) focused on gas flow in the fracture and treated FR2_G and FR3_G as two-phase gas fracture depletion with the gas influx from matrix, by applying two-phase FMB method designed for CBM wells (Clarkson et al., 2008, 2012).

kf ¼ kfi f ðpÞ ¼ kfi exp½ � ϕf ¼ ϕfi gðpÞ ¼ ϕfi exp

2.3.1. Two-phase water flowback model As discussed by Zhang and Emami-Meybodi, (2020), constant BHP model has more error than the constant rate model, because the constant BHP case has stronger pressure and saturation gradient in the fracture (i. e., the governing equation is more nonlinear). The stronger nonlinearity causes the assumption of applying pseudotime invalid for constant BHP model. Hence, we used constant rate model to analyze water-phase flowback data during FR3_W:

Bwi μwi ; wf hxf ϕfi

(2)

xf Bwi μwi ; wf hkfi

(3)

mw ¼ 2:807

bw ¼ 147:82

xf ;w ¼ 2:807

1 kfi

Z 0

t

krw ðsw Þkf pf

� � sw μw pf ϕf pf cw þ cf þ s1w

(9)

xf ;w Bwi μwi : bw wf h

(10)

2.3.2. Two-phase gas flowback model Gas flowback model is different from water due to the pressuredependent gas properties and gas influx from matrix. Based on the conceptual model depicted in Fig. 1, the gas-phase material balance equation for an infinitesimal element in the fracture is given by: � qgf ρgf �xþdx

� � � d Δx wf h ​ ϕf Sg ρgf ; qgf ρgf �x þ qgm ρgm �y¼0 ¼ dt

(11)

where x is the location in the fracture, y is the location in the matrix, Δx wf h is the control volume, qgf is the gas-phase downhole flowrate at a specific location x in the fracture, qgm is the gas-phase downhole inflow rate at a specific location y from matrix, sg is the gas saturation in the fracture, ρgf is the gas-phase density in the fracture, and ρgm is the gasphase density in the matrix. We neglected the effect of fracture pressure gradient on gas influx and assumed that gas flows uniformly into the fracture. Hence the gas influx term in eq (11) can be included into accumulation term on the right-hand side to treat as a function of time. Dividing both sides of eq (11) by Δx, the following equation can be obtained: � ∂� ∂ qgf ρgf ¼ w h ​ ϕf Sg ρgf ∂x ∂t f

� ! d τ;

(7)

Bwi μwi ; mw wf hϕfi

kfi;w ¼ 147:82

where, RNPw ¼ (ppi-ppwf)/qw is the rate normalized pseudopressure for water, qw is the water-phase flowrate from the whole fracture at surface condition, Bwi is the water formation volume factor evaluated at initial pressure, μwi is the water viscosity at initial pressure, wf is the fracture width, h is the fracture height, ϕfi is the initial fracture porosity, xf is the fracture half-length, kfi is the initial fracture permeability, and mw and bw are the slope and intercept of RNPw vs. tp,w plot, respectively. The definitions of water-phase pseudopressure pp and water-phase pseudotime tp,w are given by: Z μ Bwi p kf ðpÞkrw ðsw Þ pp ¼ wi dp; (4) kfi pb μw ðpÞBw ðpÞ tp;w ¼

� pÞ ;

Cf ðpi

where qw,j is the water-phase flowrate at the surface condition in the jth flowrate interval, tp,w,j is the water-phase pseudotime at the jth flowrate interval, and tsp,w is the water-phase superposition pseudotime. With the two-phase water flowback model, we can generate RNPw and the derivative of RNPw (DRNPw ¼ dRNPw/dln tsp,w) vs. tsp,w in log-log scale as a diagnostic plot to identify FR3_W regime (Zhang and Emami-Meybodi, 2020). Although the heterogeneity in fracture half-length and permeability will affect the transition from FR2_W to FR3W, it will not impact the identification of FR3_W, which always shows a unit-slope straight line when flow in the fracture reaches the outer boundary. RNPw vs. tsp,w in normal scale can be generated as the water-phase specialty plot to calculate xf and kfi with the slope mw and intercept bw of the straight line part for general production conditions:

Based on the flow regime analysis in Section 2.2, we adopted the two-phase water BDF model developed by Zhang and Emami-Meybodi, (2020) to analyze water-phase flowback data corresponding to FR3_W. For gas flow, we focused on FR2_G and developed a two-phase gas FMB to analyze gas-phase flowback data.

(1)

(6)

p�;

where γ is the permeability modulus. Applying the superposition principle (Dake, 1983), the constant rate flowback model can be extended to the variable rate case by replacing pseudotime in eq (1) with superposition pseudotime: � � qw;j qw;j 1 tp;w;N tp;w;j 1 tsp;w ¼ ΣNj¼1 ; (8) qw;N

2.3. Mathematical models

RNPw ¼ mw tp;w þ bw ;

γðpi

� � qgm ρgm �y¼0

(12)

p M

f , Darcy’s gas flowrate along Applying gas Equation of State ρgf ¼ ZRT

(5)

the fracture qgf ¼ wf h

dsw dpf

2h

Δx μkm ∂∂pym , gm

krg kf ∂pf μgf ∂x ,

and inflow rate from matrix qgm ¼

eq (12) can be expanded with chain rule by introducing

effective gas-phase compressibility in the fracture: � � ∂ krg kf pf ∂pf pf ∂pf ¼ ϕf ceff ; ∂x μgf Z ∂x Z ∂t

where krw is the water-phase relative permeability, sw is the average water saturation in the fracture, pf is the average fracture pressure, pb is the base pressure for integral, cw is the water compressibility, cf is the fracture compressibility. Fracture permeability kf and porosity ϕf are pressure-dependent variables (Pedrosa, 1986):

0 ceff

4

∂sg ¼ sg cgf þ sg cf þ ∂pf

∂ @ ϕf ρgf ∂pf 1

Z 0

t

1 2 km ∂pm ρ dtA; wf μgm ∂y gm

(13)

(14)

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

where μgf is the gas viscosity in the fracture, μgm is the gas viscosity in the matrix, pf is the fracture pressure, pm is the matrix pressure, km is the matrix permeability, Z is the gas deviation factor, T is the temperature, R is the ideal gas constant, M is the molar mass of gas, ceff is the effective gas-phase compressibility, and cgf is the gas compressibility in the fracture. First, we derived gas flowback model under constant rate condition and then extended it to the variable rate condition. Equation (13) is subject to the following initial and boundary conditions: � � �� krg kf pf M ∂pf �� psc M ∂pf �� pf ðx; 0Þ ¼ pi ; 2 wf h ¼ q ; ¼ 0; (15) g μgf ZRT ∂x �x¼0 RTsc ∂x �x¼xf

mg ¼ 500:2

tp;g ¼ μgf ceff

i

Z

t

0

μgf

� � krg sg f pf � � � d τ: pf ceff pf g pf

x¼0

where ηg ¼

¼

�� qg μgfi psc Zi T ∂m pf � � ; ¼ 0; 2hwf kfi Tsc ∂x �x¼xf

(24)

μgfi psc xf Zi T hwf kfi Tsc

(25)

;

where, RNPg ¼ (m(pi) – m(pwf))/qg is the rate normalized pseudopressure for gas, Vfi ¼ xf wf hϕfi is the initial pore volume of half fracture. Equation (23) can be rearranged to pressure normalized rate form: mðpi Þ

qg �¼ m pwf

1 500:2 psc Zi T tp;g qg 1 � �� þ : Vfi bpss mðpi Þ m pwf Tsc ceffi bpss

(26)

Substituting eq (22) into (26), we arrive at the FMB equation for twophase gas flow in the fracture: �! mðpi Þ m pf qg 1 1 �¼ � þ Vfi : (27) bpss Vfi bpss mðpi Þ m pwf mðpi Þ m pwf Similar to the two-phase water model, the superposition principle can be applied to obtain the analytical solution of eq (18) and eq (19) for the variable gas-phase flowrate qg: " # qg;N μgfi psc xf Zi T 1 x x2 þ 2ηg tsp;g mðpÞ ¼ mðpi Þ LT; (28) þ 2hwf kfi Tsc 3 xf 2xf2

(17)

By applying pseudopressure and pseudotime, the linearized gov­ erning equation with initial and boundary conditions are: � � ∂ 2 m pf 1 ∂m pf ; (18) ¼ ∂x ηg ∂tp;g �� � ∂m pf �� m pf ðx; 0Þ ¼ mðpi Þ; ∂x �

psc T Zi ; Tsc Vfi ceffi

bpss ¼ 26340:7

where qg is the gas-phase flowrate from the whole fracture at the stan­ dard condition, and psc and Tsc are pressure and temperature at standard condition, respectively. Considering pressure-dependent fracture permeability and porosity in eqs (6) and (7), we defined gas-phase pseudopressure m(p) and pseudotime tp,g to linearize the governing equation by evaluating pressure-dependent parameters at the average fracture pressure and saturation: Z μgfi Zi p pf kfi f ðpÞ krg ðpÞ mðpÞ ¼ dp; (16) kfi μgf ðpÞZðpÞ pb �

(23)

RNPg ¼ mg tp;g þ bpss ;

tsp;g ¼ ΣNj¼1

LT ¼

(19)

qg;j

qg;j

� 1

tp;g;N

tp;g;j

qg;N

� 1

(29)

;

∞ �X psc T xf μgfi Zi N 1 qg;j qg;j 1 Σ exp Tsc hwf kfi j¼1 ðn π Þ2 n¼1 ! � � � nπx tp;g;j 1 cos ; xf

ðnπÞ2 ηg tp;g;N x2f (30)

where qg,j is the gas-phase flowrate at the surface condition in the jth flowrate interval, tp,g,j is the gas-phase pseudotime at the jth flowrate interval, and tsp,g is the gas-phase superposition pseudotime. Equation (28) extends the two-phase gas flowback model under constant rate condition to variable rate by replacing tp and qg in equa­ tions (21), (23) and (26) with tsp,g and qg,N. Accordingly, using the mðpi Þ mðp Þ normalized cumulative production GpN ¼ Vfi mðp Þ m p f vs. pressure ð wf Þ i q normalized rate PNR ¼ mðpi Þ g;Nmðpwf Þ plot (i.e., gas-phase specialty plot or

k

fi . ðϕf ceff μgf Þi The analytical solution of eq (18) subject to eq (19) is given by Miller (1962): " ∞ X � qg μgfi psc xf Zi T 1 x x2 þ 2 ηg tp;g 1 m pf ¼ mðpi Þ þ 2 exp 2 2hwf kfi Tsc 3 xf 2xf ðn πÞ2 n¼1 ! # � � ðnπÞ2 ηg tp;g nπx cos 2 xf xf

FMB plot), xf and kfi can be obtained from the slope mFMB and intercept bFMB of the straight line, which passes through the later portion of the FMB plot:

(20) Considering the long-term approximation for boundary dominated flow (Wattenbarger et al., 1998), the volumetric average pseudo­ pressure can be obtained by taking an integration over eq (20) with respect to x: Z � 1 xf � Zi T psc qg � tp;g ; m pf ¼ m pf dx � mðpi Þ 500:2 (21) xf 0 Tsc hwf xf ϕf ceff i

xf ;g ¼

bFMB ; mFMB wf hϕfi

kfi;g ¼ 26340:7

bFMB μgfi psc xf ;g Zi T : h wf Tsc

(31) (32)

To determine xf and kfi, the FMB plot PNR vs.GpN needs average pressure in the fracture, whose calculation requires xf itself. Therefore, an iterative procedure is required as discussed in the next section. We adopted the systematic procedure proposed by Shahamat and Clarkson (2018) to conduct the FMB analysis.

where mðpf Þ is the fracture average pseudopressure in the field unit. Assuming that the average pseudopressure can be approximated by pseudo-average-pressure mðpf Þ, eq (21) can be rearranged to relate gasphase pseudotime with pseudopressure drop: !# � " Tsc hwf xf ϕf Ceff i tp;g ¼ mðpi Þ m pf (22) 500:2 Zi T psc qg

2.4. Analysis technique

Substituting x ¼ 0 and p ¼ pwf into eq (20) and taking the long-term approximation for boundary dominated flow (Wattenbarger et al., 1998), we have the RNP form solution in field unit:

An analysis technique is proposed that includes the calculation of average pressure in the fracture and a workflow to calculate xf, kfi, and γ. In this study, fracture closure is quantified by kfi and γ while length and 5

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 4. Schematic of the variation of flowing pressure and intermediate parameters in the calculation.

width are kept fixed. In other words, changes in fracture dimensions during fracture closure are considered in γ. This is reasonable because fracture closure is usually characterized by fracture conductivity, which is the product of fracture width and fracture permeability. Note that the fracture closure is directly dependent on average fracture pressure, and the average pressure is affected by production rates.

solution by Wattenbarger et al. (1998). They set the flowing pressure of inflow model to the average pressure in the fracture as an approxima­ tion, although the inflow model by Wattenbarger is derived for constant BHP condition. In this study, we modified their MBE by using the rigorous solution of gas influx under variable BHP condition. Wattenbarger et al. (1998) derived the solution of gas influx during the formation IALF under constant BHP condition in field unit: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffiffi � � h xf � km ϕm μgm ctm i ; (33) qg;inf ¼ mm ðpi Þ mm pwf 315:4 fcp T t

2.4.1. Calculation of average pressure in the fracture In order to generate the specialty plot, we calculated pseudotime based on the average fracture pressure pf at each time. Williams-Kovacs and Clarkson (2016) modified the MBE by King (1993) to calculate pf and considered gas influx from matrix using the formation linear flow

where qg,inf is the flowrate of gas influx at surface condition, mm ðpÞ ¼

Fig. 5. Workflow to evaluate xf, kfi, and γ using both water-phase and gas-phase flowback models. 6

F. Zhang and H. Emami-Meybodi

2

Rp

p 0 μg Z dp

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

is the real gas pseudopressure, ctm ¼ cgm þ cm is the matrix

Table 1 Summary of input data used in numerical simulation.

total compressibility, cgm is the gas compressibility in the matrix, and cm is the formation compressibility in the matrix. fcp is the correction factor proposed by Ibrahim and Wattenbarger (2006) and it is a function of dimensionless drawdown parameter DD. fcp ¼ 1

0:0852DD

0:0857D2D ;

mm ðpi Þ mm pwf : mm ðpi Þ

(34)

ζðτÞ ¼

� mm ðpi Þ mm pf ðτÞ ; mm ðpi Þ mm ðpb Þ

ΔζðτÞj ¼ ζðτÞjþ

psc TZ ; Tsc pf

��

;

3 � 10

Fracture compressibility, cf (psi Matrix compressibility, cm (psi

5

1 1

300 1 � 10

)

2 � 10

)

1 � 10

Reservoir thickness, h (ft)

300

Horizontal wellbore length, Lw (ft)

10000

Fracture width, wf (ft)

0.02

Initial fracture water saturation, swi (%)

70

Number of fractures Fracture spacing, Lf (ft)

Matrix water saturation, swm (%)

6

1.03

0.3

Initial water density, ρw (lb/ft3)

Fracture half-length, xf (ft)

5

0

Initial water viscosity, μw (cp)

Fracture permeability modulus, γ (psi

5

125 80

Initial water formation volume factor, Bwi (bbl/STB)

(37)

6

50

Matrix permeability, km (md)

62 1

500 )

3 � 10

4

workflow follows nine steps: 1. Guess an initial value for xf based on the field experience and assume γ ¼ 10 3 psi 1, which is an upper bound based on experimental data (Ozkan et al., 2010), 2. Calculate average pressure in the fracture pf using eq (40) and the Newton-Raphson technique, 3. Evaluate tp,w and pp for water-phase model and m(p) for gas-phase flow model at each pf , 4. Create two-phase diagnostic plot of water, tsp,w vs. DRNPw in log-log scale, to identify water-phase BDF regime (FR3_W), as well as a specialty plot, tsp,w vs. RNPw in normal scale, to calculate xf,w and kfi,w from eqs (9) and (10), 5. Generate FMB plot GpN vs. PNR in normal scale and pass a straight line through the later portion to calculate xf,g and kfi,g from eqs (31) and (32). The gas-phase diagnostic plot shown in Fig. 3b can assist to identify gas-phase fracture depletion flow regime (FR2_G and FR3_G) with the half-slope straight line, 6. Update xf,new ¼ (xf,w þ xf,g)/2, 7. Repeat step 2 to 6 until the relative error between xf,new and xf,old becomes less than an acceptable tolerance εxf , 8. Update γ with another guess in the range of 10 4 to 10 3 psi 1, and repeat step 2 to 7 until the relative error between kfi,w and kfi,g be­ comes less than an acceptable tolerance εkf , 9. Report xf, kfi, and γ.

(38)

(41)

pf

0.65 )

Matrix porosity, ϕm (%)

where Gp is the cumulative gas production at the surface condition, Bg is the gas formation volume factor, and Sg is the average gas saturation in the fracture that is calculated using the simplified MBE proposed by King (1993): � Wfi 5:615Wp Bw � ��; (42) sg ¼ 1 sw ¼ 1 2Vfi 1 cf pi pf � Bw ¼ Bwi 1 þ cw pi

1

Fracture initial permeability, kfi (md)

where the variations of ζðτÞ and ΔζðτÞj are treated by stepwise functions as shown in Fig. 4b and c, respectively. With the calculated gas influx, the average pressure in the fracture can be solved iteratively from the modified MBE: � �� 2Vfi sgi 2Vfi 1 cf pi pf sg Gp ¼ þ Gg;inf ; (40) Bgi Bg Bg ¼

160

Initial fracture porosity, ϕfi (%)

(39)

ζðτÞj ;

5000

Water compressibility, cw (psi

(35)

h xf qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffiffi km ð1 swm Þϕm μgm ctm i ½mm ðpi Þ 0:3154T Z t ΔζðτÞj mm ðpb Þ � Σ Nj¼01 qffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffi dτ; 0 t τj fcp;j

Initial pressure, pi (psi) Specific gravity, SGg

In our conceptual model, the flowing pressure of the formation IALF model is the average pressure in the fracture, which is not a constant but changing with time. The variation of flowing pressure can be treated by a stepwise function as shown in Fig. 4a. With the constant matrix water saturation swm, flowrate qg,inf and cumulative production Gg,inf of gas influx at surface condition can be obtained using Duhamel’s principle € (Ozisik, 1993): h xf qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffi qg;inf ¼ km ð1 swm Þϕm μgm ctm i ½mm ðpi Þ 315:4 T ΔζðτÞj mm ðpb Þ �Σ Nj¼01 qffiffiffiffiffiffiffiffiffiffiffiffiffi�ffiffi; (36) fcp;j t τj Gg;inf ¼

Value

Reservoir temperature, T (℉)

�

DD ¼

Parameter

(43)

3. Results and discussion

where Wfi is the initial water storage in the fracture at the surface con­ dition, and Wp is the cumulative water production at the surface condition.

In this section, first, a numerical simulation case is generated based on the reservoir and fracture properties in Marcellus Shale to test the accuracy of the proposed flowback models, average pressure calculation method, and the workflow. Next, the proposed models and workflow are applied to an MFHW in Marcellus Shale to characterize hydraulic frac­ tures. The analysis results are compared with long-term production data analysis by history matching with the IHS Harmony (2019) to have a quantitative understanding of temporal changes in fracture properties during the production life of a well.

2.4.2. Workflow for xf, kfi, and γ calculation To obtain xf and kfi using the developed flowback models, the average pressure in the fracture is calculated based on MBE, which re­ quires xf to be known. Thus, an iterative procedure is needed. Here, we present a workflow that calculates fracture properties and evaluates fracture closure using water-phase and gas-phase flowback data. Fig. 5 shows the workflow to estimate fracture attributes iteratively. The 7

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 6. Comparison of average fracture pressure calculated from material balance equation and numerical simulation for (a) Case 1 under constant gas rate and (b) Case 2 under constant water rate.

Fig. 7. Two-phase diagnostic plot of water phase for (a) Case 1 with a constant gas rate and (b) Case 2 with a constant water rate.

3.1. Model validation

regimes depicted in Figs. 2 and 3. The permeability modulus in the fracture is set to be 3 � 10 4 psi 1 according to the laboratory mea­ surements (Ozkan et al., 2010). The input data used in the numerical simulation are given in Table 1. We simulated two cases with different production conditions so that both constant rate and variable rate models can be validated. In Case 1, gas is produced at a constant flowrate of 625 Mscf/day while in Case 2 water is produced at a constant flowrate of 125 bbl/day. First, we validated the modified MBE to calculate average pressure in

We constructed a 2-D two-phase numerical model using IMEX—CMG to generate synthetic data. Relative permeability curve to water and gas in the fracture was assumed to be straight lines. We chose a logarithmic grid size distribution along the fracture and matrix to accurately simu­ late flow near fracture and wellbore. Considering that gas influx usually becomes significant during late flowback period, we conducted the simulation for 200 days, which is long enough to observe the flow

Fig. 8. Water-phase specialty plot for (a) Case 1 with a constant gas rate and (b) Case 2 with a constant water rate. 8

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Table 2 A summary of input data used in numerical simulation and calculated properties using flowback models. Case #

1 2

Set values

Water-phase flowback model

Gas-phase flowback model

Calculated values

Relative error

Calculated values

Relative error

xf (ft)

kfi (md)

xf (ft)

kfi (md)

xf (%)

kfi (%)

xf (ft)

kfi (md)

xf (%)

kfi (%)

500 500

300 300

499.6 499.8

328.0 328.5

0.08 0.04

9.33 9.48

500.3 498.9

319.3 328.5

0.06 0.23

6.44 9.49

Fig. 9. Two-phase diagnostic plot of gas phase for (a) Case 1 with a constant gas rate and (b) Case 2 with a constant water rate.

the fracture with the known xf as an input parameter. A comparison between the analytical approach and numerical simulation is shown in Fig. 6. The results show a perfect match for Case 1 and Case 2, which validate the accuracy of the proposed MBE method. Here, the accuracy of the two-phase flowback model is tested against the production data obtained from the numerical simulation by considering γ as a known parameter. Two-phase water flowback model is validated first. Figs. 7 and 8 show the two-phase diagnostic plot for water and specialty plot for Case 1 and Case 2, respectively. Super­ position pseudotime was used in Case 1 due to the variable flowrate of water. The diagnostic plots for both numerical cases show a half-slope straight line and a unit-slope straight line, which indicate the IALF regime (FR2_W) and the fracture BDF regime (FR3_W), respectively. By selecting the deviation point of the unit-slope straight lines, the start of water-phase BDF regime (FR3_W) was identified to be about 3.41 and 3.5 days for Case 1 and 2, respectively. Within the selected data interval, the slope and intercept of specialty plot are extracted to calculate xf,w and kfi,w. The calculated results and relative error in Table 2 show that the estimated fracture attributes are very close to the set values and the

relative errors are all less than 10%, which validate the accuracy of water-phase flowback model. Here, the gas-phase flowback model is tested against numerical simulation results. We generated the gas-phase diagnostic plot presented by Zhang and Emami-Meybodi (2018) to study the gas-phase flow re­ gimes during flowback in Fig. 9, in which the pseudo variables are evaluated at the average pressure in the distance of investigation of matrix obtained from the numerical simulation. The diagnostic plots for both numerical cases show a clear quarter-slope straight line and a half-slope straight line, which indicate the bilinear flow regime (FR1_G) and the formation IALF (FR2_G), respectively. By selecting the deviation point, the start of gas-phase fracture depletion regime (FR2_G and FR3_G) was identified to be about 1.8 day for Case 1 and 1.9 day for Case 2. Then, we generated the gas-phase specialty plot, as shown in Fig. 10, to conduct FMB analysis. The slope and intercept of the straight line that passes through the late portion of the specialty plot are extracted to calculate xf,g and kfi,g. The calculated results and relative error given in Table 2 show that the evaluated fracture properties are very close to the set values and the relative errors are all less than 10%.

Fig. 10. Gas-phase FMB plot for (a) Case 1 with a constant gas rate and (b) Case 2 with a constant water rate. 9

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

fracture properties are very close to the set values. With the calculated kfi and γ, we can predict the fracture permeability deterioration with the change of average fracture pressure due to fracture closure for the longterm production. Accordingly, the proposed workflow is able to predict fracture properties and evaluate fracture closure. 3.2. Field application We applied the flowback models and workflow to an MFHW drilled in Marcellus Shale in Pennsylvania, USA, to estimate fracture properties and compared the results with long-term production data analysis to characterize the fracture closure during production. The purpose of analyzing this field case is to verify the practical use of the flowback model and to quantitatively evaluate the temporal changes of fracture half-length and permeability. Table 5 shows the real field data from the Marcellus MFHW with 110 hydraulic fractures. Considering that only 40%–60% of the fracture stages are found to contribute the production in practical stimulation job, we set the efficiency of hydraulic fracturing to be 60% for this case. Fig. 12 shows production history of gas flowrate, water flowrate, and bottomhole pressure for this well during flowback and long-term pro­ duction period. After stimulation, the well was flowed back for 24 days and then continuously produced for 249 days without shut-in. In the first 4 days early time flowback period, water flowrate is dominant and re­ mains almost constant, and gas flowrate gradually increases. During the late time flowback, gas production is significant because of the gas influx from matrix, whereas water production quickly declines due to the lack of water contribution. When it comes to long-term production, gas is dominantly produced with negligible water flowrate. We applied the proposed workflow in Fig. 5 to analyze flowback data. By starting the workflow with an initial guess of fracture halflength xf ¼ 700 ft and permeability modulus γ ¼ 1 � 10 3 psi 1, the average pressure in the fracture, two-phase diagnostic plot of water, and

Fig. 11. Calculated average pressure in the fracture during iterations of workflow. Table 3 Validation results of the workflow with the known γ. Iterate xf with the known γ Iteration #

Input xf, ft

Output xf, w, ft

Output xf,g, ft

(xf,w þ xf,g)/ 2, ft

Error of xf, %

1 2 3 4 5

700 488.9 540.0 496.2 498.8

516.5 575.4 506.1 495.7 497.7

461.3 504.7 486.3 501.9 500.7

488.9 540.0 496.2 498.8 499.2

30.16 10.46 8.12 0.52 0.09

The applicability of workflow depicted in Fig. 5 is tested by calcu­ lating xf, kfi, and γ iteratively and comparing them against the set values in Case 1. The proposed workflow updates xf and γ separately by checking the convergence of calculated xf and kfi from the water-phase and gas-phase flowback data. Here the workflow is tested by setting γ as a known parameter to see if the correct xf can be determined. An initial value of xf ¼ 700 ft is chosen to calculate average pressure and then generate specialty plots for water-phase and gas-phase flowback model. The fracture attributes are calculated from the slope and inter­ cept of specialty plots. We compared the average of xf calculated from water and gas flowback data with xf,guess to check the convergence. After five iterations, the relative error meets the acceptable tolerance of 0.2%. The calculated average pressure for each iteration are shown in Fig. 11. As given in Table 3, the estimated xf is very close to the set value after five iterations. Now, we test the proposed workflow by setting the known xf as input to iterate γ. The initial guess of γ starts from 1 � 10 3 psi 1, which is the upper limit of the possible range (10 4–10 3 psi 1). Following the workflow depicted in Fig. 5, the calculated kfi,w and kfi,g are compared to check the tolerance. In the first iteration, the relative error of kfi doesn’t meet the acceptable tolerance 2.7%, thus another γ value is picked as the new input for the second iteration. Similarly, iteration continues and for γ ¼ 3 � 10 4 psi 1 the tolerance is met. Table 4 shows that the calculated

Table 5 Real Field data from a Marcellus MFHW. Parameter

Value

Initial pressure, pi (psi)

6580

Reservoir temperature, T (℉)

140

Initial reservoir pressure, pir (psi) Specific gravity, SGg

Water compressibility, cw (psi

1

4480 0.6

)

2.4 � 10

Initial fracture porosity, ϕfi (%) Matrix porosity, ϕm (%)

7.75

Matrix permeability, km (md)

Fracture compressibility, cf (psi Matrix compressibility, cm (psi

6

60

1 1

1 � 10 )

7 � 10

)

1 � 10

Reservoir thickness, h (ft)

80

Fracture width, wf (ft)

0.04

Initial water formation volume factor, Bwi (bbl/STB)

1.03

Initial fracture water saturation, swi (%)

98

Matrix water saturation, swm (%)

28

Initial water density, ρw (lb/ft3)

62

Initial water viscosity, μw (cp)

5 5 6

0.54

Table 4 Validation results of the workflow with the known xf. Iterate γ with the known xf Iteration #

Input γ, psi

1 2 3 4

1� 7� 5� 3�

10 10 10 10

3 4 4 4

1

Output xf,w, ft

Output xf,g, ft

xf error, %

Output kfi,w, md

Output kfi,g, md

kfi error, %

501.2 500.3 499.8 499.6

503.1 501.9 501.0 500.3

0.39 0.30 0.25 0.14

336.0 333.2 331.4 328.0

326.2 323.5 321.7 319.3

2.92 2.92 2.93 2.64

10

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

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

Fig. 13. Plots during flowback analysis after the convergence of xf and kfi: (a) Average pressure and permeability in the fracture during flowback period. (b) Twophase diagnostic plot of water where the pseudo variables are calculated based on the average pressure from plot (a). (c) Specialty plot of two-phase water flowback model. (d) Specialty plot of two-phase gas flowback model.

specialty plots for water and gas are generated during each iteration. The tolerances of the relative error between xf,new and xf,guess, as well as kfi,w and kfi,g are both set to be 10%. Fig. 13 presents the essential plots of flowback analysis after the convergence of both xf and kf. Fig. 13a shows the average fracture pressure obtained from the modified material bal­ ance equation and the average fracture permeability history during

flowback which is determined based on the iterated γ and kfi. Fig. 13b provides the log-log plot of DRNPw vs. tsp,w to identify FR3_W regime (Zhang and Emami-Meybodi, 2020). The half slope straight line in the diagnostic plot shows the FR2_W regime in the early time flowback period. A dominant unit-slope straight line is identified after 2.75 days flowback on the diagnostic plot, which shows the boundary dominated flow in the fracture. It can be seen that the flow regimes sequence for this well is consistent with our conceptual model in Fig. 2, which verifies the practical use of the developed models and proposed workflow. With the identified flow regime, water and gas specialty plots are generated, Fig. 13c and d, to calculate fracture properties based on the slope and intercept of the straight line. The flowback analysis results are listed in Table 6. We also analyzed the long-term production data to obtain the

Table 6 Analysis results using the flowback model and history matching. Fracture Properties

Flowback analysis

Long-term production data analysis

xf, ft kfi, md γ, psi 1 Average kf, md

873 1227.5 8 � 10 –

489 – – 46.4

4

11

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Fig. 14. History match of (a) gas flowrate and (b) bottomhole pressure during long-term production using the Horizontal-Multifractured-SRV model in IHS har­ mony software.

fracture properties in the late time. We use history matching approach to analyze the long-term production data with the single-phase gas Horizontal-Multifractured-SRV model in IHS Harmony software. Fig. 14 shows a good match of gas flowrate and bottomhole pressure history. Although the history match result is not unique, it is still deemed good enough to yield accurate results as long as the model parameters are in acceptable ranges (Clarkson et al., 2019). Table 6 summarizes the analysis results for the field case using the new flowback model and history matching approach. By comparing the flowback analysis with long-term production data analysis results, we observed a significant decrease in xf from 873 ft to 489 ft and in kf from 1227.5 md to 46.4 md during the production, which is mainly due to the fracture closure. Fig. 14b shows that the bottomhole pressure declines slowly after 3 months of production, which is around 1200 psia. Applying the obtained kfi and γ from flowback analysis, we estimated the ultimate kf to be 16.6 md when average fracture pressure reaches 1200 psia. The ultimate kf 16.6 md is slightly below the average kf 46.4 md obtained from history matching of long-term production data, which is reasonable according to eq (6) because the estimated ultimate fracture pressure 1200 psia is lower than the average fracture pressure during long-term production. Furthermore, the predicted ultimate kf 16.6 md falls into the permeability range between the fracture without proppant and with one layer of glass beads, which is obtained from laboratory experiments by Tan et al. (2018). Accordingly, the decrease in fracture permeability and fracture half-length is mainly caused by the removal of proppants and fracture closure.

the two-phase flowback models can closely predict fracture properties. We applied the workflow to flowback data from a multi-fractured hor­ izontal gas well in Marcellus Shale and utilized IHS Harmony software to analyze long-term production data. With the estimated fracture prop­ erties from flowback and long-term production data, the fracture closure was quantitatively evaluated. The Marcellus field application shows that the fracture half-length is reduced to half and the fracture permeability decreased by over an order of magnitude after nine months of production. We assumed a uniform gas influx along the fracture in the gas flowback model, which is valid when the pressure gradient along the fracture is small. If the system has a substantial pressure gradient in the fracture, we suggest considering non-equivalent gas influx in the flow­ back model. Another assumption made in this study is negligible water influx from matrix, which is valid when the water saturation is small, or the capillary pressure impedes water flow in the matrix. However, when the water leakoff during hydraulic fracturing is significant, the absorbed water in the shale matrix has a significant impact on flowback, so the two-phase flow is likely to appear in the near fracture region. Hence further study is needed to include the effect of water influx from matrix into the flowback models, as well as gas desorption and diffusion for shale gas reservoir (Xu et al., 2016; Zhang et al., 2018; Yang et al., 2017; Clarkson et al., 2016, 2017; Hamdi et al., 2018; Behmanesh et al., 2018a). Further, the modeling of two-phase bilinear flow in FR1_G needs further investigations, which is outside of the scope of this study. Declaration of competing interest

4. Summary and conclusions

The authors declare no conflict of interest.

We developed multiphase flowback models for both gas and water boundary dominated flow (BDF) by considering gas influx from matrix under both constant and variable-rate conditions. The flowback models consider pressure-dependent fracture permeability. We also developed a two-phase flowing material balance (FMB) equation to calculate average pressure in the fracture by using the variable BHP solution of gas inflow equation in the matrix. We proposed a workflow to reconcile waterphase and gas-phase flowback models and evaluate fracture closure through estimation of initial fracture permeability, fracture half-length, and fracture permeability modulus. We validated the flowback models and workflow against the results obtained from numerical simulations. The validation results reveal that

Acknowledgments The authors acknowledge the financial support from the Penn State Institutes of Energy and the Environment for the PSIEE Seed Grant and the Penn State College of Earth and Mineral Sciences for the 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 Modeling Group, Inc. (CMG) and IHS Markit.

Nomenclatures ABT

after-breakthrough 12

F. Zhang and H. Emami-Meybodi

BBT BDF CBM bFMB bpss bw Bg Bw ceff cf cgf cgm cm ctm cw DD DDA DRNPw DRP EGP fcp FMB Gg;inf Gp GpN GWR h IALF kf kfi km krg krw lf LGP mðpÞ mðpf Þ mðpf Þ mFMB mg mm ðpÞ mw M MBE MFHWs p pb pf pf pi pp ppi ppwf psc pwf PNR qg qgf qg;inf qw R RNPg RNPw

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

before-breakthrough boundary dominated flow coalbed methane the intercept of GPnorm vs. PNR plot the intercept of tsp;g vs. RNPg plot the intercept of tsp;w vs. RNPw plot gas formation factor, ft 3 =scf water formation factor, bbl=STB effective gas compressibility, psi 1 fracture compressibility, psi 1 gas compressibility in the fracture, psi 1 gas compressibility in the matrix, psi 1 matrix compressibility, psi 1 matrix total compressibility, psi 1 water compressibility, psi 1 dimensionless drawdown parameter, dynamic drainage area derivative of water-phase rate normalized pseudopressure, psia=STB dynamic relative permeability early gas production correction factor, flowing material balance cumulative gas production of gas influx at surface condition, scf cumulative gas production, scf normalized cumulative gas production, scf gas water ratio, Mscf/STB fracture/formation height, ft infinite acting linear flow fracture permeability, md initial fracture permeability, md matrix permeability, md gas relative permeability, water relative permeability, fracture spacing, ft late gas production gas-phase pseudopressure, psia gas-phase pseudo average pressure in the fracture, psia average gas-phase pseudo pressure in the fracture, psia the slope of GPnorm vs. PNR plot the slope of tsp;g vs. RNPg plot real gas pseudopressure, psia2 =cp the slope of tsp;w vs. RNPw plot mass of gas, lb/lbmol material balance equation multi-fractured horizontal wells pressure, psia base pressure, psia fracture pressure, psia average pressure in the fracture, psia initial fracture pressure, psia water-phase pseudopressure, psia initial water pseudopressure, psia bottomhole water pseudopressure, psia pressure at standard condition, psia bottom-hole pressure, psia pressure normalized rate, psia⋅day=Mscf the gas-phase surface flowrate, Mscf=day the gas-phase downhole flowrate, Mscf=day the flowrate of gas influx at surface condition, Mscf=day the water-phase surface flowrate from the whole fracture, STB=day ideal gas constant, psia⋅ft3/(lbmol⋅R) gas-phase rate normalized pseudopressure, psia⋅day=Mscf water-phase rate normalized pseudopressure, psia⋅day=STB 13

F. Zhang and H. Emami-Meybodi

RTA sg sg sw swm sw t tp;g tp;w tsp;g tsp;w T Tsc Vfi wf Wfi Wp x xf y Z

ε ϕf ϕm

ρw ρgf ρgm μw μgf μgm

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

rate transient analysis gas saturation in the fracture, average gas saturation in the fracture, water saturation in the fracture, water saturation in the matrix, average water saturation in the fracture, time, day gas-phase pseudotime, day water-phase pseudotime, psia⋅day=cp gas superposition pseudotime, day water superposition pseudotime, psia⋅day=cp temperature, R temperature at standard condition, R initial pore volume of half fracture, ft3 hydraulic fracture width, ft initial water storage in fracture, scf cumulative water production, STB location in the fracture, ft fracture half-length, ft location in the matrix, ft gas deviation factor, error tolerance, fracture porosity, matrix porosity, water density, lb/ft3 gas-phase density in the fracture, lb/ft3 gas-phase density in the matrix, lb/ft3 water viscosity, cp gas viscosity in the fracture, cp gas viscosity in the matrix, cp

γ

fracture permeability modulus, psi

Subscripts i j f m g w D

initial the jth flowrate interval fracture matrix gas water dimensionless

1

References

concept. In: SPE Western Regional Meeting, SPE-185724. Society of Petroleum Engineers. Clarkson, C., 2013. Modeling two-phase flowback of multifractured horizontal wells completed in shale. SPE J 18 (4), 795–812. 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., Jordan, C.L., Gierhart, R.R., Seidle, J.P., 2008. Production data analysis of coalbed-methane wells. SPE Reserv. Eval. Eng. 11, 311–325. Clarkson, C.R., Qanbari, F., 2016. 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., 2016. A semi-analytical method for forecasting wells completed in low permeability, undersaturated CBM reservoirs. J. Nat. Gas Sci. Eng. 30, 19–27. Clarkson, C.R., Qanbari, F., Williams-Kovacs, J.D., 2016. Semi-analytical model for matching flowback and early-time production of multi-fractured horizontal tight oil wells. J. Unconv. Oil Gas Resour. 15, 134–145. Clarkson, C., Williams-Kovacs, J., 2013. A new method for modeling multi-phase flowback of multi-fractured horizontal tight oil wells to determine hydraulic fracture properties. In: SPE Annual Technical Conference and Exhibition, SPE-166214. Society of Petroleum Engineers. 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: Unconventional Resources Technology Conference. URTEC-2019-85.

Abbasi, M., Dehghanpour, H., Hawkes, R.V., 2012. Flowback analysis for fracture characterization. In: SPE Canadian Unconventional Resources Conference, SPE162661. Society of Petroleum Engineers. 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, SPE-168982. Society of Petroleum Engineers. Agarwal, R.G., Gardner, D.C., Kleinsteiber, S.W., Fussell, D.D., 1999. Analyzing well production data using combined type curve and decline curve analysis concepts. SPE Reserv. Eval. Eng. 2 (5), 478–486. Al-Kobaisi, M., Ozkan, E., Kazemi, H., 2006. A hybrid numerical/analytical model of a finite-conductivity vertical fracture intercepted by a horizontal well. SPE Reserv. Eval. Eng 9 (4), 345–355. Behmanesh, H., Hamdi, H., Clarkson, C.R., Thompson, J., Anderson, D., Chin, A., 2018. Transient linear flow analysis of multi-fractured horizontal wells considering threephase flow and pressure-dependent rock properties. In: Unconventional Resources Technology Conference. URTEC-2884255. Behmanesh, H., Mattar, L., Thompson, J.M., Anderson, D.M., Nakaska, D.W., Clarkson, C. R., 2018. Treatment of rate-transient analysis during boundary-dominated flow. SPE J. 23 (1), 145-1,165. Clarkson, C., Jordan, C., Ilk, D., Blasingame, T., 2012. Rate-transient analysis of 2-phase (gasþ water) CBM wells. J. Nat. Gas Sci. Eng. 8, 106–120. Clarkson, C., Qanbari, F., Williams-Kovacs, J., Zanganeh, B., 2017. Fracture propagation, leakoff and flowback modeling for tight oil wells using the dynamic drainage area

14

F. Zhang and H. Emami-Meybodi

Journal of Petroleum Science and Engineering xxx (xxxx) xxx

Crafton, J.W., 2008. Modeling flowback behavior or flowback equals“ slowback”. In: SPE Shale Gas Production Conference, SPE-119894. Society of Petroleum Engineers. Dake, L.P., 1983. Fundamentals of Reservoir Engineering. Elsevier. Ezulike, D., Dehghanpour, H., Hawkes, R., 2013. Understanding flowback as a transient 2-phase displacement process: an extension of the linear dual-porosity model. In: SPE Unconventional Resources Conference Canada, SPE-167164. Society of Petroleum Engineers. Ezulike, O.D., Dehghanpour, H., 2014. Modelling flowback as a transient two-phase depletion process. J. Nat. Gas Sci. Eng. 19, 258–278. Fu, Y., Dehghanpour, H., Motealleh, S., Lopez, C.M., Hawkes, R., 2019. Evaluating Fracture Volume Loss during Flowback and its Relationship to Choke Size: Fastback versus Slowback. SPE Prod. Oper 34 (3), 615–624. Ghanbari, E., Abbasi, M.A., Dehghanpour, H., Bearinger, D., 2013. Flowback volumetric and chemical analysis for evaluating load recovery and its impact on early-time production. In: SPE Unconventional Resources Conference Canada, SPE-167165. Society of Petroleum Engineers. Hamdi, H., Behmanesh, H., Clarkson, C.R., 2018. A semi-analytical approach for analysis of the transient linear flow regime in tight reservoirs under three-phase flow conditions. J. Nat. Gas Sci. Eng. 54, 283–296. Hossain, S., Ezulike, O., Dehghanpour, H., 2019. Estimating residual fracture pore volume by analyzing post-flowback water production: an Eagle Ford black-oil case. In: Unconventional Resources Technology Conference. URTEC-2019-989. Ibrahim, M., Wattenbarger, R.A., 2006. Rate dependence of transient linear flow in tight gas wells. J. Can. Pet. Technol. 45, 18–20. IHS Harmony, 2019. HarmonyTM & Harmony EnterpriseTM. Well Performance Software. https://ihsmarkit.com/products/harmony-oil-well-performance-software.html. 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, 18(4). Society of Petroleum Engineers, SPE-135607. 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. 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 hydraulicallyfractured shale gas wells using a hybrid numerical/analytical model. Int. J. Coal Geol. 182, 14–31. Jones, M.R., Blasingame, T.A., 2019. A direct method for short-term forecasting of multiphase production rates using flowback data. In: Unconventional Resources Technology Conference. URTEC-2019-178. Kanfar, M., Alkouh, A.B., Wattenbarger, R.A., 2014. Bilinear flow in shale gas wells: fact or fiction?. In: SPE/CSUR Unconventional Resources Conference–Canada, SPE171669. 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. Kurtoglu, B., Salman, A., Kazemi, H., 2015. Production forecasting using flow back data. In: SPE Middle East Unconventional Resources Conference and Exhibition, SPE172922. Society of Petroleum Engineers. Miller, F.G., 1962. Theory of unsteady-state influx of water in linear reservoirs. J. Inst. Pet. 48, 365–379. Nobakht, M., Clarkson, C.R., 2012. Analysis of production data in shale gas reservoirs: rigorous corrections for fluid and flow properties. J. Nat. Gas Sci. Eng. 8, 85–98. Nobakht, M., Clarkson, C.R., 2012. A new analytical method for analyzing linear flow in tight/shale gas reservoirs: constant-rate boundary condition. SPE Reserv. Eval. Eng. 15, 51–59.

Nobakht, M., Clarkson, C.R., 2012. A new analytical method for analyzing linear flow in tight/shale gas reservoirs: constant pressure production. SPE Reserv. Eval. Eng. 15, 370–384. € Ozisik, M.N., 1993. Heat Conduction. John Wiley & Sons. 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, SPE-134830. Society of Petroleum Engineers. Palacio, J., Blasingame, T., 1993. Decline curve analysis using type curves—analysis of gas well production data. Low Permeability Reservoirs Symposium, 26-28 April, Denver, Colorado SPE-25909. Pedrosa JR., O.A., 1986. Pressure transient response in stress-sensitive formations. In: SPE California Regional Meeting, SPE-15115. Society of Petroleum Engineers. Shahamat, M.S., Clarkson, C.R., 2018. Multiwell, multiphase flowing material balance. SPE Reserv. Eval. Eng. 21, 445–461. Shahamat, M.S., Mattar, L., Aguilera, R., 2015. A physics-based method to forecast production from tight and shale petroleum reservoirs by use of succession of pseudosteady states. SPE Reserv. Eval. Eng 18, 508–522. Sun, Z., Shi, J., Zhang, T., Wu, K., Miao, Y., Feng, D., Sun, F., Han, S., Wang, S., Hou, C., 2018. The modified gas-water two phase version flowing material balance equation for low permeability CBM reservoirs. J. Pet. Sci. Eng. 165, 726–735. Tan, Y., Pan, Z., Liu, J., Feng, X.-T., Connell, L.D., 2018. Laboratory study of proppant on shale fracture permeability and compressibility. Fuel 222, 83–97. Wattenbarger, R.A., Alkouh, A.B., 2013. New advances in shale reservoir analysis using flowback data. In: SPE Eastern Regional Meeting, SPE-165721. Society of Petroleum Engineers. 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, SPE-39931. Society of Petroleum Engineers. Williams-Kovacs, J., Clarkson, C., Zanganeh, B., 2015. Case studies in quantitative flowback analysis. In: SPE/CSUR Unconventional Resources Conference, SPE175983. Society of Petroleum Engineers. Williams-Kovacs, J.D., Clarkson, C.R., 2016. A modified approach for modeling twophase 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, 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. Yang, Y., 2017. Mathematical Development for Flowback Rate Transient Analysis. MSc Thesis. Pennsylvania State University. Zhang, F., Emami-Meybodi, H., 2018. Flowback fracture closure of multifractured horizontal wells in shale gas reservoirs. In: SPE/AAPG Eastern Regional Meeting, SPE-191817. Society of Petroleum Engineers. Zhang, Y., Ehlig-Economides, C., 2014. Accounting for remaining injected fracturing fluid in shale gas wells. In: Proceedings of the 2nd Unconventional Resources Technology Conference. URTEC-1892994. Zhang, F., Emami-Meybodi, H., 2020. Multiphase flowback rate-transient analysis of shale gas reservoirs. Int. J. Coal Geol. 217 (103315). Zhang, Z., Yuan, B., Ghanizadeh, A., Clarkson, C.R., Williams-Kovacs, J., 2018. Rigorous estimation of the initial conditions of flowback using a coupled hydraulic fracture/ dynamic drainage area leakoff model constrained by laboratory geomechanical data. In: Unconventional Resources Technology Conference. URTEC-2901771.

15