THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING r;...., PREDICTION:' OF G.AS STORAGE RESERVOIR.EJ AVIG - ",;...[,, *-. s-*. < *s,F + s. -.r t Keith H. Coats April, 1959 IP-367

Doctoral Committee: Assistant Professor M. Rasin Tek, Chairman Associate Professor John W. Carr, III Professor Donald L. Katz Professor Robert R. White Professor Chia-Shun Yih

ACKNOWLEDGY- ENTS The author wishes to express his appreciation to the following individuals and organizations for their contributions to the research which was the basis of this dissertation: Assistant Professor M. Rasin Tek9 chairman of the doctoral committee, for his sincere interest, useful advice, and most generous donation of time to discussion of the problems which arose in the course of the research. Professor Donald Lo Katz for his invaluable suggestions concerning important and appropriate problems for consideration in the research, for his aid in securing gas field data, and for his advice regarding the application of the research results in practical reservoir problemso The other members of the doctoral committee for their advice and encouragement. Assistant Professor Bernard Ao Galler of the University of Michigan Math Department and Mrs. Shirley Callahan of General Motors for their aid in securing use of the General Motors' IBM 704 computer and in the processing of calculations on that machine. The Michigan Consolidated Gas Company and the Natural Gas Storage Company of Illinois for their donations of gas storage reservoir data. The National Science Foundation and MGA fellowship programs for their financial assistance. Glen C0 Smith for valuable assistance in many phases of the calculations o The University of Michigan Industry Program personnel for their efficient and accurate rendering of the dissertation in its final formo ii

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS..0.................ooooooooooo0oooooooooooo0o ii LIST OF TABLES..a.000000000000...............a................. a 0 vi LIST OF FIGURES oooo.ooooo...........oooooooooooooooo o.......... viii NOMENCLATURE..o............o ooooooooooooooooooooOOOOOOOOOOO.O xi I. INTRODUCTIE ON. O OaOOoooooooooo o 00000000000006060660066 1 I. INTRODUCTION............................ 0 a 6 6 0 0 0 0 d 0 a a 0 c o o a o 1 II. LITERATURE REVIEW..,,...o.... o Oo.....0000.....000 0..... 5 III o FUNDAMENTALS OF GAS STORAGE RESERVOIR PERFORMANCE.....0... 7 A. Mathematical Relationsooooooooooooooooooo......o.... 8 1. Pressure-Explicit Equation for Gas Storage Reservoir Situated on Aquifer.......00.0...006... 19 Bo Application of the Developed Equations to the Prediction of Gas Storage Reservoir Performance...... 23 1o Predicted Volume Variation of a Gas Storage Reservoir for Idealized Pressure Cycles.,........ 24 2. Applications of Equations in an Actual Case Studyoooo00000000 oooooo000000ooo0000 oooooooooooooooooo 33 C, Determination of Effective Reservoir Parameter Values from Gas Field Pressure-Production Data.....,,,.,... 43 1. Case I - Gas Storage Reservoir-Aquifer Systems Which Satisfy the Assumptions in Section III-Ao.. 44 2o Case 2 -'Aquifer StorageT Gas Reservoirs. o..... 52 3. Errors in Predicted Water Influx and Pressure Drop Values Resulting from Errors in Reservoir Parameter ValLues0.000 0000 00...................00. 69 Do Effect of Reservoir Geometry on the Calculation of Water Encroachment into the Gas Field.00.......0...06 74 iii

TABLE OF CONTENTS (CONT'D) Page Io Method of Treating Vertical Pressure Penetration into Aqufer.........4............ o......... 76 a. Diffusivity Equation Including Vertical Pressure Distribution.,....................... 76 b. Equation, Governing Water Influx into Gas.Bubb.le * 79 c. Arialytical Solution of the Developed Equations for the Vertical Pressure Penetration Flow Model.................a..a... e80 2. Numerical Solution of Diffusivity Equation in Elliptic Coordinates............................ 83 a. The Diffusivity Equation in Elliptic Coordinates.................................... 85 b. Equations Describing Liquid Flow Across an Elliptic Boundary................... O.*. 89 c. Nutherical Method of Solving the Diffusivity Equation in Elliptic Coordinates......... 92 d. Comparison Between Elliptic and Radial Flow.. 102 e, Error Analysis 0 0........... 0.................... 112 f, Stability Analysis 00 0 —.................... 118 IV. INDIVIDUAL PROBLEMS ENCOUNTERED IN THE PREDICTION OF GAS STORAGE RESERVOIR BEHAVIOR. 0..O................... 121 A. Treatment of Pressure inter:ference Between Two or More Gas Fields Situated o.n a Common. Aquifer....... 121 1. Existing Theory Relating to Interference Between O0l Reservoirs................. 122 2. Adaptation of Existing Oil Reservoir interference Theory to the Treatment of Gas Field Interference 125 30 Application of t;he Presented Gas Field Interference Theory to an Actual Case Study........... 129 B Prediction of Gas Storage Reservoir Perforjance When Initial Aquifer Pressure is a Non-uniform Function of Radius.......................................! 135 1, Approx:imate Solution to the Diffusivity Equation for the Case of a NonrLuE.form Initial Aquifer Pressure Dtstribueti o 0 0.00............ 136 iv

TABLE OF CONTENTS (CONT'D) Page 2. Application to Aquifer Storage Reeservoir E - Method A.................................. 139 3. Application to Aquifer Storage Reservoir E - Method B..............................*..... 150 C. Treatment of the Moving Boundary Problem Encountered in Aquifer Storage Operations............................. 155 1. Prediction of Aquifer Storage Reservoir Behavior with rb Assumed Constant..................... 160 2. Development of an Analytical Approach to the Moving Boundary Problem and Application to a Case Where rb is Known as a Function of Time.............. 173 3. Application of Moving Boundary Method to the Case Where rb is Unknown as a Function of Time............. 186 V. CONCLUSIONS AND SU MMARY............... 194 APPENDIX I. IBM 650 IT COMPITER PROGRAM.................. 201 A. Program Listing and Information Necessary for Compiling*.***................................... ^.,,. 202 B. Required Input or Data for Program.............. 205 C. Output or Results from Program.................. 206 APPENDIX II. SOLUTION OF DIFFUSIVITY EQUATION GOVERNING PRESSURE DISTRIBUTION IN VERTICAL PRESSURE PENETRATION FLOW MODEL........................4.. 207 APPENDIX III. ROUND-OFF AND TRUNCATION ERRORS ARISING IN THE NUMERICAL SOLUTION OF THE DIFFUSIVITY EQUATION IN ELLITPTIC COORDINATES.................... 212 REFERENCES........4.444,,*4** 4.*444444....... 4,,,,, -^, - -*+..... 217 v

LIST OF TABLES Table Page I Gas Field A Data.................................... 36 II Pressure-Volume Behavior of Idealized Gas Storage Reservoir.................................... 47 III Effect of K1 and K on DEV Value of Idealized Gas Reservoir-Aquifer system.......................... 51 IV Values of PtD and Related Derivatives vs t......... 58 V Summary of Field B Actual and Predicted Performance, November 30, 1957 - July 5, 1958............ 60 VI K3 x K2 and DEV Tabulated as a Function of K2 for Field B......................................... 65 VII tDE as a Function of Dimensionless Time tDE........104 VIII Qt/QtE as a Function of Dimensionless Time tDE...... 109 IX Individual Effects of Au, ZA, and AtDE on the Truncation Error in 1tDE......................115 X Summary of Actual and Predicted Performance of Gas Fields C and D for Period 1944-1956................. 133 XI Summary of Actual and Predicted Performance of Aquifer Storage Reservoir E for Period September 30, 1955 - May 4, 1956.............................e 141 XII P2(l tD) as a Function of Dimensionless Time tD for Aqui1te Storage Reservoir E........ 144 XIII DEV as a Function of Assumed "Gas Loss" and K3 Values for Aquifer Storage Reservoir E........ 149 XIV Summary of Reservoir E Data for Period September 1, 1956 - December 15, 1956............................ 153 XV Summary of Reservor E Observed and Predicted Performance, December 15, 1956 - March 6, 1957,_...... 156 XVI Summary of Reservoir B Observed and Predicted Performance, July 5, 1958 - November 29, 1958........ 162 vi

LIST OF TABLES (CONT D) Table Page XVII Summary of Calculated Performance of Aquifer Storage Reservoir with and without Gas Withdrawal During the Gas Bubble Growth Stage....o o o. 169 XVIII DEV values for Reservoir B, Calculated by Moving Boundary Method.O.O.O... o................ 180 XIX Reservoir B Pressures and Pore Volumes Predicted by Moving Boundary Method, November 30, 1957 - July 5, 1958 e a n o o o o o o e a a a c e.o a e a... O. a o. 183 XX Field B Pressures and Pore Volumes Predicted by Moving Boundary Method of Section IV-C-3.......,....... 188 vii

LIST OF FIGURES Figure Page 1 Typical Gas Storage Reservoir-Aquifer System Configurations.,............. o................. g 2 Representation of Assumed Flow Model............ 10 3 Approximation of Pressure Drop at Gas Bubble Boundary by a Step Function........................... 17 4 Effect of Pressure-Time Cycle on Gas Reservoir Volume........................................... 25 5 Minimum Gas Reservoir Volume vs Time................. 25 6 Effect of Pressure Cycle on Gas Reservoir Volume Variation.oo..ao.o....o^......................... * * 28 7 Effect of Pressure Cycle on Quantity of Gas Stored... 29 8 Effect of Pressure-Time Cycle on Gas Reservoir Volume...................... a......... 30 9 Minimum Gas Reservoir Volume vs Time................. 30 10 Effect of Pressure Cycle on Gas Reservoir Volume Variation.,,,..,,........,.,.,,,......,........ 32 11 Gas Reservoir Pressure-Time Cycles..,........... 34 12 Gas Field Volume Variation with Time................. 34 13 Plan Sketch of Gas Field A Boundary.................. 39 14 Comparison Between Observed and Predicted Gas Storage Reservoir Well-Head Pressures....,....,........... 41 15 Comparison Between Observed Data and Predicted Gas Storage Reservoir Volume Varittion............... 42 16 Effect of Assumed K1 and K2 Values on Agreement Between Predicted and Observed Gas Bubble Volume....... 49 17'Observed' and Predicted Gas Bubble Pore Volume for Idealized Reservoir-Aquifer System................. 50 18 DEV vs K3 for Aquifer Storage Reservoir B.......... 61

LIST OF FIGURES (CONTID) Figure Page 19 Field B Observed and Predicted Pressures vs Time, November 30, 1957 - July 5, 1958..................... 63 20 Field B Observed and Predicted PreSsures vs Time, November 30, 1957 - July 5, 1958.................... 64 21 Field B Actual and Predicted Pore Volumes vs Time, November 30, 1957 - July 5, 1958..................... 67 22 Field B Actual and Predicted Pore Volumes vs Time, November 30, 1957 - July 5, 1958............... 68 23 Vertical Pressure Penetration Flow Model............. 78 24 4eometric Relationships Between Cartesian and Elliptic Planar Coordinates,.............. 88 25 QtD and QtDE vs Dimensionless Time........ 105 26 Comparison Between Elliptic and Radial Flow.......... 110 27 Field C and D Predicted and Actual Pore Volume Ratios vs Time.....,....................... 130 28 Field C and D Predicted and Actual Pressure Ratios vs Time...................*. *. e 0. 0 0... d.. 0. 131 29 Pressure Distribution in Aquifer Underlying Reservoir E, September 30, 1955,, e,...,.................... 142 30 Effect of Assumed and "Gas Loss" Values on DEV Value for Aquifer {orage Reservoir E............ 146 31 Comparison Between Aquifer Storage Reservoir E Observed and Predicted Pressures for Period September 30, 1955 - May 4, 1956......................... 147 32 A. vs Dimensionless Time for Aquifer Storage Reservoir 33 Aj vs Time for Aquifer Storage Reservoir E, September 1, 1956 - December 15, 1956.................... 154 34 DEV vs K3 for Aquifer Storage Reservoir E, December 15, 1956 - March 6, 1957.,.o........................ 157

LIST OF FIGURES (CONT'D) Figure Page 35 Comparison Between Aquifer Storage Reservoir E Observed and Predicted Pressures, December 15, 1956 - March 6, 1957..6#0...o e......... 158 36 Comparison Between Aquifer Storage Reservoir B Predicted and Actual Pressures, July 5, 1958 - November 29, 1958.............................. 164 37 Comparison Between Aquifer Storage Reservoir B Predicted and Actual Pore Volumes, July 5, 1958 - November 29, 1958................... 165 38 Effect of Gas Withdrawal During Growth on Pore Volume and Gas Content of Aquifer Storage Reservoir. 172 39 Feilationship Between rb and y in Aquifer Storage Reservoir B......o0.......................... 177 40 DEV vs K3 for Aquifer Storage Field B............. 181 41 Comparison Between Field B Actual Pressures and Pressures Predicted by Moving Boundary Method, November 30, 1957 - July 5, 1958..........eh...... 184 42 Comparison Between Field B Actual Pore Volumes and Pore Volumes Predicted by Moving Boundary Method, November 30, 1957 - July 5, 1958..*................. 185 43 Comparison Between Field B Observed Pressures and Pressures Predicted by Moving Boundary Method of Section IV-C-3....,..........,..... 190 44 Comparison Between Field B Actual Pore Volumes and Pore Volumes Predicted by Moving Boundary Method of Section IV-C-:3....,..... o*. 191 x

NOMENCLATURE AJ pressure drop, at time t = jAt, due to initial non-uniform pressure distribution, psia c sum of aquifer water and formation compressibilities, vol/ vol-psia d distance between centers of oil or gas fields, feet DEV average fractional deviation between predicted and observed pressures f foci of ellipse in x-y plane are at x = + f, y=O fmn (Au)2 [cosh2 (mAu) - cos2 (nAv)] AtDE h thickness of aquifer formation, or of flow model, feet JO, J1 Bessel functions of the first kind and of Oth and 1st order, respectively k permeability, millidarcies or ft-lbm/sec2-psi k/t mobility, ft2/sec-psi kR ratio of permeability in vertical direction to permeability in horizontal direction, dimensionless K1 jhrpcrb2, ft3/psi K2 k/icpcrb2, 1/seconds Kt - (k/icpc) (r cp/6o0RT)2/3 K3 Ct/2ihkAt, psi/ft3 K constant determined from gas compressibility data, z = l+KP, l/psia nt gas in reservoir at time t, lb. moles ni or nj nt at t = iit or jAt, respectively N total number of time increments included in period for which performance of reservoir is to be calculated P pressure, psia xi

NOMENCLATURE (CONT'D) Pt pressure at time t, psia Pj or Pi Pt at t = Ajt or iAt, respectively P pressure drop, Po - P, psia PO initial pressure, psia PtD dimensionless Qre sure drop quantity, tabulated vs tD in the literature 1) 2), at dimensionless time tD PaiAt observed gas bubble pressure at time t = iAt, psia &Qt cumulative influx of water into circular gas field at time t, ft3 QtE cumulative influx of water into elliptic gas field at time t, ft3 QtD dimensionless production guantity for radial flow case, tabulated in the literature l)(2) vs tD, at dimensionless time tD QtDE dimensionless production quantity for elliptic flow case q rate of water movement, ft3/second 9q average rate of water influx into gas bubble during time increment from (i-l)At to iAt, ft3/second Aqi qi+l qi r radius, feet rb gas reservoir radius, feet re aquifer exterior radius, feet rD dimensionless radius = r/rb rDe dimensionless exterior radius, re/rb R gas constant, 10.73 psia-ft3/lbomole-~R, or (in elliptic flow section) (u)2 AV Sj slope of predicted pressure-time curve, psi/sec. T temperature, ~R

NOMENCMLATURE (CONTI' D) t time, seconds tD dimensionless time = (k/ipcrb2)t = K2t tDE dimensionless time for elliptic flow, (k/1cpcf2)t TD dimensionless ttime for moving boundary method, TD (k/pc)fJ dt o rb2(t) At time increment, seconds AtD dimensionless time increment u,v,w elliptic cylinder coordinates (u,v dimensionless, w in feet) ub value of u on interior ellipse (gas field), dimensionless ue value of u on exterior ellipse (aquifer exterior boundary), dimensionless Uo(XnrD) Yo(XnrDe) Jo(XnrD) - Jo(XnrDe) Yo(XnrD) V velocity vector, ft/sec. Vx, VyP Vz velocity vector components in x, y, and z directions, ft/sec. Vt gas field pore volume at time t, ft3 Vi or Vj Vt at t = iAt, or jAt, respectively, ft3 VaiAt actual gas field pore volume at time t = iAt, ft3 V0 initial pore volume, ft3 AVi 2Vi - Vi_1 - Vi+l x, y, z cartesian coordinates, ft Y0o Y1 Bessel functions of second kind, and of Oth and 1st order, respectively z cartesian coordinate, ft, or gas compressibility factor, dimensionless p density, lbs mass/ft3 c@ ~ porosity, fraction xiii

I. INTRODUCTION Natural gas producing reservoirs often occur in close proximity to large water bearing formations called aquifers. When a gas reservoir is adjacent to an aquifer, the gas-water interface becomes mobile under the effect of pressure gradients resulting from gas production. The rate of water movement is a function of the aquifer formation and fluid properties and the pressure-time relationship at this interface. The purpose of the research herein described is the development and application of general methods of predicting the performance (e.g., volume-time and pressure-time relationships) of gas storage reservoirs situated adjacent to aquifers. The storage of natural gas in underground reservoirs is referred to as either "depleted-field storage" or "aquifer storage". The former term designates gas storage operations in oil or gas fields which are converted to gas storage use after a period of natural depletion. The latter term denotes storage operations in a gas reservoir which was initiated and grown by injection of gas into an aquifer formation originally containing no gas. Gas is usually injected into a storage reservior during the summer months when consumer demand is low and is withdrawn during the winter months when the demand is high. This cyclic gas injection and withdrawal results in cyclic periods of high and low gas reservoir pressures which in turn result in aquifer water movement out of and into the gas field. This water movement results in a variable gas reservoir volume and it becomes very desirable to be able to predict the gas field volume-time and pressure-time behavior for a postulated gas injectionproduction schedule. -1

-2The prediction of the effect of aquifer water movement on the volume and pressure variations of an adjacent gas storage reservoir can provide information valuable in the study of several reservoir engineering problems. The prediction of the pressure behavior is important in relation to the deliverability potential of the gas storage reservoir. For a fixed number of wells, the deliverability (or gas production rate) of a storage field is a function of the field pressure. A prediction of the pressure during the high demand winter months thus allows calculation of the deliverability potential for this period and comparison of this potential with the estimated demand. The economics of gas storage operations are directly influenced by the influx of aquifer water since a shrinking storage reservoir results in increasing pressures for the same quantity of gas in storage. Material balance and reserve or recovery calculations obviously require knowledge of the reservoir volume-time relationship. Various other reservoir engineering calculations such as interpretation of well interference data, evaluation of physical characteristics of porous media, water coning and pressure maintenance studies are typical problems where variations of reservoir volume due to edge or bottom-water encroachment becomes important. The available literature relating to the performance of gas storage reservoirs situated adjacent to aquifers consists primarily of the papers by Van Everdingen and Hurst(l) and Chatas(2) which present dimensionless, tabular solutions to the partial differential equation governing unsteady state liquid flow through porous media. A paper by Katz, Tek, and Coats(3) describes the results of the initial stages of

-3the research reported in this thesis. Certain design aspects of the development of gas storage fields are treated by Katz, Vary, and Elenbaas.(4) The general purpose of this research has been the development of fundamental theory relating to the performance of gas storage fieldaquifer systems and evaluation of this theory by application to actual gas field studies. Several sub-problems were studied in order to attain this general purpose. Thus the specific objectives of this research were the development and application of the following theory and methods of calculation to actual gas storage reservoirs: 1) A general digital computer method to predict the pressure and volume behavior of a gas storage reservoir from any postulated gas injection-production schedule. This method should employ the assumptions of a constant gas bubble radius, a constant initial pressure throughout the aquifer, circular geometry of the gas field, negligible vertical flow effects, and absence of interfering, neighboring gas fields. 2) A method of determination of the "effective" values of reservoir physical constants, i.e., those values which result in best agreement between the calculated and actual field performance during the initial period for which data is available from the field. 3) A gas field-aquifer system model which will account for the effect of vertical pressure penetration into the aquifer.

-44) A numerical method of solving the unsteady state diffusion equation in elliptic coordinates. Use of this method should make possible certain conclusions concerning the effect of the gas storage reservoir areal geometry on the reservoir performance. 5) A digital computer method of treating pressure interference between two or more gas storage reservoirs situated on a common aquifer. 6) A predictive method which will take into account any nonuniform initial pressure distribution throughout the aquifer, 7) A method of predicting the performance of an "aquifer storage"t reservoir which accounts for the fact that the gas bubble radius is not constant but increases continuously from 0 feet as the gas bubble is grown.

II. LITERATURE REVIEW With the exception of Reference 3, the literature contains no material relating directly to methods of predicting the performance of gas storage reservoirs situated on aquifers. Schilthuis in a 1935 article(5) described calculations, relating to unsteady state aquifer motion, which reproduced quite well the observed pressure performance of the East Texas oil field. As mentioned above, Van Everdingen and Hurst(l) and Chatas(2) presented tables of dimensionless quantities obtained by solution of the partial differential equation governing isothermal, radial, unsteady state liquid flow through porous media (hereafter designated as the diffusivity equation). Their solution was obtained for the two cases of "constant terminal rate" (the rate of water movement across the gas-water interface is assumed constant) and "constant terminal pressure" (the pressure at the gas-water interface is assumed constant). They also indicated a method of combining the superposition principle with their solution to treat a variable rate or variable pressure case. Neither Reference 1 nor Reference 2 includes any application of the developed theory to an actual field study. Mortada(6) presented a method of treating interference between two or more oil fields situated on a common aquifer which requires use of the solutions presented by Van Everdingen and Hurst or Chatas. He presented illustrative calculations but did not apply the method to an actual field study. Warren(7) presented a numerical solution to the partial differential equation governing linear, unsteady state gas flow through a porous medium. He considered his results applicable to the -5

-6case of a linear gas storage reservoir bounded on all sides by an impermeable rock. Books by Muskat(8), Carslaw and Jaeger(9), Churchill(lO), and Sneddon(ll) contain analytical solutions to the diffusivity equation for a variety of boundary conditions which cover a wide spectrum of practicality in reservoir engineering problems. Muskat and Woods(l2) showed that, in the case of oil reservoirs, a minimum deviation between the actual or known oil in place and the values calculated from the material balance equation (including an aquifer water influx term) is a poor basis for choice of reservoir physical constants. Their conclusion was that if core, logging, or other data were available to fix the values of all but one of the reservoir parameters, then agreement between observed and calculated field behavior serves as a good basis for selection of this one remaining unknown parameter.

III. FUNDAMENTALS CF GAS STORAGE RESERVOIR PERFORMANCE The performance of gas storage fields can be investigated by either a mathematical or an empirical approach. The latter approach, which would consist of determining empirical relationships between the performance and certain operational variables of a specific storage reservoir, might prove useful in predicting the performance of the reservoir concerned but would very likely fail to be useful in other cases. The mathematical approach, however, has the potential of providing a general method capable of predicting the performance of many reservoirs. From an academic as well as a practical point of view, the achievement of a general predictive method is certainly preferable to the development of a myriad of empirical correlations, each of which is developed for and applicable to only one field. However, mathematical analysis of storage reservoir behavior requires certain assumptions concerning the reservoir geometry and formation and fluid properties. The generality of the methods obtained by the mathematical approach will therefore be restricted insofar as the storage fields fail to conform to these assumptions. The present section (III) of this thesis treats the storage reservoirs which satisfy, in general, the flow model chosen and assumptions given in Section III-A below. Realizing that the methods developed from these assumptions will not be applicable to all cases, the author has treated in Section IV certain important storage cases which definitely fail to satisfy certain of the assumptions in Section III-A. -7

-8A. Mathematical Relations Calculation of aquifer water movement and its effect on the performance of a gas storage reservoir requires adoption of a flow model, certain assumptions concerning aquifer formation and fluid properties, and formulation and solution of the appropriate differential equation. Three typical gas field-aquifer system configurations are shown in Figure 1 where rb denotes the gas bubble radius, h the aquifer thickness, and re the aquifer exterior radius. The aquifer exterior boundary is assumed at infinity if the boundary is so far removed from the gas bubble that its effect on the gas bubble performance will be negligible for any practical length of time. The flow model adopted for the present study approximates the systems shown in Figures la and lb and is shown in Figure 2; a model which will approximate the system of Figure lc will be discussed in Section III-D-1 of this thesis. Figure 2a shows the assumed position of the gas bubble relative to the aquifer and Figure 2b indicates the radial nature of the aquifer water flow. The following assumptions are made concerning the flow model and aquifer formation and fluid properties: 1) The rate of water influx into the gas bubble is equal to the rate of water movement in the negative r direction across the vertical dotted line boundary at r = rb indicated in Figure 2a, 2) Uniform aquifer formation permeability and porosity. 3) Negligible capillary, gravity, and vertical flow effects. 4) Single phase, isothermal, radial, unsteady state aquifer water flow.

-9GROUND LEVEL WELLS IMPERMEABLE SHALE LAYERS l GAS BUBBLE AQUIFER r_._re or O) la. Edge-Water Drive System. UBBLE = E L L~~ELS h\ AQUIFER r %rre,or O I IMPERMEABLE SHALE { b lb. Bottom-Water Drive System with Definite Aquifer Thickness. GAS BUBBLE iMPERMEABLE SHALE AQUIFER - -- r-re,or CO - - __X rb_ -- h large DISCONTINUOUS SHALE STREAKS la. Bottom-Water Drive System with Large Aquifer Thickness. Figure l. Typical Gas Storage Reservoir-Aquifer System Configurations. Figure 1

-10GAS BUBBLE AQUIFER h I, AQUIFER i, /I r -— re or C IMPERMEABLE SHALE LAYERS 2a r rb =re or a) 2b Figure 2. Representation of Assumed Flow Model.

-115) rb, the gas bubble radius, is constant with respect to time. 6) Absence of pressure interference effects from neighboring gas (or oil) fields. The above assumptions and described flow model now allow formulation of the partial differential equation governing the aquifer water flow. This equation is derived by combination of the continuity Equation (III-1), which is obtained by a mass balance about an infinitesimal volume element, with Darcy's flow Equation (III-2) and an equation of state (III-3) for the aquifer liquid.(1,8) oaX(pv,) t y(pVy) + ) ~ aC (P 1) d;/(f -o apt (III-1) "P = — No ofe / * c (P-Pa)] ~or c 5X//-3) c = sum of aquifer water and formation compressibilities, vol/vol-psi e = base of natural logarithm k/. = aquifer formation mobility, ft2/sec-psi P = aquifer water pressure, psi Pb = arbitrary pressure base, psi t -= time, seconds V = velocity vector, ft/sec

-12x,y,z = cartesian coordinates, ft p = aquifer water density, lbs/ft3 po = aquifer water density at P = Pb, lbs/ft3 cp = aquifer formation porosity, dimensionless Substitution of V and p from (III-2) and (III-3) into (III-1) yields Equation (III-4) below. 2P 1P +lr -p (III-4) vP'E t At (iii-4) Equation (III-4) can now be expressed in cylindrical coordinates(8) with the terms d2P/6g2 (@ = polar angle) and 62p/az2 set equal to 0 since the flow is assumed radial with negligible vertical flow effects. Thus (III-4) can be rewritten as in (III-5) a P # I ap _ 9P (III-5) where P = Po - P, pressure drop (Po is the initial aquifer pressure, psia), rD = r/rb, dimensionless radius, and tD = kt/Cpcrb2, dimensionless time. Equation (III-5) has been solved by Van Everdingen and Hurst(l) for the two cases of'constant terminal pressure' and'constant terminal rate'; the details of the solution are therefore omitted here. The initial and boundary conditions employed in the constant pressure and constant rate cases are given below in Equations (III-6) and (III-7), respectively.

-15(II-6) hrmit PCP t ~ P5), 8 = /(0 Thus the initial condition in both cases specifies a uniform pressure, Po, throughout the aquifer at time 0. The constant pressure case specifies that the pressure is decreased by unity at the inner aquifer boundary at time 0 and held constant at that value for all time. The constant rate case assumes a constant rate of aquifer water movement across the inner aquifer boundary, r = rb. The condition limit P(rDtD) = 0 is r 0)- equivalent to the assumption of an infinite aquifer. Van Everdingen and Hurst also presented some solutions for the case of limited aquifers. Their solution for the constant pressure case consists of dimensionless production quantities, QtD3 given in tables as a function of dimensionless time tD. The QtD values are employed in the calculation of Qt as shown in Equation (III-8), Q = 2rhc A a P (III-8) where Qt = cumulative influx of aquifer water into the gas bubble at time t, ft3

LP=Po-(P)rli; (P)r D= = constant pressure maintained at edge of gas bubble (r=rb), psia QtD = dimensionless production quantity at dimensionless time tD = kt PPcrb The constant rate solution consists of dimensionless pressure drop quantities related to the pressure drop in the gas bubble as shown in Equation (I-9).9) Pk t - hk 5ZD II-9 Pt = pressure at gas bubble boundary r = rb, at time t, psia Po = initial aquifer pressure, psia q specified constant rate of water influx into the gas bubble, ft3/sec ktD = dimensionless pressure drop quantity, tabulated as a function of dimensionless time by Van Everdingen and Hurst(1l), at dimensionless time tD kt Lcpcrb2 The restrictions of constant gas bubble pressure or constant water influx rate severely limit the usefulness of Equations (III-8) and (III-9). Duhamel's super-position principle(lo0) however, can be employed in conjunction with these two equations to yield solutions valid for a varying gas bubble pressure or water influx rate.(l) This principle states that if ~l(rDPtD) is the solution to Equation (III-5) for a condition P(l,tD) = 1, if Yl(rDs) is the La Place transform of Pl(rDtD), if P2(rDtD) is the solution to Equation (III-5) for a condition P(1 tD)=f(tD),

