THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING MULTIPHASE FLOW OF. IMMISCIBLE FLUIDS IN POROUS MEDIA Robert L.Gorring June, 1962 IP-574

^ I I's v -Z

ACKNOWLEDGMENTS The author wishes to express his appreciation to the following individuals and organizations for their significant contributions to the successful completion of this dissertation. Professor Donald L. Katz whose generous donation of time to the discussion of the theoretical details was of great help in the fundamental aspects of the thesis. Professors K. H. Coats, W. P. Graebel, V. L. Streeter and M. R. Tek for their contributions of time and discussion. The Michigan Gas Association which provided financial support for the author during the period 1957-1962. The personnel of the Industry Program for their competent production of the thesis in its final form. ii

ABSTRACT The topics treated in this thesis cover three main areas: (1) Phase distribution during flow of two immiscible fluids in a porous medium; (2) Stability of macroscopic interfaces separating two immiscible phases in motion in a porous material; (3) Residual equilibrium saturation of a non-wetting phase resulting from frontal displacement of that phase by another fluid immiscible with it. The first main problem was that of formulating and verifying the basic equations of two phase flow. The principles of dimensional analysis were applied to the equations in order to reveal the basic scaling parameters. An approximate solution to the problem of linear displacement of gas by water was derived. Experimental data were taken for constant rate of water injection into the base of a vertical sand column initially containing air and water. These data agreed well with the approximate solution, A theoretical solution to the problem of the steady state capillary end effect was obtained by means of an integration of the basic equation of multiphase flow. It was found that the saturation distribution near the outflow face of a core could be predicted quite well. As an illustration of the use of the basic flow equations, the problem of countercurrent capillary imbibition was treated theoretically. The second major topic considered is that of the stability of the displacement fronts which occur when one fluid expels another from a porous material. The rigorous solution of the problem in which the iii

two fluids are moving at an initially constant velocity, and are separated by a discrete interface, was obtained. The analogous problem, in which the fluids are accelerated from rest, was also solved, The effect of viscous damping on the propagation of instabilities was treated fundamentally by means of rigorous formulations of the work-energy and Bernoulli equations for porous media. The effect of capillary pressure gradients on the stability of quasi-steady flows was treated from a theoretical viewpoint. The effect of capillary imbibition was determined quantitatively insofar as it tends to stabilize the flow by destroying the discrete nature of the interfaces. Helmholtz flows in adjacent layers of porous media were also treated with respect to interface stability. The problem of the stability of radial, quasi-steady flow was solved for the cases of thick and thin storage reservoirs. Finally, the relation of the influence of model boundaries on instability propagation was examined theoretically. The third main topic of this work is that of the magnitude of the non-wetting phase residual equilibrium saturation. The experiments performed showed that the microscopic fluid-fluid interface moves across a given pore unstably and that a pore is either completely full of wetting phase fluid or else almost completely empty. It was also shown that the residual is due to a blocking process in which a zone of non-wetting phase becomes isolatedo An empirical correlation of available data yielded the conclusion that the non-wetting phase residual saturation is a function of porosity oily, provided that the initial wetting phase saturation is sufficiently low. Experimental work also showed that the residual gas saturation after water flooding is independent of the flooding rate over a wide range of velocities. iv

TABLE OF CONTENTS Page ACKNOWLEDGMENTS.................................................'.ii ABSTRACT................. o................ eeiii LIST OF APPENDICES... e...........,..<..........vii LIST OF TABLES................................... viii LIST OF FIGURES........... e o ** *-* e**@**ix NOMENCLATURE...................................................... xiii I. INTRODUCTION...............o...................D........... 1 A. Saturation Distributions During Immiscible Displacement... 1 B. Stability of Interfaces................................. 5 C. Residual Equilibrium Saturation........ o.............. 6 II. SATURATION DISTRIBUTION DURING MULTIPHASE FLOW IN POROUS MEDIA..............o.............o o.. *.....oo*...... 9 A. Formulation and Dimensional Analysis of Basic Equations................................................. 9 B. Saturation Distribution During Constant Rate..... Displacement.............................................. 2 C. Steady State Capillary End Effect......................... 6 D. Countercurrent Capillary Imbibition...................... 76 III. MACROSCOPIC STABILITY OF INTERFACES DURING DISPLACEMENT PROCESSES IN POROUS MEDIA.................................... 93 A. Stability of Plane Interfaces in Linear, Quasi-Steady Flow-Effects of Inertial Damping......................... 93 B. Viscous Damping of Unstable Fingers....................... 124 C. Effect of Capillary Pressure Gradients on the Hydrodynamic Stability of Quasi-Steady Flows........................... 1 D. Stability of Plane Interfaces in Accelerated Flow......... 148 E. Stability of Stratified Flows in Porous Media............. 151 F. Stability of Radial, Quasi-Steady Flows, with Application to Thick and Thin Gas Storage Reservoirs............ 162 G. Effect of Capillary Imbibition Between Fingers on the Stability of Quasi-Steady Flows...................... 173 H. Effect of Model Depth on Perturbation Growth............. 18 V!

TABLE OF CONTENTS (Continued) Page IV. RESIDUAL EQUILIBRIUM SATURATION IN POROUS MEDIAo ooooo 188 A, Description of Microscopic Flow Processes in Porous Media,, e ee..o.. oo.................. o..... o.......... 191 Bo Theoretical Aspects of Microscopic Motion................. 205 C. Prediction of Residual Equilibrium Saturation............ 207 BIBLIOGRAPHY.....e.oo oooooooouooooooooe 225 APPENDICES,.. o........ o.... o.o...o........... 233 Ao Analytical Determination of the Derivative S....o... 234 zt 2 Bo Residual Equilibrium Saturation - Tabular Data from Literature................................................ 237 C. Mounting of Consolidated Coreso........ o............. 238 D. Mounting of Unconsolidated Materials.................... 242 Eo Residual Equilibrium Saturation - Tabular Data for Figures 43 and 44.....o...........................2 244 vi

LIST OF APPENDICES Page APPENDIX as 23~ A. ANALYTICAL DETERMINATION OF THE DERIVATIVE........... 2 z B. RESIDUAL EQUILIBRIUM SATURATION - TABULAR DATA FROM LITERATURE................................7................. 237 C. MOUNTING OF CONSOLIDATED CORES............................. 238 D. MOUNTING OF UNCONSOLIDATED MATERIALS......................... 242 E. RESIDUAL EQUILIBRIUM SATURATION - TABULAR DATA FOR FIGURES 43 AND 44............................................ vii

LIST OF TABLES TABLE Page I PROPERTIES OF SAND AND FLUIDS USED IN FRONTAL DISPLACEMENTS............................................. 3 II PROPERTIES OF MATERIALS USED BY DOMBROWSKI (31)........... 67 III CHARACTERISTICS OF ASSUMED WATER-BLOCKED GAS FIELDS....... 73 IV CRITICAL VELOCITIES FOR UNSTABLE DISPLACEMENT IN POROUS MEDIA.............................................. 181 V DATA TAKEN FROM LITERATURE FOR RESIDUAL EQUILIBRIUM SATURATION-WETTING PHASE DISPLACING NON-WETTING PHASE..... 237 VI PROPERTIES OF POROUS MATERIALS USED IN RATE EXPERIMENTS... 214 VII COMPARISON OF RESIDUALS COMPUTED BY WEIGHT DIFFERENCE AND PRESSURIZATION.................................. 217 VIII WEIGHT OF CORES AFTER WATER BREAKTHROUGH................. 218 viii

LIST OF FIGURES Figure Page l(a) Saturation Front During Displacement of A by B from a Porous Medium Initially Saturated with A and Irreducible B.................................... 4 l(b) Capillary Pressure Curve.................................. 4 l(c) Relative Permeability Experiment..................... 4 2 Development of Unstable Fingers in a Porous Medium...... 7 3 Relative Permeability vs. Saturation.......1............. 11 4 Leverett Function vs. Wetting Phase Saturation - Unconsolidated Sands..1................................ 14 5 Surface in Space Generated When Saturation S' is Plotted as a Function of Distance and Time S' = f(z,t)...,.,,.... 236 6 Saturation Distribution During Displacement of Air by Water u = 0.00490 ft3............................... 37 ft2hr 7 Saturation Distribution During Displacement of Air by Water u =.0126 ft3............................. 38 ft 2hr 8 Saturation Distribution During Displacement of Air by Water u - 0.441 ft3..............,.......... 39 ft2hr 9 Saturation Distribution During Displacement of Air by Water u = 0734 ft....................... 40 ft hr 10 Saturation Distribution during Displacement of Air by Water u - e157 ft3...................... 41 ft2hr 11 Saturation Distribution During Displacement of Air by Water u =.157 ft3.........*................ s*e*S 42 ft2hr 12 Saturation Distribution During Displacement of Air by Water u =.157 ft3.................4.............. ft2hr ix

LIST OF FIGURES (Continued) Figure Page 13 Saturation Distribution During Displacement of Air by Water u = 2239 ft3 0. o 0 0 e o....... e...a o 44 ft2hr 14 Saturation Distribution During Displacement of Air by Water u =.285 ft3.. 0 0 0 0.. o o o o....... 45 ft2hr 15 Relative Permeability vs. Saturation for 150-200 Mesh Sand in Low Saturation Region o.o...o.............46 16 Saturation Distribution During Displacement of Natural Gas by Water; u = o0216 ft3.......................oo 48 ft hr 17 Saturation Distribution During Displacement of Natural Gas by Water - Capillary Imbibition...................... 49 18 Saturation Gradients During Displacement of Air by Water (t = 1.00 hr in Figure 10)..................58 19 Saturation Gradients During Displacement of Air by Water (t = 2 00 hrs. in Figure 10).......5............. 20 Liquid Saturation Adjacent to Outflow Face of Porous Bed - 100/115 Mesh Glass............................. o 69 21 Liquid Saturation Adjacent to Outflow Face of Porous Bed 60/65 Mesh Glass o...................., 70 22 Liquid Saturation Adjacent to Outflow Face of Porous Bed - 35/42 Mesh Glass.............................. 71 23 Productivity Ratios vs. Reservoir Radius for Various Cases of Water Blocking............................ 75 24 Height Above Base Plane vso Fractional Water Saturation at Various Times for Typical Gravity Drainage Process....... o....... 78 25 Fractional Wetting Phase Saturation vso Distance During Typical Countercurrent Capillary Imbibition....... 81 26(a) Distribution of Water and Gas During Linear Countercurrent Capillary Imbibition Into Porous Medium.......... 89 x

LIST OF FIGURES (Continued) Figure Page 26(b) Rate of Production of Gas per Square Foot of Producing 90 Face During Linear Countercurrent Capillary Imbibition Into Porous Medium........................................ 26(c) Total Gas Production per Square Foot of Producing Face During Linear Countercurrent Capillary Imbibition Into Porous Medium....................................... 91 27 Schematic Diagram of Linear Displacement of Fluid 2 by Fluid 1 at Start of Fingering Process................. 97 28 Differential Element of Porous Media..................... 98 29 Experiment to Determine the Pseudo-Viscosity p' of a Porous Medium Saturated with Liquid of Viscosity......... 117 30 Stratified Flow of Fluids in Porous Media................. 1 31(a) Argand Diagram Showing Components of Instability Index n.. 157 31(b) Zones of Stability for Stratified Flow in Porous Media.... 157 32 Top and Side View of Disc-Shaped Reservoir of Finite Height. Radial and Longitudinal Perturbations are illustrated....................................... 164 33 Single Gas Finger Penetrating Water-Saturated Porous Formation............................................ 175 34 Schematic Diagram of Linear Displacement of Fluid 2 by Fluid 1 in the Case of Finite Boundaries................. 185 55 Capillary Condensation in Thin Tubes...................... 194 36 Schematic Diagram Illustrating the Formation of Residual Saturation by the Blocking Process...................... 194 37 Micromodel Packed with 6 mm. dia. Glass Beads-Simple Cubic Packing.................................... 202 38 Invading Interface and Pendular Rings During Displacement of Gas by Water.................................... 203 39 Residual Equilibrium Non-Wetting Phase Saturation vs. Porosity for Unconsolidated Porous Media and Various Sandstones............................. 211 xi:

LIST OF FIGURES (Continued) Figure Page 40a Schematic Diagram of Mounted Consolidated Core and Preparation Procedureso.. o...........ooooooooooooooo 240 40b Schematic Diagram of Mounted Consolidated Core and Preparation Procedures.......o....................o oo 241 41 Schematic Diagram of Mounted Unconsolidated Core....o...... 243 42 Ruska Pump with Mounted Core. o.......................... 215 43 Effect of Flooding Rate on Non-Wetting Phase Residual Equilibrium Saturationo. o oo... o... o o oo... o.eoo o oo oo 219 44 Effect of Flooding Rate on Non-Wetting Phase Residual Equilibrium Saturation o ooo......................... o.oo.o 220 45 Effect of Initial Water Content on Residual Oil After Waterflooding................o....... o.................... 221 xii

