Transient pressure behavior for unconventional gas wells with finite-conductivity fractures

Fuel 266 (2020) 117119

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Transient pressure behavior for unconventional gas wells with finiteconductivity fractures

T

Kien H. Nguyena, Miao Zhanga, , Luis F. Ayalab ⁎

a b

New Mexico Institute of Mining and Technology, United States The Pennsylvania State University, United States

ABSTRACT

Horizontal wells and hydraulic fracturing technologies are the only way that makes commercial production from the low and ultra-low permeability unconventional plays possible. The development of analytical and semi-analytical tools to analyze transient behavior of multi-fractured horizontal gas wells have been largely biased towards solving the simplified problem of one-dimensional (1D) Cartesian gas flow in matrix domain, which are only applicable to the production from infiniteconductivity fractures. This study presents a novel semi-analytical method for transient behavior analysis of horizontal gas wells producing from finite-conductivity fractures. Unlike the available transient analysis tools for finite-conductivity fractures which are developed based on empirical and approximate solution, the proposed solution is rigorously derived by solving the governing gas diffusivity equations written for gas flow in two contiguous regions of matrix and fracture simultaneously. Nonlinear, pressure-dependent gas properties remaining in pseudopressure-based diffusivity equations—namely viscosity and compressibility—are rigorously and straightforwardly captured in both fracture and matrix systems. The validity of proposed solution is verified against commercial numerical simulator (COMSOL Multiphysics®) for a series of synthetic cases considering constant and variable -rate production constraints. The case studies presented in this paper also showcase the applicability of proposed solution to the analysis of transient behavior of multi-fractured horizontal gas wells in shale plays under the following two state-of-the-art interpretations: 1) production from finite-conductivity hydraulic fractures in a single-porosity matrix; and 2) production from infinite-conductivity hydraulic fractures in a naturally-fractured matrix.

1. Introduction Unconventional gas reservoirs, including shale gas, tight gas and associated gas from tight oil, constitute a significant source of natural gas production in the United States. In the ultra-low permeability environment found in these unconventional plays, multi-transverse hydraulic fractures in horizontal wells have become a standard production strategy. Accurate and reliable analysis and forecast of the pressure and/or rate transient behavior of these multi-fractured horizontal wells (MFHWs) are important components of the assessment of economic feasibility of tight and shale resources exploitation. Pressure transient analysis (PTA), rate transient analysis (RTA) and production data analysis (PDA) are among the most popular techniques to analyze and predict well performance as well as estimation of reservoir characteristics. The development of these inverse analysis tools heavily relies on the availability of analytical solutions to the governing fluid flow equations (diffusivity equations) applied to the initial and boundary conditions of the system of interest. These analytical tool have been traditionally available for single-phase oil (slightly compressible) reservoirs, for which closed-form, analytical solutions can be derived for the linear partial differential equations (PDE) of liquid flow (e.g. [14,9]; Chen and Raghavan, 1997; [7]). However, the PDE of gas ⁎

(compressible) flow presents strong nonlinearities due to the presence of pressure-dependent gas properties, and thus it generally lacks rigorous and exact analytical solutions. In gas applications, pseudo-pressure [2] is routinely used to group nonlinear pressure-dependent terms with pressure, but its use always leads to a residual presence of pressure-dependent variables in the accumulation term of the gas flow equation that hinder their analytical treatment. As a result, classical PTA or RTA tools that have been developed for vertical gas wells completed in conventional reservoirs normally apply the following approximations: 1) at early-time of production (infinite-acting conditions), the remaining nonlinear gas coefficient in the diffusivity equation is assumed to be constant and equal to its initial value; 2) at latetime of production (after reservoir boundaries are felt), pseudo-time function [13] is deployed to evaluate the nonlinear gas diffusivity term at a volume-averaged condition of the reservoir. Both of these approximations have been observed to not directly applicable to analysis of MFHWs’ performance in unconventional tight or shale reservoirs. Neglecting the nonlinearity of gas diffusivity in early-transient infiniteacting periods has shown to lead to erroneous results [18,21,34]. Empirical, approximate, and semi-analytical methods have been proposed to improve the predictions by capturing nonlinear gas diffusivity during early-transient periods of MFHWs [18,21,24,32,34,35]; but they are all

Corresponding author. E-mail address: [email protected] (M. Zhang).

https://doi.org/10.1016/j.fuel.2020.117119 Received 21 September 2019; Received in revised form 7 December 2019; Accepted 15 January 2020 Available online 27 January 2020 0016-2361/ © 2020 Elsevier Ltd. All rights reserved.

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

Nomenclature Aw CfD cϕ cg ct fa G h km kF kf LF Lf MW M mi mD nt nc p pi Qs

QFd QmD R SG T T tD wF wf x

cross-sectional area, ft2 dimensionless fracture conductivity rock compressibility, psi−1 gas compressibility, psi−1 total compressibility, psi−1 analytical adjustment for gas flow, dimensionless Free-space Green’s function reservoir thickness, ft matrix permeability, md or µd hydraulic fracture permeability, md or µd natural fracture permeability, md or µd hydraulic fracture half length, ft half of the hydraulic fracture spacing in naturally-fractured matrix, ft molecular weight, lbmass/lbmol normalized real gas pseudopressure, psi normalized real gas pseudopressure at initial condition, psi dimensionless pseudopressure number of discretized time element number of discretized space element pressure, psia initial reservoir pressure, psia source strength, lbmass/day/ft3

(A)

Closed boundary

y Z

dimensionless source strength in fracture dimensionless source strength in matrix universal constant gas specific gravity reservoir temperature, °F or °R time, day dimensionless time hydraulic fracture width, ft natural fracture width, ft linear distance under Cartesian system in matrix domain, ft linear distance under Cartesian system in fracture domain, ft gas compressibility factor

Greek λ µg ρg ϕm ϕf ϕF δ ηD

Hydraulic fracture

gas viscosity-compressibility ratio, dimensionless gas phase viscosity, cp gas phase molar density at reservoir conditions, lbmass/ft3 matrix porosity natural fracture porosity hydraulic fracture porosity Dirac delta function dimensionless fracture storage capacity

Matrix

LF

Horizontal well

(B)

Closed boundary

Hydraulic fracture

Matrix Natural fracture

Lf

Horizontal well

Fig. 1. (A) Schematic of a multi-fractured horizontal well with single-porosity matrix; (B) schematic of a multi-fractured horizontal well with naturally-fractured matrix. 2

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