-15and if Y2 rD s) is the La Place transform of P2(rDtD)9 then Y2( ~P5) = 5'S7)y(tr:5) (III-10) where f(s) is the La Place transform of f(tD). QtD' the dimensionless function tabulated by Van Everdingen and Hurst, is defined as(l) Employing the property of the La Place transformation [where [f(t)] = e-stf(t)dt = the La Place transform of f(t)] 0o one finds that Q(s). the La Place transform of QtD' is =;)(5) _-' _ } Now, (QtD)2) the dimensionless production quantity corresponding to the'time varying reservoir pressure' solution, P2(rDtD)I is (Q) SOt[ - and the transform of (QtD)2 is ( Q(s)x - - 5 I; &~ Y2 (TP, )g, Since Y2(rDs) = sf(s) Yl(rDs) (Equation III-10), 2 -f~S~)C~Y,(r ~j4 ((;(s - - f (S g Ye (b, O, I){US; cfcS ~ S 9 ) eg -

-16and (QtD) is found, by application of the convolution theorem in La Place transform theory, to be (Qt3) ) -If(r;7 (,-r) C, r (III-11) Now, Qt, the cumulative influx into the gas bubble at time t, is related to QtD as shown in Equation (III-12). 9tz T Ct t (III-12) Qt for the varying reservoir pressure case is then obtained by insertion of (QtD)2, the dimensionless solution for the varying pressure case, into Equation (III-11) to yield = 2 m'c'4 (fr') ar (t- r) d f (111i-13) t3 t~. If, now, the pressure drop at the gas bubble boundary r = rb, P(l,tD), is some function of time and is approximated by a step function as shown in Figure 3, then f(T) [f(T) = P(1,T)] in each of the integrals of Equation (III-13) can be taken as a constant equal to P(l,tDi) + P(l,tDi_l) 2

-17Ca 0 a. 3'n1 4. JQJ _ /,'D,'!,, _ 0 tD tD2 tD3 tD4 tD, DIMENSIONLESS TIME Figure 3. Approximation of Pressure Drop at Gas Bubble Boundary by a Step Function.

-18Equation (III-13) then becomes (III-14) +; (;tp,' -t-t">,) 2 (0 D where Pi P(1, tDi) = P - P(l,tDi) P - Pi tD =.iAtD tD = jAtD t = jAt Rearrangement of Equation (III-14) yields (III-15), a result which appears in Reference 1. 4 2 _ _ L, -t~f 2 (III-15) - ~2rrh~~~l,4 c~,~~ r++*i,,-,... 2 2 where p i P i-i'- Pi+l where AP Pi-l P.-; P-1 = Po psi Pi = pressure at gas bubble boundary r = rb at time t = iAtt, psia QtD-tDi = QJtD-i&tD = Q(j-i)AtD = dimensionless production quantity evaluated at dimensionless time tD = (j-i)AtD AtD = kAt:pcrb2 At = chosen time increment, seconds t = jAt, seconds

-19A similar application of the superposition principle to the constant rate solution, Equation (III-9), yields Equation (III-16) as the varying rate case solution. (l) A gr = s SE " B: -P~jA i t (111-16) where Aqi q+l - qi'= average rate of water influx into gas bubble during time increment from (i - 1) At to iAt, ft3/sec. 1. Pressure-Explicit Equation for Gas Storage Reservoir Situated on an Aquifer The gas law can be combined with a simple material balance and Equation (III-15) or (III-16) to give an equation explicit in the gas bubble pressure. The material balance Equation (III-17) simply expresses the fact that the gas storage reservoir pore volume is equal to t1he original volume minus the cumulative influx of aquifer water. ~' -c -t (III-17) Vt = gas bubble pore volume at time t, ft3 Vo = original gas bubble pore volume, ft3 The gas compressibility factor z can be substituted from Equation (III-18) into the gas law (II-19) to give Egqation (r1X.20) by rearrangement.

-20= - / + AkdP (III-18) f = nt R T (III-19) V, = ktRR((/< (III-20) K = a constant determined from gas compressibility data, l/psia z = gas compressibility factor R = universal gas constant, 10.73 psia-ft3/lb.mole - OR T = gas reservoir formation temperature, OR Pt = static gas storage reservoir pressure, psia nt = lb. moles gas in place in storage reservoir at time t The representation of z as a linear function of pressure is in general valid for pressures less than 700 psia. For pressures greater than 700 psia, additional terms should be present on the right side of Equation (III-18). In most practical cases, however, z can be approximated as a linear function of pressure over the limited pressure range encountered in the gas storage operations, regardless of whether the average pressure is above or below 700 psia. Equations (III-15), (III-17), and (III-20) are now three equations in the three unknowns Pt, Vt, and Qt. The role played by Pt in Equation (III-15) can be seen more clearly by removing the last term from the summation as shown in (III-15a) below.

-21-:2j12 Qt= u 7wh S C is g 29~j~~- z Pi @J - ~' 5Dt ({III- 15a > where Pj = Pt and t = jAt Eliminating Vt and Qt from Equations (III-15a), (III-17), and (III-20), one obtains (III-21). ___Ky _ _S__+ (III-21) 1t<5 k, -.;- —, I~___ __ —:2Equation (III-21) is an explicit relation for Pt, the gas bubble pressure at time t = jAt, provided the calculations for Pt are begun at t = At and continued in order t = 2At, 3At,..., ( - l)At. For if Pt has been calculated at all time steps prior to t = jAt, then all the pressures appearing on the right side of (III-21) are known at time step j (at t = jAt). The water influx rate terms, q4, in Equation (III-16) can be combined with (III-20) to yield a pressure explicit equation. Since q' is (as defined above) the average rate of water influx into the gas bubble during the time interval from (i-l)At to iAt, q! can be defined as e- VL- VL *; - = at where Vi is the gas bubble pore volume at time t = int. Substituting this expression for qj into (III-16) and removing the

-22last term from the summation, one obtains (III-16a), below. &, At 27h A a (jh i) 4 t (III-16a) A (2 V _v_ -)-fPtJ Eliminating Vt (= Vj) from Equations (III-20) and (III-16a), one obtains (III-22). Pt (-Igo + #k7 (III-22) where = ([ - -,i-Vi,). C a.t) s at s' *(2 V.., -., -RTC)J - o 7= 13P T, t.V n,,A 71<3 /,/2T7Akat PO = original aquifer pressure, psia nj = lb. moles gas in storage reservoir at time t = jAt t = jAt All the pressures on the right side of (III-22) are known at time step j if calculations are started at t = At and continued at t = 2At, 35At,..o, (j-l)At. Equations (III-21) and (III-22) are useful only when the gas in place quantities, ni, are known or specified and the gas bubble pressure-time relationship is desired. When the gas storage reservoir pressure is specified as a function of time, then the two Equations (III-21) and (III-22) must be rearranged to allow calculation of the unknown gas in place quantities.

-23(III-23) and (III-24), below, are such arranged versions of (III-21) and (III-22), respectively. L =j-l a k, = Pa Q ( III-23) My_ -P.~ ~5icFr,-,] 8 Dt (O-( (Z)At (III-24) k3 +9 L -( B. Application of the Developed Equations to the Prediction of Gas Storage Reservoir Performance Application of the equations developed in Section III-A is now made in order to (1) determine the effect of the shape of the gas field pressure-time cycle on the gas reservoir volume behavior, and (2) evaluate the applicability of the equations to actual gas storage reservoir studies. The first application should yield answers to some of the following important questions concerning the pressure cycle-volume variation relationship: For the same number of gas field pound-days above and below discovery or initial pressure and at the end of one or several complete cycles, does the gas reservoir volume return to its original or discovery value? Is there a definite phase difference between the time of change in the reservoir pressure and the time corresponding to the response in the aquifer movement? How does shape of the gas field pressure-time curve affect the water movement? Does the time at which a gas reservoir is converted from production to storage operations appreciably affect the size or

-24subsequent size-variation of the gas reservoir? Answers to these and other important questions can best be obtained by simulating some idealized cases by digital computations based on the equations of Section III-A. The second application of the developed equations involves comparison of the predicted performance of an actual gas field with available observed data. An IBM 650 digital computer program has been written to facilitate the calculations involving use of Equation (III-21). This program is described and listed in Appendix I. 1. Predicted Volume Variation of a Gas Storage Reservoir for Idealized Pressure Cycles Equations (III-15) and (III-17) have been employed in calculating on an IBM 650 computer the gas storage reservoir volume behavior effected by nine different pressure cycles. The pressure cycles are specified in a manner which will facilitate answers to certain of the above questions concerning pressure cycle-volume behavior relationships. The pressure cycles are intended to represent the static gas bubble pressure, i.e., the gas (or water) pressure at the boundary r = rb. Calculations for pressure cycles 1, 2, and 3 (see Figure 4) allow certain conclusions concerning the effect of amplitude and duration of the pressure drawdown portion of the pressure cycle on the reservoir volume behavior. Each of the three pressure cycles has an equal number of pound-days above and below the reference P/Po = 1.0 and has a period of one year.

-25o 1.4 1.2 1.0 I L0.8 --—. —-----—. _P-Ol 1000 PSIA L 0.6 V0= 2 x 107 ft3 0.4 c 3.3 x 10-6 1/ PSI 4= 0.1 1.05 k//" = 320 Millidarcys /cp h = 20 ft o 1.00 rb=1780 ft 0.95 o 0.90 0.80 0.80 — Cycle2 0.75 | I L I I --- ""Cycle 3 1 0.75 0 1 2 3 4 5 6 7 8 9 10 11 12 13 TIME IN MONTHS Fig. 4 Effect of pressure-time cycle on gas reservoir volume o 0.90 0.88 ~.o- O-O-OOO- - pp o 0.88[ — o -o-43-o -O-O-O- -O-O 0.. o-0o.0.0.m0__0__0- 0-0-0-0-0-0-0-0-0-0 0.86, 0.84 a Calculated points, cycle 1 = ~ 0.82P~ L 0 Calculated points, cycle 2 0 I I i I Co Calculated points, cycle 3 0 2 4 6 8 10 12 14 16 18 20 TIME IN YEARS Minimum gas reservoir volume vs time Figure 5

Pressure cycle 1 represents a large amplitude, short duration pressure drawdown followed by a mild, extended gas injection period. Cycle 3, on the other hand, represents a low-amplitude extended duration pressure drawdown followed by a relatively large-amplitude short duration gas injection period. Cycle 2 represents the intermediate case. The formation and fluid properties used in calculations for pressure cycles 1, 2, and 3 are noted on Figure 4. The aquifer fluid is assumed to be water of one centipoise viscosity and compressibility noted on Figure 4. Calculated volume ratio curves (Figure 4) clearly indicate that pressure cycle 1 results in the maximum gas reservoir shrinkage. Another result of interest is that the gas reservoir volume does not assume its original value after one complete pressure cycle but exceeds that value by 2 - 4%, although the net pound-days after one cycle is the same above and below discovery pressure ratio for all three cases. Perhaps the most interesting observation is that although each of the three pressure cycles has as many pound-days above the original pressure as below, the corresponding calculated V/VO curves show that the net gas reservoir volume is less than its discovery value for 10-ll months of the 12-month cycle. This fact is of practical interest because of possibility of water coning and direct proportionality of amount of gas stored to gas reservoir volume. If the pressure cycle curves in Figure 4 are reversed about the axis P/Po = 1.0 to

-27represent a gas injection period followed by a gas withdrawal period, the corresponding V/Vo curves will simply be symmetrical about the axis V/Vo = 1.0. Equation (III-15) shows this since reversal of the pressure cycle curves will simply change the algebraic sign of each APi in Equation (III-15). Thus, if the gas reservoir is converted to storage operations upon discovery by injecting gas at the beginning of the cycle and withdrawing gas at the end of the cycle, net gas reservoir volume will be greater than the original value for 10-11 months of the 12month cycle. Figure 7 compares the number of moles of gas stored in the reservoir for the two pressure cycles shown. Yearly minimum gas reservoir volumes corresponding to pressure cycles 1, 2, and 3, are plotted versus time in Figure 5. The pressure cycles were repeated for 20 years to allow calculation of points plotted in the latter figure. The points or curves clearly indicate growth of the gas reservoir with time. This increase in the average reservoir size is also shown on Figure 6 where curves of the relative volume ratio V/Vo of the gas field were plotted versus time corresponding to pressure cycles 1, 2, and 3, repeated for 20 years. This gas reservoir growth will later be seen to be dependent upon the fact that the pressure cycle curve represents an initial period of gas withdrawal followed by gas injection. Pressure cycles 4, 5, and 6, shown in Figure 8, again have the same number of pound-days above and below the reference

1.080 _ cycle 1 cycle 2 1.060 ~ cycle 3E 1.040 1.020 > 1.000-' I 0.980 o 1 i' 0~ Iti I 1' 0. 960LI I I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0.9:E 0940' 6' \!II L 7 4 41 X~~ 0.920 0.900 \ I p 0. 880 j, /~ I,5,I~~~~~~~~~~~~~~~~~~~~~~~~i Il IiIIaI!t 0.860 0.840 0.820 0 2 3 4 5 16 17 18 19 20 TIME IN YEARS Effect of pressure cycle on gas reservoir volume variation Figure 6

I~ I6 I I I 1 I I I 1 I II 62IIr I I1 I I I I Iv'o 6.6, b {~b O b O bO b bhdb O c,.6~ I I P ~II 1I II I. I~~ ~ ~ ~ I *I I I I I I I II "6. I 1 I I _1 II I' r II I~ ~ ~ I II I C3.8- ~ I I I I II I I 1' I 1 I I I I I I I I [, 1 ( I I I;g I I i I I I I i I I I o)5.6~'" " I I I I Ij I I I I, I. I I5.2 I! I I I I I I I I I II tI//! ~I I I I I Ili I! I r\)! j I I I I I4 8 I i I I I I *: 5.8 I I I I I I I I I I I I 14. I' 0210121 1 I I I I I i [ I I I I I I I 024681012TIMEINMTIIS I I''II! I I I I I [ I i- - I I I I I I,,, -I I i I 32 I I I I /11HIi d I I I 1 r I 1~, I i, [ I Cas i [1M [i [OI ] [ 111' /'!!111 I~~ I.0 I 2 3 4 5 6 7 10 I I 12 13P 14 15 16 1 8 92 TIM IN YEARS''''' ~ I I I I~~Efec of presr cyl on quntt of ga s mtored O I I I I I I I I I 7 I III I r I I I I I I ~~~~~I I [ - - - I I 36 I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~!1 ~ ~ ~ II I! I! I] ~ IIj I!i~ il 1! Iti II- -— "~'d ~brt~rj TIMI~~~~~~I II i!~il F'i~~~w(as 3 o,, i

-301.5 0 0 0.51 _ J VO= 2 x 1O- ft3 1.05 050.5=,~IP 1000 PSIA L 0.5 _ Vo 2 x 107 ft5 c = 3.3 x 10-e I /PSI 0.81, 0 >.1 k/L = 320 Millidarcys /cp 0.90 0Cc p C ycle 6 o.75 0 2 4 5 6 7 8 9 10 II 12 13 TIME IN MONTHS Minig. 8 Effect of pressure-time cycle on gas reservoir volume F0.85 0.83 51"00 o08oI o.81o'. 0. 0.0-0 0- 0r 0m4:aC 0.79 Calculated points, cycle 4 ~(~,e~'0"0' 03 Calculated points, cycle 5 Una o7~ 0 Calculated points, cycle 6 0 2 4 6 8 I0 12 14 16 18 20 TIME IN YEARS Minimum gas reservoir volume vs time Figure 9

-31P/Po = 1.0. Upper and lower pressure limits are the same for all cycles, and the latter two-thirds of each cycle is identical. Cycle 4 shows steep pressure drawdown to the minimum cycle pressure and gradual buildup to maximum cycle pressure. Cycle 6 represents gradual pressure drawdown followed by steep pressure buildup. Cycle 5 shows the intermediate case. Formation and fluid properties used in the calculations are the same as those employed for pressure cycles 1, 2, and 3, and are again given in Figure 8. The calculated V/Vo curves of Figure 8 corresponding to cycles 4, 5, and 6, show that cycle 6 gives maximum water encroachment into the gas field, and cycle 4 causes least amount of water encroachment. Again, results show failure of the gas reservoir volume to assume its original value at the end of the cycle, the decreased volume for 10-11 months of the 12-month cycle, and growth of the yearly average gas reservoir volume with time (Figures 9 and 10). Figure 12 shows the effect of the time of conversion from production to storage operations and of injection of gas as opposed to withdrawal of gas during the initial portion of the storage cycle. Figure 10 gives pressure cycles 7, 8, and 9, showing initiation of storage operations after a four-year production period, initiation of storage operations immediately upon discovery with a production-injection cycle, and initiation of storage operations upon discovery with an injectionproduction cycle, respectively. The periodic gas field pressure

1.080 - Cycle 4 1.060 -— <~ —- Cycle 5 1.060 [04 ~ Cycle 6 1.040 1.020 1.000 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0.980' 0.960 i it ii II~~~~~~~~~~~~~~~~~~~~~~~~~~~I 0.940'I I'9M I';~~ ~~~'I Ii II I I t 0.900'- I " 0.880 0.860 - -, Q820 II 0.800. 0.780 1 2 3 4 5 Y 617 18 19 ~ IMEIN YEARS Effect of pressure cycle on gas reservoir volume variation Figure 10

-35in order to allow mathematical treatment of aquifer water movement, they are seldom realized in an actual gas reservoiraquifer system. Actual aquifer formations have varying degrees of non-uniformity and gas field boundaries rarely conform closely to circular geometry. Thus the predictive accuracy of the equations of Section III-A, based on assumptions of uniform aquifer sand porosity and permeability and radial aquifer water flow, might well be questioned. The only practical method of evaluation of the usefulness of these equations is comparison of the observed behavior of an actual gas field with that predicted from the equations. The Section III-A Equations (III-20) and (III-21) have been applied in predicting the performance of a Michigan gas field ("Field A") which was discovered in 1941 and produced for six years prior to conversion to storage operations. The pressure and gas production-injection data from Field A and other data used in the calculations are listed in Table 1. The monthly wellhead pressures are shut-in pressures and are intended to represent an average uniform gas field pressure. Figure 13 is a sketch of the gas field boundary and shows the non-circular geometry of the gas sand. The gas injection-production data from Table I and Equation (III-21) were employed in calculating the predicted gas bubble pressures, Pt. These pressures were then employed in Equation (III-20) to obtain the predicted volume ratio Vt/Vo.

1.6 Monthly Cycle 7 Cycle 8 Cycle 9' 1.2- Increments =' 0.8 1r' 0.4 - J I r=y () I 2 3 4 5 0 1 0 1 TIME IN YEARS Fig. II Gas Reservoir Pressure-Time Cycles 1.010 1.00, i ).9995 0.'.990~- ~p0 =440 PSIA.985 h 2\ Vo= 5 x 107 ft3 s.980 c 3.3 x 10-6 1 PSI: ~:~0=0.1.975- k/PL - 160 Millidarcys /cp 970 h= 20 ft.970.965 —-. Cycle 7 ~.960 -— o —Cycle 8.955 --— Cycle 9.950.9465 I I I * 0 1 2 3 4 5 6 7 8 9 10 TIME IN YEARS Gas field volume variation with time Figure 11

-36Table 1 GAS FIELD A DATA Pressure Pressure Base = 15.025 psia Base = 15.025 psia Cumulative Gas Pwellhead, psig Cumulative Gas Pwellhead, psi6 Date Production Mcf at End of Month Date Production Mcf at End of Month Initial 438 1945 Aug. 6,141,991 2>4 1941 Sept. 152,222 Sept. 6,180,187 500 Oct. 305,413 Oct. 6,548,855 Nov. 504,643 Nov. 6,510,325 -79 Dec. 783,967 Dec. 6,679,127 1942 Jan. 1,043,989 419 1946 Jan. 6,849,875 275 Feb. 1,281 223 410 Feb. 7,001,380 Mar. 1,517,966 402 Mar. 7,149,025 Apr. 1,708,094 397 Apr. 7,290,974 289 May 1,814,486 394 May 7,438,005 2559 June 2,001,371 591 June 7,573,120 255 July 2,199,592 589 July 7,689,770 lt8 Aug. 2,401,853 Aug. 7,821,790 247 Sept. 2,679,791 Sept. 7,927,062 2599 Oct. 2,964,208 Oct. 8,007,537 239 Nov. 3,237,656 361 Nov. 8,010,657 241 Dec. 3,388,508 359 Dec. 8,013,657 243 1943 Jan. 3,541,680 1947 Jan. 8,015,657 240 Feb. 3,671,511 349 Feb. 8,017,657 237 Mar. 3,804,659 349 Mar. 8,019,657 236 Apr. 3,976,316 336 Apr. 8,021,657 236 May 4,187,552 May 8.022,657 245 June 4,372,740 328 June 8,023,657 2459 July 4,511,098 325 July 8,024,657 249 Aug. 4,652,784 318 Aug. 8,025,657 252 Sept. 4,803,284 315 Sept. 8,026,657 251 Oct. 4,991,379 313 Nov. 5,171,583 311 Dec. 5,306,024 08 Base = 14735 psia 1944 Jan. 5,417,306 304Base 14735 psia Feb. 5,534,524 298 Mar. 5,662,943 289 Net Gas* Pwellhead, psij Apr. 5,732,452 289 Date Stored, Mcf at End of Month May 5,842,514 294 June 5,950,056 290 Junely 6,022,81 280 1947 Sept. 20,898 1 Aug. 6,032,458 286 Oct.,95025 Sept. 6,047,028 288 Nov. 130,547 252 Oct. 6,076,636 90Dec. 146,190 252 No 6115,46 2198 Jan. 156,194 255 Nov. 6,115,460Feb. 200,560 254 Dec. 6,137,991 MaFeb. 265,049 5 1945 Jan. 6,138,991 280 Mar. 6,o49 258.5 4Feb5. 6,139,991 Apr. 540,097 75.0 Febar. 6,149,991 May 618,588 275.8 Mar. 6,141,991 June 723,296 279.2 Apr. 6,141,991 July 943,431 287.5 JuMayne 6,141,991 Aug. 1,127, 794 2955 July 6,141,991 Sept. 1,211,634 296.2 Oct. 1,265,670 297.7 Nov. 1,294,407 C98.6 Dec. 1,296,557 298.8

-37Table 1 (CONT'D) GAS FIELD A DATA Pressure Pressure Base = 14.735 psia Base = 14.735 psia Net Gas* Pwellhead, psig Net Gas* Pwellhead, psig Date Stored, Mcf at End of Month Date Stored, Mcf at End of Month 1949 Jan. 1,250,493 297.3 Feb. 1,148,212 293.4 1953 Jan. 10,821,652 521.0 Mar. 1,106,796 291.8 Feb. 8,969,990 475.0 Apr. 1,079,054 291.1 Mar. 7,821,885 441.6 May 1,094,647 291.9 Apr. 7,725,988 440,0 June 1,116,795 293.9 May 8,516,559 454.5 July 1,149,263 295.1 June 8,558,477 465.4 Aug. 1,186,690 297.1 July 11,222,572 533.1 Sept. 1,214,451 297.8 Aug. 15,574,474 594.0 Oct. 1,171,066 296.1 Sept. 14,554,598 616.6 Nov. 1,093,240 292.3 Oct. 15,176,916 627.0 Dec. 1,0o11,018 289.4 Nov. 14,321,056 605.5 1950 Jan. 927,795 286.9 Dec. 12,264,255 549.4 Feb. 854,447 283.5 1954 Jan. 9,381,780 476.5 Mar. 768,585 280.6 Feb. 7,816,040 433.1 Apr. 686,403 277.9 Mar. 6,923,033 407.6 May 823,811 284.6 Apr. 7,229,710 417.2 June 1,005,054 292.6 May 8,425,373 44.8 July 1,114,710 296.3 June 9,943,134 494.4 Aug. 1,151,136 299.5 July 11,417,86 51.2 Sept. 3,316,323 360.5 Aug. 13,298,400 577.6 Oct. 5,802,859 436.2 Sept. 14,383,519 6.2 Nov. 5,999,257 444.4 Oct. 14,309,735 598.0 Dec. 5,786,152 4534.4 Nov. 12,699,370 558.5 1951 Jan. 3,555,830 367.5 Dec. 9,132,326 467.0 Feb. 2,730,345 338.3 1955 Jan. 6,519,597 396.2 Mar. 2,509,993 331.7 Feb. 4,771,604 332.4 Apr. 3,771,416 369.1 Mar. 3,809,960 333.0 May 5,546,526 424.7 Apr. 5,735,442 529.0 June 8,551,702 504.5 May 3,695,422 329.0 July 10,037,144 543.3 June 6,922,929 409.8 ~Aug. 12~,145,655 596.1 July 10,426,744 497.7 Aug. 12,143,655 596.1Au. Yo 6 Sept. 13,356,328 622.4 Aug. 1,208,806 565.2 Oct. 13,985,542 634.4 Sept. 13,816,648 581.4 Nov. 12,237,623 584. Oct. 15,815,026 580.8 Novc. 12,257,625 5849.2 Nov. 10,923,240 520.0 Dec. 11,158,958 549.8 Dec. 422,551 451.4 1952 Jan. 8,869,235 492.4 J.,1,51 5. Feb. 7,832,494 463.4 F95. 22,5 375.5 Mar. 7,675,820 457.4 MFeb. 2,923,306 312'.8 Apr. 8,201,524 470.2 Mar. 1866,177 275.0 May 10,217,490 524.8 Apr. 1,978,250 27.7.4 June 12,459,120 582.5 JMay 2,333,641 276.6 July 14,909,489 641.9 June 3,646,166 329.4 Aug. 15,801,603 661.0July 5,434,80 38.4 Sept. 16,178,353 664.9 Aug. 7,592,756 434.5 Oct. 15,525,408 644,7 Sept. 9,501,879 502.5 Nov. 14,671,019 621.1 Oct. 9,451,169 483.2 Dec. 13,581,502 587.7 Nov. 10,140,934 497.8 Dec. 8,693,064 466.8 1957 Jan. 5,640,184 392.4 Feb. 4,336,604 356.8

TABE I (CONT'D) GAS FIELD A DATA Pressure Base 14.735 psia Net Gas* Pwe.lh.eai, psig Date Stored, M f at, Endd of Morn!th 1957 Mar. 4,862,678 360.9 Apr. 5, 987,497 389 7 May 7,381,146 427.5 June 9,291,480 476.6 July 12,303,629 546.8 Aug. 15,526, 498 622~ 8 Sept. 15,582, 23 625 2 Oct. 15,309,240 616o0 Nov. 14,353,834 595~ 7 Dec. 12,104,393 542.0 Depth of Gas Field Below Surface 1100 ft Gas Field Initial Volume 4.6 x!08 f,3 Aquifer Thickness, ft 20 Formation (Gas Field and Aquifer) Porosity.20 Gas Gravity 0.70 Total Aquifer Compressiblf-Lit-y Vol/Vol-Psi 7 x 10 Permeability of Formation, md 956 Water Viscosity, cp Forration. Temperature, OF 7C Gas Bubble Radius, ft 3020

-39W WATER GAS Plan Sketch of Gas Field A Boundary. Figure 13

The observed pressures of Table I, converted to psia and corrected for depth, were employed in Equation (II-20) to calculate the actual gas field volume ratio. The calculations employed reservoir constants noted in Table I and were programmed on the IBM 650 digital computer. Monthly time increments were used in the calculation. Figure 14 shows the agreement between the predicted and observed gas field pressure ratios, P/Poo Figure 15 permits a comparison between observed (calculated from observed pressures) and predicted volume ratios, Vt/Vc, each plotted as a function of time. Calculations were performed for several different assumed values of rDe, the ratio of the aquifer exterior radius, r5e to interior radius, rb. Results for rDe = 20 are shown in Figure 14, and results for rDe = 10 and rDe = 20 are shown in Figure 15. The value of rDe is assumed to be'approximately' 20 because the dimensionless production quantities, QtDj given by Van Everdingen and Hurst(l) are not tabulated for rDe > 10 for their'constant terminal pressure' case of unsteady state radial liquid flow. Thus the values of QtD employed in Equation (III21) for an rDe > 10 had to -be extrapolated from the QtD values given for values of rDe < 10o Equation (III-21) was also solved for the case of an infinite aquifer, but results were in poorer agreement with data than were the results for a finite aquifer with rDt = 20.

1.500 Field A Estimated original pore volume of gas sand V0= 4.6 x 100 ft3 I 400 Estimated Itk for aquifer= 0.956 darcy /cp 1.300 * Points calculated from field data A Points predicted by theory based on finite aquifer RL120 1.200 Z 1.100 a. 1- l.00C Ca 1 0.9000.700t- ~ h, ~ aA 1 Points calculated from observed data 0. 6000.5001 t - o 2 3 4 5 6 7 8 9 1 0 1 1 2 1 3 1 4 1 5 1 TIME IN YEARS Comparison Between Observed and Predicted Gas Storage Reservoir Well Head Pressures Figure 14

1.12 PoField A Estinreted original volune of gas sand Vo 4.6 x 1083 ft I. Pit Predite d for R =10 Estimated for aq uifer.956 d arc n cp ve d a Ht I _ o o oI 1.00', —0.J vuo0 E.94 O 00000000000000 000T I _ - - ~~.uA Calculated Data f A.88 aA,A, n\ o r Predicted for R 20o(estimated) a Points calculated from field data o Points predicted by theory, based on finite aquifer R 10 a Points predicted by theory, based on finite aquifer R. 20 (estimated)beyond values given by Van Everdingen and Hurst 0 1 2 3 4 5 6 7 8 9 IO 1I 12 13 14 15 16 17 TIME IN YEARS Comparison Between Observed data and Predicted Gas Storage Reservoir Volume Variations Figure 15

C. Determination of Effective Reservoir Parameter Values from Gas Field Pressure-Production Data The parameters appearing in the equations of Section III are the aquifer formation thickness h, porosity cp, and mobility k/t the aquifer compressibility c, and the gas bubble radius rb. The'effective' values of these parameters are defined as those values which result in closest agreement between the predicted and actual gas reservoir performance over the period of time for which pressure-production-injection data is available. Examination of the expression relating tD to t (page 12 ) and the previously presented Equations (III-21) and (III-22) shows that the reservoir parameters appear only in the combinations K1 = thmcrb2 (in Equation III-21), K2 = k/cpcrb2, (in expression for tD) and K3 = [i/2ithkAt (in Equation III-22). Thus the values of K1 or K3 and K2, not the individual values of h, k/ cp, c, and rbp are unique for a gas reservoiraquifer system. One could, for example, obtain mathematically identical Pt values from Equation (III-21) by employing infinitely many different combinations of h, p,, c, rb, and k/t values provided one kept constant the values of K1 and K2. Therefore the determination of the effective reservoir parameter values reduces to the evaluation of K1 and K2 [if Equation (III-21) is employed] or K3 and K2 [if Equation (III-22) is employed]. The condition that the assumptions on pages 8 and 9 are satisfied by a specific reservoir-aquifer system is, of course, implicit in the claim of uniqueness of the K1 (or K3) and K2 values for that system. Difficulty might be expected in attempts to determine the effective

-44K1 or K3 and K2 values for an'aquifer storage' gas reservoir since the assumption of constant rb is certainly not satisfied in this case. Thus the method of K1 or K3 and K2 determination is discussed below first in relation to reservoir-aquifer systems which, in general, satisfy the assumptions on pages 8 and 9, and second in relation to'aquifer storage' gas reservoirs. Section III-C-3 below presents an analysis of the errors in the calculated Pt or Vt values caused by errors in the values of the reservoir parameters. 1. Case 1 - Gas Storage Reservoir-Aquifer Systems Which Satisfy the Assumptions in Section III-A The phrase'satisfaction of assumptions' is not to be taken literally in the present context since no actual gas field-aquifer system will conform perfectly to the assumptions made in the development of the equations of Section III-A. No precise statement can be made concerning the extent to which a given reservoir-aquifer system can deviate from the assumptions and yet be amenable to analysis by the method described below. Certain extreme cases such as aquifer storage reservoirs and reservoirs undergoing interference effects can certainly be excluded from the category of systems amenable to treatment by the method discussed below. The term DEV, defined in Equation (1II-25) below, serves as a measure of the agreement between the actual and predicted pressure behavior of a gas storage reservoir; DEV is the average fractional deviation between the observed and predicted

-45gas bubble pressure over the time period for which field pressure data is available. Its value wills of course, decrease with increasing agreement between these pressures. DE L/ t=l ( t (III-25) PiAt = Pt = predicted gas bubble pressure at time t - iAt, calculated from Equation (III-21) or (III-22), psia PaiAt = observed gas bubble pressure from field data at time t = iAt, psia At = time increment chosen for calculations, seconds N = number of time increments comprising time period for which a comparison between predicted and observed pressures is desired The effective K1 (or K3 ) and K2 values for a reservoir-aquifer system can now be determined by analysing the results of the following sequence of calculations: 1) assume a value of K2 2) assume a value of K1 (or K3) 3) calculate Piat (i = 1, 2, 3,. o. N) from Equation (111-21) or (III-22) 4) calculate DEV from the observed gas bubble pressures and the pressures calculated in step 3) 5) assume a different value of K1 (or K3) 6) repeat steps 3), 4), and 5) until a plot of the DEV values calculated in step 4) versus K1 (or K3) shows a minimum DEV value -47

-467) choose a different K2 value and repeat steps 2) 6) The effective reservoir parameter values are now the K1 (or K3) and K2 values corresponding to the minimum DEV of the several minimum DEV values referred to in step 6). The above-outlined calculations are illustrated by application to an ideal gas field-aquifer system having the following characteristics: infinite aquifer K1 = thcpcrb = 10 ft3/psia K2 = k/cpcrb2 = 10-5 seconds-1 Vo = initial gas bubble pore volume = 220,000 ft3 A time increment of one month was chosen (At = one month = 2.628 x 106 seconds) and the calculations were performed for 13 time increments (N = 13). The gas field pressure at t = At, 2At,... 13At was specified as shown in Table II. The pore volumes of the gas bubble were calculated on the IBM 650 computer from Equations (III-15) and (III-17) for K1 = 10, K2 10-5 and are listed in Table II. These pore volumes were then considered as reflecting the actual or observed behavior of the gas field and the effective values of K1 and K2 were considered unknown. Since the pressures were specified and the pore volumes were calculated, DEV in this case was defined as in Equation (III-25a) below.. t(III-25a) /3

-47TABLE II PRESSURE-VOLUME BEHAVIOR OF IDEALIZED GAS STORAGE RESERVOIR Specified Gas "Observed" Pore Volume Reservoir Pressure, (Calculated for K1 = 10, Month Psia K = 105) Mcf 0 500.22000 1 460.21395 2 420.19755 3 380.17296 4 430.15437 5 480.15178 6 530.16066 7 600.18291 8 6oo.20944 9 54o.22213 10 480.21765 11 440.20281 12 500.19606 13 60oo.21382

-48Viot = predicted gas bubble pore volume at time t = iAt, ft3 Vai&t = actual gas bubble pore volume at time t = iAt, ft3 For each of several assumed values of K2, the predicted pore volumes Viot and the corresponding DEV values were calculated for several assumed K1 values. These calculations were performed in the manner outlined in steps 1) - 7), above. The calculated DEV values corresponding to the various combinations of K2 and K1 values are listed in Table III and DEV is plotted versus K1 with K2 as a parameter in Figure 16. Figure 17 shows the observed pore volumes [calculated from Equations (III-15) and (III-17) with K1 = 10, K2 = 10-5] and predicted pore volumes [calculated from (III-15) and (III-17) for K1 =.0229, K2 = 10-2, and for K1 = 60, K2 - 10-6] plotted versus time. Figure 16 reveals that, for a fixed K2 value, DEV is sensitive to the assumed K1 value. However, DEV is less sensitive to the assumed value of K2 provided that for each K2 value the gas field pore volumes and DEV are calculated for the optimum K1 value, i.e., the K1 value corresponding to points a, b, c, and d noted on Figure 16. If one performed the calculations indicated in steps 1) - 7) by fixing a K1 value and then calculating the pore volumes and DEV values for a sufficient number of K2 values to obtain a minimum on a plot of DEV versus K2, then the results would show that DEV is very sensitive to the assumed K2 value for a fixed K1 but is relatively insensitive to the

