Calculating immobile gas saturations

JOUI

Journal of Petroleum Science and Engineering 15 (1996) 33-43

Calculating immobile gas saturations Eric S. Carlson, Philip W. Johnson SPE, The Uniuersiry of Alabama,

Tuscaloosa, AL 35487-0207,

USA

Received 11 October 1994; accepted 19 September 1995

Abstract We present a simple, fast, analytical expression which can be used to estimate the reservoir pressure as a function of gas saturation, for gas saturations between zero and the critical gas saturation. Use of the relation also makes it extremely easy to assess the local recovery factor as a function of pressure from the bubble point down to the pressure at which the critical gas saturation is reached.

1. Introduction The threshold saturation where free gas begins to flow (during solution gas drive production for example), is called the critical gas saturation. Unfortunately, estimating the pressure at which the reservoir reaches this critical gas saturation has previously involved long, tedious, material balance calculations. Collecting and correlating the data for the material balance method is enough to discourage many operators from attempting to make the estimate. We have developed a method consisting of one explicit equation for determining the pressure at the critical gas saturation, or any lesser gas saturation. Charts are included which facilitate estimation of the required parameters. We hope that the material presented here will promote better reservoir management techniques by oil field operators. 2. Background For an undersaturated oil reservoir (the pressure is above the bubble point), the examination of recovery 0920-4105/96/$15.00 0 1996 Elsevier Science B.V. All rights reserved SSDI 0920-4105(95)00039-9

versus pressure indicates that incrementally more oil can be produced per unit pressure drop at pressures immediately below the bubble point than at pressures above the bubble point. This change in reservoir efficiency is due to the tremendous amount of energy supplied by free gas evolving from the oil. The apparent improvement of reservoir efficiency, however, is short-lived. As soon as the critical gas saturation is reached, free gas begins to flow, and the pressure declines rapidly. Free gas formation caused by pressure depletion usually has detrimental effects on ultimate oil recovery, although it may briefly improve production. An increase in gas saturation causes the oil permeability to drop, while the loss of gas (from the oil) makes the oil shrink and causes the oil viscosity to increase. These changes in the oil and formation properties make it progressively harder for the oil to flow through the reservoir. Under the best of circumstances, these changes reduce oil mobility and consequently the efficiency of any secondary recovery project, even if the project is implemented prior to free gas flow.

34

ES. Cat-km,

P. W. Johnson / Journul

ofPrtrolrum

In solution gas drive and gas cap drive fields it is common for the pressure to drop below the bubble point. As long as little free gas flows, the change in ultimate recovery may be negligible (for instance, if the extra production data leads to a significant improvement in the reservoir characterization, the secondary recovery design can be improved). It is wise to limit free gas production. When free gas flows, the most important source of natural pressure maintenance is lost. The effects of the decrease of oil phase viscosity, increase in oil shrinkage, and reduction of oil permeability become very pronounced. After the inevitable drop in reservoir pressure, high saturations of free gas will remain in the reservoir at low pressure. For all practical purposes, a volume of liquid equal to the volume occupied by the free gas must be injected before a strong response to a secondary recovery program occurs. Therefore, a secondary recovery project will not only be less efficient because of lower oil mobility, but also will give no incremental oil production for a long time after initiation. The deferred response to injection may cause the project to be uneconomical. Operators can gain from a method for estimating the pressure at which the critical gas saturation occurs, and by using the method to minimize the flow of free gas in the reservoir.

3. Traditional

Science und Enginrerirrg

4. Equations

