A Mathematical Model Water Movement about

Bottom-Water-Drive Reservoirs
KEITH H. COATS * THE U. OF MICHIGAN
JUNIOR MEMBER AlME ANN ARBOR, MICH.

ABSTRACT horizontal interface between reservoir fluid and

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

aquifer water and with a significant depth of aquifer
This paper presents the development and solution
underlying the reservoir. In these cases, bottom-
of a mathematical model for aquifer water move-
water drive will occur, and a three-dimensional
ment about bottom-water-drive reservoirs. Pressure
model accounting for the pressure gradient and
gradients in the vertical direction due to water flow
water flow in the vertical direction should be
are taken into account. A vertical permeability
employed. This paper treats such a model in detail
equal to a fraction of the horizontal permeability is
- from the description of the model through for-
~lso included in the model. The solution is given
mulation of the governing partial differential equation
In the form of a dimensionless pressure-drop quantity
to solution of the equation and preparation of tables
tabulated as a function of dimensionless time. This
giving dimensionless pressure drop as a function of
quantity can be used in given equations to compute
dimensionless time. The model rigorously accounts
reservoir pressure from a known water-influx rate,
for the practical case of a vertical permeability
to predict water-influx rate (or cumulative amount>
equal to some fraction of the horizontal permeability.
from a reservoir-pressure schedule or to predict
The pressure-drop values can be used in given
gas reservoir pressure and pore-volume performance
equations to predict reservoir pressure from a known
from a given gas-in-place schedule. The model is
water-influx rate or to predict water-influx rate (or
applied in example problems to gas-storage reser-
cumulative amount) when the reservoir pressure is
voirs, and the difference between reservoir per-
known.
for.mances predicted by the thick sand model of
thzs paper and the horizontal, radial-flow model is The inclusion of gravity in this analysis is actually
shown to be appreciable. trivial since gravity has virtually no effect on the
flow of a homogeneous, slightly compressible fluid
INTRODUCTION in a fixed-boundary system subject to the boundary
conditions imposed in this study. Thus, if the
The calculation of aquifer water movement into or
acceleration of gravity is set equal to zero in the
?ut ?f oil and p-s reservoirs situated on aquifers
following equations, the final result is unchanged.
IS Important In pressure maintenance studies
The pressure distribution is altered by inclusion of
material-balance and well-flooding calculations:
gravity in the analysis, but only by the time-constant
In gas storage operations, a knowledge of the water
hydrostatic head.
movement is especially important in predicting
The equations developed are applied in an example
pressure and pore-volume behavior. Throughout this
case study to predict the pressure and pore-volume
paper the term "pore volume" denotes volume
behavior of a gas storage reservoir. The prediction
occupied by the reservoir fluid, while the term "flow
of reservoir performance based on the bottom-water-
model" refers to the idealized or mathematical
drive model is shown to differ significantly from
representation of water flow in the reservoir-aquifer
that based on van Everdingen and Hurst's horizontal-
system.
flow model.
. The prediction of water movement requires selec-
tIOn of ~ flow model for the reservoir-aquifer system. DESCRIPTION OF FLOW MODEL
A phYSIcally reasonable flow model treated in detail
to dat.e is the radial-flow model considered by van The edge-water-drive flow model treated by van
Everdmgen and Hurst.1 In many cases the reservoir Everdingen and Hurst 1 is shown in Fig. la. The
is situated on top of the aquifer with a continuous aquifer thickness h is small in relation to reservoir
radius rb. water invades or recedes from the field
Original manuscript received in Society of Petroleum Engineers at the latter's edges, and only horizontal radial
office july 27. 1961. Revised manuscript received jan 17
1962. Paper presented at 36th Annual Fall Meeting of SPE' flow is considered as shown in Fig. lb. The bottom-
0'1. 8-11. 1961. in Dallas. • water-drive reservoir-aquifer system treated herein
References given at end of paper.
is sketched in Fig. 2a and 2b. Here the aquifer
*Presently associated with jersey Production Research Co
Tulsa. .,
44 /'1/,VSOCIETY OF PETROLEUM ENGINEERS JOURNAL
thickness h is appreciable in relation to Tb. water and
flows into and out of the reservoir across a roughly
horizontal reservoir fluid-water interface, and flow P(TD. y, tD) = Pi(Y) - P(TD. Y, t D), (5)
components in the vertical direction exist. The
aquifer is considered as a right circular cylinder allows Eq. 1 to be written
of height h and exterior radius Te , with upper and
lower faces impermeable except for that portion
1
(T <,.,,) of the upper face intersected by the reservoir. + (6)
The aquifer formation is considered to have constant, TD
but unequal, permeabilities in the horizontal and
vertical directions. The case of an average vertical The pressure Pi (Y) is the initial aquifer pressure
permeability equal to a fraction of the average which is assumed to be constant except for the
horizontal permeability is a practical one in aquifers vertical variation due to gravity. That is,
riddled with thin, discontinuous shale streaks. This
fraction may be taken as 1.0, of course, for applica- g
tions of this thick sand model to aquifers considered Pi (Y) = Po + P gc Z • • • • • • • •• (7)
homogeneous.

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

where Po is a constant equal to the initial aquifer
MATHEMATICAL CONSIDERATIONS
pressure at the horizontal plane of the reservoir,
The partial differential equation governing un- Z = 0 (see Fig. 2b).
steady-state flow of a slightly compressible fluid in Eq. 6 is solved here for the case of an infinite
the geometry shown in Fig. 2b appears as 2 aquifer (i.e., Te = 00) and for the "constant rate
case" 1 wherein the rate of water flow across the
reservoir-aquifer interface (z = 0, T < Tb) is specified.
1 ap p.¢c The basic solution is obtained for a constant rate
+ - - + kR (1)
T a,. k of water influx, while the general solution for an
arbitrary time-dependent rate is obtained by applica-
where k R is the ratio of vertical effective permea- tion of Duhamel's superposition principle 3 to the
bility k V to horizontal permeability k. Definition of basic solution. The velocity of water flow vertically
the new variables, into the reservoir is given by Darcy's law as
= TITh>'
u= ~ ~)
TD • • • • • • • • • • • • • • (2)
(tzP
P. ~
- p gc z = 0 .....
(8)
........... (3)

If this velocity is considered constant over the

tD =ktlp.¢ CTb2*, • • • • • • • • • •• (4) area of the reservoir (0 < T< Tb), then the volumetric
rate of water influx, e w is given by
..,hen the units given in the Nomenclature are used, tD "" 6.33kt/llcP cr b2 ,

RESERVOIR
(OIL OR GAS)\,--+_~ 1TTb
2 kv (dP
11 az -
g )
P gc z = 0' • • • (9)

\. which is equivalent to

;i'
/

AQUIFER~
FIG. 1a - EDGE-WATER-DRIVE FLOW SYSTEM.

FIG. 2a - BOfTOM-WATER-DRIVE FLOW SYSTEM.

\RESERVOIR .--r--~/ RESERVOIR

\ §- -

Z=h~~-=~---+-=~--=~--~
1/ '
Z=Q ~--""-Z---..JA--QUIFER-I-~-~l..------'
I -
----
FIG. 1 b - IDEALIZED FLOW MODEL FOR EDGE- FIG. 2b - IDEALIZED FLOW MODEL FOR BOfTOM-
WATER-DRIVE SYSTEM. WATER-DRIVE SYSTEM.

MARCH, 1962 45
 TTTb k YkR. ~ap)