10 9' 0 6 K 2 Ka IO 6 0 0 w b 0 C 0 pn QOIO 0.015 Q020 Q025 0.10 0.15 0.20 I 1.2 1.3 1.4 L5 10 60 O 80 9 Effect off Assumed K1 and K(2 Values on Agreement Between Predicted and CObserved Gas Bubble Volumes.Fiue1

-50230 /I! \, 2 10I~~~~~~~~~ / 2w~~~i\, /0 0, 180 I 170 \ fr —* —-— A Kx0.O229 K IO 160, \ /00 / 2 150 0 2 4 6 8 10 12 NUMBER OF ONE - MONTH TIME INCREMENTS'Observed' and Predicted Gas Bubble Pore Volumes for Idealized Reservoir-Aquifer System. Figure 17

-51'TABLE III EFFECT OF K1 AND K2 ON DEV VALUE FOR IDEALIZED GAS RESERVOIR - AQUIFER SYSTEM 2 2 K2 =1 /thcrb DEV 10-2 2.292 x 10-2 2.71% 10-3 1.886 x lo-l 2.12% 10-4 1.407 1.37% 10-5 1.000 x 10 0 10 -6 5.987 x 10 2.33% 10~7 1.748 x lo2 9.59%

assumed KEm value provided K2 is taken equal to its optimum value for the assumed K1. Thus, in short, DEV is relatively sensitive to the assumed value of K1 (or K2) when K2 (or K1) is fixed but is relatively insensitive to the fixed value of K2 (or K1) provided the optimum K1 (or K2) is chosen for each of the latter fixed values. 2. Case 2 -'Aquifer Storage' Gas Reservoirs The validity of the method of K1 (or K3) and K2 determination discussed in Section III-C-1 above depends upon'randomness' in the errors in the Vauit and Viot or PaiAt and PiAt values. That is, differences between the observed and predicted performance of the gas field should be due to some non-homogeniety in the aquifer formation, perhaps slight deviation of the field boundary from a circular shape, slight variation in rb, and errors in the reservoir pressure and gas flow rate measurements. Differences due to these factors can be reasonably expected to be much more random than the difference between PaiAt and Pint in the case of an'aquifer storage' gas reservoir where, while PiAt is calculated from equations derived from the assumption of a constant rb, rb actually increase from a value of 0 to values of several thousand feet as the gas bubble is grown. Analyses involving the Section III-A equations will now be made to (1) determine the effect of an increasing rb on the actual or observed performance of an aquifer storage gas

-53reservoir,(2) determine the manner in which K3 and K2 must be chosen when employing the method discussed in Section III-C -1, above, in order to produce the same trend in the predicted field performance as the increasing rb produces in the actual field performance. The results of these analyses are checked by applying the method of Section III-C-1 to an actual aquifer storage reservoir. The variable rate Equation (III-16a) is employed in these analyses and the assumption is made that the term (2Vi - Vi+l - Vi-l) in this equation is constant and less than 0. [The following analyses would yield the same results if this term were allowed to take on varying negative values; however the analyses are much simpler if (2Vi - Vi+l - Vil) is assumed constant.] This is equivalent to asserting that the second derivative of the pore volume with respect to time is greater than O. This latter assertion is in general valid for an aquifer storage reservoir during the initial growth stage; Figure 21 shows that the actual pore volume-time curve for aquifer storage reservoir B conforms very well to this assumption during the first 31 weeks of the gas bubble growth. The following analysis shows that the effect of a growing rb on the performance of an aquifer storage reservoir is a less rapid increase of gas bubble pressure with time than would occur in a reservoir of constant radius. The last term of Equation (III-16a) can be included in the summation term to

yield L- -/ where AVi = 2Vi - Vi+l - Vil; a negative constant, from the above assumption. The effect of an increasing rb on Pj can be determined by differentiating both sides of (III-26) with respect to rb. s ar, i - v dtp drbt j-)t (III-27) Now, since tD = K.2t = k/Acperb 2At 9 dt _ 2 kt t Jt S.+'c1f r3: and (III-27) becomes Thtrm d~hD/dt isL = po site i (III-28) The term tD dPtD/dtD is positive for all values of tD12. aPj/6rb is therefore negative since AVi is negative and K2 and K3 are positive. Thus as rb increases, Pj will assume lower values than those corresponding to the case of a constant rb. Nowy since the actual gas bubble pressures reflect the water movement away from a growing gas bubble and the predicted pressures are obtained by calculations in which rb is assumed constant, one might expect the actual pressures to assume smaller values (smaller than the predicted pressures) as time increases, since aPj/ rb < 0. Thus application of the method of Section III-C-1, above, to an aquifer storage

-55reservoir should yield increased agreement between the reservoir observed and predicted performance for combinations of K3 and K2 values which effect a smaller slope of the predicted pressure-time curve. The statement that the predicted pressure-time curve slope will be decreased by increasing K2 and decreasing K3 will now be proven. The slope, Sj, of the predicted pressure-time curve can be obtained by differentiating both sides of (III-26) with respect to time and noting that d/at = K2 a/atD (since tD - K2t). cj = at = -3k L (J (III-29) The total differential of Sj with respect to K3 and K2 is S/j= a- S kci (III-30) or, dealing with finite increments, as- = k; 13 A kt (III-31) The terms aSj/aK3 and aSj/aK2 can be obtained by differentiating (III-29) to yield \SJ = (-b3K2V L(cltD, i-;) D ) 3 (AC t'24 < E ( t W Kkv,> (to Dot $ apt ~ 1<3 f dViP / dt) "

-56or, upon rearrangement, + (5 (t3 a -.-LtJ The reader might interject at this point that a negative ASj (decreased slope) is easily obtained by simply decreasing K3 (negative 6K3/K3) or K2 (negative AK2/K2). However, the answer is not that simple since if K3 is decreased then K2 must be increased, and vice-versa9 to maintain close agreement between the predicted and actual gas bubble pressures. For suppose that for a certain K3 and a certain K2, a close agreement is achieved between actual and predicted pressures and it is desired to change K3 and K2 so as to achieve closer agreement. Then by decreasing K3 and leaving K2 unchanged, or vice-versa, one will simply effect decreased values of all the predicted pressures9 as is obvious from Equation (III-26). But, in general, one desires to change the slope of the predicted pressuretime curves at the same time keeping the general magnitude of the predicted pressures the same. Actually, if K3 is decreased by a certain factor then K2 must be increased by a larger factor in order to maintain close agreement between the predicted and actual pressures. This fact can be seen from Figure 18 and Table VI, which present the results of application of the method of Section ITIC-1 to an actual aquifer storage reservoir.

-57These results will be discussed below. (Table VI shows that as K3 is decreased, the product K3 x K2 is increased; that is, K2 is increased in greater proportion than that in which K3 is decreased.) Equation (III-33) indicates that a negative ASj (decreased slope) can be obtained by decreasing K3. Since the product -K2K35AVi in (III-33) is positive (AVi < 0) the term in the square brackets must be negative for LSj to be negative. Approximate values of the terms dPtD,/dtD and d/dtD (tD dEtD/dtD) are tabulated in Table IV. These values show that if tD > 90, the latter term is negligible compared to the former. Since these two terms are the coefficients of AK3/K3 and AK2/K2 in (III-33), the sign of ASj is primarily dependent upon the sign of AK3/K3. Thus if K3 is decreased, ASj will be negative; a decreased predicted pressure-time curve slope is therefore obtained by decreasing K3 and increasing K2. The agreement between predicted and actual aquifer storage reservoir performance (during initial growth stage) should therefore increase as K3 is chosen smaller and K2 larger. The method described in Section III-C-1 has been applied to aquifer storage reservoir B in order to check the above conclusions concerning the effects of an increasing rb and the manner of K3 and K2 variation on the field observed and predicted performance. Field B is an aquifer storage reservoir which was initiated in 1957 by injection of natural gas into an

-58TABLE IV VALUES OF PtD AND RELATED DERIVATIVES VS t -tD D dPt d( ddPt PpD t d VD D tD tD ptD dtD dtD dtD 8 1o556 9 1,604 o048.432 i0 1.651 o047 0470 ~038 40 2,282 50 2,388,olo6.53 60o 2,476.oo88 528 -,0002 8o 2,615 90 2,762.0057 513 100 2,723 00051,510 -.0003 400 3 406 500 35516.0011.550 6oo 360o8.0092.5515.000015 1000* 3 860* 00005* o 5* 0* * For tD > 1000o, PtD 1/2 [~8090 + en tD] (reference 1) d PtD ~dt 2t 1.0005 for tD = 1000 dtD t d d PtD ~tD=.5 and. (tD t dtD dtD dtD

-59aquifer formation having an initial uniform pressure of 1060 psiao This field is an excellent example of a gas storage reservoir which does not satisfy the assumption of a constant rb. since Field B's radius increased from 0 feet on November 30, 1957 to approximately 2600 feet in early summer of 1958. Available pressure-injection data for the period November 30, 1957 - July 5, 1958 were employed in the determination of the K3 and K values which would result in closest agreement between the observed and predicted gas bubble pressures. The predicted pressures were calculated from Equation (III-22), with the required nt values obtained from the data tabulated in Table V. These predicted pressures and the nt data were employed in Equation (III-20) to calculate the predicted pore volumes. The actual pore volumes were calculated by inserting Pt and nt from the data into Equation (III-20). A measure of the agreement between the actual and predicted field B performance was provided by the value of DEV, defined in Equation (III-25) (with N = 31). The calculations of the predicted pressures were performed with a time increment of 7 days and with various combinations of K3 and K2 values chosen as outlined in the method of Section III-C-1. Table V lists the observed field B data and the actual and predicted pressures and pore volumes (for K3 = 2.3 x 10-5 ~ = 1000) for the period November 30, 1957 - July 5, 1958 (31 weeks). Figure 18 shows DEV plotted versus K3 with K2 as a parameter. These curves show that the agreement between the

-60o TABLE V SUMMARY OF FIELD B ACTUAL AND PREDICTED PERFORMANCE FOR PERIOD NOVEMBER 30, 1957 - JULY 5, 1958 MMcf In-Situ Gas* Reservoir Gas Bubble Date @ 14. 65 Psia & 60cF Pressure* Psia Pore VQlume*, MMcf 1957 Novo 30 0 1060 d O Dec, 7 5.092d 1075 d 1074.3p.059d o059p Dec,.14 8,831d 1075.4d 1071. lp.1025d.103p Dec. 21 15.747d 1088o2d 1079.7p.180d,182p Dec, 28 20o722d 1084.9d 1075.2p.238d.241p 1958 Jan. 4 21. 697d 07;8,o8d 1064o6p,251d.255p Jan. 1i 29.254d 1093.2d 1080,9p.335d o 338p Tan_ 18 39 o080d 1098.5d 108709P.442d.448p Jar.n 25 50,330d 1100.7d 1092o3p,568d.574p 1958 Feb. 1 69 512d 1206.9d 1111 2p.702d.776p Feb 8 1.09 747d 1222,4d 1157,5p 1.09 d 1.17 P Feb. 15 10O.990d 1103.9d 1089o9p L 25 d 1.27 P Feb, 22 -11O990d 1088,5d 1070 9p 1.27 d 1.30 p 1958 Mar, 1 146.886d 119353d 1136.2p 1.50 d 1.60 p Mar. 8 191 o112d 1204.6d 1165 9p 1.93 d 2.01 p Mar. 15 235.7 d 1202,8d 1174.9p 2.39 d 2.46 p Mar. 22 254o715d 166.8d 11.40o2p 2.68 d 2.76 p Mar, 29 310 o842d 1203.6d 1181.1p 35.15 d 3522 p 1958 Apr. 5 316o410d 1202,4d 1189.1p 3.67 d 3572 p Apr. 12 403,439d 1.191,4d 1182o5p 4.14 d 4.18 p Apr. 19 459.657d 1201.7d 1196,3p 4,67 d 4.69 p Apr, 26 481,326d 1145 d 1165.1p 5,1.8 d 5.07 p 1958 May 3 497.766d 1.!43.2d 1143,5p 5,37 d 5.37 p May 1.0 558 301d 1181.2d 1178.0p 5.79 d 5,81 p May 1.7 617.53 d 11934 94o7P 6,32 d 6.31 p May 24 657o215d..180o8d 1186.4p 6.82 d 6,78 p May 31. 714o131d 1194,7d 1196.9p 7.30 d 7.29 p 1958 Julre 7 784.039d.20o8,8d 1214.p 7.90 d 7.86 p June. 14 850o387d 12070 d 1221.o5p 8,59 d 8,46 p J'une 21 925,87 d 1190o3d 12335.p 9.51 d 9.1 p June 28 1.02442 d 238d 1 9d 1256,6p 10.02 d 9084 p 1958 JL.?y 5 1.116565d 1242o7d 1267o2p 10o88 d 10.62 p * d dernotes guantity obtained from field data. * p denotes quantity predicte+ for K35 = 2,5 x 10-5 O = 1000,

2.1 k 4 2.0 K2 KI c r.73x K' 2r hk At 1.9 At 7 7 Days - 0.606 x 10 Sec s. b K2 0. I K2 1000 c d 1.6 0.2 0.3 0.4 0.5 0.6 0.7 K3 x 10 Dev vs. K3 for Aquifer Storage Reservoir B. Figure 18

-62predicted and observed field B pressures improves as K3 is decreased and K2 is increased. Thus application of the method of Section III-C-1 to field B does not yield a unique pair of K3 and K2 values corresponding to closest agreement between predicted and actual performance but yields increasingly better agreement as K3 is decreased and K2 increased. Table VI lists the K3, K2, K3 x K2, and DEV values corresponding to the points as b, c$ and d noted on Figure 18. It is of interest to compare the K3 and K2 values in Table VI with the following approximate K3 and K2 calculated from available data relating to h, (p, c, k/4, and rb: tk:= 3a2 = 9.55 xl6 /O r p/f(A at = o.cog X /O)6 k2 = k9 X - OjX/O-S sec's-1 The product K3 x K2. tabulated in Table VI, varies widely with K2 (or K3). If K2 x K3 were equal to a constant A, then the calculations involved in the method of Section III-C-1 would be greatly reduced since for every value of K2 assumed, DEV would have to be calculated for only one K3 value, namely K3 -A/ Figures 19 and 20 compare the actual and predicted field B pressures as a function of time for pairs of K3 and K2 values corresponding to points d and b of Figure 18. The actual pressures are seen to fall below the predicted pressures as time increases; that is, the'average' slope of the predicted pressure-time curve becomes greater than the'average' slope

1280 P 1260 - - OBSERVED PRESSURES -0 —0-0 PREDICTED PRESSURES 1240 - - K3' 2.3 x 10 1220 -C1000 I/. 1200 w:C 1180 I w~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ w V) 1 16 0 P 0)lID \/ r,,o. 0 1140 - > ~ ~ ~~~ ~~~~~~~~~~~~~~~~~~~~> w Ld ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-J uJ u" 1080 Io \I Cn 1120 W w 1100 1080 PI'P 1060 1040 0 2 4 6 8 10 12 14 16 Is 20 22 24 26 28 30 32 jNUBER OF 7-DAY TINE INCREMENTS Field B Observed and Predicted Pressures vs. Time, November 50, 1957-July 5, 1958. Figure 19

1280 OBSERVED PRESSURES 1260 / 0-0-0 PREDICTED PRESSURES 1240 C K,~K5.55 x 10 5 K~a.001 1220. 0,C 1200 Co 0, I 4w ~ 160 Co IC) \ C0) 0~~~~~~~~~~~ o 1140 un Ic I O~~~~~~~~~~~~~~~~~~~~~~0-0 w 0 c)r 1120 wx I 5 1140 1100 ~~~~~~/ Ix~~~~~~~~~~~~ 1100O 1080 0\ I / 1060 0 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 jNUMBER OF 7-DAY TIME INCREMENTS Field B Cbserved and Predicted Pressures vs. Time, November 30, 1957-July 5, 1958. Figure 20

-65TABLE VI K3 x K2 AND DEV TABULATED AS A FLUNTCTION OF K2 FOR FIELD B Point on Fig. 18 K2 K3 K2 x K3 DEV, % a 1.73 x 10-4.680 x 10-4 1.18 x 10-8 1.837 ~~b.001.555 x 10-4 5.55 x 10-8 1.782 C.1.382 x 10-4 3.82 x 10-6 1.702 ooo.230o x o4 023-4 1633 d 1000.250 x 10.025 1.655

-66of the actual pressure-time curve, All the pairs of K3 and K2 values employed resulted in an average predicted curve slope greater than the actual curve slope, although the former slope approached the latter more closely as K3 was decreased and K2 was increased. Figures 21 and 22 show the actual and predicted pore volumes as a function of time for combinations of K3 and K2 values corresponding to points d and b of Figure 18. Thus the results of application of the method of Section III-C-1 to field B confirm the earlier conclusions (obtained above by analysis of the predictive equations) that the effect of an increasing gas bubble radius is a less rapid increase of the field pressure with time (relative to a gas field of constant radius) and that the predicted pressures of an aquifer storage reservoir increase more rapidly with time as K3 is decreased and K2 is increased. The important conclusion to be obtained from the results reported in Figures 18 - 21 and Table VI is that application of the method described in Section III-C-1 to aquifer storage reservior B does not yield a unique pair of K3 and K2 values but yields increasingly better agreement between predicted and actual field performance as K3 is decreased and K2 is increased. The results plotted in Figure 18 show, however, that if the value of K2 is fixed, application of the method of Section III-C-1 yields a unique K3 value corresponding to a minimum DEV or best agreement between actual and predicted gas bubble performance.

I I 10 ACTUAL PORE VOLUMES 9 0-0-0 PREDICTED PORE VOLUMES -5 K 3 2.3 x 10 8 K 1000 2 2 w 7 -J 0 > 6 w I 0, 4~~~~~~~~~~~~~~~~~~~~~~~~~~~~4 ap 4 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~In IAJ Cn 3 a - 2 / 0 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 j, NUMBER OF 7-DAY TIME INCREMENTS Figure 21. Field B Actual and Predicted Pore Volumes vs. Time, November 30, 1957-July 5, 1958.

I0 = * * ACTUAL PORE VOLUMES 9 0-0-0 PREDICTED PORE VOLUMES K3= 5.55 X 105 I 8 K:.001 2 - 7 I 6 or 5 03 o /D 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 j,NUMBER OF 7-DAY TIME INCREMENTS Figure 22. Field B Actual and Predicted Pore Volume vs. Time, November 30, 1957-July 5, 1958.

-693. Errors in Predicted Water Influx and Pressure Drop Values Resulting from Errors in Reservoir Parameter Values An ideal reservoir-aquifer system, satisfying all the assumptions on pages 8 and 9, is considered in this section. Errors associated with use of the variable pressure Equation (III-15) are first considered, The correct or actual values of the parameters K1 (= -hcpcrb2) and K2 (=k/cpcrb2) for the ideal system are denoted by (Kl)c and (K2),. (Qt)c denotes the cumulative water influx calculated by inserting the (K1),, (K2),, and pressure schedule APi (considered known) values into the previously presented Equation (III-15); QLL - I -/'t a 2pt =-6)'D (III15) where (j-i)AtD = (j-i)K2At and t = jAt. If (Qt)e denotes the erroneous Qt value calculated by inserting erroneous K1 and K2 values into Equation (III-15), then the fractional error in Qt, EQt, is ( Qt)c (Qt)C Now, since Qt is directly proportional to K1, the fractional error in Qt caused by an erroneous K1 value, (Kl)c + aK19 is simply Eq ) =(III-34) where (EQt)K1 denotes the fractional error in Qt caused by use of an erroneous K1 value and the correct K2 value (K2)c.

-70For small errors in K2q the fractional error in Qt caused by an error in K2 can be obtained from the truncated Taylor series form ((?t),- (t)cQ (-3 Q t a 1K (III-35) where the erroneous K2 value = (K2)c + K2, aQt/aK2 can be determined by differentiating both sides of Equation (III-15) with respect to K2 to yield 94;)tqa C - it. -L)atp (III-56) One thus obtains I( Pt)Q 4ap)I (111-37) where (EQt)K2 denotes the fractional error in Qt caused by use of an erroneous K2 value and the correct K1 value (K1)c, Since tD dQtD/dtD ranges from 0j7QtD* for small tD to O.9QtD* for large tD~ the percentage error in Qt effected by an error in K2 is 07 to 0O9 the percentage error in K2. The fractional errors in Qt. effected by erroneous h, k/7. and pcrb2 values can be determined, by proceeding in the manner employed above, as ( <;t~h = (h)) (III-58) *Thiese figures are calculatedi from tables of QtD vs tD given in References 1 and 2.

-71o -- (k/, ~(E PZp (tg jt)At (YM) (E,)*';O by;,(ji)^t (06? (III-39) (III-40) The Equations (TIII-38)- (iii-40) reveal that an x per cent error in h causes an x per cent error in Qt, an x per cent error in k/i causes an 0.7x (for small tD) to O.9x (for large tD) per cent in Qt, and an x per cent error in pcrb2 effects an 0.3x (for small tD) to O0lx (for large tD) per cent error in Qt. Thus an error in h causes the greatest error in Qt while the same percentage error in k/~ causes a slightly smaller error in Qt and the same error in cpcrb2 effects a significantly smaller error in Qt. An approximate upper limit ( approximate' since all second and higher derivatives in the Taylor series are ignored) on the total fractional error in Qt caused by errors in K1 and K2 is Ev I (Eqtk |4I (f)I - (I')4- 0 9 (III-41) The corresponding approximate upper limit on the total fractional error in Qt caused by errors in h, k/p, and Ocrb2 is

-72_Qt A(k/) (IIi-42) Errors associated with use of the variable rate Equation (III-26) are now considered. The correct or actual value of K3 (= i/2nhkl t) is denoted by (K3)c and (Pj) (- (Po - Pj)0c denotes the pressure drop calculated by inserting the (K3)c, (K2)c, and growth schedule AVi values (considered specified or known) into Equation (III-26). y = TP -q' = k3 =AV, (j1i )^t (111-26) avL = 2 V -V,/ -~L If (Pj)e denotes the erroneous Pj value calculated by inserting erroneous K3 and K2 values into (III-26), then the fractional error in Ppj Epj, is Now, since Pj is directly proportional to K3, the fractional error in Pj caused by an erroneous K3 value, (K3)c + AK3, is simply =E k3 K3 (III-43) where (E )K denotes the fractional error in Pj caused by use of an erroneous K3 value and the correct K2 value (K2)c. For small errors in K2, the truncated Taylor series form similar to that given in (TIII35) can be employed to

-73relate (Pj)e to (Pj)e, The term aPj/aK2 required in this Taylor series form can be obtained from Equation (III-26) as k 3 - B ( d ks i.),t (III-44) 1<2 d.:o Thus tA —--- (III-45) where (E:Pj)K2 denotes the fractional error in "Pj caused by use of an erroneous K value and the correct K3 value (K3)c; the erroneous K2 value = (K2) + AK2' Now, since tD dPtDdtD ranges from 0.5 PtD* for small tD (tD.01) to less than.187 PtD for large tD (tD dPtjdtD =.187 PtD for tD = 100, tD dPtDjdtD <.187 PtD for tD > 100), the percentage error in Pj effected by an error in K2 is 0.5 to.19 (or less) times the percentage error in K2. Thus an approximate ('approximate' since all second and higher derivatives in the Taylor series are ignored) upper limit on the total error in Pj caused by errors in K3 and K2 is z I'd k3j., I _ d _/+ _O___ (III-46) (k.,) 4 -~j(k.), * This figure and the.187 PtD following were calculated from tables of PtD vs tD given in Reference 1.

-74The fractional errors in Pj caused by erroneous h, k/in and perb2 values can be easily determined as (III-47) (EtJ = jt V 2 ( f - (~,u) (III-48) Fill'~.~, L- - V.[ _,i (I I-49) Thus an x per cent error in h causes a -x per cent error in Pj, an x per cent error in k/i causes a -. 5x per cent (for small tD) to -x per cent (for very large tD) error in Pj, and an x per cent error in cprb2 effects a -.5x per cent (for small tD) or smaller (for large tD) error in Pj. D. Effect of Reservoir Geometry on the Calculation of Water Encroachment Into the Gas Field Water influx calculations reported in the literature invariably employ the two assumptions of planar (i~eo two-dimensional) liquid flow and circular geometry of the oil or gas reservoir. These assumptions are perhaps more necedsary ian formulating governing mathematical equations which can be solved. analytically Than rhey are justified by the geometry

-75or characteristics of actual reservoirs. For example, relatively few reservoirs conform closely to circular geometry. However, no solutions to the diffusivity equations analogous to Van Everdingen and Hurst's tabular solutions for the radial flow case(l), have been presented in the literature for the case of a square or elliptic sink (or source) surrounded by a large liquid bearing porous medium. Mathematicians have in fact been unable to solve analytically the diffusivity equation, in one time and two space coordinates, in a planar region bounded internally and externally by concentric rectangles or in a region bounded internally by a rectangle and unbounded externally. Perhaps the most practical and important question is not "Do actual reservoirs conform closely to the assumptions of planar flow or circular geometry?" but is "Does the actual performance of a reservoir which does not conform to these assumptions differ appreciably from the performance predicted by equations derived from the assumptions?". The calculations described in Section III-D-2, below, represent an attempt to answer this latter question with respect to the assumption of circular, reservoir geometry. The performance of an elliptically shaped reservoir, determined by numerical solution of the diffusivity equation in elliptic coordinates, is compared to the performance calculated by approximating the elliptic reservoir boundary as an'equivalent' (e.g., equal area or equal perimeter) circle and employing the radial flow equations. Section III-D-1 presents an analytical solution to the diffusivity equation governing three dimensional liquid flow in an aquifer, i.e., governing the case where the aquifer formation is very thick and vertical flow effects may be significant.

-761. Method of Treating Vertical Pressure Penetration into Aquifer The previously presented Equation (III-5) governs the radial aquifer water flow in the reservoir-aquifer systems sketched in Figures la and lb. Water flow in the neighborhood of a gas bubble situated on top of a very thick aquifer formation, as sketched in Figure lc, however, will not be strictly radial, or horizontal, and will be governed by an equation accounting for a vertical as well as a horizontal pressure distribution. a. Diffusivity Equation Including Vertical Pressure Distribution The diffusivity equation in cylindrical coordinates, with the term 62P/6z2 included to account for the vertical pressure distribution, appears as follows: __P RaP h~ M? _ 9'LP f& r/ R9f U (111-50) where k = permeability in horizontal direction, ft-lbm/ sec2-psi kR = permeability in vertical direction/permeability in horizontal direction z = vertical distance coordinate, ft. Definition of the new variables Akt i, mS

-77P(~s2St) -b tP~() ^(bDyt, ) (III-51) Po(y) = initial pressure with vertical gradient due to gravity alone, psia transforms Equation (III-50) into (III-52) below. + = --- P (III-52) The flow model represented by Equation (III-50) is shown in Figure 23. The vertical distance variable z is set equal to 0 at the upper impermeable plane covering the aquifer and set equal to h, the aquifer thickness, at the lower impermeable plane. The assumptions implicit in Equation (III-52) and made in formulating the boundary conditions, below, necessary for its solution, are 1) constant value of the gas bubble radius, rb 2) uniform aquifer formation directional permeabilities in the horizontal and vertical directions of k and kz ft-lbm/sec2-psi, respectively 3) uniform aquifer formation porosity 4) isothermal, unsteady state aquifer water flow 5) uniform rate of water influx into gas bubble over circular plane portion z = 0, 0 < rD < 1 (the rate may vary with time)

-78UPPER IMPERMEABLE SHALE STREAK GAS BUBBLE Z ~ 0b | -DISCONTINUOUS AQUIFER SHALE STREAKS r, re or Oa Z *h LOWER IMPERMEABLE SHALE STREAK Vertical Pressure Penetration Flow Model. Figure 23

-796) initial aquifer pressure is a function of z only, in accordance with the pressure variance with depth induced by gravity. b. Equation Governing Water Influx into Gas Bubble Darcy's flow law relates the aquifer water velocity across the horizontal boundary z = 0, 0 < rD < 1, to the vertical pressure gradient as shown in Equation (III-53). V - Ad (P- t)] (III-53) where kz = aquifer formation permeability in vertical direction The volumetric rate of water influx into the gas bubble, q, is then _ l (Pe- i (III-54) Now, from Equation (III-51), a P-yt) r A ) - and, employing the relations _P - P - b ke 7 one finds that 3P-(rDYtD) = g7 ( -

a8oor _P__y__tD)_ = __ _ I' P ________ (III-55) Substitution of) from (III-55) into (III-54) yields or ( 9 PSD,;~tD)) - __ 1 ___ (111-56) Thus Equation (III-56) relates the volumetric rate of water influx into the gas bubble to the z-direction gradient of the dependent pressure variable in Equation (III-52). c, Analytical Solution of the Developed Equations for. the Vertical Pressure Penetration Flow Model An analytical solution has been obtained to Equation (IIIS52) for the boundary and initial conditions given in Equations (III-57), (III-58), and (III-59) below., = (III-57) I (111-58) P()g,, i) a) = o (III-59)

The condition (III-57) expresses the fact that the plane z = h is impermeable to water flow. (III-58) represents the condition that the upper limiting plane of the aquifer formation, z = 0, is impermeable with the exception of the circular portion of that plane covered by the gas bubble (see Figure 23). The gradient ( aP t))= is set equal to -1 at this circular portion of the plane z = 0; Duhamel's superposition principle can be employed to transform the solution to Equation (III-52) satisfying condition (II1-58) to a solution satisfying the condition 3~( 9r,,y t) g(tp) ox 1 x @ (III-60) where g(tD) is any semi-continuous function of dimensionless time, tD. The initial condition (III-57) simply equates the initial pressure P(rD,y,O) to the uniform (with exception of the gravity variation with z) pressure Po(Y). The analytical solution to Equation (III-52) for the above boundary and initial conditions is given in Equation (III-61); P(~,~~~i - o t#X (; tt1' - -e (III-61) rn where M = =m - nT/p

-82The reader is referred to Appendix II for the details of the development of this equation. The presented solution, P(rD OtD)g is related to the aquifer pressure at the upper limiting plane z = 0 by the previously given Equation (III-51). Duhamel's principle, discussed in Section III of this thesis, can be employed in conjunction with Equation (III-61), above, to yield the gas bubble pressure corresponding to a variable water influx rate into the gas bubble:.V) If) ) -Trico k,1 Al oItO) (III-62) where Po(O) = the initial uniform pressure throughout the plane z = O0 psi P(lO,tD) = the gas bubble pressure at the circular boundary r = rb, z = 0O psia i i+l i q' = average rate of water influx into gas bubble during the time increment from (i-l)At to iAt, ft3/sec. tD = kt/Acpcrb 2 = dimensionless time t = jAt, seconds At = time increment chosen for purpose of approximating variable water influx rate, q9 by a step function, seconds (lOtD) - given by Equation (III-61) forrD 1

-83Numerical evalutation of the term P(lOtD) has not been performed owing to the lack of data from a gas field conforming to the vertical pressure penetration flow model described above. 2. Numerical Solution of Diffusivity Equation in Elliptic Coordinates The studies reported in the literature of water encroachment into oil and gas reservoirs invariably involve approximation of the reservoir boundary by a circle. This approximation allows use of the solutions to the radial flow equations presented in the literature.(1,2) Nevertheless, the areal boundaries of a majority of the oil and gas fields can be better approximated by ellipses of appropriate eccentricity (> O) than by circles. Thus the question arises as to whether a significant difference exists between the performance of an elliptically shaped reservoir, calculated from solutions to the diffusivity equation in elliptic coordinates, and the performance calculated by approximating the reservoir boundary as a circle and employing the radial flow equations. The present treatment of elliptic flow has been performed to answer this question. The literature presents neither analytical nor numerical solutions to the diffusivity equation in elliptic coordinates. Several texts(14t15a16) discuss the Mathieu functions, which are involved in analytical solutions to this equation, but the authors do not provide solutions useful in applications of