limited to the most simplified scenario of one-dimensional Cartesian (linear) flow in matrix domain, and only apply to the production infinite-conductivity hydraulic fractures where pressure drop and transient within fractures are neglected. For the horizontal gas wells with finite-conductivity fractures, most recently, Kanfar and Clarkson [16] proposed empirical corrections to bilinear liquid flow solution in order to capture the retaining nonlinear gas diffusivity in both fracture and matrix flow equations. However, their approach is currently empirical and strictly applicable to bilinear flow conditions only, represented by ¼ slope in diagnostic plot. This study proposes a rigorous, semi-analytical solution applicable to the analysis of transient pressure behavior for horizontal gas wells with finite-conductivity fractures. The early-time fracture Cartesian (linear) flow regime, bilinear flow regime, late-time matrix linear flow regime, as well as the transitions between these flow regimes are all rigorously modelled. The proposed solution for finite-conductivity fractures is developed based on extending the Green’s function-based solution proposed by Zhang and Ayala [34] for gas horizontal wells with infinite-conductivity fractures. Before Zhang and Ayala [34], the Green-function-based integral solution for gas reservoirs have been successfully applied to analyze the well test data of vertical gas wells [4,5,20,26,27] and unfractured horizonal gas wells [36]. While these aforementioned works are all focused on solving one governing gas flow equation written for single domain (matrix) only, this present work takes one step further to analyze the behavior of gas horizontal wells with finite-conductivity fractures in which two domains of matrix and fracture have to be solved simultaneously.

1D matrix flow (black streamlines) and 1D HF flow (white streamlines) are solved simultaneously in Fig. 1A. Moreover, we also show the proposed finite-conductivity fracture solution can be used to model the horizontal gas wells in naturally-fractured reservoirs under the context of ‘triple-porosity (dual-fracture) model’ as depicted in Fig. 1B. In this model, naturally-fractured matrix is modelled using dual porosity approach, in which rectilinear prism matrix is separated by an orthorhombic continuum of natural (micro) fractures [31,17,37]. Besides matrix and natural fracture domains, hydraulic(macro) fractures are considered as third domain with distinct properties, hence “triple-porosity (dual-fracture) model” [3,12] or “multi-porosity model” [37,30]. In this manuscript, we honor the classical triple-porosity model by assuming hydraulic fractures producing to the well while they are fed by the natural fractures only. Consequently, the matrix feeds the natural fractures only. Therefore, the inter-porosity flow is sequential from one medium to the other, as shown by the red, black and white streamlines in Fig. 1B. Al-Ahmadi and Wattenbarger [3] first employed the porosity model to analyze the production data of hydraulic-fractured gas shales with existing natural fractures (NF) in which the three inter-porosity flows (matrix to NF to HF to HW) are all assumed to be under transient (unsteady-state) period. Tivayanonda et al. [29] later showed that these three flow regions can be simplified to the consideration of only matrix to NF to HF by assuming infinitely conductive hydraulic fractures—which can be straightforwardly modelled using proposed solution. Notably, in both Al-Ahmadi and Wattenbarger [3] and Tivayanonda et al. [29], pseudo-pressure-based gas flow equations are analytically solved without the consideration of nonlinearity effect caused by remaining pressure-dependent gas diffusivity term. Kanfar and Clarkson [16] improved Tivayanonda et al.’s solution by proposing empirical corrections to account for pressure-dependent gas diffusivity but strictly for rate-normalized-pressure versus fourth-root-time plot and thus limited to bilinear flow duration only.

2. Horizontal gas wells with finite-conductivity fractures Zhang and Ayala [34] proposed a Green-function-based semi-analytical solution to simulate the pressure transient behavior of gas horizontal wells (HW) with multi-transverse hydraulic fractures (HF) and single-porosity matrix shown as the Fig. 1A below. The hydraulic fractures in Zhang and Ayala [34] are assumed to have infinite conductivity and uniform pressure distribution, and thus the gas flow system can be characterized as one-dimensional from matrix to infiniteconductivity HF only (represented by the black streamlines in Fig. 1A). This present work extends the solution to finite-conductivity HFs where

3. Semi-analytical solution for gas flow in unconventional reservoirs with Finite-conductivity fractures The semi-analytical solution proposed in this study can be applied to transient gas flow in unconventional reservoirs for the case of gas HW

y No flow boundary

Plane source (Fracture)

Matrix Lf

Q ( y ,t )

p ( x → −∞ ; t ) = p i

p ( x → +∞; t ) = pi

x Intersection of HF - HW(Fig. 1A) or NF - HF(Fig. 1B)

Fig. 2. Flow regions and boundary conditions of interest in this study.

3

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

with finite-conductivity HF (Fig. 1A) or gas HF with finite-conductivity NF (Fig. 1B). The simplified system of interest under isotropic and homogeneous conditions, for which the proposed solution is derived, consists flow regions and boundary conditions shown in Fig. 2. Additional assumptions behind the representation in Fig. 2 are explained in detail for each of the applications of interest (Fig. 1A or B) below. For the application to multi-fractured horizontal gas wells with single-porosity matrix and finite-conductivity HFs (Fig. 1A), the representative flow model to be studied is single-phase gas flow from matrix to finite-conductivity HF to horizontal wells, which represents the basic building block of the sequence in Fig. 1A. The following assumptions are employed: 1) fully penetrating horizontal well produces at a rate-specified constraint (either constant or variable specification); 2) hydraulic fractures are uniformly distributed; 3) fracture width is small compared with fracture length, which allows us to not only consider 1D Cartesian flow within the fracture but treat fractures as plane sources with uniform flux densities along height (h) and width (wF); 4) flow beyond fracture tip is neglected (no flow at top-boundary); 5) matrix flow is infinite-acting (i.e., HF interference not considered). For the application to multi-fractured horizontal gas wells with naturally-fractured matrix (Fig. 1B), the representative flow model to be studied is single-phase gas flow from matrix to finite-conductivity NF to infinite-conductivity HFs which is the building block of the sequence in Fig. 1B. Additional assumptions include: 1) fully penetrating horizontal well produces at rate-specified constraint (either constant or variable specification); 2) both HFs and NFs are uniformly distributed; 3) natural fracture width (wf) is small compared with fracture length Lf (which is referred as HF spacing in many literatures such as Al-Ahmadi and Wattenbarger [3], Tivayanonda et al. [29] and Kanfar and Clarkson [16]; 4) no flow boundary is considered at top-boundary where symmetry line is present (dash line shown in Fig. 1B); 5) matrix flow is infinite-acting (i.e., any interference between NFs is neglected).

