Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx 1

Contents lists available at ScienceDirect

Journal of Unconventional Oil and Gas Resources journal homepage: www.elsevier.com/locate/juogr 5 6

4

Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method

7

Pichit Vardcharragosad ⇑, Luis Felipe Ayala 1

8

The Pennsylvania State University, University Park, PA, United States

3

9 1 4 1 2 12 13 14 15 16 17 18 19 20 21 22 23

a r t i c l e

i n f o

Article history: Received 5 August 2014 Accepted 18 May 2015 Available online xxxx Keywords: Gas reservoirs Significant transient flow A density-based method Rate-time prediction Tight gas Shale gas

a b s t r a c t Early transient flow corresponds to the period before the effect fluid depletion has reached the nearest reservoir no-flow boundary. Production from unconventional reservoirs tends to exhibit extended periods of early transient flow because of their low permeabilities. Massive flow areas are generated, typically through the creation of multiple fractures in horizontal wells, to feasibly produce hydrocarbons from these formations at economic rates. The presence of these fractures leads to a series of non-radial flow regimes, which may continuously change before reservoir no-flow boundaries are reached, with linear flow being one of the dominant regimes. One of the significant challenges in this area has been devising a proper production analysis technique applicable to the analysis of early transient flow data. Progress has been made in the area through the use of the concept of the region of influence, which accounts for the portion of reservoir volume responsible for early transient production. In this study, we propose to implement a density-based approach to analyze early transient production data. In the density-based approach, rate-time responses of gas reservoir system are predicted by rescaling the responses of liquid system with depletion driven variables. The density-based technique has previously proven applicable to boundary-dominated radial-flow, and has been extended to analyze boundary-dominated linear-flow behavior. In this work, we show that early transient flow behaviors can be analyzed using the  calcukb density-based method that incorporates region of influence concept into rescaling variables,  lations. A density-based procedure is proposed to analyze early transient production data and its applicability is verified using simulated rate-time data. Results show that the proposed method can effectively predict Contacted Gas In-Place and the fracture half-length and square root of permeability product. The density-based methodology provides an alternative and reliable means to model and analyze data from gas reservoirs exhibiting extended early transient production. Ó 2015 Elsevier Ltd. All rights reserved.

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

49 50

Introduction

51

Reliable rate-time forecasting for reservoirs in decline is critical in the evaluation the economic feasibility of any hydrocarbon asset. For unconventional assets exhibiting extended periods of early transient flow, classical early-time linear-flow solutions for constant bottomhole pressure production conditions have proven particularly useful in the analysis when deployed using a number of modifications or corrections. Although rarely expected during early-time production, the study of early-time radial-flow

52 53 54 55 56 57 58

⇑ Corresponding author at: PTTEP, Energy Complex Building A, 555/1 Vibhavadi Rangsit Road, Chatuchak, Bangkok 10900, Thailand. Tel.: +66 (0) 2537 4000, mobile: +66 (0) 8 1849 9780. E-mail addresses: [email protected], [email protected] (P. Vardcharragosad), [email protected] (L.F. Ayala). 1 103A Hosler Building, University Park, PA 16802, United States. Tel.: +1 814 8654053.

solutions can be useful in exploring the effect of flow geometry on these early-time rate-time predictions. These classical solutions originally derived for liquid systems for both linear and radial flow geometries are summarized below along with traditional methods and more recent developments aimed at extending their applicability to gas analysis.

59

Early-time liquid solutions – linear flow

65

The constant-pwf early transient flow solution for the classical linear flow system depicted in Fig. 1 (i.e., production from an infinite conductivity fracture extending to both ends of the drainage reservoir boundaries) is given by, in the Laplace domain and for a slightly compressible fluid (van Everdingen and Hurst, 1949).

66

60 61 62 63 64

67 68 69 70

71

D ¼ q

2 1 pffiffi p s

ð1Þ

http://dx.doi.org/10.1016/j.juogr.2015.05.003 2213-3976/Ó 2015 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

73

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 2

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

Nomenclature Roman A Ac Ar Bx cg ct Cr fc Gp h k kg MW p q r re r inf rw rq R s S t T v xe xf X Y

74 75

76 78

79 80

81 83

Z entire reservoir area, L2 , ft2 cross-sectional flow area around wellbore or hydraulic fractures, L2 , ft2 reservoir area inside ROI, L2 , ft2 formation volume factor of x, reservoir volume/volume at standard conditions gas compressibility, L  t 2 =M, psi1 total compressibility, L  t 2 =M, psi1 coefficient used to calculate radius of influence, dimensionless correction factor, dimensionless cumulative gas production, L3 , SCF reservoir thickness, L, ft permeability, L2 , md gas permeability (effective permeability with respect to vapor phase), L2 , md molecular weight, M/mole, lb/lbmol pressure, M=L  t 2 , psia volumetric flow rate, L3 =t, ft3/day radius/ distance from wellbore, L, ft external (outer reservoir) radius, L, ft radius of influence, L, ft wellbore (inner reservoir) radius, L, ft wellbore-to-initial density ratio, dimensionless universal gas constant, M  L2 =t 2 -T-mole, 10.73 ft3-psi/R-lbmol Laplace variable saturation, dimensionless time, t, day or hr temperature, T, R Darcy velocity, L=t, ft/day outer reservoir boundary in x-direction, L, ft fracture half-length, L, ft multiplier used to transform time to dimensionless time, 1=t, day1 multiplier used to transform rate to dimensionless rate, t=L3 , SCFD1

where s represents Laplace variable. Eq. (1) can also be written in the real time domain as (Wattenbarger et al., 1998):

qD ¼

2

p

1 pffiffiffiffiffiffiffiffi ptD

qosc Bo lo   a1 2pkh pi  pwf

Top View

Reservoir boundary Equipotenal Line

Greek unit conversion constant, 6:328  103 (ft2-cp)/ (psi-md-day) c k, dimensionless time-average c viscosity–compressibility ratio, dimensionless fluid viscosity, M=L  t, cp Boltzmann variable, dimensionless hydraulic diffusivity, L2 =t, ft2/day fluid density, M=L3 , lbm/ft3 porosity available to gas phase, dimensionless pseudopressure, M=L  t 3 , psi2/cp or normalized pseudopressure. M=L  t 2 , psi

a1  b k

l n

g q / w

Subscript D g i o sc wf

dimensionless variable gas initial conditions oil standard conditions (14.7 psia and 60 F) bottomhole flowing conditions

Superscript  density-based definition 0 normalized definition Abbreviation BDF boundary dominated flow BSCF billion standard cubic feet CGIP contacted gas in place GIP gas in place ODE ordinary differential equation OGIP original gas in place PDE partial differential equation ROI region of influence

tD ¼

a1 kt /lo ct x2f

ð4Þ

85

ð2Þ

For these equations, the classical liquid definitions of the dimensionless variables are used:

qD ¼

gas compressibility factor, dimensionless

ð3Þ

Top View

Stream Line Wellbore Infinite Conducvity Fracture

Fig. 1. Graphical representation of linear and radial flow systems.

Early-time liquid solution – radial flow

86

The constant-pwf early-transient solution for the classical radial flow system depicted in Fig. 1 is given by, in the Laplace domain and for a slightly compressible fluid (van Everdingen and Hurst, 1949):

87

pffiffi s K D ¼ pffiffi 1 pffiffi q s  K0 s

ð5Þ

where K 0 and K 1 are modified Bessel functions of the second kind, order 0 and 1, respectively. This solution can be inverted into real time domain using numerical method, e.g., the Stehfest algorithm (Stehfest, 1970). Eq. (5) can also be written in an infinite integral form in the real time domain (van Everdingen and Hurst, 1949). In addition, Dresner (1983) provided the approximate heat transfer solution for this scenario, which can be rewritten as the approximate solution of flow in porous media in the following form:



E1 4t1D



1   ffi qD 2  exp  1 4t D

88 89

90 92 93 94 95 96 97 98 99 100

101

ð6Þ

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

103

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx 104 105 106 107

Jacob and Lohman (1952), Earlougher (1977), Uraiet and Raghavan (1980), and Ehlig-Economides and Ramey (1981) used the following alternate real-time, line-source approximation for Eq. (5) for constant-pwf infinite-acting radial systems:

108





110

1 1 1 ffi  Ei  qD 2 4t D

111

where the logarithmic approximation was shown to be correct

ffi

1 ðln t D þ 0:80907Þ 2 11

ð7Þ

4

112

within 0.1% for t D > 5  10 , only 1% for t D > 8  10 , and 2% for

113

t D > 5  103 . The radial-flow definition of qD , remains the same as the one given in Eq. (3), while the radial-flow definition of t D is given by:

114 115

116 118 119 120 121 122 123 124 125 126

a1 kt tD ¼ /lo ct r2w

ð8Þ

In the context of unconventional production analysis, early transient radial solutions do not get much attention. Pure radial flow is not expected in these tight reservoirs typically completed with massive hydraulic fractures – especially at early times. We include these liquid solutions in this study for comparison purposes and to demonstrate the generality of the density-based technique and its applicability to early-time analysis for both flow geometries.

127

Early-time gas solutions: traditional and ‘corrected’ pseudo-functions

128

For gas systems, pseudo-functions have customarily been used to linearize the gas diffusivity equation. When the linearization is successful, the same qD liquid solutions above (Eqs. (2) and (7)) would still apply if the dimensionless definitions for qD and tD are re-expressed in term of pseudopressure ðwÞ (Al-Hussainy et al., 1966) and pseudotime ðta Þ (Agrawal, 1979; Fraim and Wattenbarger, 1987) functions, respectively, where assuming c t ffi cg :

129 130 131 132 133 134 135

136

q T psc  gsc  qD ¼ a1 pT sc kh wi  wwf

a1 kt a ðlinear casesÞ /lgi cgi x2f a1 kta ðradial casesÞ tD ¼ /lgi cgi r 2w

tD ¼

138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160

ðlinear and radial casesÞ

ð9Þ ð10Þ ð11Þ