-84interest here. Analytical methods are required to obtain general solutions which can be tabulated in dimensionless form as a function of dimensionless time. However, the analytical approach presents greater difficulties than those encountered in the numerical treatment. The developments and calculations described below are therefore performed only to determine whether there exists any need for general solutions to the elliptic flow diffusivity equation. For if the water flow quantities calculated from the elliptic equation can be duplicated closely by the quantities calculated (for an equivalent circle) from the radial flow equation, then analytical treatment of the former equation would be unnecessary. The results obtained below are plotted in Figure 26 and show that significant error ( 7% in Qt for small tD) can be introduced into the influx quantity, Qtp by approximating the elliptic reservoir as an equal area circle and employing the radial flow equations. Qt is the cumulative water influx (ft3) into a circular reservoir of area equal to that encompassed by the elliptic boundary considered. QtE is the influx into the elliptic reservoirs calculated by numerical solution of the elliptic flow diffusivity equation, as described below. In the development given below, the flow model and assumptions concerning the aquifer formation and fluid properties are identical to the flow model described and assumptions given at the beginning of Section III, with the following exceptions:

-85Radial flow assumptions in Section Corresponding elliptic flow III assumptions The flow model is a thin circular The flow model is a thin disk having concentric circles as elliptical disk having coninterior and exterior boundaries focal ellipses as interior and exterior boundaries rb does not change with time The position of the interior ellipse (gas bubble areal boundary) does not change with time The diffusivity equation in elliptic coordinates is developed below and the equations describing liquid flow across an elliptic boundary are given. A numerical method of solution is then presented and employed to obtain dimensionless water influx quantities corresponding to the dimensionless quantities given by Van Everdingen and Hurst(l) in their "constant terminal pressure" radial flow case. The tabulated elliptic flow results are finally compared to Van Everdingen and Hurst's solution. Error and stability analyses of the numerical method are presented at the end of this section. a. The diffusivity Equation in Elliptic Coordinates The general form of the diffusivity equation governing unsteady state liquid flow through a porous medium is given by the previously presented Equation (III-4).' -- (III-4) I. ~P

-86The terms Pq k/,9,9 c. and t have been defined in Section I'I. The general expression for the three-dimensional Laplacian of a dependent variable in orthogonal curvilinear coordinates u9 vg and w. is(3.7)'7 a (VW 9 I2'- + r:) (a;)7 (III-63) where V (V/(;V) + (V)@; V (III-65) t -@:;7J-f(-4-C. ) (III-66) x,y. z - cartesian coordinates The derivatives required in Equations (III-64)9 (III-65), and (III-66) can be obtained from the following relations between cartesian and elliptic coordinates x = f cosh u cos v (III-67) y = f sinh u sin v (III-68) z = W (III-69) where the foci of the ellipse are at x = + f,~ y - 0, and u, v, and w are the elliptic cylinder coordinates. The geometric relationships between the cartesian planar

-87coordinates x and y and the elliptic planar coordinates u and v are shown in Figure 24. The vertical-direction elliptic cylinder coordinate is denoted by w. Differentiating Equations (III-67), (III-68), and (III-69) and inserting the resulting derivatives into (III-64), (III-65), and (III-66), one obtains be_ f cosh Q - cos v P = O' (III-70) Insertion of these expressions for a, B, and y into (III-63) then yields V'-Z~. PU iC-C ~J (nIII-71) This expression for the Laplacian of P(u,v,w) can be simplified by making use of the assumption (page 8 ) that the flow has no vertical component. This assumption is equivalent to the equality aP)~u~r~ =_ d- ( III-72) where y is the liquid specific weight, pounds force/cubic foot, since there is no liquid velocity component in the vertical direction. Since aP/6w is equal to a constant &2P/aw2 is equal to 0 and Equation (III-71) simplifies to

Geometric Relationships Between Cartesian and Elliptic Planar Coordinates. Figure 24

-89V P(U, V, I)-fs,o L - a S; (III-73) Substituting V2P from (III-73) into (III-4), one obtains the diffusivity equation in elliptic coordinates AgLoC o p (III-74) b. Equations Describing Liquid Flow Across an Elliptic Boundary Darcy's flow equation relates the superficial fluid velocity in a porous medium to the pressure gradient in the following manner, V = - 7 (III-75) where V is the velocity vector, feet/second, y is the liquid specific weight, and w is the vertical direction space coordinate, feet. The gradient of the dependent variable P(u,vw) is expressed in the elliptic coordinates u, v, and w as vpIP = P /aP.)'l:- — 7 a/ Fa (III-76) where i} j, and k are the unit vectors in the u, v, and w directions, respectively. Substitution of the previously determined expressions for a, h, and y into (III-76) yields.& Vcosh5-)a5 - vCOS [i + f -+ 6 ( III-77) PP~~~~ +

Equation (III-78) is obtained by inserting VP, from (III-77),I into (III-75) At any given point on an ellipse u - constant, the term aP/3v j is proportional to the velocity vector component in the v direction, or in the direction tangent to the ellipse at that point. This term can therefore be deleted from (III-78) since only the flow across an elliptic boundary is being considered here. Also, the term (P/Jaw + y/144) k in Equation (III-78) can be set equal to 0 since aP/6w - Y/144, from Equation (III-72). Making these simplifications in (III-78), one obtains V(v t) - / PJUC (III-79) where C = a constant V = V(vt) = fluid velocity in positive u direction across the elliptic boundary u = C, feet/second The volumetric rate of liquid flow across an infinitesimal area element dA at the elliptic boundary is simply VdA or a = _ kS (III-80) d A r~ ap (iii-80o) where dA = hds, feet2 s arec Lengtthi on ellipse u = C, feet

-91q = q(t) = volumetric liquid flow rate across entire elliptic boundary u = C, feet3/second h = thickness of aquifer formation, or height of elliptic cylinder flow model, feet ds is given on an ellipse by Us = 3dv = f coVs -eos vt (II-81) and substitution of hf cosh2u - cos2v for dA in Equation (III-80) yields f _ ( P)" sd V (III-82) q is now obtained by integrating dq over the ellipse from v = 0 to v = 2ir. v-2rr Ae) -/ -A/ v (III-83) However, because of obvious flow symmetry about the x and y axes (see Figure 24) the right side of Equation (III-83) can be multiplied by 4 and the integration carried out from v = 0 to v = t/2. A XV (III-84) The cumulative volume of water which has flowed across the elliptic boundary up to time t, QtEi is given by A ~ J7 = a jjt (AII-85) 9tE~~~~gce~~t ~~

-92c. Numerical Method of Solving the Diffusivity Equation in. Elliptic Coordinates The first step in any numerical treatment of the diffusivity equation is the replacement of the derivatives in the equation by proper finite difference forms. The resulting difference equation is called'explicit' or'implicitt depending upon the time index attached to the difference forms of the spatial derivativesi (1718) For example, consider the equation P P(x, t) _ P(K~ t) ( III-86) The explicit difference form of (III-86) is I -2.c *PmP PL, L / -Pm, (111-87) where the spatial derivative 52p/8x2 is replaced by a difference form of time index i. Thus the unknown value of P at time t = (i+l)At can be determined explicitly in terms of known values of P at time t = iAt. An implicit difference representation of (III-86) is p-,,L. 2 PM& P, P,+, -Pv - ( I -88) Here the derivative 82P/8x2 has been replaced by a difference form of time index i+l and the unknown values of P at time t = (i+l)At are bound together by a system of simultaneous equations. Thus the expllcit form seems preferable to the

-93implicit in the respect that the former allows explicit calculation of the unknown P values while the latter form usually requires an iterative procedure to calculate P at time t = (i+l)At. However, the utility of the explicit form is severely restricted by the stability requirement(19) that 6 t ( ~2x) If this restriction is violated then small errors arising at the beginning of the calculations grow and eventually become so large that the results are meaningless. This condition is commonly referred to as'computational instability'. Thus if Ax is made small to reduce truncation error (i.e., error resulting from approximation of derivatives by finite difference forms) then At must be accordingly reduced, often to such small values that an impractical number of time steps must be employed to cover the desired time period. Bruce et al(l8) report that use of an explicit difference representation of the diffusivity equation governing unsteady state gas flow would entail 20,000 time steps to remove half the gas from the system they considered. The implicit form, (III-88), on the other hand, is stable for any size time increment.(l9) The disadvantage associated with the restriction on the size of At usually outweighs the disadvantages posed by the greater difficulty of solving the implicit difference equation. Thus the majority of the numerical treatments of the diffusivity equation reported in the literature employ the implicit form.

The numerical procedure presented here is a combination of the alternating-direction implicit method developed by Peaceman and Rachford(20) and a method given by Richtmyer.(21) Definition of the new variables P(ugv t) Po - P(u9vgt) tDE = kt/Lcpf2 = dimensionless time for elliptic flow simplifies Equation (III-74) to' + 7=Cus P —-- S-) lPts ( III-89) The alternating-direction method involves replacing, at odd numbered time steps, the derivative a2p/6u2 by a difference form composed of the unknown values of P at time t = (2i+l)At and the derivative 62p/6v2 by a difference form composed of known P values at time t = 2iAt. At even numbered time steps, 62-f/u2 is replaced by a difference term involving known values of P at time t = (2i+l)At and ~2P/~v2 is replaced by a difference term involving unknown P values at time t = (2i+2)At. Thus the difference representation of (III-89) at odd numbered time steps [i.e., t = (2i+l)At] is 2n2- t, +SI, 2 anP,}; tPSM2 n-I 2. h2 2 I. - +_ ( v)' (III-90) stir.)EP ti,1

-95and the representation of (III-89) at even numbered time steps [i.e., t = (2i+2)At] is ~t:~n/2i+, -2,,n. +.F7, Pl Lt, ~.nt t -'2 mniL, ~.,. (III-91) where Pmn,i - the value of P at u = mu, v = rv and t = iAt. Equations (III-90) and (III-91) can be rearranged to (III-92) and (III-93), respectively, in which all the unknown values of P are on the left side of the equations, P.aL t(fv,+2)- tm nr)2+/ =R,2C % m2nL; (III-92) P^,nits~(m+lf 2 R)3R( eL P)mj (III-93 ) where,_'-,, C, (c +sh,,; -, os. nv),~,- T n -,it;(frnati2v +r ( mert gie by.., ~ ~ L), TRhtmyer(21) i or svg pli cit mdif en e t The non-iternative or explicit method given by Richtmyer(21) for solving implicit difference equations systems involving one spatial variable is easily extended

-96to the case of Equations (III-92) and (III-93) involving the two spatial variables u and v. Consider, for the moment, Equation (III-92). Let us assume a relationship of the following kind between Pm,n,2i+l and Pm+ln,2i+l. FM naLE/ - Crnn Fm ) n L 2I + Sm)n (III-94) Cmn and Dmqn are supposed to be functions of m and n, i.e., functions of u and v in the u-v planar grid. Then -/no; L; h/ P)m) Y))InLt +* IPD )r (III-95) and substitution of this expression for Pm-l,n,2i+l into Equation (III-92) yields Pm+I, i+l -u +L m+I (III-96) L.+n; I fmsn+2 I- 2. InlComparison of (III-96) with (III-94) shows that Cm~= - +Mn -+2 -Crmrl.r and (III-97) Similarly, assuming that P- naL+; 4tnn Prr7Jntl lai +z Frnn (III-98) one finds and _ ~~~and R ~~ f~ n(III-99) F R IFm,- +3b nR w:+ X rIn i-2\ LR'i1RErJ -i

4 -97The recursive relationships (III-97) and (III-99) allow explicit calculation of all Cmn and Dmn at odd numbered time steps [Am,n,2i is known at odd numbered time steps when P at t = (2i+l)At is being determined] and Em n and Fmn at even numbered time steps (Bmn,,2i+l is known at these time steps) if starting values of these coefficients can be obtained. After the values of these coefficients are known, the pressures Pm,n,2i+l and Pm,n,2i+2 can be calculated explicitly from Equations (III-94) and (III-98), respectively. Attention is now given to the boundary conditions since these conditions must be employed to obtain the starting values of Cmns Dmgnp Em n~ and Fmn. Since the liquid flow is obviously symmetrical about the x and y axes, as seen from Figure 24, only the quadrant of the plane from v = 0 to v = </2 need be considered. The boundaries of the flow model in this quadrant are characterized by the following m and n index values: on'the gas field ellipse: u: Ub and m = Mo = ub/Au on the aquifer exterior boundary: u = ue and m = Me - (Ue - 1/2 auyAu on the positive x axis, x > a (see Figure 24): v = 0 and n = No = O/Av = 0 on the positive y axis, y > b (see Figure 24): v = </2 and n = <2 = Ne Av

-98Because of the flow symmetry about the x and y axes, the gradient of the pressure in the v direction is 0 at v = 0 and v = it/2. If the aquifer exterior boundary (u = ue) is assumed closed (i.e., no flow across exterior boundary), the initial pressure is assumed equal to Po, and if the pressure is dropped by unity at u = ub for t > 09 then the initial and boundary conditions become (9~ ap(Lj)_.. = o d P(uVt) ) = o C? 0., C eP(re)s in P(dieneo) When expressed in difference form these conditions are PM) DJ o = - -Iv; (III-100):/l) - e. (III-101) F14enoL= (III102) (111-102) P~,mno - 0 (III-103) 7S~,rxL I (III-104) Since PMe.ni = PMe+l nyi the effective position of the exterior boundary is at u = (Me + 1/2)hu; this is the reason that Me was chosen above as (ue -1/2BuI rather than as Ue/AU..

-99The starting values of Cmn and Dm~n are now obtained in the following manner. Letting m = Mo + 1 in Equation (III-92) and inserting the value PMon,2i+l = 1 [from Equation (III-104) ] one obtains eP PtcrXon/;tl i cF)N.I/ 4 1) e (III-105) Comparison of (III-105) with (III-94) yields (III-106) Equations (III-97), in conjunction with (111-106), now allow explicit calculation of Cmn and Dm, n for m = Mo + 2, Mo + 3..., Me, and for n = 01, 2,...9 Ne. At the exterior boundary, u = ue, we have from Equations (III-102) and (III-94), Pe =VI Ctjjn +2,n, L+j +S or VP)1e nihetI = 1'c (III-107) After the value of PMeln2i+l is calculated from (III-107), the values of Pm,n,2i+lD for m - Me - 1, Me - 2,..., Mo + 1 (and for n = 0, 1, 2,. s, Ne), can be explicitly calculated from Equation (III-94). The starting Emn and Fmnn values can be obtained by letting n = O in (III-93) and setting Pm,-1,2i+2 equal

-100to Pmo,2i+2 [from Equation (III-100)1] to obtain —, OR — 2 L- -gwhich, by comparison with (III-98), yields Ems ="Ae, A (III-108) The starting values Em,0 and Fmto and the recursive relationships (III-99) now allow calculation of Em,n and Fmpn for m = Mo + 1, Mo + 2,..., Me, and for n = 0, 1 2,..., Ne. At the upper boundary v = i/2, we have from Equations (III-101), (III-98), and (III-93), and Pm,Ner* =+ P*F,,ee o I.;* + z * 2...2 Eliminating PmNe+l,2i+2 from these two equations, one obtains F )me A E (III-109) After PmgNe,2i+2 is calculated from (III-109), the values of Pm,n,2i+2 for all m and n = N e - e 2..., 0, can be calculated from (III-98). At each time step the calculated P values can be employed in the determination of a quantity, QtDE' defined __D Jo- (rtbE)dtD (III-0llO) q*~E 0

O101where 9 PCUC v t/ aK PO a b Simpson's rule of numerical integration can be employed to yield eVt~ --- z( — +( 3 2 Vg~v= o + ( =V)_ (111-112) or, expressing the derivatives oaP/u by differences, since ( =o and Oi) grati f rom the boundary condition ( -III104). The values of q at successive time steps, determined by inserting the calculated values of PMo+ly ni into (III-113) can then be employed to determine atDE by a Simpson's integration as (111i-114) + Yf (3ftp) - - - f,' (J-)4tpE)'-7 The term QtDE or QjAtDE is the dimensionless production quantity, for elliptic flow, corresponding to the dimensionless production quantity QtDl for radial flow, tabulated by Van Everdingen and Hurst.(l)

-102d. Comparison Between Elliptic and Radial Flow The above described numerical method has been employed to calculate the water influx into an elliptic cylinder sink having an eccentricity of 0.925 or a ratio of 2.63 between the major and minor axes. The exterior boundary of the flow model was specified as a confocal ellipse encompassing an area equal to 101.3 times the area of the inner ellipse. The boundary conditions employed in the solution are given in Equations (III-100) - (II-o104). The variable u ranged from a value of 0.4 (corresponding to an eccentricity of 0.925) on the interior ellipse to 2.6 on the exterior boundary. Since the area of an ellipse is itab, where a - semi-major axis = f(cosh(u)) b = semi-minor axis = f(sinh(u)), the ratio between the areas of the exterior and interior ellipses is cosh(2,6)sinhk2.6 or 101.3 (The ratio of cosh(. sinh. the area included within the exterior boundary of a flow model to that included within the interior boundary is designated hereafter as the l area ratio'.) The variable v ranged from 0 to T/2, covering the first quadrant, (see Figure 24), and dimensionless times tDE, was varied from 0 to 426 in increments of the following values:

t.DE AtDE 0-1.01 1-9.05 9-26.1 26-106 5 106-426 1 A Au increment of.08 and a Av of.098175 were employed in the numerical solution. The calculations were programmed in the FORTRAN code and carried out by an IBM 704 digital computer. The calculated values of QtDE [see Equation (III-110)] are listed in Table VII and plotted in Figure 25 as a function of dimensionless time. The latter figure shows that. for all practical purposes, a steady state or asymptotic value of QtDE is attained at tDE = 200. Figure 25 also shows QtD' the dimensionless influx quantity tabulated in the literature(lp2) for the constant pressure radial flow case, plotted versus time. These QtD values were calculated(l) for a circular cylinder flow model having an exterior radius 10 times the interior radius; that is, the area ratio for the flow model is 100. The difference between the 100 and 101.3 area ratios for the radial and elliptic flow models is taken into account in the comparison developed below. The boundary conditions employed by Van Everdingen and Hurst(l) in obtaining the QtD values are identical to those used here in solving

-1go4TABLE VII QtDE AS A FUNCTION OF DIMENSIONLESS TIME tDE kt DE =Icpcf2 Av =.098175 AtDE = given on page 103 of Thesis Au = o08 tLQ~tDE * AQt\ tDE QtDE (00 x - ) tDE QtDE oo 0 0 90 32.58599!1.47584 4,4 95 32.92836.5 1.29198 1.9 100 33 22521 1 1. 98924,99 110 33 70573 2 35,16138 120 34.06693 3 4 4,18009.51 130 34 33840 4 5 11626.50 140 34,54244 6 6o83881,43 150 34.69-578.33 8 8 42793 16o 34,81102 10 9.91663.36 170 34.89761 15 13,27534 180 34,96267.33 20 16,18515.31 190 35 o1154 25 18,70798 200 35.04821 30 20,89547 220 35.09638 35 22.79213 240 35.12344 40 24.43662 30 260 35 13860.33 45 25,86248 280 35514701 50 27,09876.30 300 35,15162 55 28L17068 320 35515404 60 29,10009 340 35015523.33 65 29.90592 360 35.15574 70 30 60 6462 380 355 15587 75 31.21043 400 35.15587.33 80 31 073569.31 420 35015587.33 85 32.1 o9112 100 x = maximum % error in QtDE' calculated from Equation (III-124)

-105-* — * RADIAL FLOW, FROM VAN EVERDINGEN a HURSTI -A- ---- ELLIPTIC FLOW, FROM NUMERICAL SOLUTION OF DIFFUSSIVITY EQUATION t D: RADIAL FLOW DIMENSIONLESS TIME t kt t DE * cfT, ELLIPTIC FLOW DIMENSIONLESS TIME 55 QtD 50 45 40 DED ~3 5 - |-.. — -- 30 ~~~o ~ ~ ~ ~ ~ ~ ~ ~~ ~o t I * Ot * 20 40 80 120160 2 0.2. 4.6.8 1.0 1.2 1.4 1.6 1.8 2.0 to or tDE 0 40 80 120 160 0tD and QtDE vs Dimensi0nless Time. toFigure 25t DE QtD and QtDg v8 Dimensionless Time. Figure 25

-106the elliptic flow diffusivity equations namely an initial pressure drop of O0 a closed exterior boundary, and a pressure drop of 1 for all time at the interior boundary. In order to compare the radial flow with elliptic flow, one must establish some basis of comparison. This basis cannot be the distance between the exterior and interior flow model boundaries since the elliptic model has no single characteristic dimension analogous to the radius of the circular model. The comparison has therefore been made on a basis of equal areas encompassed by the exterior ellipse and the exterior circle and equal areas included within the interior ellipse and the interior circle. Thus each flow model contains the same volume (or mass) of water. The equality between the areas included within the interior boundaries yields the relationship rrrj 7 f'cost a b sin bu where rb = radius of inner boundary of circular cylinder flow model ub = value of u on inner ellipse of elliptic cylinder flow model. Since ub =.4, one obtains ErL_: f ~' o-A-osC oh (,,' - o.-r.f (III-115) Equation (1II-115) is employed in relating tDE to tD as to =, -.=,..4a,2S= f,"j 2t':,- (III-116)

-107Thus corresponding QtD and QtDE values should be taken at tD = 2.245 tDE rather than at tD = tDE. The actual cumulative water influx across the inner elliptic boundary is related to QtDE in the following manner. (tjE= t#C'1P( tp ( -III-117) QtE = total volume of water influx across inner elliptic boundary at time t, ft3 AP - pressure drop at inner elliptic boundary, psi h = thickness or height of elliptic cylinder flow model, ft. f foci of ellipse are at x = + f, y = 0 (see Figure 24) The actual volume of water influx into the inner circular sink of the radial flow model is related to QtD as(l) Qt=.a h cr j Ci;MPL2t; (III-118) Qt = total volume of water influx into circular sink of radius rb, at time t, ft3 h = thickness or height of circular cylinder flow model, ft. Thus a comparison between the actual water flow into the elliptic sink and that calculated by approximating the

ellipse as an equal area circle and employing the radial flow equation is afforded by the ratio _t - 1_ ctD 0,d O9 X (III-119) where tD = 2.245tDE" The relation Qt/QtE = 1 would denote exact duplication of the elliptic flow results by radial flow quantities calculated for an inner circular boundary including an area equal to that included by the inner elliptic boundary The ratio Qt/QtE is listed in Table VIII and plotted in Figure 26 as a function of dimensionless time tDE. This plot is actually a band rather than a line because of errors in the calculated QtDE values and in the QtD values. Error in the latter quantity arises from the fact that the literature does not contain QtD values corresponding to an area ratio of 101.3, Thus available QtD for a ratio of 100 were employed in calculating the points plotted as circles in Figure 26; a band of sufficient width to include the difference between the QtD for ratios of 100 and 101.3 was then drawn through the points. The determination of the various errors which are included within this band is discussed below, Figure 26 shows that application of radial flow calculations to the elliptic cylinder flow model described above results in. errors of the order of 7% in Qt for small dimensionless time tDE. The error decreases as dimensionless time increases and approaches 0 (i.e., Qt/QtE approaches

-109TABLE VIII Qt/QtE AS A FUNCTION OF DIMENSIONLESS TIME tDE Qt =.699 QtD QtE &QtDE tD = 2,245 tDE * _100 x ~*,; QtDE/QtRE * tDE QtDE tD ~7QtDQltQNeg. Pos. x 100 D __ Error Error QtD 1 1.9892 2,245 2.636.928 -.185 +.99 3 4.1800 6.75 5.598.936 -.185 +.51 6 6,8388 13o5 9.213.944 -.185 +.43 10 9.9166 22.45 13.416.948 -.185 +.36.04 20 16.1852 45 22,09.955 -.185 +.31.09 50 2700988 112,1 37.31.965 -.185 +.30.28 80 31 7357 180 44.21.975 -.185 +.31.46 110 3357057 247 47,12.98 -.185 + o32.62 140 34.5424 314.3 48,45.982 -.185 +.33.84 170 34.8976 382 49.03.984 -.185 +.33.99 200 3550482 450 49~30.985 -.185 +.33 1017 230 35.1119 516 49.5.986 -.185 +.33 1.31 400 35 1559 900 50 o5*.999 -.185 +.33 0 The listed QtD values were calculated(1) for an area ratio of 100; tD = (QtD)o) 101 3 (Qt) 100 ZQtDE was calculated from Equation (III-124)o ***This QtD value corresponds to an area ratio of 101.3,

Qt OtEm 1.0 1. 00 Qt Qto 0.99 0.''tE i~ t DE I I I 1ll 0.98 0. 97 - ~;~tD Corresponding to correct area ratio of 101.3 employed Y10a 0.96 to calculate this point. 0.95 0.94 ERROR BAND DUE TO ERROR IN OtDE 0 CALCULATED POINTS O. 92 0 20 40 60 80 100 120 140 160 180 200 220 240 260'400 t DE DIMENSiONLESS TIO ME tOIN t c re 26. C flow. Figure 26. Comparison Between Elliptic and Radiul Flow,

-1111.0) for large time. This approach of Qt/QtE to 1 as tD or tDE grows large is a good check on the accuracy of the calculated QtDE values, since for large time the pressure drop approaches a steady-state uniform value of 1 and the total expansion of the (equal) volumes of water in both flow models (i.e,0 the cumulative influx Qt or QtE into the inner sinks) should be identical. The circle plotted on Figure 26 at the far right (tDE = 400) is the only point in this figure calculated from a QtD value corresponding to an area ratio of 101.3 (the other points are calculated from QtD values corresponding to an area ratio of 100). The reason for this fact is that the large-time or asymptotic QtD value can be exactly determined for any given area ratio as(l) (rP/rb)- IThus (QtD) asymptotic equals 50.15 for an area ratio (re/rb)2 of 101.3 and employment of this value in the ratio Qt/QtE ='QtJQtDE yields the plotted (Figure 26) value of 0.999. From Figure 26 the conclusion can be drawn that the error incurred by applying radial flow calculations to elliptic gas storage reservoirs is inversely proportional to the magnitude of the tD values employed in the calculations (for tD > 1). For example, suppose an elliptically shaped reservoir (with a finite, impervious exterior boundary) is approximated by an equal area circle and

-112Equation (III-15a) is employed to calculate the water influx, Qt, for 20 time steps. Then the error in Qt will be larger if the range of tD values, for the time period concerned, is 0 - 20, with a AtD of 1, than if the range is 0 - 1000, with a AtD of 50. This fact is evidenced by the curve plotted in Figure 26, since the error in the calculated Qt value is larger for small tD than for large tD. The fact that the error in Qt ranges approximately from 7% to 2% for intermediate tD from 2.245 (tDE = 1) to 224.5 (tDE = 100) indicates that a need exists for analytical solution of the elliptic flow diffusivity equation for larger (i.e., more practical) aquifer/reservoir volume ratios. e. Error Analysis Error in the ratio Qt/QtE =.699 QtYQtDE is caused by errors in the QtD values employed and in the QtDE values calculated as described above. Thus if AQtDE denotes the error in QtDE and AQtD the error in QtD, the error in the ratio Qt/QtE will be Thus the percentage error in Qt/QtE is (st/4a, /00 = Dt )/00 (III-120) or the difference between the percentage errors in QtD and QtDE' The error components AQtDitD and AQtDEE/ tDE have

-113been calculated, as described below, and are listed separately in Table VIII. Errors in the calculated DE values arise from truncation and computer round-off error. Truncation error is introduced by approximation of the differential Equation (III-74) by the difference Equations (III-90) and (III-91). Round-off error is caused by the fact that the computer carries only a finite number of digits. An upper limit on the round-off error is estimated in Appendix III as.00185 QtDE. A truncation error analysis based on a method given by Scarborough(22) is also described in Appendix III. This analysis yields the relation tE Ci ~6 -~ 2 +14) Bat f Catpe (III-121) where QtDE QtDE denotes the truncation error in QtDE and A, B, and C are numbers which may vary with tDE. If (QtE* - QtDE) is denoted by (AatDE)Tr and if iA29 BAv, and CAtDE are denoted by (AQtDE)Aup (EQtDE)Av9 and (EQ4tDE)ZtDE, respectively, then Equation (III-121) becomes ( t(Rt)Trun kn= ( gDE p q (III-122) The value of (SQtDE)Au can be determined by calculating QtDE for two different sets of increment values, each set having identical Av increments and identical AtDE increments but different Au increments. For example, if Qi (i = 1,2) denotes the value of QtDE calculated with a set

-114of increments (Aui, AV, AtDE), then the value of A can be obtained from (III-121) as The value of (AQtDE )u; the truncation error in QtDE due to Au, is then 1 -. B. - z,= This method of determining the error terms on the right side of Equation (III-122) has been employed to obtain the results listed in Table IX. These results show that, for tDE > 1, the error in QtDE due to Av and AtDE is small compared to that due to Au. The combined errors due to Av and AtDE are ignored here since they cause an error in Qt/QtE of less absolute value than the radius of the circles used to plot Qt/QtE in Figure 26. Thus Equation (III-122) is rewritten (8&t3e) -r (f = 14(au) (III-123) where A(Au)2 is always positive (see Table IX) and therefore causes only a negative error in Qt/QtE, as shown by Equation (III-120) o The total error in QtDE is equal to the sum of the round-off and truncation errors. Thus if QtDE denotes the true solution and QtDE is. as before, the machine solution, then iI4 _7 ~DAd /87492 (III-124)

TABLE IX INDIVIDUAL EFFECTS OF Au, Av, AND AtDE ON THE TRUNCATION ERROR IN QtDE Effect of Au Av = o098175 AtDE = given on page 103 QtDE tDE Au =.08 Au =.133333 (QtDE)u ( ) 1 1.98924 1.96054.0162 5 5.99779 5.96793.0168 9 9.18374 9.15236 o0176 26 19.17078 19.12799.0240 50 27.09876 27.03897.0336 10o6 33052980 33044528.0474 276 35 14570 35505145.0529 426 35.15587 35006323.0520 Effect of Av Au =..133333 AtDE = given on page 103 QtDE tDE Av =.0715 v = 098175 (ADE) = B(.098175) (Z —QtDE)Av = B 1 1. 96063.0 96054.000337 5 5.96797 5o96793 o000147 9 9,15240 94!5236 o 000147 26 19.12804 19.12799 o 000184 50 27.03904 27.03897.000258 106 33.44539 33.44528.000405 276 35.05122 35005145 -.000846 426 35.06275 35506323 -.00177

-116TABLE IX (CONT'D) Effect of AtDE Au = 0.1 Av =.0o98175 AtDE given AtDE = twice tDE on page 103 values on page 103 ('tDE)AtDE = C (AtDE) page 103 1 1,97852 1.97859 -~00007 5 5 98696 5.98684 o00012 9 9.18149 9.18128.00021 26 19.48793 19,48742.00051 50 28,24847 28.24800 o00047 o106 3625795 36.2567.00125 276 38,82800 38,82292,00508 426 38.86234 38~85430.oo804

-117where Q~tDE is QtDE - QtDE and +.00185 QtDE is the maximum round-off errors as mentioned above and estimated in Appendix III. nQtDE has been calculated from (I11-124) for a range of tDE values and is listed in Table VII as a percentage of QtDE [the maximum value of the right side of (III-124) was employed]. Attention is now given to the error caused in Qt/QtE by use of erroneous QtD values. The term Qt/QtE is claimed here to be a comparison between elliptic and radial flow for equal areas included by the inner boundaries and for equal areas included by the outer boundaries of the respective flow models. The equality of the areas enclosed by the inner ellipse and the inner circle is assured by the relationship (III-115). However, the QtDE values were calculated for an external boundary area 101.3 times the inner area while QtD was obtained from the literature(l) for an external circular area 100 times the inner area. If (QtD)101.3 and (QtD)100 denote the values calculated for the subscripted area ratios, then (;tp),0 3 - t)lOI 9at (III-125) where AQtD is the error of interest here. The maximum error, (~atD)max ~ as well as the maximum DAtr(/tD)100 occurs in the asymptotic value of QtD. This fact is apparent from tables of QtD versus tD given in the literature(l) for various area ratios9 as well as from the observation that for

-118sufficiently small tD, QtD will be independent of the potition of the exterior boundary. This maximum AQtD can be exactly determined since the asymptotic value of QtD is the following explicit function(l) of the area ratio: (;)csy/mpfottic 2 The term (re/rb)2 is the area ratio, i.e., the ratio of the area encompassed by the exterior circle to that encompassed by the inner circular boundary. Thus, from (III-125), /0/3 -/ /00 - / =(a = 0.6 and _ &Zt3j 1 - _ /3 / This tD (Q-tD)100 ratio is 0 for dimensionless time tD = 15(1) and increases to 0.0131 at tD = 500, where the asymptotic (QtD)100 value of 49.5 is reached. The manner in which'AQt&ttD [= AQtJ(QtD)loO] varies from 0 to.0131 has been estimated from tables of QtD vs tD, for various area ratios, given in the literature.(l) The AQtD/vtD values$ estimated as described immediately above, and AQtDE/QtDE, calculated from Equation (III-124), have been employed in (III-120) to yield the maximum percentage error in Qt/QtE. This calculated error is plotted as a band about Qt/QtE in Figure 26. f. Stability Analysis The numerical method of calculation employing Equations (III-92) and (III-93) proved to be stable and

-1.19convergent for all combinations of Au, Av, and AtDE values used. The question remains, however, as to whether the author was merely fortunate in avoiding combinations which would result in instability or whether the method is actually stable for any combination of increment values. The method is shown by the following analysis to be stable for all increment values. If the solution for Pm n i is assumed to be of the form Pn,*7, L( rDmn,, e (II 126) then the condition for stability of the numerical method of solution described above is that IIDm,n,2i+2/Dmn,n2i be < 1 for all m, n. and i.(21) This method of stability analysis is strictly valid only for difference or differential equations with constant coefficients. However, as pointed out by Mickley et al.(19) the method can be applied to variable coefficient equations by treating the variable coefficients as fixed constants in the determination of the form of IDm,n,2i+2/Dm,n,2il and then letting them take on their most adverse values when determinining whether this latter ratio is < 1 for all m. n, and i. Substitution of Pmn,i from (III-126) into (III-90) yields Pm. Wn)_ — n. / =e (II 127)

-120Similarly, substituting Pm)n i from (III-126) into (III-91), one obtains Jmtn,2Lt fm, n 2 / caS j 28 =.0 ) - (III-128) Dm)n2L + C m n 2 R (/ - cos a v) One can readily obtain from Equations (III-127) and (III-128) the desired ratio Sm. n) 2 if v fm n -2 R (1 - to) frm/ r, -2(-cor sj1U ) D,,,, 2 -. + 2!<-CSjh/, In (III-129) PM^n} 2L arm n + X(/-gosiaC) Am, n 2 (I- aSlz Now, by inspection, it can be seen that the absolute value of the right side of (III-129) is < 1 for all m, n, j,,y fmn, and R provided fm n and R are both > O. Since fm n is [(au)2/atDE] [cosh2(n~A) - cos2(nAv)] and R is (Au/Av)2, both are greater than 0 and Dm,n,2i+2/IDmn,2i is less than 1 in absolute value for all Au,, Av, and AtDE. Thus the numerical method developed above, employing Equations (III-90) and (III-91), is stable for any combination of Au, Av, and AtDE increment values.

IV. LDI7 VIDUAL PROBLEMS ENCOUNTERED IN THE PREDICTITON OF GAS STORAGE RESERVOIR BEHATVIOR The equations of Sec,,fon 11-A would allow a completely general method of predicting gas storage reservoir performance if all reservoirs satisfied the assumptions on pages 8 and 9 of this thesis. Several cases arise, however, in which certain of these assumptions are completely invalid. For example, the assumption of a constant gas bubble radius is certainly not satisfied during the initial period of growth of an aquifer storage reservoir since the radius increases from a value of O at the time of initiation of the gas bubble to values of several thousand feet. Also, the solutions to the diffusivity equation presented in Section III-A are valid neither for cases in which the initial aquifer pressure dIstribution is a non-uniform function of radius nor for cases involving two or more gas fields situated on a common aquifer. Methods are developed below to treat the cases of gas field interference, nonuniform aquifer initial pressure distribution, and aquifer storage; the methods are then evaluated by application to actual gas storage reservoirs. Ao Treatment of Pressure Interference Between Two or More Gas Fields Situated on a Common Aquifer The absence of pressure interference effects from neighboring reservoirs was one of the assumptions empl.oyed in developing the Section 11I1-A equations governing gas storage reservoir behavior. In cases where +his assumption is not valid, a modification is required in the previously described and ill2ustrated (Section II1MB) method of application of these equ-ations to the prediction of gas field performance. Mortada(6) described such a m odif cation which allows prediction of the performance of several oil fields situated or. a common aquifer but he did not apply his method in an actual -121

-122case study. His interference theory is discussed below and an adaptation of his method, applicable to interfering gas fields, is then presented. This adapted method is finally applied in the performance prediction of two interfering fields and a comparison is made between their predicted and actual performances. The validity of both the adapted method presented below and Mortadats method depends upon satisfaction of all the assumptions listed in Section III-A excepting the one specifying absence of interference effects. 1. Existing Theory Relating to Interference Between Oil Reservoirs M. Mortada(6) presented a method of aquifer water influx calculation which accounts for interference between two or more oil fields situated on a common aquifer. His calculations involve primarily Equation (IV-i), below, which is obtained by application of Duhamel's superposition principle to Van Everdingen and Hurst's "constant terminal rate case" solution to the previously presented diffusivity Equation (III-5). Equation (IV-l) is identical with Equation (III-16) except for the retention in (IV-1) of the independent variable rD in the aquifer pressure P(rD t) and in the dimensionless pressure drop term P(rD,tD). I, g-P(5Xt) 27r~k~o ~+I ) J (IV-1) PO = initial uniform aquifer pressure, psia P(DrD t) = aquifer pressure at rDxrb feet from center of oil field and at time t, psia rb = oil field radius, feet P(rD,(ji)AtD)= dimensionless pressure drop term, obtained as solution to Equation (III-5), at dimensionless radius rD and dimensionless time (j-i)YtD