where µg is gas viscosity, ct is total compressibility and equal to the summation of gas compressibility (cg) and pore compressibility (cφ). Because gas viscosity and compressibility are strong functions of pressure, has been shown to have significant impact on the recovery performance of natural gas reservoirs, including late-time behavior of conventional gas wells exhibiting radial-cylindrical flow [13,22,33] and early-transient behavior of MFHWs exhibiting linear or linear flow [24,32]. In Zhang and Ayala [34], a semi-analytical solution was proposed to solve gas 1D Cartesian flow problem (Eq.3) in which is captured in a rigorous and direct manner by means of an adjustment factor fa . The solution to Eq. (3) at fracture location (x = 0) is expressed as: where mD* is the analytical solution of Eq. (3) when is assumed to be constant and equal to 1, and which represents the classical Duhamel’s principle applied to the transient behavior in an infinite, 1D Cartesian system, as shown below:

mfD =

fa (tD) =

x

m g)

t

+ Qm ( x

x 0)

QmD ( ) 4 (tD

)

d

(6)

tD

d

+

1 ( , )

0

1

mD ( , ) G d tD

(7)

The free-space Green function for one-dimensional Cartesian flow, G (xD , tD | , ) , is the system response of location xD at time tD to an instantaneous source of unit strength happening at location and time . The expression at fracture location (xD = 0) and is written as (Cole et al., 2011):

Given the coordinate system and representation shown in Fig. 2, the governing diffusivity equation for 1D linear gas flow in matrix along x direction can be expressed as:

(

tD 0

QmD is the dimensionless form of the plane source strength of infinite-conductivity fracture, which can be calculated in terms of specified flow rate as shown in Appendix A. fa is the adjustment factor that quantifies the deviation of system response due to the nonlinear effect of in the matrix domain:

4. Matrix linear flow solution

km p = g µg x

(5)

mfD = mfD + fa

G (tD; , ) =

1 4 (tD

)

exp

2

4(tD

)

(8)

This Green-function-based semi-analytical solution (Eqs. (5) to (8)) is derived for 1D Cartesian flow along x direction under infinite-acting

(1)

where Qm* is source strength in mass per unit volume per unit time, δ is Dirac delta function, and x0 gives location of continuous source with strength Qm*. Eq. (1) can be partially linearized using normalized gas pseudo-pressure [2,22]:

µgi Zi

m=

pi

p 0

p dp µg Z

(2)

The dimensionless form of Eq. (1), using the dimensionless group proposed by Zhang and Ayala [34] and provided given in Appendix A, can be written as: 2m D xD2

=

1 mD + QmD (xD tD

xD0)

where mD is the dimensionless pseudopressure m. dependent viscosity-compressibility ratio of gas:

=

(3) is the pressure-

µgi cti µg (p) ct (p)

(4)

Fig. 3. Schematic of the coupled 1D Cartesian flow regions simulated in this study.

4

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

conditions, thus is limited to the cases in which infinite-conductivity fracture is present. In this study, we propose the extension of this semianalytical solution to analyze finite-conductivity fractures by dividing the fracture into N elements as depicted in Fig. 3. For each fracture element, the flow from matrix domain into fracture is assumed to be one-dimensional and along x-direction only. Thus, the response of i-th fracture element at n-th time step is simply given by applying Eq. (5) to each element: (n ) (n ) (n ) mfD , i = mfD, i + f a, i ; i = 1, 2, 3, ...,N

D

tD

=

0

n

QmD, i ( ) 4 (tD

)

d

(j ) QmD ,i

j=1

(tD, n

tD, j )

(n ) QfD ,i =

(11)

(

(n ) QfD ,i =

QfD =

1 D

mD tD

1

(n ) QfD ,i =

1 (n ) D i

(n ) mfD ,i

(n 1) mfD ,i

tD

; (16)

Qs(,ni ) RTµgi Zi Lf2 wf MWpi kf mi

(17)

(n ) QmD ,i

CfD

(18)

where CfD is the dimensionless fracture conductivity [8]:

CfD =

wf kf L f km

(19)

Rate-specified production condition at the interaction/connection interface (point of origin in Fig. 3) is treated as boundary condition of Eq.16. Plus the no-flow boundary at the top tip of fracture, 2N primary unknowns (mfD,i and QmD,i) can be solved simultaneously at each time step after coupling Eq. (9) and (16). , Iterations are implemented to rigorously account for the dependency of the nonlinear λ parameter on pressure. In this study, a simple Newton-Raphson procedure is deployed to obtain the final solutions. Appendix C provides a detailed flow chart for this solution method. It is worth reiterating once more that the proposed semi-analytical solution embraces the assumption that the matrix flow is 1D and along x direction. It is considered to be a good approximation when the reservoir permeability is smaller than fracture permeability by at least several magnitudes, which is a typical condition found in unconventional tight and ultra-tight formations. A more rigorous approach is to consider matrix flow under 2D Cartesian flow geometry, a model that was employed by Cinco-Ley et al. [8] to solve for the transient behavior of fractured vertical oil wells. Cinco-Ley et al. [8] coupled 2D Greenfunction-based solution of matrix flow with 1D finite-difference model of fracture flow, which are both developed for single-phase liquid system, to solve for primary unknowns at fracture elements. For the purpose of this study, we assume 1D matrix flow so that the semianalytical solution of 1D gas flow proposed for infinite-conductivity fractures by Zhang and Ayala [34] is readily applicable. The extension to 2D gas flow in matrix domain, however, requires further modifications of adjustment factor fa (Eq.7) to rigorously capture the nonlinearity effect in a 2D flow geometry, and will be further discussed in our future work.

(12)

The right-hand-side term (RHS) accounts for expansion of fluid inside fracture (fracture storage term). However, since the fracture’s volume is significantly smaller than matrix and its contribution to production rate could be considered negligible, many studies [15,34] are able to treat the RHS of Eq. (12) as zero. In Appendix B of this paper, we compare our modeling results which consider and do not consider the right-hand-side accumulation term, and it is demonstrated that both approaches differ by an extremely small amount. For simplicity, subscript ‘f’ denotes fracture in the derivations provided in this section, which refers to HF in Fig. 1A but NF in Fig. 1B. Substituting pseudo-pressure (Eq. (2)) and the dimensionless group provided in Appendix A in Eq. (12) leads to: 2m D yD2

)2

(n ) is calculated using Eq. (6) for pressure value of i-th fracture i element at n-th time-step. wf is the width of the finite-conductivity fracture. The relationship between two plane source terms QmD and QfD can be straightforwardly obtained by combining Eq. (11) and Eq. (16):

f g)

t

(n ) (n ) 2mfD , i + mfD, i