ew =
P. aYy=o ...... (10)
Since this solution varies with radius TD, the
Thus, the boundary condition at z = 0 for the basic question arises as to what value of TD between 0
(constant rate) solution is and 1 should be chosen for numerical evaluation
of P. Rather than choosing a single value of TD.
the solution Eq. 16 was integrated over the radius
to obtain the "areal mean" dimensionless reservoir-
pressure drop.

(11)

where ew is considered constant. The initial con-

dition and other boundary conditions corresponding or
to impermeable lower boundaries, finiteness of
P at T = 0 and equilibrium at T = 00 are 00 h 2 (x) e- x2tv
P(tD) 2 fo [eoth Mx-
x2 Mx
P (TD, y, 0) = 0, (12)

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

m = -(x 2 + a2m)tv~
f.~PJ
2x e
. • (18)
00

=0 •.•...•.. (13) .1 2 dx •
x2
\ YJy = hlrb y;;:;-
M m = 1 +a,n

P was numerically integrated on an IBM 704

0, . • • . . . • • • (14) digital computer for several values of the parameter
M, and the results are listed in Table 1. Although
and the integrations were carried out to a dimensionless
time of 1,600, the results can be approximated
lim quite closely * for tD > 10 by
P(TD' y, t D) = finite . . (15)
TD~ 0
1
The assumption of a constant influx velocity u P = A + 4M In tD • • • • • • • • • • (19)
(Eq. 8) over the reservoir area T < Tb is actually
open to little objection since, if u were considered where A is a constant dependent upon M as shown
a function of radius TD, the solution obtained would in Table 2 and Fig. 3. Values of P, for tD < 10 and
show the same time dependence but would differ by for values of M not listed in Table 1, can be found
a multiplicative constant. Since the solution (Eq. by interpolating between curves of P vs tD as shown
16) contains the multiplicative constant P./TTTb k {kR, in Fig. 4. Cross plots ofP vs M with tD as parameter
the appearance of another constant is immaterial would also serve this purpose.
for two reasons: (1) the factors kip., Tb and VkR
are not generally known exactly; and (2) the con- WORKING EQUATIONS
stant P./TTTbk...,;'kR will be chosen in a practical case
Combination of Eqs. 5, 7, 16 and 17 yields
by matching predicted pressures with available
field data.
The solution to Eq. 6 for the conditions of Eqs. . •• (20)
11 through 15 is derived in the Appendix and appears
as which gives reservoir pressure P as a function of
"" hex) time for a constant rate of water infl ux e W" The
E fo x [coth Mx constant TT has been absorbed into other unit con-
version factors to give the constant .0502; units of
the variables present are given in the Nomenclature.
e - x2tv 2x m = 00 :-(x 2 + a 2m )tv] Application of Duhamel's superposition principle 3
.1 e 2 J O(TDX) dx to the constant-rate-case solution (Eq. 20) gives the
Mx M m=l x2 + a m
"variable-rate-case" solution
(16)
i=j-l
where y has been set equal to zero, .0502 f1
Pj = Po - k TbVJi:R I tie; Pj-i • • (21)
i=O

E
Since !'!.ej is defined as ei + 1 - ei and ej as the
*The error between values of ji tabulate<! in Table 1 and calculated from Eq. 19 is
.E is .0502 eu'p./ r"k {kii when units given in the Nomenclature are used. less than 1 per cent at tv = 10 and decreases rapidly with increasing time.

46 SOCIETY OF PETROLEUM ENGINEERS JOllRNAL

 TABLE 1 - DIMENSIONLESS PRESSURE DROP VS
DIMENSIONLESS TIME FOR THICK SAND MODEL TABLE 2 - DEPENDENCE OF A UPON M
P Values for Different Values of M M A
tD M=.05 M=.1 M=.3 -'1=.5 M=.7 M=.9 M=I.0
.05 6.5644
.1 1.433 .741 .329 .287 .285 .285 .285 .10 3.3065
.2 2.675 1.362 .536 .412 .383 .377 .377 .30 1.1846
.3 3.649 1.849 .699 .509 .454 .436 .436 .sO 0.8010
.4 4.464 2.257 .835 .591 .512 .483 .483 .70 0.6622
.5 5.166 2.608 .952 .661 .562 .522 .521 .90 0.6000
1 7.693 3.871 1.373 .914 .743 .663 .648
1.00 0.5911
2 10.590 5.319 1.856 1.204 .950 .824 .792
3 12.442 6.245 2.164 1.389 1.082 .927 .885
4 13.789 6.919 2.389 1.523 1.178 1.001 .952 average rate of water influx during the time increment
5 14.848 7.448 2.565 1.629 1.254 1.060 1.005
6 15.722 7.885 2.711 1.717 1.316 1.109 1.049 from (i - I)At to iAt, i?i is simply
7 16.465 8.257 2.835 1.791 1.369 1.150 1.086
8 17.111 8.580 2.942 1.856 1.416 1.186 1.118 e; = (Vi -1 - Vi)/ At • • • . • (22)
9 17.684 8.866 3.038 1.913 1.456 1.218 1.147
10 18.198 9.123 3.123 1.964 1.493 1.246 1.173
12 19.089 9.569 3.272 2.053 1.557 1.296 1.217 and

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

14 19.846 9.947 3.398 2.129 1.611 1.338 1.255
16 20.503 10.276 3.508 2.195 1.658 1.374 1.288 A"ei = e;+1 - ei = (2Vi - V;-1 - Vi + 1 )/M
18 21.083 10.566 3.604 2.253 1.699 1.407 1.317
20 21.603 10.826 3.691 2.305 1.736 1.435 1.343 == A Vi/At.
24 22.505 11.277 3.841 2.395 1.801 1.486 1.388
28 23.269 11.659 3.969 2.471 1.855 1.528 1.426
32 23.931 11.990 4.079 2.538 1.903 1.565 1.459 Therefore, Eq. 21 can be written as
36 24.516 12.282 4.177
4.264
2.596
2.648
1.944
1.982
1.597 1.489
1.626 1.515 .0502 fl i=j-l
40 25.039 12.544 I
45 25.625 12.837 4.361 2.707 2.024 1.659 1.544 Pj = Po - rb k YkRAt ;=0
50 26.149 13.099 4.449 2.759 2.061 1.688 1.570
13.336 4.528 2.807 2.095 1.714 1.594
(23)
55 26.623
60 27.057 13.553 4.600 2.850 2.126 1.738 1.616 and gives reservoir pressure P as a function of
70 27.824 13.937 4.728 2.927 2.181 1.781 1.654 time for an arbitrary, time-variant water-influx
80 28.490 14.269 4.839 2.994 2.228 1.818 1.687 rate or, equivalently, reservoir pore volume variation.
90 29.077 14.563 4.937 3.052 2.270 1.851 1.717
100 29.603 14.826 5.024 3.105 2.308 1.880 1.743 The solution for the case of constant reservoir
120 30.512 15.281 5.176 3.196 2.373 1.930 1.788 pressure; while not obtained directly here, can be
140 31.282 15.665 5.304 3.272 2.428 1.973 1.827 approximated by re-arranging Eq. 23 as
160 31.948 15.998 5.414 3.339 2.475 2.010 1.860
180 32.536 16.292 5.513 3.398 2.517 2.043 1.890
200 33.062 15.556 5.601 3.451 2.555 2.072 1.916 . = .05 <£"l( [i=[;2 -p
- up
A A 11
220 33.538 16.794 5.680 3.498 2.589 2.099 1.940 Po - Pj
. 1'[,
k./
vkRo.t
A ~
;=0
o.Vj j-i
240 33.973 17.011 5.753 3.542 2.620 2.123 1.961
260
280
300
34.373
34.743
35.088
17.211
17.396
17.568
5.819
5.881
5.938
3.582
3.619
3.653
2.649
2.675
2.700
2.145 1.981
2.165 2.000
2.185 2.017
+ (2Vj _ 1 - Vj-2)i\ - Vj PI]'
330 35.564 17.806 6.018 3.701 2.734 2.211 2.041
360 35.999 18.024 6.090 3.744 2.765 2.235 2.063 and
390 36.399 18.224 6.157 3.784 2.793 2.257 2.083
420 36.769 18.409 6.219 3.821 2.820 2.278 2.101
450 37.114 18.581 6.276 3.856 2.844 2.297 2.119
480 37.436 18.742 6.330 3.888 2.867 2.315 2.135
510 37.739 18.894 6.380 3.918 2.889 2.332 2.150
540 38.025 19.037 6.428 3.947 2.909 2.348 2.164
570 38.295 19.172 6.473 3.974 2.929 2.363 2.178
600 38.551 19.300 6.515 4.000 2.947 2.377 2.190 20
650 38.951 19.500 6.582 4.040 2.976 2.399 2.210
700 39.322 19.685 6.644 4.077 3.002 2.420 2.229 10
750 39.667 19.858 6.702 4.111 3.027 2.439 2.246
6.0 "
4.0 \.
800 39.989 20.019 6.755 4.143 3.050 2.457 2.262
850 40.292 20.171 6.806 4.174 3.071 2.474 2.277

"'-
900 40.578 20.313 6.854 4.202 3.092 2.490 2.292
950 40.840 20.449 6.899 4.229 3.111 2.505 2.305 2.0 - --
1000 41.105 20.577 6.941 4.255 3.129 2.519 2.318
1050 41.349 20.699 6.982 4.279 3.147 2.532 2.330 A 1.0
1100 41.581 20.815 7.021 4.303 3.163 2.545 2.342
1150 41.803 20.926 7.058 4.325 3.179 2.558 2.353 .6
1200 42.016 21.032 7.093 4.346 3.194 2.570 ·2.364 .4
1250 42.220 21.134 7.127 4.367 3.209 2.581 2.374
1300 42.416 21.233 7.160 4.386 3.223 2.592 2.384 .2
1350 42.605 21.327 7.191 4.405 3.237 2.602 2.393
l400 42.787 21.418 7.222 4.423 3.250 2.612 2.402 .1
1450 42.962 21.505 7.251 4.441 3.262 2.622 2.411 o 2 4 6 8 10
1500 43.132 21.590 7.279 4.458 3.274 2.631 2.419 M
1550 43.296 21.672 7.306 4.474 3.286 2.641 2.428
1600 43.454 21.752 7.33291 4.490 3.300 2.649 2.436 FIG. 3 - PLOT OF A VS M FOR TlflCK SAND MODEL.

MARCH, 1962 47
 40r-----.-----.-----.---~.-----. Interpolation in Table 1 at M = .3 then yields the
first three columns of Table 3. From the data given,

30~--~-----+----~----~--~
.0502(-80,000)(1)
-7.1
3000 (.310) ..{.IT
Ia... 20 1-------+----
where e w is negative because the bubble is being
grown and water is moving away from the reservoir.
Eq. 20 now becomes

P = Po + 7.1 P = 1080 + 7.1 P,

which allows calculation of the last column in Table

tD 3. The required gas-injection schedule could be
FIG. 4 - 15 AS A FUNCTION OF tn AND M.
calculated from the gas equation of state =PV / n zRT
where P is reservoir pressure, V is reservoir pore
volume, n is gas in place and z is the compressibility
(24) factor (the reservoir pore volume Vo at "zero time",

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

when uniform pressure of 1,080 exists, would have
to be known).
Eq. 24 gives the reservoir pore volume V as a Eqs. 20, 23 and 24 are the basic "working equa-
function of time when the reservoir pressure P is tions" allowing calculation of reservoir pressure
held at a constant value less than Po by IIp. The or volume from knowledge of the water-influx rate
cumulative volume of water influx into the reservoir or reservoir pressure. In general, however, neither
is, of course, the influx rate nor reservoir pressure is known in
(25) advance; rather, a fluid-in-place (oil or gas) schedule
is known or specified and an estimate of the reservoir
EXAMPLE PROBLEM 1 performance (pressure and/or pore volume vs time)
is desired. In this case, a relationship is needed
A gas storage reservoir has been created by in-
between fluid in place, pore volume V and reservoir
jection of gas into an aquifer formation. The field
pressure p, i.e., a material-balance equation. In the
has been grown to a radius Tb of 3,000 ft and has
case of gas reservoirs, this material balance is
been shut in for a period of time sufficient to allow
exceptionally simple and will be used in conjunction
the reservoir and aquifer pressure to reach an approx-
with Eq. 23 to yield a pressure-explicit equation for
imately uniform value of 1,080 psia. Estimate the
use in solution of the problem just stated. For a
reservoir pressure as a function of time which must
gas reservoir,
be maintained to grow the gas bubble at a constant
rate of 80 Mcf of pore volume/D. The aquifer for-
mation is 550-ft thick, core data indicate a permea- V·1 = n'RT(.z.)*
1 P j' (26)
bility ratio k R of .37, and water-pumping tests
indicate an effective aquifer horizontal permeability
of .310 darcies. Other available data are ¢ = .17, and since z can be represented as a linear function
c = 7 x 10-6 l/psi, f.1. = 1 cpo of pressure over normal operating pressure ranges,

SOLUTION Zj = a + bPi
Since the rate of water movement is specified to
be constant and bottom-water drive exists, the and
constant-rate Eq. 20 will be employed. The value
Mis
M = hlTbVkR = 550/3000..[:37 = 0.3,
*The assumption is implied here that the reservoir pressure p is essentially uniform
so that PI in Eqs. 26 and 23 are in the same pressure.
so that P values corresponding to M = .3 will be
read from Table 1. If pressures are calculated at TABLE 3
30-day intervals, then at the end of i 30-day periods, Time Reservoir
(months) 'n P Pressure p (psi)
6.33k /).t
tD =i f.1. ¢ CTb2
0 0 0 1080
5.5 2.638 1098.7
2 11 3.197 1102.7
6.33 (.310) (30) 3.532 1105.1
3 16.5
1(. 17)(7xlO-<i)(3000)2 4 22 3.766 1106.8
5 27.5 3.953 1108
5.5 i. 6 33.0 4.103 1109.1

c18 ~OC.IETY OF PETROLEUM ENGINEERS JO{IRNAL

Eq. 23 can be written from various tests and sources; Vo = 4.4 X 108 cu ft,
Po = 466 psia, Tb = 1,880 ft, k = .0296 darcies,
p. = 1.2 cp, RT = 5,570 psia cu ft/lb mole, cp = 0.2,
Pi = Po - T~Z5~t1t [i1:2 t1V i P i- i c = 7 x 10-6 l/psi, h = 565 ft, k R = 1.0, a = .998,
and b = - .00016. Estimate the pressure and pore-

+ (2Vj _ 1 - Vj - 2 ) PI - Vi pJ, .. (28)

volume behavior of the reservoir using the thick
sand model, and compare this behavior to that
predicted by the horizontal, radial-flow model.
and elimination of Vj between Eq. 27 and Eq. 28 SOLUTION
yields a quadratic in Pi which gives
The production schedule plotted in Fig. 5 consists
Pi = C j + D j ••••••••••••• (29) of monthly values of cumulative gas produced over
a period of 254 months. A time increment t1t of one
where month, or 30.4 days, will be chosen therefore. From
the data and Eq. 27,

Cj * = tL fpo + FPIRTbni - F (i~~-2 t1Vi Pi- i

\ ,=0 no = Vo/RT (b + a/po)

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

= 4.4 x 10 8 /5570(-.00016 + .998/466)
+ (2Vj _ 1 - Vi - 2 )J ·
)P1 = 39.9 x 10 6 lb mole.

The gas in place at time jt1t then

Di = CJ + FPl RTnja,
IS

nj = no - Gpj Phase /RTbase

and
= no - G pj [15.025/10.73 (520)]

where Gpj is the cumulative standard cubic feet of

gas produced at the end of j months. An IBM 704
Eq. 29 allows direct calculation of gas reservoir
computer program was written to accept the gas-
pressure at successi;ve times of t1t, ·2&, 3t1t, . . . ,
production schedule, the afore-mentioned data,
for any given gas-in-place schedule Tlj. Coats, Tek
and a table of P values vs dimensionless time tD
and Katz 4 derived an equation similar to Eq. 29
for M = h/Tb ...jki = 565/1880 (1) = .3 to calculate
from the constant-rate-case solution given by van
the gas-in-place nj; and solve Eq. 29 for the predicted
Everdingen and Hurst! for horizontal, radial aquifer
pressure Pi" These pressure and the corresponding
water flow.
pore volumes, calculated from Eq. 27, are plotted
as the solid curves in Figs. 6 and 7.
EXAMPLE PROBLEM 2 The variable-pressure-case solution for the
A newly discovered gas reservoir is to be con- horizontal, radial-flow model1 IS
verted to storage use after a period of production.
The projected schedule of cumulative gas production
is plotted in Fig. 5. The following data are available

*For j '"" 1,
i ... - l
1: /lV, Pj -
~

i == O~
where !'!.pi =Pi-l - Pi+l' t1po = Po - hand Qi-i
, ... 0
is the dimensionless influx quantity Q at tD= (jzz)
t1t D tabulated by van Everdingen and Hurst. Com-
bination of this equation with Figs. 25 and 27 gives
12 a pressure-explicit equation 4 similar to Eq. 29.
Solution of this equation for the same data previously
8 --
J
I~
given, the same gas-in-place schedule and for an
LL.
U 4 L infinite (in radial extent) aquifer gave the dashed
u
V curves shown in Figs. 6 and 7.
C/)

~
~
0 DISCUSSION AND CONCLUSIONS
0-
4 --
~{ Eq. 16 or Eq. 20 is, to the author's knowledge,
0 r the only solution available to the diffusivity equation
8 governing aquifer water movement about a bottom-
water-drive reservoir. While this solution is valid
for radially infinite aquifers, other solutions can
40 80 120 160 200 240 280 be obtained by use of finite Hankel transforms
TIME - MONTHS (see Appendix) for aquifers of various degrees of
FIG. 5 - CUMULATIVE GAS PRODUCTION AS A finiteness. Example Problem 2 shows that significant
FUNCTION OF TIME. differences may arise between field performances

MARCH, 1962 49
predicted by the thick sand and horizontal radial- gratefully acknowledges the permission of R. C. F.
flow models. Figs. 6 and 7 show quantitatively the Bartels, director of the U. of Michigan Computing
differences between reservoir pressures and pore Center, to carry out calculations on the University's
volumes calculated from the two models. IBM 704 computer, and the American Gas Assn.
In recent years increasing emphasis has been Project No. 31 for its financial assistance.
placed on the "resistance-curve" technique as
opposed to the flow-model approach. The latter NOMENCLATURE
approach by necessity involves various idealizations
a = constant in equation z = a + bp, dimen-
pertinent to reservoir and aquifer geometry and sionless
aquifer homogeneity. The objection thus arises
that most practical cases violate significantly one b = constant in equation z = a + bp, Ilpsi
or more of these idealizations. The resistance- c = compressibility of aquifer water and for-
curve method meets this objection by requiring no mation, l/psi
assumptions concerning aquifer geometry and homo- E = f!.w IlI1T7'b k VfR
geneity but, rather, involves the d~ermination of a e w = rate of water influx, cu ft/day
resistance-curve analogous to the P vs tD curve of ei ,. average rate of water influx from time
this paper or the Qt vs tD function of van Everdingen (i -l)~t to i~t, cu ftl day
and Hurst l directly from field data. The paper by
~ei = ei+l - ei
Hutchinson and SikoraS is an example of this method.

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

It is the author's opinion that the resistance- gc = conversion constant, 32.17 ft-lb massl
curve method will, in the long run, replace the Ib force-second 2
model approach with resultant increased ease h = aquifer thickness, ft
and accuracy in calculation of water movement. /0. h Bessel functions of first kind, of order 0
=
However, this method has been of little value to 1, respectively
date in studies of edge-water-drive gas fields
k = aquifer formation permeability in hori-
carried out at the U. of Michigan, and requires
zontal direction, darcies
development far beyond its present state. Various
degrees of success in matching field behavior have k R = ratio of vertical-to-horizontal aquifer
been achieved by use of the horiZDntal radial-flow permea bi lity
model, reflecting the various degrees to which the kv = aquifer formation permeability in vertical
edge-water-drive fields studied satisfy the idealiza- direction, darcies
tions involved in that model. At the present time, M = parameter, hlTb{kR
flow models provide a useful tool in prediction of n i = gas in place in reservoir at time j~t, Ib
reservoir performance and, in the future, should mole
remain useful as standards for comparison to p = pressure, psia
resistance curves obtained from field data by a
reliable method. For example, conclusions pertinent
Pi(Y) = initial aquifer pressure, psia
to aquifer geometry andlor extrapolation (in time) Pi = reservoir pressure at time i~t, psia
of resistance curves might be accomplished by Po = initial aquifer (and reservoir) pressure
comparing these curves to those corresponding to at reservoir depth, psia
various models. P =: dimensionless pressure-drop function

ACKNOWLEDGMENT 15; = value of P at tD = i~tD

Suggestions and incentive from D. L. Katz, ~P i = Pi - 1 - Pi + 1

chairman of the Dept. of Chemical and Metallurgical T = radius, ft
Engineering, greatly aided this study. The author
1.0~~----r----r----~--,---~~--,

« 800
V') :>.2-iIRADIAlIFlOW ~ODEl .96 ~.-+­
700 - -
ror-~
a.. ---- - -

w ~ 0 .92
~ 600 --- -
.~
V')
1< 0 V .88
V')
~ 500 c---- - -
nO
Vo .84 I
---r-
a.. d ,,
~
-~+-
~'~ -+
o 400 - - - - r------c 1---
.80 -
«
w oOne "'-- RADIAL FLOW
I 300 -- .76 MODEL -
.....
..... ~I THI~K SA~D MO?El
w 200 .72
~ 0 40 80 120 160 200 240 280 0 40 80 120 160 200 240 280
TIME - MONTHS TIME - MONTHS
FIG. 6 - WELLHEAD PRESSURE AS A FUNCTION OF FIG. 7 -- RESERVOIR-PORE-VOLUME RATIO AS A
TIME. FUNCTION OF TIME.

50 SOCIETY OF PETROLEUM ENGINEERS JOURNAL

 rv = dimensionless radius, rlrb so that multiplication of both sides of Eq. 6 by
rb = radius of reservoir, ft rv] 0 (xTD) dT V and integration from zero to 00 yields
r e = exterior radius of aq uifer, ft
R = gas constant, 10.73 psia-cu ft/lb mole-'R cPu 2 au. (32)
ay2 -x U = atv •••••
t = time, days
IJ.t = time increment, days Solution of Eq. 32 by separation of variables,
tv = dimensionless time, 6.33 ktlp.¢CTb2 U(x, y, tv) = Y(y) (t W, yields e
T = reservoir temperature, oR
u .. velocity of aquifer water flow, ft/day U = A cosh xy + B sinh xy
V = reservoir pore volume, cu ft
2
+ e- 'A t D (C cos y.~ - x2 Y
Vi = reservoir pore volume at time i/);.t; cu ft, + D sin .J'A2 _x 2 y) (33)
V-I = V 0
= 2Vj - Vi-l - Vi+l; /);.Vo = Vo - VI
/);.Vi In order that Eq. 11 be satisfied,
Vo = initial reservoir pore volume, cu ft
f(x)

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

We = cumulative
water influx, cu ft
y .. dimensionless vertical distance, zlrbykR
z = vertical distance co-ordinate, ft, or gas where f(x) is the Hankel transform of f(TD), so that
compressibility factor, dimensionless
_'A 2 t
p = density, lb mass/cu ft Bx+Dv'A2-x 2 e D=f(x)
¢ = aquifer formation porosity
00
am = mTT/M = fo TDf(T W] O(xTD)dT V
p. = aquifer water viscosity, cp
1
REFERENCES =-f E] O(XTV)T VdT D = - E x
. . (34)
1. van Everdingen, A. F. and Hurst, W.: "The Appli- o
cation of the LaPlace Transformation to Flow
Problems in Reservoirs", Trans., AIME (1949) Vol. and D must equal zero and B must be - Eft (x)/x 2
186, 305. if Eq. 34 is to hold for all time. Thus, U is now
2. Muskat, M.: The Flow 0/ Homogeneous Fluids given by
Through Porous Media, J. W. Edwards, Inc., Ann
Arbor, Mich. (1946). ElI (x)
3. Churchill, R. V.: Modern Operational Mathematics
U = A cosh xy sinh xy
in Engineering, McGraw-Hill Publishing Co., Inc.,
N. Y.
4. Coats, K. H., Tek, M. R. and Katz, D. L.: "Method
for Predicting the Behavior of Mutually Interfering
Gas Reservoirs Adjacent to a Common Aquifer", In order that Eq. 13 be satisfied,
Trans., AIME (1959) Vol. 216, 247.
5. Hutchinson, T. S. and Sikora, V. J.: "A Generalized au o
Water-Drive Analysis", Trans., AIME (1959) Vol.
216, 169. ay y=M
6. Sneddon, I. N.: Fourier Trans/orms, McGraw-Hill ElI (x) ,-=---=-
Publishing Co., Inc., N. Y. (1951). = A sinhMx -.::..~..:... cosh Mx -C J'A2 _x 2
x
APPENDIX 2
. sin (-J'A 2 -x 2 M) e -'A t D •
The solution to Eq. 6 for conditions of Eq s. 11
through 15 is obtained by use of the infinite Hankel which requires
transform 6 defined by
00
ElI (x)
A x coth Mx • • . • • • • . . • (35)
2

cPP
Following Sneddon,6 the Hankel transform of --2 + mTT/M == am. m = 0, 1, 2, ...
iJrv
. • (36)

Thus, U is now give~ by

f
00
TV - - 2(a
P +-1-
2
ap) ]0 (xTD)dTD =-x 2 U(X,y,tD),
U
Eft (x) {COSh [x(M -y)]}
0. iJrv TD iJrD x2 sinh Mx
• • • • 0 • 0 • • • • • • • • • •• (31)

MARCH, 1962 51
 111=00

+ "" C mcosamYe -,\2 mD

.. t •••• (37)
m=O
where A'; = x 2 + a'; and the summation is imposed
in order to satisfy the initial condition of Eq. 12.
This initial condition is U(x, Y, 0) ~ 0 or

111=00
Ell (x) cosh [x(M=y)] Since the Hankel inversion integral is 6
I C m cos amY
m=O x2 sinh Mx
00
. (38) P(TD, y, tD) = fox U(x, y, tD)J 0 (TDx)dx, • • (40)
The terms cos a~ form an orthogonal set over the
interval (0, M), so that multiplication of both sides the final solution for P is
of Eq. 38 by cos (a~)dy and integration from 0
to M yields 00 It (x) {COSh [x(M -y)]
P(TD' y, tD) = E io- - X
.h
SID
M
X Mx
2Eh (x)
2 2

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

and
2x _e_-_<..,.,,_+_a;:<"m_)_t_D cos amY}
M x2-+a 2
m

which becomes identical to Eq. 16 when y is set

The final solution for U now appears. equal to zero. ***

52 SOCIETY OF PETROLEUM ENGINEERS JOURNAL

 Further Discussion of Paper Published in
Society 0/ Petroleum Engineers Journal, March, 1962

A Mathematical Model for Water Movement about

Bottom-Water-Drive Reservoirs
KEITH H. COATS JERSEY PRODUCT/ON RESEARCH CO.
JUNIOR MEMBER AIME TULSA, OKLA.

(Published on Page 44)

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

DISCUSSION
J. E. WARREN
MEMBER AIME I GULF RESEARCH & DEVELOPMENT CO.
PITTSBURGH, PA.

The mathematical problem considered by the with the basic assumption made by Hartsock and
author l can be given another physical interpretation Warren 2 in a prior steady-state study of fractured
which is of some practical significance. The alter- systems. Transient S*-values based on Coats' re-
native physical problem involves the approximate sults and steady-state values computed by the
behavior of a single well which is producing at a method of Hartsock and Warren are compared in
constant rate through an axially symmetric, hori- Fig. D-l; the small, systematic deviations from the
zontal fracture of infinite flow capacity which is ideal correlation are the result of the nonuniform
located at the center, top or bottom of a uniformly distribution of the flux over the fracture surface.
thick, horizontal, homogeneous, anisotropic reser- Because of the mathematical equivalence of the
voir of infinite lateral extent which contains a two physical problems, the following supplementary
single, slightly compressible fluid. Using Coats' conclusions can be drawn.
asymptotic results, the wellbore pressure in the l. By including the pseudo-skin resistance S*,
fractured system is given by the following equation. the performance of a reservoir with a bottom-water

Pw = Pl - 162.5 q J1. B {IO/.00633 k t\ FOR ALL STEADY-STATE GASES, FLOW

k h \¢ /1 C TW2 / CAPACITY OF FRACTURE. 125,000 kh
FRACTURE IN CENTER OF FORMATION
+ .351 + .87 S*}; t.? 1580 ¢ /1 C Tj/k
6.5
. . . . . . . . . . . . . . . . . (1)**
S* _ SMA(M/2) -In (TI1Tw) - .4045, center,
where
-t 2MA(M) - In (Tlirw) - .4045, top or
bottom,
~
Z
6.0

III
A(M) geometric constant (numerical values iii 5.5
z
are given in Table 2 of the subject '"
It:
paper), .
~

Ul
I 5.0
M h yklkv /TI'
k permeability in the horizontal direc-
=
tion, and
4.5
k v = permeability in the vertical direction.
It is obvious that Eq. 1 differs from the usual
equation for radial flow into a wellbore by the
geometric parameter S*; this result is in accord
4.5 5.0 5.5 6.0 6.5
lReferences given at end of paper.
-S*(STEADY-STATEI
**Standard AIME nomenclature is used unless otherwise in-
dicated; units are psi, STB/D, cp, md, ft, reservoir bbl/STB, FIG. D-l- COMPARISON OF TRANSIENT AND STEADY-
(psi)-1 and days. STATE RESULTS.

DECEMBER, 1962 367

drive can be asymptotically approximated by utiliz- production tests will indicate erroneously high PI's
ing conventional methods which assume radial in- if the duration of the production period is less than
flux; the values of A(M) given by Coats can be used 1580 ¢ J1 c r 2lk.
f
for limited aquifers if the external radius of the
aquifer is at least four times as large as the radius REFERENCES
of the region initially occupied by oil, r0' and t ?:
1580 ¢ J1 c r021k. 1. Coats, K. H.: "A Mathematical Model for Water Move-
ment About Bottom-Water-Drive Reservoirs", Soc. Pet.
2. If the radius of a high - capacity fracture is Eng. Jour. (March, 1962) 44.
less than one - fourth of the drainage radius, the 2. Hartsock, J. H. and Warren, J. E.: "The Effect of
steady-state values of S* controls the performance Horizontal Hydraulic Fracturing on Well Performance",
of the well for t.? 1580 ¢ J1 c r/lk; post-fracturing Jour. Pet. Tech. (Oct., 1961) 1050.

AUTHOR'S REPLY TO J. E. WARREN

The equivalence pointed out between the bottom- ally gas-storage reservoirs where cycling occurs),
drive reservoir and fractured well problems is of the production (influx) rate may vary considerably
significant interest and is an outstanding example with time which requires that superposition be

Downloaded from http://onepetro.org/spejournal/article-pdf/2/01/44/2152509/spe-160-pa.pdf by guest on 22 July 2022

of consolidation of distinct research results into a applied to Eq. 1. The result of this superposition
more meaningful whole. Only one observation is is that dimensionless pressure-drop values at small
made here. Warren's Conclusion 1 states that times are required and these values cannot be
bottom-drive reservoir performance can be asymp- approximated by the asymptotic expression em-
totically approximated by conventional methods ployed in Eq. 1. Recourse must then be made to an
involving radial influx if the S* factor is included. equation such as Eq. 21 of the subject paper, and
This conclusion applies with more meaning to the P values at dimensionless time less than 10 must
fractured well problem where production rate is be obtained from the tables.
held constant. In the case of a reservoir (especi-
***

368 SOCIETY OF PETROLEUM ENGINEERS JOURNAL