-123t t= ime = jAt, seconds. Mortada equates the pressure drop in each oil field to the sum of a number of individual pressure drop components, where each component represents the effect of one of the interfering oil fields on the pressure drop in the field considered. Thus, in the case of two interfering oil fields denoted by field R and field T, PtJ - PRc PPR ) P' A k (IV-2) P P o -to - (IV-3) Pr(jt) r (P),at)= Pr7 -t) 7t4 (tIV3) PR( jt) = total pressure drop in field R at time t = jZt, psia PT(jat) = total pressure drop in field T at time t - jAt, psia PO = initial pressure in reservoir-aquifer system, psia PRR(jZNt) = pressure drop in field R caused by R's production, psia PRT(jnt) = pressure drop in field R caused by T's production, psia PTT(jt) = pressure drop in field T caused by T's production, psia PTR(jZt) = pressure drop in field T caused by R's production, psia

-124If n interfering oil fields are involved, then the pressure drop in field i is 4? (J a it) > Ha F gjD t) (IV-4) rn- I where Pim(jAt) is the pressure drop in field i caused by field m's production. The terms PR (jAm) and PT(jlt) are calculated as shown in Equations (IV-5) and (IV-6), below. (Iv-5) X - Q( P rLn (J - pZa tj.)l klt where AtD = ( 2 and (rb)R is the field R radius in feet AtDT krt 2 and (rb)T is the field T radius in feet tDT = 2 b Ts qRi = the average rate of water influx into field R during the time increment from (i-l)At to itt, ft 3/second i = the average rate of water influx into field T during the time increment from (i-l)At to iAt, fto3/second P(1 ji)tD) = PtD, the dimensionless pressure drop quantity tabulated as a function of dimensionless time tD by Van Everdingen and Hurst(l) and Chatas(2)0 r - distance between centers of fields R and T divided DRT by field T radius, (rb)T, dimensionless P(rDT (j-i)ttD) = dimensionless pressure drop quantity presented graphically by Mortada(6) as a function of rD and tD

-125The first s': nation of the right side of (IV-5) is PRR(jAt) and the seod suratn s RT(J) Te expression for PT(jAt) is simmilarly, P (jat) = Lf ( r, - %re) v((( Ladt) &- lt -- ( Pa i( Ge6; (1,,:v-6) tR 2. Adaptation of Existing Oil Reservoir Interference Theory to the Treatment of Gas Field Interference Mortada assumes that the rates, qR and qT, of aquifer water influx can be equated to t;he rates of oil production in the corresponding fields or to these latter rates corrected for the expansion of the oil field content. In the case of gas fields, the relationship between the rate of aquifer water influx and the gas production or injection rate is, of course, not a direct proportionality and further attention must therefore be given to the terms qR and qTo (Fields R and T are now assumed to be interfering gas fieL ds) The aquifer water influx rate qi in Equation (V-i) can be related to the gas field pore volume Vi as outlined in Section III, page 21. of this thesis, to obtain. Prl t)+- =3?P t) Z trP v1 (IV-7) where P(rDt) aquifer pressure, at radius rDxrb feet from center of gas field, at time t - jAt, psia rt = gas:field radius, feet i = gas field pore volume at time t = itt, ft)3 K3 L psi/fto3 At - -t2me 4increment, seconds