where

The governing diffusivity equation of 1D gas flow within a finiteconductivity fracture along y axis shown in Fig. 3 is written as:

kf p + Qf = µg y

(15)

i = 1, 2, 3, ...,N

5. Fracture linear flow solution

g

m cti

f km cfti

( yD

(10)

Qs(,ji) is the strength of the line-source parallel to the x-axis in mass per unit area per unit time. In this study, it is assumed to be a stepwise function, which is constant for a given fracture element i and time interval j. Eq. (9) provides a system of N equations where unknowns are mfD,i and QmD,i of each element—a total of 2 N unknowns. The simultaneous solution of the set of 2 N unknowns, corresponding to conditions along the N fracture elements, is achieved by coupling the integral solution for matrix flow (Eq. (9)) with fracture flow model discussed in the subsequent section. Notably, under the assumption of 1D Cartesian flow within fracture domain, the discretized fracture domain is a one-dimensional system along the y- direction (Fig. 3). Therefore, the strength of the line plane source Qs is equal to the summation of normal flux densities of fluid flow towards both sided fracture element, as proven analytically by Pecher [23].

y

kf

(n ) mfD ,i+1

Qs(,ji) RTµgi Zi Lf MWpi km mi

is the dimensionless fracture storage capacity [8]:

Assuming a uniform discretization of fracture domain as well as piecewise constant of nonlinearity during every time interval, applying finite-difference scheme to Eq. (15) results in:

(9)

(j ) QmD , i is the dimensionless source strength of i-th fracture element at j-th time step: (j ) QmD ,i =

=

(14)

wf MWpi kf mi

D

where

mfD(n, i)

Qs RTµgi Zi Lf2

QfD =

(13)

where is the pressure-dependent viscosity-compressibility ratio defined in Eq. (6). QfD is the source strength in dimensionless form:

5

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

6. Benchmark case: Validation against classical analytical solution

term of dimensionless variables given in Eq. (15) and (19) for two scenarios, for which Cinco-Ley and Samaniego [9] provided closedform, analytical solutions: one assumes steady state flow within fracture domain while neglecting fracture storage capacity, and the other accounts for fracture storage capacity by including the RHS of Eq. (12). Comparison results are displayed in Figs. 4 and 5 for with and without fracture storage effect, respectively. Bilinear, transition and late-time linear flow signatures are identified from ¼ slope and ½ slope in Fig. 4 where fracture steady-state flow is assumed. Fracture linear flow period (early-time ½ slope) is observed in Fig. 5 where storage capacity is considered. In both figures, excellent matches between the proposed solution and Cinco-ley and Samaniego [9] inverse-Laplace solution verify the validity of our solution when applies to gas bilinear flow under the assumption used in classical well test analysis of gas wells that the nonlinearity of λ is neglected.

The proposed solution is first benchmarked against the analytical solution developed by Cinco-Ley and Samaniego [9] for single-phase vertical oil wells producing from finite conductivity hydraulic fracture in an infinitely large domain. Under the same assumptions of one-dimensional flow in both fracture and formation as this study, Cinco-Ley and Samaniego [9] solved liquid flow diffusivity equations in both fracture and reservoir domain using Laplace transform. In the classical well test analysis for conventional gas reservoirs, and the remainder nonlinear viscosity-compressibility ratio λ (Eq. (6)) in the RHS of the pseudopressure-based gas diffusivity equation is assumed to be constant throughout the test (pressure and thus time-independent) and evaluated at initial condition (λ = 1). Therefore, the classical method to apply Cinco-Ley and Samaniego’s liquid solution to compressible, gas flow, the following dimensionless functions are employed, which essentially replaces pressure in the liquid solution by real gas pseudopressure [11]: i

pD =

7. Application to multi-fractured horizontal gas wells: validation against numerical simulations In this section, synthetic cases with real reservoir fluid properties are used to test the capability of proposed semi-analytical solution to rigorously capture the transient behavior of unconventional gas wells under the effect of pressure-dependent, nonlinear term λ. For validation purpose, results obtained by using the proposed method are compared with simulation data generated using a finely-gridded, finite-element commercial simulator (COMSOL Multiphysics®). In both proposed solution and COMSOL, gas fluid properties are calculated using Dranchuck and Abou-Kassem [10] for Z-factor, Lee et al. [19] for gas viscosity, Abou-Kassem et al. [1] for gas isothermal compressibility, and Sutton [28] for pseudo-critical property calculations.

wf

(20)

i qD

qD is the dimensionless flow rate: qD =

2psc Tqsc (21)

Tsc h i km

where

=

0

p 2p dp µg Z

=

2pi µgi Zi

m.

This extension of using Eqs. (20) and (21) in gas well analysis is approximate because the nonlinearity of the retaining viscosity-compressibility ratio λ in Eq. (3) and (13) has been neglected. Neglecting this nonlinearity term is generally considered sufficiently accurate for the analysis of vertical gas wells in conventional reservoirs, which is Cinco-Ley and Samaniego [9] studied, but it is shown to lead to erroneous predictions when applied to unconventional MFHWs as discussed in subsequent case studies. To generate the results using our proposed semi-analytical solution following the same simplification/approximation, λ is taken as constant and equal to unity in both Eq. (3) and (13). As a result, the adjustment factor fa is simply dropped in Eq. (9). Comparisons are conducted in

7.1. Application #1: Horizontal gas wells with finite-conductivity hydraulic fractures Case 1 investigates the transient pressure responses of a horizontal gas well with finite-conductivity bi-wing hydraulic fractures and singleporosity matrix (Fig. 1A). It is a tight gas sand with high reservoir pressure and temperature conditions (pi = 16,200 psia and T = 760 R), which was originally reported by Rushing et al. [25] and later studied in Zhang and Ayala [34] under the assumption of infinite-conductivity

10000

Dimensionless Pressure PD

1000 100 10 1

From top to bottom CfD = 0.1 1 10 100 1000

1/2 slope

1/4 slope

0.1

Cinco-Ley and Samaniego-V (1981)

0.01

Proposed Solution with λ = 1 0.001 1.E-6

1.E-5

1.E-4

1.E-3

1.E-2

1.E-1

1.E+0

1.E+1

1.E+2

1.E+3

1.E+4

1.E+5

Dimensionless Time tD Fig. 4. Validation against Cinco-Ley and Samaniego [9], without fracture storage effect (steady-state flow).

6

1.E+6

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

10000

From top to bottom ηD = 1000 100 20

Dimensionless Pressure PD

1000

1/2 slope

100

10

1/4 slope 1

Cinco-Ley and Samaniego-V (1981)

0.1

Proposed Solution with λ = 1

1/2 slope 0.01 1.E-6

1.E-5

1.E-4

1.E-3

1.E-2

1.E-1

1.E+0

1.E+1

1.E+2

1.E+3

1.E+4

1.E+5

1.E+6

Dimensionless Time t D Fig. 5. Validation against Cinco-Ley and Samaniego [9], with fracture storage effect and CfD = 1.

7.1.1. Bottomhole pressure response of constant rate specification Fig. 2 compares the well bottomhole pressure versus time predicted by our proposed semi-analytical solution versus COMSOL’s numerically simulated results. Results under three constant-rate specifications are the displayed Fig. 6 considering following specifications: qsp = 1, 2 and 3 MMSCFD. This specified rate is the constraint at the HW-HF intersection (Fig. 2), which is treated as a boundary condition in this study. Solid lines represent results generating using COMSOL, and proposed solution is shown as dot line. For comparison purposes, dash lines are the predictions provided by classical solution that does not consider the pressure-dependency of λ term by simply assuming it equal to one. Excellent matches can be observed between straight and dot lines, indicating that that the proposed solution successfully predict the transient pressure behavior of gas wells with finite hydraulic conductivity by rigorously capturing nonlinearity effect due to pressure-dependent λ term. The impact of nonlinear λ term on BHP responses are also clearly observed. The difference between proposed solution and classical

Table 1 Reservoir and fluid properties, single-porosity matrix case. Initial pressure, pi Temperature, T Specific gravity, SG Matrix permeability, km Matrix porosity, m Thickness, h Hydraulic fracture half length, LF Hydraulic fracture permeability, kF Hydraulic fracture porosity, φF Hydraulic fracture width, wF Pore compressibility, cφ

16,200 psia 300 F 0.633 9 µd 0.0661 200 ft 300 ft 400 md 0.3 0.01 ft 1 × 10−6 psi−1

HFs. Input reservoir and fluid properties are displayed in Table 1. Comparisons are conducted in terms of: 1) bottomhole pressure (BHP) responses; 2) pressure distributions along hydraulic fracture. Both constant rate and variable rate production constraints are tested in this case. 18000 16000