15 (19961 33-43

for analysis

As shown in Appendix A, the pressure drops from the bubble point pressure, the gas saturation can be approximated (in the worst case, with about 6% error) by:

S,(P) =(I

-s-I&

[(:j’_‘Ij

(1)

or:

PF,)

l

=Pb

A,

I/(U .- I)

1-o M

+ 1

(2) I

where 0 I S, I .SgC;S, is the original water saturation; p is the current reservoir pressure; and pb is the bubble point pressure. The term cx is given by:

Bg dR, a=PCp=PFdp

(3)

where cPC designates the compressibility due to a phase change (Ramey, 1964); BP is the gas formation volume factor; B,, is the oil formation volume factor; and R, is the gas in solution. The units must be consistent because 01 is dimensionless. After the gas saturation and the pressure are specified, it is easy to estimate the incremental local oil recovery factor (relative to the initial pressure, p,) using:

methodology RF( P)

The classical method for analyzing reservoir performance prior to the onset of free gas flow is compactly summarized in Appendix B. Details and examples are presented elsewhere (Craft et al., 1991) and will not be repeated here. The method involves the simultaneous solution of two related equations, one that evaluates the oil material balance and one that determines the cumulative gas-oil ratio. The equations are strongly coupled because the oil material balance depends on the cumulative gas-oil ratio, and the ratio depends on the oil production. This interdependence usually leads to iteration. A quick inspection of the equations in Appendix B indicates that a single evaluation of the equations is complicated, and the difficulties are compounded by iteration. The following section summarizes a much simpler method for performance evaluation.

z$[l

-q(p,-p)]

(4)

0

where RF is the incremental recovery factor. This equation includes rock and connate water expansion effects which are negligible compared to gas expansion below the bubble point, but are important above the bubble point pressure. The recovery factor includes the oil recovered above the bubble point. Eq. 4 is rigorously accurate for all regions of the reservoir where water influx has not taken place. We must emphasize that use of Eq. 4 has not been possible without the definite relation between S, and pressure provided by Eq. I.

ES. Carlson, P. W. Johnson/Journal

5. Assumptions

of Petroleum Science and Engineering 15 (1996) 33-43

for analysis

reservoir as a whole does not meet these criteria, however, the equation may still be applied to local regions of the reservoir which adhere to the assumptions.

In deriving Eq. 1, we have assumed that: 1. The free gas is immobile, 2. S, is constant, 3. The pressure decreases continuously with time (dp/dt < 01, 4. The pressure dependent (Y is constant over the interval from pb to p, and 5. Below the bubble point the formation and water compressibilities are negligible, when compared to gas expansion. Other assumptions not relevant to the application of Eq. 1 are used in the derivation presented in Appendix A. The significance of the first assumption is that the equation applies only when the free gas is immobile, or nearly immobile. The equation becomes invalid when there is a mobile gas phase. The equation also does not apply to regions of the reservoir where water influx occurs, or over intervals of time where pressure has remained constant. Even when the

6. Estimation

of (Y

The effects of assuming that (Y is constant are evaluated in Section 7 (Examples). The term 01 is not constant over the pressure interval, since B, , B,, R, and p are pressure dependent. Evaluating the fluid properties at an average pressure between pb and an estimated p(S,) still proves to give very accurate results for p(S,). In practice, if fluid property data are available, the calculation is both simple and accurate. The derivative term, dR,/dp, can be estimated from: dR,_ dp - (Pa”,)

-

=

Rd P2) - Rs( PI) P2 -PI

100 deg. F

:

4 0’ API

0.2

IF,,,,,,,, 0

I,, 2000

35

, , , , , , , , , , , , , , , , , , , , , , , , , , ,, ; , , , , 4000

Pressure Fig. 1. The product of the phase change compressibility

6000

6000

(psia) and the pressure for a 0.6 gravity gas.

36

E.S. Carlson,

P. W. Johnson/

Journal

of Petroleum

(preferably where P,,, should lie somewhere halfway) between p, and p2. If fluid property measurements are not available, Figs. 1-4 present estimates of (Y as it varies with pressure for temperatures of 1OO”F, 175°F and 250”F, and for various oil gravities. Each graph represents a different gas gravity. These curves are based on the Sutton correlations for critical pressure and temperature versus gas gravity (Craft et al., 1991), the Hall-Yarborough method for z-factor calculation (Dake, 1978), Standing’s correlation for gas in solution, and Standing’s correlation for oil formation volume factor. These correlations give reasonable estimates of properties, but in extreme cases the charts may differ from actual field measurements. They should be used in the absence of more reliable data or locally accurate correlations (Glaso and @stein, 1980, or Vazquez and Beggs, 1980). The graphs are very useful for comparative purposes. The curves show that (Y varies between 0.3 and 1.6 in the validity range of the correlations. For

Science and Engineering

15 (1996) 33-43

lower gravity oils the assumption of a constant cx appears to be a good one. The greatest rate of change of the term occurs for a 0.9 gravity gas and a 60” API crude between 500 psia and 2000 psia. We use this most extreme case below to verify the validity of the method.

7. Examples Example 1 - Fluid properties known. Table 1 gives the fluid property behavior of the Canyon Reef reservoir of the Kelley-Snider Field, Texas, as presented in section 5.4 of Craft et al., 1991. Through application of classical reservoir engineering techniques, the gas saturation at a pressure of 1400 psig was estimated to be 10.2%. Using the values in the table near the bubble point this gives:

= 0.5682

1.4 6 0” API A

0.70

--.-..

Fig. 2. The product of the phase change compressibility

250 deg. F

6000

4000

Pressure

Gravity Gas 100 deg. F

6000

(psia) and the pressure for a 0.7 gravity gas

ES. Carlson, P. W. Johnson/ Journal of Petroleum Science und Engineering 15 (1996) 33-43

Using this value for o and assuming tion of 10.2% in Eq. 2 gives: ~(0.102)

a gas satura-

= (1725 + 14.7) [(z!!L)(

1

l-0.56*\ 0.568

order to apply the charts, estimates of the oil gravity, gas gravity and reservoir temperature must be available. To apply Eq. 1, the bubble point pressure also must be known. For this example, assume that the oil gravity is 60” API, temperature is lOO”F, the gas gravity is 0.9 and the bubble point pressure is 1800 psia. From Fig. 4 the value of 01 at the bubble point is estimated to be 0.91 and, if the critical gas saturation is 0.06 and water saturation is 0.25, then a first estimate for the pressure at the critical gas saturation is:

1

-I/

+1

(1 - 0.568)

=

1405 psia

This value is slightly different from the value calculated with the material balance method (1390 psig versus 1400 psig, less than 1% difference). If the properties are evaluated half way between 1400 and 1725, cx = 0.601 and psgC = 1419 psia, or 1404 psig (a difference of 0.3%). Although (Y varies over the pressure and properties range of this problem, our proposed technique nearly reproduced the results of the much more involved material balance method. Furthermore, the material balance method is a numerical approximation, subject to numerical errors. Example 2 Fluid properties unknown. In

0.2

jr

;,,,,I,,,,

0

i ,I 2000

I I,

31

~i(O.06) = 1800[(

1 y-:625 )( ’ s.g;‘)

+ 1]-“(1-o’9’1

= 1649 psia A better value for OLcan be found by using 1725 psia, the average between the bubble point pressure and our first estimate. The revised value for (Y is

I I I I I I I I I I I I I I i I I I I I I I I I I I I I I 4000

Pressure Fig. 3. The product of the phase change compressibility

6000

8000

(psia) and the pressure for a 0.8 gravity gas.

E.S. Carlson, P. W. Johnsnn/Journal

38 Table I Sample fluid property 1991) Pressure

(psig)

:bl,STB)

1725 1700 1600 1500 1400

1.4509 1.4468

1.4303 1.4139

1.3978

data for Examples

;CF/sTB)

0.00141 0.00151 0.00162 0.00174

885 876 842 807 772

for the pressure

at the

P,(OW = ,,,,[i

1

15 (I 996) 33-43

In Appendix A, we outline an accurate sequential implementation of our method. The charts indicate that the fluid specifications for this example represent the conditions which cause the most drastic rate of change for (Y and thus represent the data which least strictly adhere to the original assumptions. For a critical gas saturation of 20% (far higher than would ever be encountered in the field), using the value of OLat the bubble point pressure gives a result which deviates from the solution generated by the very accurate sequential implementation by less than 6%. Using an OLevaluated at the average pressure gives an estimated pressure which deviates by less than 1% from the more accurate result. For fluid properties and pressure regions where (Y does not change rapidly (and for typical values of the critical gas saturation), using the value of (Y at the bubble point gives a very good answer. Example 3 - Estimation of S, at a specified pressure. When the pressure is specified, no iteration is necessary, because (Y is evaluated at the known average between the bubble point pressure

I and 3 (Craft et al.,

% (bbl/SCF)

0.94, and our final estimate critical gas saturation is:

of Petroleum Science and Engineering

y.;25)(l-(&y) + 1](-“(‘-“.y4)

= 1654psia Note that there is little difference between the results from the first iteration and those from the second. Clearly, a third iteration is not necessary.

1.6

. . . . . . 250 deg. F

o.,i~,,,,,,,,,,,,,,,.,,,,,,,,,,,,,,,,,,,,,,,,,,, 0

2000

6000

4000

Pressure

Fig. 4. The product of the phase change compressibility

6000

(psia)

and the pressure for a 0.9 gravity gas

ES. Carlson, P. W. Johnson/

Journal of Petroleum Science and Engineering

and the specified pressure. For the parameters summarized in Table 1, the calculation using Eq. 2 and the estimate of CY from Example 1 indicate a gas saturation of 10.4% when the reservoir pressure is at 1400 psig. Using Eq. 3 gives an estimated recovery factor, relative to the oil in place at the bubble point pressure, of: RF( 14OOpsig) =l-

1 - 0.20

39

approach, is that local behavior can still be analyzed. If one point in the reservoir adheres to the model restrictions, a reasonable estimate of pressure versus gas saturation can be made at that point, so the current method may be particularly useful in numerical simulations.

9. Summary

(1 - 0.104 - 0.20)

15 (1996) 33-43

and conclusions

1.4509 = o.097 1.3978

The rock and water compressibilities have been neglected (assumed to be zero) in this recovery factor estimate, since the reference pressure is the bubble point. Example 4 - Average oil properties. From the figures it is apparent that for a typical oil (35” API) and a typical gas (0.6 specific gravity) 01 has a value of about 0.5. So for an average oil Eq. 2 can be reduced to:

This form is rarely as accurate as is Eq. 1, but it is handy for those in a hurry, and it requires a minimum of background information. For comparison, Eq. 5 yields an answer that is low by about 3% when used on the problem in Example 1, and low by about 6% when used on the extreme case cited in Example 2.

8. Discussion The foregoing examples indicate that with good data our method is easy to use and accurate. In actual reservoirs, variability of rock properties, pressure gradients, and many other factors are going to influence performance. It must be emphasized that the method presented here requires an immobile gas phase. In most reservoirs, a small amount of free gas will flow before the entire reservoir reaches the critical gas saturation. Any reservoir characteristic which hinders the accuracy of this method also will hinder the traditional material balance methods. One major strength of the present method, over the material balance

A procedure is given for estimating the pressure at which a reservoir with an immobile gas phase reaches a specified gas saturation (or for estimating the saturation at a given pressure). Also, with the pressure versus gas saturation relation provided, it is easy to estimate the local recovery factor in portions of the reservoir with immobile gas and constant water saturations. Although we have not presented a contrasting calculation here, we believe that the new method provides a simpler way to estimate reservoir performance by hand than the currently applied material balance techniques. Based on the foregoing discussion and examples, we offer the following conclusions: The method presented herein provides a very quick way to estimate the pressure at which a specified gas saturation will occur in the reservoir, provided that the gas is immobile; Estimates developed using this method agree closely with those from traditional material balance techniques; and The short-cut technique presented in the main text (Examples 1 and 2) gives results that are, at worst, within 6% of those from the very accurate method presented in Appendix A. If the short-cut method is iterated one time the error is reduced to less than l%, but additional iteration provides little further benefit.

10. Nomenclature

Bg Bll Cg cPC

cf

Gas formation volume factor, RB/scf [res. m3/ stock-tank m3] Oil formation volume factor, RB/STB [res. m3/ stock-tank m3] Gas compressibility, psi-’ [Pa- ’ ] Phase change compressibility, psi- ’ [Pa- ’ ] Formation compressibility, psi-’ [Pa-‘]

of Petroleum Science and Engineering 15 (1996) 33-43

E.S. Carlson, P. W. Johnson/Journal

40

L-

k r” P Ph

Psgc

Water compressibility, psiPermeability, mD [m* I Relative permeability with dimensionless Pressure, psia [Pa] Bubble point pressure, psia Pressure at the critical gas

’ [Pa- ’ I

and further expanded

V . pgug + R,V . pouO + p”z~,,. VR, respect

to oil,

ahe)

R, 5 s,c S” SW t

\

A.1. Part I. Derication equation The basic equation (Crichlow, 1977) is:

of an analytical expresgas saturation to pres-

of the ordinary

at

(A41

1

The oil flow Eq. A2 may be applied to simplify Eq. A4 to: V . pgulr + ~0~0 . W

(A9 If it is noted that R, is strictly a function of pressure and it is assumed that Darcy flow applies for the oil, the result is: V

pg~lg

+

p,,uo

.

VR,

=v.pgLlg+p,,($3P)

j$b) aR\

a(P,+s,)

=Appendix A. Development sion relating the critical sure

at

+R ~(Po+%)

psia

Recovery factor, dimensionless Solution gas oil ratio, SCF/STB Gas saturation, dimensionless Critical gas saturation, dimensionless Oil saturation, dimensionless Water saturation, dimensionless Time, s Gas velocity vector, ft/s [m/s] Oil velocity vector, ft/s [m/s] Gas compressibility factor, dimensionless pcpc, dimensionless Porosity, dimensionless Oil viscosity, CP [Pa s] Gas density, lbm/ft’ [kg/m”] Oil density, lbm/ft” [kg/m31

+ P”+%,

at

=-i

[Pa] saturation,

Pal RF

to give:

!

at

+ P&Suat 1

(‘46)

If it is assumed that pressure gradients are small, then the square of the pressure gradients is negligible, so that:

differential

(A? The right-hand

of gas flow in porous media

side of Eq. A7 can be expanded

to:

V . (pgug + R,P,~,) =-

a( R,P&%)

~(P,+%J at

and the basic equation

v. (PO%)

of oil flow is:

a(P&S,, = -

at

!

The gas equation

(‘41)

at

-

1

I

may be expanded

(A21 to give:

V . pgtia + V. RspOcO

rz-

~(p,+%) at

_ ~(RSPOWO) at

(A31

After noting that density and porosity are functions of pressure, this side can be further adjusted to give:

ES. Carlson, P. W. Johnson / Joumal of Petroleum Science and Engineering

Application of the definitions compressibilities gives

for rock and gas

A.2. Part 2. Solution of the ODE If g(p)

-v .PgUg= P,G$

+ P&J

g(p)

(AlO) Definition of compressibility from vapor to liquid as: PO dRS cpc = - -=_pp dp

due to phase change

Rs dRS

(All)

Ro dp

41

15 (1996) 33-43

is defined as:

= ew( /“Cc, + cf - cpC)dp)

(AIT)

then the general solution to this ordinary equation (Boyce and DiPrima, 1977) is:

sg=

differential

%(a - lPdP)c,,dP)

(Al81

where a is an arbitrary constant of integration. At the bubble point, pb, the gas saturation is zero. Using this condition, and defining f(p) so that:

results in: f(p)

-v .PgUg = P&Z

= /“g( p)c,,dp

(‘419)

the constant is found to be: + P,$[ &(cg + cf) + sot,,]

;

a=f(p,)

(A20)

(A12)

so that:

Prior to reaching critical gas saturation, the velocity of the gas is zero. When this is the case, Eq. Al2 reduces to:

S,=(l-S,)

O=~+[S,(c,+c,)+Socp,]~

A.3. Part 3. Approximate

(Al3)

or, noting that So + Sa + S, = 1:

f(

Pb) -f(P)

(A211

g(p)

analytical solution

Calculation of S, depends on the ease with which g(p) and f(p) can be calculated. For all pressures:

O=~+[S,(c,+c,)+(l-Ss-S,)c,,]% (A14) If S, is fixed or otherwise not dependent on time (when water influx is present in the reservoir, in the regions of the reservoir where the water saturations increase there is strong time dependence), and the pressure is not constant Cap/at Z O>, then the time dependence of the above relation can be factored out to give

1 cg=---P

1 dz

and cpc = CY(p>/p. The term CY(p) is pressure (psia), temperature (“F), gas ity, and oil gravity (API) which can with the Standing correlation (Craft et R, as: a(p)

=pcpC =p--

-(l-S,)c,,=~+S,(c,+c,-c,,)

a function of specific gravbe estimated al., 1991) for

Bg(P) = 1.204 -R,(P) Bo( P) (A231

(AW differential

Bg dRs Ro dp

o=~+[Ss(c&,+cr)+(l-S&,-S&,]

This first-order ten as:

(A22)

2 dp

equation can be rewrit-

(Ale)

Applying the definition of g(p), and assuming that cf is much smaller than the other two compressibilities, results in: ln[ g( p)] = jpi

- t 5

a(P) - Pdp

(A=)

42

ES. Carlson, P. W. Johnson/Journal

which can be simplified

of Petroleum Science und Engineering 15 (19%) 33-43

to give:

A.4. Part 4. Stepwise solution ,for greater accuracy

14g( PII = 14 P>- 14 z) - I”-

a(P) dp

P

(A25)

This leads to: ln[ g( p)] = ln( 4)

The approximate solution may be incrementally applied for better accuracy using small changes in saturation. This results in: 1 - o( pi_ ,) + 1

- /‘%dp,or

a( Pi-l) a(P) -1”pdp

g( p) = fexp The term f(p)

.f( P)

(A261

can be evaluated

= / p[ Fexp(

A&E

as:

-i”Fdp)]dp

(A27)

If o(p) and z are approximately constant over the pressure range of interest (where z and a(p) have been evaluated at some point on the interval between p,, and pCsgCl,then: P &T(P) =-

I-m (A28)

z

and: p’-a

01

dp=--

(A29)

1 - CY z,,, This results in: S,(P)

=(’

-SW)

;‘,t”i-_R[jy---

11

(A30) If z(p) S,(P)

is approximately =(I

equal to zaye, then: [( Fj’pu-

-SW)&

1]

=Pb

or, if z(p)

p(Sg)=p,

!

&y,T

1-

o!

u

l

I)

1

(A321

is approximately

&

z,,, -+I d P)

I/(u-

1-o -+1 w 01

where VS,, is the critical gas saturation the number of steps, n. At the final step, ~1, the pressure is: P, = PSgc

divided

by

( A36)

The pressure at the initial step is pb and z(p) is approximately equal to z,,, since the pressure increments are small. This solution has an advantage over direct numerical solution of the ODE, when pressure as a function of saturation is the goal. The functions of the ODE are all defined with respect to pressure so that pressure has to be incremented. A direct numerical solution will involve selection of an incremental pressure, which is applied to calculate a new saturation. The pressure is incremented until the critical gas saturation is reached or exceeded. When the approximate solution presented here is used, the saturation is incremented and the new pressures are calculated. The accuracy of psgC can be established before the interation through selection of the number of increments.

(A31)

From Eq. A3 1, the pressure, PCS,), at which a given gas saturation is reached is approximately:

P&)

(A351

n

equal to z,,,, then:

A.5. Part 5. Solution for the special case, CY= I When (IY(p) equals 1, the foregoing approximate solution is undefined. When CY(p) approaches 1 either from the positive or negative directions, the approximate solution becomes: 5

(A371

I/Cap 1) (A33) I

where t approaches infinity as CK(p> approaches 1. Assuming that 5 approaches infinity as integral num-

ES. Carlson, P. W. Johnson/ Journal of Petroleum Science and Engineering 15 (1996) 33-43

bers, the binomial theorem can be applied to expand the previous expression to give: u2

1 5 l-u-

u3

(‘438)

=l-U+l-_j+.. 5i

i

which of course is the definition of ePU. Therefore, for the case where OL(p) approaches 1:

(

4%) =hexp

-

&-

1

w

This analysis can be easily verified by substituting 1 for a(p) in the original differential equation.

Appendix

B. The material

balance method

For the sake of comparison, we present here the equations required for analysis of reservoir performance by the material balance method. Essentially, current material balance methods involve the simultaneous solution of the equations: dR, dp

=

(Rs - RP) dRF RF

(B-1)

dp

and: dRF -= dp

RF Bt -Bt(PtJ)

dB, -dP

RF Bo+Bg(R,-Rs) _BdR, g dp

+(Rp-R,)% i

oil production. Eq. B2 assumes no net water influx and neglects water and formation compressibility. The two ordinary differential equations are solved simultaneously using numerical techniques (e.g. Euler or Runge-Kutta methods), and are subject to the constraint that S, < S,,. The gas saturation can be given by: 1 - SWi

(A-39)

(B-2)

where B, is the two-phase formation volume factor given by B, = B, + ( Rsi - R,) B, and R, is the ratio of the cumulative gas production and the cumulative

43

Sg= [l-RF(p)] -{I

(l-cf)(Pi-P)

- [1 + (cw +Cf)(pi-P)]‘wi}

B, BOi (B-3)

for any pressure. Although we have never seen the classical subbubble point pressure recovery analysis presented as the foregoing system of ODE’s, the foregoing is the most compact explanation for the iterative technique. Clearly the method presented in the main text is much less complicated.

References Boyce, W.E. and DiPrima, R.C., 1977. Elementary Differential Equations and Boundary Value Problems. Wiley, New York, N.Y., p. 19. Craft, B.C., Hawkins, M. and Terry, R.E., 1991. Applied Petroleum Reservoir Engineering. Prentice Hall, Englewood Cliffs, N.J., pp. 12-45. Crichlow, H.B., 1977. Modern Reservoir Engineering - A Simulation Approach. Prentice Hall, Englewood Cliffs, NJ., pp. 47-52. Dake, L.P., 1978. Fundamentals of Reservoir Engineering. Elsevier, Amsterdam, pp. 19-20. Glaso and Bistein, 1980. Generalized pressure volume temperature correlations. J. Pet. Technol. (June), 32(5): 785-795. Ramey, H.J., Jr., 1964. Rapid methods for estimating reservoir compressibilities. J. Pet. Technol. (April), 16(4): 447-454. Vaaquez, M. and Beggs, H.D., 1980. Correlations for fluid physical property prediction. J. Pet. Technol. (June), 32(6): 96% 970.