-126The total pressure drop in gas field R is now determined by first calculating the pressure change in R which would occur if no other interfering gas field T were present (i.e., ignoring the presence of field T entirely) and then adding to that calculated pressure drop the pressure change in R caused by T's gas production (or injection). Thus the pressure drop in field R is p~ + = (P).)R.' (JAf) = 7zt) + %ut) (Iv-8) and PRR(jAt) is first calculated, ignoring field T, as -< 2 I, - P00(vt) = z3 (A - L-/'V1i)-PC/>JE L at, (IV-9) i= " Similarly, PTT(jAt) is calculated, ignoring field R, as PTT(jdt) =RJ( (~V;;-Vri-,-YT:.~r)P - -VI(/gJ- (IV-10) where V(or i) = uncorrected' gas field R (or T) pore volume at time t = iLt, ft 3. PR(jat) (or PT(jat))= gas field R (or T) pressure at time t = jAt, psia (Po)R (or (Po)T) = initial gas field R (or T) pressure, psia. The terms VTi and V are'uncorrected' pore volumes in the respect that they are calculated from the Equations (III-20a and III-20b), below, ignoring the interference effects. VKr< = n[iRfrf/ 4+ ( k (III-20a) Vr = nTr; +/ R gr( 7 (III-20bo) ~r~ ~ ~.A t)'l-~~..

-127where i (or Ti) pound oles gasin place in field R (or T) at time t =~t, PRR'it) (or PTT(,Pt)) - (Po)R - PR(iAt)[or (Po)T - PTT(iAt) Elmi natior of Vof between Equations (IV-9) and (IIl-20a) and of V'. between Equations (IV-10) and (ITI-20b) yields the following equations PR 4 t) 2 (IV-.) FTrf(ja4t) - E (IV-12) where L 1 1 - - K8 and K9T are same as K8R and K9R except all'R' subscripts are changed to ITt. Equations (IV-11) and (IV-l2) are explicit expressions a'llow-n.g calculation of PRR(jt) and PTT(jAt) if the cal. culations are begun. with j = I and continued w-+.h j 2,3.... etc., since then for a given value of j, all the terms in K8 and K9 are knowno After all the PRR(jAt) arnd PTT(,At) are calculated (j - 1,2,3, 2 N whiere N is the to+tal number of time inhrements for which the predicted gas field performance is required) the terms VI and VT (i.1 2 3.., N) are all known (as calculated from Equations (Iii-20a) and cfield R due to Ts production ) and ) and PRjt) (pressure change in. field T due to RtS production) are then calculated from. tshe equations w~pre f th tofain~mbrQftf~en@rmn~~sfo wzi-, RT p~eije

-128-'Pr j t) /t j'( V( - Kl - (JIL-),)) ( V-13) Proq ( j t)= k3 (2VW,-VRL,-VRL't) P(rfrR) - i) tDR (IV-14) = C where, letting d = the distance in feet between the centers of gas fields R and T, rDRT = d/(b ) rD = /(rb)R (rb)R(or (rb)T) = gas field R (or T) radius, ft. AtDR (or AtDT) = dimensionless time increment for field R (or T). The total pressure drops in fields R and T are then given by (lbzt) = (/jR,- (Jt) =- — jIt) (from Equation IV-11) (iv-8) +-PR-T(jt) (from Equation IV-13) — OJ = (?P)T ~ - S tPr ) = P-r(J'Dt) ( from Equation IV-12) + fr (j&t) (from Equation IV-14) The correct or predicted gas field pore volumes are finally calculated from the relations V (ji t) nj rT (/< + ) (IV-16) VTn j r) = e' l T (eok + ee) )(IV-17) In order to employ the equations presented here, one must know or specify the values of d, (Po)R, (Po)T, (nRi, nTi, i = 1,2,3,.N~, N), K3, at, KR2T = K K (where....,rN)3 K2T ~(~)'~H Wpc(rb)

-l29z L- + KP - gas eornpresso l 1t, ty factor), and RT (R = 0,3 ft -3 psia/# mole - R, -= gas:field formation temperature, ~R). 53~. Application of the Presented Gas Field Interference Theory to an Actual Case Study Michigan gas fields C and D were discovered in 1944 with a distance of approximately four miles separating their centers, These fields were originally considered as'constant volume' fields; as the fields were produced, however, the field pore volumes, calcu-.ated from the data, decreased significantly. This pore volume reducrtion, the production of water from wells near the areal boundary of the gas sand, and the approx4mately equal values of initial pressure and depth below the surface for the two fields all point toward the existence of a common aquifer underlying both reservoirs. The presence of interference effects between the two fields is evidenced by the failure of many attempts to reproduce the actual vol ume-time and pressure-time behavior of the fields by the equations and method of call.culation given. in Sections III-A and III-Bo The 12 year period from. 1944 to 1956 was chosen for the purpose of application of the gas field interference theory presented above to the gas reservoirs C and D. A time increment of one mon.th and Equat.ions (IV!1-1) through (IV-17) were employed in the cal culation of the predicted volume and pressure behavior of each of th.e two fie.lds, abe predicted and actual field performance are plo-tted inr Figures 27 and 28 and the values of the reservoir constants employed in. the calculations are noted on these figures~

1.175 I.150 1.125 - --- WL.UES OBTAINED FROM DATA 0 0 0 PREDICTED VALUES USING INTERFERENCE METHOD 1.100 A A A PREDICTED VALUES WITHOUT USING INTERFERENCE METHOD 00 0 I.075 - I.050 - / -.E'- V. A AAAAA 1.025 - FIELD C 0 0.A A AA..AA AA AA-AAAA A^ 1.000 0 Ao-AAA AA AaAA 0.975 a- - -.~00A00 00 A — AAA A u f ~ 0 AAAo A A AA:A A 0.950 _AAAAAA AAAAAA4A aa 0.925 Vo)C- 4..26 x 10e t ( Vo)D- 1.43 x 108o i w 4 4 kJ (Ki)c- -h c(rb)- 4.56x I0 ( Ki)D- h c (rb[D- 0.647 x 10: 1.000 h 6 Oft. > O165.. 0.4 65 t2/S/c.-psi ( 0.975 4 - At - I MONTH 2.628 x 10 Sec.. w 0.950 0.925 - 0.900A 000 0.9500, A ~~0.73~75 AA~~ 00000FIELD00 AA "AAAt 2 i t i,,% t,,A 0 5 AA 1 z50 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 105 lO o 10 1 120 125 130 135 140 145 0.7523 0 0.700:oo -4o NUMBER OF ONE MONTH TIME INCREMENTS FIELD C AND D PREDICTED AND ACTUAL PORE VOLUME RATIOS VS TIME Figure 27

1.7 1.6 -- WVALUES OBTAINED FROM DATA 1.5 o o o o PREDICTED VALUES USING INTERFERENCE METHOD 1.4 1.3!.2 1.1 Z 1.0 FIELD C 0.9 0.8 *0.o 0.7 00 ( P)C= 425.7 psio WELLHEAD 00.6 (Po)=D 431.7 psia WELLHEAD ~~~~~~~~w ~~(KI)c=. h c(rb) - 4.56 x 104 0.5 (Ki t h # c (rb) 0.6470 -10:w ~h ~60t. ft ~~~~~~O~~~~~~~~~~~~~~~~~~~ 0.4 x I0 ft/Sec-psi;f 1.0 ~~~~o L; d~~~~~~~At = 1 MONTH 2.628 x 106 Sec. 0.9 d 0.8 0.7 00 0.6 000000000000000000000000000 0 0.5 0.4 4 4 4 4 0.3 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 NUMBER OF ONE MONTH TIME INCREMENTS FIELD C AND D PREDICTED AND ACTUAL PRESSURE RATIOS VS TIME Figure 28

Table X lists tnr.e:field data ard pred-icted faield pressure a;,-., vo.l umen ratios for the period 1944-1956. Th. e cail.'.atorsas errployirg Equa+.olns (IV-Il) through (rV-17) were quite I ergtlhy arnd were theref-ore programmed in the FORTRAN code and perfocrmed by t-,he IBM 704 diga, -...om.pu.....er. An interesting observat.ion s ir.terj'ectd..e.re co:r-.cernit..g the necessity for taking the watler moveme:..t, i. to aco'unt. a.: a s, one might' question whether the act-u.al. fie.ld perfo:r.arce mig.rc+, bhe predicted to sufficient accuracy by smly' ass-.mig th-e:fed t;o he of con.stant vollme, An excel.lent; comparison betoween. t.he'co:rjns'ta..t., volumne and t;he twat;er infiux' metthods of predi-t+ion is giver. I.n Figure 27, since the'constan.t Vollme I perforima.r-ce is represe,.ted in this figure by the s't*raight. line VT/o' C. Sulppose or..e assulied. field C to have a constant vo'lume equal to * its original (discovervy) value of 4.2 x 108 ft 3 Tehen the error i ncurred in <:.h-e predic-.+ed pressure at 136 months after dis-ove>ry would. be approxim.ae>;.;.a- 6.o2% since an error of l6o22 would be incurred in the predicted volu'Lrre (see Figure 27) and pressure is approxirnatrel rv erver iv prcpornct - or.a —. to volumne (PTV' z~nRT) o However, Figur.re 27 shows t:ha th.>.e differe r.L.:e between the pore volume predic.ted tby the wa+;er inrifl1ux: (interf.ereenr: Ce) method and the act,tual poPre -vo um.e a+ this time is ornl c,/ 2n'6 a:rd th.-s s the error in t;he predicted pressure 136 m._onth.,,s after d.iscovrery would be approxim.,ately 2o6%. Ar error of approxi r..ate7.y 8 7% nr. t, pre dicted pressure at t = 136 mont;hs world be incurred by ass.rumin. a constant volume only over a period of 5 tmont'hs since Fl —gur. e 27 shows that the actual pore volumn charnges tby 807% f.rom t- 131 to 13(

-133Table X SUMMARY OF ACTUAL AND PREDICTED PERFORMANCE OF GAS FIELDS C AND D FOR PERIOD 1944-1956 Initial Field C Pressure = 425.7 psia, static wellhead K3 = ~/2~hk~t = 2.52 x 10-4 psia/ft3 Initial Field D Pressure = 431.7 psia~ static wellhead (Vo)C = 4.26 x 108 ft3 = IK2L k/~olrb~ (Vo)D 1.43 x 108 ft3 ~ = 1.65 x 10-8 i/seconds (K23D = k/~a(rb)~ = 1.163 x 10-7 I/seconds mt = 1 month = 2.628 x 106 seconds Cumulative SCF Wellhead Static Pressure Ratio Pore Volume Ratio Gas Production* P/Po* (End of Month) V/Vo* (End of Month) From Data Field C Field D Field C 1944 Initial 0 1 1 1 1 Dec. 114,7874.995d.992p.983P.997d.9999e.9999P.9997e.9997P 1945 Jan. 277j4254.988d.982p.966p.992d.9996e.9996p.999e.999P Feb. 410,8874.974d.973P.947p.998d.999e.999P.998e.998P Mar. 572,5614.963d.963p -.931p.999d.999e.9991o.997e.997P April 727,9884.941d.953P.917p 1.012d.998e.998p.996e.996p May 884,9484.941d.943p -.900p.999d.998e.998p -.994e.994p June 1,026,1644.937d.933P.900d.887p.993d.997e.997P.944d.993e.993P July 1,187,7894.923d.923p.883p.997d.996e.996p.991e.991p Aug. 1,342,6984.909d.913p.886d.880p 1.001d.996e.996p.983d.99Oe.990P Sept. 1,470,0~14.904d.905p.877P.996d.995e.995P.988e.989P Oct. 1,580,3004.897d.898p.873d.874p.996d.994e.994p.988d.987e.987p Nov. 1,680,2704.892d.892p -.867p.993d.993e.993P -.986e.986p Dec. 1,786,2664.880d.885p -.85910.999d.992e.992P -.984e.984p 1946 Jan. 1,903,4164.871d.878p -.845p 1.000d.991e.991p -.983e.983p Feb. 2,019,3484.869d.871p -.834p 993d.990e o991p -.981e.981p Mar. 2,137,7574.864d.863p -.825p.989d.990e.99Op -.979e.979P April 2,261,0814.852d.855P -.815p.993d.989e.989p -.977e.978p May 2,390,1164.85Od.846p.805p.984d.988e.989p -.975e.976p June 2,509,9074.838d.839P.803d.799P.989d.986e.988p.968d.973e.974p July 2,645,3874.826d.830P -.792p.992d.985e.987p -.971e.972p Aug. 2,795,3354.817d.820p.799d.786p.990d.984e.986p.954d.969e.970p Sept. 2,915,1904.810d.812p -.783p.998d.983e.986p -.967e.969p Oct. 3,027,3984.805d.805p -.775P.984d.982e.985p -.965e.967p Nov. 3,141,8364.791d.797P -.767p.992d.981e.984p -.963e.965p Dec. 3,258,8824.775d.789p -.756p 1.003d.979e o983p -.961e.963p 1947 Jan. 3,599, 0614.770d.780p.740p.996d.978e.982p -.958e 961p Feb. 3,533,6384.765d.771p.724p.989d.977e o981p -.956e 959P Mar. 3,676,6954.753d.761p.710p.992d.975e.980p ~.953e 957P April 3,796,2884.742d -753P.697P.996d.974e.980p -.950e 954p May 3,931,7924.746d.744p.688p.976d.973e.979P -.947e 952p June 4,072,1484.734d.735P.678p.980d.971e.978p -.944e 950p July 4,204,8274.725d.725p -.670p.978d.97Oe.977P.941e.948p Aug. 4,347,9904.709d.716p.664p.987d.968e.977P -.938e.945p Sept. 4,507,5984.706d.704p.660d.659P.973d.966e.976p.941d.935e.943p Oct. 4,753,0884 -.687p -.654p -.965e.975P -.933e.941p Nov. 4,984,1914 -.671p -.645p -.963e o974p.930e.9391o Dec. 5~219,8744.654p.636p -.961e.973P -.927e.937P 1948 Jan. 5,445, 3424 - ~ 638p.628p - ~ 959e o ~ 972P.924e ~ 934p Feb. 5,601, 9934 -.627p.621p -.957e.971p -.921e.932p Mar. 5,672,8754 -.622p.615p -.955e.970p -.917e.930p April 752,002b -.678p.609p -.954e.969p -.914e.928p May 1,066,749b.709d.701p.603p.968d.953e.969p -.911e.925P June 1,495,225b.731d.731p.596p.969d.953e.969p -.908e.923p July 2,227,027b.783d.783p.592p.969d.953e.969p -.905e.921p Aug. 2,742,427b.827d.818p.589p.959d.954e.970p.9(y2e.919p Sept. 3,045,194b.844d.838p i86d 585p 965d.954e 972p.915d.917p..... 899e Oct. 3,043,155b.847d.837P -.581p.961d.955e.973P -.896e.915p Nov. 2,774,029b.831d.817p -.578p.957d.955e.974p -.893e.913p Dec. 2,495,443b.808d.797P -.574p.962d.954e.975P -.890e.911p 1949 Jan. 2,171,643b.784d.774p.570P.962d.953e.976p -.887e.909p Feb. 1,857,267b.763d.752P.566p.960d.952e.976p -.884e.907p Mar. 1,676,078b.747d.739P - 564p 965d.951e.977P -.881e.905p April 1,721,350b.746d.742p - 561p 971d.950e.976p -.878e.904p MaY 2,336,595b.783d.786p - 557P 979d.949e.976p -.875e.90~p June 2,957,088b.825d.829p - 554p 982d.950e.976p -.872e.899P July 3~646~471b.879d.876p - 552p 973d.951e.977P -.869e.897P Aug. 4,324,981b.930d.921p - 550p 968d.952e.978p -.866e.895P Sept. 4,846,694b.967d.955P.542d.548p 966d.954e.98Op.904d.863e.893P Oct. 4,429,520b.920d.924p -.546p.988d.955e.983p -.860e.891p Nov. 4,397,883b.931d.921p -..544p.972d.955e.984p -.857e.888p Dec. 3,345,222b.863d.847p.542p.968d.955e.987p -.854e.886p 1950 Jan. 2,417,884b.793d.782p.539P.973d.953e.98910 -.851e.885p Feb. 1,952,432b.753d.749P.537P.984d.951e.99Op.848e.883p Mar. 1,294,869b.706d.7(Y2p.535P.985d.949e.991p.845e.880p April 1,400,996b.70~d.710p.533P 1.00Qd.947e.990p -.842e.878p May 2,937,336b.835d.820p.532p.967d.947e.986p -.839e.876p June 5,481,165b 1.009d.994p.531p.968d.949e.983p.836e.873P July 8,008,379b 1.166d 1.158p.530P.977d.954e.985p -.834e.870p Aug. 11,714,392b 1.391d 1.385p.530p.983d.962e.989p -.831e o867p Sept. 1R,867,813b 1.469d 1.444p.53Op.977d.969e.996p -.828e o864p Oct. 12,780,643b 1.441d 1.431p.53Op.995d.975e 1.003p -.825e.861p Nov 10,629,085b 1.313d 1.295P.530P.994d.978e 1.009p -.822e.858p Dec. 7,754,323b 1.110d 1.113p.530p 1.017d.977e 1.014p -.820e.856p 1951 Jan. 6,061,630b 1.006d 1.003P.529p 1.01Rd.975e 1.016p -.817e.853P Feb. 4,099,974b.874d.873P.529p 1.017d.972e 1.019P -.814e.851p

-134Table X ~oNT'o~ SUMMARY OF ACTUAL AND PREDICTED PERFORMANCE OF GAS FIELDS C AND D FOR PERIOD 1944-1956 ~umulative SCF Wellhead Static Pressure Ratio Pore Volume Ratio Gas Production* V/Vo* (End of Month) P/Po* (End of Month) Month From Data Field C Field D Field ~ Field D Mar. 2,857,796b.786d.789P - o528p 1.024d.968e 1.021p -.811e.849p April 2,915,068b.79Od.794p -.528p 1.025d.965e 1.019p -.809e.846p May 6,328, 30lb 1.033d 1.026p -.528p 1.00~d.966e 1.009P -.806e.844p June 9, 379, 874b 1.236d 1.219P -.529p ~ 993d ~ 971e 1. 009P -.803e.84Op July 11,103,548b 1.335d 1. 319P - ~ 530p ~ 999d.976e 1. 013p -.801e ~ 837P Aug. 1}, 236, 984b 1. 457d 1. 440p -.532p 1. 004d.983e 1. Ollp - ~ 798e ~ 833p Sept. 14,468,642b 1. 522d 1. 504p - ~ 532p 1. 009d ~ 989e 1. 023p - ~ 795e.829p Oct. 14,456, 835b 1. 506d 1. 495P - ~ 533P 1. 021d.994e 1.030p - ~ 793e.825P Nov. 12,930, 613b 1. 403d 1. 401p - ~ 534p 1. 034d ~ 997e 1.055P - ~ 79Oe.82 lp De c. 11,612,164b 1.326d 1.32 Op - ~ 535P 1. 034d ~ 999e 1.040p - ~ 788e.819p 1952 Jan. 9,306,955b 1.175d 1.178p ~ 535d ~ 537P 1. 048d.998e i. 044p.819d ~ 785e.816p Feb. 8, 384,215b 1.117d 1.121p ~ 519d ~ 538p 1. 050d ~ 997e 1. 046p.847d ~ 783e.815p Mar. 8,249, 887b 1.109d 1.112p ~ 539d.54010 i. 050d ~ 997e 1. 046p.814d.780e.8!3p Apri 1 8,732,541b 1.141d 1.142p.542d.541p I. 046d ~ 997e 1.046p.810d ~ 778e.811p May 10, 788, 594b 1. 272d 1.266p ~ 547d.544p I. 039d ~ 999e 1.045p.802d ~ 775e.807p June 12,962,585b 1.40~d 1.39OP ~ 557d ~ 551p 1. 037d 1. 003e 1. 046p ~ 792d ~ 773e.802p July 15,331, 883b 1. 534d 1. 520p.590d.570p 1. 038d 1.009e 1. 049p.769d.771e.796p Aug. 17,165,763b 1.625d 1. 614p.648d.606p 1. 046d 1.015e 1.054P.735d.770e.790p Sept. 17, 581,084b 1. 631d 1. 628p.708d ~ 653P 1. 059d 1. 021e 1. 061p.720d ~ 769e ~ 785p Oct. 17, 200, 676b 1. 606d 1o 611p.748d.689P!o 070d 1.025e!o 066p.716d ~ 769e ~ 783P Nov. 16,333,966b 1.52 6d 1. 548p.723d.690p 1. 090d 1. 028e 1. 071p.744d ~ 769e.782 p Dec 14,570, 85 lb 1. 432d 1. 446p.715d.691p 1. 089d 1.03Oe 1. 076p ~ 752d ~ 768e.781p 1953 Jan. 11,525,577b 1.253d 1.270P.711d.691p 1.099d 1.029e 1.081p.757d.767e.780p Feb 8,279,304b 1. 059d 1. 077P.710d.691p 1.107d 1.02 6e 1. 087p ~ 759d.766e.780p Mar. 6,437,259b.946d.965p.710d.692p 1.114d 1. 022e 1. 089p.759d.765e.780p Apr. 6,368, 997b.949d ~ 963p.710d ~ 693p 1.104d 1. 019e 1. 086p ~ 759d.764e ~ 778p May 6,812,015b.984d ~ 992p.711d.696p 1. 094d 1. O17e 1. 083p.758d.763e.776p June 7,274,647b 1. 019d 1. 022p.712d.699P 1o084d loO16e 1.081p.757d.762e.772p July 10, 02 O, 604b 1.19Od 1.19OP.714d.702p 1. 075d 1. O18e 1. 075P ~ 755d ~ 761e ~ 769p Aug. 12,692,443b 1.349d 1. 344p.714d.705p 1. 070d 1.O22e 1.075P.755d.760e.765p Sept. 14,014,979b 1.416d 1.414p.716d.708p 1. 076d 1. 027e 1. 078p.753d.759e.762p Oct. 14,432,099b 1.430d 1.432p.717d.711p 1. 084d 1. 031e 1. 082p.752d.758e.759P Nov. 13,519,224b 1.370d 1.376p.718d.713p 1.093d 1.034e 1.087p.751d.757e.757P Dec. 10, 999, 458b 1.217d 1.229p.72Od.714p 1.105d 1. 034e 1. 093P ~ 749d ~ 756e.756p i954 Jan. 7,909,681b 1. 030d 1. O44p.721d.715p 1.116d 1.032e 1. 099P.748d.755e.755P Feb. 6, 849, 802b.965d.981p.722d o715p lo120d 1o029e 1 o!OOp.747d.753e.754p Mar 5,774, 258b.904d.915p.723d.716p 1.117d 1.026e 1.102p.746d.752e.75~P Apr. 5,995,596b.923d ~ 931p.725d.718p 1.109d 1. 024 e 1. 099P.744d ~ 751e.752p May 7, 389,775b 1. 012d 1. 019p.725d.72Op 1.102d 1. O23e 1. O92p.744d.750e.75Op June 8, 630, 622b 1. 0968 1~ 095P.726d ~ 723p 1. 088d 1. 024e 1. 090p ~ 743d ~ 749e ~ 747p July 10, 2 98, 292b 1.194d 1.194p ~ 727d.726p 1. 088d 1.02 6e 1. 087p.742d.748e.743P Aug. 11,837, 297b 1.281d 1.282p ~ 729d ~ 729p 1. 089d 1. O29e 1. O87p ~ 739d ~ 747e.74Op Sept. i3,453,985b 1.367d 1.371p.729d.732p 1.092d 1.033e 1. 088p.740d.746e.737P Oct. 13,270,064b 1. 350d 1.356p.732d.734p 1.098d 1.037e 1. 093P.737d.745e.735P Nov. 12,211,907b 1.282d 1.293P.733d.736p 1.108d 1.039e 1.097P.736d.743e.733p Dec. 9, 266, 347b 1.102d 1.119p.734d.737P 1.124d 1.038e 1.104p.735d.742e.732p 1955 Jan. 6,746,445b.950d ~ 968p ~ 734d ~ 737P 1.134d 1. 035e 1.110p.734d.741e ~ 731p Yeb. 4,278, 932b.797d o812p.736d.737P 1.142d 1. 031e 1.119p.733d.740e.732p Mar. 3,251,797b.733d.748p.736d.737P 1.145d 1. 026e 1.121p.732d.739e.732p Apr. 2,820, 878b.708d.721p.738d.738p 1.144d 1. 022e 1.120p.731d.738e.731p May 2,820, 449b.716d.723p.739d.739P 1.129d 1. 019e 1.117p.729d.737e.729p June 2,820,192b.722d.725p.740d.741p 1.119d 1.016e 1.114p.728d.736e.727p July 2,820, 037b.725d.727p.747d.744p 1.114d 1.014e 1.111p.722d.735e.724p Aug. 4,818,371b.869d.86Op.744d ~ 737P 1. 085d 1. 013e 1. 097P ~ 725d ~ 734e ~ 72 lp Sept. 7, 616, 379b 1.045d 1. 038p.744d.750p 1. 078d 1.015e 1.086p.725d.733e.719p Oat. 8, 939, 003b 1.127d 1.116p ~ 745d.752p 1. 075d 1. 019e 1. O86p.724d ~ 732e.716p Nov 8, 516, 400b 1. 094d 1 o 088p o747d o754p 1.084d 1.021e 1.09Op.722d.731e.715p Dec. 6,005,686b.929d.930p.747d.754p 1.102d 1.021e 1.101p.722d.730e.714p 1956 Jan. 2,999, 464b ~ 730d ~ 735P.748d ~ 754p 1.125d 1.017e 1.116p.72 ld.728e.714p Feb. 1,503,130b.627d.637P.748d.755P 1.144d 1.012e 1.125p.721d.727e.714p Mar. 345,285b.546d.559P.748d.755P 1.163d 1. 007e 1.134p.720d.726e.714p Apr. - 164,269b.515d.526p.75Od.756p 1.161d 1.002e 1.136p.719d.725e.713p May - 367,194b.504d.513p.750d.757P 1.157d.998e 1.135P.718d.724e.712p June 1,521,965b.648d.648p.753d.759P 1.106d.996e 1.105p.716d.723e.710p July 4,004,281b.833d.818p.754d.761p 1.064d.997e 1.085p.715d.722e.708p Aug. 6,294,072b.991d.967p.754d.763P!.047d.999e 1.076p.715d..721e.706p Sept. 6, 345,981b.987d.968p.755d.765p 1.055d 1. 002e 1.079P.713d.720e.704p N 6, 343,531b.982d ~ 966p ~ 756d ~ 766p 1. 061d 1.003e 1.0SOp.713d ~ 719e ~ 703p ~ot. v. 7,105,189b 1o023d 1.014p.756d.767p 1.067d 1.005e 1.079P.713d.718e.702p Dec. 6,860, 494b 1. 007d.996p ~ 757d ~ 768p 1. 070d 1. 006e 1.082p.712d ~ 717e.701p ~ a denotes cumulative gas production, SCF ~ 60~F & 15o025 psia b denotes cumulative gas injected, starting April, 1948, SCF ~ 60~F & 14.735 psia d denotes quantity obtained from field data e denotes quantity predicted by Equations (III-20a), (III-2063, (IV-11), and (IV-12), which ignore interference effects p denotes quantity predicted from interference theory - from Equations (l-V-8), (ZV-15), (~-163, a~d (~-~73

mronths. Trie percentage error in predicted field D pressures, effected by use of the constant; volume method, would, approach 29% at 1.45 months after discovery since the actual field volume differs from the original volume by 29% at this time. Figure 27 shows that use of the'nterference method, described above, produces considerably better agreement between. the actual and predicted performances of fields C and D than that achieved by ignoring interference effects, The difference between the two field C predicted pore volume curves represents the decrease in field C's pressure (or increase in pore volume) caused by gas production in field D, The crossing of the two field D predicted pore volume curves at approximately 110 months after discovery is due to the gas injection and subsequent pressurizing of field C during the time period following t = 50 months. B, Prediction of Gas Storage Reservoir Performance When. Initial Aquifer Pressure is a Non-Uniform Function of Radius Th.he aqun'fer pressure at the time of discovery or initiation of a gas reservoir is, in, general, uniform or constant throughout the aquifer. In some cases, however, it is necessary to predict the performance of a gas bubble sit+uated on an aquifer when the initial aquifer pressure is a non-uniform function of radius, A typical case is a gas field from which little or poor data were obtained during a period immediately following the field discovery or 7initiatiaon, In such case thLJe predictive calculations must specify as 0 time some date wh1ich marks the end of the period during which poor field data were recorded, Another trypical case is that of an aquifer storage field, since if the predictive calculations specify the reservoir initiation date as O time

-136then these calculations must take into account the growing radius. If, however, time is considered O at some date after the radius has attained a relatively constant value, then the calculations are much simpler since they can employ equations derived from the assumption of a constant gas bubble radius. These equations must, of course, take into account the effect of a non-uniform initial aquifer pressure distribution, since the gas injection during the growth of the aquifer storage field will result in a developed aquifer pressure distribution at the selected date corresponding to 0 time. The literature contains no solutions to the diffusivity equation corresponding to the case of a source or sink (gas field) of finite radius surrounded by an infinite aquifer in which the initial pressure distribution is a non-uniform function of radius. An approximate solution for this case is developed below and the solution is then employed in predicting the performance of aquifer storage reservoir Eo 1, Approximate Solution to the Diffusivity Equation for the Case of a Non-Uniform Initial Aquifer Pressure Distribution An approximate solution to the previously presented diffusivity Equation (III-5) is presented here for the initial condition (IV-18) and boundary conditions (IV-19) and (IV-20). d {(5>oD) e / pP6(D! to) =; P~it)(III-5) /,';7" Py,r'.$t)- 0 (Iv-19)

-157r The initial condition (IV-i8) equates the initial aquifer pressure P(r DO) to P - frD), where Po is the constant aquifer pressure value at large radial distances from the center of the gas bubble. Condition (IV-19) is equivalent to the assumption of an infinite aquifer and (IV-20) specifies a variable water flow rate across the gas field boundary rD = r/rb = 1i. If the pressure drop P(rDtD) is set equal to a sum of two pressure drop components PfrDgWt) = Pf~ tD+ P;~tp) ~(IV-21) where Pl(rDtD) satisfies Equation (III-5) and conditions (IV-22) and P2(rD,tD) satisfies Equation (II1-5) and conditions (IV-23), then the sum of P1 and P2 will satisfy Equation (III-5) and conditions (IV-18), (IV-19), and (IV-20), as desired. l (IV-22) (;4(rD', D) = O The sum P1 + P2 represents only an approximate solution to Equation (iiI-5) for the conditions (IV-18), (IV-L9), and (IV-20) in the respect that P2(rDtD) ) does not approach 0 at increasing radial distances from the gas field center but is identically equal to O at rD > rD>

-138Solutions to the diffusivity equation (III-5) for boundary conditions similar to (IV-23) have been developed in the literature(9); the steps involved in the development of the solution P2(rDtD ) are therefore omitted here and the final solution is given in Equation (IV-24), (rv, >t an e 0 ( An 11) 5(Iv-24) n=n and Xn are the roots of the equation w oher ee)4(rn) Jo (An tc) qI (), (IV-25) obtained in Section III-A by application of Duhamel's superposition A (ATt) = A?. Z (2 ku - (-A/ )k (IV-26) where g(t) is related to the pore volumes V by where g(tD) is related to the pore volumes Vi by 9((t) ^ t (itD) 5 roat (Ad.- VA and tD = jatD = jK2At K _ K3 2chknt Setting rD = 1 in (IV-24) and (IV-26), one obtains the following expression for the gas bubble pressure drop at time tD, P(1,tD)'

-13 9 o=j-I o- P/_ = 37 ( -Ik -ti,) P-Q,, ihyZ)i4 e t/ n) (IV-27) Denoting the infinite summation 2 [ne 0 ) (or r.j n=, n-, Uo(An)) by Aj, and employing the previously presented Equation (III-20) to eliminate the term Vj from (IV-27), one obtains the pressureexplicit equation Pt P(2 t+ (IV-28) where Pt Pj _ P(l,tD) = gas bubble pressure at time t = jAt or tD = jAtD, psia, K6 = K6 + Aj, and K6 and K7 are defined on page 22 of this thesis. One fact of interest is readily apparent from Equation (IV-27), The error incurred by neglecting the term P2(1,tD).'ne n e(,n) will decrease as the value of tD, for a fixed value of time t, increases, since a larger tD results in a smaller value of e ntD and thus in a smaller value of P2(1 tD) at time t. Therefore, the error incurred by neglecting P2(1 tD) will be smaller as the aquifer formation mobility k/4 or the factor l/cpcrb are taken larger, since tD K2t k/l)(l/cpcrb.)t 20 Application to Aquifer Storage Reservoir E - Method A Gas field E is an. aquifer storage reservoir which was initiated in 1953 by injection of gas into an aquifer formation. original.l.y containing no gas~ Subsequent gas leakage from reservoir

-140E to shallower formations was indicated by increased gas vent rates from wells completed to small gas pockets in these shallow formations. Calculations performed in 1954 and later years indicated that the effective reservoir E gas content was significantly less than the cumulative inventory gas injected into the storage field. September 30, 1955 was specified as 0 time in the calculations described below since the gas bubble E radius had attained a relatively constant value at that date and the gas leakage had been significantly reduced in 1955 by the introduction of a gas recirculation system to remove gas from shallow formations and return it to the reservoir E. The calculations were made with the assumption of a 0 gas leakage rate or, in other words, with the assumption that the gas leakage out of the field E was replaced in equal quantity by injection of the shallow formation gas. The initial (September 30, 1955) aquifer pressure distribution was obtained from: observation well pressures and is plotted in Figure 29. The reservoir E pressure and gas in place data for the period September 30, 1955 - May 4, 1956 are listed in Table XI. The gas in place data listed in this table are inventory values which represent the net total cumulative gas injected into the reservoir, Calculations performed by the personnel of the operating company have indicated that the effective reservoir content is 4 - 12 billion standard cubic feet less than the inventory volumes listed in Table XI, The following K2 and K3 values were calculated from available data relating to the reservoir-aquifer system physical constants: = - k/.crb = O.3,X/o-g" //see. k3 p /27rArhk = / ift /O- ps'a/f7,3

Table XI SUMMARY OF ACTUAL AND PREDICTED PERFORMANCE OF AQUIFER STORAGE RESERVOIR E FOR PERIOD SEPTEMBER 30, 1955 - MAY 4, 1956 K = -.0001565 1/psia (where z = 1 + KP) Po = 715 psia T = Gas Field Formation Temperature = 70~F K2 = k/qccrb2 =60.35 x 10-5 i/seconds At =.606 x 10 seconds =1 week Predicted Pressures, Psia Inventory Observed Aj O Aj = Gas in Place, Mcf Reservoir K3 =.43 x 10-5 K3 =.46 x 10-5 Date (@ 60~F & 14.65 psia) Pressure, Psia "Gas Loss" = 12 M3cf "Gas Loss" = 12 M3cf 1955 Sept. 30 21,718,753 731.5 Oct. 7 21,921,162 735 726.4 724.2 Oct. 14 22,079,299 736.2 729.8 727.7 Oct. 21 22,271,891 738 733.4 731.7 Oct. 28 22,498,054 738.2 737.2 736.1 1955 Nov. 4 22,728,245 740 740 739.3 Nov. 11 22,806,731 736 735.8 735.3 Nov. 18 22,894,276 732.5 734.1 733.7 Nov. 25 23,155,528 739 740.1 740.1 1955 Dec. 2 23,020,604 728.5 727.9 727.8 Dec. 9 23,123,118 733 730.3 730.2 Dec. 16 22,905,135 721 718.4 718.1 Dec. 23 22,911,443 720 720 719.6 Dec. 30 23,211,928 734 732.2 732.2 1956 Jan. 6 23,404,623 737.4 735.2 735.6 Jan. 13 23,606,872 738.5 737.8 738.4 Jan. 20 23,843,593 742.5 741.2 742 Jan. 27 23,878,525 736 735.8 736.7 1956 Feb. 3 23,996,994 735.5 735.9 736.8 Feb. 10 24,274,831 741 742 743.1 Feb. 17 24,403,826 738 740.3 741.5 Feb. 24 24,594,969 739.6 741.9 743.2 1956 Mar. 2 24,787,774 741.2 743.1 744.5 Mar. 9 23,548,830 692 690.4 690.3 Mar. 16 23,483,865 704 700.1 699.7 Mar. 23 23,614,564 718 712.1 711.8 Mar. 30 23,971,217 733 727.6 727.9 1956 Apr. 6 24,289,732 742 736.1 736.8 Apr. 13 24,416,308 739 735 736 Apr. 20 24,499,116 735 733 734.1 Apr. 27 24,683,076 735 736 737.2 May. 4 24,869,274 738 738 739.2

735 730 o W 725 (0 U) w o. _ 720- 0 W IU 715 710 705 L I l 70~~~~~ I ~~~~I I I I I 2 3 4 5 6 7 8 9 10 r rD - rb Pressure Distribution in Aquifer Underlying Reservoir E. September 30, 1955 Figure 29

-143A time increment At of one week or 0.606 x 106 seconds was employed in the calculation of the above K3 value and in the calculations described below. The value of Po for reservoir E is 715 psia. The term DEV, defined as in Equation (III-25), serves as a measure of agreement between the field E observed and predicted pressures. The initial pressure distribution curve given in Figure 29 extrapolates to a value of O at an rD vaule of approximately 9 - 11. An rDe value of 10 was chosen for the purpose of calculating the value of the term ZQ anee A l"t/o(h) as a function of dimensionless time tD. The calculations were programmed in the IBM 704 FORTRAN code and performed by the 704 computer; the results are listed in Table XII and plotted in Figure 32. The values of E ahe U0 (xo for tD = K26t, 2K2At, 3K2at,..., 31K2Lt, where K2 = 0.35 x 10-5/ second and At = 0.60o6 x 106 seconds, were obtained by interpolation between the values listed in Table XII. The resulting interpolated values were then employed as the Aj terms in Equation (IV-28), where A3 -Zane g (Xh) For each of several assumed'gas loss' values, where'gas loss' = inventory gas content (from Table XI) - effective gas content, the gas in place values given in Table XI were reduced by the assumed'gas loss'. The resulting effective gas in place values were employed

-144TABLE XII P2(1 tD) AS A FUNCTION OF DIMENSIONLESS TIME tD FOR AQUIFER STORAGE RESERVOIR E Assumed rDe = 10 00 2 tD p2(1,tD) = n=1 an ent DUo(Xn)= Aj 0 - 16 2 - 4.954 4 - 3.518 6 - 2.762 8 2.280 10 - 1.937 12 - 1.674 14 - 1.463 16 - 1.286 18 - 1.134 20 1.003 22.888 24.787 26.698 28.618 30 0549 32.487 34 -.432 36.383 38 -.340 4o.301 42.267 44.237 46.- 210 48 -.187 50 - o166 52 -.147 54 - 130 56.116 58.102 60.091 62. 081

-145in Equation (TV-28) to calculate the predicted pressures, Pj (J =1, 2, 3,...., 31), for sufficiently many assumed K3 values to fix a minimum on a plot of DEV vs K3" A constant K2 value of 0.35 x 10 was employed throughout these calculations, which were performed by the IBM 650 computer, The resulting DEV values are listed in Table XIII as a function of K3 for assumed'gas loss' values of 11, 12, and 13 billion standard cubic feet. DEV is plotted vs K3 in Figure 30 with the'gas loss' value as a parameter. The reservoir E predicted pressures, for assumed'gas loss' and K3 values of 12 M3 cf and 0.43 x 10-5 psia/ft.3, respectively, are listed in Table XI and are compared graphically to the observed pressures in Figure 31, The predicted reservoir E pressures were also calculated with the term P2(!,tD) ignored, i,e., with all the Aj terms in Equation (IV-28) set equal to 0o This neglect of P2(l1tD) reduces (IV-28) to the predictive Equation ( II-22), which was developed for the case of a uniform (constant) initial pressure throughout the aquifer. A'gas loss' of 11 billion standard cubic feet was assumed and the predicted pressures Pj (j = 1, 2, 3,...., 31) were calculated from (IV-28) for sufficiently many assumed K3 values to fix a minimum on a plot of DEV vs K3o The resulting DEV values are listed in Table XIII and are plotted vs K3 in Figure 30. The predicted pressures obtained are graphically compared with the observed reservoir E pressures in Figure 31 and are 1listed in Table XIo

-1460.0039 38 GAS LOSS"= 11 x O10 f 3< 37 A _ i ~-, ~jA 0 36 aI 334 S GAS LOSS a 12 x 10 scf 35 o 34 IGAS LOSS =13 x 109scf 0.0032 0.41 x 10- 0.42 0.43 0.44 0.45 0.46 0.47 K3 Effect of Assumed K3 and "Gas Loss" Values on DEV Value for Aquifer Storage Reservoir E. Figure 30

Aau'A 740 730 \\!~~~~~~~~~~~~~~~~~~~~~~I. 720 I W 0 —-O PREDICTED PRESSURE \ I a.~ Aj t b -A J 710 K3=O.43 x10 > GAS LOSS =12x10 scf / a: OBSE RVED co PRESSURE Ir) -- PREDICTED PRESSURE I 70010 Aja 0 )-~~~~~~~~~~~~~aK3:0.46 x10/ 0- W GAS LOSS9 II x109 scf / 690 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 j NUMBER OF I-WEEK TIME INCREMENTS Figure 31. Comparison Between Aquifer Storage Reservoir E Observed and Predicted Pressures for Period September 50, 1955 - May 5, 1956.

To 16 01, o 0 5~~~~~~ CALCULATED PON 4 a.Uj~~~~4 56 40 4 24 0 s 16 LES TME to 0 DIMENSION-ES io~~aes Tmefor Aquifer Storage V Dilnlnsine Tm Fjigure 32

-149TABLE XIII DEV AS A FUNCTION OF ASSUMED "GAS LOSS" AND K3 VALUES FOR AQUIFER STORAGE RESERVOIR E 2 K2 = k/cpcrb = 0.35 x 10 5 1/seconds K = - 0001565/Psia At = 1 week =.606 x 106 seconds T = Gas Field Formation Temperature = 70~F Po = 715 Psia Results for Aj 0 Results for A = 0 Assumed "Gas Loss" = 11 MJcf (@ 60F & 14,.65 psia) Assumed "Gas Loss" Assumed Assumed M3cf ( 60~:F & 14.65 psia) K3 Value DEV K3 Value DEV 11.41 x 10-5.00365.45 x l0-5.00385 11.42 o00345.46 x 103.00382 11 e43 o00336.47 x 10-5.00383 11.44.00330 11.45.00336 12.42 x 10-5.00332 12 *43.00322 12.44.00324 13.41 x 10o5.00345 13.42.00333 13.43.00329 13.44.00336

-150The results listed in Table XIII show that inclusion in Equation (IV-28) of the pressure drop term accounting for the aquifer initial pressure distribution, P2(itD), results in nearly a 16% decrease in the average fractional deviation (DEV) between the observed and predicted pressures (from 0.00382 to 0i00322 for an assumed'gas loss' of 12 billion standard cubic feet). The definite minimums in the curves of DEV vs K3 in Figure 30 show that a minimum DEV value serves as a good criteria for choice of the K3 value when K2 is fixed. The above described Method A of treating cases involving a non-uniform initial aquifer pressure distribution may prove impractical in several respects. First, the value of rD must be e known or assumed and this value is, in general, unknown. Second, the calculation of the term P2(1,tD), required in application of Method A, is very lengthy since the roots Xn of Equation (IV-25) m 2 t must first be calculated, and the an,Uo(Xn), and e n D values must then be determined and inserted into Equation (IV-24) in order to determine P2(l tD) The practicality of Method A is further restricted by the requirement that the initial aquifer pressure distribution be known. Gas storage reservoir-aquifer systems may or may not have a sufficient number of aquifer observation wells to allow determination of this initial distribution. 3. Application to Aquifer Storage Reservoir E - Method B The alternative Method B now presented requires neither qualitative nor quantitative knowledge of the initial pressure distiribution. and requires no specification of the value of rD.

-151This method does require krnowledge of the K2, K3, and Po (where P0 is the constant aquifer pressure approaced at large radial distances from the gas reservoir) valueso The K2 ( = - ) and K3 ( = i ) values can be obtained from available data relating to the aquifer formation and water properties or from calculations similar to those described in Section III-C of this thesis. The terms Aj, calculated above as A. = D ne U n=I can be calculated directly from field data by employing a rearranged form of Equation (IV-27) 6 = j: Paj ~ -,k52Z(20-iAt,-V J)fpr(. 4' zt) (IV-3o) where Paj = observed gas reservoir pressure at time t - jZt, psia Vaj = actual gas reservoir pore volume at time t = iAt, ft.3. Suppose the predicted performance of a gas field is desired for a time period of (N-M) time increments, where M and N are both integers and N > Me Let M be such that 0 time is specified at M time incremernts prior to the beginning of the period for which the predicted performance is desired0 If field data is available for the period from t = to t = MIt then Aj can be calculated from Equation (IV-30) with j =, 2, 3, 0 0, M. Employing these Aj values for j < M, Aj values equal to 0 for j > M, and the projected gas in place schedule for the period from t - (M+i)At to t - Nat, one could then calculate Pj from Equation (IVP28) with j = I, 2, 3, O.a, No However, the P values calculated for j = 1, 2, 3,...., M, will simply be the Paj values employed in ca2l.cuating Aj from Equation (IV-30)o Thus there is no need to calculate the terms Ai (from Equation IV-30) and Pi

-152(from IV-28) for j = 1, 2, 3,...., M; the identical predicted pressures Pj (j = M+l, M+2,...., N) will be obtained from (IV-28) by simply setting Pj = Paj and Vi = Vai for j or i < M. setting Aj = 0 for j > M, and beginning the calculations with j = M+1 in (IV-28). The approximation Aj = 0 for j > M will be valid in direct proportion to the value of tD = MAtD since The method B outlined immediately above has been applied to gas storage reservoir E. The validity of the claim that the Aj values calculated from field data and Equation (IV-30) will decrease with increasing time was first tested. Pressure-injection data for the period September 1, 1956 - December 15, 1956, a time increment of three days, and K2, K3, and'gas loss' values of.35 x 10-5 1. 08 x 10-5 and 12 M3 cf, respectively, were employed in calculating Aj from Equation (IV-30)o The resulting Aj values are listed in Table XIV with the gas field E data for the period September 1, 1956 - December 15, 1956. Figure 33 shows the decrease of Aj with increasing time. Method B has been applied in calculating the predicted performance of aquifer storage reservoir E for the period December 18, 1956 - March 6, 1957. September 1, 1956 (36 three-day time increments prior to December 18, 1956) was specified as 0 time and the Pj and ni values for j or i < 35 were taken from the data listed in Table XIV. The predicted Pj values were then calculated from Equation (IV-28) (with Aj = O, j. 36) for j = 36, 37,....= 62 (March 6, 1957). A three-day time increment and K2 and'gas loss'

-153Table XIV SUMMRY OF RESERVOIR E DATA FOR PERIOD SEPTEMBER 1, 1956 - DECEMBER 15, 1956 Aj Mcf Inventory Calculated from (IV-30) Gas in Place K2 =.35 x 10-5 Observed (@ 14.65 psia & 60~F) K3 = 1.o8 x 10-5 Reservoir (Includes Known 2.9475 "Gas Loss" = 12 M3cf Date Pressure, Psia M3cf Vent Loss) At =.25925 x 106 sec. 1956 Sept. 1 735 23,904,133 Sept. 4 735 23,929,791 - 16.2 Sept. 7 734.2 23,946,195 - 14.0 Sept. 10 734.5 23,961,023 - 16.o Sept. 13 735 23,981,894 - 16.2 Sept. 16 734.2 23,999,853 - 13.0 Sept. 19 734 24,016,902 - 13.4 Sept. 22 734 24,029,576 - 14.3 Sept. 25 734 24,073,009 - 9.8 Sept. 28 734 24,104,968 - 10.3 1956 Oct. 1 734 24,149,612 - 8.0 Oct. 4 735.2 24,202,403 - 9.6 Oct. 7 735.5 24,254,215 - 7.9 Oct. 10 736 24,299,321 - 9.1 Oct. 13 737 24,339,011 - 11.7 Oct. 16 736 24,363,984 - 8.9 Oct. 19 735.4 24,393,082 -8. Oct. 22 735.3 24,426,468 - 8.1 Oct. 25 735.2 24,471,546 - 6.1 Oct. 28 735 24,515,877 - 5.2 Oct. 31 735.2 24,564,983 - 5.0 1956 Nov. 3 735.7 24,616,976 - 5.4 Nov. 6 736.2 24,664,850 - 6.3 Nov. 9 736 24,690,894 - 7.7 Nov. 12 735.8 24,735,508 - 5.1 Nov. 15 736.5 24,777,252 - 7.6 Nov. 18 736 24,807,763 - 6.3 Nov. 21 736 24,824,767 - 9.3 Nov. 24 733 24,831,932 - 1.7 Nov. 27 733.5 24,838,476 - 9.2 Nov. 30 733.5 24,879,548 - 4.4 1956 Dec. 3 734 24,928,980 - 4.1 Dec. 6 735.5 24,987,209 - 6.1 Dec. 9 734 25,001,096 - 4.2 Dec. 12 734.5 25,001,751 - 11.1 Dec. 15 732.5 25,021,261 - 2.2 Dec. 18 733 25,052,917 Dec. 21 734.5 25,109,917 Dec. 24 735 25,152,600 Dec. 27 735.5 25,192,432 Dec. 30 735.5 25,233,545

-20 -18 X-16 | ~\ o o 0Aj VALUES CALCULATED FROM eq. IZ -30 0 I 0 -14 0 0 - -12 ~ I: -10 - Ow~~O H W -8 G ~ ~ ~ 4 4..J 0I 0I I jNUMBER O 3-DAY TIE INCREMENTS 4 - -2 en 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 j,NUMBER OF 3-DAY TIME INCREMENTS Aj vs. Time for Aquifer Storage Reservoir E, September 1, 1956 - December 15, 1956. Figure 33

-155values of ~35 x 10 5 and 12 M3 cf, respectively, were employed in the calculations0 The December 18, 1956 - March 6, 1957 gas in place data used in the calculations are listed in Table XV, along with the observed and predicted gas bubble pressures, The value of K3 was varied until the value of DEV, defined in Equation (III-25) (with N = 62-35 = 27), passed through a minimum. The DEV values obtained are plotted versus K3 in Figure 34 and the predicted reservoir E pressures are compared to the observed reservoir pressures in Figure 35. The reservoir E predicted pressures for the December 18, 1956 - March 6, 1957 period were also calculated, ignoring the initial pressure distribution, from Equation (III-22). December 15, 1956 was specified as 0 time and a three-day time increment, and K2 and'gas loss' values of.35 x 10-5 and 12 M3 cf, respectively, were used in the calculations. The K3 value was varied until DEV, defined in (III-25) (with n = 27), passed through a minimum, DEV is plotted versus K3 in Figure 34 and the predicted reservoir pressures calculated from (III-22) are listed in Table XV and plotted in. Figure 35o Figure 34 shows that use of method B yields significantly better agreement between the observed and predicted pressures than that obtained by ignoring the initial pressure distribution. C, Treatment of the Moving Boundary Problem Encountered in Aquifer Storage Operations As pointed out previously (Section I), the distinguishing characteristic of an aquifer storage' gas field is the growth of the field radius or

-156Table XV SUMMARY OF RESERVOIR E OBSERVED AND PREDICTED PERFORMANCE, DECEMBER 15, 1956 - MARCH 6, 1957 Predicted Pressure, Psia "Gas Loss" = 12 M cf K2 =.35 x 10-5 6 At = 3 days =.25925 x 10 sec. Mcf Inventory Po = 715 Psia Gas in Place (@ 60~F & 14.65 psia) Observed Ignoring Initial (Includes Known Reservoir From Method B Pressure Distribution Date 2.9475 Vent Loss) Pressure, Psia (K3 =.63 x 10-5) (K3 =.55 x 10-5) 1956 Dec. 15 25,021,261 732.5 Dec. 18 25,052,917 733 728.5 725.7 Dec. 21 25,109,917 734.5 727.8 724.2 Dec. 24 25,152,600 735 727.1 723,2 Dec. 27 25,192,342 735.5 726.6 722.7 Dec. 30 25,233,545 735.5 726.3 722.4 1957 Jan. 2 25,164,139 731.5 722.1 718.5 Jan. 5 25,170,410 728.5 721.9 718.5 Jan. 8 25,211,675 733 722.9 719.6 Jan. 11 24,951,832 723 712.4 709.7 Jan. 14 24,409,203 713 694.1 692.7 Jan. 17 23,641,655 687.5 671.1 671.7 Jan. 20 23,746,900 688.3 688.5 689 Jan. 23 23,829,446 702.5 698.5 698.2 Jan. 26 23,606,157 700.5 692.5 692.8 Jan. 29 23,464,763 698 692 692.5 1957 Feb. 1 23,581,566 702 701.9 702 Feb. 4 23,813,816 718 713.1 712.4 Feb. 7 24,075,535 728.5 722.0 720.5 Feb. 10 24,129,902 729.5 720.6 718.9 Feb. 13 24,219,643 731 721.8 719.8 Feb. 16 24,307,362 733.5 722.8 720.7 Feb. 19 24,326,681 734.5 721.3 719.3 Feb. 22 24,405,428 733 722.7 720.6 Feb. 25 24,466,676 734 723.0 720.8 Feb. 28 24,503,658 733.5 722.4 720.3 1957 Mar. 3 24,560,432 732.8 722.9 720.7 Mar. 6 24,637,685 733 724.1 721.9

1.50p 0 x i.40 a DEV VALUES CALCULATED IGNORING INITIAL PRESSURE DISTRIBUTION a? ~ 1.30 - DEV VALUES CALCULATED C~ 1'3 0 EMPLOYING METHOD B N 1- ~ | 1.20 II 3 iW| "GAS LOSS"= 12 M cf K =.35 x I0- 1.10 K 35X0o o | t 3 DAYS =.25925XlO6 SECONDS 0 — _ P-Po= 715 PSIA 1.00 I I I I I 0.40 0.45 O. 50 0.55 0.60 0.65 0.70 K 3X 105 Figure 34. DEV vs. K3 for Aquifer Storage Reservoir E, Period December 15, 1956 - March 6, 1957.

740 730 0. &A ~00.00.0 0 "%t 001 0 720 A ~710 g3 zz- OBSERVED RESERVOIR.. 700 PRESSURE CL 7000 i O -0-0 PREDICTED PRESSURES > 690 oCALCULATED BY METHOD B W co cn LU PREDICTED PRESSURES W CALCULATED IGNORING A INITIAL PRESSURE o 670 DISTRIBUTION 660 I 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 j,NUMBER OF 3-DAY TIME INCREMENTS Figure 55. Comparison Between Aquifer Storage Reservoir B Observed and Predicted Pressures, December. 15, 1956 - March 6, 1957.

-159volume w:ith.. increasin!.g time. The radius grows from a value of 0 at the time of init"iation of the storage reservoir to values of several thousand feet. This increasing rb (gas bubble radius) Complicates the solution of the diffusivity Equation. (III-5) since the term K2 -k/4pcrb which appears in tD (tD = K2t) is no l.onger constant but is a function of time, The solution of Equation. (iTI-5 for this case is common.ly ca ied a!moving boundary problemT' A general method of solving analytically a true moving boundary probl.e.s* has not been devised to date. The difficulty of solution lies in the fact that the diffusivity or governing differential equation cannot be sol.ved urntil t;he time variation of the movin.g boundary position is known., but thi-~s boundary posItiton is n.ot known as a function of time until the differential equation has been solved. The literature(9) does contain analytical solutions to the d-ffusivity equation for cases involving moving line or moving point sources when the moving line or point position is known as a function of time. A numerical. procedure(3) has also been presented for the solution. of a liquid1i quid disp acemen.t problem in a square, where one liquid is injected at one corner and tPhe oth.er liquid is removed at an equal rate at the diagonally opposite corner. The moving source solutions are not applicable to the case considered here: the displacement or flow of a liquid away from a moving circular boundary at which. a time-variant pressure is maintained. The numerical procedure can simulate only finite aquifers, entails a tremendous volum.e of calcu.lations for an aquifer of reasonable extent (eog. a ratio of 1.0 or greater between. the aquifer exterior and interior radii), and is incapable of yieldi ng general solutions applicablie to any gas field. * A "true movin.g boundary problem.: is one in. which the movement of the boundary is not known a priori as a function of time~

The problem of predicting aquifer storage reservoir performance (during the initial growth stage of the gas bubble) is considered below from three points of view. First, the agreement between predicted and actual performance which can be obtained by ignoring the variation of rb with time is determined. Second, an analytical approach to the moving boundary problem is developed and the resulting method is applied in a case study where the moving boundary position is known as a function of time. Third, the developed analytical approach to the moving boundary problem is employed in predicting the performance of an aquifer storage reservoir when the moving boundary position is not known a priori as a function of time. The reader should recall that, as pointed out in Section IV-B, above, the difficulty presented by the growing rb of an aquifer storage reservoir can be avoided entirely if the gas bubble is already grown and a prediction of future performance is desired, For, in this case, 0 time can be specified at a date when rb has attained a relatively stable value and method B of Section IV-B can be employed to account for the non-uniform pressure distribution throughout the equifer at that date. 1. Prediction of Aquifer Storage Reservoir Behavior with rb Assumed Constant Calculations described in Section III-C of this thesis yielded the predicted performance of aquifer storage reservoir B for the period November 30, 1957 - July 5, 1958. The radius rb was assumed constant in these calculations, i.e. the term K2 = k/Icp~crb was assumed to be independent of time, The calculations showed that no unique combination of K3 and K2 values resulted in best

-.161agreemen+t'-oetween predicted and act. ual performance but that this agreement.'. improved as K3 was decreased and K2 increased~ However, comparison between, the predicted and actual. pore vo.umes, plotted in Figures 2- and 22, shows that the agreement obtained is satisfactory for engineering purposes. The observation. was made in Sectioln I1>.C that the actual pressures in field B were increasing less rapdl'y than the predicted pressures; thus the question arises as -to whether the discrepancy between the actual performance of field B and the predicted performance (predicted with the assumption of a constant rb), for a time period starting from July 5, 1958, will become sufficiently large to invalidate the predicted performance for engineering purposes~ In order to answer this question, the predicted performance of field B has been calculated from Equation 11>.L-22) (which was derived from the assumption of a constant rb) and compared to the actual performance over the period July 5, 1958 - November 29, 1958. A time irrement of one week and K3 and K2 values of 23 x 10 5 and 1000, respe tivel.y, (the same values employed in ca.c.ul.ating the results plotted in Figure 19) were used in the cal.cu-.atio.ns... Gas in p-lace ard pressure data for the July 5, 1.958 - November 29, 1958 period are l.isted ir Table XVIo The gas in place data from Tables V anud XVI were empoyed in calculating the predicted pressures from tB..e variablesra+te Equation (III-22). These predicted pressu-res and: the gas in. place data were ther employed irn Equat on'T:1T,20) to obtain. the predicted pore volumes. The observed pressures Tabt.les V and Efb) and gas in palce data were employed in

-162TABLE XVI SUMMARY OF PESERVOIR B OBSERVED AND PREDICTED PERFORMANCE, JULY 5, 1958 - NOVEMBER 29, 1958 MMcf Gas in Place'* Reservoir Gas Bubble Pore Date @ 60~F & 14.65 Psia Pressure,* Psia Volume* MMcf 1958 July 5 1116.565d 1242o7d 1267.2p 10.88d 10.62p July 12 1212.612d 1247 d 1276.6p 11.76d 11.53p July 19 1300.974d 1242.6d 1278.5P 12.68d 12.24p July 26 1394.177d 1244.6d 1282.7p 13.55d 13.06p 1958 Aug. 2 1492.338d 1244.ld 1288.4p 14.52d 13.90p Aug. 9 1581.520d 1239.7d 1288.0p 15.42d 14.74p Aug. 16 1979.101d 1241.3d 1292.Op 16.38d 15.59P Aug. 23 1783.678d 1241.7d 1298.2p 17.30d 16.46p Aug. 30 1861.472d 1231.4d 1291 p 18.30d 17.30p 1958 Sept. 6 1994.792d 1250 d 1309.5p 19.28d 18.21p Sept. 13 2132.184d 1258.4d 1324.9p 20.35d 19.18p Sept. 20 2257.872d 1258.4d 1331.9P 21.64d 20.18p Sept. 27 2378.619d 1256.6d 1335.6p 22.80d 21.19p 1958 Oct. 4 2509.062d 1259.1d 1342.0p 24.06d 22.21p Oct. 11 2642.448d 1259 9d 1348. Op 25.26d 23.26p Oct. 18 2758.985d 1254.2d 1347.3P 26.42d 24.31p Oct. 25 2877.581d 1249.6d 1347.6p 27.81d 25.34p 1958 Nov. 1 2997.308d 1251.4d 1348.2p 28.89d 26.38p Nov. 8 3127.890d 1254.9d 1352.0p 30 o09d 27.44p Nov. 15 3233.416d 1253.5d 1347.9P 31.06d 28.47p Nov. 22 3362.744d 1255.3d 1351.3P 32.31d 29.51p Nov. 29 3419.751d 1228.6d 1334.6p 33.74d 30.49p "d" denotes quantity obtained from field data. * "p" denotes quantity predicted for K5 = 2.3 x 10-5, K2 = 1000

Equation (Ii-,20) to obtain the actual pore vQlumes. All the calculated quantities are tabulated in Table XVI; Figures 36 and 37 compare the actual and predicted pressures and actual and predicted pore volumes, respectively. The predicted pressures and pore volumes for the November 30, 1957 - July 5, 1958 period are identical to the predicted values tabulated in Table V and plotted in Figures 19 and 21. The value of DEV, calculated from Equation (III-25) (with N = 52 since there are 52 weeks between November 30, 1957 (initiation of reservoir B) and November 29, 1958), is 3.27%. That is, if one takes an average pressure of 1200 psia for field B over the first 52 weeks of operations then the average difference between the predicted (assuming a constant rb) and observed pressure is.0327 x 1200 or 39.2 psia. Figure 36 shows that a difference of 106 psia occurs between the predicted and observed field B pressures 52 weeks after initiation of the gas bubble, This error of 803% in the actual 1228 psia pressure invalidates the predicted pressure curve for engineering purposes. However, the predicted pore volume curve, shown in Figure 37, allows a useful estimation of the bubble size as a function. of time. For example, suppose the operators of storage field B, after 2 or 3 weeks of gas injection, estimated the injection schedule for the first year of operations and waited to know the period of time required to grow the bubble to a pore volume of 30 MMcf, The predict i've method of calculation employed here would have given them an answer of 51L weeks, which Figure 37

-1641360 O...-.-0 —'0-0O 1320 O 1280 0o0o~0-0~-00 W. c) 1200 o X- |.. * -*- OBSERVED PRESSURES | cr cU a 1160 - - Tr | -0- PREDICTED PRESSURES oa co 1120 I: W 0 Z: 1080 1060 PSIA = INITIAL AQUIFER PRESSURE 1040,. 1040o- I I I I I I I I I I 31 33 35 37 39 41 43 45 47 49 51 53 j,NUMBER OF 7-DAY TIME INCREMENTS (j=O AT NOV.30,1957) Figure 36. Comparison Between Aquifer Storage Reservoir B Predicted and Actual Pressures, July 5, 1958 - November 29, 1958.

-16534 *O$- PORE VOLUMES CALCULATED FROM DATA 32 -0-0-0- PREDICTED PORE VOLUMES -5 30 K3= 2.3 x 10 / 0 K = 1000 0 28 / 26 / o 0 r / -= ~ 24 0 22) w ~~~~/ 0 D 22~ J /O IJ~~ 0 20 C 0D /0O 0) 0' Z 16 0 ~ /~~ 14 2 / ~~~~~/10 /9 ~ ~ ~ 31 33 35 37 39 41 43 45 47 49 51 53 j,NUMBER OF 7-DAY TIME INCREMENTS (jxO AT NOV.30,1957) Figure 57. Comparison Between Aquifer Storage Reservoir B Predicted. and Actual Pore Volumes, July 5, 1958 - November 29, 1958