Bottomhole pressure (psi)

14000

1 Mmscf/d

12000 10000

2 Mmscf/d

8000 6000 4000

COMSOL Proposed Solution

2000

3 Mmscf/d

Classical solution 0

0

50

100

150

200

250

300

Time (days) Fig. 6. Bottomhole pressure response for three different flow rate specifications.

7

350

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

Rate-normalized pseudopressure(psi/(scf/d))

0.01

1/2 slope

1/4 slope

0.001

0.0001 0.01

0.1

1

10

100

1000

Time (days) Fig. 7. Pseudo-pressure normalized rate vs time solved by proposed solution (qsp = 1MMSCFD).

16000 Bottomhole pressure (psi)

4

COMSOL

14000

Proposed Solution

3.5

Classical Solution

3

12000

2.5

10000

2

8000

1.5

6000 4000

1

2000

0.5

0

Specified flow rate (MMscf/d)

18000

0

0

50

100

150 200 Time (days)

250

300

350

Fig. 8. Bottomhole pressure response under variable-rate specification.

solution is relatively small at low rate. Increased rate will amplify the difference, and the classical solution fails to provide reasonable prediction after 250 days at the highest flow rate (3 MMSCFD) specification. This observation is expected as classical solution assumes the viscosity-compressibility product remains at initial condition (λ = 1).

In other words, the more the BHP drops from initial pressure, i.e., the higher the initial drawdown, the less accurate the classical solution becomes. These comparisons demonstrate that the pressure-dependent viscosity and compressibility product can have a very significant impact in transient behavior of multi-fractured gas wells, especially when the

18000 16000

5 days

Pressure (psi)

14000 55 days

12000

100 days

10000 8000

305 days

6000 4000

COMSOL

2000

Proposed Solution

0 0

50

100 150 200 Distance from wellbore, y (ft)

250

300

Fig. 9. Pressure distribution along fracture for constant-rate specification qsp = 3 MMSCFD.

8

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

18000 16000

5 days

14000 55 days

Pressure (psi)

12000

100 days

10000

305 days

8000 6000

4000 COMSOL

2000

Proposed Solution

0 0

50

100 150 200 Distance from wellbore, y (ft)

250

300

Fig. 10. Pressure distribution along fracture for variable-specification.

wells are experiencing significant pressure drawdown at higher rate specifications. Same observations have been made by Zhang and Ayala [34] for production from infinite-conductivity HFs in infinite acting reservoir, and have been now further confirmed for the case of production from finite conductivity HFs. Fig. 7 displays pseudopressure normalized rate versus time for the constant production rate of 1 MMSCFD. The rate-normalized pressure is (m m ) calculated as: i wf . This plot shows clear bilinear flow regime (one-

7.2. Application #2: Multi-fractured horizontal gas wells with naturallyfractured matrix This section presents the examples of proposed solution applied to a naturally-fractured matrix under the context of a triple porosity model. There is a total of three (sequential) flow regions (as displayed in Fig. 1B): flow from matrix to natural fracture, from natural fractures to hydraulic fracture, and from hydraulic fracture to wellbore. To apply the proposed semi-analytical solution to the analysis, the hydraulic fracture is assumed to have infinite conductivity: pressure drops along HF is negligible and pressure in HF is thus uniform and identical to well sand-face pressure. The proposed semi-analytical method, which is has been successfully applied to finite-conductivity HF and single-porosity matrix model in Application #1 in previous section, is thus applied to solve the coupled 1D flow from matrix to natural fracture and 1D flow from natural fracture to hydraulic fracture under a triple porosity model. This case study is based on a synthetic gas shale case investigated by Kanfar and Clarkson [16] in which they studied the horizontal well production performance under constant-pressure specification using numerical simulation and an empirical method. Input reservoir and fluid properties are summarized Table 2. For the purpose of this paper, transient pressure behavior of horizontal gas wells under constant- and variable- rate specifications are investigated. The specified rate corresponds to their value at each intersection location of infinite-conductivity HF and finite-conductivity NF shown in Fig. 2. Under the assumption of 1) both hydraulic and natural fractures are identical and uniformly-spaced, and 2) no significant pressure losses along the horizontal wellbore and HF, multiplying the specified rate with the number

qsp

fourth slope) from 0.1 day up to a few days. After a few hundred days, the late-time matrix liner flow regime dominates (half slope). 7.1.2. Bottomhole pressure response of variable rate specification Fig. 8 presents the well BHP comparisons under variable flow rate constraint between numerical simulation results and the proposed solution. The well is scheduled to start producing at high rate value of 3 MMSCFD, which then drops to 1.5 MMSCD after 150 days and finally shuts-in after 300 days. Throughout the 365 days of production, BHP responses predicted by proposed solution (dot line) matches perfectly with simulation data of COMSOL (solid line). The good agreement between proposed solution and simulated results verifies the validity of proposed semi-analytical approach in solving the coupled nonlinear gas diffusivity equations of matrix and fracture domain under variable-rate production constraint. In addition, predictions given by classical solution, which wrongful assumes viscosity-compressibility at initial condition during production, shows a very clear disparity from simulation and proposed solution. 7.1.3. Pressure distribution along hydraulic fracture Figs. 9 and 10 present the comparisons of the pressure distributions along the fracture solved by proposed solution and COMSOL for constant- and variable- rate, respectively. Four different time steps are selected: 5, 55, 100 and 305 days. In both figures, good consistency is observed between the proposed solution and simulated data generated using by COMSOL for all time steps tested, suggesting that the proposed semi-analytical method can accurately capture flow dynamics within fracture domain. Notably, the numerical model built in COMSOL for all cases presented in this paper follows the schematic shown in Fig. 2 in which matrix flow is modelled under 2D Cartesian geometry. The good matches in term of fracture pressure distributions between propose solution in which 1D matrix flow is assumed and COMSOL model that uses rigorous 2D flow model confirms again the applicability of the use of 1D matrix flow model in the proposed coupled solution.