pffiffi Based on observations of early-time qgsc performance vs. t (Eq. (2) for linear flow) in field data from unconventional gas reservoirs, Wattenbarger et al. (1998) suggested that tD could be defined as a function of actual time t during the early-transient flow period, with pseudotime ta only needed for boundary-dominated flow pffiffi analysis. The linear relationship between qgsc and t was then used by the authors to evaluate flow properties such as permeability or fracture half-length using production rate-time data. The estimated time when gas rate starts deviating from this trend (i.e., the so-called time to the end of the half-slope line) is also be used to estimate distance from the fracture to no-flow outer boundary and reserves. Anderson and Mattar (2007) first argued that the use of the traditional pseudotime ta definition may be ill-suited for production data analysis of reservoirs with extended periods of infinite acting flow. The traditional pseudotime ta , routinely used for boundary-dominated production data analysis based on reservoir average pressure over the entire reservoir drainage volume, fails to recognize that during early transient flow only a portion of (and not the entire) reservoir volume is responsible for production. Consequently, Anderson and Mattar (2007) implemented a modified version of the pseudotime function (i.e., a ‘corrected’

3

pseudo-time) based in the estimation of an average pressure inside a region of influence (ROI) during the calculation of the pseudotime integral. Ibrahim and Wattenbarger (2006) showed that the use of the pffiffi qgsc vs. t relationship in early-transient gas linear analysis, as originally proposed by Wattenbarger et al. (1998), suffered from a pressure drawdown dependency. They proposed an empirical, numerically-derived factor, f c , to correct for this drawdown dependency. Later, Nobakht and Clarkson (2012a) showed that this drawdown dependency could be rigorously calculated by evaluating t D for linear flow (Eq. (10)) as a function of the ‘corrected’ pseudotime variable proposed by Anderson and Mattar (2007). In their work, they derived a key insight: that average pressure inside the ROI remains constant throughout early transient flow period during gas linear flow under a constant-pwf production specification. As a result, the drawdown dependency previously observed by Ibrahim and Wattenbarger (2006) is shown as arising because of the use of actual time instead of ‘corrected’ pseudo-time in the pffiffi qgsc vs. t relationship. Because average pressure remains constant inside the ROI, they also analytically established a linear relationship between t and corrected ta variables through the square-root of the ratio between the lg ct product at initial and average pressure conditions. This ratio becomes the new correction factor ðf c Þ allowing users to analyze properties of gas reservoirs pffiffi using the convenient qgsc and t (actual time) relationship. More recently, Qanbari and Clarkson (2013, 2014) and Chen and Raghavan (2013) presented an alternate semi-analytical method for the prediction of gas rate-time responses during early transient linear flow without pseudotime or ROI calculations. Their work is based on a transformation of the original linear-flow partial differential diffusivity equation (PDE form) into a non-linear ordinary differential equation (ODE). This technique makes solutions for early-transient flow for infinite-acting conditions directly possible, for the cases considered by the authors, while sidestepping ROI and pseudo-time computations. Qanbari and Clarkson (2013) also presented a detailed procedure on how to solve the resulting non-linear ODE via integration and how to iteratively handle non-linearities. A description of those details was not presented in the Chen and Raghavan (2013) study. The Qanbari and Clarkson (2013) methodology directly solves the following linear-flow PDE:

  @ @wD 1 @wD ¼ @r D @rD gD @tD

ð12Þ

where:

wD ¼

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200

201 203 204

205

w0i  w0 w0i  w0wf

pffiffiffiffiffi rD ¼ r= Ac tD ¼ a1

161

gi

Ac gD ¼ g=gi

t

ð13Þ ð14Þ ð15Þ ð16Þ

207

and r represents distance from wellbore, Ac represents cross-sectional flow area around the wellbore or hydraulic fracture. In this methodology, the diffusivity equation is not linearized using the pseudotime function and Eq. (12) remains non-linear because of the pressure-dependent gD hydraulic diffusivity coefficient defined by Qanbari and Clarkson (2013) as:

208

gðpÞ ¼

a a 0 ¼ ð@b=@pÞ b

ð17Þ

where

k aðpÞ ¼ lg B g

209 210 211 212 213

214 216 217

218

ð18Þ

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

220

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 4

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

222

bðpÞ ¼ /=Bg

223

and where the following version of the pseudopressure function, normalized by the variable a evaluated at initial condition is used:

224

225

Z

2 w ðpÞ ¼ ai 0

227 228 229 230

231

ð19Þ

p

a  dp

ð20Þ

p0

Qanbari and Clarkson (2013) show that, for an infinitely large gas reservoir produced under constant pwf condition, a solution to the diffusivity PDE in Eq. (12) is given by

233

 R n ^  exp 2 0 gn d^n  dn D  wD ¼ 1  R R n ^n   1 exp 2 0 g d^n  dn 0

234

where n represents the well-known Botzmann variable, defined as:

Rn 0

235 237

ð21Þ

D

n¼

rD 2t 0:5 D

ð22Þ

247

For a gas system, Eq. (21) is solved iteratively because of the pressure-dependency built into the dimensionless hydraulic diffusivity ðgD Þ. The accuracy of the resulting semi-analytical solution is subject to the numerical integration scheme applied to evaluate the finite and infinite integrals, which is analogous to a grid discretization used by a numerical simulator. Once wD is determined, dimensionless flow rate ðqD Þ predictions can be derived by coupling Darcy’s flow equation with @wD =@n solved from Boltzmann’s transformed diffusivity equation and wD solution from Eq. (21). The qD solution is given by

250

pffiffiffiffiffiffiffiffi 1 ¼ 2pf c pt D qD

251

where dimensionless flow rate variable is defined as

254

qgsc lgi Bgi   qD ¼ 2pa1 kpffiffiffiffiffi Ac w0i  w0wf

255

and correction factor ðf c Þ is given by

238 239 240 241 242 243 244 245 246

248

252

ð23Þ

1

256 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283

pffiffiffiffi 

1 p @wD ¼ 2 fc @n

 ¼ n¼0

pffiffiffiffi

p

ð24Þ

Z 0

1

n

gD

dwD

ð25Þ

The iterative scheme required for Eq. (21) can be initiated by evaluating wD ðnÞ with gD of 1.0. For any subsequent iterations, gD ðnÞ is evaluated using the available wD ðnÞ profile and iterations continue until no significant improvement results from further trials. While Qanbari and Clarkson (2013) and Chen and Raghavan (2013) restricted their methodology to the analysis of infinite-acting linear-flow responses, Appendix A presents and verifies the extension of this method for the analysis of early-transient radial flow, derived by this study for comparison purposes. In this study, a density-based methodology is proposed and validated to analyze early transient production data. The technique was first used to analyze unsteady state flow of natural gas reservoirs without the use of pseudofunctions by Ye and Ayala (2012). In density-based analysis, the effect of reservoir depletion on fluid properties is accounted for through the use of dimensionless vari These variables are then used to rescale the solutions k  b. ables  readily available for slightly compressible fluid (oil) systems; and the rescaled solutions are shown to successfully model unsteady state flow behavior of natural gas systems. The density-based technique has proven applicable for boundary-dominated radial flow analysis under constant-rate and constant-flowing-pressure constraints (Ye and Ayala, 2012; Ayala and Ye, 2013a) and variable rate/flowing pressure production (Ayala and Zhang, 2013b). The analysis method has also been extended to analyze gas reservoir with significant rock compressibility (Zhang and Ayala, 2013),

systems with significant slippage and desorption effects (Vardcharragosad and Ayala, 2013), and systems with long-term linear flow production (Vardcharragosad and Ayala, 2014). In this work, the density-based approach is further extended to analyze unsteady state flow in gas reservoirs with significant early transient flow. The analysis is validated and demonstrated for both linear and radial flow geometries.

284

Region of influence (ROI)

291

The region of influence (ROI) concept, as originally proposed by Anderson and Mattar (2007), is based on the notion that two gas reservoirs with different sizes but identical reservoir and fluid properties and produced under the same production conditions should yield identical reservoir responses during early transient (or infinite acting) flow. A graphical demonstration of this concept is illustrated in Fig. 2. In the figure, reservoirs A and B, with identical reservoir and fluid properties but different sizes, are produced under the same production specification condition. At any time t during early transient flow, the ROI of reservoirs A and B represented by regions inside red-dash circles are identical, regardless their total available reservoir volume and actual location of their nearest reservoir no-flow boundary. Average pressure inside the highlighted ROI for these reservoirs will also be identical; but the average pressure within the entire available reservoir volumes (entire square) will be different. Consequently, the pseudotime function (Agrawal, 1979; Fraim and Wattenbarger, 1987) of reservoirs A and B at time t will be identical if based on average pressure inside the ROI, but traditional pseudotime based on average pressure within the entire drainage area will be different. The modified pseudotime calculation based on average pressure inside the ROI is called the ‘‘corrected pseudotime’’ function. For practical purposes, Anderson and Mattar (2007) suggest invoking the traditional material balance equation to evaluate average pressure inside the ROI domain. For a dry gas natural gas reservoir undergoing volumetric depletion, the average pressure inside ROI would be evaluated as:

292

   pi Gp p ¼ 1 CGIP Z Zi

ð26Þ

where CGIP (Contacted Gas In-Place) denotes amount of GIP (Gas In-Place) inside the ROI. For a reservoir with uniform thickness and porosity, the CGIP can be calculated from

CGIP ¼

Ar h/qi

qsc

ð27Þ

where Ar represents the reservoir area inside the ROI, which changes with time as ROI grows, and / represent the porosity available to gas phase. While contacted reservoir volume may be determined from the portion of reservoir which experiences pressure

Fig. 2. Graphical representation of the region of influence (ROI) in radial flow.

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

285 286 287 288 289 290

293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318

319 321 322 323 324

325 327 328 329 330 331

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx 332 333 334 335 336 337 338 339 340

341

343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375

depletion beyond a certain threshold level, for practical purposes the authors suggest the use of the traditional radius of investigation ðrinf , in this work) formulation for the determination of Ar and CGIP values. Because rinf would increases with time, the ROI will extent to the near-wellbore vicinity at very early time and expand as production continues. It will be bound by the location of reservoir no-flow boundaries during boundary dominated flow (BDF) period, at which point Ar ¼ A and CGIP = OGIP. For a radial geometry (Energy Resources Conservation Board, 1975; Lee, 1982):

r inf ¼ C r 

a1 kt /lgi cti

!0:5 ð28Þ