shows is + 2~u weeks in error ~ Thus t,.e ass'imption of a rons'tan:.:-' rb and use of the predictive equat.ions devel.oped in Secton. _:A provides a useful means of estimatirg the size of ar agu-fle, r stoc age reservoir as a function of tim.e The discussion. of feld B in Se't io,-j T7 i.>-'2 and 4 - t-io section has attributed any gross diufferen-e bet.wee:,. ihe, a- ua. ar:. predicted performance to the in.creasi.g bubbie radius, However. the observed difference, shown in Figures 36 and. 37, 7ta:r be at t:,r. i ted to two other factorso T'he first is leakage of gas t;ro'gS: f.a:_]sts or permeable streaks in the fie.ld B cap rook, IL._e ffectFe of gas leakage would be that of decreasing ti'e observed pressAres or -~ r creasing the pore volume, since P 1./', by:.e gas law) be.'ow +he pressure corresponding to no leakage. Sirnce F`9uire 36) shows t,.:-at the actual pressure does fall below thJe predict-ed values, t::he possibility of gas leakage must be considered, T7Ie second> fac-t.;or is vertical pressure penetration` into'the aquifer forma;io:ro he predictive calculations assumned a circu-lar dis:k f.low'rodel hoI.ded on the top and'bottom r.y impervious pl.ares,. Howeter, a psress u. re test made in fall of 1958 iLndi.caied fluid pressure com.-nx.. -.aT onyr between the Mt. Simon and next lower formationr, wi. had prL-I viously been considered separated ky an imperrea l..e slhaLe s';rea.< The effect of pressure comLmulnir;at;ion. to lower strat:a wo;..d.:- a. increased pore volume growth rat'e sir.;nce the wat-er c re...eede I:.. the vertical as well as horizorn.tal directions~ Sinc::e thle a-.t.ual. growth rate exceeds the predicted rat e, as shown r: F i~g.re 37 the possibi.ity of vert^ical pressure penetration.w-s: a:lso >e

-167considered, The possibility of gas leakage and vertical pressure penetration is discussed further in Section IV-C-3, below. An example is now presented of how the water influx calculations can provide information useful in the specification of gas injection and withdrawal during the growth of an aquifer storage reservoir. The Section III-A equations, derived from the assumption of a constant rb, are employed in this example. The question considered is a very practical and important one concerning the effect of gas withdrawal during the initial growth stage of the gas bubble. The question is: what difference might be expected between the pore volumes and gas content of two aquifer storage reservoirs, separately grown in identical aquifer formations, at the end of a given number of years during which the first reservoir was grown under a constant bubble pressure and the second was grown under the same constant bubble pressure except during an annual 3 to 4 month gas withdrawal period? This question is of practical importance because storage field operators would like to compare the economics of withdrawing and selling gas from a reservoir as it is grown with the economics of withdrawing and selling gas only after the field is grown to the desired size. Calculations have been performed to determine the aquifer storage reservoir performance effected by two different pressure schedules, The following reservoir properties were assumed:

-168infinite aquifer Po = 700 psia K3 = 2~27 x 10-5 psia/ft.3 K2 = ~89 x 10-3 seconds-1 along with a time increment of two weeks in calculating gas in place quantities from Equation (III-24) (when the pressure was specified) and in calculating pressures from Equation (III-22) (when the gas in place was specified). The first case considered was the growth of the reservoir under a constant bubble pressure of 900 psiao The bubble pressure was increased in four 50 psia increments from 700 to 900 psia and was maintained at the latter value for 186 weeks, The predicted gas in place quantities were calculated from Equation (III-24) and were then employed in Equation (III20o) to yield the pore volumes. The specified pressure schedule and the calculated gas in place quantities and pore volumes are listed in Table XVII and plotted in Figure 38. The second case treated was the growth of the reservoir at a constant bubble pressure of 900 psia, except during the last 12 weeks of each successive 48 week period when 20% of the current gas in place was withdrawn. Equations-(III-22), (III-24), and (III-20) were again employed to calculate the pressure (when the gas in place was specified), the gas in place (when the pressure was specified) and the pore volumes, The calculated pore volumes and specified or calculated pressures and gas in place quantities and pore volwumes are listed in Table XVII and plotted in Figure 38.

-169Table XVII SUMMARY OF CALCULATED PERFORMANCE OF AQUIFER STORAGE RESERVOIR WITH AND WITHOUT GAS WITHDRAWAL DURING THE GAS BUBBLE GROWTH STAGE K3 = 4/2Ahkt = 2.27 x 10-5 psia/ft3 K2 = k/icpcrb =.89 x 10-3 secondsAt = 2 weeks Case 1: Gas Bubble Grown at Constant 900 psia Pressure Case 2: 20% of Current Gas in Place Withdrawn During Last 12 Weeks of Each 48 Week Period j, MMMcf Gas in Number of Gas Bubble Place*@ 14.65 MMcf Gas Bubble Two-Week Pressure, Psia Psia & 60OF Pore Volume Time Increments Case 1 Case 2 Case 1 Case 2 Case 1 Case 2 0 700 700 0o 0 0 1 750s 750s.034.034.605.605 2 800s 800s.108.108 1.76 1.76 3 850s 850s.226.226 3.44 3.44 4 900s 900s.395.395 5.63 5.63 5 900s 900s.541.541 7.71 7.71 6.682.682 9.73 9.73 7.821.821 11.69 11.69 8.956.956 13.63 13.63 9 1.o9 1.o9 15.53 15.53 10 1.22 17.41 17.41 11 1.35 1.35 19.26 19.26 12 1.48 1.48 21.10 21.10 13 1.61 1.61 22.92 22.92 14 1.74 1.74 24.73 24.73 15 1.86 1.86 26.52 26.52 16 1.99 1.99 28.30 28.30 17 2.11 2.11 30.07 30.07 18 2.23 2.23 31.83 31.83 19 848.2 2.36 2.16s 33.58 32.96 20 808.1 2.48 2.09s 35.33 33.64 21 775.6 2.60 2.01s 37.06 33.99 22 748.5 2.72 1.94s 38.78 34.09 23 725.3 2.84 1.86s 40.50 33.96 24 705.1 2.96 1.79s 42.21 33.66 25 900s 3.0o8 2.51 43.91 35.77 26 3.20 2.65 45.61 37.71 27 3.32 2.78 47.30 39.58 28 3.44 2.91 48.99 41.40 29 3.56 3.03 50.67 43.19 30 3.67 3.16 52.34 44.96 31 3.79 3.28 54.01 46.71 32 3.91 3.40 55.68 48.45 33 4.02 3.52 57.34 50.17 34 4.14 3.64 58.99 51.88 35 4.26 3.76 60.64 53.58

-170Table XVII (CONT'D) SUMMARY OF CALCULATED PERFORMANCE OF AQUIFER STORAGE RESERVOIR WITH AND WITHOUT GAS WITHDRAWAL DURING THE GAS BUBBLE GROWTH STAGE jH MMMcf Gas in Number of Gas Bubble Place @ 14.65 MMcf Gas Bubble Two-Week Pressure* Psia Psia & 60~F Pore Volume Time Increments Case 1 Case 2 Case 1 Case 2 Case 1 Case 2 36 900s goos 4.37 3.88 62.29 55.28 37 900s 900s 4.49 4.oo00 63.93 56.96 38 900s 900s 4.60 4.12 65.57 58.64 39 4.72 4.23 67.21 60.32 40 4.83 4.35 68.84 61.98 41 4.95 4.47 70.47 63.64 42 5.o6 4.58 72.09 65.30 43 860.7 5.17 4.43s 73.71 66.47 44 826.3 5.29 4.28s 75.33 67.27 45 795.6 5.40 4.12s 76.95 67.75 46 767.7 5.51 3.97s 78.56 67.95 47 742.0 5.63 3.82s 80.17 67.90 48 718.0 5.74 3.67s 81.77 67.63 49 900s 5.85 4.89 83.38 69.62 50 5.96 5.02 84.98 71.46 51 6.o8 5.14 86.57 73.23 52 6.19 5.26 88.17 74.97 53 6.30o 5.38 89.76 76.68 54 6.41 5.50 91.35 78.37 55 6.52 5.62 92.94 80.o4 56 6.63 5.73 94.52 81.71 57 6.74 5.85 96.10 83.36 58 6.86 5.97 97.68 85.oo 59 6.97 6.o8 99.26 86.64 60 7.08 6.20 100.84 88.27 61 7.19 6.31 102.41 89.90 62 7.30 6.42 103.98 91.52 63 7.41 6.54 105.55 93.13 64 7.52 6.65 107.12 94.74 65 7.63 6.76 108.68 96.35 66 7.74 6.87 110.24 97.95 67 865.1 7.85 6.65s 111.80 99.13 68 833.2 7.96 6.42s 113.36 99.95 69 803.6 8.07 6.19s 114.92 100.47 70 775.8 8.17 5.96s 116.48 100.70 71 749.6 8.28 5.83s 118.03 100.68 72 724.7 8.39 5.50s 119.58 100.43 73 g00s 8.50 7.18 121.13 102.36 74 8.61 7.31 122.68 104.14 75 8.72 7.43 124.23 105.86 76 8.83 7.55 125.77 107.55 77 8.94 7.66 127.31 109.21 78 9.04 7.78 128.85 110.85

-171Table XVII (CONT'D) SUMMARY OF CALCULATED PERFORMANCE OF AQUIFER STORAGE RESERVOIR WITH AND WITHOUT GAS WITHDRAWAL DURING THE GAS BUBBLE GROWTH STAGE j, MMcf Gas in Number of Gas Bubble Place*@ 14.65 MMcf Gas Bubble Two-Week Pressure, Psia Psia & 60~F Pore Volume Time Increments Case 1 Case 2 Case 1 Case 2 Case 1 Case 2 79 900s 900s 9.15 7.89 130.39 112.48 80 9.26 8.01 131.93 114.10 81 9.37 8.12 133.47 115.71 82 9.47 8.23 135.00 117.32 83 9.58 8.35 136.54 118.92 84 9.69 8.46 138.07 120.51 85 9.80 8.57 139.60 122.09 86 9.90 8.68 141.13 123.68 87 10.01 8.79 142.66 125.25 88 10.12 8.90 144.18 126.83 89 10.23 9.01 145.71 128.40 90 10.33 9.12 147.23 129.96 91 867.4 10.44 8.82s 148.75 131.96 92 836.8 10.55 8.51s 150.28 131.96 93 807.9 10.65 8.21s 151.80 132.49 94 780.3 10.76 7.90s 153.31 132.74 95 753.9 10.87 7.60s 154.83 132.73 96 728.4 10.97 7.30s 156.35 132.48 * "s", following number, denotes specified quantity * Every number not followed by "s" was calculated from Equation (III-20), (III-22), or (III-24).

-172CASE I. GAS BUBBLE GROWN AT CONSTANT BUBBLE PRESSURE OF 900 Psia CASE2. GAS BUBBLE PRESSURE MAINTAINED AT 900 Psia EXCEPT DURING LAST 12 WEEKS OF EACH 48 WEEK PERIOD WHEN 20 % OF GAS IN PLACE IN PLACE IS WITHDRAWN. -S 3 K3 a 2.27 x 10 Polo/ft. 2 vh k At K2 k 0. 89 x I16SIc'. /*+ c rb2 At a 2 WEEKS 00000 GAS BUBBLE PORE O00 VOLUME 140 PO * 700 Pio. 000000000 0o00 AA"" 0 CASE I 120 - 00000 0,, A A CASE 2 w 1200 ADA 10 oOOOOOAAA o 0000AAAA w 800 000000lA 1i 00 A&~~~00 A00 o AAAAAAAA.2~~. >*~~~~~~ O A,, o O 0 G AAA.IPA 10 0000000 0,AA AAAA 0 R~ ~oo00 AAA- CA00S A2 w 60 - 00 00 00000 AAAAAA 0j w 000AAAAAAA c 6 A0oo AA AA CD oO AAAA A m 40 -O 0 AAA AAA ~~~~~OOo AAA. ov Ao U) if~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ mD.oOO.(~AAAO 4 20 m o 00 o w~~~~ I90 900 oC0 GAS IN PLACE U) ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~o 10 0000 W O~~~~~~~~0000 o~~~B5O~~~~~~ *oooooo~ ~~ Ia~~~~~~ I - CASEI * I I * I A A* —-CASE 2 LJW 00000 A A A, AE oBOO I I I 6 O00 ~~~~~~~~~~h 00 A 0.1 4 \*1aa 000 AAA OAAA 50 - 0 2 3 4 50 0 0 8 F(DOETTED LINE COINCIDENT WITH SOLID LINE EXCEPT WHEN SHOWN SEPARATELY) VOLUME GAS CONTENTOFAQUIFERSTORAGERES RGAS BUBBLE,,\ I'~'~ ~FgPRESSURE 850 CASE I w CASE 2 800 ~ ~\ \ 2 11 \i'I\\ 0 10 20 30 40 50 60 70 SO 90 100 jNUMBER OF TWO-WEEK TIME INCREMENTS EFFECT OF GAS WITHDRAWAL DURING GROWTH ON PORE VOLUME a GAS CONTENT OF AQUIFER STORAGE RESERVOIR.n~~~~~~Fgr ~ 8

-173The results listed in Table XVII show a gas in place difference between cases 1 and 2 of 1.21 M3 cf (10.33 - 9.12 M3 cf) at the end of 180 weeks of bubble growth. However, the cumulative gas withdrawn during the three withdrawal periods prior to the 180th week was 2.72 M3 cf. Similarly, the difference between gas in place values at the end of 132 weeks is.87 M3 cf although a cumulative 1.45 M3 cf of gas was withdrawn (in case 2) prior to the 132nd week. Thus one can conclude (at least for the cases 1 and 2 considered here) that, at the end of a given gas injection season, the difference between the reservoir gas content obtained by continuous growth and the content obtained by interrupted injection of gas will be less than the cumulative gas withdrawn prior to the end of the given injection season. 2. Development of an Analytical Approach to the Moving Boundary Problem and Application to a Case Where rb is Known as a Function of Time As mentioned above, analytical solutions to the moving boundary problem, for cases where the position of the moving boundary is not known a priori as a function of time, have not been presented to date. A method is presented here which can be used in cases where rb is known as a function of time. The usefulness of this method in cases where rb is unknown a priori as a function of time is considered in Section IV-C-3 below, The dimensionless time variable in the previously presented diffusivity Equation (III-5); raJ - (II"-5)

can be redefined in a manner which will take into account any arbitrary variation of rb with time. tD, as defined previously in this thesis, is related to t by k t and dtD = k/2cprb dt. Now let rD be defined as Ar = dr /) (Iv-31) where q is the dummy variable of integration, so that kdr = d t = c/tp (IV-32) Then Equation (III-5) becomes ZP - Pp,-p (IV-33) with TD defined by Equation (IV-31). The variation of rb with time causes no difficulty with the boundary conditions Equation (111-6) or (III-7) since the conditions P(l,tD) = 1 or (P/~rD)(lltD) 1 spec the pressure and the pressure gradient at rD = r/rb = 1 and this rD value is independent of the time varying rb value. Since the form of the diffusivity equation and its boundary conditions are unchanged by substituting TD for tD, the tables of QtD and PtD vs tD, obtained by solution of the diffusivity equation for the conditions (III-6) and (III-7) and tabulated in the literature (12), can be interpreted as QTD and P vs TD and can be applied in cases where rb varies with time. The terms Q(j-i)AtD and P(j-i)AtD or QjAtD-iztD and PjatD-iAtD in Equations

-175(III-21) and (III-22) must now be interpreted as QT D and P where — ID - TDi &Q -T the dimensionless production quantity tabulated in TDj-TDi (12) the literature'' vs. tD, for tD = ~j - - and P = the dimensionless pressure drop quantity tabulated -TD T D.2 in the literature (12) vs. tD, for tD = The integral fj t dq can be evaluated approximately as XJL stJ~g _ i -it -[ tl / __/ (Iv-34) where J is obtained by setting rb equal to some assumed func0o r(q) 2 tion of time for small time. For example, if rb is assumed proportion2 al to Sjt for small time, then where A is some constant. The value of A can be determined if rb is known at ti m At as A 2 rb2(t)ime. aIf rb is assumed proportional to tx for small time then x must be less than 1 for the integral w1t dq to exist. In the present discussion, the value of rb2 at t = jLt (j = 1, 2, 3,..., N) is assumed to be known from field data. a. Application to Field B The above outlined method of treating cases where rb varies with time has been applied in determining the effective reservoir parameter values for aquifer storage reservoir B, Section III-C-2 describes the attempt to determine the field B

parameter values by calculations which employed equations derived from the assumption of a constant rb. The result of this attempt was the conclusion that no unique combination of K3 and K2 values corresponded to best agreement between the actual and predicted performance of field B (over the November 30, 1957 - July 5, 1958 period) but that this agreement continuously improved as K3 was decreased and K2 was increased. The hypothesis was made that field B's nonconformance to the constant rb assumption might be responsible for the failure to obtain a unique, optimum combination of K3 and K2 values, Thus application of the above outlined moving boundary theory should provide a test of this hypothesis. The values of rb required in Equation (III-34), above, were obtained from field B's pressure-production data in the following manner. The relationship between rb and y, where y is the vertical distance from the apex of the field B cap rock to the water level, was obtained from a contour map and is given in Equation (IV-35) below, i = 00Y (Iv-35) This relationship between rb and y is shown in the idealized vertical cross-section of reservoir B, sketched in Figure 39, Since the volume of a right circular cone is 1/3 Tr2h or, in the terms employed here, 1/3 trb2y, the field B pore volume is related to rb by the equation = Trr / (Iv-36) or r("t)= ( Vt / (Iv-37)

-177GROUND LEVEL WELLS 200__ GAS WATER L.. r e IMPERMEABLE CAP ROCK Figure 39. Relationship Between rb and y (Distance of Water Level from Apex of Cap Rock) in Aquifer Storage Reservoir B.

-178Substituting Vt from Equation (III-20) into (IV-37), one obtains 2 the following relationship between rb and nt and Pt. ( ~O + I ( 7t ) (IV-38) Equation (IV-38) now allows the calculation of rb from available pressure-production data (Pt and nt)o Insertion of rb from (IV-38) into (IV-31) yields,7) kjt)' jd/tl t ZK- it, (Iv-39) The value of TD(t) for t = At (= one week =.606 x 106 seconds) was obtained by assuming the field B pore volume directly proportional to time for small time. Thus V, t- r(k M-) = where C is some constant, or tI(k< p)- /9 t where A is some constant to be determined from field data, The data listed in Table V give = 3 opf K obt oe P6t = /07 D 7ps The value of K was obtained as -.000152 from compressibility data for the gas injected into field Bo Thus n't (k~ -A);092'/O /) /~&t - -, 060/ -2 -iaLt *,g x x/0" c2/7/x/0Y It A.rnoks/ps/, -s;- ecnds

-179and nt(K + 1/Pt) =- 171 x 10 4t for 0 < t < Ato Inserting this expression for nt(K + 1/Pt) into Equation (IV-39) and integrating from t = 0 to t A- t, one obtains __5_Dt = X x,7X / (IV-40) The values of TD = TD (j 1, 2, 3,..., N, where N was j (t jAt) taken as 31) are given by t%, __ 2 Dn,.,(Er+C/ is 1/3 +n^(g+T Z7 (IV-41)' )j + where _ / and nj and Pj are the actual (from data) pound moles gas in place and gas bubble pressure, respectively, at time t = jAt. The predicted field B pressures, for the 31 week period from November 30, 1957 to July 5, 1958, have been calculated from Equation (III-22). The nt (gas in place) quantities required in this equation were obtained from the data listed in Table V. The term P (ji)AtD or _PjAtD _itD appearing in (III-22) was interpreted as P, defined above; TD. was calculated from Equation (IV-41) D. D. and th~ prelssure and gas in place data listed in Table V. The predicted pressures Pj (j = 1,,, 3,...., 31) were calculated for various combination of K3 (= 1/2~ThkIt) and K2 values, chosen as outlined in Section III-C-1. These pressures were then employed in. calculating, for each combination of K3 and K2 values, the value of DEV, defined in Equation (III-25) (with N = 31)o DEV is listed in Table ZXIII and plotted vs K3 with K2 as a parameter in Figure 40o The predicted pressures and pore volumes, calculated

-180TABLE XVIII DEV VALUES FOR RESERVOIR B, CALCULATED BY MOVING BOUNDARY METHOD (SECTION IV-C-2) K2 = /4qpi(jTp/60ORT)2/ K3 = i/2ithkat At = 7 days =.606 x 106 seconds K2 K3 DEV (in %) 7 x 10-6 328 x 10-3 1o627 7 x 10o-6 335 x 10-3 1.593 7 x 1o~6 o341 x o10-3 1.568 7 x lo~6.348 x lo-3 1.583 7 x 10-5.162 x 10-3 1o549 7 x 10-5 o165 x 10-3 1.515 7 x 10-5,169 x 10-3 1,512 7 x 10-5.172 x 10-3 1.536 7 x 10-4.965 x 10-4 1.531 7 x 10-4.984 x 10-4 1.482 7 x 10-4.100 x 10-3 1o453 7 x 10-4.102 x 10-3 1465 7 x 10-3.676 x 10o4 1,532 7 x 10-3.689 x 10-4 1.482 7 x 103.703 x 10'4 1o472 7 x 10-3.717 x 10-4 1,480

.0170 0168 -6. 0166 X K2:7 10 0 DEV VALUES CALCULATED FROM MOVING As. 0164 BOUNDARY THEORY. _4.0162 o." --. 0160 2: 7x 1 -5 CIL 2 4' _. 0158 K 2= 7x 10.0156 o.0 152.0150 0.6 1.0 1.4 1.8 2.2 2.6 3.0 3.4 3.8 4 F 4 KS x 10 =:.......x'" DEV vs K3 for Aquifer Storage Field B. Figure 40

-182- 4 1;-4 for K3 = 1.0036 x 10- and K2 = 7 x 10 are listed in Table XIX and are compared to the actual pressures and volumes in Figures 41 and 42. Figure 40 shows that the following unique combination of K3 and K values corresponds to best agreement between the field B actual and predicted performance: I/g k3 - /,2 - i;= 7x//o sa 76o0r J (IV-42) Available field data yield h = 94' k/4 = 200 - 500 millidarcys /centipoise cp =.16 RT = 5687 psi - ft.3/# mole while a reasonable value of c, the sum of the aquifer water and formation compressibilities,3) is 7 x 10-6 vol/vol.-psia. Inserting the above values of c, Tp, and RT into (IV-42) and considering h and k/4 as unknowns, one obtains h = 92.5' k/4 = 385 millidarcys/centipoise Thus use of the developed moving boundary theory in conjunction with the method of Section III-C-1 yields unique field B parameter values, which in turn yield formation permeability and thickness values in excellent agreement with those estimated from core analysis and field data. This fact, along with the fact that the method of Section IIIC-1 yielded no unique combination of K3 and K2 values when equations derived from the assumption of a constant rb were employed, indicate

Table XIX RESERVOIR B PRESSURES AND PORE VOLUMES PREDICTED BY MOVING BOUNDARY METHOD (SECTION IV-C-2), NOVEMBER 30, 1957 - JULY 5, 1958 Reservoir Pressure, Reservoir Pore Volume, psia MMcf Date Observed Predicted* Actual Predicted* 1957 Nov. 30 106o 0 Dec. 7 1075 1078.9.059.059 Dec. 14 1075.4 1075.7.103.103 Dec. 21 1088.2 1084.2.180.181 Dec. 28 1084.9 1080o.5.238.239 1958 Jan. 4 1078.8 1070.7.251.253 Jan. 11 1093.2 1084.5.333.336 Jan. 18 1098.5 1091.4.442.446 Jan. 25 1100.7 1096.0.568.571 1958 Feb. 1 1206.9 1112.5.702.775 Feb. 8 1222.4 1151.2 1.09 1.174 Feb. 15 1103.9 1102.4 1.25 1.251 Feb. 22 1088.5 1087.3 1.27 1.272 1958 Mar. 1 1193.3 1135.1 1.50 1.598 Mar. 8 1204.6 1158.6 1.93 2.028 Mar. 15 1202.8 1167.7 2.39 2.478 Mar. 22 1166.8 1144.3 2.68 2.745 Mar. 29 1203. 6 1176.3 3.15 3. 239 1958 Apr. 5 1202.4 1182.8 3.67 3.741 Apr. 12 1191.4 1179.3 4.14 4.191 Apr. 19 1201.7 1191.3 4.67 4.717 Apr. 26 1145 1169.3 5.18 5.052 1958 May 3 1143.2 1155.0 5.37 5.303 May 10 1181.2 1181.8 5.79 5.785 May 17 1193.4 1193.2 6.32 6.322 May 24 1180.8 1186.6 6.82 6.776 May 31 1194.7 1195.4 7.30 7.297 1958 June 7 1208.8 1208.7 7.90 7.903 June 14 1207 1214.5 8.59 8.522 June 21 1190.3 1223.8 9.51 9.192 June 28 1238.9 1242.1 10.02 9.987 1958 July 5 1242.7 1250.0 10.88 10.799 * Predicted for K2 =(k/ p~pc(mcp/600RT)2/3 = 7 x 10-4 K3 = ~/2nhkAt = 1.0036 x 104 At = 7 days =.606 x 106 seconds

1280 1260 A 1240 a' 1220 1200 1160 _ 1120 _ CD 1100 1080 PPE -OA R M H 1060 1000 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 NUMBER OF 7 - DAY TIME INCREMENTS Comparison Between Field B Actual Pressures and Pressures Predicted by Moving Boundary Method, November 30, 1957 - July 5, 1958. Figure 41

IO ACTUAL PORE VOLUME 9 - -0- -0-.-0-PORE VOLUME PREDICTED BY _8 MOVING BOUNDARY METHOD w 7 7 0 w 6 0 5 W Cn I Figure 42 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 NUMBER OF ONE WEEK TIME INCREMENTS Comparison Between Field B Actual Pore Volumes and Pore Volumes Predicted by Moving Boundary Method, November 30, 1957 - July 5, 1958. Figure 42

the validity and usefulness of the moving boundary method developed. The usefulness of this method is further indicated by the fact that the minimum DEV value obtained by the calculations described in this section was 1.45% (see Figure 40) whereas the minimum DEV obtained in Section III-C-2 ('constant rb' treatment) was 1.63%. That is, the moving boundary treatment of field B presented in this section increased the agreement between predicted and actual performance 11% above the agreement obtained in Section III-C-2. Figure 41 also shows the field B pressures (listed in Table V) predicted by assuming a constant rb, as described in Section III-C-2. This figure shows the dampening effect on the rate of pressure increase effected by taking the moving boundary into account. That is, the'moving boundary pressures' are initially somewhat higher, on the average, than the'constant rb pressures' but are considerably lower after 31 weeks of bubble growth. 3. Application of Moving Boundary Method to the Case Where rb_ is Unknown as a Function of Time The analytical approach to the moving boundary problem, developed in Section IV-C-2, above, is employed here in treatment of a case where rb is unknown a priori as a function of time. Description and application is given of a trial and error method in which an assumed position of the moving boundary at time step j is employed in calculating Pj and Pj is then used to establish a more accurate position of the boundary. The gas in place schedule is considered specified and the pressures and pore volumes are considered unknown in the developments below,

-187The method developed in Section IV-C-2, above, involves replacement of the term P(j-i)LAt in Equations (III-16a) or (III-22) by P T where TD is given by Equation (IV-41)o Thus Equa-TD - Di tion (III-22) becomes (Iv-43) j ~ ~ 2 where Kt' and K7 are the same as K6 and K7 given on page 22, with and P( ) replaced by P and P, respective— Dj-DD j-i'Dij i ly. Equation (IV-43) allows an iterative procedure in which 1) a Pj value is assumed and employed in calculating TDj from (IV-41), 2) a new Pj is then calculated from (IV-43), 3) the new P. is employed to yield a new TD, from (IV-41), 4) etc., until successive values of Pj, calculated from (IV-43), are sufficiently close to one another. This iterative procedure has proven convergent in the application to field B, described below. The field B gas in place data tabulated in Tables V and XVI and K2, K, K3 values of 7 x 10-46 x respectively (the optimum combination found in Section IV-C-2, above), have been employed in calculating Pj by iterative solution of (IV-43) for j = 1,2,3,..., 52. These Pj were then employed in Equation (III-20) to yield the predicted pore volumes and in Equation (III-25) to yield DEV. A TD value of 3~78 x 105 K2, developed in Section IV-C-2, was assumed. The calculations were programmed in the FORTRAN code and carried out by an IBM 704 computer. The calculated pressures and pore volumes are listed in Table XX. Figures 43 and 44 compare these predicted quantities to the actual field performance. Figure 43 also shows the

-188TABLE XX FIELD B PRESSURES AND PORE VOLUMES PREDICTED BY MOVING BOUNDARY METHOD OF SECTION IV-C-3 Reservoir Reservoir Pore Date Pressure* Psia Volume* MMcf 1957 Nov. 30 o1060d Dec. 7 1075d 1078p.059d.059p Dec. 14 1075..4d 1075.6p.1025d.1025p Dec. 21 1088.2d 1084. lp.180d.181p Dec. 28 1084.9d 1080.5p.238d. 239P 1958 Jan. 4 1078.8d 1070.7p.251d.253P Jan. 11 1093.2d 1084.4p.333d.336p Jan. 18 1098.5d 1091.3p.442d.446p Jan. 25 1100.7d 1095.9P.568d.571p 1958 Feb. 1 1206.9d 1112.lp.702d.77~p Feb. 8 1222.4d 1150.2p 1.09d 1 18p Feb. 15 1103.9d lOl101.9p 1.25d 1.25p Feb. 22 1088.5d 1086.9p 1.27d 1.27p 1958 Mar. 1 1193.3d 1134.5p 1.50d 1.60p Mar. 8 1204.6d 1157.7P 1.93d 2.03p Mar. 15 1202.8d 1166.8p 2.39d 2.48p Mar. 22 1166 e8d 1143.6p 2.68d 2.75p Mar. 29 1203.6d 1175.5p 3.15d 3.24p 1958 Apr. 5 1202.4d 1181.9p 3.67d 3.74p Apr. 12 1191.4d 1178.5p 4.14d 4.19p Apr. 19 1201.7d 1190.3P 4.67d 4.72p Apr. 26 1145d 1168.5p 5.18d 5 o06p 1958 May. 3 1143.2d 1154.3p 5.37d 5,31p May 10 1181.2d 1181. p 5.79d 5.79p May 17 1193.4d 1192.4p 6,32d 6.33p May 24 1180.8d 1185.8p 6.82d 6.78p May 31 1194.7d 1194.6p 7030d 7.30p 1958 June -7 1208.8d 1207.9P 7.90d 7-91p June 14 1207d 1213.6p 8.59d 8.53P June 21 1190.3d 1223.0p 9.51d 9.20p June 28 1238.9d 1241,3p 10.02d 9.99p 1958 July 5 1242.7d 1249.3p 1o.88d 1O.81p July 12 1247d 1257.Op 11.76d 11.65p July 19 1242.6d 1259.2p 12.68d 12.47p July 26 1244.6d 1263.6p 13.55d 13.531p 1958 Aug. 2 1244old 1269.1p 14.52d 14.17p Aug. 9 1269.7d 1269.7p 15.42d 15.00p Aug. 16 1241 3d 1274.op 16.38d 15.86p Aug. 23 1241.7d 1279.7p 17.30d 16.76p Aug 30 1231. 4d 1274. 7P 18.30d 17o 57P

-189TABLE xX (CONTID) FIELD B PRESSURES AND PORE VOLUMES PREDICTED BY MOVING BOUNDARY METHOD OF SECTION IV-C-3 Reservoir Reservoir Pore Date Pressure* Psia Volume* MMcf 1958 Sept. 6 1050d 1290p 19.28d 18.55P Sept. 13 1258.4d 1302.7P 20.35d 19.59P Sept. 20 1258.4d 1308.0p 21.64d 20.64p Sept~ 27 1256o6d 1311e0p 22.80d 21.68p 1958 Oct. 4 1259o1d 1316.5p 24.06d 22,75P Oct. 11 1259,9d 1321,8p 25o26d 23.85p Oct. 18 1254.2d 1321.6p 26,42d 24090p Oct. 25 1249.6d 1322.6p 27.81d 25.95p 1958 Nov. 1 125io 4d 13240op 28.89d 26o99P Nov. 8 1254,9d 1327 9p 30 09d 28 o06p Nov. 15 125305d 132553P 31.06d 29.o8p Nov. 22 125503d 1329lp 32 o31d 30o14p Nov. 29 1228.6d 1316.1p 33.74d 31.03p * "d"t denotes number obtained from field data. "p" denotes predicted quantity K = 7x 10-4 K3 = 1.0036 x 10-4

-1901380 1360 1340 _ -— e DATA 0-O-O PREDICTED BY MOVING BOUNDARY METHOD AA 1320 OF SECTION IV- C - 3 A'OOO i A A A PREDICTED IN SECTION IV-C-I, A 0 1300 WITH rb ASSUMED CONSTANT A 0 1280 -AA A,o,,oo'~~ 1260AO 1240 0 1220 - 1200o g 1180 1120 1100 1080 Comparison Between Field B Observed Pressures and Pressures Predicted by Moving Boundary Method of Section IV-C-3. Figure 43

-.19140 $2 32 - - -*- CALCULATED FROM DATA oO o o 0 PREDICTED BY MOVING BOUNDAR METHOD OF 28 SECTION IV-C-3 0 24 O 4 0 i, NUMBE OF 7- DAY TIME 16;3 12 I, 0 4 8 12 16 20 24 28 32 36 40 44 48 52 j, NUMBER OF 7- DAY TIME INCREMENTS Comparison Between Field B Actual Pore Volumes and Pore Volumes Predicted by Moving Boundary Method of Section IV-c..,. Figure 44

-192predicted pressures calculated with the assumption of a constant rb, as described in Section IV-C-i1 This figure reveals that the pressures calculated from Equation (IV-43) are, in general, in considerably better agreement with the observed pressures than are those calculated with the assumption of a constant rb. As mentioned previously (Section IV-C-1), a DEV value of 3.27% was obtained for the latter calculated pressures. However, a DEV of 2.47% was calculated by inserting the pressures from Equation (IV-43) into Equation (III-25). Thus the moving boundary method described here resulted in a 24.4% (100 x (3.27 - 2.47)/ 3.27) decrease in the average difference between predicted and observed pressure. The calculations described in Section IV-C-1 yielded the field B performance predicted with the assumption of a constant gas field radius. The results (Figure 36 or 43) show that the predicted pressures steadily increased after the 31st time increment (July 5, 1958) while the observed field pressures remained relatively constant. This systematic divergence between the predicted and observed field performances was discussed at the end of Section IV-C-1. Three possible causes of this divergence were cited: 1) a growing gas bubble radius, 2) gas leakage through the reservoir cap rock, 3) vertical pressure penetration into lower formations. The results obtained in this section (see Figure 43) reveal that the divergence between predicted and actual field performances was caused partly, but not wholly, by failure to take into account the growing gas bubble radius. After the 31st time increment, the pressures predicted here by the moving boundary method increase above