Table 2 Reservoir and fluid properties, naturally-fractured matrix case. Initial reservoir pressure, pi Temperature, T Gas gravity, SG Matrix permeability, km Matrix porosity, φm Hydraulic fracture permeability, kF Natural fracture permeability, kf Natural fracture porosity, φf Hydraulic fracture spacing, 2Lf Natural fracture width, wf Reservoir thickness, h Pore compressibility, cφ

9

4000 psia 200 F 0.7 10 nd 0.07 Infinite 0.05 md to 0.5 md 1 500 ft 0.01 ft 300 ft 1 × 10−6 psi−1

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

4000

HF pressure (psi)

3500 3000 2500 2000 1500

From top to bottom kF = 0.5 md 0.1 md 0.05 md

1000

500

COMSOL

Proposed Solution Classical Solution

0 0

200

400

600

800

1000 1200 Time (days)

1400

1600

1800

2000

Fig. 11. HF Pressure responses under different natural fracture permeability.

4500

4500 COMSOL Proposed Solution

3500 HF Pressure (psi)

4000 3500

Classical Solution

3000

3000

2500

2500

2000

2000

1500

1500

1000

1000

500

500

0 0

200

400

600

800

1000

1200

1400

1600

1800

Specified rate (scf/d)

4000

0 2000

Time (days) Fig. 12. HF pressure response under variable rate conditions.

4500

25 days

4000 3500

Pressure (psi)

3000 1000 days

2500

500 days

2000 1500 1000

COMSOL 500

Proposed Solution

2000 days 0 0

20

40

60

80

100

120

140

160

180

Distance from HF, y (ft) Fig. 13. Pressure distribution along natural fracture for constant rate case with kf = 0.05md. 10

200

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

4500 50 days 4000

Pressure (psi)

3500 3000

1000 days

2500 1450 days 2000 1500

200 days COMSOL

1000

Proposed Solution

500 0

20

40

60

80

100

120

140

160

180

200

Distance from HF, y (ft) Fig. 14. Pressure distribution along natural fracture for variable-rate case with kf = 0.1md.

of HF-NF intersections equals the total well production rate. Pressure within the HF is assumed the same as that of the horizontal well.

7.2.3. Pressure distribution along natural fracture Pressure distributions along the natural fractures at different time steps are displayed in Figs. 13 and 14 for constant- and variable-rate cases, respectively. Compared to the hydraulic fractures pressure distributions shown in Figs. 9 and 10, pressure transient within natural fractures is confined to a much smaller area because of the significantly lower permeability in natural fractures (kf ≪ kF), which also leads to a higher pressure drawdown around the intersection of HF-NF networks. Despite the more significant pressure drawdown, the proposed solution is hereby shown to successfully predict the entire pressure profiles for all scenario tested in both figures by matching with COMSOL results. Notably, given the steeper pressure changes in near HF-NF intersection area, more refined grids in fracture 1D flow model are required to obtain the following match results.

7.2.1. Bottomhole pressure response of constant rate specification Fig. 11 shows the comparison between numerical simulation and the proposed solution for the pressure of infinite-conductivity hydraulic fractures producing at the same flow rate of 2.3 MSCFD but for different natural fracture permeability kf of 0.05, 0.1 and 0.5md. As expected, pressure of HF drops faster under higher natural fracture permeabilities. The predictions given by proposed solution continues to show good consistency with numerically simulated results (COMSOL) for all case tested, further validating the applicability and reliability of the proposed solution of analyzing the transient behavior of multi-fractured horizontal gas wells with naturally-fractured matrix. The classical solution significantly overestimates the pressure drop for all cases, with the worst offender being the scenario that experiences the most dramatic pressure decline (kf = 0.05 md) because of the inherent assumption of viscosity-compressibility remaining at initial condition throughout the production (λ = 1).

8. Concluding remarks This work proposes a straightforward semi-analytical solution for transient pressure behavior of multi-fractured horizontal gas wells producing from finite-conductivity fractures. The proposed solution is developed based on coupling the solutions of two contiguous flow domains—matrix and fracture—both of which exhibit one-dimensional Cartesian flow. The gas flow equation within matrix is solved using Green and source function method, which was proposed in our previous work for infinite-conductivity fractures [34]. Gas flow within finiteconductivity fracture is modelled using finite-difference scheme. The two domain solutions are later coupled to solve for pressure and fluxes along fracture simultaneously. We present the application examples of proposed solution under two different interpretations: 1) finite-conductivity hydraulic fractures in a single-porosity matrix; and 2) infiniteconductivity hydraulic fractures in a naturally-fractured matrix. Validations are provided by comparing against commercial numerical simulator COMSOL for production scenarios of both constant and variable-rate specifications in both application examples. Excellent matches are found in all case tested. One of the most important contribution and novelty of the present work is the effect of retaining nonlinearity of gas viscosity-compressibility product in transient pressure behavior of gas horizontal wells