NOMENCLATURE a radius reservoir a constant defined in Equation (142) a,b,c,d defined by Equation (166) and Equation (167) a1 rate of acceleration in Section III D a2 defined in Equation (25) A constant defined in Equations (95) ff A negative saturation gradient at z b perturbation imposed on displacement front, cf. (Equation (94) and (181) ) B amplitude of perturbation C1cc2 principal interfacial curvatures Cn constants Dn constants in Equation (28) D hydraulic diffusivity for countercurrent imbibition (cf. Equation (63) ) f denotes "function of..." F1 = 1 + rgw] -1 K [I rw g FD = frictional force exerted by pore walls on flowing fluid Fr 7 Froude number for stratified flow in porous media gi denotes "function of..." g acceleration of gravity I modified Bessel function of first kind i -l *Note: Nomenclature conforms to AIME Code (Trans. AIME. (Petr.), 207, (1956). xiii

NOMENCLATURE (Continued) J Leverett J function, cf. Equation (4) K modified Bessel function of second kind k wave number of perturbation = X k* reduced wave number, defined by Equation (178 b) k angular wave number, defined following Equation (180) L characteristic dimension of porous material length of sand pack & half width of viscous finger shown in Figure 33 il',2 cf. Figure 34 m wave frequency, defined following Equation (180) n stability index (cf. Equation (94)) or normal derivative designation p pressure of total flow p pressure of undisturbed flow p' pressure of disturbance only Pi. pressure in phase i when capillary forces are present c Pc capillary pressure q,2 square of perturbation velocity, only, in porous media -2 q square of total flow velocity in porous media &Q flow rate of unblocked reservoir &QB- flow rate of water blocked reservoir r radius r1, r2 principal radii of curvature rB well bore radius xiv

NOMENCLATURE (Continued) rb field extorior radius Rt capillary tube radius S' - S rw S normalized saturation, ~ or surface element Si - St rg rw S' saturation (absolute) S1 breakthrough saturation, cf. p. 167 of (107) S' residual gas saturation* gr S' residual water saturation wr S' inlet face saturation during countercurrent capillary imbibi0 tion S' - S' - S' gr wr S - SO - wr t time ui velocity of undisturbed flow in i direction u' perturbation velocity (only) in i direction u microscopic point velocity during flow in porous media u, v, w total volumetric flow velocities in x, y and z directions u, v', w' perturbation volumetric velocities in x, y and z or r, Q and z directions U, V, W normalized volumetric velocities in principal directions V uniform velocity of undisturbed flows W, volume element Vc critical velocity for instability Vc critical velocity for instability considering capillary annihilation S' is, in this work, taken to mean the"ater saturation which exists when gr the porous medium is at residual equilibrium gas saturation. xv

NOMENCLATURE (Continued) x,y,z principal directions in Cartesian coordinates X,Y,Z principal normalized directions in Cartesian coordinates X,Y,Z body force components in principal directions z. distance from inflow face of porous medium or height above datum plane Zo depth of fluid penetration in displacement experiment, cf. Figure 1 (a) xvi

NOMENCLATURE (Continued) GREEK SYMBOLS cx fraction of displaced phase flooded out (cf. Figure 33) a defined in Equation (107) ff _ defined in Equation (136) ff 5x. differential length element in i direction' desaturating potential, cf. Equation (44) ff Gab',cd phase angles defined in Equation (168) ff 0 contact angle or angular direction in Sec. III F 0 wave angle ~K permeability (1 Darcy = 1.06 x 10-11 ft2) Kh hydraulic conductivity for countercurrent imbibition (cf. Equation (63) ) X wavelength of perturbation M. viscosity v kinematic viscosity (p ):rho normalized distance ( ) zo it 53.1416 p density aii normal stress component on i face in i direction a surface tension; o* effective bulk interfacial tension in porous media ~T ~ shear stress 0 porosity (D velocity potential of total flow in porous medium nQ ~ body force potential (equals + g z in gravitational field) Xo vorticity X source strength xvii

NOMENCLATURE (Continued) SUBSCRIPTS 1 phase 1 (usually displacing phase) 2 phase 2 (usually displaced phase) g gas aJ ~ imaginary component i summation index m maximum or naturally occurring o denotes value of quantity at inflow face (z = o) R real component r relative t total U undisturbed w water z denotes value of quantity at leading edge of saturation front (z Zo) xviii

I. INTRODUCTION The object of the engineering research undertaken in the area of multiphase flow in porous media is to predict the outcome of complex flow experiments on the basis of simply measured basic properties of the porous media and the fluids involved. It is necessary that the properties so measured be representative of the system to which the calculations are applied and that these measurements lead to expressions which can be used to formulate the problem mathematically. The research presented here is designed to elucidate the phenomena which occur when one fluid displaces another which is immiscible with it from a porous material. These processes present a formidable degree of complexity even under the best conditions. The number of possible situations that can arise is large. Consequently, an approach which purports to achieve any reasonable degree of generality must necessarily resolve itself into an attack on certain separate problems. A. Saturation Distributions During Immiscible Displacement The determination of the fluid phase distribution as a function of time and distance constitutes the solution of problems in two-phase flow in porous media. When one fluid displaces another which is immiscible with it from a porous matrix, the fraction of the pore space filled with the displacing fluid will vary with time and distance. This fraction is called the saturation, and must be understood in a macroscopic sense, i.e. as the ratio of the volume of displacing fluid to the total pore volume when averaged over a large number of pores. Commercially important examples of this displacement phenomenon are: (1) recovery of underground -1

-2oil by sweeping with water, (2) movement of natural gas and water in a storage reservoir, (3) simultaneous movement of two immiscible phases in packed absorption and distillation towers, (4) dewatering of filter cakes in chemical operations. As the displacing fluid moves into the saturated porous body, an interface between the immiscible fluids is formed. This interface (or front) is visualized as the distance in which the displacing phase saturation changes from its maximum to its minimum value. This saturation front is illustrated in Figure l(a). Ahead of the front, the pores will be occupied by essentially all displaced phase fluid (A), while far behind the front, the pores are nearly entirely filled with displacing phase (B). These terminal saturation values are not 0 and 1.00, however (saturation defined as fraction of pore space filled with whichever fluid is chosen as the reference fluid, in this case fluid B). In this work the displacing phase, B, is usually water, the displaced phase A being gas. The minimum permissible wetting phase saturation is usually called irreducible minimum (water) saturation, S' and corresponds to the connate water wr' saturation found in reservoirs. It represents the least fraction of the void space that can be filled with wetting phase. The minimum irreducible non-wetting phase saturation is called residual equilibrium (gas) saturation, S', and represents the maximum wetting phase saturation that can be gr~ attained by frontal displacement. The wetting property is defined as the ability of a fluid phase to be spontaneously imbibed into the porous medium. Since the porous media under consideration often are cores taken from wells, the term "core," is used interchangeably with experimental laboratory porous media. A core into which water is spontaneously drawn with resultant expulsion of, say, oil is termed water wet. It would be

-3desirable, ideally, to be able to predict the saturation distribution curve of Figure l(a) from the viscosities, densities and surface tension of the fluids and the properties of the porous media, such as porosity and permeability. However, at present, data on residual equilibrium gas saturation, irreducible minimum water saturation, capillary pressure, and relative permeability are needed for the purpose of formulating the problem. The capillary pressure concept is illustrated in Figure l(b) and may bethought of as the height corresponding to a given liquid saturation when the liquid has been allowed to imbibe vertically to equilibrium. The relative permeability concept is illustrated in Figure 1(c). Fluids A and B are entering and leaving the core at steady state and generating a pressure drop. The pressure drop depends on the saturation in the core at a given time. The relative permeability thus defines the influence of one phase on the flow of the other, and is equal to the ratio of phase permeability at saturation S' to that at S' = 1.00. When a non-wetting phase is displacing a wetting phase, and the system has reached steady state, i.e. the wetting phase fluid in the core is not moving, there is always observed to be an increase in the wetting phase saturation near the outflow face of the core. This saturation increase is called the capillary end effect, and is due to the retention of the water at the outflow caused by the saturation discontinuity at that point. In dewatering of filter cakes, and in performing laboratory experiments on petroleum cores, it is desirable to minimize this effect, or at least to have some idea of its magnitude.

-4Sgr B \ PHASE A A ENTERING B PHASE B LEAVING zo Z -* DISTANCE FROM INFLOW FACE Figure l(a). Saturation Front During Displacement of A by B from a Porous Medium Initially Saturated with A and Irreducible B. 1-~ AP HEIGHT A i —-- A ABOVE 4 I I BASE B,B Figure l(c). Relative Permeability Experiment. 0 S (frc) 1.00 5 B ) Figure l(b). Capillary Pressure Curve.

-5An important problem in reservoir mechanics is that of countercurrent capillary imbibition into a porous formation, In general, a wetting phase will spontaneously imbibe into a porous body. If an experiment is set up in which one end of a core is open and the other closed, the initially contained non-wetting phase must be expelled countercurrent to the imbibing phase. It frequently occurs thatthis is the most important means of recovering petroleum from reservoirs in which the oil is contained on long horizontal streaks which cannot be produced by frontal displacement. In many practical applications, it is desirable to know the saturation space-time history of such a system. B. Stability of Interfaces The problem of stability in two phase displacement problems is the determination of whether small disturbances on a (macroscopic, discrete) interface separating the fluids will grow or decay with time. In these stability problems, the macroscopic interface is assumed to be sharp, i.e. without significant saturation gradients on either side. It is important, at this point, to distinguish between the macroscopic and microscopic approaches. The microscopic approach involves examination of the actual process going on inside a given pore. All multiphase flows in porous media have their microscopic aspects irrespective of whether the gross flow is macroscopically stable or not. The macroscopic approach involves the examination and averaging of many hundreds of pores, and results in mean quantities. The problem of stability is considered as a macroscopic problem insofar as it is treated in this thesis. The growth of the perturbations in the macroscopic unstable case results in the formation of long fingers of displacing fluid moving at relatively high

-6velocity through the formation and bypassing large quantities of the inplace fluid. Photographs of fingers in porous media are shown in Figure 2. The light areas in Figure 2 are water fingers moving into a glass powder pack which was initially saturated with oil of the same refractive index as the glass. The simplest problem in stability analysis is that arising when two immiscible fluids are found flowing at constant velocity and separated by a plane interface on either side of which the saturations are constant. Small perturbations are then assumed to appear on this surface, giving it a corrugated shape. Growth of these perturbations results in finger formation. Decay results in steady, frontal displacement, as previously described. It is desirable to find a relation between such parameters as permeability, porosity density, surface tension, etc., and the tendency of the porous medium to form unstable fingers. Energy relations must also be satisfied in the growth of fingers, so that the effect of viscous damping on the kinetics of growth must also be considered. Finally, if fingers of wetting phase fluid are formed, capillary imbibition must take place between the fingers due to the sharp saturation gradients generated. This imbibition will tend to annihilate the fingers by filling the space between them. These problems, together with others which are related to stability, are solved in this thesis. C. Residual Equilibrium Saturation When a wetting phase moves into a porous formation containing a non-wetting phase, a certain amount of the non-wetting phase is left behind and is non-recoverable by the frontal displacement process. The problem to to be solved here is that of the determination of the final saturation of

-7N ~-W ^-. ^ N 3 13%- W. 2- 3% p 1 S _ H i N: W. 9. 5 N - 34%; W. = Io0% I N = Wi. 12% N 52%; W. 650% Figure 2. Development of Unstable Fingers when Water Displaces Oil from a Porous Medium (from van Meurs, Trans. AIME, (Petr. ), 213, 103, (1958). Np = oil production in pore volumes W. = cLmulative water injection in pore volumes. )

-8the non-wetting phase under equilibrium conditions, i.e., the quantity S' in Figure l(a). For example, when water displaces gas from a water gr drive reservoir or in aquifer gas storage projects, this phenomenon occurs. The gas is physically trapped in the interior of the pores by the advancing water. The ideal case would be to be able to predict the residual gas saturation for natural porous media from measurements of the viscosities, densities and surface tension of the fluids and porosity, permeability and other simple properties of the rock. All the available evidence, both empirical and theoretical, shows that, in the microscopic sense, the wetting phase can remain in a stable position only up to a certain critical point as it enters the pore. The interface then jumps through the pore to a new position of stability by bypassing some of the non-wetting phase. The magnitude of this amount of bypassed fluid is of prime importance in the recovery of gas and oil. The quantitative description of this unstable interface jump appears to be the key to an understanding of residual equilibrium saturation. Experimental values of residual equilibrium saturation on naturally occurring rock formations are useful for field applications, but give little insight from a fundamental point of view because of the unspecified pore geometry. Consequently theoretical and experimental consideration of the microscopic phenomena occurring in ideal packings of spheres is necessary. In particular, the effects of displacement rate, surface tension and contact angle have to be reduced to a quantitative basis if an adequate basic understanding of residual equilibrium saturation is to be achieved.

II. SATURATION DISTRIBUTION DURING MULTIPHASE FLOW IN POROUS MEDIA The four subsections of Section II include a formulation of the basic equations followed by the application fo these equations to the problems of frontal displacement, capillary end-effect, and countercurrent capillary imbibition. The general procedure in each of the subsections will be to include a literature survey, theoretical development, experimental results, if obtained, and a final discussion. A. Formulation and Dimensional Analysis of Basic Equations This section describes the development of the mathematical statement of the two-phase flow problem starting from the simple concepts of Darcy's Law, capillary pressure and relative permeability. 1. Single Phase Flow in Porous Media The general theory of single phase flow in porous materials is highly developed, and may be found in extended treatments (85), (108). General background on the problems of the oil and gas industries may be found in (61) and (95). An excellent general introduction to the theoretical aspects of flow in porous media is given in (15), which must be consulted for the most elementary background. 2. Relative Permeability The fundamental concept in two-phase flow in porous media is that of relative permeability. It appears to have been first enunciated by Wyckoff and Botset (137). The relative permeability to a given phase is defined to be the ratio of the permeability at the specified saturation to the permeability at 100% saturation of the phase in question. -9

-10Physically, it is thus a measure of the reduction of the permeability to one phase due to the presence of the other. The quantitative expression for the permeability is given by Darcy's Law ^^-*(Vp~pV^) (1) u = - ( p + p V7 ) The dynamic relative permeability equations for two phase flow, taking the saturation into account, are written for gas and water by means of a natural extension of Darcy's Law. u = - (Vg+ pgV) (2a) UW = - (Vp + pVQ) (2b) Flw The above equations are simply speculations about the nature of two phase flow, and serve as a convenient method of correlating experimental data. An example of relative permeability curves is given in Figure 3. The original experiments in (137) demonstrated that relative permeabilities were functions of saturation only, to a good first order of approximation. Additional evidence bearing on this point is given in (107). Also presented (107) is a thorough discussion of the various experimental methods used in relative permeability measurements. A recent experimental investigation by Sandberg et. al. (106) provides experimental evidence that relative permeability is, to a first approximation, a function of saturation only. When the saturation distributions due to capillary end effects were properly taken into account, it was shown that, within limits, the relative permeability is not a function of flow rate or viscosities. In the absence of conclusive proof, and for

-111.0, 1.0.9 - WYCKOFF BOTSET (ALL UNCONSOLIDATED SANDS) 0 150-200 MESH SAND.8 ^ - ~ ~~ ^ ~ ~~~~ ~^ - ~.6.6 0\/ wr ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ a..3.3 0.2.3.4.5.6.7.8.9 1.0 WATER SATURATION, Sw Figure 5. Relative Permeability versus Saturation.

-12purposes of analytical formulations, it will be assumed in this thesis that relative permeability is a function of saturation only. Additional comments regarding the uniqueness of the relative permeability functions are given by Scheidegger (107), p. 156, and Richardson (102). The latest development in techniques for measuring relative permeability was given by Johnson et.al. (56). It is called the relative injectivity method, and involves measuring the pressure drop and production history across a core during a simple frontal displacement. The method is theoretically rigorous and gives individual relative permeabilities rather than their ratio. The method is far simpler than the others, and gives results in agreement with them. The data shown in Figure 3 were taken by this method and are seen to be in agreement with those taken by Wyckoff and Botset by the steady state method. Further work on relative permeability is given in (26), (39), (43), (96) and(102). 3. Capillary Pressure Capillary pressure is the second basic concept in multiphase flow theory. The primary work in this field is by Leverett (73). Capillary pressure was defined to be the difference between the pressure in the nonwetting phase and that in the wetting phase P = Pg - Pw (3) Leverett was also able to correlate the capillary pressure curves for a wide variety of unconsolidated porous materials by means of a single function. J(S') = () f1/2 Pc or(~) ~~~~~~~~~(4)

-13where J(S') is a single valued function of the wetting phase saturation. The J function is slightly different when the wetting phase saturation is increasing (imbibition) than when it is decreasing (drainage). The static capillary pressure can also be interpreted as the height above a given datum to which the wetting phase fluid will spontaneously rise to produce a given saturation. Experimental work which establishes an important conclusion regarding two phase flow was done by Brown (9). In particular, he showed that the capillary pressure vs. saturation curve was the same whether this pressure is measured with both phases moving or with both phases stationary. This means that the relation between capillary pressure and saturation will not change in two phase flow situations in which the saturations are changing. Data are presented which show that each lithologic type of consolidated rock will possess a nearly unique J function describing the capillary pressure. A diagram of a capillary pressure curve is given in Figure 4, together with the average drainageimbibition curves of Leverett (73). Additional treatments of capillary pressure are given in (5), (12), (48), (58), (122) and (139). 4. Wettability Wettability is the third basic concept necessary to an understanding of multiphase flow. The classical approach to wettability is to define it as being some function of the three phase contact angle G. This is the angle, measured through the wetting phase, which the interface between the two fluids will assume at the point of contact with the solid. Wagner (127) elucidates this point of view, and gives a discussion of it, together with diagrams. The quantitative definition of wettability, according to this point of view (129), is that it is numerically equal to

-141.0 - 1 0.9 08 DATA OF LEVERETT (7), Avg. X 150-200 MESH SAND 0.7 0 03 X ~~ I^ I x \ -BIO O. OZlOb 0.6 uF 4. Lv~-t Fucto vDRAINAGE 03 02 0 01 0.2 0. 004 0X5 6 0.7 0.8 09 1.0 1.1 1.2.3 1.4 d X X ~l X / IMBIBITION 0 OJ 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 Sc, WETTING PHASE SATURATION Figure 4. Leverett Function vs. Wetting Phase Saturation - Unconsolidated Sands.

-15the cosine of 0. Denek&s et.al. (30) illustrate another concept of wettability. Relative wettability is here (30) to be given as the ratio of the initial rate of capillary imbibition of the given fluid to the rate at which water is initially imbibed into the same core. A proposal which purports to give a quantitative interpretation of wettability has been outlined by Calhoun (16). The suggestion is that the function of 9 (f(Q)L) of Equation (4) be replaced by cos 0 for the purpose of showing quantitatively the effect of contact angle on capillary, pressure. The last, and by far the most basic treatment of the effect of contact angle on capillary pressure has been given by Melrose (80). He exposes the u.se of cos 0 in Equation (4) as an unjustifiable attempt to avoid- consideration of the contact angle as a boundary condition dependent only on the three phase surface free energies of the gas, liquid and solid phases. Using the concept of the ideal soil, Melrose derives expressions for the drainage and imbibition capillary pressures, which are in good agreement with experimental data. He then goes on to show that the effect of contact angle cannot be included in the capillary pressure function by so simple a procedure as multiplying by cos 0. ItI-was also demonstrated that the microscopic approach to interface movement yields an equation which must be solved numerically. Additional information may be found in (2), (39), (63) and (81). 5. Formulation of Basic Equations The basic equations of two phase flow in porous media were first formulated by Muskat (85), p. 302, without the capillary pressure term. The complete formulation.was achieved by Rapoport (97) and (98), who included frictional, gravitational, and capillaric effects.* Rapoport (98) also showed how the boundary conditions could be properly specified. *The fractional flow equation (Equation (9)) was formulated by Leverett (73).

-166. Scaling Scaling of multiphase flow experiments has often been used in lieu of adequate analytic or numerical solutions to the basic equations. The theory of scaling has been developed for the purpose of aiding in the application of laboratory experimental data to field operating conditions. The most extensive treatment of scaling of two phase flow experiments has been given by Geertsma, et.al. (41). Problems in cold and hot water displacement of oil, and miscible flooding are treated in detail. Important similarity groups are derived for each case and the conditions under which scaling can be achieved, are discussed. It is shown that exact scaling is very difficult. Other treatments of scaling are given in (27), (70), (74), (90), (97), and (98). The problem to be treated is that of obtaining a mathematical statement of the multiphase flow phenomenon. The procedure will be to use the quantitative statements of the foundations of the theory as outlined previously to formulate the governing differential equation and boundary conditions. This formulation will then be used as the basis for a dimensional analysis. 7. Theoretical Development. The general problem of the movement of two immiscible phases in a porous matrix may be stated by writing the following system of equations: Subscript g indicates gas phase, subscript w indicates liquid phase. Equations of Motion; c.f. Equations (2).

-17Equations of Continuity: - V. (Pg) - ( S a(Ps~) Equations of State: Equations of State: Pg = g(Pg) (6) Pw = g(pw) The general solution of the two phase flow problem requires the simultaneous solution of this system of equations for arbitrary boundary conditions. However, certain simplifying assumptions, as applied to the above stated system of equations, will yield a single partial differential equation. Equation (3) is assumed to apply to the system. Under the assumption that both phases are incompressible, it is noted that the total flow ut equals the sum of the gas and liquid flows. ut = Uw + g (7) This latter assumption will hold very well for liquids. For a gasliquid system, the assumption is still good provided the pressure drops are not a significant fraction of the total pressure. An example of this would be a natural gas storage reservoir in which the pressure drawdowns do not exceed 25%,/ of the ambient pressure. The gas is thus behaving as if incompressible. In line with the evidence presented previously, relative permeabilities are considered to be single valued functions of saturation,. Viscosity is assumed constant.

-18Proceeding to the formulation of the equations, we obtain from Equation (2) -w Uw = - 7w - Pwv2 (8a) k g=- Vp - p VQ (8b) K Krg g g subtracting these equations, substituting for u from Equation (7), and g using Equation (3), we calculate uw as ut + ziK rgb/l [VP -A p (g (9) 1+ Krgkw/ Krw kg For the case of both phases incompressible, the Equations (5) of continuity reduce to the single equation: d s -V () (10) 6t W Substituting Equation (9) in Equation (10), and expanding; ~+ u +V -- +w. ~ LL- " X dS' u _d+s +V as- s +w + a ( -t dS L ax ay z g x r dS (11) ay ( dPc S a ( rgF dP cS ) g d (KrgF1)a = z which is the fundamental equation to be solved. where; A P - p; = P + V - Xi + Yj + Zk; Z =-g. and the positive z direc ng -u Te li r s m is c i and the positive z direction is upward. The linear system is considered

-19to be vertical, with water injected at the bottom, at rate ut. The absolute saturation. S' is the dependent variable, and is a function of the independent variables x, y, z and time t. Examination of Equation (11) reveals that it is, in general, non-linear. The coefficients of the spatial derivatives are, except in special cases, functions of S. They are not always representable in simple analytic form. It seems highly probable that the solution of Equation (11) is beyond the power of presently known, exact, closed form mathematical techniques. If capillarity is neglected, Equation (11) reduces to the well known Buckley-Leverett equation, whose ramifications have been described elsewhere (107). In order to examine the effects of the various parameters in a qualitative sense, Equation (11) is subjected to the following normalizing substitutions, which were suggested by Rapoport (90) and Perkins (98), x y z u v w X = Y Z =; U V v, W tut where' = S' - S' gr wr The result is, when S = wr AS' a:dS L x Y J 2 dS 3a (12) gr\wr

-20where f (S) = rg dP C F1 PgUtL dS f3(S) = opg~ d(KrgFl) gUt dS vas _'2s a as as VS -2 + 62 + 2 It can be shown, for linear geometry in which the z coordinate is in the same orientation (upward or downward) as the gravitational effects, that the scaling problem reduces to a question of only two dimensionless functionals. Reducing Equation (11) to linear geometry and using only the substitutions Z = and T tut L ~ls' the result is OS dfl 6S 2 S df?2 (S 2 6T'dS 6Z + f __+ dS (0z) - (13) -~ ~dS az "2z s a = where fl = F1 - pgt (rgF1) f K = rF1 Pc 1gutL dS the functionals fl and f2 are measures of the frictional-gravitational, dP and capillaric effects, respectively. The derivative d may be obtained dS from experimental data, or from correlations such as that given in Equation (4) and Figure 4. The expression of rw -Krg F1 and Pc in terms of the normalized saturation S was a concept arrived at independently by the author. This is a logical concept, since these quantities are not defined outside the limits Swr to Sgr. Rapoport (98)was the first to formulate Equation (11). The present analysis has been brief in view of the previous work (98), which must be consulted for further detail.

-218. Discussion A qualification must be made to the basic Equation (11). Equation (3) is based on the assumption that there is a finite pressure jump across the gas-liquid interface when that interface is curved. For purposes of a theoretical analysis of the mechanics of gas-liquid flow, it is assumed that the capillary pressure can be treated in the manner implied by Equation (3). It is almost certain that the existence of hysteretic effects involved with the wettability characteristics of some sands relative to some two-phase systems must compel modifications of the present approach (i.e., that the contact angle is constant). For example, contact angles are known to be a function of whether or not the wetting phase is advancing or receding across the solid matrix (88). The equations, as derived, however, will be expected to be valid for the movement of gas and water in reasonably clean, non-reactive porous media. Generalized scaling analyses were done by Rapoport (98), and Geertsma et.al. (41). Each of these authors suggests at least four different dimensionless groups which must be either equal or proportional in model and prototype in order that scaling be achieved. The initial and boundary conditions must also be the same. The author feels that, while the two analyses are internally consistent', and correct as far as they go, they are far too restrictive from a practical point of view. Equation (12) shows that, if the three functions of S, dSl' f (S) dS 2 3 are the same, then the differential equation is the same, and there is nothing further to be said about the differential equation itself, as far as the scalling aspects are concerned. Fl' f2 and f are functions 1' 2 5

-22of S only. Fl, f2(S) and f3(S) can be termed the dimensionless scaling functions. The general scaling problem for multiphase flow may therefore be stated as; 1. Model and prototype geometry must be similar. 2. Initial, boundary and operational conditions (sequence of events) must be the same. 3. dF f (S) and f (S) must be the same. dS 2 3 dF1 It must be emphasized that,~ Y f2 and f3 are single valued functions of dS S alone, obtainable from simple experimental measurements which can be taken in a core laboratory. The general procedure for conducting a df1 scaled linear experiment is simply to make the functions - and f dS 2 equal in model and prototype as well as to satisfy conditions 1 and 2 above. Exception is therefore taken to the previous work in this area (41), (98) in that the author feels that it is easier to work with two groups than it is to work with four or more. It is one thing to say that the capillary pressure function and the viscosity ratio and the relative permeability must be the same or proportional. It is quite another to say that all that is necessary is for the function f2, for example, to be the same in model and prototype. In this sense, it is felt that the scaling procedure suggested above is far less restrictive than those previously suggested (41), (98), and that it may simplify considerably the planning of model work. Examination of Equation (11) yields the total number of parameters affecting the results of multiphase flow experiments. They are: uw, ug K, ut, p, L, P,, Krw and Krg' Other parameters such as pore size, pore size distribution, and contact angle do not enter directly into the picture. They will auto

-23matically affect the parameters listed above in such a way that the equivalence of the groups of Equation (12)ff. will produce scaling. B. Saturation Distribution During Constant Rate Displacement 1. Previous Solutions of Two Phase Flow Equations As has already been shown, there are three major forces acting on fluids in mulitphase flow in porous media - frictional, capillaric and gravitational forces. A rigorous formulation, on the basis of these three forces is given in Section IIA of this thesis, and yields Equation (11). The first attempt to solve the problem was made by Buckley and Leverett (15), who included only the frictional forces. A detailed exposition of the two phase flow solution from the Buckley-Leverett viewpoint is given in (107) p. 163. Although the method which Buckley and Leverett used to solve the basic equation was not strictly correct, they did successfully obtain the right answer, and the theory which bears their name has become the foundation of all non-capillary multiphase flow theory for porous media. The next step was the inclusion of the capillary pressure effects, which make the partial differential equation of saturation as a function of space-time a second order equation. The general equation of two phase flow is usually non-linear in all the space derivative terms. Recognizing this, Philip (93) proposed a numerical solution. The numerical solution is laborious, but is claimed to be an improvement over solutions (65), (92), proposed up to that time. The solution is for linear geometry. The next step in order of complexity was the proposal of a generalized numerical solution for multidimensional two phase flow with

-24the inclusion of gravity and capillarity. This was accomplished by Douglas, et.alo (32), who were able to observe some qualitative agreement between their numerical calculations and the experimental data. The most elaborate computer scheme treating this problem was given by West (133), who not only considered complex geometry, but took into account the variation of density and viscosity with pressure. Other theoretical work related to the solution of the basic two phase flow equations is given in (17), (18), (33), (57), (72) and (79). 2. Experimental Techniques Reported The obtaining of good experimental data in two phase flow through porous media is very difficult because of the problem of measuring dynamic saturation distributions accurately. One of the finest pieces of experimental work that has come to the author's attention is that of Leverett, et.alo (74). Model studies were done for the case of water coning in multiphase flow, and for the displacement of oil from sands of non-homogenous permeability. Saturation distributions were obtained by running the same experiment on a number of identical models, and stopping the flow at different times, Samples, which were mechanically removed from the sand pack at these times were then analyzed to determine the saturation. Other experimental techniques are outlined in (21) and (64).. The measurement of saturation can be done in many different ways. These include electrical resistivity (3), (75), radioactive tracers (25), (59), (130), X-ray techniques (7), (71), (83), magnetic susceptibility (134), and matched refractive indices (51), (116). It must be noted that the electrical resistivity method must be viewed with

-25reservation for dynamic flow studies, since Dunlap, et.al. (34) have reported that it takes 5 to 25 days for most cores to reach equilibrium with respect to electrical resistivity saturation measurements. 3. Experimental Two Phase Displacements Reported Determination of saturation distribution during frontal displacement of natural gas by water was made on laboratory scale models by Geffen, et.al. (42). The conclusions were that the same behavior will occur in the core as will occur in the field. Saturation distributions during oil-water displacements from a consolidated alundum core were obtained by Levine (75) using the electrical resistivity technique. However, his results were made less useful than expected due to a significant non-homogeneity of the core. Point values of the flow potential were also measured. Saturation distribution data during oi.-ater displacements in a sand pack were taken by Perkins (89). These data illustrate the effect of capillary forces on the flooding behavior of short laboratory cores. The information obtained here indicates that both inlet and outlet capillary end effects have an influence on the saturation history. Terwilliger, et.al. (121) presented saturation distribution data for the case of gravity drainage performance. The basic concepts related to the theory of the stabilized zone were also presented here. The general form of the saturation distribution in linear, two phase flow will now be discussed. Equation (12) and experimental evidence indicate that for sufficiently high injection rates ut, the BuckleyLeverett theory (13) applies. In those flows which can be adequately described by the Buckley-Leverett theory, the water saturation drops slowly from the inflow face, reaches the breakthrough saturation (107),

-26and then drops suddenly to zero, forming a sharp saturation shock. At low rates ut, the gravitational and capillaric effects assume importance, and the saturation front is more diffuse. Examples of this behavior can be seen in Figures 1 and 2 of (89), It must be emphasized that Equation (11) is not a rate equation in the strict sense of the word. It simply describes the variation of saturation with time and distance for a given injection rate ut, Equation (11) does not describe the pressure phenomena which give rise to the flow, but describes the consequences of a given total flow in terms of multiphase behavior. Background information on the saturation profile problem has been given in the introduction. In the next section, we will consider the description of the saturation front in linear, one-directional, constant rate, upward flow of water into an unconsolidated sand column which initially contains air and a small water saturation. The shock front is defined as the region in which the water saturation is changing rapidly in comparison to its rate of change (with respect to distance) in other sections of the porous material. For low enough rates, the shape of the shock front is close to that at equilibrium with respect to capillary imbibition. For rates between the latter and those flow situations in which capillary pressure is negligible, (Buckley-Leverett flows), the saturation profile is a function of rate. It is these intermediate conditions which will be considered here, since they provide a reasonable basis on which to obtain a meaningful comparison between prediction and experiment, i.e., a comparison in which all three (frictional, capillaric, gravitational) forces are acting.

-27The practical implication of the determination of saturation profiles during frontal displacement is easily seen when considering the transition zone between gas and water in gas storage reservoirs. Gas bearing strata may generally be considered of the order of magnitude of 10 to 100 feet thick. Leverett (73), on the basis of laboratory measurements, concluded that the transition zones can also be of the same order of magnitude in thickness. The transition zone is here defined to be the section of porous structure in which the saturation changes from a high to low value. It can be considered to be equivalent to the thickness of the saturation shock front. An example of the shock fronts to be investigated is given in Figure l(a). 4. Theoretical Development The physical problem to be treated here is that of water being injected into the bottom of a linear, vertical porous core having a low initial water saturation with the remainder of the pores filled with gas. The general method of solution outlined below is analogous to the Von Karman-Polhausen (109) integral technique for the approximate solutions to problems in boundary layer flow. The technique is described in detail on pages 255-260 of Knudsen and Katz (66). Normalization of the saturation by the definition S' Swr S wr S I- Swr gr wr in Equation (9) and Equation (11) gives the equation, -, aS = a F I u rFtl - + rgFlg + (14) Fu^ og ^g g

-28The expression for the volumetric velocity of water at any point is, Uw [ut + r —g (a - g) F1 (15) Consideration of experimental data and the general physical characteristics of the displacement process in porous media yields the following boundary and auxilliary conditions which must be satisfied by any solution, These conditions are, l, S (z, o) 0 2, S (o, t) = 1.00 3. S (zo,t) = 0 4. l(Zo t) = -A1 5. o < S (z, t) < 1.00 6s 6. - (z, t) < 0 7. (Z t) A1 - 6t ~~ at lim ns 8. (o, t) = 0 Boundary condition 1 expresses the initial condition that the (normalized) saturation is everywhere zero at the beginning of the displacement. Boundary condition 2 expresses the experimentally observed fact that the inlet saturation instantaneously rises to lo00 and remains there throughout the process, This is true primarily in liquid-gas displacements in relatively long cores. Boundary condition 3 states that the saturation drops off to zero at the point of farthest penetration

-29(zo) of the liquid. An illustration of a typical displacement process is given in Figure l(a), Boundary condition 4 implies that the saturation gradient at the leading edge of the displacement front (zo) is constant throughout the process. Conditions 5 and 6 imply that the saturation must always lie between 0 and,OO and that the saturation gradient must always be negative. The former condition is evident from the physical nature of the process and the latter is derived from the experimental observation that saturations always decrease from the inlet to the leading edge of the front Boundary condition 7 is derivable from the definition of the partial derivative, This derivation is given in Appendix A, Condition 8 expresses the restriction, derived mainly from experiment, that all spatial derivatives of the saturation are zero at the inlet after a large number of pore volumes of displacing fluid have been injected, This is illustrated in Figure 6 of (42), All the physical experience heretofore gained with displacements of this type indicates that the displacing phase usually forms a saturation shock front. This means that the water invaded zone has a definite length z0 beyond which the saturation corresponds to that initially given., The saturation is observed to drop off fairly rapidly near zo, This behavior is described by the statement of the boundary conditions, as given above, It must be emphazised that the integral method, as applied in the present work, depends on the assumption that a function which is required to satisfy nearly all of the critical properties of the desired solution will consequently approximate all of the other properties of the solution. If this assumption is to be realized, a reasonable functional form is then developed into an approximate solution by means of the satisfaction of all the boundary conditions and the partial integral of the differential equation.

-30The wetting phase saturation is assumed to occur in the form of a series, (16) S(z t) = C + C1 () gl(t) + C ( )2 + * C z (Z)gn +... S(z, t) =Co+ i + ~ 2 n zO The procedure used is to calculate the constants C and functions gi(t) from the boundary conditions 1 through 8. The series must be truncated so that the number of constants in it be equal to the number of numerical boundary conditions available. The more boundary conditions the function is required to satisfy, the more likely it will be that the approximate solution will cme close tothe exact solution. The function satisfying all of the conditions 1 to 8 is, S = - Al o (17) z where = Equation (17) was obtained by truncating Equation (16) after the first two terms. We now have a saturation functions S(z, t) which satisfies all the boundary conditions associated with the problem. The satisfaction of the integral relation derived from the basic differential equation will now be examined. The result of the integration of the basic differential equation will be the determination of the relation between z and t o (zo is a function of t only). This zO (t), together with the computation of A1 will provideI the quantities for use in predicting the saturation profiles according to Equation (17). Integrating both sides of Equation (14) from 0 to zo gives; f - dzPc dz = lut - g kg dS (18)

-31the integral on the left is, by Leibnitz' rule, o oS f z dzo as 0 at Udz = M S0 a Sdz - AS' S(zot) dt (19) 0t z the last term on the right of Equation (19) is zero by virtue of boundary condition number 3. Evaluating the left side of Equation (19) Zo a 1 t 0 0at =0 a f zo [so(l - lo ) ] dt= a A Z2 +1 But, this expression is also equal to the expression on the right side of Equation (18). z ZO AS' f I dz =- - F lUt. KKrg.Fl (20) o at L lg ig dS o Noting that rw = 0 and Krg = 1 at zO, and Krw = 1 and Krg = 0 at z = 0, and using the expression for F1 in Equation (11) ff., we obtain from Equation (20) Zo'0 at -dz =ut (21) Making use of Equations (19) ff. and (21) gives Alz 2NST utt (22) Alzo+ 1 0 From which Al1uit Altt 2 AA1u'utt + ut Zo = (23) 2A 6S' 1

-32Turning now to the evaluation A1, we place Equation (15) in the fractional flow form, uw =r(alPcgs ] Utg L dS Iz Ut [ g1 + a-( The stabilized zone concept of Terwilliger, etoalo (121), and Jones-Parra at (57) will now be used to evaluate the quantity S at z, ioe,, at the leading edge of the shock fronto The basic feature of the stabilized zone method is that w is assumed to be a linear function of S from S = 0 ut to the saturation corresponding to the Buckley-Leverett breakthrough saturation, designated S'o Experimental evidence supporting the nearlinearity ofW with S is given in (75), and (79) For a complete ut discussion of the determination of S', see Figure 3 of (57), or (107), page 165. The stabilized zone method is not applicable to saturations above the cutoff, or Buckley-Leverett breakthrough saturation, the solution becoming unbounded in this region. It cannot, therefore, be used to give saturation distributions over the complete range of space and time. It is desired to solve Equations (24) for the saturation gradient Hi at the point zo, i.eo, where S is zero, It will be assumed, uw by analogy with the stabilized zone method that - is a linear function ut of S only in the immediate neighborhood of zo, i.e., near S = 0. (Stabilized zone theory assumes its linearity over the whole saturation range.) In addition, it is assumed that his linear function of S is obtained by the method of the stabilized zone, i.e., by the Welge (132) tangent construction, which yields

-33uw U- a2S ut Solving Equation (24) for az as dP, c ~W azcrC n cr 1 ( 26 )_ ~ r [r ^__ 7 (26) aS gUt Uw Kr rg APg Krg g Ut [l+rw[g - + utg "rgUw rw rW At low saturation gw > > 1 and g = 1.00, and >ut ~w" >> 1 grwg Ut rw Making all these substitutions in (26) and simplifying gives dPc z ut [ d w S g] (27) S z U I w + h Kg Ut Krw Ut dPc The value of d- near zo (ioe,, near S = 0) can be read from the static capillary imbibition curve. The relative permeability K is expanded in a power series rw in S = S + D o Dn +. (28) Kw D1S + D2S2 + oo DnS + which reduces to K = DS for small S. rw 1 Substituting this and Equation (25) into Equation (27) gives dPc az K K dS Zo 1 (29) a zo ut a2 w + ApgK A1 0 1 D1 Ut For very small values of ut, as dPc (30) az ApgdS

-34which is the capillary equilibrium curve. The equations giving the solution are (17), (23) and (29). Al and zo are determined from (23) and (29), respectively, and are then put into Equation (17) to give the required space-time distribution. The problem which has been solved is that of the saturation distributions arising when water is injected at constant rate into the base of a vertical, linear core initially containig connate water and air. 5. Experimental Results The purpose of this section is to present experimental data taken on the physical system corresponding to the problem outlined above. These data will then be compared to the predictions developed in that section. Experimental data were taken for water displacing air upward from a vertical glass tube which initially was filled with sand, air, and enough water to give a saturation corresponding to the irreducible minimum water saturation. A glass tube (1.15 inches inside diameter) was packed uniformly with 150-200 mesh clean Ottawa sand. The pack was formed by allowing the sand to settle through the water at the rate of 30-35 gm. per minute, with constant vibration of the tube. The vertical pack was then blown down to 10-15% water saturation by flowing air downward through it at 20-25 psi pressure drop. The packs were 3 to 4 ft. in height. Water was then introduced at the bottom of the vertical core at a constant rate. This rate was controlled by a Ruska pump having a range of 1 cc/hr to 224 cc/hr to a precision of 1%. The properties of the sand pack and fluids used are given in Table I. The capillary pressure curve was obtained by placing the bottom of a sand pack (containing 10-15% water) in a beaker

-35of water, and letting the water imbibe to equilibrium, The pack was then divided into cylindrical sections, each of which was immediately sealed to prevent evaporation. The sampling operation for the most part consisted of pushing the pack from the tube in the form of a slug and cutting off cylindrical sections (whose diameter was that of the tube) as they emerged from the tube. The pack was held vertical during this operation, and the total sampling time was 3-4 minutes, In some cases, samples were removed from the tube by displacing them with a flat steel strip, The samples from the saturation transition zone were taken first in all cases, The saturation distribution in adjacent layers was not disturbed during the sampling procedure on a specific element. The samples were then weighed, muffled at 175-1800C for 20-25 hours, and then weighed again, The average saturation in each element was then calculated from the known pore volume of the dry sand and the total volume of water that had been evaporated from the pore volume, TABLE I Properties of Sand and Fluids Used in Frontal Displacements Sand grain size 150-200 mesh Permeability 11o7 darcies Porosity, fraction voids.409 Water density 62~4 lb m/ft3 Water viscosity (68cF) 1.00 cp Air density.0749 lb m/ft3 Air viscosity o0184 cp

-36The sand was strongly water-wet. The capillary pressure data are plotted in Figure 4. The saturation distribution curves under the dynamic flow conditions were measured in the same way as was the saturation distribution under capillary equilibrium, the pump being shut off at a specific time after the start of flooding and the samples taken immediately. A new sand pack was used for each run. The results are shown in Figures 6 to 14. Figure 10 shows a profile development with time. In the series of profiles shown in Figure 10, the rate of injections was kept constant, and the distribution taken at various times. Figures 6 to 14, except for Figure 10, show profiles for different rates, but at only one stage da-w in the flood. The value of D = - =.025 at very low saturations was 1 dS measured by the relative injectivity method (56). The results are presented in Figure 15, and indicate that the curve does have a finite slope at S = 0. The value of a = 1.00 was determined by the Welge tangent construction method (c.f. Scheidegger (107), p. 166). The experimental technique permitted a precision of the order of 3-4%. The solid lines in Figures 6 to 14 represent the predictions indicated by Equation (17), Equation (23) and Equation (29)~ The relative permeability for the 150-200 mesh sand used in these experiments was measured by the relative injectivity method of Johnson, et.al. (56). The data obtained are plotted in Figure 3 for comparison with the average curve for unconsolidated sands determined by Wyckoff and Botset (137). The average deviation of the 185 data points in Figures 6 to 14 from the predicted curves is 3*9% in spite of the fact that small disagreements tend to cause very large errors in the

0.9 0.8 0.6 ~~~~~~__ _~~~~ (n z I- O 0. W. 05 I 0.4 z 0 THEORY - I I I 0 EXPERIMENT - I ~ ~ t 132 HOURS u =0.00490 ft/ft-hr 02 INITIAL WATER SATURATION, <\ 0.1 0 0 0.2 0.4 0.6 08 1.0 1.2 14 1.6 1.8 20 22 24 26 28 30 z(ft),VERTICAL DISTANCE FROM INFLOW FACE Figure 6. Saturation Distribution During Displacement of Air by Water from 150-200 Mesh Sand.

0.9 0.8..... 0.7 ~ 3 0.6 z 0 THEORY 0.4 0 EXPERIMENT t 70.2 hrs 0.6 ~ ~ ~ ~~0.0126 ft~/ft -hr _ \ ir~~~~~~~~ 03 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 40 4.2 0.2.. 0.1 ff r 0 02 04 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 24 26 28 3.0 3.2 34 3.6 3.8 40 42 z(ft), VERTICAL DISTANCE FROM INFLOW FACE Figure 7. Saturation Distribution During Displacement of Air by Water from 150-200 Mesh Sand.

0.9 0"~ ^ 0 0 0 ^^ ^0 0.8 0 0.7 ~ 0.6' \ z 0 cr Q5 C,) 4 \ I Q4 ~ ~ - ~ THEORY - ~ ~ ~ ~ ~ ~ ~ J 0 EXPERIMENT g t =8.37 hrs \ u:0.0441 ft/ft-hr 03 0.2 0.1 0 _. __I____ ______!~lLSIjTO_ ____________ _____ INITIAL SATURATION _ 0 0.I 0.2 0.3 04 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 z(ft),VERTICAL DISTANCE FROM INFLOW FACE Figure 8, Saturation Distribution During Displacement of Air by Water from 150-200 Mesh Sand.,

Q9 0.8 0 I 0 0.7 0 0.6 -U~ 0 z E Q5 C,) w 0.4 THEORY -__. 0 EXPERIMENT 0 0 z t =1.88 hr~0 I. - u =0.0735 ft3/ft-hr 0.3 LL 0 Q2 INITIAL SATURATION 0.1 0 0.1 02 0Q3 0.4 0.5 0.6 0.7 0.8 09 z(ft),VERTICAL DISTANCE FROM INFLOW FACE Figure 9. Saturation Distribution During Displacement of Air by Water from 150-200 Mesh Sand

Q9 El ~0.9~0 0 o 0.8 7A N ~ X ^~~ 0 X 0.6::: THEORY a I \0 EXPERIMENT Mc u =0.157 ft/ft-hr 0 0.4 z 0 01 0.2 0.3 04 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 Q 0.2 0 1 0.2 0.3 4 0.5 0.6 07 0.8 0.9 1.0 1.1 1.2 1.3 14 1.5 1.6 1.7 1.8 1.9 2.0 2.1 22 2.3 25 z(ft) VERTICAL DISTANCE FROM INFLOW FACE Figure 10. Saturation Distributions During Displacement of Air by Water from 150-200 Mesh Sand.

0.9 0.8 0 o 07 ________ 0\ 0.6 0 1 05 0.4 THEORY 0 EXPERIMENT z u =0.157 ft/ftE-hr ir i t = 1.92 hrs. 0.2 0.1 0~- __ _ INITIAL SATURATION 0 0.1 0.2 0.3 0.4 05 0.6 0.7 0.8 0.9 1.0 1.1 1.2 Z(ft),VERTICAL DISTANCE FROM INFLOW FACE Figure 11. Saturation Distributions During Displacement of Air by Water from 150-200 Mesh Sand.

0.9 ---- C 0 08 - - -- - - - - - 0.8 0~~~~~~~~~~~~~~~~~~~~~~~~ \ 0 O~~~ 0.6 ~ ~~~~~~~~~~ z 0 ^ 05 ^ ~ ~~~~~ THEORY \ IW 0 EXPERIMENT 0.4 ~ ~ ~ - u=0.157 f /ft-hr s? t 6.09 hrs. \ ~f=~ \J z Sr 0.3 ~~~~~~~~~- ~ QLL 0.2 0 INITI L SATURATI N ___ ___ ________1MSAUAN___ __ ______ ____ ____ ___ \ _ 0.1 ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ 0 0.2 04 0.6 0O 1.0 1.2 1.4 1.6 18 20 22 24 2.6 28 3.0 3.2 3.4 3.6 3.8 40 z(ft), VERTICAL DISTANCE FROM INFLOW FACE Figure 12. Saturation Distributions During Displacement of Air by Water from 150-200 Mesh Sand.

0.9 ---- 0 0 0 0.8 0 0 0 0.7 06 05 n,C,! 0.4 ~ ~ ~ ~ ~ THEORY 0 EXPERIMENT Z t =4.45 hrs 0.3 00.12 INITIAL SATURATION 0 0.2 0.4 06 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 24 2.6 2.8 30 3.2 34 3.6 3.8 4.0 z(ft), VERTICAL DISTANCE FROM INFLOW FACE Figure 13. Saturation Distribution During Displacement of Air by Water from 150-200 Mesh Sand.

0.8 ~~~~~~ o ~~~ 0 0 0.7 ~ ~~~~~~~ 0.6 0.5 4P THEORY C,, Q3 I-. 02 ~ INITIAL SATURATION 0 02 04 0.6 OB 1.0 1.2 1.4 16 1.8 20 22 24 2.6 2B 30 32 34 3.6 3.8 4.0 4.2 44 46 z(ft), VERTICAL DISTANCE FROM INFLOW FACE Figure 14. Saturation Distribution During Displacement of Air by Water from 150-200 Mesh Sand.

0.007.,~~ 0 0.006 Q005 ~~..~ ~~~~ 150-200 MESH SAND LOW SATURATION REGION w I-.U 0.004 0: I —0.003 W 0 004 0.08 0.120.025.28 0.32S 0.002 / w / 0.001 ~ / ^0 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 0.36 0.40 FRACTIONAL WATER SATURATION,S' Figure 15. Relative Permeability versus Saturation.

_47portion of the curve where the slopes are large. The prediction curves of Figures 6 to 14 were calculated as follows. The saturation profiles were calculated from Equation (17) using the computed values of Al and zo. A1 was computed from Equation (29) with D1 =.025 obtained from Figure 15, dPc and a = 1.00 obtained by the Welge tangent construction technique. dS at S = 0 was obtained from Figure 4. The depth of penetration z was calculated from Equation (23) from a knowledge of ut, t and Al. The calculation procedure implied in Equations (17), (23) and (29) was used in an attempt to predict the results of some experiments performed by Geffen, et.al. (42) on Nellie Bly sandstone (K = 590 md. (air), and 0 =.265). Natural gas was displaced upward from a vertical core first by a constant rate injection of water and then by simple cocurrent capillary imbibition. The calculation procedure was the same as that used for the unconsolidated sands. The comparison of prediction with experiment is shown in Figures 16 and 17. The equations are seen to predict the experimental saturation profiles satisfactorily. The volumetric velocities for the imbibition displacement were determined by calculating the area under the saturation vs. distance experimental curve. 6. Discussion The solid lines of Figures 6 to 14 are seen to deviate from the experimental data by about the same order of magnitude as the precision of the data themselves (3-4%). It must be emphasized that the theory developed permits semi-theoretical predictions to be made when Srw, Srg, capillary pressure and relative permeability data are known. These solutions do not involve curve fitting or the introduction of arbitrary constants. Since the agreement of theory with the data is as good as

0.9 0.8 07 I m r a 0 \, \ 0.6 O 0.5 0.6 \ ~ ~~~ ~\ ~~~~ 01.1~~ 0 gig \ur 15t.25 hrs I =3.1hrs N G t=5.05 hrs UN 0THEORY — J I o co __ INITIAL SATURATION 0.2 ~ 0 00 EXPERIMENT DATA FROM FIG.6 0F(42) u =0.026 f?/ft2-hr 0.1 0 0.1 0.2 03 04 0.5 0.6 07 0.8 0.9 1.0 1.1 1.2 z (ft),VERTICAL DISTANCE FROM INFLOW FACE Figure i6. Saturation Distribution During Displacement of Natural Gas by Water from Nellie Bly Sandstone.

0.9 "~~ THEORY Q8 * ~ O^I EXPERIMENT DATA FROM FIG. 9 OF(42) 0.7 0. z Q2 0.0.4r 1 0.34~~\-~ l _______I I_______ _______INITIAL \ SATURATON_____ 0.2 0.1 0 0.1 0.2 0.3 04 0.5 0.6 07 0.8 0.9 1.0 z(ft), VERTICAL DISTANCE FROM INFLOW FACE Figure 17. Saturation Distribution During Displacement of Natural Gas by Water by Capillary Imbibition from Nellie Bly Sandstone.

-50the precision of the data, there would be nothing to be gained by curve fitting in any case. Inasmuch as agreement of some experiments with theory does not necessarily establish the theory, the validity of the results will now be discussed. The major point at issue as far as physical theory is concerned is that of the effects of the capillary pressure, and whether or not the capillary forces cause flow in exactly the manner given in Equation (3). The development of the basic equation is, after all, only the result of a formal procedure. It involves the supposition that there is a finite pressure jump (Equation (3) ) across a curved interface separating two immiscible fluids. One author goes so far as to say that no such pressure difference exists (54). In any event, it must be pointed out that the exact nature of capillary pressure and capillary rise is still not precisely known, as evidenced by a recent publication in the area of physical chemistry (114). For purposes of multiphase flow in porous materials, however, it is considered that the agreement of the theory with experiment (Figures 6 to 14) appears to substantiate the validity of the method. Another point of discussion is related to the validity of Equation (2) to properly represent the volumetric velocities during multiphase flow. These equations state that the velocities in multiphase flow can be derived from a potential in the same way as in single phase flow, except for the relative permeability function. These relative permeability equations are not theoretically derived, and may or may not represent the what actually goes on in porous media. Work by Irmay (55) shows experimental data indicating that Equations (2) are

-51valid provided both phases are continuous The equivalence of the static and dynamic capillary pressures was discussed previously. The equivalence of dynamic and static relative permeability curves is also a problem. Relative permeability data are sometimes taken with the fluids flowing at steady rates and the saturation in the core unchanging. The question arises as to whether these static curves can be used to describe the situation occurring when the saturations are varying considerably with space and time. Johnson, etoal. (56) compared relative permeabilities measured under steady state conditions with those measured by the dynamic relative injectivity method. The curves were identical. It may be concluded from this that steady relative permeability curves may be used in calculating the behavior of unsteady state systems. Some final comments with regard to the physical theory must be made, The development of Section IIA presupposes a gas-solid-liquid contact angle which is the same as the static contact angle and which is independent of rate. The difference between advancing, static and receding contact angles has been previously documented (1), (19), and (80). Changes in the value of contact angle will generate changes in the capillary pressure, as shown by Equation (4), This does not happen with the silica-air-water system as is shown in Table I of (80), This is the reason that sand-air-water was chosen for the experimental work of this thesis. It must be emphasized, however, that the contact angle for other systems may be a function of rate of liquid advance, or at the very least a function of whether or not the liquid is advancing. The equations developed in Section IIA must, therefore, be used with care.

-52In particular, the static capillary pressure curves must be modified to properly allow for the change in contact angle with advancing or receding wetting phase. A large amount of data related to contact angle is given by Melrose (80). Kinetic effects with respect to wettability may also occur in flows in which a wetting phase is displacing a more viscous non-wetting phase. Putting aside all question of hydrodynamic stability, we may imagine that the rate at which the wetting phase actually wets the solid surface is a function of the viscosity of the in-place fluid as well as of the simple adsorption kinetics. Alternately, the time taken for the wetting fluid to move past the viscous oil to the solid surface. may be of the same order of magnitude as the time taken for the flood front to move past the element of surface in question. The flood water may then behave as if it were the non-wetting phase when it actually would wet the solid if given enough time. A problem of considerable proportions in oil reservoirs is that of wettability reversal. In physical terms, a core which is initially strongly water wet can be made oil wet by placing it in contact with crude oils containing certain components. These components appear to be adsorbed on the rock surface. The longer the surface is in contact with these oils, the more oil wet it becomes and vice versa with respect to water. In water flooding of oil by water, there may be a change, or even reversal of wettability near the flood front, depending on the kinetics of change of contact angle. This change in wettability did not enter into the frontal displacement experiments done in this thesis. However, such changes must always be consered as a possibility when

-53examining other multiphase displacements, especially with two-liquid phases. Quantitative data on the wettability reversal phenomenon are given by Denekas, Mattax and Davies (30)0 and Wagner (127). The approximate method was devised because oi the very considerable mathematical difficulties involved in the exa ct solution of Equation (14) even when the simple linear geometry is assumed. The success of the method in predicting the experimental data, together with its simplicity, are its principal justifications since it is not entirely rigorous. In view of the approximate nature of the solution presented, the validity of the results and the engineering precedent supporting the technique will now be examined. The following discussion (extending to the bottom of p. 56) illustrates other uses of the integral technique and may be omitted by those familiar with these applications. Much attention has been turned in the last fifty years to methods for the solution of the (non-linear) Navier-Stokes equations. Schlichting (109) has summarized many of the results of the boundary layer approach to the problem of viscous flow. The essence of the boundary layer approach is that the effects of the viscosity are confined to a liquid layer of definite thickness which lies in the immediate neighborhood of the solid bounding surface. Approximate methods of solving the boundary layer equations have been devised by many workers and have been summarized by Schlichting (109), page 238 ff., who introduces the problem in the following manner. "The general problem involving the flow of fluid around a body of arbitrary shape, which is particularly important in practical applications, cannot be solved with the aid of the analytical methods developed so far...."

"It is, therefore, important to devise approximate methods which would in such cases quickly lead to an answer, even if their accuracy were to be inferior to that of the numerical methods. Following Th. Von Karman and K. Polhausen, it is possible to devise such simplified methods if it is agreed to satisfy the differential equations of boundary layer flow only in the average, and over the boundary layer thickness rather than to try to satisfy the boundary conditions for every individual fluid particle. Such a mean value function can be obtained from the momentum theorem. This, in turn, is obtained as an integral of the equations of motion over the boundary layer thickness." The momentum equation referred to by Schlichting is the partial integral of the governing differential boundary layer equation, integrated from the solid wall to the edge of the boundary layer. Its physical meaning expresses the fact that the shearing stress at the wall is equivalent to the loss of momentum in the boundary layer. The analogy between this Von Karman-Polhausen method and the solution given in Equations (17), (23) and (29) may now be detailed as follows. The analog of the boundary layer thickness is the depth of fluid penetration into the porous material (the zo of Figure l(a) ), The governing differential Equation (14) is integrated from z = 0 to z = z0 (Equation (22) ) in the same way that the governing boundary layer differential equation is integrated from the solid wall (y = O) to the edge of the boundary layer (y = 6) to yield the total momentum balance. The Von Karman-Polhausen method is completed by expanding the velocity in a power series of similarity parameters and computing the constants in the series from the boundary conditions (cf. Schlichting

-55(109), po 244). Of course, the necessary auxiliary (compatibility) conditions are satisfied by the form of the function chosen. The shape factor parameter is obtained from the partial integral equation (momentum equation), cfo Schlichting, page 246-8. The similarity between the Von Karman integral technique and that given in this thesis is complete, except for the calculation of the function z (t) (Equation (23) ) in the present solution in place of the shape factor function. Also, the saturation profile does not have the similarity characteristic possessed by the boundary layer velocity distribution. Schlichting makes the following comments regarding the agreement between the exact and approximate methods. "It is seen that the approximate method leads to satisfactory results in the case of a flat plate at zero incidence, and the extraordinary simplicity of the calculation is quite remarkable, compared with the complexity of the exact solution." "The agreement between the approximate and the exact values (of the solutions to the two-dimensional boundary layer equations) is here also completely satisfactory" (Schlichting, pp. 251-2). "No general criterion regarding the admissibility of the approximation has been given so far, and it seems that it will be difficult to obtain, Judging by the above and similar calculations, as well as by experimental results, it appears certain that Polhausen's approximate method leads to very satisfactory results.... " (Schlichting (109), P. 253). The Von Karman-Polhausen method has also been outlined by Knudsen and Katz (66), p. 258.

-56Another application of the integral method was given by Squire (45) for the solution of the problem of natural convection from a vertical flat plate. Partial integrals of both the momentum and energy equations are obtained and evaluated over the boundary layer thickness. Expressions for temperature and velocity distributions satisfying all the boundary conditions are written. Substitution of these expressions in the partial integral equations gives results for the velocities,, boundary layer thickness and Nusselt number which are quite close to those computed by the exact method. A recent application of the Von Karman integral technique has been achieved by Friedlander (40) for diffusion from spheres. The partial integral of the diffusional transport equation evaluated over the boundary layer is obtained. It expresses the fact that the (steady state) amount of material diffusing from the surface at any point must be balanced by the convective transport through the boundary layer at that point., An expression for the concentrations, satisfying all the boundary conditions except one, is substituted into the over-all balance equation to obtain the solution. Agreement with experimental data is good. Vetrov, Petrenko and Todes (126) used the integral technique to obtain an approximate solution to the problem of non-steady gas motion in porous media. It must be emphasized that the approximate method is dependent upon a reasonably accurate knowledge of the slope of the capillary pressure and relative permeability curves at the saturation occurring at the leading edges of the displacement front (cf. Equation (29) ff.). When this saturation is the irreducible minimum wetting phase saturation,

-57the determination of these slopes can become quite difficult if the proper experimental data are not at hand. Examination of Equation (17) As will show that the slope aZ at zo must not be 0 or infinite. This requires a careful examination of the capillary pressure and relative permeability curves at the saturation at z o The approximate method can also be used for linear floods in which the saturation at the leading edge is not equal to the irreducible minimum wetting phase saturation. This involves the evaluation of dPc Uw the quantities dS-, and Krw at the saturation specified. Figure 11 dS Ug shows that excellent agreement is obtained between the experimental data and approximate solution even though the sand pack was at zero absolute saturation (S') at the beginning of the run. In the final analysis, the approximate method presented here must be considered as a semitheoretical model which successfully predicts the outcome of flow experiments to an acceptable degree of accuracy. The success of any semitheoretical model must be judged not only by the ability to predict absolute values of the desired quantities, but also the slopes of the curves. Consequently, experimental data shown in Figure 10 were used to form the derivative aZ and these were compared to the gradient predicted by Equation (17). The horizontal lines in Figures 18 and 19 represent the experimentally measured saturation gradients. The theoretically predicted saturation gradients (solid lines) are seen to lie quite close to the measured gradients in Figure 18. Figures 18 and w illustrate the beost and worst agreement obtained between predicted and measured gradients. Most of the curves (Figures 6 to 14) show agree

3.0 2.5 FROM 150 -200 MESH SAND DATA FROM FIG. 10 JN | WITH t 1.0Ohr o20. 1 THEORETICAL GRADIENT (EQ.17) o.L.. z 4 O ___ 02 000Q7.9 CD Z i Q5 05 0 0.1 02 0.3 0.4 0.5 Q6 0.7 0.8 09 1.0 z(t), VERTICAL DISTANCE FROM INFLOW FACE Figure 18. Saturation Gradients During Displacement of Air by Water.

3.0 2 _5 _ | FROM 150-200 MESH SAND ~2.5~ (DATA FROM FIG. 10 WITH t=2.00hrs) N 2.0 Z/ THEORETICAL GRADIENT,(EQ. 17) 1 1.5 z I- 1.0 C,) 0.5 VRIA 0 0.2 0.4 0.6 0.8 1.0 1.2 14 1.6 1.8 z(ft), VERTICAL DISTANCE FROM INFLOW FACE Figure 19. Saturation Gradients During Displacement of Air by Water.

-60ment intermediate between these two extremes. The conclusion can be drawn that the data stand up well to a first differentiation and that the integral technique predicts the saturation gradients satisfactorily. In summarizing the support for the validity of the approximate solution to the multiphase flow equations, the following are mentioned. 1. The approximate method agrees with experimental data on an unconsolidated sand to within 5.9%. 2. The approximate method successfully predicts the saturation gradients occurring in displacement of air by water from an unconsolidated sand. 3. The approximate solution predicts saturation distributions during displacement of natural gas by water from a consolidated sandstone. 4. The approximate method is an adaptation of an accepted engineering technique which has been successfully used by Von Karman and Polhausen (109), Squire (45), Friedlander (40) and others. The approximate method is applicable to the calculation of saturation distributions in natural gas and oil reservoirs when the fluids may be assumed to be moving linearly, for example, in a bottom water drive. The method can also be applied to the movement of fluids in filter cakes and packed beds. C. Steady State Capillary End Effect It has been observed for many years (31) that anomalous saturation behavior frequently occurs in two-phase flow experiments in short cores. Specifically, water will build up at the outflow face of

-61the core, resulting in average residual water saturations much higher than for long cores, Increases in the desaturating gas velocity are observed to decrease this residual water. This behavior is due to the discontinuity in capillary pressure at the outflow face of the core and is called capillary end effect. The water buildup at the outflow face causes a drastic reduction in the permeability to gas. This means that very much higher pressure drops are required to flow a given quantity of gas through the core. The problem is visualized physically as a wetting phase such as water being driven out of a porous core by a nonwetting fluid such as air. At equilibrium, the water will be immobile and the air will flow through the core at a constant mass velocity. The water saturation will build up at the outflow face because of the capillary forces which tend to draw the water into the body of the core rather than letting it flow out. This outflow saturation buildup is always present even in dynamic two-phase flow experiments. Additional discussions are found in Leverett (73), Kyte (70), and Dombrowski (31), Kocatas and Cornell (67), Brownell and Katz (11), Brownell (10), and Buzinov (14). The problem then arises concerning methods to predict the magnitude of the end effect from simple measurements on the fluid and rock properties. The saturation will be highest at the outflow face, dropping off to the irreducible minimum at large distance from the outflow face. The forces acting to desaturate the core are gravity and friction. The force holding the fluid in the core is capillarity. The problemrof the end effect with both phases in motion and the outlet saturation corresponding to residual equilibrium non-wetting phase was treated by Richardson (102).

-62Silverblatt (112), Nelson (86), and Dombrowski (31) have presented data for filter cakes. Haruni (49) has presented data for centrifuge cakes, Nenninger (87) presented a theory for the unsteady dewatering of filter cakes, Of course, the work presented in this paper is applicable to prediction of moisture content in filter cakes or any porous material, provided the irreducible minimum wetting phase saturation is known. In the remainder of this section a theory of end effect will be developed and compared with some data of Dombrowski (31). 1o Theoretical Development The starting point for the quantitative treatment of the end effect will be Equation (15) with the positive z direction this time being taken as downward, at equilibrium uw = 0, ut = ugo The final equation is F)dPc u + dP Apg1 ] (34) g 1.9 dz To obtain the equation in more explicit form, we know that (Cof, Equation (4) ) P = f(G) ( /2J(S) c K where J(S.) is Leverett, J function, Leverett found that for a wide variety of unconsolidated materials, Equation (4) would describe the experimental data. J(S) is apparently a universal curve, From Leverett's paper we can represent the J function by the following curve: o 110 J(s) =.200 + 0110 S (5

-63where S is the normalized saturationo S' is.070 for the materials under consideration, S is the irreducible gas saturation or the residual gas gas saturation, It is estimated at.90 on the basis of the experimental. data of Dombrowskio The definition of normalized S is necessary so that the capillary pressure data can be fitted to an equation of the form (35). These manipulations serve to fulfil the conditions that J(S) becomes very large at S = 0 (i.e. at S' =.070 rather than S' = 0), for drainage. With the above definitions' we can write the following equation, noting that S is to be the dependent and z the independent variable. dPc dPc dS -.110 Of(Q) ()l1/2 dS o ( ) -- (36) dz dS dz S2 K dz From the data of Wyckoff and Botset, (137) it is seen that for saturations S' up to 75% the relative permeability to gas, for a wide range of unconsolidated materials, can be represented by rg = 1 - 12 S (37) Substituting all these results into Equation (34), we get: ug + -.llOaf(0) (1/2 1 + (38) Multiplying through by Ug and rearranging, g+ tpg --.ollf(Q) (j (39) (1-1.25S) K sg2 dz

Let -.110 f( ) ()/ =2 A K Apg = B lgUg' K Equation (39) is an ordinary differential equation in and z which can be solved by separation of variables. AdS dz =.. 25. + B 1-1.25S and z = —E dS 1+E 1.25 dS (40) S2(1-1.25FS) S(1-1.25FS) where A B E =; F =; D=B+C D D Integrating Equation (40) z = E +1 1.25 (F-l)ln(1-1FS ) + ln(const.) (41), S Using the condition that S = S1, at z = 0 and taking antilogs, the final functional relation between S and z is t (1 1 iS 1-1.25FS1 The computation of the end effect for radial coordinates in the absence of gravitational effects proceeds as follows:

-65Krg dP Ug -ag dr Going through the same development as before, the analog of Equation (39) is obtained as, (1 9 ADS) -.110af((3 t)1/2 1 dS ~_ -= -.,loaf(eS)2 dr (1-1.25S)K 2r Integrating this equation and using the boundary condition S = S at r = rB, the final saturation distribution is r i 1 -- + 1,25 In + rB -C L S where A and C are the same as in Equation (39) ff. Calculation of A, B and C for various cases of commercial importance will readily show that all three of the viscous, gravitational and capillaric forces are usually significant. It is often the case in the desaturation of a porous body that the pressure drop is specified rather than the volumetric velocity (u) of the desaturating phase. The method of handling the problem is then as follows: The core of length L is divided into n' differential elements. The ith element is of width zi and has a pressure drop across it of -p dependent on the saturation of the ith element S. mp = E mi i -1 _Uigi i Krgi

-66n X = _ ui1gti i=l S^rgi But all mass passing the ith element must pass all the others so that u is a constant. Then, u _.~. n' z PgA Zi i=l Kirgi Passing to the limit u = (43) L dz Pog f r o Krg The procedure in the case that ~p only is specified will be to estimate u and calculate the saturation distribution from Equation (42). g u may then be computed from Equation (43) recognizing that k is a known g g function of S. The process is repeated until the original estimate agrees with the final calculation from Equation (42). Of course, if the desaturant mass flow rate is given directly, this procedure is not necessary and Equation (42) may be used directly. The following table shows some of the data on the systems used by Dombrowski (51)

-67TABLE II Properties of Materials Used by Dombrowski (31) Sphere size Permeability Porosity Bed Surface lb-f) Contact* _(sq.ft.) Length(ft) tension ft angle() 100/115 mesh 9.17 x 10-11.367.339 2.44 x 10-3 25~ 60/65 mesh 4.03 x 10-10.356.623 2.44 x 10-3 25 35/42 mesh 1.03 x 10-9 367.984 2o44 x 10-3 25~ Data from Newcombe etoal. (88) for air, mineral oil, and glass beads. Dombrowski (31), presents some data indicating the saturation distributions in packed beds before the pressure differential is great enough to cause gas flow. To visualize this situation, we can consider an analogy with liquid rise in a glass tube due to capillarity. If the free liquid surface and the liquid surface inside the tube are at the same pressure, the liquid will rise to a height equivalent to the capillary pressure. However, if a small pressure is put on the gas in the tube, the liquid level will fallo This is what is happening in the porous beds, except that there will be a saturation distribution. Setting ug = 0 in Equation (38) and integrating A l z = -(- - 1)+ constant (44) B s In accordance with the discussion in the preceding paragraph, it will be assumed that the heiht above the free liquid surface is depressed by the imposed pressure, The fluid will then have to assume a saturation distribution such that the gravitational and capillaric

-68forces are in equilibrium. Since at high saturations, the capillary pressure varies negligibly with saturation, it will be assumed that the height at which S = 1 is found will be depressed in proportion to the applied pressure. This depression can be easily calculated from the original capillary pressure curve. The constant can thus be evaluated. 2. Comparison of Theory with Experiment The data of Dombrowski (31), p.1216, are shown in Figures 20, 21 and 22. The solid lines represent the theory as given in Equation (42) and (44). The outflow face is at z = 0 and the saturation there is taken from the experimental data. The parameters of the curves are given as the quantity equations (42) and (44). The outflow face is at z = 0 and the saturation there is taken from the experimental data. The parameters of the curves are given as the quantity +k Ap lb-f (' =p+ L ft3 This ~' may be thought of as the combined desaturating potential consisting of gravitational and pressure factors. Since a higher pressure drop across a core will cause increased desaturation, higher t' values imply lower saturations. This is illustrated by the data. The experimental data points are rarely more than 20% from the theoretical. In most cases, the integrated average experimental saturations agree with theory to within about 10%. The curves for which the outflow face saturation is 1.00 correspond to the -situation prior to breakthrough.

28 ~'____II 24 THEORY, EQ. 42 AND 44 X 0 DATA OF DOMBROWSKI (31) 20 100/115 MESH GLASS E 0 12 6 0 12 1:i60.3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 09 1.0 FRACTION LIQUID SATURATION, SA Figure 20. Liquid Saturation Adjacent to Outflow Face of Porous Bed.

14 _~~~ 0 12 E 10 ~ ~\ \ ~ (\~- ~ THEORY, EQ. 42 AND 44 X 0 0X DATA OF DOMBROWSKI w () 0 U ^~~~~~~~~~~~~~~~~~~60/65 MESH GLASS 3:8 L1+ 0 L.L 0 o ______\_____ \______ ^^ ________ 0 w 6 LL. ~ y I 0 -89.4 Z 105.6 U)n 04 5=163.5 2 0 0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10 FRACTION LIQUID SATURATION, S Figure 21. Liquid Saturation Adjacent to Outflow Face of Porous Bed.

6' THEORY, EQ.42 AND 44 E \ X E 0 DATA OF DOMBROWSKI (31) C 4 l \ \, L 35/42 MESH GLASS 4 0..L0 116.5 0 \ z 2 \49 \X C = 198.0 ~~ ^~ ~ ~ ~____F L~ S,E 0 0.1 0.2 0.3 04 0.5 0.6 0.7 0.8 FRACTION LIQUID SATURATION,S' Figure 22. Liquid Saturation Adjacent to Outflow Face of Porous Bed.

-723. Discussion The application of this work to gas reservoir engineering is immediate, and based on practical field experience. It is well known that operators in the field report that the ability of a well to produce gas under a specified pressure drawdown is greatly reduced when a well has been watered out. The water around the well bore is usually very difficult to eliminate. Since the well bore represents a discontinuity in the porous body and hence in the capillary pressure, end effectswill be present. This work should assist in explaining the occurrence and persistence of water blocking at the well bore. The previous theory on capillary and effect suggests that the problem may have commercial significance for gas reservoirs. As preceding theoretical and experimental work has been seen to demonstrate, the increased saturation at the outflow face (well bore) may cause significant decreases in the permeability to gas. A hypothesis similar to this has been recently advanced by Ribe (100). He makes certain calculations showing that the performance of an oil reservoir blocked by water at the well bore is greatly different from that of a normal reservoir. It is suggested in this thesis that the saturation build-up due to end effect near the well bore may significantly affect the performance of commercial gas storage reservoirs. In order to test this hypothesis, the following model was set up. Pressures at the well bore and aquifer inner radius are assumed known. The radius of the well bore and the gas reservoir are known. A zone of known water saturation and radius is assumed to surround

-73the well bore. The system is at steady state with none of the radii changing significantly. Calculations will be made comparing the production rate from the reservoir described above to that of a reservoir which has no zone of saturation near the well bore. The following table shows the data used in making the calculations. TABLE III Characteristics of Assumed Water-Blocked Gas Fields Symbol Quantity Case A Case B rB Well bore radius (ft).300.300 rS Radius of saturated zone (ft) 30 3 Permeability in unsaturated..one k k Krg Permeability in saturated zone.1 K. 333 rb Reservoir radius (ft) 102 to 106 12 to 106 The geometry is considered to be radial. The theoretical solution of the problem is quite easy - being simply that of radial steady state flow through two layers of unequal permeability. This method of calculation is given in:McCabe and Smith (78), pp. 422-4. Relative permeability curves show that the water saturation in the saturated zone is 79% for case A and 40% for case B. It is estimated that cases A and B are roughly equivalent to the limiting cases of *Actual saturation distributions around well bores will, in general, not necessarily be the same as this assumed distribution.

-74commercial importance0 It is assumed that the "unsaturated zone," ioeo, from rS to rb, contains connate water. In general, the volumetric flow rate per foot of pay thickness is Q _ in (45) L R where R is the resistance to flow of the reservoir and is defined as (for the unsaturated reservoir) R (rB - rb) (46) KAL where _ = viscosity of gas, and AL = logarithmic mean areao It can be immediately seen that, provided the relative permeability curves do not change, the flow ratio between the blocked and unblocked reservoirs is independent of the absolute permeability and the viscosity. Calculations have been made using the data of Table III and the results presented in Figure 23. Results are shown as QB/QU vso rB where QB = flow rate under given total Zp in blocked reservoir Qu = flow rate under same total 4p in unblocked reservoir rb = radius of reservoir. The term blocked reservoir is taken to mean a reservoir having a zone of higher saturation around the well boreo The curves of Figure 23 show that the very small amount of water around the well bore will drastically change the ability of a gas storage reservoir to flow gas under a given pressure drop as compared to flow through dry sando In making the calculations, it was assumed that the water near the well bore remained immobile In the case that this water does move, the curves of Figure 20 give the initial productivity ratio0

1.01 1 -11 — 1 1~ — - -- 0.8 3 CASE B - -- 0 0.4 0 0.2 ~ CASE A =~- - 2 3 4 5 6 10 10 10 10 10 RESERVOIR RADIUS, re (ft) Figure 23. Productivity Ratios vs. Reservoir Radius for Various Cases of Water Blocking.

-76The computation of the liquid content of unconsolidated sands when capillary end effects are present is seen, from Equation (34), to involve only the following factors. 1. Minimum irreducible liquid saturation S'. rw 2. Relative permeability to the non-wetting phase. 3. Capillary pressure. 4. Absolute permeability. 5. Fluid densities. 6. Viscosity of non-wetting phase. 7. Volumetric velocity of displacing phase. The viscosity and relative permeability of the wetting phase are not involved, as might be guessed from the fact that this phase is not moving. In general, the saturation at the outflow face must also be known before the calculation can be made. The outflow face saturations used in the Figures 20 to 22 were taken from the curves of Dombrowski (31). In the absence of any theoretical method of predicting these outflow face saturations, it is suggested that they be obtained from empirical correlations. D. Countercurrent Capillary Imbibition A concept which is closely related to capillary pressure is that of drainage and imbibition. These two terms are usually reserved for phenomena in which only the capillaric and gravitational forces are acting. The capillary forces tend to imbibe the wetting phase into the

-77porous medium while the gravitational forces tend to drain the fluid out, or at least to keep it segregated at the bottom of the element. The drainage phenomenon results when the wetting phase fluid has been raised vertically to a level higher than it would rise by capillary action alone (cf. Fig. l(b) ). When the external forces are released, some of the wetting phase drains out until the fluids are at capillary equilibrium. Saturation distributions during a typical drainage experiment are shown in Figure 24. Terwilliger, et.al. (121), give a theoretical and experimental study of gravity drainage. Their analytical technique involved finite difference calculations. Also presented in this paper (121) is an account of the stabilized zone technique, which is an alternate method of calculating saturation distributions during dynamic multiphase flow. The stabilized zone technique was used in the calculation of the saturation distributions in Section IIB. The imbibition phenomenon results when the wetting phase has not yet reached the height to which it will normally rise by capillary action. This imbibition process may be either cocurrent or countercurrent depending on whether the top of the porous element is open or closed. Further information on drainage and imbibition may be found in (16), (30), (44), (46), (77), (115), and page 489 of (61). When a wetting phase is applied to one end of a sand pack having a low initial saturation, the wetting phase will spontaneously imbibe into the pack even though no external dynamic forces have been applied. The model to be considered is a long, horizontal linear porous medium filled initially with gas and connate water. One end of the porous element is closed and the other is suddenly placed in contact with water (wetting

-78t =O w w 0 CAPILLARY EQUILIBRIUM I I I I I..,, I 0 1.0 FRACTIONAL WATER SATURATION, S Figure 24. Height Above Base Plane vs. Fractional Water Saturation at Various Times for Typical Gravity Drainage Process.

-79phase). Water will spontaneously flow into the sand while the gas is expelled countercurrentlyo The practical importance of examining the problem of horizontal, countercurrent capillary imbibition lies in the fact that it may be an important method of production of petroleum and natural gas from underground reservoirs. In these reservoirs, there may be isolated regions of high permeability embedded in regions of low permeability. As the aquifer water rises by, say, bottom water drive, these relatively impermeable regions can be by-passed partially or totallyo The only way in which the gas or oil contained in these pockets can be produced is by countercurrent capillary imbibitiono Specifically, oil and gas can be trapped in long beds in the reservoir, which are open at only one end and otherwise surrounded by impermeable shale streaks, Under such condition, direct frontal displacement by water, of the kind described in Section IIB, is virtually impossible. Under these circumstances, the aquifer water will rise to the open face of the permeable oil or gas bearing zone by means of bottom water drive, thus allowing a countercurrent capillary imbibition to occur, Since this may be the primary production mechanism in many cases of practical importance, it is desirable to be able to determine the rate and efficiency of oil or gas production. 1o Theoretical Development The model to be treated is that in which a semi-infinite, horizontal porous medium, initially containing gas and connate water, is touched at its open end by an infinite supply of watero The effects of gravity are neglected and the linear porous element is assumed bounded by impermeable barriers at all points except the inflow faceo Water is

-80imbibed and gas is expelled countercurrently. Under conditions where the basic differential equation has non-linear coefficients, a technique analogous to that given in Section IIB can be applied. The total velocity ut at all points is zero by virtue of the assumption that the pressure drops are small relative to the total pressure (pseudo-incompressible flow) and that the porous medium is closed at one end. Because of the countercurrent nature of the flow, w = - Ug and Equation (7) yields ut = 0. Using ut = O0 neglecting gravity and making the normalizing definition, s= wr = wr (61) S' -S' S' o wr where S' is the saturation occurring at the inflow face of the porous medium. The governing equation is obtained from Equation (11) as,,- - - (- g-Fl ap) - ( 1,rgFl dPca) (62) at az p.g z az kg dS 6z Note that the saturation is normalized with respect to the saturation S' attained at the inlet face (z = 0) of the porous element. 0o The determination of the saturation distribution is to be made according to the following approximate solution.. A typical saturation profile at a specific time during the course of a countercurrent imbibition process is given in Figure 25. Two separate solutions, S1 for z < zo and SIi for z > zo are assumed, as shown in Figure 25. These solutions are then matched at the point of farthest water penetration

-811.00So so wI oI I[ z 0 I — w Figure 25. Fractional Wetting Phase Saturation vs. Distance During Typical Countercurrent Capillary Imbibition.

-82(Zo) by means of the boundary conditions. zo is a function of time only. The following boundary conditions are specified. 1. S (z, 0) = 0 (initial condition) 2. SI (O, t) = 1.00 3. S (z,t) = S (Zt) = 0 4 aSI(ot) - II(zt) = 6Z 6z ~2sI 5. ~(zot) = 0 az2 Boundary condition No. 3 expresses the continuity of saturation at zo. Boundary condition No. 4 expresses the continuity of flux at z. This may be seen by noting that the flux at any point is given by Equation (9) as, (ut = 0, gravity neglected), rgF1 dP aS 6S u rgl c= (63),w dS az haz 1rgFl dPc where the group (rgF ) may be thought of as a hydraulic conductivity [g dS Kh relating the flow rate to the saturation gradient. Boundary condition No. 4 is thus inferred from the fact that a saturation gradient is the only cause of flow (cf. Equation (63) ), and that the normalized saturation is always zero in region II (cf. Figure 25). Boundary condition No. 5 is derived as follows. Designating KK F1 dP r ( c) as some hydraulic diffusivity D (S), a form of the diffusion dS equation is obtained from Equation (62) as,

-83At z z, Equation (64) gives as (D - as a D asi (65) at Zo az az z z z zz Z 2 The right side of this last relation is simply D 2 because of boundary condition No 4, The left side of Equation (65) is zero by virtue of the total differential relation dS as + 6s dZo dt z t z dt z and boundary conditions No. 3 and No. 4o Boundary condition No. 5 is thus established. The development of boundary condition No 5 was suggested by Ko H. Coats. This additional boundary condition was apparently overlooked by Vetrov, et.al, (126) in a similar analysis. Since SI is zero for all z > z, the subscript I will be dropped from S, and a single solution will be treated. I The saturation is now expanded into a power series of the normalized distance -= L. S = a + al +a2. +.oo an n +... (66) Calculating the constants a to a3 in the (truncated) series from boundary conditions 2 to 5, there results S = (1 - )3; (67) The functional relationship between z0 and t is now established by means of the partial integral of the differential Equation (64)o

-84zo a - a f - dz = f (D -) dz o at o az a This becomes, using Leibnitz' rule, 0o d Zo dz - - f Sdz + S( t) zo - D dt o dt z substituting for S from Equation (67) gives 1 dzo 3Do 4 dt zo From which z = [24 Dt 2 (68) where Do and Dz indicate the values of the hydraulic diffusivity o (cf. Equation (63) ) at z = 0 and z = zo, respectively. The total production of the in-place fluid from the porous medium at time t is, zo Total production = J f S'dz - S'r z 0r o Substituting from Equation (61) and Equation (66) gives, - 1 1 3Dot ___ D_ ft5 Total production = S' ( 2 = 1.22 S (Dot)2 (69) ~~~2 ~ft2 (69) The rate of production at any time t is 1 1 sD3o - Do, ft5 Production rate = - ( )2 =.61 X AS' (-) 2 (70) 2 2t t ft hr.

-85The total (cumulative) production and the rate of production are given in actual cubic feet of non-wetting phase at the pressure of the operation. The depth of wetting phase penetration, saturation distribution, rate of production and cumulative production are given by Equations (67), (68), (69) and (70), respectively, with Do defined by, Do -Kr /1 dP\c^ (71) o'r L ds Jz=o The procedure in using the method is as follows. Do is calculated at S' (saturation at z =0) from Equation (71) using experimental determinations of relative permeability and capillary pressure. Some estimate of the saturation S' at z = 0 (inlet saturation) must be available, The functional relation between z and t is computed from Equation (68). Finally, cumulative production and rates of production are computed from Equations (69) and (70) and a knowledge of the cross sectional area perpendicular to flow. In the case of a porous medium having the special property that, bD = rg1 dPc = constant = Tc (c.fo Equation (62)) (72) ns I g dS the governing equation is linear, ite,, as 2S ~ = DC ~6^ (73) at c^

-86The solution for the normalized saturation is then, for a semiinfinite system, S - erfc \ (74) The analogous problem in heat transfer has been worked out previously (22), page 132 ff. The cumulative non-wetting phase efflux is given by, (cf. Equation (6) on page 133 of (22) ). Total production = ~ (4ct) / 1 (ct) (75) and the production rate is, Production rate = ~6 ( 1) = 5 (c )/ (76) 2 at 2. Discussion The question of the value of the saturation S' at the inflow o face (z = 0) of the linear element is important in capillary imbibition. The value of Do in Equation (68) is dependent on S', and this in turn 0 affects Equations (67), (69) and (70). There does not seem to be any way of calculating this inlet saturation a priori from presently known theory. Co-current capillary imbibition experiments (48) indicate that the inlet saturation will be somewhere near the residual equilibrium non-wetting phase saturation. In the absence of other information, a value of S equivalent to about 90% of the residual equilibrium nonwetting phase saturation is indicated,

-87It is seen that Equations (75) and (76) are the same as Equations (69) and (70) except for a small difference in the numerical factor. Handy (48) observed experimentally that cumulative production varied as the square root of time for co-current capillary imbibition in the case where gravitational effects are negligible. The solution given in Equations (67) and (68) is an approximate one in the same sense as the solution for the saturation distributions in Section IIB. Hence the argument supporting its validity will be the same as that given in the discussion of Section IIB. No data on countercurrent capillary imbibition have been found in the literature. However, the experimental determination that cumulative production is proportional to the square root of time (48) for co-current imbibition is encouraging, as is the square root relationship resulting from the linear diffusion model (cf. Equation (75) ). In order to test the ability of the approximate solution to give a correct prediction of saturation distribution, production rate and cumulative production, the following example calculation will be made. Assuming a constant hydraulic diffusivity Dc (cf. Equation (72) ), the exact expressions for saturation distribution, rate, and cumulative production will be computed from Equations (74), (75) and (76). The approximate solution, embodied in Equations (67), (68), (69) and (70) will then be used to compute the same quantities for comparison. 3. Example Calculation Natural gas of viscosity.01 cp. is assumed trapped in a thin, lamina-like permeable streak which is closed on all sides except for the face adjacent to the water flooded zone. The relative permeability to

-88gas and water for the porous medium have been measured to be Krg = 1- S and Krw = S (S is normalized saturation) The connate water content (S' ) and residual equilibrium saturation (S' ) have been measured as wr gr o10 and.90, respectively. Permeability and porosity are found to be 250 mde and 25%, respectivelyo Aquifer brine viscosity is 1i00 cpo The capillary pressure of the porous medium has been found to be P lb-f d c 100 lb-f P, = 100 (1 - In S) ft2; - = -S ft The hydraulic diffusivity Dc (cfo Equation (73) ) can be simplified by noting that, to a very good approximation, D KrgF1 Pc - Kw dPc 2 1-g idS - 1W;dS over the whole saturation range except for a very small region at very high water saturation. The saturation at the inlet face is taken to be St = 81 (or 90% of the residual equilibirum gas saturation) The hydraulic diffusivity is in this example a constant and is computed as D = 634 x 10 o The saturation distributions, rates and cumusec lative productions were computed by both the exact and the approximate techniques. The results are presented in Figures 26a, 26b and 26co The total cumulative gas production (Figure 26c) and the production rate (Figure 26b) according to the approximate solution are seen to be 7~8% higher than the same quantities predicted by the exact solutiono Figure 26a shows that the approximate solution not only agrees well with the form of the saturation distribution predicted by the exact method, but *Porous media will, in general, not have even approximately constant hydraulic diffusivities. The example chosen is, therefore, illustrative only.

0.9 0.8 0.7 k =250md, =0.25 \ \ ~EXACT SOLUTION Eq.(74) -c Q6 ~\ ~~~ ---— APPROXIMATE SOLUTION Eq.(67) z 0 Q5 c, i~ \ — \-: 0.4 i 0 25 50 75 100 125 150 175 200 225 250 275 300 z (f t ), DISTANCE FROM INFLOW FACE Figure 26(a). Distribution of Water and Gas During Linear Countercurrent Capillary Imbibition Into Porous Medium.

-I 10 __ ___________ -^ _______ ^',,_______ _ _ _... PRODUCTION RATE'C _1 PER SQ.FT. OF PRODUCING FACE ~^^b~~~'DURING LINEAR COUNTER CURRENT Z ______ ______ CAPILLARY IMBIBITION -- " k=250md =25% D ^ ~I | I~~ I~ I~ I I ( | |__EXACT SOLUTION (EQ. 74) 0 ---- APPROXIMATE SOLUTION (EQ.67) 0- 10 oI L........ ~ ~l~~ILI1 ~" C.!, < ~ ~~~~~ —-~ ~~I___~ ~~ __1 1 10 2 10 10 10 10 TIME,( DAYS) Figure 26(b). Rate of Production of Gas per Square Foot of Producing Face During Linear Countercurrent Capillary Imbibition Into Porous Medium.

, I I ii I iI I I I I I I I II i II t I I I I I ~~ TOTAEL PRODUCTION I 1 1 CL PR SQ. FT. OF PRODUCING FACE ~~ ~~ ~ ~ -'-,^^^~~- 1 1 1 l l l>DURING LINEAR COUNTER CURRENT.0.. < ~ Q ~~- CAPILLARY IMBIBITION O, k=250md ~25% " ~ ~ ~ ~~~~~~"~^ ^ ~ ~ ~~ | -~EXACT SOLUTION (EQ.74) ~' H ______ ___ ___- ^______ ___ ___ ~ l- -APPROXIMATE SOWUTION (EQ.67) 10 10 10 101 TIME (DAYS) Figure 26(c). Total Gas Production per Square Foot of Producing Face During Linear Countercurrent Capillary Imbibition Into Porous Medium.,.... Imbibition Into Porous Medium.

-92differs from it only by a few percent. The approximate solution thus predicts saturation distribution, total production and production rate with surprising accuracy in spite of its non-rigorous character. Summary (Sections II A-D) The preceding material serves to indicate the theoretical considerations involved in the scaling of laboratory flow experiments. The scaling parameters were derived from the basic partial differential equations governing the flow of two immiscible fluids in porous media. The results show that there are three dimensionless functions which must be scaled in two and three dimensional experiments and two functions for linear experiments. Experimental saturation distribution data were taken for the case of linear constant rate displacement of gas by water. An approximate solution to the equations of multiphase movement yielded good agreement with these data and with others taken from the literature. A theoretical solution to the problem of capillary end effect was obtained. This solution is applicable when the wetting phase is immobile and only the non-wetting phase is moving. Comparison of this theory with data from the literature yielded satisfactory agreement. An approximate solution to the problem of countercurrent capillary imbibition into a linear medium was obtained by means of an integral technique.

III. MACROSCOPIC STABILITY OF INTERFACES DURING DISPLACEMENT PROCESSES IN POROUS MEDIA A, Stability of Plane Interfaces in Linear, Quasi-Steady Flow-Effects of Inertial. Damping When a less viscous fluid, such as gas, displaces a more viscous fluid, such as oil or water, from a porous material, long fingers, of macroscopic dimensions relative to the pore size, may form under certain conditions. These fingers may bypass large amounts of the in-place fluid, with resulting technical and economic disadvantage. In natural gas storage reservoir development, the water (aquifer) which is initially found saturating the porous structure is displaced by the injection of high pressure gaso The shape of the interface between the gas and the water, and the rate at which gas can be injected are important considerationso In particular, if the gas moves into the aquifer unstably, there is the possibility that some of the long fingers may be pinched off and become trapped as permanently unrecoverable gas, The study of viscous * fingering is, therefore, expected to give insight into the details of the fluid movements involved in fluid-fluid displacement, and particularly into the mechanics of the natural gas reservoir development process, The recovery of oil by frontal displacement is also an example of a physical situation in which fingering can occur. The less viscous flood water may enter the porous medium (initially saturated with oil. and connate water) in the form of fingers, and break through to the producing wells before any significant amount of oil has been recovered, It is desired to be able to predict the occurrence and the form. of the viscous fingering from general. considerations relating to the fluid and *The term "viscous fingering" was first used by Engelberts and Klinkenberg (36) and refers to the unstable fingering resulting from the displacement of a fluid by one of lesser viscosity. -93

-94porous medium properties, such. as porosity, permeability, interfacial tension, densities and viscosities, etc, The physical model to be assumed in the following treatment is that in which one fluid is being linearly displaced from, an isotropic, homogeneous porous medium at constant volumetric velocity by another immiscible fluid which is assumed separated from the first by a discrete interfaceo Saturation distributions, as treated previously in Section II, are not considered. The porous medium is assumed 100% saturated with the displacing fluid up to the separating interface and 100% saturated with the displaced fluid beyond the interface The fluids are immiscibleo It is worthwhile to distinguish between the viscous fingering described above, and some of the other phenomena which can give rise to unequal fluid flow in porous media. They are1. Finger formation due to permeability stratification. 2. Streamers of displacing fluid resulting from, the potential field - as in the reaching out of water toward the corner wells in a five-spot flooding arrangement' or in the upward coning of water toward a single producing well. 3o Unequal flows due to gravity segregation of fluidso 4o Mechanical instability which can occur when a heavier fluid is temporarily found above a lighter oneo The kind of instability to be considered here must be carefully distin. guished from unequalflows caused by the items listed above, The theory and phenomenology of viscous fingering (hydrodynamic instability) to be considered here has been discussed by Chuoke etoalo (23), Saffman and Taylor (105), and Van Meurs (125)o Photographs of the processes are

-95given in these references, and one from Van Meurs (125), is reproduced here as Figure 2o [Other papers on this subject are (20), (21), (29), (36), (104), (108), and (124)o] The analysis employed in the papers mentioned involved the use of complex variables in the manner of Rayleigh and Kelvin, i.e., first order perturbation theory. The final step in the development was to equate the pressures on the interface between the two fluids, These pressures were obtained by straightforward computation from Darcy's Law, with interfacial tension taken into account. Also, these analyses (23) and (105) apply only in the case of fluids moving at a sufficiently low velocity. It is to be demonstrated in the following work, that the nonsteady and inertial effects involved in the finger acceleration may both affect the flow, and that both of these may be included in the theoretical analysis. The kind of instability to be treated in the whole of Section III is designated by the general term "hydrodynamic instability" and is to be carefully distinguished from the other forms of non-uniform behavior listed above. The term hydrodynamic instability is used, since the fingers are due solely to the fluid dynamic effects and are not conditioned by non-uniformities in either the porous medium or the applied forces, The work of this section differs from that of Chuoke etoal. (23), and Saffman and Taylor (105) in that inertial and unsteady effects were taken into account, The analysis presented here yields the result that the initial growth rate of the instabilities is less than previously calculated in (23) and (105).

-961. Theoretical Development The physical model to be considered is that shown in Figure 27, Fluid 1 of p1 and p1 is displacing fluid 2 of P2 and 12 at volumetric velocity V. The fluids are immiscible and separated by a discrete interface on either side of which are no saturation gradients. The porous medium is assumed homogeneous and isotropic. The first step in the solution is to formulate a complete expression for the pressure in a porous material under unsteady conditions Consider the flow of a viscous, incompressible fluid through the differential volume element shown in Figure 28. The velocities to be used will be the integrated average (Darcy) velocities, which are related to the true point velocities in the volume element by ui _ 1 uid~ (77) where i signifies any principal direction, and ui is the microscopic point velocity within the pore. J indicates integration over the void space (pore volume) b6x6ybz of the total volume element 6x6ybz. The total u. volume element bxbybz contains both fluid and solid. The quantity is recognized as the average pore velocity. The pore spaces are considered to be completely filled with fluid, i.e., at 100% saturation. The velocity u. is known as the volumetric velocity (fluid volume per unit area per unit time), and is a statistical average. It is assumed that the porous medium is so constituted that the small differential volume elements used in the analysis do not approach the pore size. It is only under the preceding limitations that it is possible to speak of

-97P2, 2 FLUID 2 Lx-'1~~~~~~ -x FLUID I PI. 11 V Figure 27. Schematic Diagram of Linear Displaeement of Fluid 2 by Fluid 1 at Start of Fingering Process.

-98yy lv w x x zFi Figure 28. Differential Element of Porous Media.

'99' a point volumetric velocity uio One of the basic premises of flow through porous media is that the volumetric velocity is derivable from a potential, U. - V (78) The result is experimentally based, and has its origins in the work of Darcy (28), It has only been verified in the case of steady flow. As has been shown by Streeter (p. 234 of (118) ), microscopic inertial effects in the steady state flow of a homogeneous fluid through porous media can be shown to be negligible if Equation (78) holdso In addition, when the fluid is incompressible, V ui = 0 (79) and the Laplace equation holds 2 O = O (80) The velocity potential D has the same units as the pressure (the force potential concept of fluid motion is implicit in Equation (78)). In the steady state, the force exerted by the walls of the pores on the fluid contained in the pore space of the volume element of Figure 28 is calculated as Fi - 0 6xbybz (81) This Darcy viscous reaction is calculated in the x direction, say, by multiplying the potential drop - by 5x and the area of application 86ybz. The quantity FD can also be thought of as the momentum flux 1

-100between the fluid and the walls of the pores contacted by that fluid. The admission of negligible inertial effects on a microscopic scale, i,e,. use of Equation (78), by no means justifies elimination of the macroscopic inertial effects. The distinction between the approaches is that on a microscopic scale, the actual processes occurring inside individual pores are examined whereas, in the macroscopic approach, the processes are treated as the result of averaging over a very large number of pores. Microscopic inertial effects are related to the acceleration and deceleration of fluid particles as they pursue tortuous paths through the pores. Equation (78) implies that these effects are negligibleo Macroscopic inertial effects are related to gradients in, and time rates of change of, the volumetric velocities. To clarify this point, it must be noticed that the width of unstable fingers is usually many orders of magnitude larger than the average pore diameter. Equation (78) implies that the microscopic (point).inertial effects are negligible But the fingers, being of much greater than microscopic dimensions, are thought to generate appreciable inertial effects on a macroscopic scale due to their acceleration, The situation is analogous to that in which a rigid bar is accelerated. While there is no acceleration of the internal particles of the bar relative to one another (microscopic scale), the acceleration of the bar as a whole will give rise to inertial effects on a macroscopic scale. From the point of view of the pressure in porous systems, it may be stated that microscopic inertial effects are negligible because there is essentially complete pressure recovery as a fluid particle accelerates and then decelerates to its initial velocity if the gross flow is steady. These microscopic effects are averaged to zero nAll topics in this thesis are treated from a macroscopic viewpoint unless otherwise stated.

-101over a large number of pores due to the isotropic statistical properties of the porous material. When the total flow is being accelerated, the gross inertial effect must, of course, be the sum of all the microscopic inertial effects. The continuity of this gross acceleration is what causes the macroscopic inertial effects to influence the pressure. In the light of the preceding discussion, we now turn to a formulation of the momentum balance for fluids moving in the unsteady state in porous materials. The balance is taken.on the fluid contained in the volume element of Figure 28. In addition to the body force, normal pressure force, and stresses, the viscous reaction of the pore walls on the fluid is considered. These forces are equivalent to impulse per unit.time and are equated to the net momentum flux for the fluid element. At any given time, the influx of momentum into the differential element in the u direction of Figure 28 is, at the point x, y, z. Influx= = (pu) 5yz + (pu) XSz + H (pu)Bxy (82) These momentum flux terms were found by multiplying the average pore velocity, e.g., u/0 by the mass flow rate e.g., pu by the cross sectional area, e.g., bybz..across which the fluid was flowing. In formulating these equations, it is important to keep in mind that the fluid volume is 5x6y65z, that its mass is pO5xSy5z, and that the area normal to flow in the x direction, say is fSy5z. Also the summation notation with i and j is used throughout. Here the subscript i is the free index, indicating which component is being considered, and j is the dummy index indicating repeated operations. The momentum balance is being made on a control

-102volume comprising the fluid in the pore space of the differential element of Figure 28. The outflux of momentum in the u direction at the point x?. y., z is; (u+ x 5x) Outflux = () (pu+ (Pu 5x) y5z (83) 6v (+ T ~ y) (pu + ( y) 6xz +. (pu + y) W 6y (w+ z) + ~ (pu +(pu) z) P5x5y The time rate of change of momentum is a.(p xy Lt To show that the Darcy velocity, ui can actually represent a physically real momentum, consider a very small volume element dV (much smaller than the pore size) in which the average velocity is ui. The momentum of the fluid in this volume element is, (in the i direction) puidV. The total momentum in a differential (macroscopic) element is then, from Equation (77) J puid~ = puij where ~ = 56xy6z The developments leading to the momentum influx and efflux terms are thus seen to be correct since thevolumetric velocity pu has been shown to correctly represent a momentum term.

-103The next step in the development of the unsteady momentum balance for porous media is to consider the shear stresses which may 6u. be set up due to velocity gradients —. Graton and Fraser (47) and axj Coberly andMarshall (24) have observed velocity gradients in flow through tubes packed with porous materials. These velocity gradients suggest the presence of macroscopic internal stresses and, in particular, shear stresses. The relation between the velocity gradients and the shear stresses which may occur is not known with any certainty. Also, the velocity distributions observed in flow through porous media in tubes may be due to inhomogeneities in the packing. The flatness of the velocity profiles obtained experimentally (24), indicates that the constant relating the shear stress to the velocity gradient in porous media is much smaller than the corresponding constant (viscosity) in homogeneous fluids. The stress-deformation relation for a fluid of viscosity i. flowing through a porous medium of permeability K will be assumed to be represented by the (Newtonian) relation.i = (. + __) (84) ~~ wherexj where Tij = the stress in the j direction of the face perpendicular to the i direction, and f1 = a pseudo viscosity which is expected to be a function of \1 and K. The forces exerted on the fluid at x and x + bx on a face 6ybz perpendicular to these directions are po6ysz and p - dx)#6ybz, leaving ox

a net momentum interchange of a- o x6y6z (85) ax The body force on the fluid mass pb5x6ybz is Xp06xy5 z (86) Formulation of the shear and normal stresses as on p. 216-218 of Streeter (118) yields, for the momentum flux due to shear stress in the x direction on the face normal to the y axis -TX 0 bxbybz (87) by The last force acting on the fluid in the differential element is the drag from the walls of the pores surrounding that element. This is the Darcy viscous reaction. It will be designated F', with the prime Di denoting the unsteady state. The associated momentum interchange will be F' Di ~ bxCybz (88) The momentum balance on the fluid in the element of Figure 28 may be expressed verbally by saying that the net flux of momentum out of the element in the x direction is equal to the negative time rate of change of the momentum content of the fluid plus the net momentum interchange occasioned by.the forces acting on the fluid. The balance is formulated from Equations (82) through (88) with second order terms neglected as,

-105Outflux- Influx - time rate of change + external forces, or, I Fau au au au- FD 1 up ~ 1 a + + ~ V W + X + XaLt ax 6y azJap p x -x p ox (89) + _X + zx p by p az and likewise for the other principal directions. Proceeding according to Streeter (118) pp. 217-22c, the desired equations are obtained as +1 ari au i - i = V Iopi + xi (90) where v' = = pseudo kinematic viscosity. P The form of the momentum balance Equation (90) can be clarified by recognizing that the left side of Equation (90) represents the total momentum change per unit volume per unit time and the right side represents the corresponding unit impulse. When macroscopic velocity gradients are absent and the flow is steady (i.e., ui and a__ are zero) Equation (90) is seen to reduce at axj to Darcy's law, as required. (Note: under these assumptions FD = FDi = 6xSyz and 2ui = 0). The quantity F' represents the force exerted on the fluid by the walls of the pores, The Equation (90) is the momentum-impulse equation for porous materials and is analogous to the Navier-Stokes equations for homogeneous fluids. In fact, the Equation (90) represents spatial mean macroscopic Navier-Stokes equations for porous media in much the same way that time-mean equations are developed to represent the macroscopic aspects of turbulent flow,

-106It is now possible to proceed to the complete formulation of the momentum balance for the unsteady case. The object of this development is to obtain results which can be applied to unstable finger motion, which is a special example of unsteady flow. The integration of the basic momentum equations depends on the possibility of finding a suitable expression for F. (the force exerted by pore walls on fluid in the unsteady state). The steady state value (FDi) is given in Equation (81). Philip (91) has recently shown that, unless varying potentials of large magnitude are very suddenly applied; or are periodic with very small period, the microscopic velocity field set up in the pores of a porous medium are functions only of the instantaneously applied potential, and not of the past history of the applied potential. This means that, as long as attention is confined to a single phase flowing in a porous medium, the volumetric velocity ui at any point is a function only of the potential gradient at that point, and not of the unsteady character of the flow. The result of Philip's (91) dem onstra-tion is that the velocity u. and the velocity field in ui do not lag behind the applied potential except under conditions which almost never happen in real displacement processes. Since the microscopic velocity fields are set up instantaneously, it is reasonable to suppose that the drag at the pore walls will also be instantaneously a function of applied potential. It is therefore concluded that the pore wall drag is, within limits, the same in the.unsteady state as in the steady state, i.e., pu ~Fi gx~y~ z A(91) D. D. 1 1

-107The physical content of Equation (91) is that the pore wall drag is a function of the instantaneous volumetric velocity only, and is independent of the unsteady nature of the flow. From Equation (78), it follows that the flow is irrotational, ioe., Vu. = 0. It is also being assumed that the pseudo viscosity CL' is zero, Substituting this last relation and Equation (91) into Equation (90) yields the simplified form of the impulse momentum equation. for porous media, 1 6 i 6ui1 I 6P - ui i - + uj ax = - ~i + Xi - (92) L 6t J xj pP x - pK Integrating these equations in the manner of Streeter (118), po 23 gives, for fixed volume element, 2.1 l( _O + gz + p+ C = 0 (93) 2 X at PK P where q = - V 0 and the bars indicate quantities taken with respect to fixed volume elements C = arbitrary constant. Equation (93) can be called the Bernoulli-Darcy equation, and is the result which is desired for purposes of future applications. It is seen to include the usual terms of the Bernoulli equation plus the effect of of the Darcy viscous reactionso

-108The stability of the system illustrated in Figure 27 is now computed. Formulations of the simpler forms of the stability equations have been done by Chuoke et.al.,(23), and Saffman and Taylor (105). A summary only will be given here. Fluid 1, of viscosity K1 is assumed to be displacing fluid 2, of viscosity 12 (immiscible fluids), with velocity V. Referring to Figure 27, assume that the co-ordinate axis is moving along in the position of the undisturbed interface with velocity Vo The governing equations are, v2 The (three dimensional) perturbation is, nt k k b = Be cosx cosk-y (94) 12 12 i.e., periodic in x and y. Positive n will indicate an unstable system. 2. Boundary Conditions 1. Perturbations vanish as z -+ oo 2. Pressures equal on interface except for surface tension. ab 3 =w = w on z = b at 1 2 where primes refer to the perturbation velocity only, and ) refers to the velocity potential of the perturbation only, since we have eliminated the uniform flow by the expedient of moving with the velocity Vo The solutions are, nt kz k k t — A e e osx cos y (95a) nt kz k k = A e e cos —x cos — y (9b) 1,2 12

-109where nB A = k Assumptions are that the perturbations are small and that the quadratic terms of the substantial derivative can be neglected in developing boundary condition 3. (c.f. p. 189 of Chuoke et.al. (23) ). In order to convert the Bernoulli-Darcy Equation (93) to the moving coordinate system, the following transformations are used, = - Vz 2 (6)2 ()2+ (z + v)2 6x 6y 6z These transformations are used in Equation (93) to give 1 pq' pw'V 1 pV p _ _Vz _ = -=- - p~ -- - ----, -..... +- -- pgz + C (96 2 20 t K K where,2 ()2 + (_s)2 + ()2 2 q- ==total velocity of perturbation only. Suppose that the value of the pressure on the interface (z = 0) is po in the absence of the perturbation. Then C = p + 1/2 P and Equation (96) becomes p =+- pw +e_ a,V + I- pgz + po (97) 2 e ert K K All quantities without bars refer to the perturbation only, which is the center of interest,

-110The expression for the pressure in Equation (97) is now linearized by recognizing that q will be much smaller than V and that 1 pq. pwL'-V the term 2 will be negligible compared to. The term in q is deleted from Equation (97) based on the assumption that the second order terms of the perturbation are insignificant. The term in w'V represents the inertial cross-linking effect. The pressure balance at the interface is, cf. Chuoke, et.al. (23) Pi - P2 = (c + c2) (98) where c and c2 = interfacial curvatures with relation to phases 1 and 2, respectively. Making use of Equations (94) and (95) in Equation (98), and the fact that ekz = ekb = 1 on z = b, there results, n2(plp2) (P2-Pl)nVk (A2 C1 r^2P,,t + n ( ) - -- 1)V+(p2-Pl)g k+ak3= 0 ^ 9K2 Ki LK2 K1 (99) The validity of Equations (97) and (99) may be inferred by a comparison with the corresponding equations derived by Chuoke, etoal, (23) in their simpler analysis. Equation (97) is the same as their Equation (4) except for the terms 1 pq'2 _ +V aP 2 6 - t Equation (99) is the same as their equation (12) except for the terms n^, (pa -pl)nVk n- (P + P2)' (2.-)n which result from the terms mentioned above,

-111Solving for n from Equation (99), I 2 "i' (P2-Pl)Vk]2 - (1.2 111 )- (P2 - Pl)Vk] + 4(+ (P+P2) ( l)V(pp2 K;2 Ki (P2 +) + K1 )VMpPP2) (-)(p^) % n = 2 (P1 + P2) ~~ (:too)(100) Instability occurs only when n is positive, and this will happen only when (2 - 1) V + (p2 - P)g- k2 > 0 (101) K2 I1 This turns out to be the same stability criterion previously derived by Saffman and Taylor (105) and Chuoke, et.al. (23). The actual value for n is different however (cf. Equations (12) to (15) of Chuoke, et.al. (23)). Chuoke, et.al. (23) calculate [2 "1 v + (P2- l)gk - ] n=~. (102) (-;. + 1) K2 K The critical velocity and wavelength of critical stability also turn out to be the same as previously derived (23), (105), i.e.,

-112X = 2n 2 and m 3 Xc (103) I92'l K2 Kl where Vc is defined by ^2 C1g ( l- ) Vc + (P2 - Pl) g -rk= O (104a) 2'l C The wavelength of maximum instability is found from Equation (100) by 6n solving = 0. The result is 6k (104b) 2(p p1)2V2 4(p2-P1)V 4(124a(pl+P2) 4(Pl+P2) L[2 1l 2 ^\ ~^~ #4 +^ f { +(2 1)v 2(P ~.... +(p2-Pl)g + 2 (-+ ) 2] Km * - 24(P + P2) a and 2 2\T2 em " k m This wavelength of maximum instability is not the same as that found by Chuoke, et.al. (23), who computed Note: For explanation of the meaning of X or Xc, see footnote at bottom of p. 190 of paper by Chuoke, etal. (23).

-113(-2 +I) (p - P ) g km: (104c) af3 It is seen that deletion of the terms in (P2-Pl) of Equation (104b) reduces Equation (104b) to Equation (104c). These terms in (P21) arise from the retention of the cross-linked intertial term (of form PW~) in the Bernoulli-Darcy equation, Equation (97). Thus, the specific difference between the wave number of maximum instability (k ) as computed in this thesis (Equation (104b) ) and by Chuoke, et.al. (23) results from inertial effects only. This means that the present analysis does not indicate a change in the naturally occuring finger width (X ) when the term w is not significant with respect to the others in Equation (97). The unsteady term - a results in a change in the instability index (n of Equation (100) ) but does not affect.m3. Discussion The previous solutions (23), (105) did not include the effect of the time rate of change of the velocity potential (-) or the cross6t linked interial'effects (Pw ) in the analysis. For brevity, the change in the value of n (Equation (100) ) resulting from the inclusion of the a6s at term in Equation (97) will be called "unsteady damping." The change in n resulting from inclusion of p- in Equation (97) is designated "inertial effect."

-114While it may happen that the simpler derivations in (23) and (105) will hold for certain cases, it is seen that the present development, taking proper account of the total momentum balance (Equation (90) ), is more complete. The expressions for n, (Equation (101) ), and the wavelength of maximum instability, X, (naturally occurring wavelength) are the onlybases on which it is possible to choose between the two derivations, and this must, apparently, be done experimentally. It is suggested that critical experiments measuring the rate of early growth of the instabilities, and the tip to tip distance (X ) of the fingers may m decide the matter, since these quantities represent the only points at which this derivation differs from the others. The present analysis ultimately rests on two assumptions. One is that the flow is potential - a reasonable one for porous media. The other assumption is that the pore wall drag on the fluid, represented by the term F' in Equation (91) depends only on the instantaneous velocity Di ui and not on the unsteady character of the flow. This second assumption can be clarified by an analogy. A successful method of computing the trajectory of a sphere starting from rest and accelerating to terminal velocity in a fluid is to assume that the drag coefficient at any instantaneous velocity is the same as the steady state drag coefficient at that velocity. Of course, the fact that the sphere is accelerating through that velocity rather than being constantly at that velocity may have some effect on the drag, depending on the magnitude of the unsteady effects. The effect, however, is small. In the same way, the past history of the fluid acceleration through the porous medium is here

-115assumed to have negligible effect on the drag which is exerted at a specified volumetric velocity. The demonstration of the existence of macroscopic shear stresses being produced on fluids contained in porous media by the action of Darcy velocity gradients has yet to be achieved conclusively. Also, the form of the relationship, i.e., whether Newtonian or not (macroscopic scale) (cfo Equation (84) ) must be determined, The phenomenon is complicated because of the action of these shearing forces on the pore surfaces as well as on the fluid. It is thought that only unconsolidated materials will be able to allow any significant amount of shear stress to occur, since it is only in these materials that the fluid is continuous to any great extent, The continuous nature of the solid and the associated interconnecting channel structure of the pore space seems to preclude macroscopic shear stresses from being generated in consolidated materials. The shear stresses referred to, and defined by Equation (84) are to be thought of as macroscopic shear stresses generated by gradients in the (macroscopic) Darcy velocities. They represent the summation of the large number of microscopic shear stresses occurring in the particular fluid element under consideration, The gradient, for example, is ay defined as, au = i 1a d (105) 6y F e ay These microscopic shear stresses are generated at all times, whether or not there are macroscopic velocity gradients occurring,

-116This distinction between microscopic (pore scale) and macroscopic shear stress may be clarified by comparison to an analogous situation which occurs in turbulence. The Prandtl approach to turbulence (cf. 298303 and 368-371 of Rouse (103) ) leads to the conclusion that the mean macroscopic shear stress is equal to the eddy viscosity multiplied by the gradient of the mean velocity. The mean velocity in turbulence is the time mean of the microscopic fluctuations. These small fluctuations also lead to microscopic velocity gradients which must lead to microscopic shearing stresses. The Prandtl approach thus leads to a mean macroscopic, integrated average description of the turbulent phenomena. The approach to flow in porous media which has been made in this thesis is analogous to the Prandtl approach to turbulence except that the mean quantities are space-mean, rather than time-mean, The microscopic scale in the porous medium is the pore scale, whereas the microscopic scale in the Prandtl approach is defined by the amplitude of a typical turbulent fluctuation (Prandtl mixing length). Thus, in both turbulence and flow in porous media, the macroscopic quantities are summations of the directed character of the microscopic quantities which always underlie the basic flow. Since there are, at present, no values of p' available in the literature, an experiment by means of which it might be determined is illustrated in Figure 29. The plate is drawn along the surface of the liquid which is saturating the porous medium and also extending a certain distance above it, Anemometers placed at the downstream end of the model are used to measure the velocity distribution in both the clear section and the packed section..' is determined by a comparison with the theoretical solution for a plate being drawn over the surface of a homogeneous

-117PLATE CLEAR LIQUID / \ POROUS MEDIUM VELOCITY SATURATED WITH I GRADIENTS SAME LIQUID Figure 29. Experiment to Determine the Pseudo-Viscosity of a Porous Medium Saturated with Liquid of Viscosity A.

-l18fluid of viscosity resting above a homogeneous fluid of viscosity which is immiscible with it. The exact nature of the stress-strain relation in porous media would have to be determined by fitting the experimental data to various Newtonian or non-Newtonian theoretical models. To examine the conditions under which the value of n in Equation (100) may differ significantly from the n of the work in (23) and (105), the following example calculation is given. 4. Example Calculation An oil of 10 cp viscosity is being displaced from a porous ft 3 medium (K = 100 Darcies) by water (p = 1 cp) at a rate of 1.0 ft.5 ft. sec. Assume interfacial tension = 20 cm and 2 =.25. Density of the oil is lb 37.6 -. Compare the index of instability n for the simple case of Equaft tion (102) and the more complex case of Equation (100). Chuoke, et.al. (23) have shown that the necessary condition for instability in the simpler system considered by them is that the velocity V be greater than the critical velocity V defined by V (P Pl) (106) v ((o6) c 112 41 K K 2 1 provided wavelengths greater than X are present (cf. Equation (103) ). Reference to Equation (101) shows that the critical velocity for the more complicated case considered here is the same as given for the simpler case. For both cases, V is calculated to be 1.28 x 10- -, assuming c sec. K2 =K = 100 Darcies.

-119The wavelength of maximum instability is computed from the following equation, for both cases,* (cf. Equation (15) of Chuoke, et. alo (23) ) Am 2X j3 2 = 00927 ft 2 Il) (V,-V) 2 1 ~ok - 6780 ft-l The valuesof n calculated from Equations (102) and (100) are 3560 sec and 3360 sec1, respectively. Physically, this means that the initial rate of propagation of the fingers is smaller when the more complex analysis is taken into account (i.e., the fingers are being damped). The inclusion of time rate of change of the velocity potential, as derived from the complete specification of the inertial effects, thus tends to damp the initial growth rate. The criteria for stability and critical wavelengths are unchanged. Examination of Equation (100) reveals that the damping effect (shown by the lower n from Equation (100) will, in general, be greater as the sum of the fluid densities and the displacement velocity become greater. The quantitative criterion establishing the range of variables over which these inertial damping effects are important is determined as follows. The expression for n in Equation (100) is expanded in series by means of the expression given in Pierce (94) No. 817. The result is Inertial damping negligible,

-12042 c1 (p2 Pl)Vk K2 Kl.,.. j (icy) (P2 8 - 16 (107) 2(P1 + P2) where 1(p+ {1) 2 + (p2 pP)g k-ak - L2+~1 J^ (p2-Pl)Vk L 2 K, It is seen that the expression in Equation(107) reduces to that of Equation (102) (i.e., simplified model) whena — is small relative to unity. (i.e., quadratic terms and higher are negligible). The criterion may be set up such that the value ofcsatisfies the inequality, a <.10 (108) To give the maximum difference between the two values of the amplification factor n as 2.5%. In the main difference between the work presented here and that previously done (23) and (105), is that the initial rate of growth of the fingers, for a given situation, is predicted to be smaller., The quantitative measure of this difference is given by the respective expressions for n in Equations (100) and (102). The criterion by which it can be estimated whether or not the morecomplex solution differs from the simpler one in any given problem is given in Equation (108). Examination of Equation (107) shows that the term (p + p2), arising from unsteady (at) damping always stabilizes the flow, i.e., makes n smaller. The term (p2 -P9) stabilizes only when p2 (displaced

-121fluid density) is larger than p1, and makes the system more unstable otherwise. The assumption that the flow is potential is basic to the analysis. Therefore, some additional justification of its validity is considered necessary. The formation of a rotational boundary layer is usually associated with flow parallel to a solid boundary. For small perturbation amplitudes, the flow considered here is predominantly normal to the surface of separation of the two fluids. In the treatment of unsteady boundary layers in homogeneous fluids, the first approximation to the flow for small times is usually taken to be the potential solution. This is illustrated for the case of objects started impulsively from rest in Schlichting (109), page 207~ Immediately after the start of the motion, the flow in the whole fluid space is irrotational with the exception, of a very thin layer near the solid boundary. This comes about because, in the very early stages, not enough time has elapsed to allow the formation of any significant boundary layer. It can easily be shown that irrotationality is usually formed by diffusion of vorticity from solid boundaries at which there is no slipo It is one of the basic postulates of boundary layer theory that viscous flows are irrotational at sufficient distances from solid walls. In the problem solved in this section, there are no solid boundaries present. The fluids extend infinitely in the directions perpendicular to the flow, and semi-infinitely in the directions normal to the separating interfaceo There is thus no solid boundary which can cause vorticity diffusion, Since the viscous fingering problem is characterized as flow which is essentially normal to the interface, in its early stages,

-122and without the influence of solid boundaries, the potential flow assumption is indicated. In any case, the initial flow can always be assumed to be potential, since the form of the perturbation can be arbitrary. This amounts to imposing a potential flow perturbation on Equation (92) and computing the resulting stability, Experimental data for flow in homogeneous fluids can also be advanced to support the assumption of an irrotational perturbation, Taylor (120), in computing the stability of interfaces between homogeneous liquids accelerated normal to their planes, assumed that the flow was potential even though the analysis was supposed to apply to real liquids. Lewis (76), in a subsequent verification of this theory, found that the instabilities were propagated in accordance with potential theory as long as the height of the perturbations did not exceed 40% of the distance between fingers. Air-water, water-carbon tetrachloride, benzene-water and air-glycerine pairs were accelerated at up to about 50 g in a channel of 1" 1" cross section 2 - 2 by. Since the resulting fingers were observed to be three dimensional, it could be concluded that the flow is potential in the initial stages even though the fluids were viscous. Since flows in porous media are considerably slower than those observed by Lewis (76), it is presumably safe to conclude that the initial stages of instability in porous media will also be irrotational It remains to demonstrate the persistence of irrotationality in the flow field if the fluid is irrotational in the state of rest. The equations governing vorticity diffusion have been developed on page 210 of Rouse (103) for a homogeneous fluid. Retaining the Laplacian

-123in Equation (90) and performing the appropriate manipulations, the equations for diffusion of macroscopic vorticity in porous media can be derived. For example, if - u ~v = - - -L --- Darcy vorticity (109) az x by then, - 6 - v - ( 11 Dt X x Y + y Zy+ v z P (1o) likewise for the other vorticity components: If the pseudo-kinematic viscosity v' is zero, then the right-hand side of Equation (110) is zero in regions of the flow field where flow is irrotational, As long as no non-slip conditions are imposed in the region in question, there is no way that vorticity can be introducedo In the problem of viscous fingering, the fluid is irrotational (uniform flow) far behind the separating interface. Since the substantial derivative (moving with the fluid) of the vorticity is exactly zero far behind the interface, there is no way that the vorticity itself can change, and the total. flow field must be potential at all times and places. Finally, it is to be shown that no net vorticity is diffused into the macroscopic flow from the microscopic velocity gradients~ The total vorticity Xo in a small fluid element 65xSyz (containing many pores) must z be the sum. of all of the microscopic vorticitieso 1 au av z (y f ( x)~U ~

-124But, from Equation (109) ff., -z must be zero. This means that the net amount of spin of the fluid particles in the pores in element bxy6bz is zero. Since the microscopic vorticities are certainly not zero at all points, the conclusion can be drawn that the statistical properties of the porous media are such that w is zero. That this is physically reasonable can be shown by the following example. In the very slow flow of a viscous fluid past a sphere, the vorticity is, in general, finite at any given point. However, from considerations of symmetry, the total vorticity in the flow field is zero. A given fluid spin at a point on one side of a plane of symmetry is balanced by an equal and opposite spin at the point's mirror image. Since the porous media can be visualized as an isotropic assembly of particles it is to be expected that the total vorticity in a fluid element will be zero by virtue of the symmetry of statistical properties. The results of Equation (110) ff. are that the flow will always remain potential if it is potential far behind the interface, and that no vorticity can be diffused into the flow field from the pore walls. B. Viscous Damping of Unstable Fingers The previous section treated the problem of the effect of unsteady and inertial damping of instabilities generated in displacement processes in porous media. In this section, the problem of the effect of viscous damping on finger formation and growth will be considered. The developments of Chuoke, et.al. (23), Saffman and Taylor (105), and the previous section of this thesis consider stability from the point of view of the momentum balance only. Since the flows under consideration are dissipative in nature, it is reasonable to consider the possible influence

-125of the energy effects on stability, The object is to establish, on the basis of energy considerations, the effect of the pseudo-viscosity vJ' (defined in Equation (84) ) on the stability of the interface in the system shown in Figure 27. The description of this system is given following Equation (93). The minimum unstable wavelength will be determined by computing the necessary relation between the time rate of change of the potential plus kinetic energy and the various rates of energy dissipation within and on the boundary of the fluid body. The core of the analysis is that the Darcy velocity gradients bui do, in fact, give rise to macroscopic shear stresses on the inside and on the surface of the volume element shown in Figure 28. The velocity field obtained for the potential flow of Section IIIA will be taken to provide a uniformly correct description of the flow for purposes of computing the energy terms. This amounts to assuming a viscous, irrotational flow on a macroscopic scale in the areas in front of and behind the interface. The velocity field is given rigorously for small perturbations by the potential solution, and the energy terms (for example, the dissipation function) follow from it according to known principles. Of course, the dissipation due to gross angular deformation can only be construed in a macroscopic sense, i.e., as limited to the Darcy approach to porous media. 1. Theoretical Development The generalized work-energy equations for fluid motion in porous media will now be derived, The general impulse-momentum equations for the unsteady case were derived in Section IIIA, and are (cf. Equation (89) )

-126- a Jj p u. ax + Xi i (111) at L J 6x.j p xi p6 K where T. is the deviatoric stress tensor, related to the velocity gradient 13 according to Equation (84). The velocity ui represents the total perturbed velocity. To transform the left side of Equation (111) into the appropriate energy term, note that, for example, Du au 1 a 2 22 au u aw = + - (uvV+w ) + v (- — ) + w ( Dt at 2 ax by ax az ax or, since the flow is considered irrotational Du aui 1 q2 1 1 l Dt at 2 axi Then, Equation (111) becomes 1 _+ I - - + X. + 1 - (112) at 2 xi p P ax PK J The total perturbed flow is the sum of the unperturbed (uj) flow and the disturbance flow u' j u = u. + u' (113) a J J where u' is the disturbance velocity, and u. is, in this case, simply J J the uniform velocity V. The total pressure is the sum of the unperturbed pressure p and the disturbance pressure p',

-127p p + p (1.14) Substitution of Equation (113) and Equation (114) into Equation (112) gives, 1 [ 6u: 1. q 6u 1. 6P 6u IF 1i + ~ \u _ _ + __p- -X +t 2 6x' J ax ax J p xi.xi (115) 4ui 1 YP- 4ui + X p K i P axi PK But, from Darcy's law for the undisturbed flow 1 aP.ui i P ax. PK 1 Also, 6uj 0 - since u is a constant (V) ax 1 Since the perturbation velocity is negligible compared with gross velocity, I aq' -2 a is deleted from the equation. 2 ax: The final equation for the disturbance only is9 _au: -+ u! ip.. 6uJ1 T + — I __^ + UJ __^ - - + - ~L- - "(.i.i6) at ax-. p axi P axi Pp Multiplying each of Equation (116) by its respective perturbation velocity u' gives

-1281 q i +. (uu) - + ui I 2 at j axi - ujuj xi (117) 2 6U.1 6u!!Iq 1 ra(uip') auil 1 a 1 u1 1iL P- axi ax p ax Ti ju Tij a pK bu: Using the incompressibility condition 1. = 0 to eliminate the axi third term on the left and the second term on the right, and integrating with respect to a moving body of fluid gives _1 1 a I + XU3 d2 1 a2i (u =-i -- - f q I dJI + - dai - (u p')dr 1 a 1 6u:! I 2 P+: P 71 ax: PK(1 \\^^^ ^^7~f~; p qxi] + Converting the appropriate volume integrals to surface integrals by means of Gauss' theorem, 1 2 1 1 1. P'-dS f qTdV + f niVuJuldS = - f niu!p'dS + - T ijunjdS Where n indicates the appropriate normal components of the quantities in question. The integrals in Equation (119) represent, in order of appearance

-129(1) Time rate of change of kinetic energy of the system. (2) Rate at which energy is being convected into the system, (Cross-linked energy effect). (3) Rate of working of external forces on system. (4) Rate of working of surface stresses on system. (5) Rate of energy dissipation in the volume of the system due to macroscopic angular deformation in the flow (dissipation function). (6) Rate of (Darcy) energy dissipation due to viscous drag of the fluid on the pore walls. The quantity uj is the uniform velocity V in the z direction. Converting the expressions in the surface integrals of Equation (120) to equivalent forms gives, P aLf q'dV + P f Vw1(utx+vay+wiaz )dS = - pT(utx+v+w z)dS 20 6z, 0, g 6n rn 6n g 6n dn 6n F 6x 6y T + f (UTxx xy xz +n + (U'T +V'T +W'T ) + (U'T +VT +WT S xx xy xz an yx yy yzn zx zy zz an - Ti - i df - f qI'dd (120) 1, 6x. PK It will now be shown that - and - are negligible compared to - under an an an the conditions of the present problem. The surface in space defined by Equation (94) is, nt k k z - Be cos - x cos - y =0 (121) \/2 N/

-130Ox The quantity - is the rate of change of x with respect to arc length an along the normal to the curve of the perturbation (cf. Kaplan (60), pp. 107-9). The unit normal to a curve is, in general n =cos i + cos i + os k where cos a, cos A, cos 7 = direction cosines - -~ -- = unit vectors in x, y and z directions. i, j, k, A vector which is normal to the surface defined by Equation (121) is, k nt k k - k nt k. k ^ = - - Be s in —cos —-y i - -Be cos -x sin-y J + k E2 2 2 2 2 E2 E2 ai + bj + k where F (x,y,z) denotes the function defining the surface in space. The direction cosines are then calculated as, 3x a cos a = / 2-~ an V a2+b +1 by b cos = = an A2+b2+ az 1 cos y =-= -"n - - a+b +1 The assumption involved in the first order stability theory is that the maximum height of the perturbation is not significant compared to the wavelength x, i.e.,

-13122 b max < < X max or 2T f2 Bent <<1 but, 2 k Therefore, the inequality is kBent < < 1 (123) It is easily seen that B2e2nsin' k: cos-k y + — Be eo.s k sink y < kB2e 2 4> 42-2 2 4 (124) From Equation (123), Equation (124) and the definitions of a and b in Equation (122), it can be concluded that, a2 + b2 < <1 and a << 1, b << 1. Therefore ax O; 0;aZ 1 (125) an an an The results of Equation (125) substituted in Equation (120)give the final form of the work-energy perturbation equation as,

-132a q2d + a f VwdS + J p'w'dS - f (u'T z+v'Tzy+W' Tz)dS 2 at zx zy zzS S + Tilj i d + -;P q'db = 0 (126) lv 6xi 8x PK V It is necessary to be quite clear that these energy terms apply to the perturbation only. Continuity of stress in the z direction at the interface between the fluids requires that, xz1 xzz zz1 xz2 yz2 zz2 1 2 z=b.Txz.i+ -, + T a = T- + T- + o C1+C2) | (127) But 6wi v zz1 -- P + 2 1 =- P T1 ZZ2 = - P2 + 2L2 aZ - P2 + TZZ2 using Equation (84) in Equation (127) gives i -,, 1i, I I Iu p + p ( + 2 + (. +.- ) (,28) I I Ix + (z 6 By o x 6y (128),, O2' u 2, / 2 2 =P2 + P2 + 22 + 2 + - + -- )+ (C 1+C2) 22 2u t2;y 62 6x 6Z Z-b To.compute p on z = b, note that,, ap b P=b z=o6 az + ---- (129) Pz= = Pz=O + z z=o

-133Also, 6p 6p? P ( + p, + () + b+ ---- (130) z=b z=o az az Z=o Combining Equation (129) and Equation (130) gives PT =p -p b + --- (131) z=b z=b z=o az z=o But -- Vz p = - _ - pgz and p is given by Equation (97) derived previously. Combining these results in Equation (131) yields P'=b - p + + + + higher order terms (132) z=b 6 t K z=b Solving Equation (128) successively for p!, and p2, and using Equation (132) in each case gives, P C1~ _ P2w 2V + P2 602 + 4202 Pl = (pl 2)gz -'2V + (^-^1 Vz2 6 K2 +a(2'w2 _i6w I (v2 u2 w2 u2 + 2( - + + + ++ ) 2 Pwl +2, aZ aZ ax ay 6x!z av alv aW a l - ^~J+..... + a + c +) (i33a) -x +y 6x 8z p - (~-l 2)Vz -(pl-P)gz -p + l + 4_'_! 2 K K2 K t 2 KI + 2(n, ~7+'2i- z) + 1(x -+ — + -~^+ — ) awlZ ax ay 6x z v1,'u2 zw2 1 u2 * - I.L2( + + —- + -) - (1+c2) (13lb) ox ay ax az

-134where D designates a perturbation velocity potential is the porosity. The velocity potentials and the perturbation to be used in this problem are given in Equation (95) and Equation (94), respectively. Substituting Equations (95a) and (95b) into the equation for work-energy (Equation (126) ), performing the integrations, and equating the form of Equation (126) for phase 1 to the equivalent form for phase 2 (the left sides of each form of Equation (126) being equal to zero) gives, n2(ol+P2) n(pP)- V(p2-Pl)nk - 2( f2-1)(+lt2)nk + n( -+ - \ )V + (P1-P2)g] k - ak3 = 0 (134) L ^2 i1 The surface integrals in Equation (126) are evaluated across 1 wavelength each in the x and y directions, the boundaries falling at the nodes of the perturbation. Volumetric integrals are evaluated from z = o to + o. It can be seen that the stability equation developed by the work-energy approach is identical with that developed by impulsemomentum (Equation (99) ) except for the term in (41 + p).

-135The expression for n is, [(12+ 1) _ (P2- 1)Vk _ 2(-2-1) (1 + 2)k2 K:2 1i 1 + J 2(P1 + P2) (135) I(ll2+1 (P2-P1)Vk 2222)(fL1+42+)k2]2+ (12)((2 + 2 + -~) +P (2 (+tZk K2 K1 ~1 1 2 KvK +(p2-P1)g]k c!~*~} 2(P1 + P2) In the case of viscous damping, the stability criterion is not changed. Examination of Equation (134) shows that Equations (101), (103) and (o10) hold whether the flow is damped in the manner described in this section or not. As previously discussed, only the value of n and the wavelength of maximum instability are changed. A series expansion n gives 2( ) - (P2 -P1)Vk - 2(2-1) )k 21 1 2 12 1 3 n ~a +CX - a+ -- 2(p1 + p2) L 2 -8 + - (136) where V 2 K1 1(prI, ) { L(- )V (p2-p1)g k ~k} a ~ " =:" ~

-136Examination of Equation (136) shows that the viscous effects contributed by (41 + 42) always damp the flow. 2. Discussion The boundary condition imposed on the solution at the interface is given in Equation (127) and states that the stresses in z direction must be equal. The perfectly rigorous boundary condition at the interface is that all the stress components be balanced, not just the normal component of the stress, i.e., that all of the following equations be satisfied. 1 a 2 +2 2 (1 7a) 1XXl an TyX1 an ZX1 an XX2 an YX2 an ZX2 an T a y T X+a -x+ T 1xy +yy an + Tzy an TXy2 E Y2 (137b) T x +y yz a _ x y zz (17c) 1 2 2 These equations represent the equality of the stress components on the fluid-fluid interface. From Equations (95) and (84), + av 2 nt kz k k Txy Tyx = (yu+ ) = -2 (Ak2ent e sinkx cos y) T T - c (1+ =) ^XH(Ak e e sin-x cos-y) ~xz =rzx = z Tx)'' XZ ZX aZ aX'[2 12 z =,(u+v w) = t,,(Ak ent e cos x siny) ax = - a -) + a -k n xz zz an2 Decomposing the normal stress into its vector components,

-137where an = total stress (normal plus deviatoric) in the i direction. It has been shown that ax y << a an an an for first order theory. Then ~ x n "zz The question arises as to whether the term a ora may, in some cases xan ysn be of the same order of magnitude as the deviatoric part of the stress component in the z direction, a. Since ax, say, is largely of order p~ (neglecting gravity and higher order terms), and the deviatoric component of r is of order i-zl, a quantitative estimate of their z relative magnitudes can be determined. From Darcy's law, p = - V (neglecting gravity) ~K.~ ~ ~x Then, from Equations(121), (122) ff., and (95) P- is of order on uV 2 2nt 2k 2k kBe cos -- x cos y k4 --— is of order 21 A k k 1 A kentek zcosyx cosy' y 1 "21"2 ~2

-138The ratio between azzand the deviatoric part of the normal stress T T ishenAk (T T ) is then, (using n = from Equation (95) ). n z B ax on kiVb =^ (138),,w',sn az Equation (138) shows, that regardless of the value of,V the ratio of K'~n ax to the deviatoric normal stress can be made as small as desired by making b (height of perturbation) sufficiently small (n is not dependent on the amplitude of b). Since the theory specifies that any instability however small will magnify under certain conditions, it is seen that the size of b presents no restrictions. Since the x and y components of stress can thus be made negligible without affecting the conditions of stability, it is considered sufficient to satisfy only the normal (z direction) stress continuity condition at the interface. Evaluation of Equation (126) shows that T.- d - (U'T + V'T + W'T ) dS (139) Tij 6Ox g zx zy zz!u! T. I is recognized as the macroscopic dissipation function, repre1~1ax senting the energy dissipated by macroscopic angular deformation. The surface integral represents the rate of working of the viscous stresses on the fluid boundaries. These quantities are equal and opposite for potential flow. Equation (126) then states that the change in kinetic energy content is equal to the kinetic energy flux plus the sum of the

^139work done on the surface minus the volum.etric Darcy dissipation. The result is quantitatively identical to the momentumrl approach except for the terms in' o The postulation of the stress deform.ation relation of Equation (84) is admittedly heuristic. In fact, the relation may ultimately turn out to be non-Newtonian, However, it is argued that the general structure of the solution should be given by Equation (134) It should be realized that the hypothesis contained in Equation (84) and (91) (namely, that the stress effects can be separated into the integrated mr!icroscopic (,>) and the gross macroscopic ( ui ) parts) are not meant to describe the exact details of the viscous damping so much as the net effects. gTaylor on page 590 of(O5)has given a succinct s atemernt, of this kind of macroscopic approach in speaking of his assumption that vorticity is transferred unchanged by turbulence in homogeneous flow; he said, "This is clearly a very drastic assumption, because it. is certainly untrue in details, though it may be true when considered in relation t, the effects produced on the mean motio"n," C Effect of Capillary Pressure Gradients on t:he Hydrodynaramic Stabil.ity of Quasi-Steady Flows The computation of the stability of quasi-steady flows, as presented in Section IIIA, did not take into account the possible effects of saturation gradients near the fluid-f.luid:interface The results obhtained are strictly correct only if no saturation gradent s are present, and the interface is sharpo Experimental verification of the simpler form of the stability equations has been obtained on Hele-Shaw cells.

-140Chuoke and co-workers (23) showed that the critical velocity (cf. Equation (104) ) was much higher than predicted from the sharp interface theory. That is, real systemsof porous media were much more stable than the simplified analysis predicted. In two phase flow in real porous media (as distinguished from the Hele-Shaw work), saturation gradients will exist near the transition region between the two phases. It is suggested that these saturation gradients give rise to capillary pressure effects in such a manner as to tend to cause stabilization of the flows, as observed in practice. To handle the problem analytically, it is convenient to assume that the saturations are invariant except in the immediate neighborhood of the interface. This saturation transition region has been treated experimentally and theoretically in Section IIB of this thesis. The object of the analysis is to incorporate the effect of these saturation gradients on hydrodynamic stability insofar as the saturation gradients give rise to (capillary) pressure gradients, which in turn will affect the pressure balance. 1. Theoretical Development The physical model assumed is exactly the same as that described in Section IIIA and illustrated in Figure 27. The saturation gradients behind the interface are the same as those illustrated schematically in Figure la. The fluid-fluid "interface" is to be taken arbitrarily near the point of farthest advance of the displacing phase, i.e., the point z in Figure 1 a. 0

-141The argument is as followso Capillary pressure is a known function of So But S is, from the work of Section TIB, a known function of Zo Therefore the capillary pressure gradient must be a known function of z according to the equation aPC dPc as ~^ = ~^ ~~ _ d~~~~~~~c (140) az dS az The strict definition of the mean interface will be that it is located at the point zo - 65 ioe., at a distance 6 behind the point of farthest advance of the displacing phase zo, Phase 1 will be defined as that which is located behind the interface z = b, cfo Figure 27; phase 2 being defined as that which is located ahead of z - bo The pressure in phases 1 and 2 will be changed due to the saturation gradients near the interface~ The effect of these capillary pressure gradients will now be examined in detail, The present analysis is not going to include the effects of (a) unsteady damping (b) inertial damping (c) viscous damping. This simplified approach is considered to be clearer and more illustrative of the phenomenon to be examined. Direct integration of Darcy's law gives, in the absence of capillary action K1 1 K1 P- - " - Plgz + K + pg (141a) P 12Vz g 202 + (4b = - ~ - p-gz +- + P (14lb) 2 K2 where po -= pressure on the interface in the absence of the perturbation.

-142Experimental measurements indicate that capillary effects will cause a reduction in the wetting phase pressure and an increase in the non-wetting phase pressure. These pressure changes occur over and above the hydrodynamic pressure gradients occurring in the non-capillary flow. Accordingly, the actual pressures acting in the non-wetting and wetting phases will be written as = dPc as lc- dc- z) (142a) Plc = P (a dS az z neglecting higher order terms dPc as P2 P= 2 + ( a + dS z) (42b) c z where subscript c indicates the total pressure acting under the influence of capillary pressure effects. a is a constant to be determined. The conditions which must be satisfied at the interface (z = z ), in the o absence of the disturbance are; Pl + P2 P 2 0 P2 - PI = Pc co co 0 where p.c = capillary pressure existing at the saturation occurring at the interface. p0 = mean pressure at the interface.

-143Solving these last two equations gives C ~ o c P = PO- 2 (143b) c Comparison of Equations (143) with Equations (142)(using Equations (141) on z = o and in the absence of the perturbation), shows that - PC a = - and that..the. resulting expressions for the phase pressures under the influence of capillarity are, from Equations (141) PCo dPc S lVz P-l1 - - z Plgz + (144a) c = o 2 dS az o 1 (144a),c0 dPc S __Vz Y2(b P -2 Po + 2 + z- o2gz + (144b) C K2 22 Then continuity at the interface requires that, P - P2 =- PC + ac (145) c c The displacing phase is here assumed to be the wetting phase. Substitution of Equations (144) into Equation (145) yields the value of n as,.(-)v ~1 dPc 6a 2l 41)v + (p,-p,)g - - -k0 (14) n( 2 + 2) - ) V + (PF -p)g dS a +azk I = (146) K2 K1 1 2 1 1 dS 6z z^i The stability criterion is _ v+ (p -p)g dPc- S 5 ( —1) V + (pp)g-~ - ok0 (147) "2 K1 dS ~z z ( "- o0

-144The critical velocity is given from,u2 il1 dPc as ( - )V +(P2-p1 dS a)gz o (148) K -2 c 2- 1 d zo Provided the perturbation contains wavelengths greater than, X 2 ~ ~ 2 (149) - 1 (V-V) K2 K ) The wavelength of maximum instability is $3 Xc from n = O - k 2. Discussion For purposes of making the previous calculations, it was assumed that the real, system (with saturation gradients) is replaced by the following model, The porous medium behind the interface was considered to be filled with a homogeneous fluid of viscosity u. The pressure in this phase is, however, assumed to be reduced by the amount indicated by the capillary pressure gradients, It is also assumed that the saturation distribution underlying the process near the interface is changed negligibly from the steady state by the disturbance, The position of the interface z = b can be made arbitrarily close to zo (cf. Figure la)by making the quantity 5 arbitratily small (the mean interface is at zo - 5). It is then possible to compute the dPc (6S quantity dS (;) of Equation (147) for the purpose of quantitatively determining the effect of capillary pressure gradients on flow stability. Within the limits of the stated assumptions, this quantity can be calculated approximately as, cf. Equation (29)

-145dP C aZs V a + ]PK (150) dS az z l D1 VO It will be noted that Equation (146) predicts that capillary pressure gradients will tend to stabilize the flow in the case where the disdP, aS placing phase is wetting, since in this case, both d c and - are negative. In this case where the non-wetting phase is being injected, the definition of capillary pressure becomes P -P =P lc 2c c and Equation (146) becomes P-2 1'l (2 l dPcS 1 2 +K - 2i ( 2-1 +z k -P&-k 0 (151) where S is, as before, the wetting phase saturation. Equation (151) shows that capillary pressure gradients will also tend to stabilize flows as in which the non-wetting phase is being injected, since in this case - az dPc is positive and d is negative. dS Experimental data on systems where capillary damping is presumably taking place is almost non-existent, except for some given in a paper by Chuoke, et.al. (23). An example calculation will now be made on one of these experiments. Water (u = 1 cp) is displacing oil (u = 200 cp) horizontally from a sand pack of permeability 200 Darcies at a volumetric velocity of.115 cm/sec. The interfacial tension is 25 dynes/cm. Since the pack is unconsolidated sand, the relative permeability curve of

-146Figure 3 will be used. The Leverett J factor curve of Figure 4 is used for the capillary pressure relation. Figure 5 of Chuoke, et.al. (23) shows that the experimentally observed finger spacing was 3 to 4 cm. The stability condition can be written, from Equation (146) and Equation (150), as (gravity neglected because of horizontal flow) [ 42 1 a2 1 ] Vc + (p p,)g = 0 (152) The term D21 represents a pseudo inverse mobility. a2 is determined as KDi 2 4.0 by the stabilized zone technique (57) and D1 is.025 from Figure 3. In the absence of gravitational effects (horizontal flow in a very thin flat plate), Vc = 0 and the wavelength of maximum instability is = ___2 _1 _ ) 2 (153) K2 K1 KD1 The naturally occurring wavelength (xm) is computed from Equation (153) as.364 cm. X in the non-capillary case (a2 = 0) is computed to be m icdt.10 m Lager wavelengths are associated with inreants tend to stability It the flow appreciably. However, the degree of stabilization is not enough to account for the 3 to 4 cm wavelength actually observed in the experiment, although the correction is appreciable and in the right direction. The stabilizing effect of capillary pressure gradients can be illustrated more vividly by assuming that the above experiment is conducted.with the water displacing the oil upward at a density difference

-147lb of 15 -- The non-capillary critical velocity is then computed from ftft Equation (106) to be o658 f and the capillarity stabilized critical ft ~ day velocity is 3.36 from Equation 152, The critical velocity determiday nation is thus expected to yield a more critical test of the present theory than could be obtained from. measurements of naturally occurring wavelengths. Unfortunately, no experimental data on critical velocity are currently available, For neutrally wet materials, capillary forces are presumably not supposed to exist. Equation (146) shows that, for given fluid-rock properties and displacement rate, the neutrally wet materials should exhibit the highest degree of instability because of the absence of dPc as dS (a~z) However, for porous media classified by them as neutrally wet, Chuoke, etoalo (23) observe that the factor 2g f3 must be replaced by 30 in order to achieve an accurate correlation, This means that the neutrally wet materials do, in fact, show a greatly increased stability over the predictions of non-capillary theoryo The cause of this increased stability is unknown at present, The increases in the naturally occurring wavelengths up to factors of 20 (compare 3-4 cm from Figure 5 of Chuoke, etoalo (23) with the computed non-capillary X of.160 cm.) suggest that, although the present analysis gives corrections in the right direction, the magnitude of the stabilization is not accounted for, Also, severe difficulties are encountered in obtaining accurate enough experimental data to yield quantitative predictions, This arises from the fact that the denominator of the quantity in brackets in Equation (155) must be the difference of two numbers whose magnitudes are very nearly equal (I2 1 and a2nd ) K2 1 KD1

-148if any appreciable correction in Xm is to occur. Since the accuracy of these quantities is, at best of the order of 4% a reliably calculated increase in Xm of a factor of more than about 5 cannot be expected (cfo Equation (1553) ). However, since the volumetric velocity V of critical stability is normally of greater practical significance, the conclusion can be drawn that the present analysis will yield useful results. For example, the calculation on page 146 shows that for water displacing oil, capillary ft pressure gradients can increase the critical velocity from.658 da to day ft 3036 day @ This factor of 5 increase can be accurately computed from other experimental measurements on the porous medium (i.e., K, K~ K2' ) etc.) It can be concluded from this work that Equation (146) could represent the minimum, degree of stabilization possible in real porous media. Since the data of Chuoke, et.al, (23) indicate that the degree of stabilization is even greater than this, we must turn to some other stabilization mechanism to account for the results, This mechanism is thought to be related to annihilation of the fingers by capillary crossflow, and will be discussed later in Section IIIG. Do Stability of Plane Interfaces in Accelerated Flow Previous work in this thesis has dealt with the problem of stability of plane interfaces in quasi-steady flow (Section IIIA). The basis for those calculations was that the fluids were initially found flowing at a constant velocity V and separated by a plane interface, It is now proposed to examine the case in which the fluids are accelerated

-149continuously from a state of rest. Taylor (120) has examined this problem for homogeneous fluids. 1. Theoretical Development The model to be assumed is that shown in Figure 27 except that phase 1 is displacing phase 2 at a velocity equal to alt, where al is the acceleration. Going through exactly the same development as in Section IIIA yields as the counterparts of Equation (97) and Equation (100), respectively. p - pV t P atz- - p(g + a) z + p (154) ^ ^ at PK K where the kinetic energy term has been neglected and (155) n (P!+P2) (P -P )naltk 1.2 L1 L n(p1 ( -p a + n (!-+ 1) (2(- [ l )alt + (P2-Pl)(g+al)]kj k == O K2 Kl LK2 K1 where a is the upward acceleration. The axes are again fixed in the plane of the unperturbed interface. The inertial body force term is introduced by subtracting palz from the pressure in each phase. This downward pull is the force experienced by an object which is being accelerated upward. 2. Discussion Examination of Equation (155) shows that the acceleration tends to stabilize the flow when the acceleration is upward and to make it less stable when downward (assuming lower fluid is always the more dense under the influece of gravity). Equation (155) also shows that, unless al is significant compared to g, then the effect of acceleration is negligible

-150as far as the stability is concerned. It can thus be concluded that for vertical flow in real porous media, the effect of any usual acceleration - 2 on the stability is negligible. (Compare, say, 10 3 ft/sec with 32.2 2 ft/sec ). The stability is thus a function only of the instantaneous velocity and not of the instantaneous rate of acceleration. For flows in which gravity can be considered negligible, the rate of acceleration could be an important factor, however. Neglecting the first two terms and g in Equation (155) (unsteady and inertial damping effects) yields, 2n 1 l)] [(L2 l )at + (P2-1)a k - 0, (156).K2 Ki J KL ^2 Ki 1 P These conditions will normally hold in horizontal flows through porous media whose vertical thickness is small compared to the lateral dimensions. The consequences of Equation (156) with respect to critical velocity, critical wavelength, and naturally occurring wavelength are the same as those derived by Chuoke et.al. (23) except that a is replaced by g. The significance of gross fluid acceleration on oil-water movement will now be illustrated by an example calculation. Water of A = 1 cp is being injected into a thin oil bearing (p = 2 cp), horizontal, linear porous formation (K = 1D), at such a rate that the volumetric velocity is given by V = alt with a1 = 10-6 ft/sec2. Compute the distance which the interface has moved at the instant that the flow becomes unstable if the acceleration is assumed constant. The time at which instability wilL:.... occur is given by t = -P1 P2 = 9.83 x 10-7 sec. The distance moved in (2 ~1

-151this time is Distance = alt2 = o968 x 1-18ft. The velocity attained is V = a t = ~983 x 1012 ft/sec = 3.1 x 10-5 ft/year. It can be concluded that, in the horizontal systems described above, instability is achieved nearly instantaneously even under the influence of the largest accelerations and the largest mobility commonly found in reservoir problems. As in the (non-gravitational) case, the critical velocity is still virtually zero. Gross fluid accelerations thus have no significance on the stability of displacements occurring in underground reservoirs. E. Stability of Stratified Flows in Porous Media The previous work has dealt with the stability of flows in porous media when one fluid is displacing another frontally. Many problems in reservoir mechanics involve the parallel, or stratified flow of fluids past one another when separated by a more or less well defined interfaceo Gravitational drainage is usually the cause of this fluid stratification, as when oil is found above the water in a producing zone. Stratification may also occur when the same fluid is flowing through a multilayered porous formation which has zones of varying permeability. In either case, the difference in fluid mobility in various parts of the porous medium can generate effects which lead to unstable behavior.

-152-,1 Theoretical Development It is assumed that two fluids are flowing in the direction through a porous medium of permeability K separated by an interface at z = 0 (cfo Figure 30). Both fluid zones are assumed to extend semiinfinitely in the plus and minus z directions and infinitely in the x and y directions. The z axis is directed positively upward. The pressure drop across any section perpendicular to the zy plane is constant in each phase. The perturbation is given by i k bBenmt f2 (x+y) b = Be e (157)* Note that in this case, the coordinate system is not moving with any unperturbed interface. It is fixed at an origin which is arbitrarily assumed to be located in the plane z =0 The total pressure drop between any planes A and B (cf. Figure 30) of the two-layer system is assumed to be known and the system is at steady state. Under these conditions there will be no net cross-flow across the interface, Darcy's law states that, 1 6Pl K2 6P2 V = - l; v - -- (158) But, since the pressure drop across both layers is constant p~ -p2 ~then ~V1 2 (159) axI x -I x Note: The wave numbers of x and y could have been made non-equal, but this does not contribute to the generality of the solution.

-153A E +0 Kv —-*~VI X2 ~aV2 VII - I y FIgure30. Sriie FlwoFlisnPrusMdaI. I V, - I Y -00 Figure 30. Stratified Flow of Fluids in Porous Media.

-154If G(x, y, z, t) is the equation of the interface, then DG aG dG 6G dG - = + u- + v- + w- = 0 (160) Dt at ax ay az But G = z - b Velocity potentials satisfying the steady flow, the Laplace equations and the vanishing condition at + o are, i k (+y) = Aleintekze - (x+y= l Vlx (161a) i k et - ve e Vx = - V x (161b) ~2 2 2 2 (l6lb) Neglect of second order terms in Equation (160) leads to 6b db w = + V- (162) at the interface for each phase. Use of Equation (157) and Equation (161) in Equation (162) leads to A1 - - iB(n+-7); A2 = iB(n + Y) 12 The Bernoulli-Darcy Equations are, ignoring unsteady effects r2 1plql' j1t1 = 1 2 p + l plgz + C (163a) P2 = 1q + 22- p gz + C2 (16 2 1 K2

-155 but, 2 c^~D 2 2 Oi 2 oi v~ qbut (q-=)2 + +(-7) + (-F)2 - 2V.i + + i 6x 6y oz i In the absence of the perturbation, the pressure on the surface z 0 xsCiVix is p ='i., where o is the ambient pressure at the arbitrarily chosen origin. The constants in Equations (163) are then evaluated as, 1 2 -i Po + PiVi also a('.i) (+,,) + (ai) ax 6y az Note that C indicates the potential of the disturbance only and that 0 is the potential of the total flow. Then PlVl1 a 1~11 i x1 P P O + X px+ - Plgz (164a) 1 ^ ^ 6x K^ K pl 0 + p a 2 + 2-2 P2gz (.64b) The condition at the interface is P - P2 = (c + c) (165) 1 2 Substitution of Equations (164) into Equation (165), use of Equations (161). and computation of the curvatures in the usual manner gives,

-1562 2' \ kV1 1 2 kV21 - +i \_- - + L pl^J + - (P 2-Pl)g k +i n = (166) V PlV1 P2V2 ] i[11 2 1 L X + 2 k-l kK n is thus a complex number which is the ratio of two complex numbers. Let n be represented as a + i b n = _. - (167) c +i d The division indicated in Equation (167) yields the following complex number for n 2 _2 n =, - l cos (e - Q ) + i sin (e - Be ) (168) c +d The quantity.n is represented in the form n R + inl (169) where the real part nR represents the celerity or speed of movement of the wavy disturbance on the interface and n4 represents the amplification or damping factor. The stability of the flow depends on whether no is positive or negative. The Argand diagram representing the components of n is given in Figure 31(a). Since b is always positive, and a can be positive or negative, a + i b must always lie in quadrants I or II. Since c is always positive and d is always negative, c + i d must always lie in quadrant IV.

V II I'0 6 STABL ZONE\ -71~ ~ 1 O X ~ V~\~C~ ^ -7T^UN~~~~~~~~~~~~~S'T^iFDRBID ~' ABEDEN &Bs &~d Figure 31a. Argand Diagram Showing Figure 31b. Zones of Stability for Stratified Components of Instability Flows in Porous Media. Index n.

-158Therefore, - 2 < (Oa_ - 0_-) < - - (170) ab cd 2 where the angles e-b and 0 — are shown in Figure 31(a) and are defined by, a c Cos a COS= 2 2 — ab cd Reference to Equation (157) and Equation (168) shows that sin(0b - 0-d) must be positive for the system to be stable. With this information and Equation (170), the stability sectors for the flow under consideration can be represented graphically. This is shown in Figure 31(b). The argument leading to Equation (168) serves to define - < Oab- E <0 as a forbidden zone, i.e., a zone in which the angle -- - e-d cannot occur under the conditions of the problem. The section -2t < 60 - E < - E gives positive ni, so that the sector of instability is limited to the range -A < 9a - 9d < - A which is shaded in Figure 31(b). The existence of the forbidden zone implies that the amplification factor is unique when the system is unstable (there is only one value of (a6- - -a) which yields any given sin (a- - 0-) ). When the flow is stable, Figure 31(b) shows that the damping factor is not uniquely a function of ga - e. a~ - e1" equal to -_ is the point of critical stability. Let LB = 9 - 9O, then a c b d cos A0 - + (171a) a2+ a2 c2 + dd2 2 + d2 b c ad sin A0 = + (171b) a2 f y2 T2 2 72 2 2+ d

-159Using Equation (171) in Equation (168) gives n a c (bc - a- + (d) (= +d) ()2 + 2) (172) The instability criterion is that be - ad be negative, i.e., that 2 112 41 V k ~ LJK2) + (1~1 K12 (7 (K2 ) 2k (P1 + P 2 1 ( P ) [(2 P) >0 (173) where bc - ad has been computed from Equation (166) and Equation (167) and Equation (159) has been used to eliminate V2o The wave number of critical instability is determined by solving Equation (173) for k (with the left side equal to zero). v2 4~l2) 2 2(4 K2- 2,) (p1 -- 2 2 - 1 P2 i ) + 4p2'1 ( 2 2l U ru1o M.-1"i^ + 4 a 1-2+ -j ( P2 - P1 )g k = 2 cr 2 1 (174) The wave number of maximum instability is determined by differentiating n with respect to k. The result is, ~2 4 ~.L1K2~ + (~ - ~i)2 (~+'~ ~l~2)2 V1 (12 _. )(p+p:2) \ (2. ) (p+p 12) 20 K2 K1 21 2 2K 402 2 K1 21 + 12 j -$. + j (p2-p1)g m - 6 i(~*[12 + i ~] 1 ^rJ2C2 8J C(175)

-16oEquation (174) was obtained by means of the following simplification. For porous media c is of the order 102 or less and d is of the order 106 or more. Then _2 -2 -2 c +d d The velocity df critical stability is determined by eliminating a between Equation (174) and Equation (175). This gives stability condition. V2 < 4 1 + -2 [(Pl - P2)g 1 (176) i2 c1 1 E 1 Since k is related to V1, it is necessary to use Equation (175) and Equation (176) in conjunction to determine whether or not a given flow will be stable or not. k sis computed from the fluid and rock m properties and the velocity V which is linked to the stated pressure 1 drop through Equation (158). This value of k is then used in Equation (176) to determine whether the inequality is actually fulfilled. There are thus only certain combinations of V1 and k which will actually indicate instability according to Equation (176). From Equations(172), (166) and (167), the real part of n is kV2 rkVV __ ( —2 + 2 0 } [ 2JJ L~1{2 +2T2Jl k~1 kK2'R = - F~ 1 +

-161From which an order of magnitude analysis yields $2 AlVlk n2 - iV1 (177) From Equation (157) and Equation (168), the perturbation is given by, k k b: e'n l tei(nR + J2 x + Y) of which the real part is, -in t k 2Vlt b =e n cos [ 2V x - y (178) (1 + 12 ] 2V1 The waves then move along the interface with a velocity of 2V1 at the same time as they arebeing damped or amplified. (1 + 2) The stability criterion can be expressed in more compact form by means of the following definitions V2 - = Fr = Froude number (178a) gkK +2 [P1 P2 Kl 121L (178b) kl*2 - 2 [P~p~lK~j[~L] (178b) L12 ~]J 2l 40 [e [p P1-P2 L"2 KiJ Letting ki = k2 = k in Equation (176) and rearranging gives the stability condition as, Fr < - k1_A

-162where k* is thought of as the reduced wave number, and the Froude number is especially defined for the problem. Use of the method will now be illustrated by an example calculation. A perturbation of wave number 102ft1 is imposed as a porous formation which has permeability 10 Darcies. The formation contains oil (i = 10 cp) overlying water (I = 1 cp) in a horizontally stratified arrangement (cfo Figure30). The density difference between the oil and water is 5 lb m/ft3 and the interfacial tension is 5 dynes/cm. Calculate the pressure drop at which instability occurs when the fluids are made to flow in a direction parallel to the interface. Also calculate the celerity of these surface waves. Assume 1 = = = 10 Darcies. The velocity of critical stability is, from Equation (176), V =.10 ft/sec. The pressure drop corresponding to this critical velocity is 133 psi per foot (from Darcy's law). The celerity of the surface waves is, from Equation (178) ff. 2V1 2-1 = 0.182 ft/sec. 1 + — 2 K 1 K-2 Fo Stability of Radial, Quasi-Steady Flows, with Application to Thick and Thin Gas Storage Reservoirs Previous sections of this thesis have dealt with the stability of fluid-fluid displacements in porous media in cases where linear geometry prevailed. In practice, however, gas storage reservoirs and

-163oil fields usually are found in relatively flat, radially symmetric shapes. The question arises as to the conditions of stability which will prevail when a less viscous fluid is displacing a more viscous fluid in a radial manner. Specifically, it is desired to determine the rate at which gas can be injected into a gas storage field before the gas no longer expands in a radial disc, but begins to protrude into the confining aquifer in the form of long fingers. The solutions presented here also have application to the recovery of oil in which water is injected at the perimeter of an oil zone and serves to drive the oil toward a producing well. 1. Theoretical Development A flat, disc-shaped gas storage field composed of phase 1 within a radius r = a has a line source of phase 1, of strength a, at r = a. The surrounding aquifer, composed of phase 2 (immiscible with phase 1), is infinite. The reservoir has a finite height. The physical model is illustrated in Figure 32. There is no gross variation in radius between different planes perpendicular to the z axis. The positive z axis is in the upward direction. Both fluids are assumed to be subjected to sufficiently small pressure drops that their densities do not vary significantly. The governing equation is _2 2 2 be 1 b2~ 1 be - + -+- -1 +- =0 (179) ar2 r2 2 r Or az2 Where r, z and 0 are the radial, longitudinal and angular directions, respectively, and ~ is the velocity potential of the total disturbed flow.

00 PHASE 2 PHASE L. U1 I. — a b Figure 32. Top and Side View of Disc Shaped Reservoir of Finite Height. Radial and Longitudinal Perturbations are Illustrated.

-165The undisturbed source flow is = - aln r (180) u where X = line source of strength in ft3/sec.fto The form of the perturbation is b = Be cos m O cos kz (181) where k = perturbation wavelength in z direction, i.e., 21r kz k and m represents the wave frequency 2it e where 0 is the wave angle, i.e., the angle subtended at the origin by one angular perturbation. The angular wavelength is given by 2 - A - - = Oa k where k = angular wave number. The velocity potential of the perturbation only is, D = (182) Substituting Equation (182) into Equation (179) gives the equation governing the disturbance only,

-1662 2 2 + - - + - + 0 (183) ar2 r ae r ar az The solution to Equation (183) can be obtained by separation of variables or direct construction from the equation of the perturbation as, ( = ICl Im(kr) + C2 Km (kr) e ncos m 6 cos kz Using the condition that 1 must be finite at r = o and that lim r 2= O, we obtain, r — oo 1 = Cl I (kr) e cos m Q cos kz (184a) ~2 = C K (kr) entcos m 6 cos kz (184b) where m must be integral. Im is the modified Bessel function of the first kind and of order m. K is the modified Bessel function of the -m second kind of order m. C and C are arbitrary constants. -1 -2 The substantial derivative of the separating interface is, DG uG + u - + + w - 0 (160a) Dt at ar r aO az where G =r - a- b u, v, w = velocities in r, 0 and z directions, respectively. Neglect of second order terms yields b= u at the interface, (185) at

-167Continuity of velocity at the interface gives ul = u2 at r = a + b (186) Use of Equations(185) and (186) allows calculation of the constants in Equations (184) as -nB 1k IM+ (ka) +k Im (ka)] nB c (188) k K _m(ka) + k K (ka) The pressure balance at the interface gives p1 -p2 = Pc + C (1 +c ) on r = a+b (189) where phase 2 is wetting. The curvatures are calculated as a2b c1 = - a2 (190a) r2 + 2(2r2 (2r -) 2 r 1+ 1 62r.(a 2 ___ - r (190b) r2 = r 2 3/2 r r2 6&2 r +( )] From, cu k ap r r

-168the undisturbed pressure distribution is calculated as r p = po - in (a) Inclusion of the capillary pressure gradients and the pressure discontinuity due to surface tnesion in the manner illustrated in Section IIIF, Equation (141) through Equation (144) gives the Bernoulli-Darcy equations for radial flow as, P Pg~ Pco dPc LS r i1 a o 2 PO + 2 +(r-a) i- l ) + 1 + a l P2a (191a) P2 = P +.... (ra) n (r) + _n K2 2 2a P dP S Lr 1 a* co ~c, ^ ^2^ ^i^ p.^ cr (19lb) recognizing that, In r = in b) b a a a b for small a.

-169The value of n is obtained by substituting Equations (184) into Equations (191), using Equations (190) to compute the curvatures. The result is, from Equation (189), n[ 2n((ka) I Im(ka) LK2[k l(ka) + kK(ka)] Kl[klm+(ka) + IKm(ka)] (192).e Ira W ~1 1 dPc 3S ( I2 ki 2 2 1 c a K2 K + ) + dS dr=a In the preceding analysis, inertial damping, unsteady damping, viscous damping, and gravitational effects have all been neglected. The quantity ka will normally be of the order of 104 or higher. For these large arguments, ka -ka e K (ka) 2 e- I (ka) e -m2ka m 2 The use of these simplifications in Equation (192) gives - 1 c as o n + kl2 l) - ( _2_ 1) + a(k + k2 1 dc as(193) cr (k + a-+ - a) + --- k ~ a (193) The stability criterion is D (2h) k2 + k + 2 1) + cs I (194) al a K r=a or, alternately, the system is unstable if (-2 ( Il) + dPc as a K2 K1 dS r r=a

-170and the disturbance contains wavelengths greater than r ] 1/2 X V112 1 CD - ~b)^ ) ~K2 K1 a The wavelength of maximum instability (naturally occurring wavelength) is determined from an/at = 0 ~2 2 52 k 2 where 72 = k2 + (1 +5)k where 5 is a proportionality factor, The result is 1/2 x~m = [- (195) K2 KI Since the quantity ~-..- will rarely be more than about 10"-, it can be concluded that, n =-3 Xc Equation (195) shows that capillary pressure stabilizes the flow, since dPc as _c_ is negative for phase 1 wetting. dS 6r The preceding development treated the case where the discshaped reservoir had a finite dimension in the z direction. Since many gas storage reservoirs are very thin relative to their lateral dimensions, it is natural to investigate the stability of the system which has zero thickness0

-171The schematic diagram illustrating the physical model is the same as that shown in Figure 32, except that the side view is reduced to a line. 2T = 2 T o (196) which becomes, a2 +L ai + La = 0 (197) ar2 r2 ar r ar under the transformation u Using the same nomenclature and boundary conditions as in the solution for the thick reservoir gives, nB 1-m m nt.1 =-m a r e cos mg (198a) m 3= al+mr'ment cos m (198b) 2 m where the original perturbation is nt b = B cos mO (199) Use of Equation (189), Equation (190) and Equation (196) through Equations (199) for proper formulation of the pressure balance gives n (- + 1l) - - (. _ ) _ 2 ) dPc S. (200) k 2 K a = 0T a=dS (200) k K2 K1 a -2 K1 a dS -r dS ~rr=a

-172The system is unstable if the source strength is larger than a critical source strength defined by Mc,, 2 dPc 6S c (2 _ 1)+ (_2 =0 (201) a K K dS r r r=a and the disturbance contains wavelengths greater than.1/2 X = 27 [ ~1/,c -L(2 _- i)(-' c)J K2 K1 a 6n r The root of -k = 0 gives, = \3 Xc under the simplification previously noted.. Discussion The results obtained here indicate that gas can be injected unstably into a storage reservoir under certain conditions. Thus, calculations on underground water movement based on the existence of stable displacement of water by a gas zone having a discrete, definable radius may not be based on the correct physical model. Certainly the pressure behavior will be expected to be different in the case where the gas enters the aquifer in long fingers than in the case where frontal displacement occurs. The solution shown in Equation (184) are based on the assumption that m is integral. This is derived from the expectation that the perimeter of the storage field will contain only a whole number of perturbations. Non-integral m yields the solution, 0 = C11 (kr) e cos mg cos kz + C2I (kr) cos mg cos kz

-173and the detailed working out of the problem is the same as previously shown, The use of the results developed here will be illustrated by means of an example calculation, A gas storage field of permeability 100 md and porosity of 25% has a radius of 1000 feet and a well bottom pressure of 500 psi. The thickness is 10 feet~ The aquifer is assumed to be infinite. The gas and aquifer brine viscosities are o011 cp and 1.2 cp, respectively, The equilibrium capillary pressure gradient at low water saturations has been determined experimentally on core samples to be.2 psi/ft. Assuming that the field can be considered to be an infinitely thin disc and that the saturation gradients at the field boundaries correspond essentially to capillary equilibrium, calculate the maximum rate of gas injection which is compatible with stable reservoir development. From Equation (201), the critical source strength is computed as nc =.0122 ft3/ft.sec (actual ft3 at 500 psi). This corresponds to a rate of gas injection of about 5,300,000 standard cubic feet of gas per day into a well with 10 feet of sand exposed. Rates of injection higher than this will be expected to cause the gas to enter the aquifer in form of long fingers, Go Effect of Capillary Imbibition Between Fingers on the Stability of Quasi-Steady Flows All the work heretofore done describing the stability of fluid-fluid displacements in porous media has been confined to the purely hydrodynamic influences, These include the effects of frictional pressure drops, inertial, unsteady, and viscous damping, effect

-174of capillary pressure gradients, and the effect of acceleration. It was shown in Section IIIC that the inclusion of the capillary pressure gradient in the Bernoulli-Darcy equations (cf. Equations (144)) predicted an increased stability for displacements in porous media. However, the data taken by Chuoke, et al. (23) showed that the water wet system was even more stable than predicted by this correction. It is postulated here that the very large stabilizing effect contributed by the presence of the porous media (as distinguished from the parallel plate models) is due to the sidewise imbibition of the wetting phase fluid into the fingers, Physically, the effect can be described as follows. The discrete gas finger is presumed to start forming under the conditions outlined in Section IIIAo The sharp saturation gradients existing at the edge of these gas fingers will naturally cause high capillary pressure gradients in a direction perpendicular to the gross flow, The water will then flow laterally, filling in the fingers and thus annihilating them in the process. The net effect of this flow will be to keep the flow stable until much higher velocities than predicted by the parallel plate model are achieved. However, velocities will eventually be reached at which even the relatively high sidewise capillary flow potentials can no longer counterbalance the unstable tendencies inherent in the hydrodynamics. Then gas fingers will actually form, and progress forward faster than the laterally oriented capillary forces can fill them in, This gradual thickening of unstable fingers is illustrated during the course of experiments in Figure 2, It is proposed, in the following

-175y= l AS/y=o -W III GAS WATER Figure 33. Single Gas Finger Penetrating Water-Saturated Porous Formation.

-176development, to obtain a quantitative relationship which will predict the displacement velocities which actually must be obtained in order to overcome the tendency of the capillary forces to annihilate the fingers. 1. Theoretical Development Gas is assumed to be displacing water from a porous medium of known properties. Water is considered to be the wetting phase. The finger under consideration is assumed to be fully developed into the shape shown in Figure 33. The flow is two-dimensional with gravitational effects negligible. The gas finger is propagating in the z direction. A single unit of one wavelength X in breadth is shown in Figure 33. Not far behind the nose of the finger, the finger boundaries become essentially parallel to the z axis (cf. Figure 9 of Saffman and Taylor (105)). Water imbibes from the edge of the finger at y = Q to the center of symmetry at y = 0. y = 0 is considered as an insulated boundary since no water imbibes across it because of the vanishing of the saturation gradient there. Water imbibes into the space 0 < y < i as gas is expelled countercurrently and carried away in the main stream. In making the same approximations as those leading to Equation (77) and using Equation (4) to compute dPc/dS, the governing equation becomes _ -.0f(g)ey(o 2)1- 2 (202) as = Ky PC aS =.060 f(9), K)" S a s DC S (202) dt O-r1w dS z 0 AS' Iw for unconsolidated media. The constant.060 was obtained from a best linear fit of the KW and J(S) data from Figures 3 and 4.

-177f(9) is some function of contact angle (cfo Equation (4) ff, ) It will be assumed that f(G) must be obtained experimentally in the laboratoryo The imbibition is assumed to be linear and countercurrent in the y direction. The boundary conditions are lo S (y,0) = 0 20 1yy (Ot) = 0 30 s (et) = lo00 The solution to this problem is given in (22), page 137, and is, for the center of symmetryo c0 S(o,t) (-1)n 2 erfc (2n +l } n=0 2 Dct S is the normalized saturation. The time which is taken for the saturation at y = 0 to rise to.99 is computed as followso S =.99 when. 354, From 2:b t c this 212 t = (203) 099 - c It can reasonably be assumed that the bulk of the fingers will be filled when the saturation at the center of symmetry reaches.99, If the distance of forward travel of the nose of the finger is less than i in the time taken for the finger to fill in (t 99) then the flow must be considered stable, ie.,,

-178" Under these conditions, the finger is filling in sidewise faster than the nose of the finger can advance, The fraction of the two-dimensional system which is flooded out by gas is, 2a oa = x If the total volumetric velocity (volumetric velocity far in front of the finger) is V, then the finger velocity U, must be (cf. Saffman and Taylor (105)), U = (205) 1 -a The time taken for the finger to travel a distance & is thus, t.= (1 - (206) = V Substituting from Equation (203) and Equation (206), Equation (204) gives VaX D > (I -) (2o7) The naturally occurring wavelength X is given from previous theory in Equation (103)o (The displacing phase is designated as phase 1 (gas in this case) while the displaced phase is designated as phase 2 (water in this case))o Substituting Equation (103) and Dc (cf. Equation (202)) into Equation (207) and assuming that K1 = K2 = K gives, f(g)'2J \7 aV \*6, ^ a)(/ ( > ~ 0V~ (/,2 (208) ^2^~ J AS'2 o 06(1 - a)(v ( V )2

-179as the condition under under which the flow is stable, Saffman and Taylor (105) have shown that a = 1/2 over a wide range of conditions in parallel plate models, Inspection of Figures 4 and 5 of the paper by Chuoke, et al, (23) shows that a = 1/2 is a reasonably accurate number for real porous media~ Van Meurs (125) reports recoveries of 47 to 53 percent in the displacement of oil by water from a sand pack of dimensions 18 x 9 x 2 cmo For horizontal flow Vc O= The final equation defining the velocity of critical stability is, from Equation (208) _ /2 = f() )a(2 1)/( c 33a4ASsTr 3 04 Note = V indicates the critical velocity computed as in Section IIIA without the inclusion of the effects of capillary annihilationo Vc includes capillary annihilation, Velocities greater than V of Equation (209) indicates that the stability of fluid flows in which capillary annihilation is occurring is not dependent on the permeability, The forces governing capillary crossflow, flow under a pressure gradient, and finger spacing are apparently balanced in just such a way that permeability cancels out of the equation, 2,Discussion The stability criterion embodied in Equation (209) provides a method for predicting the stability of real porous media under conditions where the displacement takes place at an adverse viscosity ratio, It will be expected to apply in a very practical sense, to

-180systems in which viscous oil is being displaced by water. The parallel plate models, since they do not take into account the possibility of capillary annihilation, will tend to predict smaller critical displacement velocities than are actually observed. Chuoke, et al. (23) obtained agreement between theory and experiment in Hele-Shaw cells, but actually observed experimentally that real porous media showed increased stability. In the work of Chuoke, et al. displacement velocities had to be raised much higher than predicted by the parallel plate theory before instability actually occurred, In order to test the ability of Equation (209) to correctly predict the stability of displacements in real porous media, experimental data from the literature were compared with Vc from Equation (209). Before turning to the actual data, the general character of displacements at adverse viscosity ratios will be discussed. The set of experimental data which best illustrate the effect we are looking for are given in Figure 8 of the paper by De Haan (29). (Note that the lower graphs of Figures 8 and 9 of (29) should be reversed.) At very low displacement rates, Figure 8 of (29) shows that the displacement is essentially frontal and stable up to a velocity of about 10-3 cm/sec for the displacement at ~f = 25. However, between velocities of 10-3 cm/sec and about 2 x 10-2 cm/sec the breakthrough recovery drops to a low value and then levels off as velocities are increased further. This behavior is interpreted to mean that the fingers are being totally annihilated at velocities below 10-3cm/sec and are fully developed at velocities above about 2 x 10-2 cm/sec.

-181" In the intermediate velocity range, it is considered that the fingers are being only partially annihilated, Of course it is the percentage recovery at breakthrough that must always be examined when testing for stability, Even in highly unstable floods, the recovery approaches the maximum if enough pore volumes of displacing fluid are pumped through the system. In view of the approximate nature of the derivation leading to Equation (209), and the inherent complexity of the unstable displacement process it must be expected that the predictions will not agree with the data to a high degree of accuracy, Because of the difficulty of ascertaining the precise velocity at which the minimum breakthrough saturation is reached, (cfo Figure 8 of (29)), the predictions must be limited to an order of magnitude. The amount of experimental data which is available in the literature is very limited, A comparison of the predictions of Equation (209) with these experimental data is shown in Table IV below. TABLE IV CRITICAL VELOCITIES FOR UNSTABLE DISPLACEMENT IN POROUS MEDIA ft ft Displacing Fluid Displaced Fluid Vccalc(s —) Vobs(s-) Source Water (u = o768cp) Oil(u=150cpest).67 x 10-4 o85 x 10- (99)^Fig.l0 Water (u =,895cp) Oil(u - 3o58cp) o46 x 10-2 072 x 10-2 (29),Fig,8 Water (u =.895cp) Oil(u = 22.3cp) o91 x 10-3 75 x 10-3 (29),Figo8 Water (u = 0894 ) Oil(u = 00cpest),23 x 10-3 ~67 x 10-3 (27),Figo3 Water-alcohol Oil(u = 635,) o41 x 10-3,29 x 10-3 (37),Figo12 (u = o644 cp ) and Table II

-182The agreement between the calculated and observed values obtained above is not very close, and yet, when everything is taken into account it is encouraging that there should be agreement even within a factor of about 2. In the first place, the parallel-plate model of porous media predicts a zero critical velocity for the experiments performed above (gravity negligible). Velocities of 103 ft/sec are relatively high for displacements conducted under the usual field operating conditions. It is, therefore, of extreme importance to ascertain the order of magnitude of the critical displacement velocity for real porous media. To know that the critical velocity Vc lies somewhere between zero and values which are considerably higher than occurring in practice is not of much practical use. The real porous medium-fluid system will have a critical displacement velocity somewhere above zero, but the possible values vary over such a wide range that it must be considered as quite satisfactory if the theoretical value comes within a factor of 2 of the experimental value. The fact must be stressed that no claim is made for the strict accuracy of the mathematical development leading to Equation (209). At best, the theory has only been advanced to provide a means of calculating a quantity which could previously not be calculated at allo It must be considered quite encouraging that the theory is in such relatively good qualitative accordance with experiment. That this agreement between theory and experiment is not fortuitous is borne out by range of variables encountered in the experiments cited in Table IV.

-183An example calculation will now be shown to illustrate the use of the method and to indicate some qualitative conclusions. A porous medium of porosity 20% and permeability 500 md contains oil of viscosity 20 cp which is being displaced by water of viscosity 1 cp. The oil water surface tension is 25 dynes/cm. Connate water is 20% and residual equilibrium oil saturation is 20% (.. AS' =.60). The porous medium is strongly water wet. From Equation (209), the critical volumetric velocity Vc is computed to be 162 ft/day. This is a rate far in excess of any normally encountered in the field. It thus can be concluded that the theoretical generalizations obtained here indicate that most real homogeneous, isotropic porous media will be stable to immiscible displacements even at highly adverse mobility ratios, H. Effect of Model Depth on Perturbation Growth During the course of laboratory experiments in which viscous fingers occur (cf, Chuoke, et al, (23)), the influence of the boundaries of the experimental apparatus may be considerable. It is, therefore, reasonable to inquire into the quantitative relationships between such parameters as finger wavelength, fluid properties, and the proximity of the disturbance to the walls of the experimental apparatus, Probblems of this type have been treated previously by, for example, MilneThomson (82) for the case of progressive waves of homogeneous fluid.

-1841. TTheoretical Development The physical model to be examined is illustrated in Figure 34. Assume that the porous medium has thickness 1 + Q2 and that the mean position of the perturbed interface is designated z = 0. The perturbation velocities will be assumed to vanish on z = -21 and z = +22 because of the termination of the porous medium at these planes. Fluid 1 is being pumped into the lower face of the porous slab at -21 at velocity Vo Fluid 2 is leaving the upper face at i2 at velocity Vo The Laplace equations govern the flow in each of phases 1 and 2. i.e., V1 = v2d The general solution to the Laplace equation is 0 = A cosh k(f(z)) + sinh k(fz) ent cos x cos 2y cos x cos where f(z) is a linear function of z. The boundary conditions are: 1. w i(+ ) = w( ) 2b 1 1 2, at = w1 = w2 on interface, Boundary condition 1 is seen to eliminate the sinh term from the general solution and establish the functions f(z) as z - 22 and z + 21

-185z = t2 FLUID 2 z:Q0 FLUID I P, m', Z =-', Figure 34. Schematic Diagram of Linear Displacement of Fluid 2 by Fluid 1 in the Case of Finite Boundaries.

-186The velocity potentials for each phase are then, 1 = A- cosh k(z + n1) ent cos x cos (210a) 1 hk7z+1)e 122 \f2 nt k k D1 = A2 cosh k(z - 22) e cos x cos y (210b) From the boundary condition 2 and the definition of the perturbation as, nt k k b = B cos - x cos 2 y A= -nB A nB (211) k s k sinh k k sinh k2 It is desired to determine the effect of the magnitude of 21 or Q2 on the magnitude of the velocity potential. Specifically we wish to determine the extent of deviation of the solutions given o in Equations (210) (boundaries at finite distance) from those derived for the case of boundaries at infinity in Equations (95). Substituting from Equation (211) in Equations (210) and expanding the cosh termas, 1! = - cosh kz coth k1 - sinh kz e cos g x cos y (212a) 1 k 172 I c 2 2 = + ncosh kz coth k1 -sinh kzent cos x cos y (212b) 2 k _2- x cJs 72 2 It is seen that Equations (212) reduce to Equations (95) when coth kil, and coth kQ2 are near 1.00. For 2 =, coth kQ = 1.090 and for 2 = - coth ki = 1.0038.

-187" 2. Discussion The results obtained above show that, if the model boundaries are one-half a wavelength away from the interface, the solution for finite boundaries differ from the infinitely deep case by only.38%o. The implications for experimental work are clear. The point at which fingers can be allowed to form ought to be at least one-half of a wavelength from the model top or bottom in order to insure that the solutions (Equations (95)) for the infinitely deep case will apply. This point is worth considering carefully, since all of the classical solutions (23), (105) are for infinitely deep systems. However, at least some of the experimental work heretofore presented (23) shows fingers being formed exactly at the model base without any distance intervening. Summary (Sections III A-H) The impulse-momentum equations for non-steady flow of fluids through porous media were developed and used to establish the criteria for the stability of fluid-fluid displacements under conditions where unsteady and inertial effects were significant. The work-energy equations were also developed to demonstrate the influence of viscous damping. The stability criteria were also established for the cases of accelerated, stratified and radial flows as well as for the effects of model boundaries. It was shown that the stability of fluid-fluid displacements in real porous media can only be accounted for by reference to capillary annihilation of the fingers approximate predictrol methods were developed.

IV, RESIDUAL EQUILIBRIUM SATURATION IN POROUS MEDIA When a wetting fluid displaces a non-wetting fluid from a porous material, some of the non-wetting phase is left behind in the pores. In natural gas reservoirs, the aquifer water is the wetting phase and the gas is non-wetting. Also, in petroleum reservoirs, the oil is usually the non-wetting phase, and is to be displaced by water during secondary recovery operationso With the demand for oil and gas reaching higher and higher levels, and the location and exploitation of new reserves becoming more difficult and expensive, a knowledge of the practical limits of secondary recovery operations is imperative, An indication of this practical limit is given by the residual equilibrium saturation concept, Residual equilibrium saturation may be defined as the minimum fraction of the pore space that can be occupied by gas or oil after the porous medium has been flooded to completion with a very large number of pore volumes of water, It is desired to be able to predict a priori the residual equilibrium saturation from the measured properties of the fluids and porous medium. These measured properties are expected to include permeability, porosity, viscosity density, etc,, as well as any other characterization of the pore morphology that may be at hand. To date, no papers treating theory of the magnitude of the residual equilibrium saturation have appeared in the literature, In fact, very little reliable experimental data is available, The phenomenology of -188

-189the residual equilibrium saturation is, as yet, imperfectly understoodo It appears that even qualitative determinations of the parameters affecting trapped gas or oil have not been attempted. This is in large part due to the inherent complexity of the internal structure of the porous material and the fact that residual equilibrium saturation must, apparently, be treated from a microscopic viewpointo This means that the processes actually going on inside a given pore must be examined quantitatively in detail if adequate understanding is to be achieved. The microscopic nature of the gas trapping process appears to rule out exact approaches which are based on measurement of macroscopic quantities, such as permeability, The Darcy approach to porous media (28) allowed the integrated macroscopic equivalent of fluid resistance (permeability) to be established experimentally in order that the considerable mathematical difficulties inherent in the exact treatment of the microscopic fluid motions could be avoided. This macroscopic approach is quite successful as long as macroscopic problems are being treated (pressure or saturation distributions, etco), but cannot be expected to be as easily applicable to the solution of microscopic problems such as that of residual equilibrium saturation, An attempt to achieve an adequate understanding of microscopic processes in porous media from a truly phenomenological viewpoint was made by Smith (113), in connection with the theoretical treatment of capillary pressure, It was established at that time, and subsequently confirmed by Melrose (80), that the ideal soil concept of porous media is the proper point of departure necessary to

-190an adequate understanding of the microscopic mechanics of fluids in pores. Melrose (80) shows, on the basis of the success of the ideal soil model in predicting capillary phenomena, that it is superior to the capillary tube bundle theory of porous media. In addition, the hydraulic radius theory has been severely criticized from the point of view of its usefulness in treating single phase flow (cf, p. 107 of (107))o The essence of the ideal soil concept is that a porous medium is composed of idealized packings of spheres. The porosity is to be varied by varying the packing mode and the permeability by varying the sphere size. The work of Smith (113) and Melrose (80) has provided a qualitative basis on which the fundamental problem of residual equilibrium saturation can at least be stated with some degree of confidence. However, the ideal soil concept is rendered highly limited by the fact that ideal packings are assumed as a model. Real porous media, having enormously complicated microscopic internal boundaries, cannot, therefore, be represented exactly by this approach, In spite of the simplifications involved in the ideal soil concept, it provides a basis on which qualitative conclusions can be drawn. The major conclusion indicated by ideal soil theory is: The fluid-fluid interface jumps unstably across a pore from a position of critical stability. The position to be developed in this thesis is that it is the kinetics of this interface jump which determines the extent of residual equilibrium saturation.

A, Description of Microscopic Flow Processes in Porous Media The formation of a trapped globule of non-wetting phase in the pores of a porous medium is the process by which the residual equilibrium saturation is established, In order to understand the trapping process thoroughly, it is necessary to construct an adequate phenomenology of the microscopic hydrodynamics which occur, It must be emphasized that the microscopic approach involves the examination of processes taking place in individual pores rather than in an assemblage, All the previous treatments of the microscopic aspects of fluid distribution in pores have been concerned with the equilibrium hydrostatics of interfaces (1), (5), (16), (19), (21), (22), (80), (113)o These papers treat the problem of the position of interfaces between immiscible fluids in porous media from the viewpoint that the fluids are not moving and are at thermodynamic equilibrium. The problem of moving interfaces has been largely ignored in the literatureo Melrose (80) has considered microscopic interface hydrodynramics he mentions that interfaces may Jump through a pore in an unstable process. Microscopic interface behavior can most easily be discussed by separating those processes involving equilibrium from those which involve dynamicso lo Microscopic Hydrostatics in Porous Media When two immiscible fluids are at rest in any configuration whatever, the fundamental equation which relates the pressure difference across the interface with the principle radii of curvature is

-192P1 - P = aL2 ( )1 -+ () (213) 2 %2 r1 r2 where rl and r2 are the principal radii of curvature and o12 is the interfacial tension between the fluidso Equation (213) is known as the Laplace equation. Its differential form is, ~2- 2^ Y l p (214) yax y O12 where ~ represents the height of any element of the curved surface from the mean datum plane. Equation (213) describes the geometry of the interface under all conditions of static equilibrium, The second basic concept necessary to an understanding of capillary hydrostatics in porous media is that of the contact angle. The contact angle is the angle which the fluid-fluid interface subtends at the surface of the solid phase which composes the matrix of the porous material. It is a thermodynamic quantity and depends only on the following interfacial tensions (1) phase 1 to solid,' Ols (2) phase 2 to solid, o2s and (3) phase 1 to phase 2, O12, The quantities als and a2s must be considered in the immediate vicinity of the three phase line of contact only. The quantitative relationship between these interfacial energies and the contact angle is given by Young's equation 0ls - 02s cos @ =-s - (215) ~12 In real systems, contact angles must be measured experimentally, since there seems to be no way of obtaining,cra or a2s independently.

-193Also, these quantities are known to be strongly dependent on the adsorption of small amounts of extraneous components on the rock surface, as well as to the crystalline form of the solid. The third basic foundation of equilibrium hydrostatics is contained in Kelvin's theorem (123). According to Kelvin, if two similar masses of liquid, confined by identical boundaries, are at equal distances above a free liquid surface, then the vapor pressures over both of the surfaces must be the same if the system is at equilibrium. The basic ideas behind Kelvin's theorem can be outlined as follows. Referring to Figure 35(a), it is noticed that if a thin capillary tube, open at both ends, is placed in a liquid contained in a larger sealed vessel, the height to which the liquid will rise is given by Laplace's equation as 2c12 h = 12 Rt(PQ - Pg) where Rt is the capillary tube radius and pg and p. are the gas and liquid densities. It is assumed that the gas is composed of the vaporized liquid (i.e., the water-steam system, for example). r1 and r2 have been equated to Rt, which is reasonable for small diameter capillary tubes. The immobility of the system in the final state, implies that it is in physical and thermodynamic equilibrium. The fact that the vapor possesses a hydrostatic head implies that the vapor which is in equilibrium with the curved liquid surface in the capillary tube differes from the vapor pressure in equilibrium with the plane liquid surface by approximately -p h (neglecting gas g

STEAM STEAM ^~h 4 ^ PRESSUREPo ___ _ W R- WATER.,-WATER~ a b Figure 35. Capillary Condensation in Thin Tubes.

0ompressibility)0 Suppose a new experiment is carried out, as shown in Figure 35(b). A capillary tube of the same radius Rt. this time closed at the bottom end is placed in the liquid in the larger vesselo Assume that a small amount of liquid is in the bottom of the tube, and that the mean liquid level in this tube is below that corresponding to equilibrium in the open bottomed tubeo The radius of curvature of the liquid surface in the case of Figure 35(b) is still Rto However9 the gas exposed to its surface is at a pressure greater than po-p ho The system is therefore not at equilibrium, and vapor will condense in the closedbottomed tube until the liquid height in it reaches ho This phenomenon is known as capillary condensation or Kelvin condensation, A corollary of the Kelvin theorem is that fluid-fluid interfaces must have the same curvature if the vapor pressure over each of them is equal, and if the system is at equilibrium, Kelvin's theorem only treats the equilibrium situation, and not the kinetics with which it is achieved~ The foundation of all microscopic hydrostatics in porous media are thus (1) Laplace s equations (2) Young's relation, and (3) Kelvin's theorem, The Laplace equation determines the surface in space defined by the fluid-fluid interface, The contact angle is the boundary condition that must be satisfied at the three phase line of intersection between the solid and two fluids, Kelvin s theorem allows determination of the curvatures of surfaces which are in equilibriumo It must be noted that the Laplace equation is rigorously true only for equilibrium conditions,

-196The use of microscopic equilibrium hydrostatics has allowed the successful prediction of hysteresis in capillary pressure on the basis of the ideal soil model (80), (113), The success of the method in treating equilibrium capillary pressures has served to establish the ideal soil model as a reasonably good working model for the interpretation of phenomena in porous materials, 2, Microscopic Hydrodynamics in Porous Media Equations of Motion. On a microscopic scale, the fluid-fluid displacement process must be governed by the Navier-Stokes equations. Due to the presence of the pore walls and the large velocity gradients generated, an exact treatment of microscopic dynamics cannot be based on potential flow, Non-steady Capillary Condensation, The Kelvin theorem, as outlined previously indicates that mass transfer will occur between fluid-fluid interfaces when the system is not at equilibrium (57)o Specifically, the example of water displacing gas can be used as an illustration. The connate water which is found in the zone which is not flooded out is present in the form of pendular rings, Reference to Figure 38 shows the physical situation, The advancing water-gas interface is supposed hemispherically shaped and located at AB, It advances to successively assume positions 1, 2 and 3, The gas in the pore will, of course, have a certain partial pressure of water corresponding to the existing pressure and temperature. Thus, both the pendular ring and the main interface at AB are both in contact with the same gas, Kelvin's theorem then demands that the interface

curvatures be the sarme in both the pendular rings and the main interface, Since the main interface curvature decreases from position 1 to 3, then the pendular ring interface curvature must increase in exactly the same way if the equilibrium conditions are attained, The only way that the pendular ring curvature can decrease is for the pendular ring to grow. This is accomplished by capillary condensation according to Kelvin s theorem, The pendular ring grows by means of water evaporation from the main interface at AB, diffusion through the gas and condensation on the pendular ring, Under equilibrium conditions, the pendular rings grow and the main interface advances until they touch and coalesce at some point. In dynamic displacements in porous media, equilibrium conditions may not occur, If the process of evaporation-diffusion-condensation from the main interface to the pendular ring is slower than the gross motion of the main interface, then non-equilibrium conditions will occur, The curvature of the pendular ring is then not the same as that of the advancing interface, The size of the pendular ring is expected to influence the point of coalescence with the main interface, The residual equilibrium saturation is in turn affected by this point of coalescence, The conclusion can therefore be drawn that the kinetics of the capillary condensation phenomenon are a factor in determining the magnitude of the residual equilibrium saturation, Wettability and Adsorption Kinetics, An important basic difference between microscopic hydrostatics and hydrodynamics in

-198porous materials is the rate of adsorption of extraneous compounds on the solid surface, Also involved is the problem of advancingreceding contact angles and the establishment of the wettability condition, The latter problem has been previously discussed in Section IIA. The contact angle is known to be affected by the adsorption of extraneous compounds on the solid surface (30). However, the kinetics of this adsorption process has received little attention, In a qualitative sense, it must be expected that the kinetics rather than the statics of the wettability phenomenon are the primary considerations in the cases where equilibrium does not prevail, Specifically, it is necessary to know whether or not the contact angle is a function of the rate of advance, and if so, what is that function. It is also necessary to have some idea of the rates of adsorption of extraneous surface active components in relation to the rates of gross fluid flow, 35 Experimental Procedures One of the primary considerations in treating the dynamic motion of fluid-fluid interfaces through pores of a porous material is whether or not the motion is stable or not, The stability to be considered here is that of fluid motion inside a single pore - not the stability of large macroscopic fingers treated in Section III, It does not seem that the magnitude of the residual equilibrium saturation has, as yet, been attributed to the dynamics of microscopic fluid motion,

-199In order to visually observe the formation of the non-wetting phase residual, the following experiment was performed, Glass beads, 6 mm, in diameter were saturated with benzene dyed with Sudan IV red dye, The beads and glass containing vessel were thoroughly cleaned with dichromate solution before use. Distilled water was used to displace the benzene vertically upward from the bead packo During the displacement process the fluid-fluid interfaces were carefully observed. The results of this qualitative experiment were as follows, Any given pore is either entirely full of non-wetting phase (except for connate water) or it is all full of wetting phase. It can thus be said that any given pore is in either a 0 or 1 state, The movement of the fluid interface through the pore takes place in the form of a jumpo This jump is, apparently, between two positions of critical stability, When a fluid-fluid interface jumps through a given pore, the nonwetting phase is completely displaced from that pore, The residual equilibrium saturation appears to be the result of a given pore being surrounded by flooded out pores and then being isolated by a final pinching off process, The resulting residual benzene structure consisted of isolated globules of various shapes, The pores that had experienced an interface jump were completely flooded out, Those that had been surrounded and cut off by the pincer movement of water in the surrounding pores contained an irregularly shaped benzene globule. Thus, the residual equilibrium saturation appears to be the result of the surrounding of pores rather than in anything inherent in the mechanics of the unstable jump across a given pore.

-200If the residual equilibrium saturation were dependent on the interface jump itself we would expect to find a non-wetting phase globule in every pore through which an interface had jumped. Since just the opposite of this is true, the residual can be ascribed to the surrounding process, A given pore can be surrounded and blocked off when it is encircled by smaller pores. The capillary forces draw the displacing water into the smaller surrounding pores and pinch off the larger pore before the applied potential is strong enough to cause displacement in the larger poreo The various stages of this blocking-off process are illustrated in a two dimensional schematic diagram in Figure 36, Chatenever (20) observed that oil globules extending over from 1 to 10 pores were formed when oil was displaced by water from water wet micromodel bead packs, The surrounding process can thus extend over more than 1 pore diameter. Residual equilibrium saturation is thus attributable to a statistical nonhomogeneity in the porous medium which leads to the possibility that any given pore can be more or less completely surrounded by smaller pores. These nonhomogeneities cause the wetting phase to rise around the pore in question due to increased local capillary suction, and ultimately pinch off that pore, The size and shape of the trapped globule will be expected to depend on the kinetics of the interface jump, The question naturally arises as to whether a perfectly homogeneous porous medium will exhibit different behavior with respect

-201TIME 3 TIME 2 - WATER TIME I SAND Figure 36. Schematic Diagram Illustrating the Formation of Residual Saturation by the Blocking Process.

-202AIR OUT LIZ AIR -I WATER GLASS BEAD WATER IN Figure 37. Micromodel Packed with 6mm. dia. Glass Beads Simple Cubic Packing.

-203PENDULAR RING GAS PARTICLE S WATER Figure 38. Invading Interface and Pendular Rings During Displacement of Gas by Water.

-204to residual equilibrium saturation than the more randomly constituted materials. The most obvious example of a perfectly homogeneous porous medium is that formed from uniform spheres in simple cubic packing. Since all of the pores are of exactly the same size and shape, there would appear to be no chance for a residual saturation to form by means of the fluid bypassing mechanism. In order to test the hypothesis that the residual saturation is caused by local nonhomogeneities, a micromodel of a perfectly uniform porous medium was built as follows. A lucite cube of inside dimensions 24 x 24 x 24 mm. was packed with uniformly sized 6 mm. glass beads in simple cubic packing. Inlet and outlet ports were 1/64" in diameter and centered at the top and bottom of the model. Needle valves were screwed into the top and bottom of the model and used to seal it. The construction of this micromodel is illustrated in Figure 37. The pores were thus oriented in vertical columns and horizontal rows. Colored benzene was introduced into the pores of the model completely saturating the pore space. Water was then introduced into the bottom inlet port and displaced the benzene upward. The gross flow was exactly normal to the horizontal planes of beads. The interfaces were observed to jump across the pores in all cases. The benzene-water interface flooded out all the pores in the model starting from the lowest plane of pores near the base, and progressing upward. It was observed that all of the pores in any given horizontal plane were displaced before any pore in the horizontal plane next above was displaced. The residual equilibrium

-205saturation occurring in the micromodel at all flow rates between 1 cc/hr and 224 cc/hr was exactly zero. Every pore in the micromodel experienced an interface jump and hence was completely displaced. The same results were obtained with the air water system, The experimental results obtained above allow the conclusion to be reached that a perfectly uniform packing of spheres will give zero residual saturation, and that the local inhomogeneities in real systems are the cause of residual saturations therein. In a random packing of the same 6 mm. glass spheres, a residual equilibrium benzene saturation of about 2-3% was obtained. Since the packing mode was the only variable between the micromodel and the random pack, the pore size distribution and orientation must be the cause of the residual saturation, B. Theoretical Aspects of Microscopic Motion As stated in the previous section, the Navier-Stokes equations govern the flow of both the displaced and displacing phases in their microscopic aspects. Laplace's equation (Equation (213)) relates the pressure difference across the fluid-fluid interface to the curvatures, The governing equations are thus, tui -l + U I aP x.. - — i _' 1 ui x Xi i Xi il U(26) St m2 ax+ u u2 i a2 + Xi + P2 vu2 U (216b) Bt j 2 px P where subscripts 1 and 2 signify displacing and displaced phases, *For simple cubic packings of very small spheres (i.e. > lmm) this result has not been definitely established.

-206= = respectively. ui and p denote the actual microscopic point velocities and pressures in the pores. Suppose that the position of the interface is designated fi(x,yz,t)=o and that the surface in space defined by the tortuous pores is denoted G (x,y,z)=o,fi is as yet unknown. The boundary conditions are: 1 uil = ui2 = on G (x,y,z) = o 30 T =T Uil = ui2 on f (xpy"z~tj = ijl ij2 Dfi = Ui: -- on f (x,y,z,t) = o uil =Ui2 =Dt The initial conditions are: fi (x)y,zo) = o and u. = Boundary condition 1 states the non-slip condition at the pore walls, Boundary conditions 2 and 3 express the continuity of velocity and stress at the fluid-fluid interface, Boundary condition 4 is the kinematical condition required by the existence of the interface (cf. Streeter (118), po16)o It is seen that Equations (216) represent six simultaneous non-linear partial differential equations, The fact that the interface is moving due to the gross displacement caused by fluid injection means that the problem is one involving a moving boundary, Also, since the fluid-fluid interface is not constrained, the surface of separation is a free boundary,

207 The theoretical aspects of microscopic fluid motion can then be summarized as follows, The hydrodynamic problem to be solved is that involving six simultaneous non-linear partial differential equations coupled with a movring free boundary, In addition, rigorous treatment involves the consideration of wettability and diffusional kinetics, which must, of course, also be coupled to the purely mechanical aspects. Co Prediction of Residual Equilibrium Saturation lo Empirical Correlation of Data Cursory examination of the theoretical problem just outlined shows that the possibility of evolving a general theory of microscopic multiphase motion in porous media must be viewed with pessimism, Solutions to simultaneous non-linear equations involving moving free boundaries seem to be well beyond the power of any method that is presently available0 Even in the case of motion in ideal packings of spheres, the mathematical difficulties are virtually insuperable. That the interface motion has been experimentally observed to occur in a series of unstable jumps is yet another complicating factor, The main difficulty encountered in a consideration of the exact solution revolves around the boundary conditions, Muskat (85), page 184, in considering the problem of single phase flow in real porous media makes the following statements. "The myriads of interstitial and multiply connected grain surfaces are obviously beyond the power of quantitative description,"1 "And even if by some miraculous insight the detailed geometry of the almost infinite variety of

-208these intertwining labyrinths could be crystallized into a quantitative description, it would be utterly futile to attempt to find solutionsoo that would be amenable to adjustment for satisfying the boundary conditions at these surfaces," "It will thus be clear that an approach to the treatment of fluid flow through porous media by the classical method of hydrodynamics must be condemned from the start as being totally devoid of any hope of success," Illustrations of the complexity of actual porous media are given in (117)o Attempts to predict the permeability of porous materials for single phase flow based on general physical considerations have not met with success (cf. Scheidegger (108), p. 89)o Kozeny (68) in an attempt to experimentally verify the theory previously developed by him found discrepancies of -69 to +86 percent between theory and experiment. Most of the other attempts to predict permeability involved the introduction of "factors," which are nothing more than a convenient means of bringing the data into line with the predictions (128), (138). Multiphase flow is at least an order of magnitude higher in complexity than single phase flow, The failure of any attempt to obtain a general theory of single phase flow, together with the inherent greater complexity of the relations (cf, Equations (216) ffo) governing multiphase flow allows the conclusion to be drawn that the microscopic approach to residual equilibrium saturation is not likely to be successfulo Putting aside all questions of agreement of theory with experiment,

-209it would'be highly unlikely that meaningful results of any kind could be extracted from the boundary value problem of Equations (216) ffo The next simplest approach to the prediction of residual equilibrium saturation would seem to be through the use of the macroscopic measured characteristics of the porous material such as porosity, permeability. capillary pressure, relative permeabilityo It also might'be assumed that the fluid properties, such as viscosityr density, interfacial tension) etco are significant variableso This approach might be compared to that used in thermodynamics, where desired properties of materials are computed from known relations between the various state functions, For example, the effect of pressure and temperature on heat capacities can be calculated indirectly by measuring the over-all behavior of the system when other state functions are varied'Unfortunately'there are not presently available any theoretical relations connecting any of the characteristic properties of porous media (such as porosity) permeability) capillary pressure, relative permeability, etCo ). It is for this reason that the measurement of these over-all macroscopic quantities (black box technique) is not presently capable of yielding any way of theoretically predicting the microscopic residual equilibrium saturation, The formidable theoretical difficulties involved in the determination of residual equilibrium non-wetting phase saturation from general considerations indicates that an empirical correlation must be obtained. Accordingly, experimental data from the literature covering a wide range of parameters were examined. Careful analysis

-210of these data showed that, to a good first order of approximation, the non-wetting phase saturation was a function only of the porosity of the porous material. The plot of the data is shown in Figure 390 The tabular compilation of the data is given as Table V, Appendix Bo The trend of the data appears to follow the curve, Residual Equilibrium Saturation = o62 -.13 ~ (217) The average deviation of the 60 data points from this mean curve is 204% of the pore space. The average absolute percentage deviation is 9. 4o% For example, an observed residual of 33% vso a prediction of 30%) results in an error of 3% of pore space and an absolute error of 10%.o The data plotted show residual equilibrium saturation when the non-wetting phase has been displaced by the wetting phase. The data include cases in which gas is the displaced phase and also cases where oil is displaced. Water or brine is usually the displacing phase. These data cover a very wide range of porous media properties, The materials whose porosities are above 30% were mostly unconsolidated sand or glass bead packs. The materials with porosities less than 30% were sandstones, Limestone data were eliminated from the correlation because of the erratic nature of the results caused by the vuggy porosities encountered. The initial water saturation in all the data in Figure 39 corresponded to connate water or lesso The viscosity ratios (viscosity of displaced phase divided by viscosity of displacing phase) varied from about.01 for gas-water displacements to about 10 for some of the oil-water displacements0 Permeability varied in the range of about 10 md to about 100,000 mdo

-21152.. 0 50 48 46 426 0 ~ 0~0 ~~ ~ ~~ 44 0 42 No 0 40 38 ~~~~ -0 ~ ~NST~~~~ 0 0 34 0 \ 00 z 0 o ~G3 I I I I i W o o D A P t 28 o):LLJ _ C 22 20 RE SIDUAL STU RAT I ON _ 00 1 2 iz O 1 20 0 ^ ~ ~ \ RESIDUAL SATURATION ~~cn^J) I ( FOR UNCONSOLIDATED POROUS MEDIA Q Q\ QT AND VARIOUS SANDSTONES 0 \ 18 WETTING PHASE IS DISPLACING PHASE 0 SEE TABLE v, APPENDIX A. 16 14 0 12 10 0 004 0.08 0.12 0.16 0.20 0.24 028 032 0.36 040 POROSITY, (FRACTION VOIDS) Figure 39. Residual Equilibrium Non-Wetting Phase Saturation vs. Porosity for Unconsolidated Porous Media and Various Sandstones.

-212Flooding velocities were not always specified, but may be said to have varied over a wide range. Surface tensions varied in the range of 20 to 70 dynes/cm. The effect of wettability could not be included in the correlation because of the lack of availability of experimental data. All the materials on which data was taken for Figure 39 can be classed as water wet. In all of the data shown in Figure 39, the permeability, viscosities, flooding rate, and surface tension were carefully examined and found to have no systematic effect on the residual equilibrium saturation, At adverse viscosity ratios, the residual equilbrium saturation data must be interpreted with extreme care because of the possibility of viscous fingering. This situation is clearly illustrated in Figure 2. The water fingers are here seen moving through the viscous oil and breaking through the outflow face of the model. Even after 6.5 pore volumes of water have been injected, the oil recovery is only 52%. In stable floods, 6.5 pore volumes is usually much more than enough to reach the true residual saturation. At high rates and highly adverse mobility ratios, it is possible to have the viscous fingers form and persist for a long time after breakthrough, as shown in Figure 2. It may even be possible to have them persisting indefinitely, or at least until very high water-to-oil producing ratios are attained. In the latter case, the displacing water is simply moving through the water fingers from inlet to outlet of the core and not performing any displacement. The final oil saturation will not reflect the true residual under these conditions - it will simply reflect the extent of the instability

-213of the system, Thus, observation of an apparent "effect of rate and viscosity ratio" may not reflect a meaningful determination of the residual saturation at highly adverse viscosity ratioso Thus, all data in which it could reasonably be expected that viscous fingering occurred (7 > 10, and high rates) were eliminated from the correlation. Also, the data of Searcy and Stewart (110) were not used because their experiments were done by capillary imbibition rather than frontal displacement. These data are not in agreement with the correlation of Figure 39, Data taken by Morse and Yuster (84) and Richardson and Perkins (101) showed no effect of rate on residual. Also, another study (62) yielded no effect of surface tension in a series of systematic experiments, 2. Experimental Procedures In order to determine the effect of water flooding rate on the recovery of gas from porous media at very low rates, the following experiments were performed. Four consolidated and two unconsolidated porous materials were mounted and flooded at various rates with distilled water, Air was the displaced fluid in all cases9 and the cores were at room temperature and pressure except for the pressure drop through the core, After each run, the cores were desaturated by displacement with air, The properties of the porous materials used are given in Table VIo All cores and core materials were fired at 200~C for 24 hours before use, The volume of the storage chambers in the Ruska pump was 500 cCo The cores were mounted vertically and connected

-214to the pump outlet by about 1 foot of 1/8 inch steel tubing. All air was removed from the pump outlet chamber and the core inlet chamber before the beginning of each run. The cores were flooded in a vertical position with the water injected into the base, Total volume injected was read from a calibrated dial on the pump. TABLE VI PROPERTIES OF POROUS MATERIALS USED IN RATE EXPERIMENTSCore PermeDesig- Porous Core abilnation Material Dimensions Porosity Porosity ity (inches) (Vaco Imbibo ) (True Density) (md) A Sintered 1.48 x 12.02 0.243.259 85 3 Alundum B Sintered 1,98 x 11,97 0.249.243 303 Alundum C Sintered 1.97 x 11,77 0,238.246 227 Alundum D Berea 2.03 x 2.08 x 11,10 0o180 o191 114 Sandstone E Ottawa Sand, 1o 17 x 10o15 o 409 404 11700 150/200 mesh F Glass Powder 1.17 x 10.15 o 419.411 1160 The consolidated cores were mounted according to the technique outlined in Appendix CO The mounting technique for the unconsolidated materials is described in Appendix D. The water was injected into the base of the vertical cores at rates varying between 1 cc/hr and 128 cc/hr a precision of + 1% with a piston displacement pump manufactured by the Ruska Instrument Corporation of Houston, Texaso All fittings were of 1/8" steel tubing, A photograph of the Ruska pump is shown in Figure 42,

::::::j::::if f:::~:: ~::~~~~~~~~~~~~~~~~::? m! W F ~~~~~~~~~~~~~~~~~~~~;::fi.i: ti'g I e: a,, ji::: Figure 42.: Ruska Pump ith Mounted Core.:: i:' Ii ti ti i0 /~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~~~~i~ -su:r -- & dj.............:... *Z K mL1~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~iiiiiii.::i Figure 42. Ruska Puinp with Mounted Core.~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i.: —-i~:~i.iii:l.

-216The porosities of the sintered alundum cores were measured in two different ways. The first method consisted of drawing deaerated water into the bone dry core at about 0.1 mm of pressure by means of a vacuum pump. The change in weight gives the interconnected pore volume. The second method of determining the pore volume of the alundum cores consisted of calculations based on the core weight, its over-all dimensions, and the true density of alundum. The two methods checked as indicated on Table VI. Residual equilibrium saturations were obtained by weighing the core at various times after the throughput of increasing numbers of pore volumes of displacing water. The residual saturation was calculated from a knowledge of the total pore volume, the weight of the completely dry core, and the weight after displacement of the air by water. The magnitude of the residual was cross-checked by closing off the outlet end of the core and injecting a specified volume of water into it with the Ruska pump. The rise in pressure per unit volume of water injected allowed calculation of the total volume of gas in the core. The residual saturation calculated by weight difference checked quite closely with that calculated by the pressure rise procedure, as is shown by the sample data given in Table VII below. No air production after water breakthrough was observed on any of the experimental runs performed. Example data are given in Table VIII below. Apparently, the low viscosity of the air and the relatively large capillary end effect cause retention of the water until the core achieves maximum saturation, at which time breakthrough occurs.

TABLE VII COMPARISON OF RESIDUALS COMPUTED BY WEIGHT DIFFERENCE AND PRESSURIZATION Alundum Core, K = 85 md, Volumetric Velocity 257 ft3/ft2yr Volume of gas measured by weight difference = 14o95 cco Pressure rise per volume water injected =,o668 atm/cc. Volume of gas measured by pressurization = 14 9 cco Alund-um Core, K = 303 md, Volumetric Velocity 4672 ft3/ft2yr Volume of gas measured by weight difference = 2674 cc, Pressure rise per volume water injected =- 0382 atm/cco Volume of gas measured by pressurization 26,2 cco Alundum Core, K = 227 md, Volumetric Velocity 18688 ft3/ft2yr Volume of gas measured by weight difference = 26,57 cco Pressure rise per volume water injected = o0370 atm/cco Volume of gas measured by pressurization = 27,0 cc,

TABLE VIII WEIGHT OF CORES AFTER WATER BREAKTHROUGH Alundum Core (K = 85 3 md)(2056 ft3/ft yr) 150-200 Mesh Sand Core (3088 ft3/ft2yr) cc Water Injected Core Weight (gm) cc Water Injected Core Weight (gm) 113 0 2660, 62* 112, 8 846, 47 131, 1 2660,39 130,9 846 09 139o0 2660,48 138 8 845090 158,4 2660053 158 3 846,06 176 4 2660 50 176,2 846 44 * Breakthrough point,

30 25 z 1 S 20 - - - ~~ o 0 4 o D 5 < 0 ^ 13 0 _ < 10 w cr; ALUNDUM CORE, k=227 md, 2.238 II - -A |El ALUNDUM CORE,k=85.3mdd,.243 0 150-200 MESH SAND k 11700 md, ~409 0 2 3 4 5 10 10 10 10 VOALUMETRIC VELOCITY (ft2 -year) Figure 45. Effect of Flooding Rate on Non-Wetting Phase Residual Equilibrium Saturation. Residual Equilibrium Satburation.

30 25 z 0!R 20 20 ~~~El El 0 ]~ (/) Q 15 r w 10._ w 5 0 BEREA SANDSTONE, k 114 md., =.180 o ALUNDUM CORE, k= 303 md., =.249 0 GLASS POWDER, k=160 md., =.6419 10 10 10 VOLUMETRIC VELOCITY (ft3/f-yeor) Figure 44. Effect of Flooding Rate on Non-Wetting Phase Residual Equilibrium Saturation.

-22150 40 Lu L 30 rY oI 0^ 0 CORE Ia 2r 0 CORE 2a 0 CORE 3 a 20 0 20304 CORE4a ~0 -J I0 0 __ ~~~~ 0 10 20 30 40506070 INITIAL WATER CONTENT, Swi Figure 4~. Effect of Initial Water Content on Residual Oil After Faterflooding. Data of Dykstra and Parson's (35).

-222The results of the experiments are given in Figures 43 and 44 and summarized in Appendix E. The data show conclusively that the rate of displacement has no effect on the residual gas at the rates employed. The residuals obtained on the consolidated materials are somewhat lower than would be predicted by Equation (217) because the initial water saturations at the start of the floods were somewhat higher than the connate water. The connate water content of sintered alundum and Berea sandstone was found by Kyte, et al. (69) to be 15.4% and 26.00, respectively. The data appearing in the correlation of Figure 39 were all taken on cores in which the initial water saturation was at connate water or less. Dykstra and Parsons (35) obtained some data showing that initial water saturation may have an effect on the residual oil retained after water flooding. Some of their data is shown plotted in Figure 45. Cores 3a, 4a and 6a show a decreased residual with increasing initial water saturation while cores la and 2a show no such effect. The results are thus somewhat contradictory. It may be concluded, however, that the initial water content can affect the residual under certain conditions. This is, apparently, what occurred in the experiments on the consolidated cores whose residuals are given in Figures 43 and 44. These figures show, however, that, whatever the residual saturation, it is not influenced by the displacement rate. Summary As a result of the work of this thesis and also of work previously done (21), (22), (29), (80), (136), the following conclusions

-223can be drawn regarding various aspects of multiphase flow in porous media. 1o The fluid-fluid interface moves through a given pore in a highly unstable manner, The velocity of this jump is many times larger than the gross fluid velocity through the porous medium, and is controlled largely by capillary forces, 2, A pore is either completely empty, with the exception of connate water, or completely fullo The transition between these states is effected by the interface jump. 5o Any pore through which an interface jumps is completely swept of non-wetting phase fluido 40 The source of the non-wetting phase residual saturation is in a pore blocking movement, in which a given pore is surrounded by other pores which have been completely flooded. outo The given pore thus never is swept by an interface jump and retains its original non-wetting phase contento 50 The pore blocking process is due to a differential capillary suction process, in which the smaller pores surrounding a larger pore imbibe water to a greater extent. If the small local non-homogeneous character of the porous medium penrmits certain pores to be so surrounded, then these smaller pores will be flooded out first by capillary suction, and will pinch off the

-224large poreo The local pore size non-homogeneity is indicated as the cause of non-wetting residual saturation because of the zero residual saturation observed when displacements on a micromodel with spheres in perfect cubic packing were performedo 6 Empirical correlation of all available data shows that non-wetting residual is, to a first order of approximation, a function only of the porosity of the porous material~ The rate of displacement of gas by water has no effect on the residual gas saturation in the range of flow rates studied,

BIBLIOGRAPHY 1. Adam, N. K., Disc. Faraday Soc., No. 3, p. 5 (1948). 2. Amott, E., Trans. AIiME. (Petr.), 216, 156 (1959). 3. Archie, G., Trans. AIME. (Petr.), 146, 54 (1942). 4. Bass, D. M. and Crawford, P. B., Trans. AIME. (Petr.), 207, 293 (1956). 5. Bethel, F. T. and Calhoun, J. C., Trans. AIME. (Petr.), 198, 197 (1953). 6. Blair, P. M. 4th Biennial N. Texas Sect. AIME., Secondary Recovery Symposium, Wichita Falls, May 2, 1960. 7. Boyer, R. L., Morgan, F. and Muskat, M., Trans. AIME. (Petr.), 170, 15 (1947). 8. Breston, J. N. and Hughes, R. V., Trans. AIME. (Petr.), 186, 100 (1949). 9. Brown, H. W., Trans. AIME. (Petr.), 192, 67 (1951). 10. Brownell, L. E., Chem. Eng., 56, 112 (1949). 11. Brownell, L. E. and Katz, D. L., Chem. Eng. Prog., 43, 601 and 703 (1947). 12. Bruce, G. H., Oil and Gas J., 46, 223, July 26, 1947. 13. Buckley, S. E. and Leverett, M. C., Trans. AIME. (Petr.), 142, 107 (1942). 14. Buzinov, S. N., Dokl. Acad. Nauk., App. Phys. Sect., 116, 321 (1957). 15. Calhoun, J. C., Fundamentals of Reservoir Engineering, Univ. of Oklahoma Press, Norman, Okla. (1953). 16. Calhoun, J. C., Trans. AIME. (Petr.), 198, 197 (1953). 17. Cardwell, W. T., Trans. AIME. (Petr.), 216, 271 (1959). -225

-22618. Carter, R. D. and Tracy, G. W., Trans. AIME. (Petr.), 219, 195 (1960). 19. Cassie, A. B., Disc. Faraday Soc., No. 3, p. 11 (1948). 20. Chatenever, A., American Petroleum Institute Research Project 47B, 1957. 21. Chatenever, A., Trans. AIME. (Petr.), 195, 149 (1952). 22. Churchill, R. V. Operation Mathematics, McGraw-Hill Book Company, Inc., New York (1958). 23. Chuoke, R. L., Van Meurs, P. and Van der Poel, C. Trans. AIME. (Petr.), 216, 188 (1959). 24. Coberly, C. A. and Marshall, W. R., Chem. Eng. Prog., 47, 141 (1951). 25. Coomber, S. E. and Tiratsoo, E. N., Inst. of Petr., 36, 543 (1950). 26. Craig, F. F., Trans. AIME. (Petr.), 195, 325 (1952). 27. Croes, G. A. and Schwarz, N., Trans. AIME (Petr.), 204, 35 (1955). 28. Darcy, H., Les fontaines publiques de la ville de Dijon, Dalmont, Paris (1856). 29. De Haan, H., 5th World Petroleum Congress, Sect. II, Paper No. 25 (1959). 30. Denekas, M. 0., Mattax, C. C. and Davis, G. T., Trans. AIME. (Petr.), 216, 330 (1959). 31. Dombrowski, H. S. and Brownell, L. E., Ind. Eng. Chem., 46, 1207 (1954). 32. Douglas, J., Peaceman, D. W. and Rachford, H. H., Jr., Trans. AIME. (Petr.), 216, 297 (1959). 33. Douglas, J., Trans. AIME. (Petr.), 213, 96 (1958). 34. Dunlap, H. F., Bilhartz, H. L., Shuler, E. and Bailey, C. R., Trans. AIME. (Petr.), 186, 259 (1949). 35. Dykstra, H. and Parsons, L. R., Secondary Recovery of Oil in the United States, p. 160, American Petr. Inst., New York (1950).

-22736. Engelberts, W. and Klinkenberg, L., 3rd World Petroleum Congress, 2, 544 (1951). 37. Everett, J. P., Gooch, F. W., Jr. and Calhoun, J. C., Trans. AIME. (Petr.), 189, 215 (1950). 38. Fayers, F. J. and Sheldon, J. W., Trans. AIME. (Petr.), 216, 147 (1959). 39. Fatt, I. and Klikoff, W. A., Jr., Trans. AIME. (Petr.) 216, 426 (1959). 40. Friedlander, S. K., A.I.Ch.E.J., 3, 43 (1957). 41. Geertsma, J., Croes, G. A. and Schwarz, N., Trans. AIME. (Petr.) 207, 119 (1956). 42. Geffen, T. M., Parrish, D. R., Haynes, G. W. and Morse, R. A., Trans. AIME. (Petr.), 195, 29 (1952). 43. Geffen, T. M. and Owens, W. W., Trans. AIME. (Petr.), 192, 99 (1951). 44. Gillespie, T., J. Colloid Sci., 14, 123 (1959). 45. Goldstein, S., Modern Developments in Fluid Mechanics, Oxford, Clarendon Press, Vol. II, p. 641 (1938). 46. Graham, J. W. and Richardson, J. G., Trans. AIME. (Petr.), 216, 377 (1959). 47. Graton, L. C. and Fraser, H. J., J. Geol., 43, 785 (1935). 48. Handy, L. L., Trans. AIME. (Petr.), 219, 75 (1960). 49. Haruni, M., Chem. Eng. Sci., 2, 203 (1953). 50. Hawthorne, R. G., Trans. AIME. (Petr.), 219, 81 (1960). 51. Heller, J. P., Rev. Sci. Inst., 30, 1056 (1959). 52. Holmgren, C. R. and Morse, R. A., Trans. AIME. (Petr.), 192, 135 (1951). 53. Holmgren, C. R., Petroleum Technology, 11, 1, July, 1948. 54. Housel, W. S., Univ. of Michigan, Eng. Res. Inst., Rep. Ser. No. 16, p. 465.

-22855. Irmay, S., Trans. Am. Geophys. U., 35, 463 (1954). 56. Johnson, E. F., Bossier, D. P. and Naumann, V. O., Trans. AIME. (Petr.), 216, 370 (1959). 57. Jones-Parra, J. and Calhoun, J. C., Trans. AIME. (Petr.), 198, 335 (1953). 58. Jones-Parra, J., Trans. AIME. (Petr.), 198, 314 (1953). 59. Josendal, V. A., Sandiford, B. V. and Wilson, J. W., Trans. AIME. (Petr.), 195, 65 (1952). 60. Kaplan, W., Advanced Calculus, Addison-Wesley Company, Reading, Mass. (1952). 61. Katz, D. L., ed., Handbook of Natural Gas Engineering, McGrawHill Book Company, Inc., New York (1959). 62. Kennedy, H. T. and Guerrero, E. T., Trans. AIME. (Petr.), 201, 124 (1954). 63. Kinney, P., World Oil, 134, 174, March, 1952. 64. Kimbler, O., Oil and Gas J., 55, 85, December 16, 1957. 65. Klute, A., Soil Science, 73, 105 (1952). 66. Knudsen, J. G. and Katz, D. L., Fluid Dynamics and Heat Transfer, McGraw-Hill Book Company, Inc., New York (1958). 67. Kocatas, B. and Cornell, D., Ind. Eng. Chem., 46, 1219 (1954). 68. Kozeny, J., Wasserkr. u. Wasserw., 22, 86 (1927). 69. Kyte, J. R., Stanclift, R. J., Stephan, S. C. and Rapoport, L. A., Trans. AIME. (Petr.), 207, 215 (1956). 70. Kyte, J. R. and Rapoport, L. A., Trans. AIME. (Petr.), 213, 423 (1958). 71. Laird, A. D. K. and Putnam, J. A., Trans. AIME. (Petr.), 192, 275 (1951).

-22972. Lee, E. H. and Fayers, F. J., Trans. AIME. (Petr.), 216, 284 (1959). 73. Leverett, M. C., Trans. AIME. (Petr.), 142, 152 (1941). 74. Leverett, M. C., Lewis, W. B. and True, M. E., Petroleum Technology, 5, 1, January 1942. 75. Levine, J., Trans. AIME. (Petr.), 201, 21 (1954). 76. Lewis, D. J., Proc. Roy. Soc. A., 202, 81 (1950). 77. Marx, J. W., Trans. AIME.(Petr.), 207, 88 (1956). 78. McCabe, W. L. and Smith, J. C., Unit Operations of Chemical Engineering, McGraw-Hill Book Company, Inc., New York, p. 422 (1956). 79. McEwen, C. R., Trans. AIME. (Petr.), 216, 412 (1959). 80. Melrose, J. C., Am. Inst. Chem. Engrs., Preprint No. 63, Kansas City, May 17-20, 1959. 81. Michaels, A. S. and Timmins, R. S., Trans. AIME.(Petr.), 219, 150, (1960). 82. Milne-Thomson, L. M., Theoretical Hydrodynamics, 4th Ed., Macmillan Company, New York, p. 394 (1960). 83. Morgan, F., Trans. AIME., 189, 183 (1950). 84. Morse, R. A. and Yuster, S. T., Oil Weekly, 125, 36 (April 7, 1947). 85. Muskat, M., Physical Principles of Oil Production, McGraw-Hill Book Company, Inc., New York (1949). 86. Nelson, P., Chem. Eng. Prog., 53, 320 (1957). 87. Nenninger, E., A.I.Ch.E.J., 4, 305 (1958). 88. Newcombe, J., McGhee, J. and Rzasa, M. J., Trans. AIME. (Petr.), 204, 227 (1955). 89. Perkins, F. M., Trans. AIME. (Petr.), 210, 409 (1957). 90. Perkins, F. M. and Collins, R. E., Trans. AIME. (Petr.), 219, 383 (1960). 91. Philip, J. R., Austr. J. Phys., 10, 31 (1957). 92. Philip, J. R., Trans. Faraday Soc., 51, 885 (1955). 93. Philip, J. R., Soil Sci., 83, 345 (1957).

-23094. Pierce, B. 0. and Foster, R. M., A Short Table of Integrals, Ginn and Company, Boston, 4th Ed. (1956). 95. Pirson, S. J., Oil Reservoir Engineering, McGraw-Hill Book Company, Inc., New York (1958). 96. Rapoport, L. A. and Leas, W. J., Trans. AIME. (Petr.), 192, 83 (1951). 97. Rapoport, L. A., Trans. AIME. (Petr.), 198, 139 (1953). 98. Rapoport, L. A., Trans. AIME. (Petr.), 204, 143 (1955). 99. Rapoport, L. A., Carpenter, C. W., Jr. and Leas, W. J., Trans. AIME. (Petr.), 213, 113 (1958). 100. Ribe, K., Trans. AIME. (Petr.), 219, 1 (1959). 101. Richardson, J. G. and Perkins, F. M., Trans. AIME. (Petr.), 210, 114 (1957). 102. Richardson, J. G., Trans. AIME. (Petr.), 195, 187 (1952). 103. Rouse, H., Advanced Mechanics of Fluids, Wiley and Sons, Inc., New York (1959). 104. Saffman, P. G., Quart. J. Math. and Appl. Mech., 12, 146 (1959). 105. Saffman, P. G. and Taylor, G. I., Proc. Roy. Soc. A., 243, 312 (1958). 106. Sandberg, C. R., Gournay, L. S. and Sippel, R. F., Trans. AIME. (Petr.), 213, 36 (1958). 107. Scheidegger, A. E., Physics of Flow Through Porous Media, Macmillan Company (1957). 108. Scheidegger, A. E., Physics of Fluids, 3, 94 (1960). 109. Schlichting, H., Boundary Layer Theory, 4th Ed., pp. 238-265, McGrawHill Book Company, New York (1960). 110. Searcy, D. F. and Stewart, J. R., Paper presented at Univ. of Michigan, Research Conference on Flow of Natural Gas, Ann Arbor, Michigan, June 30, 1955. 111. Sheldon, J. W., Zondek, B. and Cardwell, W. T., Trans. AIME. (Petr.), 216, 290 (1959). 112. Silverblatt, C., Ind. Eng. Chem., 46, 1201 (1954). 113. Smith, W. D., Physics, 4, 184 and 425 (1933).

-231114. Srivastava, S. N., Am. J. Phys., 28, 299 (1960). 115. Stahl, R., Trans. AIME. (Petr.), 151, 138 (1943). 116. Steber, J. D. and Marsden, S. S., Paper presented at 19th Technical Conference on Petroleum, Penn State Univ. (1952). 117. Stewart, C. R., Craig, F. F. and Morse, R. A., Trans. AIME. (Petr.) 198, 93 (1953). 118. Streeter, V. L., Fluid Dynamics, McGraw-Hill Book Company, Inc., New York (1948). 119. Taylor, G. I., Quart. J. Math. and Appl. Mech., 12, 265 (1959). 120. Taylor, G. I., Proc. Roy. Soc. A., 201, 192 (1950). 121. Terwilliger, P. L., Wilsey, L. E., Hall, H. N., Budges, P. M. and Morse, R. A., Trans. AIME. (Petr.), 192, 285 (1951). 122. Thomeer, J., Trans. AIME. (Petr.), 219, 73 (1960). 123. Thomson, W., Phil. Mag. (4), 42, 448 (1871). 124. Van der Poel, C., Trans. AIME. (Petr.), 213, 103 (1958). o 125. Van Meurs, P., Trans. AIME. (Petr.), 210, 295 (1957). 126. Vetrov, B. N., J. Tech. Phys. USSR, 29, 1021 (1959). 127. Wagner, O. R. and Leach, R. 0., Trans. AIME. (Petr.), 216, 65 (1959). 128. Walas, S. M., Trans. Am. Inst. Chem. Engrs., 42, 783 (1946). 129. Warren, J. E. and Calhoun, J. C., Trans. AIME. (Petr.), 204, 22 (1955). 130. Watkins, J. W. and Mardock, E. S., Trans. AIME. (Petr.), 201, 209 (1954). 131. Welge, H., Std. Oil Co. of New Jersey Tech. Publ., p. 17 (1948). 132. Welge, H., Trans. AIME. (Petr.), 192, 91 (1952). 133. West, W., Trans. AIME. (Petr.), 198, 125 (1954). 134. Whalen, J., Trans. AIME. (Petr.), 198, 111 (1954). 135. Willman, B. T., Valleroy, V. V., Runberg, G. W., Cornelius, A. J. and Powers, L. W., J. Pet. Tech., p. 681, July, 1961.

-232136. Wilson, D., Fluid Distributions in Porous Media, Stanolind Oil and Gas Co., Technical Publication. 137. Wykoff, R. and Botset, A., Physics, 7, 325 (1936). 138. Wyllie, M. R. J. and Rose, W. D., Nature, 165, 972 (1950). 139. Yuster, S., Penn State Univ. Mineral Ind. Expt. Sta. Bull., No. 52 (1948).

APPENDICES -233

APPENDIX A ANALYTICAL DETERMINATION OF THE DERIVATIVE a| As a check on the consistency of the solution, the total differential of the saturation is written; as as dz dS = dt + - - at az dt This relation must hold at all points, and, in particular at the leading edge (z = zo) of the saturation front, i.e., dS as I + dz dtlz = ZO 6t z = ZO z dt dt 5zdt but dS =0 dt z = ZO from boundary condition number 3, Equation (15) ff. dz = dzo; Also s = -Al dt dt dz Z = ZO Zo Then as + Al d at any zo, say zo (31) [t dt The relationship between the dependent variable S and the independent variables z and t is shown in Figure 5. The function S = f (z, t) represents a surface in space. The definition of S is, at dS lim S_ lim S(zo, to + At) - S(Zo,t) (32) Atto - At - A

-235Let zo = Zo at to - dzo zo = zo + - At at to + At dt The derivative in Equation (32) is now evaluated by means of the expression for S obtained in Equation (17). Then, dzo A1(ZO + da At) S(zo, to + At) = 1 - [ Z dzo Zo +. At Using formula number 823 of Pierce (94), gives dzo Ai(Zo + d At) S(to + At) = 1 - [ o (1 - Z At +...)] Using formula number 828 of Pierce (94), gives,a ~ dzo r /~ 1_ d zo S(to + At) = -A1(Zo + dt tt)[log( ) + log(l - Z -t At)+....]+.... Using Pierce number 842 gives (t At) = -A1(Z+ dzo At) [log(Z ) 1 dzo At....... S(dt Z+ zt dt where all quadratic terms and higher have been omitted. Then, from (32) aS lir +A1 d-z At 2 d- At2 + A1 dzo -Al lmdz d2 2. At- o At(33) at A -o At zo and aS Al dzo 6- - d-t zoi as required by Equation (31).

-2361 t S I dt Figure 5. Surface in Space Generated When Saturation S' is Plotted as a Function of Distance and Time. S' = f(z,t).

APPENDIX B RESIDUAL EQUILIBRIUM SATURATION - TABULAR DATA FROM LITERATURE TABLE V DATA TAKEN FROM LITERATURE FOR RESIDUAL EQUILIBRIUM SATURATION-WETTING PHASE DISPLACING NON-WETTING PHASE Displacing Displaced Porous Medium Porosity % RES. obs. % Source Phase Phase Water Oil Torpedo Sandstone 24.9 28.1 (69) Water Oil Torpedo Sandstone 22.8 28.9 (69) Water Oil Sintered Alundum 26.6 25.3 (69) Water Oil Redwater Limestone 19.7 32.4 (69) Water Oil Berea Sandstone 19.3 36.7 (69) Water Oil Unconsolidated Sand 34 16.4 (127) Water Oil Unconsolidated Sand 34 18.1 (127) Water Oil Sintered Alundum 24 32 (70) Water Natural Gas Unconsolidated Sand 36 16 (42) Water Natural Gas Slightly Consolidated Sand 33 21 (42) Water Natural Gas Sintered Alundum 23.2 25.0 (42) Water Natural Gas Sintered Alundum 24.1 24.2 (42) Water Natural Gas Sintered Alundum 22.3 25.0 (42) Water Natural Gas Wilcox Sandstone 21.3 37.0 (42) Water Natural Gas Wilcox Sandstone 26.8 37.0 (42) Water Natural Gas Frio Sandstone 18.3 32.5 (42) Water Natural Gas Frio Sandstone 27.6 30.0 (42) Water Natural Gas Nellie Bly Sandstone 26.5 33.0 (42) Water Natural Gas Nellie Bly Sandstone 26.5 32.3 (42) Water Natural Gas Springer Sandstone 22.2 33.0 (42) Water Natural Gas Springer Sandstone 21.9 33.0 (42) Water Natural Gas Tensleep Sandstone 15.9 46.8 (42) Water Natural Gas Tensleep Sandstone 10.9 49.0 (42) Water Natural Gas Tensleep Sandstone 19.5 45.8 (42) Water Natural Gas Tensleep Sandstone 18.4 38.6 (42) Water Natural Gas Bradford Sandstone 13.5 35.0 (42) Water Natural Gas West Beaumont Sandstone 32.9 16.7 (42) Water Natural Gas West Beaumont Sandstone 31.0 18.5 (42) 29.6 23.2 (42) 25.8 29.0 (42) 22.9 32.0 (42) 25.7 30.6 (42) Water Natural Gas Gulf Coast Sandstone from 28.2 25.5 (42) Gas Producing Reservoir 27.9 25.0 (42) 22.7 34.0 (42) 22.5 31.0 (42) 32.3 19.1 (42) Water Oil Nellie Bly Sandstone 27.4 28.8 (52) Water Oil Rushford Sandstone 16.9 38 (8) Water Oil Glass Beads 35 18 (101) Water Oil Glass Beads 35 20 (101) Water Oil Glass Beads 38 12 (125) Water Oil Glass Powder 34 18 (89) Water Oil Unconsolidated Sand 32.9 19 (4) Water Oil Weiler Sandstone 20.0 38 (56) Water Oil Bartlesville Sandstone 23.9 33.4 (135) Water Oil Torpedo Sandstone 29.6 33.2 (135) Water Oil Weiler Sandstone 21.0 41.4 (135) Water Oil Torpedo Sandstone 26.2 33.4 (135) Water Oil Torpedo Sandstone 27.0 32.7 (135) Water Oil Torpedo Sandstone 26.0 24.3 (135) Water Oil Torpedo Sandstone 26.5 28.0 (135) Water Oil Cores from California Fields 11 48.0 (35) Water Oil Cores from California Fields 13 42.3 (35) Water Oil Cores from California Fields 11 43.2 (35) Water Oil Cores from California Fields 18 32.4 (35) Water Oil Cores from Franklin Mine Project 16.5 40 (84) Water Oil Cores from Franklin Mine Project 11.1 51 (84) Water Oil Cores from Franklin Mine Project 11.0 46 (84) -257

APPENDIX C MOUNTING OF CONSOLIDATED CORES The general technique for mounting of the consolidated cores was as follows: (1) Casting of end chambers and attaching to core. (2) Sealing off the remaining exposed surface of the core with nearly polymerized resin. (3) Perform the final casting operation. The end chambers were cast in paper cups in the form of thick, right circular cylinders. Conical inlet chambers were then machined into them on a lathe. The blanks for the end chambers could be held very well in the jaws of the lathe chuck. The polymer used was Scotchcast number 2 resin (Minnesota Mining and Manufacturing Company). The hardened polymer was found to be easily machinable without flaking or cracking. The end caps were machined down to a diameter slightly greater than the core to be mounted. The end cap is cemented to the core by placing the core over the inverted cap as shown in Figure 40a and placing nearly solidified polymer in the circular fillet, as shown. This operation is performed on both ends of the core. The next step is to seal off the cylindrical sides of the core. This operation is to prevent the capillary imbibition of the low viscosity polymer which is poured around the core in the final casting operation. The materials needed for the sealing operation are simply a spatula and some very tacky, almost hardened polymer. The polymer hardens on the surface of the core almost as soon as it is applied. After the end caps have been attached and the core surface sealed, the final casting is accomplished as illustrated in Figure 4Ob. -258

-239A piece of flat glass coated with silicone grease is used as a base. A glass tube, of diameter slightly larger than the end caps is then pressed firmly into the grease. The core (with end caps) is then centered in the glass tube. The resin is then poured into the annular space between the core and the glass tube to a depth of about 1 inch, and allowed to harden and seal the base. The remainder of the polymer is then poured to the indicated level and allowed to set. Fittings are then screwed into the end caps after drilling and tapping through the center of the end cap to the inlet chamber. The volume of the inlet chambers was about 5 cc.

-240CORE CLAMP FOR HOLDING CORE HAS BEEN OMITTED ANNULUS TO BE FILLED WHEN CEMENTING CAP END CAP~ //0 Figure 40a. Schematic Diagram of Mounted Consolidated Core and Preparation Procedures —Section View Showing Arrangement for Attaching End Caps to Core Before Cementing.

-241FILL TO HERE- GLASS TUBE, CORE WITH END CAPS ATTACHED AND PRE-COAT SPACE FOR RESIN - NOTE: CLAMP FOR HOLDING GLASS TUBE HAS BEEN OMITTED. GLASS PLATE Figure 40b. Schematic Diagram of Mounted Consolidated Core and Preparation Procedures —Sectional View Showing the Setup for Making the Final Casting Before Filling.

APPENDIX D MOUNTING OF UNCONSOLIDATED MATERIALS The unconsolidated materials were mounted as follows. The sand was poured into the glass containing tube with continual tapping to insure the closest possible pack. As the sand level neared the top of the tube a rubber stopper was inserted, and pressed firmly down until its base touched the top of the pack. The rubber stopper had a 3/32" hole bored through its center. The tapping of the tube continued and more sand was introduced through the hole until the sand column was flush against the bottom of the rubber stopper. This process continued until no further settling of the sand bed occurred. The outlet hole from the rubber stopper was then sealed with electrician's tape, and the rubber stopper sealed firmly in place by casting an end block of Scotchcast number 2 resin around it in a paper cup. This is illustrated in Figure 41. The tape was then removed from the outlet port. The mounting was then completed by inserting a 1/8" steel tube firmly packed with glass wool into the cylindrical hole in the rubber stopper. -242

-243STEEL TUBE GLASS WOOL / -RUBBER STOPPER HARDENED RESIN GLASS TUBE SAND Figure 41. Schematic Diagram of Mounted Unconsolidated Core.

APPENDIX E RESIDUAL EQUILIBRIUM SATURATION - TABULAR DATA FOR FIGURES 43 AND 44 Initial Water Water Displacement Residual Equilibrium Core* Saturation Rate ( ft3 Air Saturation % Rft2yr. A 24.6 257 19.1 A 26.2 514 18.7 A 27.8 1028 18.7 A 25.3 2056 19.2 A 25.1 4112 19.2 A 26.5 8224 19.6 A 24.9 16448 20.0 A 27.1 32896 21.2 B 23.9 146 19.1 B 24.2 292 18.7 B 24.1 584 19.4 B 25.8 1168 19.5 B 24.7 2336 19.5 B 26.2 4672 19.4 B 24.0 9344 17.8 B 26.9 18688 17.6 C 24.4 146 18.3 C 25.6 292 18.0 C 23.8 584 17.6 C 24.7 1168 17.8 C 25.4 2336 18.6 C 24.1 4672 19.1 C 24.0 9344 17.5 C 26.2 18688 18.5 D 36.0 105 27.3 D 33.1 210 27.1 D 34.2 420 26.2 D 34.0 840 27.0 D 35.4 1680 28.2 D 35.8 3360 27.5 D 32.4 6720 27.0 D 35.3 13440 25.9 E 11.9 386 12.8 E 11.5 772 12.8 E 10.7 1544 12.8 ~2U4

APPENDIX E (CONT'D) Initial Water Water Displacement Residual Equilibrium Core* Saturation % Rate ( ) Air Saturation % ft yr. E 11.0 3088 12.8 E 11.3 6176 11.8 E 12.4 12352 12.5 E 10.4 24704 13.2 E 11.1 49408 13.0 D 11.5 386 12.0 D 12.2 772 12.0 D 10.4 1544 11.0 D 11.9 3088 11.4 D 11.4 6176 11.3 D 10.8 12352 10.9 D 10.1 24704 11.2 D 11.7 49408 12.0 * Key to core designation given in Table VI

UNIVERSITY OF MICHIGAN 3 9015 03127 3520