In this equation, using the traditional definition of radius of investigation for a radial flow geometry, C r ¼ 2:0 (Energy Resources Conservation Board, 1975; Lee, 1982). Nobakht and Clarkson (2012a) recommended, based on numerical simulation, the use of a C r of 2.554 to determine ROI of linear flow under a constant bottomhole flowing pressure production. More recently, Behmanesh et al. (2014) implemented two analytical-based methodologies, the impulse response and transient/ boundary-dominated solution intersection methods, for r inf calculations. For constant pwf production, the impulse method yields a C r of 2.44, while the transient-BDF intersection method results in a C r of 2.264. They concluded that the rinf calculated from the impulse response method is more suitable for early transient linear flow analysis than the classical r inf originally used for linear analysis by Wattenbarger et al. (1998) and applied by Anderson and Mattar (2007). In addition, Tabatabaie et al. (2013) proposed new procedure to determine region of influence for constant rate production by imposing a unit-slope at each point in time. Using this method, region of influence can be evaluate without the use of rinf calculation. The concept of ‘‘radius of investigation’’, a classical concept in well test analysis, has been presented in different ways by different researchers. In most cases, the radius of investigation is presented in the same functional form as shown in Eq. (28) with C r values ranging between 0.379 and 4.29 (see Table 1 in Kuchuk, 2009). In other cases, it is defined in entirely different terms (Kuchuk, 2009). ROI (as used in production data analysis) and ‘radius of investigation’ (as used in well test analysis) concepts are not necessarily identical, even when the latter is used to estimate the former. In the context of this study, the ROI focuses on the determination of the portion of reservoir volume that has been influenced and, in return, can influence wellbore production at a given time. It is also important to highlight the interplay between

i

ii

5