7.2.2. HF pressure response under variable rate specification Fig. 12 demonstrate the HF pressure responses when the horizontal gas wells are subject to a variable-rate production schedule. Natural fracture permeability kf is selected to be is 0.1 md for this example. Specified rates are also shown in Fig. 12, which considers a more severe production scenario when a high rate is imposed starting from the onset of production. Under such condition, classical pseudo-pressure solution fails to provide reliable predictions due to neglecting the significant effect that pressure drawdown has on fluid properties. This predictive failure particularly worsens when the wells experiences significant drawdown, as evidenced by the attainment of negative pressures after around 700 days. The proposed semi-analytical solution, on the other hand, shows excellent matches with simulation data reported by COMSOL. This observation once again confirms the applicability of the proposed method to the analysis of transient data of unconventional gas wells using triple porosity models that rigorously capture the nonlinear, pressure-dependent, gas viscosity and compressibility effects.

11

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

with finite-conductivity fractures is captured in a rigorous, direct manner by coupling matrix and fracture solutions. It is also shown that neglecting this nonlinearity effect—which is a widely-used approximation in classical pressure transient analysis tools for gas wells producing from finite-conductivity fractures could lead to substantial prediction errors. In this paper, we focus on rate-specified well pro-

duction. However, the general finite difference form of the fracture flow solution should be applicable to any production constraint, such as bottomhole-pressure-specified. In addition, the solution can also be further extended to trilinear model in which both HF and NF have finite conductivities. These additional topics will be explored in our future work.

Appendix A Dimensionless variable group Table A-1 below provides the definitions of the dimensionless variables used in this manuscript. The term Lf refers to hydraulic fracture half length in single-porosity matrix case and half of hydraulic fracture spacing in naturally-fractured matrix case. Table A1 Dimensionless variables group. tD

km t 2 m µgi cti Lf

xD yD mD

x/Lf y/Lf

QmD

Qs RTµgi Zi Lf

mi m mi

MWpi km mi

Appendix B Effect of fracture storage term Fracture storage effects are sometimes insignificant in a fracture flow model given the extremely small fracture volume compared to the contribution of matrix. Under such approximation, the storage/accumulation term in the fracture flow equation (Eq. (12)) is dropped, resulting in:

y

g

kf p + Qf µg y

0

(B-1)

In order to verify this steady-state approximation, the original proposed solution is compared against the solution using Eq. (B-1) as fracture flow model. Fig. B-1 shows the bottomhole pressure responses for two case of Application #1: constant rate of 3 MMSCFD and variable rate example. The perfect match between the results of with and without storage term confirms the validity of the assumption of fracture storage effect being negligible in the case tested in this paper.

Fig. B-1. Proposed solution: with and without fracture storage effect.

Appendix C Flow chart of proposed semi-analytical solution Fig. C-1 below showcase the detailed flow chart of employing proposed semi-analytical method to solve coupled matrix and fracture equations simultaneously.

12

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

Start

Time step n =1

(n) (n) Initial guess of QmD ,i , m fD ,i

t D( n +1) = t D( n ) + ∆t D( n ) ; n = n + 1 (n) ∗( n ) ( n) Rmi = m fDi − (m fD , i + fa , i )

Matrix Eq: Fracture Eq:

R if =

 ∂ R m ,1 0 0  ( n)  ∂ m fD ,1  ∂ R m ,2  0 0  ∂ m (fDn),2  0  0   0 0 0    ∂ R f ,1 ∂ R f ,1 0  ∂ m ( n) ∂ m ( n) fD ,1 fD ,2   ∂R ∂ R f ,2 f ,2   ∂ m (fDn),1 ∂ m (fDn) ,2   0   0 0 

∗( n ) (n) Update mf D ,i , f a ,i using Eq.6 &7

Update λi( n ) using Eq.4

(n) (n) (n) QmD ,i = QmD ,i + ∆QmD ,i

m

(n) fD ,i

=m

(n) fD ,i

+ ∆m

m

( n) fD ,i +1

− 2m

( n) fD ,i 2

+m

(Eqs. 7, 9 &10)

(n) fD ,i −1

( ∆yD )

∂ R m ,1

0

0

( n) ∂ Q mD ,1

∂ R m ,2

0

0

0

0

0

0

0

∂ R m ,N ∂ m (fDn),N

0

( n) ∂ Q mD ,1

0

∂ R f ,N ∂ m (fDn),N

N

( n) ∂ Q mD ,2

∂ R f ,1

0

n

−

∂ R f ,2

0

( n) ∂ Q mD ,2

0

0

0

0

( ) QmD ,i

C fD

0 0

−

n

(n) i

ηD λ

0 0 0

0 0 0

∂ R m ,N ( n) ∂ Q mD ,N

0 0 0

0

n −1

m(fD),i − m(fD ,i )

1

∂R f , N ( n) ∂ Q mD ,N

∆t

(n) D

(Eq. 16)

       ∆ m ( n)  fD ,1     ∆ m (fDn),2       ∆ m f( Dn),N   =− ( n)    ∆ Q mD ,1   ∆ Q ( n)  mD ,2      ( n)    ∆ Q mD ,N     

 R m ,1    R m ,2      R m, N    R f ,1  R f , 2      R   f, N 

R < 1−7

(n) fD ,i

Y (n) ( n) S tore Qm D, i , mfD,i

t D( n ) < t D ,max

N

End Fig. C-1. Detailed flow chart of proposed semi-analytical solution.

References

[7] Brown M, Ozkan E, Raghavan R, Kazemi H. Practical solutions for pressure-transient responses of fractured horizontal wells in unconventional shale reservoirs. SPE Reservoir Eval Eng 2011;14(06):663–76. [8] Cinco-Ley H, Samaniego VF, Dominguez AN. Transient pressure behavior for a well with a finite-conductivity vertical fracture. SPE J 1978;18(04):253–64. [9] Cinco-Ley H, Samaniego VF. Transient pressure analysis for fractured wells. J Petrol Technol 1981;33(09):1749–66. [10] Dranchuk PM, Abou-Kassem H. Calculation of Z factors for natural gases using equations of state. J Can Pet Technol 1975;14(3). [11] Energy Resources and Conservation Board (ERCB). 1975. Theory and Practice of the Testing of Gas Wells, Canada. [12] Ezulike DO, Dehghanpour H. Characterizing Tight Oil Reservoirs Using Dual- and Triple-Porosity Models. Society of Petroleum Engineers, 2013. [13] Fraim ML, Wattenbarger RA. Gas reservoir decline-curve analysis using type curves with real gas pseudopressure and normalized time. SPE Form Eval 1987;2(4):671–82. [14] Gringarten AC, Bourdet DP, Landel PA, et al. A Comparison Between Different Skin