-193the observed pressures in much the same manner, but not to as great an extent, as the pressures predicted with the assumption of a constant rbo The increase is too definite and systematic, however, to be considered as a random divergence about the observed performance. Thus one or both of the two other factors of gas leakage and vertical pressure penetration, which as pointed out previously would cause a systematic increase of the predicted above the observed pressures, must be held accountable for the divergence obtained. As mentioned in Section IV-C-1, data have been obtained which indicate pressure communication between the field B aquifer formation and lower formations. Since no physical evidence exists to support the possibility of gas leakage, the leveling out of the field B pressure after July 5 (31st time increment) is considered here to indicate the vertical displacement of water from the neighborhood of the gas bubble into lower formations.

V. CONCLUSIONS AND SUMMARY (1) Solutions to the diffusivity equation, given in the literature(l), have been combined with the gas law and a material balance relation to give predictive equations which allow calculation of gas storage reservoir performance. The reservoir pressure can be calculated for any time-variant gas in place schedule or the gas in place can be calculated for any time-variant pressure schedule. Several assumptions concerning the geometry and physical characteristics of the reservoir-aquifer systems were necessary in the formulation of the predictive equations. The time necessary for manual execution of the calculations, which arise in application of these equations in an actual case study, is in general so large that use of a digital computer is almost mandatory. All of the calculations described in this thesis, involving use of these equations, were programmed in the IT compiler language and carried out by the IBM 650 digital computer or were programmed in the FORTRAN compiler code and carried out by the IBM 704 digital machine. (2) The predictive equations have been employed in calculating the performance of several gas storage reservoirs. The agreement obtained between the actual and predicted performances is considered to be, with the exception of one case, more than adequate for engineering purposes; this fact indicates that the assumptions made in developing the predictive equations are in general valid for gas storage reservoirs. (3) The utility of the developed equations and methods of calculation is evidenced by the better understanding of reservoir behavior obtained by their use. The pore volume variation of a gas storage reservoir was calculated for various different idealized pressure cycles. The -194

-195results show that a pressure cycle having an equal number of pound days above and below the initial pressure will result in a slowly growing or shrinking reservoir, according to whether the storage cycle is begun with gas production or gas injection. The results also revealed that identical gas field pressure cycles can result in greatly different yearly average reservoir volumes, depending upon whether storage operations are begun upon discovery or after a period of production and, if begun upon discovery, whether with gas injection or withdrawal. (4) Consideration has been given to the problem of determining the effective reservoir parameter values. These parameters appear in the predictive equations as simple functions of reservoir physical constants such as porosity and permeability. Their'effective' values are defined as those which result in best agreement between predicted and actual field performance over a time period for which field pressure-production data are available. A computational method of determining these effective parameters is described and applied to an actual aquifer storage reservoir. The values obtained are in excellent agreement with the values estimated from laboratory and geological data. (5) The assumption of negligible vertical flow effects, made in this research, is in general valid for oil and gas reservoirs. Exceptions to this assumption have been considered in the development of an analytical solution to the diffusivity equation for the case of a storage reservoir situated on top of a very thick aquifer formation. Dimensionless flow or pressure drop quantities, similar to those presented by Van Everdingen and Hurst for the case of negligible vertical flow effects, can be obtained by numerical evaluation of the integral form of this analytical solution.

-196(6) A numerical method of solving the diffusivity equation in elliptic coordinates has been developed and employed to calculate the performance of a storage reservoir having an elliptic areal boundary. The results indicate that significant error in the calculated reservoir performance can be incurred at small and intermediate dimensionless time values (tD < 200) by approximating the elliptic boundary as an equal area circle and employing the radial flow equations. The error decreases with increasing time. This result is of importance because water influx calculations reported in the literature, as well as those described in this thesis, invariably approximate the reservoir boundary by a circle. In most cases, however, the boundary can be more accurately approximated by an ellipse. (7) A method of treating pressure interference between two or more reservoirs situated on a common aquifer is described and applied in an actual case study. The actual performances of the two reservoirs involved are duplicated closely by the performances predicted by the interference method. The interference method has been found to be considerably more accurate than the predictive method (ignoring interference effects) developed earlier in the thesis. (8) Two methods have been given here for predicting reservoir behavior when the initial aquifer pressure in a non-uniform function of the radius. The first method requires quantitative knowledge of the initial pressure distribution and of the position of the aquifer exterior boundary. This method has been applied to a reservoir-aquifer system in which the initial pressure distribution was a known function of radius.

-197The agreement obtained between the actual and calculated field performances is significantly better than that achieved by ignoring the initial pressure distribution. The second method requires knowledge of neither the initial pressure distribution nor of the position of the aquifer exterior boundary. The method does, however, require field pressure-production data over some time period. The second method again results in better agreement between the actual and predicted performances than that obtained by ignoring the initial pressure distribution. The first method is considered impractical since the initial pressure distribution and exterior radius are seldom known in actual case studies. (9) The moving boundary problem encountered in aquifer storage operations has been described and treated. An aquifer storage reservoir is created by injecting gas into a water-bearing formation originally containing no gas. Thus its bubble radius increases from zero to values of several thousand feet as the reservoir is grown. When this increasing radius is assumed constant, the method of determining the effective reservoir parameter values is shown to yield no unique combination of values corresponding to best agreement between calculated and actual performance. However, when the moving boundary method developed is employed, a unique combination is obtained which yields aquifer permeability and thickness values in excellent agreement with those estimated from laboratory and geological data. In summary, analytical methods have been developed for the treatment of the general problem of predicting gas storage reservoir performance. The investigation indicated that several special cases had

-198to be considered separately. The predictive accuracy of the developed theory has been critically examined through application in actual field studies. The results indicate that both the accuracy and applicability are more than adequate for present engineering needs. Through the theory and the applications, a substantially better understanding of gas storage reservoir performance has been obtained. The ultimate goal of ability to predict the performance of any storage reservoir, on the other hand, has not been accomplished. Theory must be developed for several additional problems and must be evaluated by application in actual case studies. The following items are recommended for future study: (a) An analytical solution to the diffusivity equation for the case of the vertical flow model has been developed and presented in the form of an infinite integral. This integral should be numerically evaluated and the results should be applied, as described in the thesis, in calculating the performance of an actual gas storage reservoir situated on top of a very thick aquifer. This predicted performance should then be compared with the actual field performance and with that predicted by the method (given in this thesis) developed from the assumption of negligible vertical flow effects. These comparisons should allow evaluation of the usefulness of the solutions obtained for the vertical flow model case. (b) The diffusivity equation in elliptic coordinates should be solved analytically for the case of an infinite aquifer.

-199Dimensionless production and/or pressure-drop quantities should be numerically determined from the analytical solution and tabulated as a function of dimensionless time. The procedure detailed in the thesis should then be employed to compare the performance of an elliptically shaped reservoir with the performance calculated by approximating the field boundary as an equivalent circle and employing the radial flow equations. If the performances are significantly different then the dimensionless production and/or pressure-drop quantities for the elliptic flow case should be calculated and tabulated for ellipses of several eccentricities in the range 0-1. (c) The effect of relative saturation on aquifer formation permeability and thus on the performance of the gas reservoir should be investigated. During the initial growth stage of an aquifer storage reservoir, the water is pushed outward into a formation essentially 100% saturated with water. Thus the permeability to water flow during this stage should be that corresponding to 100% saturation. During the first gas production cycle, however, the water moves back into a formation having a significant gas saturation. This partial water saturation will result in a permeability to water flow lower than that corresponding to 100% water saturation. This variation of permeability with time could be treated in much the same manner as the

-200time-variation of the gas bubble radius was treated in the moving boundary analysis, described in the thesis. Comparison of the observed performance of an actual gas reservoir with the performarnes predicted by ignoring and including relative saturation effects might allow conclusions concerning the importance of these effects in aquifer water flow.

APPENDIX I IBM 650 IT COMPILER PROGRAM 201

-202The TIBM 650 IT compiler program herein described effects the solution of Equation (III-21), below, as shown in the following block diagram. (IIIn21) The output or results from the program consist of the predicted reservoir (bottom-hole, psia) pressure and the predicted gas bubble pore volume at time t -At, 2At, 3tj, oe.., NAt, where N is the total number of time increments specified0 The only required input to the program is the gas storage reservoir production-injection schedule, certain physical constants representing the aquifer fluid and formation properties, and dimensionless production quantites, QtD, available from Van Everdingen and Hurst's 1949 paper(l) or Chatas' 1953 paper(2), The latter reference is recommended in preference to the former since it contains more values of o. Values of Q are tabulated vs tD D D for finite aquifers of various re/rb ratios (where re is the radius of the aquifer exterior boundary and rb is the gas field radius) as well as for'infinite' aquifers. A, Program Listing and Information Necessary for Compiling Program Language: Perlis IT compiler for IBM 650 (1-pass, Michigan - October, 1958) Packages and Suboutines Required: Compiler Package #1, square root subroutine

-203bLock 1/nGRqTtH I E RFiD STIART,, nr)n3 -3 io V0 K,, RT, /,-<'= (nRr+ T, xOt, P c~,IP - V- < IPo z).-J L > K f= k e <-'k' -,-XX k as =(ig A RT K r - I- V, - / -<K, Il t, V; r' T(k + b kL-c4q RT ) V.,j

-204Header Card Data: nI = 11, reservation for compiler'I' variables ny = 700, reservation for compiler'Y' variables nC = 365, reservation for compiler'C' variables nS = 5, number of numbered program statements nA = 50, greater than number of absolute constants assigned by program nE = 350, reservation of machine locations for package and square root extension. Program Statements: 0: READ F 0: C357 Z 6J2832 x C351 x C351 x C352 x C353 x C354 F 0: C358 Z J5 x C357 F 0: C359 Z LL C356 x Cl x C355 RS C358 x Y350 x Y1 RD C358 x Yl F 0: C360 Z 8J x C1 x C355 D C357 x Yl F 0: Y351 Z J5 x C359 SQ20EKL C360 S C359 x C359 RQ F 0: 3KIlK2KlKI0K F 1: C359 Z OJ F 0: 2KI2KOKlKIlM2K F 2: C359 Z C359S C358 x YLIlMI2R x YL351SI2RMYL349SI2R F O: C359 Z LLC356 x CI1 x C355RMC361MC359SC358 x YL348SI1R x YIRDC358 x Y1 F 0: C360 Z L8J x CI1 x C355DC357 x YIRSC359 x C359 F 3: YL350SIlR Z J5 x C359SQ20EKC360Q F 0: 5KI1LKKKIiOK F 4: YI1 Z CI1 x C355 x C356S1JDYL350SI R F 5: TYL350SIlR TYI1 TC351 TC353 TC354 F 0: H FF

-205B. Required Input or Data for Program Compiler Variable Name Physical Quantity c361 Vo, initial pore volume of gas bubble, ft.3 I10 N, an even integer, < 348; the number of time increments desired. If calculations are to be performed for a period of 5 years and a time increment of one month is specified, then N = 512= 60. C1, C2, 0.o, CN n1, n2,..o, nN, the number of pound moles of gas in place in the bubble at time t = At, 2At,.... NAt. Y350, Y349 Po, initial gas bubble and aquifer pressure, psia ~1, Y2,...., yNJ ~ atD' Q2AtD -' Yl, Y2, *o.., YN.. QAtD, Q2~tD. a.. QNZtD dimensionless production quantities corresponding to dimensionless time, tD = kAt 2kAt kNAt Suconrs b bcr C356 K, constant in equation z = 1+KP, to be determined from compressibility data for gas in question C355 RT; R, gas constant, = 10.73 psia - ft3/OR - # mole, T = average gas bubble temperature,~R C354 c, sum of aquifer fluid and aquifer rock compressibilities, 1/psia C353 h, thickness of aquifer formation, ft. C352 c, porosity of aquifer and gas sand C351 rb, radius of gas bubble, = interior radius of aquifer, ft.

-206C. Output or Results from Program Compiler Variable Name Physical Quantity Y351, Y352,..., Y(350 + N) Gas reservoir (bottom-hole) pressure, psia, at times t = At, 2At,.0.,., NAt. Y1, Y2, *0., YN Gas bubble pore volume, ft. 3 at times t = At, 2At, 3At,.o., NAt. C351, C353, C354 rb, h, and c as fed in as data

APPENDIX II SOLUTION OF DIFFUSIVITY EQUATION GOVERNING PRESSURE DISTRIBUTION IN VERTICAL PRESSURE PENETRATION FLOW MODEL -207

-208The diffusivlty equation and the corresponding boundary conditions governing the aquifer water flow in the vertical pressure penetration model appear as in Equations (1), (2), (3), and (4) below. _ - __ (1) =tay-l ~ O(2) X = - (<) 2o i < o (3) (4) where M = h/rb fkR The condition given in Equation (3) specifies an infinite aquifer; replacement of Equation (3) by a condition simulating a finite aquifer actually simplifies somewhat the process of solution of Equation (1). Taking the Hankel transform with respect to rD of both sides of Equation (1), one obtains an equation in y, tD, and parameter of the Hankel transform, q: a'b/~,y,; 5 tl~r(l,>gg~ to) r, Y~tp) (5) where (rj,y,tD) = rDP(rDy tD) Jo(rD)drD = Hankel transform of P(3,7,tD) and the Hankel transform of ( + _ j) Equation (5) can be solved by the method of separation of variables; y is assumed equal to YG where Y is a function of y only and G is a function of tD only. Equation (5) then becomes =/-, I - (6)

-209where X2 is a constant since the left side is a function of y only and the right side is a function of tD only. For X = 0, Y(y) = A cosh ry + B sinh qy; G = 1, or only other real constant and for X 0O, Y(y) = C cos I7~ 7 y + D sinJ27-y = e 2tD where A, B, C, and D are constants to be determined from the initial and boundary conditions. Setting y = YG = A cosh ry + B sinh ry + [cos y + D sin /72D-r2y] e- 2tD, one finds that But (Oy/ay)y=0 = f(qr), where f(q) is the Hankel transform of f(rD), from Equation (3); therefore Bq*Sz,-^ = faOe ~x 2-~f(U f>2xx=-x=- (7) The equalities D = 0 B = - J1(r)/2 must hold if Equation (7) is to be satisfied for all tD and r. Inserting these values of D and B into the expression for y, one obtains Osrl~yt7,)= sinA( C)o<~s -7;'eJ (8) Now, since (yy)y M = 0, where M = h/rb fkR by definition, one finds (-9 ) - i- c1 - ff /> t (9)

-210and sin NIX277M must = 0, or 1X2-2' _= mn/M = Cm, by definition, 2 2 2 and X = + n Also, from Equation (9), Inserting the values of /x72-~ and A into Equation (8), one obtains t~n,,9} C~p) = I; Z $rJofX ~~)coSh~(T)- hAny)] fC cos(odq,ye ) 9-rti s;5 n, (,Mr) A sum of terms of the form Cm cos am ye X-tD will be required in order to satisfy the initial condition, Equation (4). Thus or TCk_ 5(I3 ovr- -- or m C 05krrw) =- s/'n h Mt The terms cos Ctmy are orthogonal over the interval (O,M), where Cm = m/M, and by multiplying both sides of Equation (11) by cos Cmy and integrating over y from 0 to M, one finds:,08 5t6M y) d = -/ (,) (COSA /}(t.-)M or, performing the indicated integrations, 6k = - a:, ~k) C - J;Ct) (11) and Co = - --- (12) Thus the expression for y(T,O,tD) now becomes n s i th) (re17) - _e _ _ (13) and since the inverse Hankel transform of g(r) (where g(%) is the Hankel

-211nsform of g(r) is() Jo(rD) d, the solution for P( is -trwsorm) = iD l~ W ) coth(rDl - tD. D) C)(r on,)t (14) __ ~.~=' z/!D *1>

APPENDIX III ROUND-OFF AND TRUNCATION ERRORS ARISING IN THE NUMERICAL SOLUTION OF THE DIFFUSIVITY EQUATION IN ELLIPTIC COORDINATES -212

-2153An estimated limit on the round-off error is obtained here by determining the maximum round-off error in the asymptotic (or large time) DE and DE assuming this maximum is not exceeded in any smaller QtDE values. For example, if a maximum round-off error of 0.3% occurs in the value of QtDE at the 600th time step, then this assumption would place an upper limit of 0.5% on the roundoff error in QtDE at any time step prior to the 600th. A maximum round-off error in the asymptotic or steady-state, machine-calculated QtDE value is determined as follows. The machine-calculated QtDE (for Au =.0o8, Av =.098175, AtDE given on page 103) reaches an asymptotic value of 35.15587 as shown in Table VIIo The actual or true asymptotic QtDE can be readily calculated from a simple material balance and can be compared to the machine value of 35.15587. For at steady-state (large dimensionless time), the pressure drop attains a uniform value of 1 throughout the elliptic cylinder flow model. Thus, since the initial pressure drop was 0, the total volumetric expansion of the water in the flow model (at steady-state) must be equal to the total original volume of water multiplied by the compressibility and the pressure drop of 1. But this volume of expansion must be equal to the asymptotic or final value of QtE, the cumulative water influx into the elliptic sink, since the water can only expand into the inner sink (the exterior boundary impervious). Thus one has 1*c.fA ra(fl e — 77q h)] = Pk ez'71. or b( (7-0e e But jraebe is the area encompassed by the exterior elliptic boundary which is

-2141012.3 times the area (taibi) within the inner elliptic boundary, which is, in turn, 7r4A bt A tf6h t6 (n lh = I if rosh 07 ) ~inh( 4) -. Y 7 f7Therefore (ok ) -/O 3. 3;rr since the asymptotic machine solution for QtDE is 35.15587, the difference between the two values due to truncation and round-off error is.0048. Since the truncation error analysis developed below places a maximum error of.06 in (QDE)asymptotic (see Table IX) due to truncation error, an upper limit on the round-off error is (35.15587 +.06) - 35.151 or.065, This latter figure represents.185* of 35.15587 and an upper limit on the round-off error in D will be assumed to bel85 % for all tDE. The following truncation error analysis is based on a method described by Scarborough(22). Inclusion in Equations (III-90) and (III-91) of the lowest order Taylor series derivatives ignored in obtaining the finite 2 2 difference forms of the derivatives a P/au, aP/6v, and aP/6tDE would result in the following term added to the right side of the equations: _(4a7ld 9P~ 4ivl)'P Y 7' P__ /2 Ma; av9 L 2 tP The total truncation error in Pm,,n i is therefore the sum of the integrals of these terms; this sum can be written 6(zd)+1 B(dVia V 4 topI where A, B. and C are numbers which may be dependent upon m, n, and i. Thus if the true solution is denoted by P m,n,i and the machine solution by P then m,n, i' /Ln1z',4m,(4 L= (X") + B&ai,74 e&tie (1)

-215We now proceed to deduce the truncation error in QtDE from this form of the truncation error in Pm i' If q denotes the value of and q denotes the value of, then -=i/,+ Q f {a1I IC /V V V + (qd b Ub o Lib (2) _ 9 (()ArDL b) 6mJe C atPE where S(I) = fr( r Jv)/ and TC r ) C v If q denotes the value of as determined by the Simpsonts integration (Equation III-112), then f t f _ pA V (3) since the numerical integration introduces an error term proportional to Av. Adding Equation (2) to Equation (3), one obtains which is the expression for the difference between the true and machine values of qo The same procedure employed in obtaining q -q can be employed to obtain 4**'9 e 5 3 v );g V ((4d Vf (C e) a4 t (5) where PE= tDE and tDE = the value of q dtDE as determined from Simpsons's integration (Equation II1-114).

-216Rewriting Equation (5), and assuming the term D(L)Av to dominate the term B ()v, one obtains Qtie due= o tu)nc&atio e+()rj (6) which is an expression for the difference between the true and machine solutions due to truncation error alone.

REFERENCES 1. Van Everdingen, A. F. and Hurst, W., "The Application of the Laplace Transformation to Flow Problems in Reservoirs", Journal of Petroleum Technology, Dec., 1949. 2. Chatas, A. T., "A Practical Treatment of Nonsteady-State Flow Problems in Reservoir Systems", The Petroleum Engineer, May, 1953. 3. Katz, D. L., Tek, M. R, and Coats, K. H., "The Effect of Unsteady State Aquifer Motion on the Size of an Adjacent Gas Storage Reservoir", Paper #1114G presented at 33rd Annual Fall Meeting of AIME, October, 1958. 4. Katz, D. L., Vary, J. A. and Elenbaas, J. R., "Design of Gas Storage Fields", Paper #1059G presented at 33rd Annual AIME Fall Meeting, October, 1958. 5. Schilthuis, R. J. and Hurst, W., "Variations in Reservoir Pressure in the East Texas Field", Petroleum Trans. A.I.M.E. Vol. 114, p 164, (1935). 6, Mortada, M., "A Practical Method for Treating Oilfield Interference in Water-Drive Reservoirs", AIME PetrolQum Transactions, Vol. 204, 1955. 7. Warren, J. E, "The Unsteady-State Behavior of Linear Gas-Storage Reservoirs", The Petroleum Engineer, November, 1956. 8. Muskat, M., The Flow of Homogeneous Fluids Through Porous Media, J. W. Edwards, Inc., 1946. 9. Carslaw, H. S. and Jaeger, J, C., Conduction of Heat in Solids, Oxford at the Clarendon Press, 1947. 10. Churchill, R. V., Modern Operational Mathematics in Engineering, McGrawHill, 1944. 11. Woods, R, W. and Muskat, M., "An Analysis of Material Balance Calculations" Petroleum Trans. Vol. 160, 1945. 12. Sneddon, I. N., Fourier Transforms, McGraw-Hill, 1951. 13. Douglas, JO, Garder, A. 0. and Peaceman, D. W., "Numerical Solution of a Two-Dimensional Moving Boundary Problem", presented at the June 19-21 Houston Meeting of the Association for Computing Machinery. -217

-21814. Whittaker, E, T. and Watson, G. No, A Course of Modern Analysis, Cambridge University Press, 1952. 15. Morse, P. M. and Feshbach, H,, Methods of Theoretical Physics, McGraw-Hill, 1953. 16. Margenau, Ho and Murphy, G. Mo, The Mathematics of Physics and Chemistry, Van Nostrand, 1956. 17. Douglas, J, and Peaceman, D. W., "Numerical Solution of Two-dimensional Heat-flow Problems", AIChE Journal, December, 1955. 18. Bruce, G. H., Peaceman, D. W., Rachford, H, H. and Rice, J. D., "Calculations of Unsteady-State Gas Flow Through Porous Media", Petroleum Transactions, AIME, Vol. 198, 1953. 19, Mickley, H. S,, Sherwood, TQ K. and Reed, C. E,, Applied Mathematics in Chemical Engineering, McGraw-Hill, 1957. 20. Peaceman, D, W. and Rachford, Ho H., J. Indo Appl. Math., 3, 28, (1955). 21. Richtmyer, R. D., Difference Methods for Initial-Value Problems, Interscience Publishers, 1957. 22. Scarborough, J. B., Numerical Mathematical Analysis, The Johns Hopkins Press, 1950. 23. Hall, H. N., "Compressibility of Reservoir Rocks", Journal of Petroleum Technology, January, 1953.

UNIVERSITY OF MICHIGAN I3111111111 1[111111111511111111 3970I 3 9015 02827 3970