rinf (Eq. (28)) and average-density-within-ROI (Eq. (26)) calculations for the purposes of production data analysis. Depending on the selected value of C r used for rinf estimations in Eq. (28) above, two regions (regions i and ii) may be created within the depleted region highlighted in Fig. 3. Fig. 3 shows density distribution profile inside a hypothetical reservoir at time t. Regions i and ii represent amount of reservoir fluid which has been produced from the reservoir volume inside and outside rinf (= ROI, if r inf is used to define ROI), respectively. From this figure, we see that depending on the resulting location of r inf at time t (controlled by the C r value used in Eq. (28), if all other properties are assumed the same), there is no guarantee that all produced reservoir fluid at time t (i.e., Gp Þ would have been produced from region i alone when implementing Eq. (26). It follows that the ‘proper’ value of C r , for the purposes of the material balance equation in Eq. (26), is the one that eliminates the potential presence of region ii during pseudotime calculations. This may not be always accomplished for all possible scenarios given the large range of C r values available to the analyst, described above, and its stated dependency on flow geometry and production specification conditions.

376

Early transient rate-time responses: a density-based approach

396

In a density-based approach, rate-time responses of a gas reservoir system producing under a constant pwf constraint are pre kb dicted by rescaling the responses of liquid systems via 

397

dimensionless variables (Ye and Ayala, 2012). The rescaled solution is calculated as:

400

     bt qgas t D ¼ k  qliq D D D

381 382 383 384 385 386 387 388 389 390 391 392 393 394 395

398 399 401

402

405

qD

qsc li ci qgsc   ðlinear and radial casesÞ ¼ a1 2pkh qi  qwf

406

407

ð30Þ

tD ¼

a1 kt /li ci x2f

ðlinear casesÞ

ð31Þ

tD ¼

a1 kt /li ci r2w

ðradial casesÞ

ð32Þ

409

where the ‘⁄’ superscript denotes the use of a density-based version of the dimensionless variable. In the transformation in Eq. (29), qgas D

410

is the predicted dimensionless gas rate and qliq D is the dimensionless rate solution of the analogous liquid system. The viscosity– kÞ, originally proposed by Ye and Ayala compressibility ratio ð (2012) and later rigorously redefined by Ayala and Zhang (2013b) and Zhang and Ayala (2014) is given by the expression:

412

¼1 b t

Fig. 3. Graphical representation of density distribution in a gas reservoir.

380

where qD and t D are the density-based definitions of dimensionless rate and time, given by:

l c l c k ¼ gi gi ¼ gi gi  l g cg 2h q qwf

re

379

404

Z

t

k  dt

411

413 414 415 416

417

ð33Þ 419

where h ¼ RT=MW and the original pseudo-pressure definition of  is the time-averaged value of Al-Hussainy et al. (1966) is used. b k, thus defined as: the viscosity–compressibility ratio 

rinf

378

ð29Þ

wwwf

rw

377

ð34Þ

0

 ¼ 1 in Eq. (29), the liquid solution is readily By setting  k¼b  definitions are valid as long as ct ffi cg kb recovered. These  (Zhang and Ayala, 2013; Vardcharragosad and Ayala, 2013). The rescaled transformation in Eq. (29) has been proven readily applicable to boundary-dominated rate analysis in linear flow (Vardcharragosad and Ayala, 2014) and radial flow (Ye and Ayala, 2012) when  k is calculated in terms of the entire reservoir drainage

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

420 421 422

423 425 426 427 428 429 430 431 432

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 6

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

445

area. For BDF radial flow, Ye and Ayala (2012) rescaled Laplace liquid radial solutions and successfully matched gas rate responses from numerical simulators. For linear BDF, the generalized liquid solution used by Wattenbarger et al. (1998) was successfully rescaled to match gas rate responses from numerical simulators (Vardcharragosad and Ayala, 2014). The rescaling has proven to remain valid for systems with significant rock compressibility effects (Zhang and Ayala, 2013) and systems with significant slippage and desorption effects (Vardcharragosad and Ayala, 2013). In this section, we show that this density-based rescaled transformation in Eq. (29) remains fully valid for early transient flow rate analysis for both linear and radial systems as long as the k is calculated in terms the ROI (i.e., ROI-based  kÞ. rescaled variable 

446

Density-based early-transient solution: linear flow

447

b  transformation stated in Eq. (29), and On the basis of the k the corresponding liquid counterpart in Eq. (2), the density-based rate-time analytical solution for infinite-acting gas linear flow is  rescaled model: kb readily written as the following 

434 435 436 437 438 439 440 441 442 443 444

448 449 450

451

qD ¼ k  453

2

p

1 qffiffiffiffiffiffiffiffiffiffiffi   pbt

463 464

where the CGIP is calculated as:

457 458 459

461

465 467 468



CGIP ¼

 ð36Þ

Ar h/qi

ð37Þ

qsc

For linear flow, the contacted area defined by r inf , or Ar , is:

469 471

   Ar ¼ 2xf  2rinf

472

where

473

475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494



r inf

a1 kt ¼ Cr  /lgi cgi

Initial reservoir pressure (psia) Bottomhole flowing pressure (psia) Reservoir temperature (F) Gas specific gravity Porosity Permeability (md) Gas saturation Thickness (ft) Fracture half-width (ft) Reservoir outer boundary (ft)

5000 14.7, 1810, and 3220 290 0.717 0.15 0.05 1.00 92 392.7 1

1.E+01

1.E+00

q Gp ¼ 1 qi CGIP

456

Value

ð35Þ

D

460

455

Input

1.E+02

where  k is the viscosity–compressibility ratio calculated via its def is its time-averaged value as inition in Eq. (33) on a ROI-basis and b defined in Eq. (34). As per our preceding discussion about ROI, average density in Eq. (33) is calculated in terms of contacted gas-in-place (CGIP) using the appropriate material balance equation within the ROI. For dry gas reservoirs under volumetric depletion, Eq. (26) shows average density within the ROI is given by:

454

Table 1 Input data for early transient solution of linear flow.

qD

433

ð38Þ

!0:5 ð39Þ

To illustrate the proposed methodology and thus verify the validity of the transformation in Eq. (29) for transient flow, rate-time profiles are generated for the hypothetical reservoir undergoing linear flow as described in Table 1. Three different levels of bottomhole (fracture) pressure conditions are considered: wide-open production (pwf ¼ 14:7 psia); intermediate drawdown (pwf ¼ 1; 810 psia); and relatively low drawdown (pwf ¼ 3; 220 psia). Analytical liquid qliq D

solutions (corresponding to in Eq. (29)) are generated using the applicable liquid model, i.e., Eq. (1) inverted into real time space or Eq. (2) for infinite-acting linear flow. The liquid response for this scenario is plotted as a black solid line in Fig. 4. A suitable numerical simulator, or any other independent methodology able to model gas non-linearities, can be used to generate the associated gas responses. Gas responses in this study were obtained from an in-house simulator and independently via the Qanbari and Clarkson (2013) and Chen and Raghavan (2013) semi-analytical methodology. In these calculations, pseudo-critical properties are estimated via Sutton (1985), Z-factors and gas isothermal compressibility via Dranchuk and Abou-Kassem (1975) and Abou-Kassem

1.E-01 1.E-04

1.E-03

Liquid Soluon

1.E-02 tD Gas Soluon, pwf = 14.7 psia

1.E-01

1.E+00

Gas Soluon, pwf = 1810 psia

Gas Soluon, pwf = 3220 psia

Fig. 4. Early-transient gas responses results vs. analytical liquid response – linear flow.

et al. (1990), and gas viscosities using Lee et al. (1966). Resulting gas rate and time data were transformed into their density-based dimensionless values using the definitions in Eqs. (30) and (31), respectively, for plotting and comparison purposes. Results are shown in the same figure (Fig. 4), represented by red markers for every pwf value. Unlike the liquid responses, shown to be insensitive to the prevailing pwf value when plotted in terms of qD vs. t D variables, the gas qD vs. t D responses exhibit the drawdown dependency reported by Ibrahim and Wattenbarger (2006). In Fig. 5, the analytical liquid solution (black-solid line) has been  using average pressure of kb rescaled (Eq. (29)) using ROI-based  ROI volume. Rescaled solutions are represented by red-dash lines in the figure. It is shown that gas rescaled solutions (red-dash lines) are able to successfully match all gas rate-time responses derived from numerical simulation and the independent semi-analytical methodology (red markers). In all these cases, ROI volumes were calculated based on the radius of influence definition in Eq. (39) with a C r of 2.554 as suggested by Nobakht and Clarkson (2012a) k values used for linear flow. Fig. 6 plots the resulting ROI-based  in each of the rescaled solutions, from where it is seen that required  k values become time-independent. According to Eq. (33), a conk corresponds to a constant stant (time-independent) ROI-based  (time-independent) average density (and average pressure) within the ROI. Fig. 7 shows such time-independent average density used to evaluate the ROI-based  k plotted in Fig. 6. These figures demon q  are constants during infinite-acting linear  , andp k; b; strate that  flow period at each pwf condition. This is consistent with the findings of Nobakht and Clarkson (2012a), who first reported and analytically proved that average pressure inside the ROI remains constant throughout early transient flow period for gas linear flow

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

 model in Eq. kb The validity of the density-based rescaled  (35) for early-transient linear flow is thus verified. The availability pffiffiffiffiffi of a rescaled qD vs. tD model for early-transient linear flow  rescaled expokb analysis is akin to the availability of successful  nential models derived and successfully verified by Ayala and Zhang (2013b) and Zhang and Ayala (2014) for BDF analysis.  ¼ const early-transient linear Further, based on the observed  k¼b behavior, it follows that the ratio between the gas and liquid linear solutions is given by

1.E+02

qD

1.E+01

1.E+00

qgas D qliq D

1.E-01 1.E-04

1.E-03

1.E-02 tD

1.E-01

1.E+00

531 532 533 534 535 536

537

ð40Þ 539

wf

551

0.60

0.50

pffiffiffi 2 1 qD ¼ k  pffiffiffiffi  pffiffiffiffiffi p p tD

0.40

which in terms of dimensionless gas cumulative production ðGPD Þ is expressed as:

Gas Soluon, pwf = 14.7 psia

Gas Soluon, pwf = 1810 psia

Gas Soluon, pwf = 3220 psia

Rescaled Sol., pwf = 14.7 psia

Rescaled Sol., pwf = 1810 psia

0.70

0.30 5.E-05

GPD 5.E-04

5.E-03 tD

5.E-02

5.E-01

¼

Z

t D 0

pffiffiffi 4 k pffiffiffiffiffi  qD  dt D ¼ pffiffiffiffi  t D

p p

ð41Þ

ð42Þ

11.0

Eq. (42) demonstrates that gas cumulative production during pffiffi early-transient linear flow is directly proportional to t. Eqs. pffiffi (37)–(39) show that CGIP is also directly proportional to t. Therefore, the Gp =CGIP ratio in the ROI material balance statement in Eq. (36) becomes time-independent, leading to the observed  (and  kÞ characteristic of early-transient linear time-independent q flow first reported by Nobakht and Clarkson (2012a). Substituting these Gp and CGIP definitions into Eq. (36) (see Appendix B), the material balance statement for early-transient linear flow becomes:

10.5

  G 2 r pffiffiffi q ¼ qi  1  p ¼ qi  1  pffiffiffiffi q k CGIP p Cr

pwf = 14.7 psia

pwf = 1810 psia

pwf = 3220 psia

Fig. 6. ROI-based  k of early transient linear flow.

11.5

Average Density (lbm/3)

530

able behind all density-based transformations prescribed by Eq. (29). The density-based gas analytical solution for infinite-acting linear flow becomes:

Liquid Soluon

0.80



10.0

9.5

9.0 5.E-05

5.E-04

pwf = 14.7 psia

5.E-03 tD pwf = 1810 psia

5.E-02

5.E-01

pwf = 3220 psia

 of early transient linear flow. Fig. 7. q

527

529

540

Fig. 5. Early-transient gas responses results, analytical liquid response, and  rescaled solution – linear flow. ROI-based  kb

526

ffi pffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k lgi cgi    ffiffiffi p ¼ k¼ ¼ const ¼  l g cg q ;q b

528

i.e., solutions are offset by a constant value equal to the square root of the viscosity–compressibility product in Eq. (33). The offset would be non-existent for fluids whose viscosity and compressibility values are independent of pressure (e.g., liquids, for which  k ¼ 1Þ. For all other fluids, this offset becomes a function of the pwf specification and the fluid viscosity–compressibility product dependency with pressure. This is again consistent with Nobakht and Clarkson (2012a) who showed that the correction factor ðf c Þ pffiffi to the liquid qgsc and t (actual time) relationship is equal to the square-root of the ratio between the lg ct product at initial and  rescaling variROI average pressure. In this work, this ratio is the k

Rescaled Sol., pwf = 3220 psia

525

7

under a constant-pwf production specification. It should be noted  is constant, as per the definition in Eq. (34), that when k ¼  is constant as well ð¼  b k – i.e., b kÞ.



ð43Þ

541 542 543 544 545 546 547 548 549 550 552 553

554 556 557 558

559 561 562 563 564 565 566 567 568 569 570

571 573

i.e., the ROI average density is constant during early-transient pffiffiffi k, the density drawdown ratio linear flow and a function of  ðr q Þ, and the C r value used in the radius of influence expression. This ROI density remains constant because the amount of mass per unit ROI volume added by the expanding control volume balances the mass of fluid leaving the ROI per unit ROI volume. The pffiffiffi  rescaling correction needed for the k density-based

574

infinite-acting linear models in Eqs. (41) and 42 is directly calculated from its definition once initial pressure, temperature, gas specific gravity and prevailing pwf is known:

581

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi u l cgi gi k ¼ u t q q  wf 2h ww 

575 576 577 578 579 580 582 583

584

ð44Þ

wf

where the ROI average density is calculated simultaneously through Eq. (43). Few successive substitutions would be needed, initiated

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

586 587 588

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 8

591 592

wf

595

relatively low drawdown (pwf ¼ 3; 220 psia), which are the values in display in Fig. 6 as originally calculated during the forward-solution.

596

Density-based early-transient solution: radial flow

597

600

 transformation stated in Eq. (29), kb On the basis of the same  and using the corresponding radial liquid counterpart in Eq. (7), the density-based rate-time analytical solution for infinite-acting gas  rescaled model: kb radial flow is readily written as the following 

603

   1 1 1 1 1 1   ¼    Ei    ffi   ln bt D þ 0:80907  qD k 2 4bt D k 2

593 594

598 599

601

604 605 606 607 608 609 610

611

613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648

5.E-01

r inf ¼ 2 

a1 kt /lgi cgi

2.E+03

ð46Þ

2.E+04 tD Gas Soluon, pwf = 14.7 psia

2.E+05

2.E+06

Gas Soluon, pwf = 1810 psia

Gas Soluon, pwf = 3220 psia

ð45Þ

!0:5 þ rw

5.E-02 2.E+02

Liquid Soluon

 is its time-averaged value defined k and b where, again, ROI-based  in Eqs. (33) and (34), respectively. Average ROI density in Eq. (33) is calculated in terms of contacted gas-in-place (CGIP) using the appropriate ROI material balance equation – e.g., Eq. (36) for dry gas reservoirs under volumetric depletion. CGIP is calculated from Eq. (37); with the following updated early transient radial definitions for Ar and r inf :

  Ar ¼ p r 2inf  r2w

qD

590

pffiffiffi pffiffiffi  ¼ qi in Eq. (44), to obtain  k. This  k calculation applied to using q the case study above yields, unsurprisingly,  k ¼ 0:349 for the k ¼ 0:549 for the wide-open production case (pwf ¼ 14:7 psia);   ¼ 0:735 for the intermediate drawdown (p ¼ 1; 810 psia); and k

Fig. 8. Early-transient gas responses results vs. analytical liquid response – radial flow.

5.E-01

qD

589

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

ð47Þ

To showcase the applicability of the density-based technique to early transient radial flow, radial liquid responses and gas rate-time responses were simulated using the same reservoir and fluid data of the linear case (Table 1) but using with a characteristic length of rw ¼ 0:25 ft in lieu of xf . The analytical infinite-acting liquid solution (black-solid line in Fig. 8) was obtained via Stehfest (1970) inversion of Eq. (5) and also via the exponential–integral solution in Eq. (7) to double check for consistency. Gas rate-time responses were independently generated from an in-house simulator and further verified with the radial extension of the Qanbari and Clarkson and Chen and Raghavan (2013) semi-analytical methodology derived by this work in Appendix A. Gas solutions are displayed in red markers in Fig. 8, after being transformed into dimensionless values using the density-based dimensionless radial definitions in Eqs. (30) and (32). It is noted that gas solutions also exhibit a drawdown dependency. In Fig. 9, the analytical liquid solution  variables in Eqs. kb (black-solid line) is rescaled by ROI-based  (33) and (34). The resulting density-based rescaled solutions are shown as red-dash lines, and they show that they closely mimic gas rate-time responses (red markers). In these examples, ROI volumes are calculated based on radius of influence (Eq. (28)) using the classical radial C r value of 2. It should be noted that the rate of decline rate in radial flow is much less pronounced than in linear flow. While linear rates declined at a 1:2 qD to t D log-cycle ratio in Fig. 5, radial rates decline at about a 1:8 qD to tD log-cycle ratio in Fig. 9 for the case under consideration. Fig. 10 plots ROI-based  k values generated by the methodology to obtain the rescaled solutions. Fig. 11 shows values of corresponding ROI average density. Both figures display a clear trend: the ROI is not subject to density depletion during infinite-acting radial flow. Rather, ROI average density, and hence ROI-based  k values, increases with time. Eq. (36) explains that this happens when CGIP increases at a faster pace than Gp as the control volume expands. This is because contacted area, and thus CGIP, increases quadratically with r inf in radial flow (Eq. (46)) compared to linearly

5.E-02 2.E+02

2.E+03

2.E+04 tD

2.E+05

2.E+06

Liquid Soluon

Gas Soluon, pwf = 14.7 psia

Gas Soluon, pwf = 1810 psia

Gas Soluon, pwf = 3220 psia

Rescaled Sol., pwf = 14.7 psia

Rescaled Sol., pwf = 1810 psia

Rescaled Sol., pwf = 3220 psia

Fig. 9. Early-transient gas responses results, analytical liquid response, and  rescaled solution – radial flow. ROI-based  kb

0.85

0.75

0.65

0.55

0.45

0.35 2.E+02

2.E+03

pwf = 14.7 psia

2.E+04 tD pwf = 1810 psia

2.E+05

2.E+06

pwf = 3220 psia

Fig. 10. ROI-based  k of early transient radial flow.

in linear flow (Eq. (37)). At the same time, radial Gp increases at a much-reduced pace compared to linear flow given the much-reduced area available for flow at the wellbore. This causes the rate of addition of mass within the expanding control volume

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

649 650 651 652

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 9

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

12.4

Average Density (lbm/3)

12.2

11.6

2.E+03

2.E+04 tD

pwf = 14.7 psia

2.E+05

pwf = 1810 psia

2.E+06

pwf = 3220 psia

 of early transient radial flow. Fig. 11. q

655 656 657 658 659 660 661 662 663

to be larger than the ability to produce fluids out of the ROI through the wellbore flow area – leading to an increasing ROI density with time. During BDF, however, drainage area density would start to decrease (i.e., deplete) for both linear and radial cases once the control volume is fixed and no longer expanding.  model in Eq. (45) kb The validity of the density-based rescaled  for early-transient radial flow is thus verified. The dislocation or offset between gas and liquid radial solutions is time-dependent, and thus gas infinite-acting radial solutions require the calculation  traces as a function of t . The q vs. t  gas type-curve is genkb of  D D D erated following the rescaled model:

664 666 667 668 669 670 671 672

qD

  1 1 1 ¼ k  Ei    2 4bt D

676 677

ð48Þ

It should be noted that for the purposes of building a type-curve, the more rigorous finite-well Eq. (5) for constant-pwf infinite-acting radial systems can be used in lieu of its line-source approximation in Eq. (48). The Stehfest (1970) inversion of Eq. (5) would be suitable for such purpose. Once qD vs. tD data is generated, dimensionless cumulative production is obtained as:

673 675

GPD ¼

Z 0

t D



qD  dtD

ð49Þ

with the ROI density trace for radial flow is calculated from (see Appendix B):

678

680 681

q ¼ qi  1  2rq 

685 686 687 688 689 690 691 692 693 694

!

GPD ðC r 

ðtD Þ0:5

2

ð50Þ

þ 1Þ  1

698

In previous section, we showed that the density-based technique successfully predicts rate-time responses from gas reservoirs exhibiting extended infinite-acting linear and radial flow behaviors. The key to the methodology is evaluating rescaling variables   within ROIs. In this section, the density-based technique is kb used to solve inverse problem, i.e., the evaluation of rate-time responses for GIP inside ROIs (CGIPs) estimations. The rate-time analysis procedure is summarized as follows:

699

Step 1: Select the flow model that best fits the rate-time early transient data. - Linear flow if the rate-time data is best described by a pffiffi straight line in a cartesian 1=qgsc vs. t plot (Eq. (41)); - Radial flow if the rate-time data is best described by a straight line in a semi-log 1=qgsc vs. t plot (Eq. (45)); Step 2: Generate the gas analytical solution (rescaled solution) for early transient flow period. - For linear flow, the rescaling entails the calculation of pffiffiffi   k using Eq. (44), and the q t a single value of 

707

responses using Eq. (41); k - For radial flow, the rescaling entails the calculation of   profiles as a function of dimensionless t  using and b  D  tD from their definitions (Eqs. (33) and (34)) and the q   Eq. (50) or Eq. (51). The qD t D gas profile can then be calculated using Eq. (45); Step 3: Determine time-multiplier ðXÞ and flow rate-multiplier ðYÞ that providing the best match against the analytical gas rate-time responses generated in Step 2. This can be done either through type-curving or straight-line matches. For linear flow, it should be noted that type-curving against a straight-line (linear flow case) yields multiple possible pffiffiffiffi Y  X pairs but with the same Y X ratio. Instead of finding individual X and Y values through type-curve matching, pffiffiffiffi the analyst may calculate the value Y X directly from the slope of the plots:

718

q ¼ qi  1 

2rq C 2r



GPD t D

D

pffiffiffiffi rq p p pffiffiffiffi pffiffi ¼ pffiffiffi  Y X  t qgsc 2 k pffiffiffi pffiffi GP 4 k ¼ pffiffiffiffi pffiffiffiffi  t rq p p  Y X

D

697

! ð51Þ

b  calculations start at least 5 log-cycles prior In this study, k the desired minimum tD for type-curve construction purposes. k-guess based on initial Calculations are initialized with a first   ¼ qi Þ. In this first time step, iterations are required pressure ðq k is attained via successive applications of the until a converged  rescaled solution (Eq. (48)) and material balance statements (Eq. (50)) using the short-time approximation GPD ¼ qD  t D for Eq. (49). The iterative process continues until there is no further k is attained. Throughout the iteration improvement on this first  ¼ for initial  k; b k which is valid for initial conditions. After the first

GP ¼ 2t qgsc

700 701 702 703 704 705 706

708 709 710 711 712 713 714 715 716 717 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734

735

ð52Þ ð53Þ

which follow from Eqs. (41) and (42) using the dimensionless variables definitions in Appendix B. Also, the multiplication of both expressions above yields:

which at later times (when r inf >> r w ) becomes:

682

684

Early transient rate-time analysis applications

696

11.8

11.2 2.E+02

654

695

12.0

11.4

653

  pair is obtained, the solution marches in small t increments kb D  from the previous time kb with new qD values evaluated using   via their definitions. kb step and new 

ð54Þ

737 738 739 740

741 743

i.e., cartesian plots of Gp =qgsc vs. time would yield a straight lines, which is the theoretical basis for Doung’s linear flow forecasting model (Duong, 2011). It is interesting to note pffiffiffi k and r q dependenthat this multiplication eliminates the 

744

cies in the linear-flow gas forecasting equation. For radial flow, straight-line analysis is also possible via the use of the logarithmic approximation in Eq. (45):

748

rq qgsc

 Y   þ 0:80907 ¼  ln bXt 2k

ð55Þ

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

745 746 747 749 750

751 753

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 10

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

758

The implementation of this equation requires iterations for  as a function of actual time. This the calculation of  k and b study has found that, for constant pwf data, an excellent estimation of the Ymultiplier may be obtained via the following approximation to the logarithmic expression above:

761

1 Y ¼ lnðtÞ þ intercept qgsc 2rq ki

754 755 756 757

759

764 765 766

767

CGIP ¼

769

C r pffiffiffiffiffiffiffiffi  pffiffiffiffi  t last p Y X 2

770

Radial Flow:

773

CGIP ¼

771

774 775 776 777 778 779 780

781

1.E+02

  ki Þ. k is explicitly evaluated at initial pressure ðk where  Step 4: Calculate CGIP using X  Y data from Step 3 and the last point of early transient (infinite-acting) data. CGIP expressions are derived in Appendix B, reproduced below. Linear Flow:

1.E-01

C 2r  t last 2Y

ð58Þ

 pffiffiffi  1 q 1 pffiffiffiffiffiffiffiffiffiffiffiffi 1 pffiffiffiffiffi sc pffiffiffiffi lgi cgi  pffiffiffiffi xf k ¼ 2p a1 qi h / Y X

1.E-02 1.E-05

785 786

787

 kh ¼

790 789 791

ð60Þ

1.E-03

1.E-02

1.E-01

1.E+00

1.E+01

1.E+02

Type-Curve (Rescaled Early Trans. Sol.) Fig. 12. Type-curve with gas analytical linear solution – Example 1 – Step 2.

1.E+02

ð59Þ



qsc 1 l cgi  2pa1 qi gi Y 1

1.E-04

tD

For radial flow, the formation capacity (kh) may be estimated from Eq. (60) as a function of the Y-multiplier (as derived in Appendix B).

784

1.E+00

ð57Þ

Step 5: For linear flow, the product of fracture half-length and square root of permeability may be estimated from Eq. (59) (as derived in Appendix B). Reservoir porosity and thickness are additional data required to do this calculation. From this, xf values may be estimated if k estimates are available.

783

1.E+01

qD

763

ð56Þ

1.E+01

qD

762

5:45  109 (day0.5/SCF). Fig. 14 shows the corresponding pffiffi straight-line cartesian 1=qgsc vs. t analysis, from whose pffiffiffiffi slope one obtains Y X ¼ 5:45  109 (day0.5/SCF).

1.E+00

1.E-01 792

Example 1. Production data analysis of early transient linear flow.

793

Step 1, 2: In this example, numerically generated data for infinite-acting linear flow conditions based on the inputs listed in Table 2 is analyzed. The type-curve rescaled solution corresponding to the conditions in Table 2 is shown in Fig. 12. For linear flow, we use C r   of 2.554. For this case, r q ¼ qi  qwf =qi ¼ 0:8979.

795 796 797 798 799 800 801 802

803

1.E-02 1.E-05

807 808

1.E-02

1.E-01

1.E+00

1.E+01

1.E+02

Rate-Time Responses

Fig. 13. Type-curve match – Example 1 – Step 3.

8.0E-07 7.0E-07

1

X ¼ 4:162  103

day

Y ¼ 8:445  108

SCFD1

6.0E-07

It is again noted that this type-curving against a straight-line yields multiple possible Y  X pairs but pffiffiffiffi pffiffiffiffi with the same Y X ratio, equal in this case to Y X ¼

806

1.E-03

Type-Curve (Rescaled Early Trans. Sol.)

Step 3: Multipliers X and Y that can match the actual rate-time against the gas analytical solution are determined in this step. Values of multiplier X and Y that yield the type-curve match in Fig. 13 are:

805

1.E-04

tD

1/qgsc (SCFD-1)

794

y = 2.65E-08x

5.0E-07 4.0E-07 3.0E-07 2.0E-07

Table 2 Input data for early transient analysis of linear flow.

1.0E-07 Input

Value

Initial reservoir pressure (psia) Bottomhole flowing pressure (psia) Reservoir temperature (F) Gas specific gravity

5000 480 290 0.717

0.0E+00 0

5

10

15

20

25

30

me0.5 (day0.5) Fig. 14. Straight-line match – Example 1 – Step 3.

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

809 810 811 812

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 11

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

814 815 816 817

818 820

821 822 823 824 825 826 827

Step 4: CGIP is estimated. The last point of data is given at t = 772.5 days with no signs of BDF (i.e., no deviation from infinite-acting behavior) as per the matches above. By substituting all values into the CGIP equation, we have:

CGIP ¼ 8:38 BSCF

1.E+00

qD

813

Step 5: In this last step, the fracture half-length and square root of permeability product is evaluated. Reservoir thickness is equal to 92 ft and reservoir porosity is equal to 0.15, which their actual values used as inputs for the numerical simulator. Density of gas at initial and standard conditions are equal to 12.46 lbm/ft3 and

1.E-02 1.E+03

3

5:488  102 lbm=ft , respectively. Gas viscosity at initial conditions is equal to 0.026 cp, gas compressibility

828 829

at initial conditions is equal to 1:418  104 psi1. Substitution these values into Eq. (59) yields

830 831

832

n o 3 5:488  10 lbm=ft pffiffiffi 1 1 n o pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi  xf k ¼  3 2p a 1 92fftg 0:15 12:46 lbm=ft rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n offi 1  0:026fcpg  1:418  104 psi

1.E-01

1.E+04

1.E+05

1.E+06

1.E+07

1.E+08

tD Type-Curve (Rescaled Early Trans. Sol.) Fig. 15. Type-curve with gas analytical radial solution – Example 2 – Step 2.

2



1 n o 0:5 5:45  10 day SCF1 9

qD

834

1.E+00

835

or

838

pffiffiffi 0:5 xf k ¼ 87:2 ft-md

836

pffiffiffi The actual xf k used to simulate production data is equal to 82.8 ft-md0.5. The accuracy of this value, as well as CGIP, is subject to the proper selection of C r -values in the radius of influence equation during the calculation of pffiffiffi  correction. the k

839 840 841 842 844 843 845

1.E-01

1.E-02 1.E+03

1.E+04

1.E+05

1.E+06

1.E+07

1.E+08

tD Type-Curve (Rescaled Early Tran. Sol.)

Rate-Time Responses

Fig. 16. Type-curve match – Example 2 – Step 3. 846 847

Example 2. Production data analysis of early transient radial flow 1.6E-06

849 850 851 852 853 854 855 856 857 858

859

Step 1, 2: In this example, numerically generated data for infinite-acting radial flow conditions based on the inputs listed in Table 3 are used. The type-curve rescaled solution corresponding to the conditions in Table 3 is shown in Fig. 15. For radial flow, we use the traditional radial C r value of 2.0. For this case,   r q ¼ qi  qwf =qi ¼ 0:8979. Step 3: Multipliers X and Y matching the actual rate-time data with the gas analytical type-curve radial solution are obtained in this step. Values of multiplier X and Y that yield the type-curve match in Fig. 16 are: 3

X ¼ 9:41  10

Y ¼ 8:89  108

861

1.5E-06 y = 9.47E-08x + 9.22E-07 1.4E-06

1/qgsc (SCFD-1)

848

1.3E-06 1.2E-06 1.1E-06 1.0E-06 9.0E-07

1

day

8.0E-07 -1.00

SCFD1

0.00

1.00

2.00

3.00 log(me)

4.00

5.00

6.00

7.00

Fig. 17. Straight-line match – Example 2 – Step 3. Table 3 Input data for early transient analysis of radial flow. Input

Value

Initial reservoir pressure (psia) Bottomhole flowing pressure (psia) Reservoir temperature (F) Gas specific gravity

5,000 480 290 0.717

Fig. 17 shows the straight-line analysis using 1=qgsc vs. logðtÞ plot approximation in Eq. (56). The slope of this 8

1

straight-line plot is equal to 9:47  10 SCFD , and ki , which can be calculated from Eq. (33), is equal to 0.5193. Substitution of these values (with r q ¼ 0:8979) into the slope expression of Eq. (56) yields

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

862 863 864 865 866 867

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 12 868 869 870 871

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

Y ¼ 8:83  108 SCFD1, which is very close to the value rigorously obtained through the type-curve match. Step 4: CGIP is calculated by substituting all values for tlast ¼ 522:9 days. We obtain:

872 874 875 876 877 878

CGIP ¼ 11:765 BSCF Step 5: In this step, formation flow capacity (kh) can be calculated from Eq. (60). Density of gas at initial, standard, and bottomhole conditions are equal to 12.46 lbm/ft3, 3

881

5:488  102 lbm=ft , and 1.272 lbm/ft3, respectively. Gas viscosity at initial conditions is equal to 0.026 cp, gas compressibility at initial conditions is equal to

882

1:418  104 psi1.

879 880

883

n o 5:488  102 1 kh ¼  0:026fcpg  1:418  104 psi 12:46 2 p a1 1

 885

!

1 n o 8 SCFD1 8:89  10

886

or

887 889

kh ¼ 4:6 md-ft The actual kh used to generate production data is exactly equal to 4.6 md-ft.

890 891

Uncited references

927

Nobakht and Clarkson (2012b,c), Nobakht et al. (2012), Zhang (2013).

928

Acknowledgement

930

The authors would like to express their gratitude to the members of the Unconventional Natural Resources Consortium (UNRC) at Penn State U. for their support of our work in the area of unconventional phenomena analysis in gas plays.

931

Appendix A: Semi-analytical early-transient solution for gas radial flow

935

Following Qanbari and Clarkson (2013) and Chen and Raghavan (2013) methodology for linear flow, a semi-analytical solution of early transient (infinite-acting) period for non-linear gas radial flow has been derived in this work for comparison purposes.

937

Gas diffusivity equation

941

The continuity equation may be written as:

@ r  ðqv Þ ¼ ð/SqÞ @t

ðA-1Þ

893

!   1 @ k @p @ / ¼ r r @r @t Bg lg Bg @r

On the basis of this study, the following conclusions are drawn:

895

(1) Density-based analysis method can successfully capture early transient behavior of gas reservoirs under constant pwf production condition by incorporating region of influence   b,  calcuconcept into depletion-driven rescaling variable, k

896 897 898 899 900 901 902

(2)

903 904 905 906 907 908

(3)

909 910 911 912 913 914

(4)

915 916 917 918 919 920 921 922 923 924 925

(5)

lations. Numerically simulated and semi-analytical results are used to verify the modified density-based solutions on both linear and radial flow geometries.  which calculated k  b, The corrected rescaling variables  from average pressure inside region of influence, are constant during early transient period of linear flow. This behavior corresponds to a constant average pressure behavior within the ROI, which is consistent with reported findings by previous works.  increase with time durFor radial flow, the ROI-based  kb ing infinite acting conditions. This finding, previously unreported, corresponds to an increasing average pressure inside the region of influence during infinite-acting conditions. This behavior has been verified through numerical and analytical arguments. Early transient behavior of gasradial flow under constant pwf production can be modeled using a newly proposed semi-analytical solution based on the extension of an existing linear-flow methodology. This solution works by using similarity variable to transform the governing PDE into non-linear ODE, which then solved by numerical integration. ODE non-linearity is handled iteratively. Contacted Gas In-Place, the product of fracture half-length and square root of permeability (for linear flow) and the formation flow capacity (kh, for radial flow) can be effectively predicted using newly proposed density-based analysis procedure.

939 940

945 946 947

948

ðA-3Þ

951 952

953 955 957

p

2 w ¼ a  dp ai p0 a a g¼ 0¼ @b=@p b k a¼ lg Bg / Bg

ðA- 4Þ ðA- 5Þ ðA- 6Þ ðA-7Þ

Dimensionless diffusivity equation

  1 @ @w 1 @wD rD D ¼ rD @r D @r D gD @tD

961

962

ðA-8Þ

where definitions of dimensionless variables are given by

wD ¼

w0i  w0 w0i  w0wf

rD ¼

r rw

959

960

Eq. (A-3) can be written in dimensionless form as

tD ¼ a1

938

956

Z

b¼

936

950

where 0

934

ðA-2Þ

By using definition of (normalized) pseudopressure function (Al-Hussainy et al., 1966), Eq. (A-2) can be re-expressed as

  1 @ @w0 1 @w0 ¼ r r @r @r g @t

933

943

If we assume 1D (r-direction) radial flow in single-phase gas reservoir and applying Darcy’s law, Eq. (A-1) can be written as

894

932

942

892

Concluding remarks

929

964 965

966

ðA-9Þ

968 969

ðA-10Þ

gi r 2w

971 972

t

ðA-11Þ

926

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

974

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 13

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

975 977 978 979 980 981

982 984

gD ¼

wD ¼ 0 at t D ¼ 0 and r D P 1

993

999

1000

wD ¼ 0 at t D > 0 and r D ! 1

ðA-15Þ

Introducing the Boltzmann’s variable ðnÞ as

rD 2t 0:5 D

1002 1003

where initial and outer boundary conditions are combined into

1004 1006

wD ¼ 0 when n ! 1

1007

and inner boundary condition is transformed into

1010

1011

1013 1015

1016 1017

1018

ðA-18Þ

1 2t0:5 D

    @ @w 2n @w n D ¼ n D @n @n gD @n

wD ¼0

wD ¼1

ðA-21Þ

0 1   Z 1 Z n ^ n ^A @wD 1 @ dwD ¼ 1 ¼ n   exp 2 dn  dn @n n¼ntd ntd n ntd gD

1029

ðA-22Þ

1030

If we multiply both sides of Eq. (A-21) by dn=n and integrate from wD ¼ 1towD ðn ¼ ntd to nÞ, we will also have

1031

1032

0 1   Z n Z n ^ n @wD 1 dwD ¼ wD  1 ¼ n   exp @2 d^nA  dn @n n¼ntd ntd n wD ¼1 ntd gD

Z

1034

1037

1040 1041

1042

ðA-25Þ 1044

r¼r w

!   pkh wi  wwf @wD n lgi Bgi @n

1045 1046

1047

ðA-26Þ 1049

n¼ntd

or

1050

1051

  1 @w qD ¼ n D 2 @n n¼ntd

ðA-27Þ

where dimensionless flow rate is defined by

qD ¼

1055

qgsc lgi Bgi qgsc T lgi Z i p   ¼ sc   2pa1 kh w0  w0 2pT sc khp w0  w0 i i wf i wf

ðA-28Þ 1057

which is a well-known constant. As mentioned by Qanbari and Clarkson (2013), qD calculated from Eq. (A-27) is inaccurate. Alternative expression must be used to achieve a more accurate qD value, which can be derived by multiplying both sides of Eq. (A-20) by dn:

  @w 2n2 d n D ¼ dw @n gD D and integrate both sides of the resulting n ¼ ntd to 1ðwD ¼ 1 to 0Þ, we will have

1058 1059 1060 1061 1062 1063

1064 1066

equation

from

1067 1068

1069

    Z wD ¼0 2 @w @w 2n n D  n D ¼ dw @n n¼1 @n n¼ntd gD D wD ¼1

1071 1072

1073

ðA-29Þ

Given that @wD =@n as n ! 1 is equal to zero. Substituting Eq. (A-29) into Eq. (A-27) yields

qD ¼ 

1053 1054

1

  Z wD ¼1 2 @w 2n n D ¼ dw @n n¼ntd gD D wD ¼0

If we multiply both sides of Eq. (A-21) by dn=n and integrate from wD ¼ 1 to 0 ðn ¼ ntd to 1Þ, we will have

Z

1036

1039

or

@n n¼n td

1024

1027

ðA-20Þ

    n@w@nD Z n ^ n¼n @wD 2n ^ log n ¼ dn @w @n ntd gD ðn D Þ

or

1026

ðA-19Þ

If we multiply both sides of Eq. (A-20) by ½dn=nðdwD =dnÞ and integrate them from n ¼ ntd to n, we will have

!     Z n ^ n ^ @wD @wD n ¼ n  exp 2 dn @n n¼n @n n¼ntd ntd gD

1025

qgsc

!

1021

1022

1035

ðA-24Þ

D

  Note that value of psc = 2pT sc  6:328  106 is equal to 711,

From Eq. (A-17)

1020

ðA-17Þ

ODE solution

1012

n

! kð2prhÞ @p ¼ lg Bg @r

qgsc ¼

and substituting r D ; @rD , and @tD with appropriate form of n, Eq. (A-8) can be transformed into ODE form as shown below:

ntd ¼

td

From Darcy’s equation,

ðA-16Þ

wD ¼ 1 when n !

R1

 ^n ^ gD dn  dn  ^n ^ n  dn d g

If we substitute definition of rD ; n, and wD into Eq. (A-25), we will have

ODE diffusivity equation

n¼

ntd

ntd

ðA-14Þ

    @ @w 2n @w n D þ n D ¼0 @n @n gD @n

1008

wD ¼ 1 

Outer boundary condition:

995

998

Rn

 Rn 1  exp 2 ntd n  Rn 1  exp 2 n

Flow rate solution

wD ¼ 1 at tD > 0 and r D ¼ 1

994

997

By dividing Eq. (A-23) with Eq. (A-22), we can write expression for wD as

ðA-13Þ

Inner boundary condition:

989

990 992

ðA-12Þ

Note that a1 is a unit conversion factor. For infinitely large reservoir produced under constant bottomhole flowing pressure constrain, initial and boundary conditions are given by Initial condition:

985

986 988

g gi

Z

wD ¼1

wD ¼0

n2

gD

1075 1076 1077

1078

dwD

ðA-30Þ

1080

which is an alternative for qD calculation. Note that qD will be positive because the minus sign on the RHS of this equation will be cancelled with negative value of the integral (n decreases as wD increase).

1081

Numerical verification

1085

Eq. (A-24) and Eq. (A-30) derived above are the basis of the semi-analytical solution methodology. For dimensionless pseudopressure:

1086

1082 1083 1084

wD

ðA-23Þ

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

1087 1088

1089

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 14

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

Rn wD ¼ 1 

1092

1095

R1 ntd

1091 1093

ntd

 Rn 1  exp 2 ntd n  Rn 1  exp 2 ntd n

 ^n  dn d gD  ^n ^ g dn  dn ^n

ðA-24Þ

D

qD ¼ 

Z

wD ¼1

wD ¼0

n2

gD

dwD

1109

    R nmax þ 1 R n ^n  ^n  dn   exp 2 d 0 gD nmax   n  6 tolerance 1      R R nmax 1 n ^n ^   exp 2 0 g dn  dn   ntd D n

1099 1100 1101 1102 1103 1104 1105

1107

1110 1111 1112 1113 1114

ðA-31Þ

1119

Type-curve ready mathematical development in this section follows the development for boundary-dominated flow presented by Ayala and Ye (2013a).Linear flow: The linear definitions of dimensionless flow rate and dimensionless time are:

1120

qD

qsc lgi cgi

q   qgsc ¼ Y  gsc ¼ rq 2pa1 kh qi  qwf

1117 1118

tD ¼

a1 k t ¼Xt /lgi cgi x2f

ðB-2Þ

where X and Y are multipliers which render flow rate and time dimensionless, and rq is the dimensionless density drawdown:

qi  qwf rq ¼ qi XY ¼

ðB-3Þ

qsc

Input

Value

Initial reservoir pressure (psia) Bottomhole flowing pressure (psia) Reservoir temperature (F) Gas specific gravity Porosity Permeability (md) Gas saturation Thickness (ft) Wellbore radius (ft) Reservoir outer boundary (ft) OGIP (BSCF)

5000 14.7, 1810, 3220 290 0.717 0.15 0.05 1.00 92 0.25 500 2.469

and

Z

1 q 1 pffiffiffiffiffiffiffiffiffiffiffiffi 1 pffiffiffiffiffi sc pffiffiffiffi lgi cgi  pffiffiffi 2p a1 qi h / xf k

ðB-5Þ

qgsc  dt

ðB-6Þ

If we substitute Eq. (B-1) and (B-2) into Eq. (B-6), we have, for constant pwf production:

rq Gp ¼ XY

Z 0

t D

qD 

 dt D

¼

r q GPD

ðB-7Þ

XY

where dimensionless gas production for the transient period is R t  GPD ¼ 0D qD  dtD . For linear flow, Contacted Gas In Place is volumetrically calculated as:

1.E+06

CGIP ¼

4xf r inf h/qi

ðB-8Þ

qsc

1.E+04 1.E-01

2 r pXY inf ;D

ðB-9Þ

where r inf ;D is expressed as a function of tD by combining definition of rinf (Eq. (38)) with definition of t D (Eq. (B-2)).

rinf ;D

1.E+00

1.E+01 1.E+02 Time (day)

1.E+03

1.E+04

Iterave Integral @ Pwf 14.7 psia

Iterave Integral @ Pwf 1810 psia

Iterave Integral @ Pwf 3220 psia

In-house Simulator @ Pwf 14.7 psia

In-house Simulator @ Pwf 1810 psia

In-house Simulator @ Pwf 3220 psia

Fig. A-1. Gas rate comparison for early transient solution verification of radial flow.

1127 1128 1129

1130 1132

1136

1140 1141 1142

1143

t

If we substitute Eq. (B-4) into Eq. (B-8), resulting equation can be written as

1.E+05

1125

1138



0

CGIP ¼

1.E+07

1124

1137

pffiffiffiffi  Y X¼

Gp ¼

1123

1134

pffiffiffi pffiffiffiffi Eq. (B-5) can be used to determine xf k once Y X is known. Recall the definition of cumulative gas production

Table A-1 Input data for early transient solution verification of radial flow.

1122

1133

ðB-4Þ

2ph/x2f qi

1121

ðB-1Þ

Based on these multiplier definitions, it follows:

Gas flow rates calculated from Eq. (A-24) are verified with results from numerical simulator. Input reservoir data is shown in Table A-1. Calculated flow rates are plotted in Fig. A-1. We can see that gas rates calculated from Eq. (A-24) and transformed using Eq. (A-28) match very well with the results from numerical

Gas Rate (SCFD)

1098

Appendix B: Early transient rate-time dimensionless equations

1116

ðA-30Þ

1106

1097

1115

For dimensionless flow rate:

These radial-flow expressions resemble the linear-flow expressions derived by Qanbari and Clarkson (2013). There are a number of important differences. First, the characteristic length is r w for pffiffiffiffiffi radial flow instead of Ac used in linear flow and inner boundary condition is specified at r D ¼ 1, instead of r D ¼ 0. Thus, lower limit of integration on the RHS of Eq. (A-24) becomes time dependent variable, n ¼ ntd ðnatrD ¼ 1Þ. Second, outer boundary condition used in wD evaluation is not practical (n ! 1 and wD ! 0Þ and has to be replaced. Undisturbed boundary condition ðn ! 1Þ is replaced by (n ! nmax Þ and nmax is defined as a location that satisfy the following condition

1096

simulator. However, there are small differences during early time which may potentially cause by discretization of the simulator. Large differences also occur at late time basically because of boundary dominated effect.

r inf C r a1 kt ¼ ¼  /lgi cgi xf xf

!0:5 ¼ C r  t 0:5 D

1146 1147

1148 1150 1151 1152 1153

1154 1156 1157 1158

1159 1161 1162 1163

1164

ðB-11Þ 1166

If we assume dry gas volumetric depletion, simplify material balance equation can be expressed as

  q Gp ¼ 1 qi CGIP

1145

ðB-12Þ

By substituting Eqs. (B-7) and (B-9) into Eq. (B-12), we will have

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

1167 1168

1169 1171 1172

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

1173 1175

prq GPD q ¼1 2 r inf ;D qi

ðB-13Þ

1177

or, after substituting Eq. (41) for flow and Eq. (B-11)

1180

q 2 r q pffiffiffi k ¼ 1  pffiffiffiffi qi p Cr

1176

1178

1181 1182

1183 1185

1186 1187 1188

1189

GPD

in infinite-acting gas linear

ðB-14Þ

In addition, the relationship between CGIP and actual time can be derived by substituting Eqs. (B-4) and (B-11) into Eq. (B-9):

pffiffiffiffiffi 2 C r pffiffi 2 pffiffiffiffi  t CGIP ¼ C Xt ¼ pXY r pY X

ðB-15Þ

Radial flow: The radial definitions of dimensionless flow rate and dimensionless time are:

qsc lgi cgi

qD ¼

q   qgsc ¼ Y  gsc rq 2pa1 kh qi  qwf

ðB-16Þ

t D ¼

a1 k t ¼Xt /lgi cgi r 2w

ðB-17Þ

1191 1192 1194 1195 1196

1197 1199 1200

1201

1202 1204

where X and Y are multipliers which convert dimensional flow rate and time to dimensionless variables. It follows that

XY ¼

qsc

ðB-18Þ

2

2p/hrw qi

It is also noted that the multiplier Y can be used calculate formation capacity kh by rearranging its definition to show:

kh ¼





qsc 1 l cgi  2pa1 qi gi Y 1

ðB-19Þ

1206

For radial flow, Contacted Gas In Place is volumetrically calculated from

1209

p r2inf  r2w h/qi CGIP ¼ qsc

1205

1207





ðB-20Þ

1211

If we substitute Eq. (B-18) into Eq. (B-20), resulting equation can be written as

1214

CGIP ¼

1210

1212

1215 1216

1217

1219

 1  2 r 1 2XY inf ;D

ðB-21Þ

where r inf ;D can be expressed as a function of t D by combining definition of rinf (Eq. (46)) with definition of t D (Eq. (B-17)):

r inf ;D

r inf C r a1 kt ¼ ¼  /lgi cgi rw rw

!0:5 þ

rw ¼ C r  t 0:5 D þ1 rw

ðB-22Þ

1221

For dry gas volumetric depletion, the material balance equation may be expressed as

1224

q Gp ¼ 1 qi CGIP

1220

1222



 ðB-23Þ

1226

By substituting the GPD definition (Eq. (B-7)) and Eq. (B-21) into Eq. (B-23), we have

1229

q G ¼ 1  2r q 2 PD qi rinf ;D  1

1230

or

1233

q GPD ¼ 1  2r q  2 qi C r  t 0:5 þ1 1 D

1234

when r inf >> rw this becomes:

1225

1227

1231

1235

ðB-24Þ

ðB-25Þ

2r q G q ¼ 1  2 PD  qi C r tD

15

ðB-26Þ

1237

The relationship between CGIP of radial flow and actual time can be derived by substituting Eq. (B-22) into Eq. (B-21):

1238



2 pffiffiffiffiffi 1  C r  Xt þ 1  1 CGIP ¼ 2XY

1240

ðB-27Þ

which, for rinf >> rw , simply becomes:

1239

1242 1243

1244

CGIP ¼

C 2r 2Y

t

ðB-28Þ

1246

References

1247

Abou-Kassem, J.H., Mattar, L., Dranchuk, P.M., 1990. Computer calculations of compressibility of natural gas. J. Can. Pet. Technol. 29 (5), 105–108. http:// dx.doi.org/10.2118/90-05-10. Agrawal, R.G., 1979. Real gas pseudo-time – A new function for pressure buildup analysis of Mhf gas wells. SPE Paper 8279 Presented in 54th SPE Annual Technical Conference and Exhibition, Las Vegas, Nevada, 23–26 September. http://dx.doi.org/10.2118/8279-MS. Al-Hussainy, R., Ramey Jr., H.J., Crawford, P.B., 1966. The flow of real gases through porous media. J. Pet. Technol. 18 (5), 624–636. http://dx.doi.org/10.2118/1243A-PA, SPE-1243-A-PA. Anderson, D.M., Mattar, L., 2007. An improved pseudo-time for gas reservoirs with significant transient flow. J. Can. Pet. Technol. 46 (7), 49–54. http://dx.doi.org/ 10.2118/07-07-05. Ayala H, L.F., Ye, P., 2013a. Density-based decline performance analysis of natural gas reservoirs using a universal type curve. J. Energy Resour. Technol. 135 (4), 042701. http://dx.doi.org/10.1115/1.4023867 (10 pages). Ayala H, L.F., Zhang, M., 2013b. Rescaled exponential and density-based decline models: extension to variable rate/pressure-drawdown conditions. J. Can. Pet. Technol. 52 (6), 433–440. http://dx.doi.org/10.2118/168223-PA. Behmanesh, H., Tabatabaie, S.H., Heidari Sureshjani, M., Clarkson, C.R., 2014. Modification of the transient linear flow distance of investigation calculation for use in hydraulic fracture property determination. Paper SPE 168981 Presented at the SPE Unconventional Resources Conference, The Woodlands, Texas, USA, 1–3 April. http://dx.doi.org/10.2118/168981-MS. Chen, C., Raghavan, R., 2013. On the liquid-flow analog to evaluate gas well producing in shales. SPE Res. Eval. Eng. J. 16 (2), 209–215. http://dx.doi.org/ 10.2118/165580-PA, SPE-165580-PA. Clarkson, C.R., 2013. Production data analysis of unconventional gas wells: review of theory and best practices. Int. J. Coal Geol. 109, 101–146. http://dx.doi.org/ 10.1016/j.coal.2012.11.016. Dranchuk, P.M., Abou-Kassem, J.H., 1975. Calculation of Z factors for natural gases using equations of state. J. Can. Pet. Technol. 14 (3), 34–36. http://dx.doi.org/ 10.2118/75-03-03. Duong, A.N., 2011. Rate-decline analysis for fracture-dominated shale reservoirs. SPE Res. Eval. Eng. 14 (3), 377–387. http://dx.doi.org/10.2118/137748-PA, SPE137748-PA. Earlougher Jr., R.C., 1977. Advances in Well Test Analysis, Monograph v. 5. Society of Petroleum Engineers of AIME, New York-Dallas. Ehlig-Economides, C.A., Ramey Jr., H.J., 1981. Pressure buildup for wells produced at a constant pressure. SPE J. 21 (1), 105–114. http://dx.doi.org/10.2118/7985-PA. Energy Resources Conservation Board, 1975. Theory and Practice of the Testing Gas Wells, third ed., ERCB-75-34. Energy Resource Conservation Board, Calgary, Alberta, Canada. Fraim, M.L., Wattenbarger, R.A., 1987. Gas reservoir decline-curve analysis using type curves with real gas pseudo pressure and normalized time. SPE Form. Eval. 2, 671–682. http://dx.doi.org/10.2118/14238-PA, SPE-14238-PA. Ibrahim, M., Wattenbarger, R.A., 2006. Rate dependence of transient linear flow in tight gas wells. J. Can. Pet. Technol. 45 (10), 18–20. http://dx.doi.org/10.2118/ 06-10-TN2. Jacob, C.E., Lohman, S.W., 1952. Nonsteady flow to a well of constant drawdown in an extensive aquifer. Trans. Am. Geophys. Union 33 (August), 559–569. Kuchuk, F.J., 2009. Radius of investigation for reserves estimation from pressure transient well tests. SPE 120515 presented at the SPE Middle East Oil and Gas Show and Conference, 15–18 March, Bahrain. http://dx.doi.org/10.2118/ 120515-MS. Lee, J., 1982. Well testing. SPE Textbook Series, vol. 1. Society of Petroleum Engineers of AIME, New York-Dallas. Lee, A.L., Gonzalez, M.H., Eakin, B.E., 1966. The viscosity of natural gases. J. Pet. Technol. 18 (8), 997–1000. http://dx.doi.org/10.2118/1340-PA, SPE-1340-PA. Nobakht, M., Clarkson, C.R., 2012a. A new analytical method for analyzing linear flow in tight/shale gas reservoirs: constant-flowing-pressure boundary condition. SPE Res. Eval. Eng. 15 (3), 370–384. http://dx.doi.org/10.2118/ 143989-PA. Nobakht, M., Clarkson, C.R., 2012b. A new analytical method for analyzing linear flow in tight/shale gas reservoirs: constant-rate boundary condition. SPE Res. Eval. Eng. 15 (1), 51–59. http://dx.doi.org/10.2118/143990-PA.

1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

JUOGR 55

No. of Pages 16, Model 5G

1 June 2015 16 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342

P. Vardcharragosad, L.F. Ayala / Journal of Unconventional Oil and Gas Resources xxx (2015) xxx–xxx

Nobakht, M., Clarkson, C.R., 2012c. Analysis of production data in shale gas reservoirs: rigorous corrections for fluid and flow properties. J. Nat. Gas Sci. Eng. 8 (March), 85–98. http://dx.doi.org/10.1016/j.jngse.2012.02.002. Nobakht, M., Clarkson, C.R., Kaviani, D., 2012. New and improved methods for performing rate-transient analysis of shale gas reservoirs. SPE Res. Eval. Eng. 15 (3), 335–350. http://dx.doi.org/10.2118/147869-PA. Qanbari, F., Clarkson, C.R., 2013. A new method for production data analysis of tight and shale gas reservoirs during transient linear flow period. J. Nat. Gas Sci. Eng. 14, 55–65. http://dx.doi.org/10.1016/j.jngse.2013.05.005. Qanbari, F., Clarkson, C.R., 2014. Production data analysis of multi-fractured horizontal wells producing from tight oil reservoirs - bounded stimulated reservoir volume. SPE 167727 presented at the SPE/EAGE European Unconventional Resources Conference and Exhibition, Vienna, Austria, 25–27 February. http://dx.doi.org/10.2118/167727-MS. Stehfest, H., 1970. Algorithm 368: numerical inversion of laplace transforms. Commun. ACM. 13, 47–49. Sutton, R.P., 1985. Compressibility factors for high-molecular-weight reservoir gases. Paper SPE 14265 presented at the SPE Annual Technical Conference and Exhibition, Las Vegas, Neveda, 22–25 September. http://dx.doi.org/10.2118/ 14265-MS. Tabatabaie, S.H., Mattar, L., Pooladi-Darvish, M., 2013. Pseudotime calculation in low permeability gas reservoirs. SPE 167185 presented at the SPE Unconventional Resources Conference, Calgary, Alberta, Canada, 5–7 November. http://dx.doi.org/10.2118/167185-MS. Uraiet, A.A., Raghavan, R., 1980. Unsteady flow to a well producing at a constant pressure. J. Pet. Tech. 32 (10), 1803–1812. http://dx.doi.org/10.2118/8768-PA. van Everdingen, A.F., Hurst, W., 1949. The application of the laplace transformation to flow problems in reservoirs. Trans. Am. Inst. Mining Metal. Pet. Eng. 186, 305–324.

Vardcharragosad, P., Ayala H, L.F., 2013. Density-based rate-time production analysis of unconventional gas reservoirs under gas slippage and desorption. Paper SPE 166377 presented at the SPE Annual Technical Conference Exhibition, New Orleans, Louisiana, USA, 30–2 October. http://dx.doi.org/10.2118/166377MS. Vardcharragosad, P., Ayala H, L.F., 2014. Boundary-dominated decline behavior of gas reservoirs under linear flow. SPE J. (in review). Wattenbarger, R.A., El-Banbi, A.H., Villegas, M.E., Maggard, J.B., 1998. Production analysis of linear flow into fracture tight gas wells. SPE Paper 39931 Presented at the 1998 SPE Rocky Mountain Reginal/Low Permeability Reservoirs Symposium and Exhibition held in Denver, Colorado. http://dx.doi.org/10. 2118/39931-MS. Ye, P., Ayala H, L.F., 2012. A density-diffusivity approach for the unsteady state analysis of natural gas reservoirs. J. Nat. Gas Sci. Eng. 7 (July 2012): 22–34. http://dx.doi.org/10.1016/j.jngse.2012.03.004. Zhang, M., 2013. Variable pressure-drop/flow-rate system analysis for gas reservoirs: a density-based approach (MS Thesis). The Pennsylvania State University, University Park, Pennsylvania, August 2013. Zhang, M. and Ayala H, L.F., 2013. Rate-time decline analysis for natural gas wells with significant rock compressibility effects. Paper SPE 166320 presented at the SPE Annual Technical Conference Exhibition, New Orleans, Louisiana, USA, 30–2 October. http://dx.doi.org/10.2118/166320-MS. Zhang, M., Ayala H, L.F., 2014. Gas-rate forecasting in boundary-dominated flow: constant-bottomhole-pressure decline analysis by use of rescaled exponential models. SPE J. 19 (3), 410–417. http://dx.doi.org/10.2118/168217-PA, SPE168217-PA.

Please cite this article in press as: Vardcharragosad, P., Ayala, L.F. Rate-time forecasting of gas reservoirs with significant transient flow: A density-based method. J. Unconventional Oil Gas Resourc. (2015), http://dx.doi.org/10.1016/j.juogr.2015.05.003

1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369