[1] Abou-Kassem JH, Mattar L, Dranchuk PM. Computer calculations of compressibility of natural gas. J Can Pet Technol 1990;29(5). [2] Al-Hussainy R, Ramey HJ, Crawford PB. the flow of real gases through porous media. J Petrol Technol 1966;18:624–36. [3] Al-Ahmadi HA, Wattenbarger RA. Triple-porosity Models: One Further Steps Towards Capturing Fractured Reservoir Heterogeneity. SPE 149054 presented at the SPE/DGS Saudi Arabia Section Technical Symposium and Exhibition held in AlKhobar, Saudi Arabia, 15-18 May 2011. [4] Barreto AB, Pires AP, Peres AM. A New Rigorous Analytical Solution for a Vertical Fractured Well in Gas Reservoirs. Paper SPE 150663 presented at SPE Latin America and Caribbean Petroleum Engineering Conference, 16-18 April, Mexico City, Mexico, 2012a. [5] Barreto AB, Peres A, Pires AP. A variable-rate solution to the nonlinear diffusivity gas equation by use of green’s-function method. SPE J 2012;18(01):57–68.

13

Fuel 266 (2020) 117119

K.H. Nguyen, et al.

[15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

[26]

and Wellbore Storage Type-Curves for Early-Time Transient Analysis. Presented at the SPE Annual Technical Conference and Exhibition, Las Vegas, Nevada, 23-26 September 1979. SPE-8205-MS, 1979. Jia P, Cheng L, Huang S, Liu H. Transient behavior of complex fracture networks. J Petrol Sci Eng 2015;132:1–17. Kanfar MS, Clarkson CR. Rate dependence of bilinear flow in unconventional gas reservoirs. SPE Reservoir Eval Eng 2018;21(01):17–30. Kazemi H, Merrill LS, Porterfield KL, Zeman PR. Numerical simulation of water-oil flow in naturally fractured reservoirs. SPE J 1976;16(6):317–26. Ibrahim MH, Wattenbarger RA. Analysis of rate dependence in transient linear flow in tight gas wells. SPE paper 100836 presented at Abu Dhabi International Petroleum Exhibition and Conference, 5-8 November, Abu Dhabi, UAE. 2006. Lee AL, Gonzalez MH, Eakin BE. The viscosity of natural gases. J Petrol Technol 1966;18(8):997–1000. Miranda FA, Barreto AB, Peres AMM. A novel uniform-flux solution based on the green’s function method for modeling the pressure-transient behavior of a restricted-entry well in anisotropic gas reservoirs. SPE J 2016;21(05):1870–82. Nobakht M, Clarkson CR. A new analytical method for analyzing linear flow in tight/shale gas reservoirs: constant-rate boundary condition. SPE Res Eval Eng 2012;15(1):51–9. Palacio JC, Blasingame TA. Decline-Curve Analysis with Type Curves - Analysis of Gas Well Production Data. SPE 25909 presented at the Low Permeability Reservoirs Symposium, 26-28 April, Denver, Colorado, USA, 1993. Pecher R. Boundary element simulation of petroleum reservoirs with hydraulically fractured wells. PhD dissertation. University of Calgary. 1999. Qanbari F, Clarkson CR. A new method for production data analysis of tight and shale gas reservoirs during transient linear flow period. J Nat Gas Sci Eng 2013;14:55–65. Rushing JA, Perego AD, Sullivan RB, Blasingame TA. Estimating Reserves in Tight Gas Sands at HP/HT Reservoir Conditions: Use and Misuse of an Arps Decline Curve Methodology. SPE paper 109625 presenetd at the 2007 SPE Annual Technical Conference and Exhibition held in Anaheim, California, U.S.A 11-14 November, 2007. Sousa EPS, Barreto Jr. AB, Peres AMM. Finite-wellbore-radius solution for gas wells by Green’ s functions. SPE J 2015;20(04):842–55.

[27] Sousa EPS, Barreto AB, Peres AMM. Analytical treatment of pressure-transient solutions for gas wells with wellbore storage and skin effects by the green’s functions method. SPE J 2016;21(05):1858–69. [28] Sutton RP. Compressibility factors for high-molecular-weight reservoir gases. SPE paper 14265 presented in the 60th SPE Annual Technical Conference and Exhibition. Las Vegas, Nevada 1985. [29] Tivayanonda V, Apiwathanasorn S, Ehlig-Economides C, Wattenbarger RA. Alternative Interpretations of Shale Gas/Oil Rate Behavior Using a Triple Porosity Model. SPE 159703 presented at the SPE Annual Technical Conference and Exhibition held in San Antonio, Texas, USA, 8-10 October 2012. [30] Wang C, Xiong Y, Huang Z, Winterfeld P, Ding D, Wu Y-S. A Multi-Porosity, MultiPhysics Model to Simulate Fluid Flow in Unconventional Reservoirs. Paper SPE 182698 presented at SPE Reservoir Simulation Conference, 20-22 February, Montgomery, Texas, USA, 2017. [31] Warren JE, Root PJ. The behavior of naturally fractured reservoirs. SPE J 1963;3(3):245–55. [32] Zhang M, Vardcharragosad P, Ayala H, L. F., The similarity theory applied to earlytransient gas flow analysis in unconventional reservoirs, J Nat Gas Sc Eng, 21, 659–668, 2014. [33] Zhang M, Ayala H, L. F., Gas-rate forecasting in boundary-dominated flow: constant-bottomhole-pressure decline analysis by use of rescaled exponential models. SPE J., 19(3), 410–417, SPE-168217-PA, 2014. [34] Zhang M, Ayala H. Variable rate and pressure integral solutions to the nonlinear gas diffusivity equation in unconventional systems. Fuel 2018;235:1100–13. [35] Zhang M, Ayala H. Application of superposition principle to variable rate/pressure production analysis of multi-fractured horizontal wells in unconventional gas reservoirs. J Natural Gas Sci Eng 2019;72:103011. [36] Vicente BJ, Pires AP, Peres AMM. A Semi-Analytical Solution for a Buildup Test for a Horizontal Well in an Anisotropic Gas Reservoir. Integral Methods in Science and Engineering, Volume 2. 1ed.: Springer International Publishing, 2017, v. 2, p. 285298. [37] Rocha AC, Murad MA, Vicente, Bruno J, Pires AP, Garcia ELM. A new model for flow in shale-gas reservoirs including natural and hydraulic fractures – application to well tests. In: 15th European Conference on the Mathematics of Oil Recovery (ECMOR XV), 2016, Amsterdam.

14