Two-Phase Flows in Porous Media

TWO-PHASE FLOWS IN POROUS MEDIA H. J. MOREL-SEYTOUX1 Department of Civil Engineering Colorado State University Fort Collins, Colorado I. INTRODUCTION

120

A. Raisond'fitre B. Subsurface Hydrology Problems C. Petroleum Production Problems II.

BASIC DEFINITIONS, CONCEPTS, AND EQUIVALENCES

A. B. C. D. E. F. III.

Petroleum Reservoir Terminology Petroleum Fluid Terminology Darcy's Law for Saturated Flow Darcy's Law for Unsaturated Flow Darcy's Law for Two-Phase Flow Capillary Pressure

UNSATURATED FLOW THEORY IN A LINEAR MEDIUM

A. Mathematical Formulation of Conservation Laws B. Imbibition into a Horizontal Slab C. Infiltration into a Vertical Column IV.

TWO-LIQUID FLOW THEORY

A. B. C. D. E. F. G.

120 121 121 123

124 125 126 128 128 130 131

131 132 134 139

Mathematical Formulation of Conservation Laws 140 The Integral Equation 142 Imbibition into a Thin Horizontal Slab 143 Displacement Process in a Horizontal Core 146 Displacement Process in a Horizontal Slab with a Closed End . . 149 Infiltration into a Vertical Column 150 Drainage in a Vertical Column 155

V. LIQUID-GAS FLOW THEORY

157

A. B. C. D. E.

Mathematical Formulation of Conservation Laws 157 The Integral Equation 158 Imbibition into a Horizontal Slab 158 Infiltration into a Vertical Column with a Closed Bottom . . . . 1 5 9 Experimental Observations of Infiltration into Vertical Columns with Closed Bottoms 174 F. Hysteresis Effect During Infiltration and Redistribution . . . .176 G. Experimental Study of the Effect of Air Movement and Compressibility in a Heterogeneous Medium 176

1

Present address: 3 rue de Paris, 38320 Eybens, France. 119

H. J. MOREL-SEYTOUX VI.

QUASI-LINEAR THEORY FOR TWO-LIQUID FLOW

A. B. C. D. VII.

Reservoir Behavior Simulation Saturation Equation for a Quasi-linear System Discretized Quasi-linear Model Typical Reservoir Conditions and Operations

APPROXIMATE TREATMENT OF THREE-PHASE FLOW

A. Water-Oil-Gas Systems B. Steamflood Displacement Process

179

179 180 182 184 193

193 195

VIII. TWO-DIMENSIONAL TWO-PHASE FLOWS

199

REFERENCES

200

I. Introduction A. RAISON D'ETRE

With the proliferation of scientific and technological journals it is becoming more and more difficult to read the complete literature in one's own field. Treatises that synthesize the current state of the art are welcome additions since they extract the essential information from hundreds of articles, permit the reader to bypass the task of actually reading most of these articles, and direct him to the really fundamental papers in his field of interest. Two more such treatises in the field of flow through porous media (Bear, 1972) and soil and water (Hillel, 1971) have recently appeared. When the books cover a broad field they run two risks: one of being at times rather shallow and the other of being partly obsolete by the time of publication. It is hoped that the present article will avoid these two hurdles because it is both narrow in scope and limited in length. This article was conceived with three purposes in mind. The first purpose was to present a review of what had already been done in the field of twophase flow of interest to the soil scientist and the hydrologist. The second purpose was to review the "petroleum industry know-how" in the area of two-phase flow in porous media in a form which would appeal to the hydrol­ ogist or soil scientist. Finally, the third purpose was to show how these particular techniques and their extensions can be used fruitfully for a better understanding and description of the movement of water and air in soils. In summary, the basic philosophy of this article gravitates around the following central themes: (1) water movement in the soil is (at least) a twophase phenomenon and may not be at times adequately described or under­ stood by the prevailing current unsaturated (one phase) approach, and (2) much of the petroleum literature know-how can be used effectively to solve two-phase flow problems. 120

Two-Phase Flows in Porous Media B. SUBSURFACE HYDROLOGY PROBLEMS

Many current refined models of the behavior of a watershed, the "rainfallrunoff" models, are not properly named. More precisely they really are "excess-rainfall-direct-runoff" models in which the infiltration and the interflow components are treated somewhat casually. It is apparent from the literature that the development of subsurface hydrologic models has been lagging behind the development of surface models. Past concern for unsaturated soil water movement seems to have been confined to the ranks of the soil scientists. Indeed the probably most often quoted work on the subject of infiltration (Philip, 1957), certainly a topic of great importance in hydrology, was published in the Soil Science Journal. According to the same author (Philip, 1969) however "their colleagues in engineering hydrology have exhibited an increasing interest in this field in recent years." Prior to 1969 with the exception of a relatively few investigators (Free and Palmer, 1940; Horton, 1940; Corey, 1957; EIrick, 1961; Wilson and Luthin, 1963; Youngs and Peck, 1964; Adrian, 1964; Peck, 1965b; Adrian and Franzini, 1966; McWhorter and Corey, 1967; Garner et αί, 1969) the presence of the second phase, air, was largely ignored. Theoretical efforts toward the incorporation of the second phase in the description of soil or ground water movement started (in the author's opinion) in 1969 (Phuc, 1969; Brustkern, 1970; Noblanc, 1970; Green et al, 1970; Brutsaert, 1970; McWhorter, 1971). Quite recently other publications (Vachaud et al., 1971, 1972; Dixon and Peterson, 1971) substantiating the possible serious consequences of neglecting air in subsurface hydrologic models have appeared. C. PETROLEUM PRODUCTION PROBLEMS

The petroleum engineer is concerned with the production of oil. Oil seldom occurs alone. Even when only oil is initially present in the reservoir (that is the ground formation containing oil) and not in contact with gas or water, it has the ability to generate a second phase: gas. Gas is liberated from solution in the oil as the pressure in the oil decreases. The development of a second phase may occur either throughout the entire reservoir or in the neigh­ borhood of the well where oil flows into tubing and through the tubing to the surface. The understanding of simultaneous flow of immiscible fluids in a porous medium is thus a prerequisite to sound operations in the practical business of production of oil and of natural gas. The petroleum industry has devoted much effort to gaining a better understanding of this subject. This knowledge can be ben ficially introduced to help resolve subsurface hydrology problems. 121

H. J. MOREL-SEYTOUX

30i

,(l) ^^^(100)

Repeated 5-Spot Single Interface 25

o

/ /(50)

20

o 1

15

l·

^ ^ ^

(0.5)

/

(2)

-^-^

I 0 hh

(,0)

\-

5l·

1 \l 0.5

f/^^

i

1

I

1.0

1.5

2.0

(Normalized) Pore Volumes

|

Injected

FIG. 1. Influence of mobility ratio on producing water-oil ratio in a repeated five-spot (Morel-Seytoux, 1966).

Cumulative

Injection, VhP

FIG. 2. Composition of produced fluid versus volume injected, five-spot, mobility ratio of 10 (Fitch and Griffin, 1964).

122

Two-Phase Flows in Porous Media CO _l GO rlOjOOOi CD T> - eooo

5 r 6000 t_ Q) - 4000

£

8. -

1>»

2000

- 1000

L *<5 υ

fJ

o

j-

800

-

600

(I) ._

-

400

m

-

200

*H± L

o

*fi

*£

σ»0

i ,


FIG. 3. Performance history, Pilot II (Goolsby and Anderson, 1964).

It is true that even petroleum engineers have attempted to circumvent the difficulties associated with two-phase flow problems by assuming sharp inter­ faces between fluids to render problems amenable to mathematical treatment. The simple model thus obtained has proved valuable in a variety of cases (e.g., reservoirs in which oil and gas are thoroughly segregated). However, in many instances the assumption of a sharp interface is too gross to yield exploitable results. Generally the problem is not to calculate where one fluid ends and where the other one begins, but rather to calculate how much there is of each fluid throughout the porous medium. In fact one of the preoccupa­ tions of the reservoir engineer is to forecast for instance the ratio of water to oil (WOR, water-oil ratio) production at the well head. The petroleum litera­ ture is replete with theoretical, experimental, or field curves such as shown in Figs. 1, 2, and 3. For a more complete discussion of the technology of petro­ leum production, the reader is referred to standard texts (e.g., Smith, 1966). The concept of saturation is thus basic to reservoir engineering. Saturation of a given phase is defined as the volumetric fraction of the total pore space occupied by the phase.

II. Basic Definitions, Concepts, and Equivalences The common denominator of the two disciplines (reservoir engineering on one hand, subsurface hydrology on the other) is evident. In both areas two immiscible fluids are present: water and oil (or oil and gas) and water and air. 123

H. J. MOREL-SEYTOUX

However there are distinct differences which probably explain why there has been little interplay between the disciplines. In hydrology the mechanisms for water and air movement in the soil are primarily natural, so that fluid flow movement is caused by forces of relatively small magnitude. In petroleum production on the contrary fluid flow movement is often induced by con­ siderable external pressure gradients. For example, under the common practice of water injection to displace oil (a procedure called water-flooding in the petroleum literature) pressure differences of several thousand pounds per square inch between an injector (well) and a producer (well) are not unusual. It is assumed that the reader is familiar with the classical hydrologic con­ cepts. For this reason, only a few terms frequently used in reservoir engineering will be defined (Calhoun, 1053). A. PETROLEUM RESERVOIR TERMINOLOGY

The oil and gas bearing zones (see Fig. 4) form what is known as the reservoir. The fraction of the reservoir which is filled with gas is known as the gas cap. The gas cap may be present at the time of discovery of the oil field

FIG. 4.

Schematic cross section of an oil reservoir. 124

Two-Phase Flows in Porous Media or may develop later as gas liberated from solution due to a drop in oil pressure proceeds updip by buoyancy. The oil zone and the gas cap are usually surrounded by aquifers or aquicludes. B. PETROLEUM FLUID TERMINOLOGY

Oil, gas, and water are compressible fluids. A given mass of oil at a given pressure and temperature, e.g., the pressure and temperature which prevail in the reservoir, occupies a different volume in the reservoir than when brought to the surface and stored under atmospheric conditions. The oil production of a field is measured at atmospheric pressure and measured in stock tank barrels (STB). The dimensionless ratio of the volume of oil occupied in the reservoir at a pressure p to the volume of oil N it occupies at atmospheric conditions (standard conditions) is known as the oil formation volume factor: Reservoir volume at pressure p B°(p) (1) N It is expressed in reservoir barrels/stock tank barrels (RB/STB). A typical curve of oil formation volume factor versus pressure is shown in Fig. 5. Gas CD H

l.50i

ω 1.40

1.30

1.00

1000 Pressure,

2000 PSIA

FIG. 5. Formation volume factor of the Big Sandy Field reservoir oil, byflashliberation at reservoir temperature of 160°F (Craft and Hawkins, 1964).

and water formation volume factors are defined similarly. The gas formation volume factor is, however, expressed in reservoir barrels per standard cubic feet (RB/SCF) or reservoir barrels per thousand cubic feet (RB/MCF). In Fig. 5 the oil formation volume factor displays a bend at a pressure labeled bubble point pressure. On the high pressure side of the bend, a volume 125

H. J. MOREL-SEYTOUX

expansion corresponds to a drop in pressure. When the pressure reaches the value corresponding to the bend, bubbles of gas evolve from the oil. The loss of solution gas from the oil more than offsets the volumetric expansion of the remaining oil phase. Further pressure drops thus bring about an oil shrinkage. A given mass of liquid hydrocarbon, which originally occupied a certain volume in the reservoir, will separate at the surface into a volume N of oil and a volume of gas evolved from solution in an amount proportional to N. The constant of proportionality is known as the solution gas oil ratio Rs, and is expressed in SCF/STB or MCF/STB. A typical solution gas-oil ratio is shown in Fig. 6.

700

1

567scf/stb.

ω 600 CO

£

500

CO

° § 3

200

co

100

At 1200 PSIA — R8 = 33 7 *

Φ

0

«A M Φ

«_ 1 » • · Si

300

0.

"5 'c

CD

1000

2000

*

3000

Pressure , PS I A

FIG. 6. Solution gas-oil ratio of the Big Sandy Field reservoir oil, by flash liberation at reservoir temperature of 160°F (Craft and Hawkins, 1964).

C. DARCY'S LAW FOR SATURATED FLOW

In hydrology the most familiar form of Darcy's law (Darcy, 1856) is written as Q = AKAh/L

(2)

where Q is the volumetric flow rate, A the bulk cross-section area of the soil column in a direction perpendicular to flow, K is the (saturated) hydraulic conductivity and Ah is the drop in head in the direction of flow over the length L. The head is defined as h = (p/pg) - z

(3)

where p is fluid pressure, p is the fluid density (or specific mass), g is the acceleration of gravity, and z is the vertical coordinate oriented positive 126

Two-Phase Flows in Porous Media downward (this convention will be used throughout the article with the exception of Section VI). Equation (2) is usually written as Q/A = v = KAh/L

(4)

and v is called the (Oarcy) fluid velocity. It is known that K is both a function of the soil characteristics (ability to transmit fluid) and of the fluid properties (density and viscosity). As long as hydrologists were concerned solely with the flow of water, it was not deemed essential to separate the various constituents making up the hydraulic con­ ductivity. Probably for this reason Darcy's law was generalized in vectorial form by the hydrologists (De Wiest, 1967) as v = - K g r a d / * = -KVh

(5)

where grad or V are symbols for the gradient (vector) operator. It was recognized later that Eq. (5) was only valid for an incompressible fluid. For a compressible fluid it has been shown fairly convincingly that Darcy's law should be generalized as

= -Kv[j(\/pg)dp-z\

(6)

In the petroleum industry (Muskat, 1949) Darcy's law was generalized early in the form y=-(k^)Vp

+ (k/ß)pgVz

(7)

where k is intrinsic permeability (a property of soil reflecting its ability to transmit fluids, no longer dependent on fluid properties) and where μ is fluid (dynamic) viscosity. Possibly just by chance, Eq. (7) is also valid for com­ pressible fluids. In fact Eq. (6) can be rewritten equivalently as v = -(K/pg)Vp

+ KWz

(8)

The vector V[fj0 (l/pg) dn], where π is a dummy variable of integration, has components d/3£[f{!0 (\/pg) dn] where ξ is any cartesian coordinate. For a compressible fluid p = p(p), i.e., density is a function of pressure and the integral is a function of p, the upper limit of integration. Applying the chain rule of differentiation one obtains d_ TrP dn 1 = d_ TrP dn 1 dp δζ [J p o p(n)g\ dp [iPo p(n)g\ δξ and applying the rule of differentiation of an integral with respect to the upper limit variable, d \ r dp~\ 1 dp (10) δξ V pg\ P(P)9 5ξ 127

H. J. MOREL-SEYTOUX

In other words,

(l/w)Vp = v[j(l/p«r)
(11)

From Eq. (11) and given the equivalence relation between (hydraulic) con­ ductivity and (intrinsic) permeability, namely K = kpg/μ

(12)

one sees that Eqs. (6) and (7) are completely equivalent. Equation (7) seems easier to write and it will be taken in this article as the general form of Darcy's law for saturated flow, valid for an incompressible or a compressible fluid. D. DARCY'S LAW FOR UNSATURATED FLOW

This particular generalization of Darcy's law was introduced by soil scientists (Richards, 1931). Because they were concerned only with the flow of water (an incompressible fluid under natural hydrologic conditions of minor pressure differences) under unsaturated conditions, Darcy's law was generalized in the form: v = -K W VA = -JK(0)VA

(13)

where X w (or Κ(θ)) is the unsaturated hydraulic conductivity which is a function of 0, or volumetric moisture content (i.e., volumetric fraction of the bulk soil volume occupied by water). E. DARCY'S LAW FOR TWO-PHASE FLOW

This generalization was introduced by the petroleum engineers (Muskat, 1949; Pirson, 1958; Smith, 1966). For phase / (water, oil, gas, or air), Darcy's law takes the form y.^-h^Vpi

+ k^PtgVz

Mi

(14)

Mi

where kTi is the relative permeability to phase / which is a function of the saturation of phase /. By definition, kTi is dimensionless and varies between 0 and 1. Saturation of phase ι is defined as the volumetric fraction of the total pore space occupied by phase /. If the phase i is water, its saturation Sw is simply related to the moisture content by the equation 0=<£S W 128

(15)

Two-Phase Flows in Porous Media where φ is porosity. Typical curves of relative permeabilities to water and oil for a water-oil system are shown in Fig. 7. Note that both relative permea­ bilities are expressed as functions of water saturation. It is convenient to use only one saturation variable since the water and oil saturations are simply related by the equation: S w + S0 = 1

(16)

Figure 8 shows relative permeabilities to water and air for a water-air system. In either case (Figs. 7 and 8) the relative permeability to water goes to zero

j\

i.O

0.8

σ α>

\

ro

\

E Φ Q.

a> tr

\

V

5 0.6

>

\\

0.4

0.2

0

0

/ K rw

0.2

0.4 0.6 Water Saturation , S w

0.8

1.0

FIG. 7. Typical relative permeabilities of Kansas-Missouri and Wyoming heavy oil sands (Holm et al, 1965).

before its saturation reaches the value zero. The (water) saturation at which the relative permeability goes to zero is termed the residual (water) saturation throughout this article. More accurately it should be called the displacement residual water saturation. It is the limiting lowest water saturation that can be achieved through attempts of displacement of water from the medium by another fluid. The water content could be reduced further by heating the soil column. In this case however the water would not leave the medium in liquid form as a continuous phase but in the form of vapor. The curve of Fig. 8 does not include this " o t h e r " relative permeability. The residual water (air) saturation is denoted S wr (5 ar ). 129

H. J. MOREL-SEYTOUX 100 90 80

<£

70

>»

5 60 oΦ b 50

*

a>

40

o> 30 a> (XL

o • Δ a

krw(Surface) krw (Subsurface) kra (Surface) kra (Subsurface)

20 10 ~0

10 20 30 4 0 50 60 70 80 90 100 Water Saturation %

FIG. 8. Relative permeabilities as a function of saturation (water-air system) (Corey, 1957). F. CAPILLARY PRESSURE

Between two immiscible fluids there exists at the interface a pressure difference. This pressure difference is called the capillary pressure; by defini­ tion, Pc = Pnw - Pw

(17)

where p nw and p w are the pressures in the nonwetting and wetting phases respectively. Petroleum engineers use this definition of capillary pressure very uniformly. A typical capillary pressure curve for a water-air system is shown in Fig. 9. By definition pc is always a positive quantity. Soil scientists being only interested in the water-air system have defined a water pressure, expressed as a height of water, relative to atmospheric pressure and thus negative and denoted often by the symbol Ψ. If pc is ex­ pressed as a water height one can then define: K = Pclp*9 Similarly, one can define an equation of equivalence as V=-hc

+ (pa-pA)/Py,g

(18) (19)

where pA is atmospheric pressure and p a is air pressure. It is assumed that the reader is familiar with the basic concepts of wettability, capillary pressure, the various modes of entrapment of a wetting or nonwetting phase, and hysteresis; in the negative the reader is referred to a text (De Wiest, 1969). 130

Two-Phase Flows in Porous Media 0.30,

r

' ' ΓΤ

I

\ \

i \

0.25

r o D Δ

^ \

1

1

Core 1 Core 2 Core 3

κ \

0.20 -

*

r

\

>\

\

V\

\

ΗΛ

0.15

v

^·

■\\

\S \

0.10

V

^

-Q^

X

\

0.05

1

1

1 0.20

1

1 1 1 1 0.40 0.60 Water Saturation

1

0.80

1

Ί

1.00

FIG. 9. Static capillary pressure as a function of saturation for sandstone cores (airwater system) (Handy, 1960).

III. Unsaturated Flow Theory in a Linear Medium A. MATHEMATICAL FORMULATION OF CONSERVATION LAWS

It is not the purpose of this section just to rederive the well-known Richards' equation and to show how to solve it by well-known and well-used methods. The equation will be rederived to make very explicit the assumptions tradi­ tionally made in the derivation of the equation and to emphasize thereby the limitations of the equation in describing certain phenomena. Finally, a new method of solving the equation will be introduced for a particular simple problem. The one-dimensional version of Eq. (13) in the vertical direction has the form i;w = X w (l - dhjdz)

(20)

where vw is the water velocity (in a Darcy sense, i.e., a volumetric flow rate per unit bulk cross-section area), Kw is the unsaturated hydraulic conductivity (a function of moisture content), hw is the water pressure head (i.e., its pressure 131

H. J. MOREL-SEYTOUX

measured as the height of a water column), and z is the vertical coordinate oriented positive downward. Logically one would write down an equation similar to Eq. (20) for air flow. Clearly if water moves in, air moves out. However, traditionally the air movement has been assumed unimportant and no Darcy's equation is written for air movement. From the principle of mass conservation a continuity equation is derived as de/dt + dvjdz = 0

(21)

where Θ is moisture content and t is time. Substitution of Eq. (20) in (21) yields

Tt^Jz [ " Χ - ( Θ ) dhJdz

+

*w(0)] = °

(22)

This equation provides one single equation for two unknowns, Θ and /zw. Again traditionally the method of circumventing this difficulty is to neglect air compressibility and viscosity by assuming that air pressure is atmospheric (or constant) and remains so at all places and at all times in the medium. In this case, then, the water head is simply the negative of the capillary pressure head hc, and therefore a function of Θ only (neglecting hysteresis). With these assumptions Eq. (22) transforms to

Traditionally a diffusivity is defined: £>(©) = - K w ( 0 ) dhjde

= K w (0) ί/Ψ/ί/θ

(24)

By its very definition, Z)(0) is always a positive quantity. Equation (23) takes the final customary form

ΙΓ-&1* Θ ) &] + Ί>Γ" 0

(25)

This equation is usually referred to as Richards'' equation or as the diffusivity equation. The problem is now to solve Eq. (25) given boundary and initial conditions. When gravity is neglected, Eq. (25) takes the form <90

d r dei

It B. IMBIBITION INTO A HORIZONTAL SLAB

In this section the problem of water penetration into a semiinfinite hori­ zontal core with a uniform initial moisture content 0 j from a plentiful water tank is considered. Let x be the horizontal coordinate. The equation describ132

Two-Phase Flows in Porous Media ing the motion is Eq. (26) with x substituted for z. The boundary condition at x = 0 for all times (t > 0) is © = Θ where Ö represents moisture content at residual air saturation. Throughout this article the tilde (~) sign over any quantity or function say ( ~ ) indicates that the quantity ( ) is evaluated at residual saturation of the other phase. For example, in a water-air system 5W means the value: (1 - Sar), Xw means XW(5W), £ra means &ra(Swr) etc. It has been recognized (Philip, 1969) that Eq. (26) can be reduced to an ordinary differential equation by a similarity transformation x = w(©)i1/2. The result­ ing ordinary differential equation, udG 2 du

d Γ de~\ du I du\

(27)

is solved numerically and its solution is denoted u0{&). The solution of Eq. (26) is then known, and is denoted x = u0(&)t1/2 (28) This form of the solution indicates that the profile stretches with time without deforming. In other words, given any two moisture contents ©j and © 2 and their locations xu and x2t at any time t and their locations at any other time τ, xu and x2x (see Fig. 10), it is always true that (29)

* l r / * l t = *2f/*2T

This result is evident from Eq. (28) by calculating the relative change of position of a given moisture content, namely /dx\ _ 1 w0(©)f 1/2 dt = I \ Τ / θ " 2 ι#0(Θ)ί 1/2

27

Moisture Content

Θ

FIG. 10. A stretching type of moisture profile evolution with time. 133

(30)

H. J. MOREL-SEYTOUX

The relative change of position of a given Θ does not depend on the value of Θ. Another way yet to look at this result is to consider the rate of propagation of a given moisture content Θ. From Eq. (28) one obtains an expression for this rate, namely (dx/dt)e = u0(Q)l2t1'2 (31) This result shows that the rate of movement of a given Θ does depend on the value of Θ and slows with time. It is interesting to observe that the stretch­ ing nature of the solution can be deduced from the fact that Eq. (26) trans­ forms to (27) by a similarity transformation, without solving (27). Similarly, from Eq. (31) one can deduce that the imbibition rate will vary as the inverse of the square root of time. C. INFILTRATION INTO A VERTICAL COLUMN

1. Philip's Perturbation Solution The equation describing the soil water movement is Eq. (25). Philip's method of solution (Philip, 1957; Swartzendruber, 1969; Childs, 1969; Philip, 1969; Hillel, 1971) is to view Eq. (25) as a combination of Eq. (26) plus a perturbation term, dKJdz. One then expects that the solution of Eq. (25) can be fruitfully put in the form z = z0 + Δζ. Naturally, no gain is achieved if Δζ is not a separable function of Θ and t. Thus the procedure is to look for a solution of Eq. (25) in the form z = u0(&)tx'2 + ul(e)g(t) + ε(Θ, t). After substitution in Eq. (25) one obtains an equation for uu one neglects the terms in ε, and one selects the function g(t) that reduces the approximate remaining equation to an ordinary differential equation. Repeating the operation a number of times, one obtains a solution in the form z(0, 0 = u0(G)t1/2 + «!(0)i + w 2 (0)i 3/2 + · · ·.

(32)

From this solution one can deduce an equation for the rate of infiltration as a function of time: / = crl/2

+ A + Dt1/2 + Et + Ft2,12 + · · ·.

(33)

where the constants C, A, . . . , depend on the soil characteristics and the boundary conditions 0 j and 0 b where 0 , is the initial uniform moisture content and 0 b is the prescribed moisture content at the soil surface. Equa­ tion (33) predicts that at small times the infiltration rate is high, which is appropriate. For large times, Eq. (33) predicts that the infiltration rate increases indefinitely. Since this result is totally contradictory to experimental evidence, this indicates that the method of solution does not converge for large times, which is clear from Eq. (32). It is accepted practice (Philip, 1969) 134

Two-Phase Flows in Porous Media now to truncate Eq. (33) after the constant term A and to interpret A to be the hydraulic conductivity at residual air saturation K. The inconsistency of the results arises from the implicit assumption of the perturbation procedure that gravity is small compared to capillarity. At large times (of the order of only 1 hr in some instances), the gravity force will actually become the dominant one. Finally it must be noted that the solution given by Eq. (32) is valid (mathematically speaking if not physically) also for a closed medium of finite depth, at least until the wetting front reaches that lower boundary. A similar remark holds for Eq. (33). 2. Moving Strained Coordinate Solution The boundary conditions are again the traditional (simple) ones of constant moisture content at infinity downward from the soil surface, denoted 0 i ? and of a given moisture content at the soil surface, denoted 0 b and not necessarily constant with time. Initially the moisture content profile was uniform at the value Θ{. a. PHYSICAL BASIS OF THE METHOD OF SOLUTION. The boundary value problem associated with Eq. (25) is very difficult due to the nonlinearity of Eq. (25). It is not possible to solve it exactly except under very special circum­ stances. Thus, in general, one seeks an approximate solution. A " g o o d " approximate solution is one that (1) yields a solution close (in some sense) to the exact solution and (2) achieves criterion (1) at minimal computational expense. There are many ways to derive approximate solutions. One way is to consult standard texts of mathematical analysis (Ames, 1964), look at the repertoire of known techniques, e.g., series expansion, successive approxima­ tions, Ritz' methods, etc., and try them all. Another approach is to gain a knowledge of the qualitative behavior of the solution from scrutiny of the experimental information and from a physical understanding of the flow phenomenon. The second approach was followed. First it can be stated (and this without need of experimental evidence) that the evolution of the moisture content profile is the combined result of a translation and of a deformation. In some cases translation dominates, in others deformation dominates. For example, when gravity is not operating as in the case of, e.g., "free water imbibition " (Morel-Seytoux, 1969) in a semiinfinite horizontal soil column, it is well known (Philip, 1969) that the behavior of the moisture content profile is of the stretching kind, a special combination of translation and deformation. This result is due to the fact that the ratio of viscous resistance to the capillary drive, though varying with moisture content, is the same at all times. When a variety of forces interact, such as the capillary drive, the force of gravity, an external pressure gradient, the viscous resistance, etc., the evolution of the 135

H. J. MOREL-SEYTOUX

profile is not likely to be of the stretching type because the relative importance of these various forces varies in time and with moisture content. However, even in this complex case, there may be a region where translation dominates and a region where deformation dominates. Close examination of experimental results (Morel-Seytoux, 1969; Swartzendruber, 1969) shows that during infiltration there is a zone, in the lower range of moisture content, where the profile propagates ultimately essentially as a translation. Experimental profiles at various times can be superposed and they match very well. On the other hand, in the region behind the front and up to the soil surface, the profile seems to evolve essentially by stretching. In this region particularly close to the soil surface, the curvature of the profile is not pronounced. From this experimental evidence, one can formally design a mathematical approximation technique to solve the infiltration problem. b. FORMAL MATHEMATICAL APPROACH. Let us define a new coordinate system Z, traveling, so to speak, with the wetting front, namely Z = z - C(0

(34)

where ζ(ί) is the coordinate of a point (yet undefined) traveling with the front. Utilizing (34), Eq. (25) transforms (Morel-Seytoux, 1972) to δθ

άζθθ

δ Γ ,^ΘΊ ÖK 0) + =

^-^äz-äzh

Λ

äzj äz °

(35)

which is the unsaturated flow equation in the new coordinate system. This equation will be simplified differently depending upon the region of interest. c. MOISTURE EQUATION IN THE FRONT REGION. Since in the front region the moisture profile propagates essentially by a translation, an observer traveling with the front does not witness much change with time of the profile. In other words δθ/dt (short-hand for 8Θ/δήΖ not to be confused with d®/dt\z) is negligible in Eq. (35). The approximate moisture equation in the front region is, therefore, D{&)

Tz[

Tz-

K{&) +

dCdQ

Έ-ζ'°

<36)

At a fixed time this equation is an ordinary differential equation in Z, and can be rewritten as

H I - * W + ;JH = 0

3Z _

136

(37)

Two-Phase Flows in Porous Media which can be integrated with respect to Z, with the result D(6) j | _ Ζ)(θ,)(^) - [Κ(θ) - K(0,)] + j t (θ - β,) = 0

(38)

At the downstream end the profile is uniform ahead of the wetting front and the second term of Eq. (38) vanishes. Equation (38) can be solved for δθ/δΖ = δθ/dz with the result, back in the original system of coordinates, as άθ

Κ(θ) - Κ(β,) - (
Finally, the equation of the moisture profile in the front region is ζ(θ, 0 = f

D(6) άθ + 2(0 J e Κ(θ) - Κ(θ{) - (άζ/Λ)(θ - 0i)

(40)

where 0m is the upper moisture content in the front region and where Θ <6m. Note that Eq. (40) involves two unknown functions of time: 0m and άζ/dt. Note also that the shape of the profile is still changing with time in the front region because both 0m and dC/dt are functions of time. d. MOISTURE EQUATION IN THE SURFACE REGION. In the surface region the moisture profile displays little curvature, and Eq. (25) can be approximated as dz/dt = K'(0) - D'(0) δθ/δζ

(41)

where the prime indicates differentiation with respect to Θ. Equation (41) can be rewritten (again neglecting second-order terms) as

d(om)idt = K\eb) - zr(0 b )(i/o

(42)

where am = dz/δθ defines the constant slope of the profile in the surface region. Equation (42) is solved by a finite-difference approximation, either explicitly or implicitly. With the explicit scheme, Eq. (42) gives om = σ ^ + [K°(ßh) - Z T ( 0 b ) / O Δί

(43)

where am° represents the slope of the profile at time t — At. e. MATCHING CONDITION. The solutions for moisture content in the two regions must match at the interface of the two regions where moisture content has the value 0 m . At that point the slopes of the two profiles must be equal. Expression of this condition yields

σ

=

ma 137

(44)

H. J. MOREL-SEYTOUX

This equation determines αζ/dt in terms of 0 m , the only remaining unknown Elimination of αζ/dt between Eq. (44) and Eq. (39) yields for the front region: r

0m

Ώ(θ) αθ

where Km is short-hand for K(6m) and similarly for Kx and Dm. f. MATERIAL BALANCE CONDITION. It remains to express the condition that the net water influx must equal the change of moisture content in the medium. The change of water content in the medium between the time t — At and t is

" ™-Κ·-(«"-*<-%)& (46) where W0 is the water content of the soil at time t — At. The net water influx is obtained from Darcy's law, namely AW = (Kh - DJam - K{) At

(47)

Equality of the right-hand sides of Eqs. (47) and (46) provides one equation in the only unknown 0 m , namely, ,

-

°-f V(0b2 2

°

-ei)+ f m

m

Ρ{θ)-άθ

...~. -

J», Κ(θ) -Κ;

r.....

..

WJ ~ K;

zΚΘ„)](Θ-ΘΛ x e j i / Θ - Θ,

-W0 = {Kb-y-

Κ^ At

(48)

Once the proper value of 0m is obtained, the infiltration rate for the period t - At to t is simply AW/At. g. REMARKS CONCERNING THE BOUNDARY CONDITIONS AT THE SOIL SURFACE.

If the boundary condition at the soil surface is one of flux, i.e., a fixed influx rate is prescribed with value less than the maximum potential infiltration rate, then the explicit scheme of Eq. (43) is adequate. The value of 0b is deduced from the condition that I=K(eh)-D(eh)/am 138

(49)

Two-Phase Flows in Porous Media However, because am depends on 0 b , Eqs. (49) and (43) must be solved simultaneously. If, on the other hand, the boundary condition at the surface is one of ponding, i.e., if the water supply rate exceeds the maximum potential infiltra­ tion rate, then one must assume that 0b = φ{\ — Sar) where φ is porosity and 5 a r is residual air saturation. This assumption is necessary when Richards' equation is used to describe the soil water movement because that equation provides no mechanism for air counterflow (Morel-Seytoux and Noblanc, 1971). Physically this assumption is incorrect, but it is the only mathematically acceptable one when Richards' equation is used. With a boundary condition of ponding at the surface, the explicit scheme of Eq. (43) should be replaced by an implicit scheme; i.e., instead of Eq. (43) the following equation should be used:

<7m = *m° + {

K"(6b) + K"(0J1

\D"(eb) + - '

2σ„

D'WJ

Δί

(50)

where 0b = φ(1 - SJ

(51)

In this case Eqs. (48) and (50) must be solved simultaneously for the two unknowns 0m and
H. J. MOREL-SEYTOUX

A. MATHEMATICAL FORMULATION OF CONSERVATION LAWS

First let us state the assumptions for the case considered. The porous medium is homogeneous and isotropic. Its porosity and permeability are constants also with respect to time. The soil column has a uniform cross section. Then Darcy's law can be written for each phase (say, water and oil):

,

w

=

. ^ ^

U

^

M

(52)

kTO op0 kT0 V = ~£ — "^- + k—p0g (53) μ0 dz μ0 A water-air system would be treated in the same way. Again provided that air compressibility is negligible, all the results obtained for the water-oil system will be applicable to the water-air system. All that is required is to substitute the subscript a for the subscript o. It is convenient to define a total velocity v as v = vw + v0

(54)

and to define the water fractional flow as Fw = vjv

(55)

Since the two fluids have been assumed incompressible, the equations of mass conservation for both fluids are φ dSJdt + dvjdz = 0

(56)

φ dSJdt + dvjdz = 0

(57)

where φ is porosity. Adding Eqs. (56) and (57) and noting that the sum of the two saturations is always 1, one obtains the result that the total velocity does not depend on z. (It may depend on t, however.) In other words, the system of Eqs. (56) and (57) is equivalent to the system φ dSJdt + dvjdz = 0 v = v(t)

(58) (59)

This system can be rewritten, using Eq. (55), as φ dSJdt + v dfjdz

=0

(60)

Of the three variables 5 W , z, and i, only two can be viewed as independent variables. Viewing Fw in Eq. (60) as a function of 5W and t, Eq. (60) can be rewritten as φ dSJdt 4- i>Fw' dSJdz = 0 (61) 140

Two-Phase Flows in Porous Media where F w ' is a short-hand notation for 3Fw/<3Sjf. Using Eqs. (52) and (53) and the capillary pressure, Eq. (17), to eliminate p w and p0, one can deduce an explicit form for the general fractional flow function. Multiplying the two sides of Eq. (52) by pjkkrw, the two-sides of Eq. (53) by —p0/kkro and adding the resulting equations, yields kkrw

kkrn

dz

dz

(62)

or μ*, F W

μ„(ΐ - Fw)

i (dp.

-\®+«*)

(63)

kL· kL· which can be solved for Fw with the result

F

-4i+£(£♦*»)]

(64)

= (ι + ^φΓ

(65)

where Δρ = p w — p0 and/ w is a function of saturation only (i.e., not a function of z or i), namely

A

Figure 11 shows typical/ w curves for a given set of relative permeabilities and different viscosity ratios. The relative permeabilities used to obtain the curves in Fig. 11 were (66) kTa = 0.3(SW*)3 kro = 0.75(1 - Sw*)2 i.o

0.8

0.6

0.4

0.2

FIG. 11. Effect of viscosity ratio on the fractional flow function, function / w 141

(67)

H. J. MOREL-SEYTOUX

where Sw* is the normalized saturation defined as ^w

=

(^w ~~ S wr )/(1 — o w r — Sor)

(68)

Within the range of variations in saturation by liquid displacement, the normalized saturation varies between 0 and 1. Figure 12 shows curves of lO.Of-

- V 40

FIG. 12. Effect of viscosity ratio on the derivative of the fractional flow function/w.

/w' = dfJdS^ for the same viscosity ratios. These curves will prove of interest later. Equation (64) can be rewritten as

K =/w{l + [kkT0^0v][P:{dSJdz) + Apg]}

(69)

where the prime, e.g., as in pc\ indicates differentiation with respect to saturation. By substitution of Eq. (69) into (61), an explicit partial differential equation for the unknown saturation is obtained. Note however that this equation contains another unknown, namely v. It is the purpose of the next section to derive an equation for v. B. THE INTEGRAL EQUATION

From Eqs. (52) and (53) an equation for v in terms of pressures can be obtained, namely , , , , \ , , , = ~ T 2 + Po9 + / w ^ +/w Δρ<7 (70) KKJl·** + KolVo) dz dz This equation is valid at all times and anywhere within the porous medium. In particular at a fixed time both sides of Eq. (70) can be integrated with respect to z between any two positions denoted 1 and 2, with the result: 142

Two-Phase Flows in Porous Media 2

v

/

2

7

f ΊΤΓ~~ί

1 ;

2

/ Λ = Poi ~ Pol + Po^(^2 - Zi) + \ /w <*pc + Apfif f / w dz (71)

=

(Pol - P o ^ l ) - (Po2 ~ P o f l ^ ) + J?/w ^Pc + Apflf J ? / w rfz

This last equation is of value because even though the pressure distribution is not known, often values of pressure are known at particular locations (e.g., on the boundaries). Thus several terms of Eq. (72) can be evaluated exactly, whereas the others have to be approximated. Of course Eq. (72) can be written in terms of water pressures instead of oil pressures, with the result (Pwi ~ Pw2) + Ρ*Φι ~ zi) ~ &P9 JiO ~/w) dz - \\{\ - / w ) dpc ί IKKJP* + KolPo)]'1 dz In the numerator the first terms, i.e., p w l — /? w2 , expresses the drive mecha­ nism due to an external pressure gradient. This term will be very high, relative to the other terms, when a fluid is injected to displace another under very high pressure. This is a situation frequently encountered in waterflooding opera­ tions. The second term expresses the gravity drive of water. It will be the dominant term under conditions of ponded infiltration for large times. For intermediate times, the third term (expressing buoyancy) will not be negligible as the water flow fraction throughout the soil column is different and smaller than 1. The fourth term expresses the capillary drive Finally, the denominator of Eq. (73) expresses the viscous resistance to flow. In the next section, the problem of imbibition into a horizontal core of infinite extent (i.e., large L) will be discussed again but from the two-phase point of view. C. IMBIBITION INTO A THIN HORIZONTAL SLAB

Figure 13 shows a graphical representation of the problem. At time zero the air saturated slab is placed in contact with the water tank. From this initial time on water imbibes into the slab. It will be assumed for simplicity that the slab extends infinitely to the right. In this case the saturation equation takes the form φ dSJdt + d(v¥w)/dx = 0

(74)

Fw =7w[l + (**„/*!, v)pe' dSJdx]

(75)

where

143

H. J. MOREL-SEYTOUX Atmospheric Pressun ■llCT . "Φ-

Water Tank

Oil Saturated Thin Slab

Atmospt Pressur

L X

(a) Sw k

(b)

FIG. 13. Water imbibition into an air-filled core, (a) System geometry; (b) later time saturation profile.

Taking location 1 at x — 0 and location 2 at x = oo in the integral equation yields V =

Jo[fr(frrw/^w +

fcra//0]

* dx

(76)

Several comments are now in order. After substitution of Eq. (75) into (74), the resulting equation has the form ,3S

d Γ

k ,YC,5SW]

(77)

where for brevity, a new function is defined by (78)

Equation (77) with the new notation takes the form +

*ir of

+

^r dx

"H s wrr° ox L/*a 144

dx J

(79)

Two-Phase

Flows in Porous

Media

Except for the second term, this equation has the same structure as Eq. (25). One can expect a solution of the stretching kind particularly since only the denominator (resistance) term of Eq. (76) varies with time. Let (80)

x = u(S)tn

where the w subscript is dropped for convenience since no confusion is possible. By expressing the differential of any function ( ) in terms of the independent variables x and t or S and t, one obtains the identity dx +

dx

dt =

dt

öS

dS +

en dt

dt

(81)

and from Eq. (80) one has dx =

du

t" dS + nu(S)f-1 dt

Ö~S

(82)

Substitution of this expression into Eq. (81) and identifying coefficients of the independent variables dS and dt yields dx

,dS du , dS

= Γ

(83)

and

d

A2

dt

d

12

dt

, dS -nw(S)i_1 — s ou

(84)

dS

Using these operators Eq. (79) transforms to -nu(S)t~l

dS du

+ f/w'i""

dS

k

2dS

du

μ3

du

E ( S SS\~[ )

I d Γ

U[

d,H (85)

which simplifies to μΛ

öS

£(S)ash

dJ=°

(86)

From Eq. (76) and using the assumed form of the solution given in Eq. (80) one obtains (87) All the integrals are not functions of time. The integral in the numerator can be readily calculated from basic properties of the soil and of the fluids. The 145

H. J. MOREL-SEYTOUX

integral in the denominator can be calculated once u(S) is known. Even though u(S) is not yet known, the functional dependence of v on time is known as v(t) = C/Rf = oc/f

(88)

where C is known and both R and n are unknown constants. Substitution of Eq. (88) in (86) yields

-φηηφΓ1

+ a / w ' r 2 M + - r-in2 n ^ \E(S) ^ k

Ö

μ3

30 SS L

1 =0

du U

(89)

This equation will reduce to an ordinary equation for the choice of n = \. Finally, the equation for S and u is d_ du

dS'

6u k

Ί dS

T-* s

(90)

This equation has to be solved numerically. The solution u(S) will depend on a = C/R and the proper value of a is the one that satisfies the compatibility equation C = rsu'(S,a)dS cc Ji-s ar /c(/c rw /// w + /cro//i0) An alternative method of solution was devised by McWhorter (1970) in which the fractional flow function was viewed as the basic unknown function rather than u(S). The numerical procedures associated with the method have been found to converge rapidly (McWhorter, 1970). D. DISPLACEMENT PROCESS IN A HORIZONTAL CORE

This case corresponds to the situation when a fluid is injected at high pressure into a core originally saturated with another fluid. In this case the pressure drop between the injection and the production faces is large. It means that in Eq. (76) the term p^gH is large compared to the capillary drive term. It seems that after the neglect ofthat capillary term in the integral equation, one could proceed as before. However, if the capillary terms are neglected not only in the v equation but also in the expressions for i;w and vQ, the saturation equation (Buckley and Leverett, 1942; Morel-Seytoux, 1969) takes the simpler form (assuming v different from zero) φ dSJdt + divfj/dx

=0

(92)

or φ dSJdt + vfj dSJdx = 0 146

(93)

Two-Phase Flows in Porous Media from which the velocity of propagation of a given saturation is deduced (Buckley and Leverett, 1942; Scheidegger, 1960; Morel-Seytoux, 1969) as (

(94)

The location of a particular saturation 5W at a time t is thus (95)

x = Xo + (fJ/Φ) f' υ(τ) dx = x0 + WI4>)W{t)

where W(t) is the cumulative volume influx of water per unit area. It follows from Eq. (95) that the shape of the saturation profile is essentially the same as that of/w' versus S w . It is clear from Fig. 12 that for a given/ w ' value there exist two saturation values. As a result, the saturation profile determined from Eq. (95) shows multiple values of saturation at a fixed x (see Fig. 14). Thus, the strict application of Eq. (95) is not physically acceptable. The problem of multiple values was recognized by Buckley and Leverett (1942) who suggested the introduction of a front (see Fig. 14) to replace the multiple-valued part of the profile. The upper value of the saturation front is denoted S B L . Its value can be obtained by drawing a tangent (Welge, 1952; De Wiest, 1969) from the point of abscissa Swi to the/ w curve (as shown on Fig. 14). From this construction it follows that the following relation must be true: /BL = / B L / ( S B L - S„)

(96)

From Eq. (95) one obtains at a fixed t dx = (f:/)W(t)dS

(97)

Substitution in Eq. (76) yields

P*gH+tiriS*"V*dPe JJ i-s

/cro//i0) o r ^(/c r w //i w + KolVo) Φ i -Sor KKvlVv,

kkTOl

(98) φ

J

where L is the length of the slab. Equation (98) is valid only prior to break­ through of water at the producing face. Though this expression may look formidable, let us note that every item involved can be evaluated in terms of the properties of the soil and of the fluids. For convenience, let us define J μ0 s

-

-ns)äs

(/crw//iw + Κ0/μ0)

In particular, let ABL = A(5BL). Let us define if = P.SH+ f" 147

" /„ dp,

(100)

H. J. MOREL-SEYTOUX

—X

1.0

T 1 1 1 1

0.8

1

0.6 0.4

j

1 1 I I 1 1 Γ

0.2 h

l·

/

/

1 1

/ /

1 ^wr

/ !

Ut-ς

y

uL^r

Swi

/

/ /

'/ 0.2

/ / /

/

/

7

i

// / / /

i l SBL

0.4 0.6 Saturation

i .. 1 .

0.8

j

-H

1

1.0

FIG. 14. Multiple-valued saturation profile and the Welge graphical construction for the front.

Then Eq. (98) takes the form ν = ΓΑρ

(101)

where Δρ is the driving mechanism (pressure) for flow and Γ the conductance is [^(0/^][ABL-^L]+L

(102)

or more simply Γ =

Kok/ßo a 148

bW{t)

(103)

Two-Phase Flows in Porous Media Since v = dW/dt one obtains a differential equation for W which is easily integrated with the result for v: krokApl

2kk~roAp t t l

f -w*

It is clear that for some times (t = 0) the resistance to flow is due to the sole presence of oil since hardly any water has as yet imbibed. Note that the total velocity v (also the imbibition rate of water in this case) is not infinite at time zero. In fact, the imbibition rate increases with time because generally f^L is greater than ABL when μ0 exceeds μ^. The initial imbibition rate will be infinite if the viscosity of the displaced fluid is assumed negligible, which is unrealistic for oil but reasonable for air. Even in the case of air Eq. (104) (now with the subscript a instead of o) may provide a better approximation than the unsaturated flow solution particu­ larly for a high noncapillary pressure drop across the slab. Equation (104) has the physical advantage over the unsaturated solution of predicting a finite imbibition rate at time zero. In the unsaturated flow solution, the length of the slab does not appear in the equation for the imbibition rate. From Eq. (104) it can be seen that the square-root of time law of decrease of v is pronounced for small cores but masked for long ones. Ironically the unsatu­ rated flow solution for a semiinfinite core provides a better approximation when the slab is short than when it is long! For a drive which is primarily of capillary nature, it may be necessary to solve Eq. (90) for u(s) and then Eq. (87), specially if a good determination of the moisture profile is wanted. If the imbibition rate is the prime quantity of interest, Eq. (104) is probably sufficiently accurate (this assertion will be substantiated in Section V).

E. DISPLACEMENT PROCESS IN A HORIZONTAL SLAB WITH A CLOSED END

In this case since v is not a function of z and is zero at the end of the core, v is zero everywhere at all times. Thus the integral equation is not needed since v is known. The fractional flow concept however is no longer applicable since division by zero is not permissible. The saturation equation can be written as φ dSJdt + ej dSJdx = 0

(105)

where ej is shorthand for deJdSJ[t and *w = - U P . ' ψ μα °χ

149

= - £(SW) ψ μΛ dx

(106)

H. J. MOREL-SEYTOUX

As before the rate of propagation of a given saturation is (dx/dt)Sw = e^'/φ

(107)

The solution of the problem is the same as for Eq. (77) with the provision that v = 0. The solution for u(s) is that of Eq. (90) with a = 0, or d Γ du L

dS~\ μΛφ dS du\ 2k du

This equation is similar in structure with Eq. (27) of the unsaturated flow theory but there is a fundamental difference. The diffusivity is a monotonic function of Θ whereas E(s) is not. This means that there exists a saturation Sb for which u(s) = 0 which is equivalent to say that ej is zero at all times for some saturation Sh. Only saturations between Swr and Sb propagate into the core. Water moves in and air moves out and consequently at the inlet (on the slab side) the saturation of water is less than (1 — Sar). Because the solution is of the form x = u(s)t1/2 the cumulative imbibition volume increases like t1J2.

F. INFILTRATION INTO A VERTICAL COLUMN

The system of equations describing this problem consists of Eqs. (60), (64), and (72) Equation (64) will be rewritten as Fw = F w + (k/μ. v)E(SJ dSJdz

(109)

F*=f*[l+(kkJw)Apg]

(110)

where:

Fw is the fractional flow function when capillarity is neglected but gravity is retained;/ w is the fractional flow function when both capillary and gravity effects are neglected. Whereas/ w was only a function of saturation, F w is a function of saturation and of v, and therefore a function of time. Figure 15 shows several F w curves corresponding to several values of v. The curves of Fig. 15 are obtained with the same relative permeability data used to obtain Figs. 11 and 12 and with a viscosity ratio /iw//ia = 100. The complete saturation equation is somewhat complicated by the presence of the capillary term. If one is primarily concerned with the infiltration rate rather than with the moisture profile, it would seem intuitively reasonable to neglect the capillary term in the saturation equation where its presence com­ plicates matters and to retain it in the integral equation where its presence 150

Two-Phase Flows in Porous Media

0.5

02

0.3

0.4 0.5 0.6 Water Saturation

0.8

0.9

1.0

FIG. 15. Effect of total velocity on the fractional flow function Fw for a viscosity ratio of 100.

introduces little complication. Under these conditions, the saturation equa­ tion has the simple form (Brustkern and Morel-Seytoux, 1970) φ dSJdt + vFJ dSJdz = 0

(111)

for which the solution is z = z0+

|VW70MT)A

(112)

giving the position z of a given saturation Sw originally located at z0 at time zero, provided Sw exceeds the Buckley-Leverett value 5 B L or 1

zf = zn +

ν(τ) άτ

(Π3)

giving the location zf of the front and where the superscripts — and + refer to conditions just ahead and just behind the front. We shall consider the case of a vertical column of depth D under conditions of constant depth of ponding at the surface H and atmospheric pressure for the air at the bottom of the column. It is assumed initially that the water saturation in the column is uniform and residual. Under these conditions the integral equation takes the form: _ PwgH + PagD + J?^ s . r /wPc' dS + Apg J?/ w dz (114) Jo[fe(ferw/^w +

151

fcra/^a)]

' ^Z

H. J. MOREL-SEYTOUX

or, employing the approach analogous to that used to derive Eq. (98), v = W(t) rD kpqt r SeL rSwr dS ; / w dpc + Apg -f/ w / " dS + -~- f /WG" Py,gH + PagD + I rSBL

G"c/5

t J»SBL

/ " dS

Φ λ-Sar ΚΚ^Ιμ* + kJva)

ψλ-5.Γ

KKJVv,

+ ^ra/Ma)

+

Mi

77HÖ-*BL] fcfc

(115) Whereas Eq. (98) could be integrated exactly to yield Eq. (104) for v, Eq. (115) cannot be integrated so easily because / depends on v and consequently also S BL which is thus a function of time. An alternative procedure is to start with an initial v denoted v0 defined by kk vo = -η:

o°Γ

Γ PwffH +

°· 2

0.4

L

1

rSwr PagD

1

+

0.6 1

1

0.8 1

1 fwPc' dS\ 1 1

(116)

1.0. f I

5.86hr

1000 S

wr

!

c

H °ar

2000

3000 139.20 hr

1

4000

5000

■

6000

-

7000

-

8000

1

214.20 hr

339.20hr Qnnn z(cm)

FIG. 16. Advance of saturation profile in semiinfinite medium (Brustkern and MorelSeytoux, 1970). 152

Two-Phase Flows in Porous Media

250

E υ

100 150 Time ( h r )

200

250

FIG. 17. Infiltration rate and cumulative infiltration versus time fr.r a semiinfinite medium (Brustkern and Morel-Seytoux, 1970).

For a ponding of 1 cm and a 3-m sand column, v0 will be of the order often times the hydraulic conductivity. Given v0, it is possible to find the evolution of the saturation profile after a small time increment At from Eqs. (112) and (113). Once the new profile is established, a new v is calculated directly from Eq. (114) since 5w(z) is now known. The procedures are then repeated for a new time step and so forth. Figure 16 shows the saturation profiles for various times down a column of semiinfinite extent. Figure 17 shows the infiltration rate and cumulative infiltration for the same case. It is interesting to note that the limiting value of the infiltration rate is below JKW, i.e., the hydraulic conductivity at residual air saturation. This is due to a counterflow of air (upward movement of air) that escapes through the surface as shown by the desaturation occurring in the soil zone close to the surface (Fig. 18) in spite of the fact there is ponding of water above the soil surface. Because v decreases 153

H. J. MOREL-SEYTOUX

with time, sooner or later the JFW curve will exceed the value 1 at high values of saturation (Fig. 15) which means that va is negative. As time increases and v continues to decrease, the boundary of the zone of air counterflow which first developed close to the surface will propagate downward. More and more of the column will experience counterflow. Though this result is based on the study of a fractional flow function in which the capillary term has been 0.880

0.890

0.90

0.910

0.920

Q930

0.940 Sw

1000

2000

3000

4000

5000

6000

7000

8000

9000 Z(cm) FIG. 18. Advance of saturation profile. Enlargement of a section of Fig. 17 (Brustkern and Morel-Seytoux, 1970).

neglected, the result is still valid because the zone of air counterflow is the very one where the capillary gradient is small (dSJdz ~ 0 and Sw is large) and in the process of infiltration (dSJdz < 0) for a given 5 W , Fw is greater than F w . Of course one can (and probably should) question the validity of the approx­ imation involved in obtaining the solution. Figure 19 shows a comparison 154

Two-Phase Flows in Porous Media between the solution obtained by the method herein presented and a solution of the nonsimplified system of Eqs. (56) and (57) by a finite-difference tech­ nique (Phuc, 1969). The approximation, at least for prediction of the infiltra­ tion rates, is obviously good.

Time ( hr )

FIG. 19. A comparison of infiltration rates as obtained by the Brustkern and Phuc procedures (Brustkern and Morel-Seytoux, 1970). G. DRAINAGE IN A VERTICAL COLUMN

It would not be realistic to discuss the infiltration of water into an air-filled column of limited depth with a closed bottom without including the effect of the compressibility of the air. However, it is interesting to discuss a similar case, that of redistribution of a water-oil mixture (Martin, 1960). In this case both fluids can be regarded as incompressible. In the case to be discussed the initial saturation profile is assumed to be uniform of value 5 w i . Both ends of the core are sealed. At time zero the core that was laying flat is brought in a vertical position. The question is: What happens? Can the evolution of the saturation profile be predicted? Again the problem can be solved if the capillary effects are neglected. Without making any approximation it can be stated that since v does not depend on z and since both vw and va are zero at the top and bottom of the column at all times, it follows that v is zero throughout the column at all times. Since v is known to be exactly zero, there is no need for the integral equation. However, the fractional flow concept cannot be used as defined 155

H. J. MOREL-SEYTOUX

previously because the division by zero is not permissible; but the limit of vFw for v = 0 exists. It will be denoted Gw and is in fact nothing more than vw. The notation Gw will be used to emphasize the similarity of the procedures with the case where Fw was defined. The equation for water saturation is φ dSJdt + GJ dSJdz = 0

(117)

from which the velocity of propagation of a given saturation is deduced: (dzW)Sw = G„'l4>

(Π8)

γ/7///////////////////// \Z (cm) (a)

j6«

\ / /

// ! M

S wl

h

M

a >

If

1 1

»

i

!

j_

\j \b

Sor

vv !

A >*L

! -1..

^

(b)

FIG. 20. Water-oil mixture redistribution in a sealed vertical column, (a) Saturation profile; (b) Gw function. 156

Two-Phase Flows in Porous Media where Gw is a function of saturation defined as G

kApgkrwkra /^w kra + kry/ μ3

^kApgknfw μ3

Figure 20 shows the function Gw and the evolution of the saturation profile. Throughout the system the water flux is downward and the oil flux is upward. A wetting front propagates upward and a drainage front propagates down­ ward. After some time the two fronts will run into each other. Then they will merge into a new front of bigger size (graphically displayed as ab on the Gw curve). Since the slope ab in Fig. 20b is slightly negative the new front will propagate upward slowly. It is a wetting front. Had the slope of ab been positive it would have been a drainage front propagating downward. The evolution of the profile after the merger of the two fronts is not as simple as before because the new front (ab) is now moving into a zone where the saturation is no longer uniform. The front will grow in size and gradually slow down to the equilibrium position when oil will sit on top of the water.

V. Liquid-Gas Flow Theory In the previous section the infiltration of water into an air-filled column was studied under the assumption that air was incompressible. For a long column open at the bottom, it seems intuitively that the assumption may not invalidate the results. For a closed short column, however, both the movement and the compressibility of air probably affect the results. For this reason, the air is no longer assumed incompressible in this section. One can either try to treat the problem exactly or view the compressibility effect as a correction term which need not be treated with great refinement. The second approach will be discussed in more detail.

A. MATHEMATICAL FORMULATION OF CONSERVATION LAWS

The equation of mass conservation for water (assumed incompressible) is still given by the usual equation: φ dSJdt + dvjdz = 0

(120)

However, the equation of mass conservation for air takes the form φ d(pa Sa)/dt + d(pa va)ldz = 0

(121)

and it is no longer true that v is only a function of time and not a function of z. 157

H. J. MOREL-SEYTOUX

B. THE INTEGRAL EQUATION

In the liquid-gas flow theory the proper form of the equation is a slight modification of Eq. (72), namely r2

v(t,z)dz

r2

,

r2

+ g fp^l-fjdz.

(122)

Even though v varies with z there may exist zones of almost uniform v separated by a sharp transition zone of steep change in v. If positions 1 and 2 refer to the extremities of a zone in which v changes little, then as an approxi­ mation one obtains for an average v in this zone: - = (Pai - Pai) + f?/w dPc + ^P9 f i/w dz V

ilWJ^ + kJriridz

{

}

C. IMBIBITION INTO A HORIZONTAL SLAB

Let us consider the problem illustrated in Fig. 13 when the right-hand side extremity of the slab is sealed and no flow can take place at this end. Utilizing the method of solution discussed in Section IV,D and integrating from the left-hand side of the slab to just ahead of the wetting front, one obtains - _PV,QH

V

+ PA - A»2 + ίι/w

\m^ Φ

dpc

K }

fids *i-Sar/c(KJPv, + KJPa)

where P A is atmospheric pressure. The only unknown is the air pressure pa2 just ahead of the front. Since the air displacement is almost complete behind the wetting front and no flow of air takes place at the end of the slab, one can assume that the air pressure ahead of the front is uniform and its value at time t can be obtained from the perfect gas law, namely (125)

**-L-w\tm Substitution in Eq. (124) yields . _ p„gH + ff / w dpe - p A W(t)^L

-

W(t)Y'

(126)

The numerator has three terms. The first two terms were previously denoted Ap, a constant when the head of water H is maintained constant, which is 158

Two-Phase Flows in Porous Media assumed for simplicity in the following discussion. Equation (126) indicates that v will take the value zero when the cumulative infiltration has reached the value W(t)/4>L = Ap/(pA + Ap) = y

(127)

which will eventually occur. The rate at which v goes to zero is of interest. Let W0 be the value of W satisfying Eq. (127). The solution for W* = W/W0 as a function of time is obtained by integration of Eq. (126) (with a result analogous to that of Adrian and Franzini (1966) who assumed a piston­ like displacement): yW*2

_ 2(1 _ y)w* + 2(1 _ y)in(_l_) = 2^1

(128)

'

\\ - W*f W0^aABL This is a transcendental equation which cannot be solved for W*. It is easy enough though to calculate values of t for various W*. Close to 1, small changes in W* will induce large changes in t, which implies that v goes to zero when W* approaches one. Of course this result was already known. In terms of W0 the expression for v takes the form v = (W0-

W)/cW((l)L - W)

(129)

where € = μΆΑ^/φΚΑρ

+ ρΑ )

(130)

From Eq. (121) one can calculate dv/dt, namely dv =-v[W2 dt

-2WW0 2

cW ((l)L-W)

+ W0L] 2

(

}

which evaluated in the neighborhood of W0 gives dv/dt = - v/cW0(
(132)

That is, v approaches zero asymptotically for large times. Once v reaches the value of zero, the situation is analogous to the problem discussed in Section IV,E, i.e., a case of capillary counterflow. D. INFILTRATION INTO A VERTICAL COLUMN WITH A CLOSED BOTTOM

The problem will be discussed from several viewpoints. First the techniques discussed in Sections IV,F and V,B will be combined to obtain a solution. Then a more refined method, which can be called the moving strained coordi­ nates method, will be discussed. Finally, results obtained by a finite-difference method are presented. 159

H. J. MOREL-SEYTOUX

1. The Brustkern Procedure In this approach the saturation equation to be solved is φ dSJdt + vFJ dSJdz = 0

(133)

and v is obtained from the integral equation (124). In other words the capil­ lary terms are neglected in the saturation equation but they are retained in the integral equation. The effect of air compressibility is reflected in the term pa2, which until counterflow occurs is given by the perfect gas law p 2

*

=

PAJ^

m4

{io

D-w(t)/
x }

The saturation profile is obtained by moving every saturation over a time step Δί by the amount Αζ = (ν/φ)Γ„Άί

(135)

and by moving the saturation front by the amount

From the definition of the fractional flow function F w [Eq. (110)] for large v the Fw curve deviates little from/ w . However for small v the F w function can become large and in particular will exceed 1 in the upper range of saturations. Typical Fw curves for several v values are shown in Fig. 15. If F w exceeds 1 it means that va is negative, that is air is flowing upward or counterflowing. A case of counterflow was already encountered in Section IV,G where v was always zero. Counterflow will occur when the primary drive mechanism for water is gravity and that for air is buoyancy. The initiation of counterflow starts at the time when the curve F w just exceeds 1 for some high values of saturation. As v decreases further, more and more saturations travel upward and the zone of counterflow reaches further downward. Shortly after the beginning of counterflow, air starts to escape at the surface and the air pressure ahead of the wetting front is no longer given by Eq. (134). A material balance coupled with the perfect gas law determines the new air pressure PAK Pa2

φϋ{\

-j—-

[£>( 1 - Swl) - j\i

+ Pc(Su)/pA)(I - v) άτ^

(137)

where 5 U denotes the saturation of water at the upper boundary on the soil side. Equation (137) is valid for all times because prior to air escape I — v = 0 160

Two-Phase Flows in Porous Media and Eq. (137) reduces to Eq. (134). Figure 21 shows how the infiltration rate is affected by the depth of the porous medium. It is interesting to note that the infiltration rate drops much below the hydraulic conductivity at residual air saturation X w , then recovers and asymptotically seems to reach a value significantly below Kw. This is understandable because the air continually escapes from the medium and the water saturation at the upper boundary is less than (1 — Sar). A look at typical kry/ curves shows that kTW decreases rapidly for small desaturations around the value S w .

Time (hr)

FIG. 21. Comparison of infiltration rates when depth D of water table varies (MorelSeytoux, 1971). Air counterflow occurs much before the wetting front reaches the bottom impervious boundary. By that time v has been zero for a long time already. The saturation profile evolution when the wetting front hits the bottom is similar to what happens in the case discussed in Section IV,G. Figure 22 shows a graphical display of the construction of the saturation profile, the reflection of the downward front at the boundary and the propagation of a new front wetting the soil from the bottom up (Brustkern and Morel-Seytoux, 1972). 2. The Nob lane Procedure The Brustkern procedure seems susceptible of generalization for more complex problems (e.g., the cases of nonuniform initial saturation, realistic hydrologic boundary conditions of intermittent rainfall, etc.). Some measure of the validity of the procedure was already obtained in Section IV,F. Neverthe less the neglect of capillary terms in the saturation equation may cast some doubt on the method. The Noblanc procedure is somewhat similar in concept in the sense that it also uses the fractional flow function idea and the integral equation. It differs 161

H. J. MOREL-SEYTOUX

much in the treatment of the saturation equation. The technique used was described already in the unsaturated flow theory discussion of Section III. Proceeding similarly (Morel-Seytoux and Noblanc, 1971) one obtains for the saturation profile in the front region dz dS

sE(S) FJ(S$S - F w + C(t)

(138)

where for brevity ε = k/νμ^ is defined.

Water Table

Z(cm) FIG. 22. 1972).

(e)

J

The development of an upward moving front (Brustkern and Morel-Seytoux,

In the soil surface region, to the contrary of what happens in the front region, the profile does not move as a translation but rather stretches with time. In this region curvature of the profile is not pronounced (Swartzendruber, 1969). Also as results have shown, the extent of this zone into the medium is limited to a small region near the surface. In fact, after counterflow of air occurs that zone disappears. For these reasons, Eq. (61), which can be equivalently rewritten as dz/dt = (ϋ/φ)¥„' 162

(139)

Two-Phase Flows in Porous Media can be approximated as dz/dt = ( i # ) [ F w ' + eE'(S) dS/dz]

(140)

or d[dz/dS] dt

= (v/)[F: + E'(S)dSldz]

(141)

where F$ is defined by F$ = (dFw'/dS)t = (d2FJdS2)t. To obtain the complete saturation profiles, Eqs. (138) and (141) must be solved in their respective regions and the two profiles matched at the interface of the two regions. Numerically the solution cannot be obtained without specifying the boundary and initial conditions. The computational procedures are the same as in the Brustkern procedure. The total velocity v is obtained from the integral equation. Given v, now Eqs. (138) and (141) must be solved. Given v note that two unknown parameters still appear in Eq. (138): 5 ζ and C(t). These are eliminated by the requirement that the profiles calculated from Eqs. (138) and (141) must match smoothly at a point of saturation Sm (matching condition) and by the requirement that the net water influx into the soil during the time interval At as calculated by the solution of Eq. (141) must equal the change of volume under the matched composite profile obtained from Eqs. (138) and (141) (balance condition). The expression of these requirements permits the reparameterization of the problem in terms of Sm alone and of the slope of the profile at S = Sm, denoted am ( = dz/dSm). Since without air counterflow the procedures are fundamentally the same as the one described in Section III, only the situation following initiation of counterflow will be discussed. As long as there is no air counterflow the saturation at the surface has the value 5 b r . The slope of the saturation profile near the surface om is negative. During the infiltration process \om\ increases with time. Then, before counterflow occurs, |<7m| decreases, which means that drying takes place. Finally, am becomes zero. When this occurs, and thereafter, the entire profile of saturation is governed by Eq. (138). The matching point has moved all the way back to the surface. The change of area (water content) takes the form • L FV'(S#S

- Swl) - FW(S) + Fw(Swi)

'

0

where A is the area at the end of the previous time step and the balance condition is

JF W (SJ + ε-ψ^ - Fw(Swi)J v At = φ ΑΑ 163

(143)

H. J. MOREL-SEYTOUX

3. Results of the Noblanc Procedure a. CASE 1. In this first situation the initial water saturation is assumed to be uniform and equal to the residual water saturation. At the upper surface it is assumed that water is ponded to a constant depth of 5 cm throughout the duration of the infiltration process. The medium is considered to be semiinfinite in extent. The results of a computer analysis of this problem (Noblanc and MorelSeytoux, 1972) are illustrated in Figs. 23 and 24. The curve of rate of infiltraSemi-infinite Medium x Occurrence of Air Counterflow

.£

5

- Hydraulic Conductivity at Residual Air Saturation

2h 50

100 Time ( hr)

150

200

FIG. 23. Rate of infiltration versus time for a semiinfinite medium (Noblanc and Morel-Seytoux, 1972).

tion versus time is shown in Fig. 23. During the first 3 hr of infiltration the rate is very large because the capillary effect is important. During the next 9 hr the rate decreases, and is of the same order of magnitude as the hydraulic conductivity at residual air saturation. It dips slightly below the hydraulic conductivity because even in a semiinfinite medium there is a zone of in­ creased air pressure behind the wetting front. As a result, even though the soil surface is still fully saturated (except for residual air saturation), the hydraulic gradient is less than unity. So to speak, the weight of the column of water is not moving it down as fast as it would in a vacuum because it is running into a cushion of air. After 12 hr of infiltration duration, air counterflow occurs. The rate of infiltration decreases further, this time because there is slight water desaturation at the soil surface with a consequent reduction in hydraulic conductivity, 164

Two-Phase Flows in Porous Media b

wr

Q4

0.5

" "

' Γ ~

Ί

0.6 I

0.7 1

f

0.8

0.9

I

30,000 sec 8.33 hr

1.0

10

200,000 sec 55.55 hr

20

3 0 0 , 0 0 0 sec 83.33 hr

4 0 0 , 0 0 0 sec I I I . I I hr 30 5 0 0 , 0 0 0 sec 138.88 hr 6 0 0 , 0 0 0 sec 166.66 hr

40

7 0 0 , 0 0 0 sec 194.44 hr

J y

1

/ Spr

J A

J J

FIG. 24. Advance of saturation profile in a semiinfinite medium (Noblanc and MorelSeytoux, 1972).

and tends to stabilize at a value lower than the hydraulic conductivity at residual air saturation. Figure 24 shows the evolution of the saturation profile versus time. Even though capillarity is taken into account, the profile is very similar to that obtained by solving Eq. (133). Except for the larger values of water saturation, the motion of the saturation profile is essentially a translation. b. CASE 2. The boundary condition at the upper soil surface, and the initial moisture condition are the same as in the first case. However, now a water table is assumed to exist at a depth D. Figure 25 shows the rate of infiltration as a function of time for different depths of the water table: 50, 150, and 300 cm. On the same figure, the cumul­ ative infiltration is shown as a function of time for the same depths. Figure 26 165

H. J. MOREL-SEYTOUX 1.5

1

1

r

1

ι

\^y

i

i

1

1

i

ix

1

1

c£>

1.0

S^ i>0.5

1

1

0.1

[1

1

1

0.3

1

1

1

1

^ D = 50 cm ^-D=l50cm ^ - D =300 cm

•E 2

1

1

0.5 0.7 Time ( hr )

I

0.9

I

I

I

I

!

^-Hydraulic Conductivity at Residual Air Saturation

LI ^

7T

0.1

1

-J 0.3

1

1 L._ . 1 0.5 0.7 Time ( hr)

1

1

1

0.9

FIG. 25. Variation in infiltration rates and in cumulative infiltration resulting from changes in water table depth (Noblanc and Morel-Seytoux, 1972).

shows saturation profiles when the water table is at a depth of 150 cm. Figure 27 displays the increase in air pressure ahead of the wetting front in the soil with time. The method for predicting the saturation profiles, and correlatively the infiltration rates, discussed here is very different from the Brustkern method. However, because the results are similar (for instance, the infiltration curves display similar dips), the reader is referred to the earlier paper (Brustkern and Morel-Seytoux, 1970) for physical interpretation of the behavior of the curves. 4. Comparison of Results of the Two Procedures In the case of a semiinfinite medium, the Brustkern procedure provides a larger rate of infiltration than in the Noblanc method; air counterflow occurs at a later time, and the rate of infiltration tends toward a constant value that 166

^wr 1

0.4

0.5 1

1

0.6

0.7

0.8

'

Γ"

,0.9

1.0

Y( ΪΙ

I 30 sec ._ Sar

4h

L

70sec.^

^Ρ^ J^y

D = 150 cm

S .= S wi

wr

FIG. 26. Progress of saturation profile with time (water table at a depth of 150 cm) (Noblanc and Morel-Seytoux, 1972). 1.060 D=50 cm

5

1.050

a

1.040

1.030 0.05 Time

0.10 ( hr)

0.15

0.20

FIG. 27. Average air pressure in soil ahead of the wetting front versus time for various depths of water table (Noblanc and Morel-Seytoux, 1972).

H. J. MOREL-SEYTOUX Semi - infinite Medium T

6 •

» 51 c o σ

4

£

6

2

Brustkern Solution Noblanc Solution Occurrence of Air Counterflow

\\

\ \ -\

t~ I

^

/

^Hydraulic Conductivity at Residual Air Saturation

x "

·

.

1 50

1

i

100 Time ( hr)

150

I

I

200

FIG. 28. Comparison of infiltration rates as obtained by the Brustkern and the Noblanc procedures for a semiinfinite medium (Noblanc and Morel-Seytoux, 1972).

is still lower than the hydraulic conductivity at residual air saturation (see Fig. 28). The same comparison is made in Fig. 29 in case of a bounded medium when the water table is at depths of 150 and 300 cm. In a sense the model of infiltration described in Section V,C,2 can be viewed as merely an improve­ ment of the Brustkern procedure rather than a new technique in itself. The same concept of fractional flow function is used, the same integral equation is used to compute the total velocity v from the knowledge of the saturation profile and the compressibility of air is treated in the same way. The improve­ ment consists in considering the influence of capillarity on the shape of the saturation profile. Even though the computed profile is an approximate one, its shape is in better agreement with the shape of experimental profiles (De Wiest, 1969). In the Brustkern procedure there exists a sharp front in the solution for the saturation profile. Figure 24 indicates that for large times (i.e., after a few hours) the profile in the wetting zone is indeed an extremely sharp front. For large times, the Brustkern procedure must be valid. In the early stages and particularly for shallow depths to an impervious boundary the profiles do not show steep wetting fronts (Fig. 26). As a result the infiltration curves are noticeably different. However, compared to the high value of hydraulic conductivity at residual air saturation, which by the currently accepted theory is supposed to be the asymptotic limit of the infiltra­ tion rate, the two curves are not very different (Fig. 29). Of course, the fact 168

Two-Phase Flows in Porous Media that the two curves are close does not mean that either is correct per se, only that the two methods are consistent. However, it has been shown earlier that the Brustkern method deviates unappreciably from a finite-diiference solution of the simultaneous unsimplified equations of flow for the water and air phases (Fig. 19). Additional results from the finite-difference solution are described in the next section. There again the ultimate infiltration capacity lies below the hydraulic conductivity at residual air saturation. Finally, it must be mentioned that the dips in the infiltration rate curves were noticed in some fairly old experiments (Free and Palmer, 1940) and in recent ones (McWhorter, 1971; Gaudet, 1972). From a hydrologic point of view this study of Noblanc confirms the value of the Brustkern procedure for predicting infiltration rates. The approxima­ tions in the procedure have little detrimental effect on the accuracy of the infiltration rates, and it has the advantage of being simpler.

0.4

0.6 Time (hr)

FIG. 29. Comparison of infiltration rates as obtained by the Brustkern and Noblanc procedures for water table depths of 150 and 300 cm (Noblanc and Morel-Seytoux, 1972). 169

H. J. MOREL-SEYTOUX

On the other hand, if the distribution of soil moisture is the information of interest, a better description is obtained by the Noblanc procedure. In particular, the procedure is still valid, and indeed simpler, for treating hori­ zontal flows in which air movement and compressibility are important; under these conditions the Brustkern procedure is probably too crude since gravity is no longer the ultimate major driving mechanism of flow.

5. The Phuc Procedure The Phuc procedure is a finite-difference approach to the solution of the problem defined by Eqs. (120), (121), Darcy's law for two immiscible fluids, and the perfect gas law (Phuc and Morel-Seytoux, 1972). It makes no assump­ tions at all, except those inherent in the finite-difference algebraic approxi­ mations of a system of partial differential equations. The computer program developed by Phuc (1969) can handle all hydrologically meaningful initial and boundary conditions, including effects of water movement, air movement, air compressibility, and hysteresis, but not of heterogeneities in soil conditions.

Time , sec

FIG. 30.

Infiltration response to successive bursts of rainfall (Morel-Seytoux, 1971). 170

Two-Phase Flows in Porous Media The versatility of this program is well illustrated in Fig. 30, which shows the infiltration response curve to two bursts of rainfall, for a soil of hydraulic conductivity of 0.55 x 10" 3 cm/sec and a depth 40 cm above an impervious boundary. The initial saturation of the soil was assumed uniform and equal to residual water saturation (0.2). Naturally the initial conditions for satura­ tion at the start of the second burst is not one of uniform saturation. Satura­ tion profiles at some important times are shown in Fig. 31 as well as the cumulative water heights for rainfall, infiltration, and ponding.

Hyetograph

13

I _L 200

_L 300

^

400

500

Time , sec

1 2 3 4

Cumulative Values

End of the First Ponding Just Before the Second Rain End of the Second Ponding Redistribution-1000 Seconds

Rainfall

300 400 Time ,sec

600

700

FIG. 31. Saturation and infiltration responses to rainfall (Morel-Seytoux, 1971).

Characteristically (see Figs. 32-35) the saturation profile is steep in the wetting zone, desaturation occurs in the soil surface zone, the air pressure builds up ahead of the wetting front, the zone of air counterflow (zone of increasing air pressure with depth) gradually extends downward and ultimately affects the entire column, and the ultimate infiltration rate is less than Ä L . 171

: Residual Air Saturation Depth : 4 9 5 cm 10 h - Rainfall

1556s 10031s 20092s 35981s ^77777,

V777777777777777777777777/

Zcm

FIG. 32. Infiltration process in the case of an impervious lower boundary (Phuc and Morel-Seytoux, 1972). 0.91 .912 Sw

S n r, Residual Air Saturation

• 1556s A 10031s o 20092 Δ 35981 Impervious Lower Boundary Depth 4 9 5 cm

FIG. 33. Enlargement of Fig. 32 showing effect of air escaping through the surface during infiltration process (Phuc and Morel-Seytoux, 1972). 172

I0 6

PA

I.I x 10

1.05x10

Pa

> Υ//////////////////7//////////////λ Z cm

FIG. 34. Air pressure curves during infiltration with a no-flow lower boundary (Phuc and Morel-Seytoux, 1972).

Ponding Occurs 1550s

12 ω 10

K w = Gravity Flow

2Rw

s |lr\— §

β

σ

6|

Impervious Lower Boundary Depth 4 9 5 cm Kw

a « ^ - ^

§ 41

•·

2 1

1

10000

1

•

··

1

t

20000

·· 1

1

30000

Time , sec

'O

Hyetograph - 10h - Rain

8

l 4 σ

£

2P

- 0 . 5 Kw

Time , sec

35600

FIG. 35. Hyetograph and corresponding infiltration curve (Phuc and Morel-Seytoux, 1972). 173

H. J. MOREL-SEYTOUX

E. EXPERIMENTAL OBSERVATIONS OF INFILTRATION INTO VERTICAL COLUMNS WITH CLOSED BOTTOMS

Several investigators have experimentally studied this problem (Free and Palmer, 1940; Peck, 1965b; Adrian and Franzini, 1966; McWhorter, 1971). Said Free and Palmer: Direct comparisons between data for the open and closed columns are possible. It took from 10 to 100 times longer to wet columns closed at the base than it took to wet columns open at the base. . . . Infiltration rates for the open columns continued to decrease as the depth of wetting increased. However, the rates for closed columns of sand, finer than the 20-30 grade decreased rapidly to a much lower level and continued to decrease, often at a very low rate, until there was a release of the confined air. Bubbles of air could then be seen escaping at the top of the column. A sharp increase in infiltration rates followed this release of air, and infiltration continued at a much higher rate until the column was entirely wet. Although a free gravitational flow of water was unable to take place against the pressure exerted by the confined air, powerful capillary forces apparently did permit a slow continued penetration of water. Unfortunately, in most experiments some characteristics of the soil are miss­ ing, prohibiting a complete comparison of calculated and observed infiltration rates. Almost complete information is fortunately available in one recent set of experiments (McWhorter, 1971) for an oil-air system. The experimental results of McWhorter (1971) are shown in Fig. 36. Translated into infiltration rate terms, the data of Fig. 36 indicate the existence of dips in the curves. The shorter the columns the more pronounced are the dips and the earlier their time of occurrence. This behavior corresponds exactly to the predictions obtained by the theory and are shown in Fig. 21. Both Free and Palmer (1940) and McWhorter (1971) also observed humps in the infiltration rates. They attributed the humps to the observable disturbance of the soil top zone as air started to escape. In a sense, from the time of air escape on, the soil column is made of two zones: a deep undisturbed layer capped by a more, porous, more pervious thin layer. The explanation is substantiated by the fact that no humps were apparent in an experiment using a sandstone, i.e. undeformable, core (McWhorter, 1971). On the other hand, Fig. 30 seems to indicate that such succession of a dip and hump in the infiltration rate curve may occur without a disturbance of the upper soil region. However, in the experiments herein discussed the disturbance of the upper soil zone was observed and seems the most likely explanation for'the hump. Figure 37 shows a comparison of observed and calculated infiltration rates 174

Two-Phase Flows in Porous Media -

o

D

a

~ -

a

1J -"»-ao

60 h

a a

O o

D-

D A

9 9 0 cm

a

o

.393cm o

°

o°

-

·

*

I

1

1

Δ

a

°

D o D o O 0 0 °D0 Λ Δ ^

30h

A

i

σ

k 401L

A

• v> Δ

Δ

-^

1 20

AA

A

^ ^ - 2 3 3 cm

X

1 8 5 cm

10 £ 0 D 5

1

0

1

1

10

1

1

20

J

30

1

1

1

40

J

L

50

1

1

60

1

1

70

80

1

90

1

1

100

1

1

1

110

1

1

120 125

Time , minutes

FIG. 36. Comparison of cumulative infiltration volume when depth of impervious boundary varies (McWhorter, 1971).

1

0.7,

Experimental Data by McWhorter

(1971)

0.6 o 185 cm Column Δ Open Column

0.5

Numerical Solution of Two Phase Flow by Sonu (1972) • 185cm Column A Open Column

Tl

0.4

0.3

1 &A 11 A \

•E 0.2

0.1

^

^

_

▲

-1

k 0 1 *

10

20

30

°

1

^

40

50

O

1

w

-r

60

1

70

1

80

1

1

90

100

Time, Minutes

FIG. 37. Comparison of observed (McWhorter, 1971) and calculated (Sonu, 1972) infiltration rate curves. 175

H. J. MOREL-SEYTOUX

by a slightly modified Brustkern procedure (Sonu, 1972). Prior to and shortly after air escape the agreement is excellent. For large times the agreement is again excellent.

F. HYSTERESIS EFFECT DURING INFILTRATION AND REDISTRIBUTION

In the unsaturated flow theory under ponding conditions at the surface there would be no need to consider hysteresis during the infiltration process. However, it has been shown that even during the infiltration process desaturation takes place in the upper soil zone and affects larger and larger portions of the column as time goes on. Thus the infiltration rate will be affected by the phenomenon of hysteresis. Figure 38 shows the effect of hysteresis on the

Without Hysteresis With Hysteresis a)

8

1o 6 σ 4 ^

2 5,000

FIG. 38.

10,000

20,000 Time, seconds

30,000

Effect of hysteresis on the infiltration rate (Phuc, 1969).

infiltration rate for a 495 cm deep medium. The infiltration rate following ponding at the surface is 2-3 % lower than when the hysteresis effect is not considered. However, at large times the two curves reattach. Thus the hysteretic effect is of minor importance in comparison to that of air movement and especially of air compressibility. Figures 39 and 40 show the effect of hysteresis on the saturation profiles during the infiltration phase and the redistribution phase.

G. EXPERIMENTAL STUDY OF THE EFFECT OF AIR MOVEMENT AND COMPRESSIBILITY IN A HETEROGENEOUS MEDIUM

In a very recent study (Vachaud et αί, 1972) it was found that "(1) the local soil air pressure can be very significantly different from the external 176

Two-Phase Flows in Porous Media 0.3

0.5

0.7

0.9 1.0

S

-

Zcm

FIG. 39. Effect of hysteresis on the saturation profile during the infiltration phase (Phuc, 1969).

atmospheric pressure, e.g., 50 mb for rain with an intensity of 3 cm/hr or 15 mb for gravity drainage and (2) the water flow depends not only on the boundary conditions in terms of water content and of water pressure, but also is strongly dependent on the pressure conditions. We conclude that the air pressure must be taken into account when determining the soil water suction, and we suggest that the flow equations must be written in terms of two-phase immiscible fluid flow." Figure 41 shows a comparison of the evolution of the moisture content profile when lateral air flow is or is not permitted under conditions of a constant infiltration flux at the top of the column. When no lateral air escape 177

H. J. MOREL-SEYTOUX 0.9

k Sw

FIG. 40. Effect of hysteresis on the saturation profile during the redistribution phase (Phuc, 1969).

is permitted, the infiltration rate from the upper layer into the middle layer is significantly reduced. Instead of proceeding totally downward only a smaller fraction of the water proceeds and the rest is reflected upward. What happens in the upper layer is similar to the situation described in Section V,D (Fig. 22) or IV,G (Fig. 20). What happens in the middle layer corresponds to the situation discussed in Sections V,D (Fig. 21) and IV,F (Fig. 16) with a variable depth of ponding due to the accumulation of water in the upper layer. In this section, except for the brief discussion of these very important experimental results, the soil was always assumed to have uniform properties. In the next section the effect of heterogeneities will be discussed. 178

Two-Phase Flows in Porous Media Water Content, cm^cm3 0.1 0.2 0.3

Water Content, cm^cm3 0.1 0.2 0.3

0.4

I20LCI0_

0.4

CI0-

FIG. 41. Constant infiltration flux for pair of experiments with or without lateral air escape (Vachaud et ai, 1972).

VI. Quasi-linear Theory for Two-Liquid Flow In all previous sections strictly one-dimensional situations for cores of uniform properties were discussed. In the case of, e.g., infiltration into the soil, the two-dimensional character of the water availability or the hetero­ geneous nature of the soil may be important. An attempt will be made to take these effects into account without actually trying to solve the actual twodimensional problem. In this chapter a water-oil system will be considered. The application of the concepts to water-air systems should be a relatively simple matter. This section and Section VII might be of particular interest to the growing number of hydrologists concerned with the migration of surface oil spills toward the groundwater. A. RESERVOIR BEHAVIOR SIMULATION

Once the basic concepts of oil recovery by displacement have been under­ stood, the next step is to apply them for prediction of actual oil-field perfor­ mances. An oil field cannot, unfortunately, be assimilated to a laboratory 179

H. J. MOREL-SEYTOUX

cylindric core, being just a little bit bigger. The exact simulation of the behavior of an actual oil field is a difficult task and one must resort to mathe­ matical models that range from rather crude to highly sophisticated. Even the most sophisticated ones, however, leave a lot to be desired. There are two basic ways of simplifying the model. One consists of reducing the configuration of the reservoir to an idealized simple shape yet retaining a fairly rigorous analysis of the displacement process. The extreme example of this approach was studied in Section IV. If the reservoir is assumed homo­ geneous and of a perfectly uniform cylindric shape, if water injection is assumed to take place at one end of the reservoir and only there and produc­ tion taken at the other end of the reservoir, then Eq. (60) of Section IV describes the performances of the reservoir. The value of such prediction will depend on the actual reservoir properties and how it is operated. The other category of simplification consists in assuming a very crude displacement process yet describing the reservoir configuration with a great deal of realism. The extreme example consists of assuming fluids of identical properties and a pistonlike displacement. Again the value of such predictions will primarily depend on how different the fluid properties are. The present section will deal with a model of the reservoir of type 1. It is a natural extension of the study of the displacement process in a strictly linear core. It differs from the latter in the sense that a touch of realism has been added. The reservoir is no longer described as a cylinder of uniform properties; it is no longer assumed homogeneous. As in an actual field, production is taken from many points in the reservoir not merely from an outlet face. The same remark applies to injection. In the displacement process the influence of gravity is included.

B. SATURATION EQUATION FOR A QUASI-LINEAR SYSTEM

A quasi-linear system is a system in which flow is on the whole one dimen­ sional, but in which absolute permeability, porosity, cross-sectional area, dip angle, and flow rate are not uniform along the axis of the flow. The equations of fluid flow are obtained by application of the law of mass conservation and Darcy's law. In particular, for water vw = —X^ldpjdx + pwg sin 0}

(144)

where Aw = kkryv/pw is the mobility of water. Production of fluids is assumed along the mean streamline (Fig. 42) in a continuous manner and the rate at which fluids are produced is proportional to the length of the interval dx along which fluids are withdrawn. The oil production linear distribution is q0(x) and has the dimensions of volume per unit length and time. 180

Two-Phase Flows in Porous Media

FIG. 42. Quasi-linear flow geometry. X is curvilinear abscissa in mean direction of flow and Θ is the dip angle (algebraic).

The oil volume influx into the streamtube through cross-section A(x) per unit of time is A(x)v0(x). It is positive if oil is effectively flowing in the direc­ tion of increasing x. Equating net oil influx to oil production and oil accumu­ lation, one obtains -d{A(x)v0(x,

t)}/dx = q0(x, t) + Α(χ)φ(χ) dSQ(x, t)/dt

(145)

A similar expression holds for water, -d(Av w )ldx = Αφ dSJdt + ? w .

(146)

For convenience, the following notation is introduced: Q = Av

(147)

Assuming v different from 0, one obtains the expression for the fractional flow function K =/w{l + Kkn A/Qßo[dpJdx

- Apg sin 0]}

(148)

with the usual notations. Neglecting capillary pressure, the expression for the fractional flow function becomes „( kA sin 0) F* =/w(l - Ko &P9 —Q—j

(149>

Substitution of vw in term of F w into Eq. (146) yields an equation for water saturation AA

dS

™ ^ ηπ ' ^ ot ox

r d® ox 181

_L. I

f

A

^

sin 0) ax

H. J. MOREL-SEYTOUX

where A is the cross-sectional area of flow, φ the porosity, Sw the water saturation, t the time, Q the total flow rate of the two phases, FJ the deriva­ tive with respect to 5W of the fractional flow function for water including the gravity term, x the distance along the mean direction of flow (see Fig. 42), #w the water production flow rate per unit length,/ w the fractional flow function for water without the gravity term, λτο the relative mobility to oil, Δρ the difference of the specific masses of water and oil, g the acceleration of gravity, k the absolute permeability, and Θ the dip angle (Fig. 42). Equation (150) is a first-order partial differential equation for 5W which must be solved to find Sw as a function of time and distance along a quasilinear porous medium. The equation reduces to the Buckley-Leverett equa­ tion (Morel-Seytoux, 1969) with appropriate simplifications. For simplicity, this equation is rewritten as P(x) dSJdt + R(x, S w , 0 dSJdx = T(x, S w , t)

(151)

The exact expressions for P(x), R(x, S w , t) and T(x, 5 W , t) can be obtained by comparison with Eq. (150). The solution of Eq. (151) can be obtained by reducing it to two ordinary differential equations: dt/P(x) = dx/R(x9 S w , 0 = dSJT(x,

S w , t)

(152)

The method of reducing partial differential equations to systems of ordinary differential equations is known as the method of characteristics. The method was used by J. C. Martin (1960) to solve a variety of problems fitting the description of quasi-linearity. First we shall concentrate upon another approach less mathematical in nature and closer to the Buckley-Leverett method of solution.

C. DlSCRETIZED QUASI-LINEAR MODEL

The solution of the general equation in saturation for quasi-linear flow P(x) dSJdt + R(x, 5 W , 0 dSJdx - T(x9 S w , t) = 0

(153)

cannot be obtained by the procedure followed to obtain the saturation profile in a homogeneous core with injection and production at the ends. The func­ tional relationship between S, x, and t can be written as Sw = Sw(x, t), or in a differential form as (dSJdt) dt + (dSJdx) dx - dSw = 0.

(154)

5 w = Sw(x, t) can be considered as the equation of a surface in the threedimensional system of coordinates x, t, S. On the surface Sw = 5w(x, t) an 182

Two-Phase Flows in Porous Media infinity of curves can be drawn, all the curves which satisfy an arbitrary additional differential equation. Contour lines, for instance, are curves on the surface that satisfies the differential equation dS = 0. However, due to the inhomogeneous character of Eq. (153), elimination of dSJdx and dSJdt between Eqs. (153) and (154) with dS = 0 to obtain the differential equation relating x and t is no longer possible. The method of characteristics involves the solution of two ordinary differ­ ential equations. The Buckley-Leverett method is consistent with this method even though it appears to involve only one ordinary differential equation. In fact, the second equation is trivial, namely dS = 0. Therefore, the BuckleyLeverett type of simplification follows immediately whenever Eq. (153) is a homogeneous equation, i.e., T(x, S, t) = 0. If T(x, S, t) = 0, then <7w + / w dQ/dx - 4>/ w ^99 d{kA sin 6)/dx = 0

(155)

Solving Eq. (155) gives (kA sin 0)/Q = C(t)

(156)

Equation (156) implies that the partial differential equation (Eq. (150)) for saturation will be homogeneous over any region in which the quantity (kA sin 9)1 Q is uniform. In other words, the equation will be homogeneous whenever the injection and production schedule and the geometrical charac­ teristics of the porous medium are so matched as not to alter the ratio of pressure forces along the main direction of flow. Therefore, if (kA sin 6)/Q is constant, the Buckley-Leverett type of calculation using fractional flow curves can be carried out. Usually this is not the case in reservoir problems of physical significance, and the fractional flow concept loses most of its usefulness because F w becomes a function not only of saturation but also of x and t. Thus we see that quasi-linear flow cannot be treated by the normal BuckleyLeverett method if the properties of the medium are represented as con­ tinuous functions of position. However, if the formation can be reasonably well represented by a series of blocks of uniform characteristics, a generalized Buckley-Leverett method using fractional flow curves can be used. If a reservoir is divided into a number of blocks of uniform properties, as indi­ cated in Fig. 43, and the quantity (kA sin 6)/Q is constant in each block, it follows that within each block a single fractional flow curve is applicable and a description of the path of a given saturation can be obtained by techniques similar to that of Buckley and Leverett. Therefore in order to use a general­ ized Buckley-Leverett method, the principal problem is to investigate what happens at boundaries between blocks. 183

H. J. MOREL-SEYTOUX

WO I

10 i

20 i

30 i

40 i

50 i

0 G I

FIG. 43. Schematic block approximation of reservoir. D. TYPICAL RESERVOIR CONDITIONS AND OPERATIONS

1. Variation in Reservoir Properties The fractional flow function Fis useful if it is only a function of saturation, that is if the quantity (kA sin 6)/Q is constant. However, the parameter (kA sin θ)/ζ) will usually vary because properties of the reservoir vary with space and time. Variations in flow rate Q will be discussed in the following sections. For the moment, our attention will be focused on the influence of variations in reservoir properties k, A, or Θ on the evolution of saturation profiles under steady-state flow conditions. To handle this problem, we assume the reservoir is either naturally divided into zones of uniform proper­ ties with sharp transition at boundaries or artificially segmented, for calcula­ tion purposes, into blocks of uniform characteristics (Fig. 43). Then, within each block the velocity at which a saturation travels is obtained from the slope of the tangent to the fractional flow curve and is uniquely and well defined. However, at the boundary the situation is unclear. Due to the discontinuity in reservoir properties, a resulting discontinuity in saturation can be intuitively expected to occur at the boundary. Once created the discontinuity could either be stationary or propagate into the adjacent blocks. Whatever occurs at the 184

Two-Phase Flows in Porous Media boundary must satisfy the law of conservation of mass and be compatible with the information deduced from the fractional flow curves on either side of the boundary. The following typical case is presented to show what may happen at a boundary. a. ENCROACHING UPDIP WATER DRIVE. Consider the case of an updip advance of a water-oil contact. This is both a simple and common incidence and we shall use this particular case as an introduction to the general approach. At some distance updip from the original water-oil contact the properties of the reservoir are assumed to change abruptly. The block in contact with the aquifer is denoted block 1. The next block updip is block 2 and therefore water moves through block 1 toward block 2. For simplicity, only two blocks are considered. At the boundary the properties of the reservoir change so that the quantity kA sin Θ has a different value on each side of the boundary. It is immaterial whether this variation is due to a change in absolute permeability, cross-sectional area, dip angle, or a combination of these because the influence of formation characteristics is impressed on the fractional flow curve through the lump parameter kA sin 0. For purposes of this example, a change in dip angle at the boundary is assumed to be such that the angle remains positive. First consider the case of an increase in kA sin θ at the crossing of the boundary in the direction of flow. That is (kA sin 0)x < (kA sin θ)2

(157)

Subscripts or superscripts 1 and 2 always refer to blocks 1 and 2. Because the density difference in a water-oil system is positive, it follows from the expres­ sion of F w that for any given saturation S /w(S)>iV(S)>iV2(S)

(158)

as shown in Fig. 44. Subscript w, referring to water, will often be omitted. Unless specifically stated otherwise, F and S will be understood to mean F w and S w . Because block 1 is in complete contact with the aquifer the saturation profile in block 1 behaves as a normal Buckley-Leverett profile with a velocity that is obtained from the proper values of F^1. When the BuckleyLeverett front for block 1 reaches the boundary between blocks 1 and 2 water starts moving into block 2. However, the saturation profile in block 2 will not be a normal Buckley-Leverett profile because only water saturations between 5WC and 5 a , shown in Fig. 44, are available to move into block 2. As the water drive displaces more and more oil from block 1 to block 2, the saturation of water on the left side of the boundary will gradually increase from 5 a to Sh, to Sc, and so forth. At the instant that any value of saturation, Sa or Sh, reaches the boundary on the left side we assume that the saturation 185

H. J. MOREL-SEYTOUX Variation in Reservoir IS

Properties

Boundary

Aquifer

Water Encroachment

^

Buckley - Leverett I Ic A Λ

Front-^j

Unstable

-j

I

Fronts

I

ΊΪ

BuckleyLeverett Front o

Evolution of Saturation Profile

Fractional Flow Curves

Abrupt Change in kA sin Θ

BL, FIG. 44.

BU

Effect of variation in reservoir properties on the saturation profile.

in block 2 is not Sa or Sh but S a ' and S b '. Therefore the saturation profile undergoes a discontinuity at the boundary. The size of the discontinuity, including zero, will be obtained by considering material balance at the boundary. Because there is neither injection or production at the boundary, the total flow rates must be the same on both sides, or Öl = Ö2 = Q

(159)

This equation does not mean that the flow rates of oil and water must be conserved on both sides. It only says that the sum of the flow rates of water and oil is conserved at the crossing of the boundary: Öw 1 + ß o 1 = Ö w 2 + Ö o 2 186

(160)

Two-Phase Flows in Porous Media or using the fractional flow function and dividing through by ß Fj(S)

+ F0\S) = FW2(S') 4- F02(S')

(161)

Because F0 = 1 — Fw, Eq. (161) is an identity no matter what 5 and S' may be, therefore Eq. (161) alone cannot give the value of the saturation dis­ continuity. Because water flowing in block 2 comes from block 1, the flow rate in block 2 cannot exceed that in block 1. Symbolically, ßw2 < ßw1

(162)

or FW2(S') < FW\S)

(163)

2

FWA(S)

First assuming the case where ß w * = ß w , then it follows that = FW2(S') and consequently the jump from 5 to S' is given by the equation FJ(S) = FW2(S')

(164)

The solution can be obtained graphically using the fractional flow diagrams, by drawing horizontal lines such as aa\ bb\ and cc' to obtain values of Given a saturation profile at an arbitrary initial time t0 (profile oabc, Fig. 44) the saturation profile at a later time t is easily constructed. The rate of advance of a given saturation (e.g., Sh) in block 1 is obtained in the usual manner. If over the time interval (t — t0) this saturation does not reach the boundary of blocks 1 and 2, its location at time t is easily obtained as OF* '(Sh) x(Sh, 0 = * ! \ h ) (t - t0) + x0(Sh, io)· ^ιΦι

(165)

If on the contrary the saturation 5 b reaches the boundary at a time t' less than f, one must proceed to calculate the saturation jump. The saturation Sh' is the solution of the equation ^ ( S O = F2(Sb')

(166)

a solution which is easily obtained graphically. Then the location of Sh' at time t is obtained by the formula x(S'h, 0 = ^

^

(f - t0) +

QF

/^h)

(t - O + x0(Sh, *o)·

(167)

Proper applications of Eqs. (165) or (167) for all values of saturation yield the complete saturation profile at time t. 187

H. J. MOREL-SEYTOUX

The displacement process in block 2 is more efficient than in block 1. This result is due to the fact that (kA sin 0)2 > (kA sin θ)1 which means physically that in block 2 the gravity force plays a relatively more important part in the flow of the fluids. In block 2 buoyancy makes the displacement more efficient, speeding up the flow of oil, slowing down the flow of water. The greater efficiency of the displacement process in block 2 should not be misunderstood to mean that the size of the Buckley-Leverett front in block 2 is necessarily greater than that in block 1, because it is not always the case. Figure 45 illustrates a case for which the size of the Buckley-Leverett front in region 2 is smaller than in block 1. The evolution with time of the saturation profile is also shown for this case. The essential difference beteeen the cases of Figs. 44 and 45 is that in one case the front immediately assumes the Buckley-Leverett value S B L 2 , Encroaching Updip Water Drive Si

Boundary

Aquifer

Encroaching Water Q

Evolution of Saturation Profile

BL2 BL,

FIG. 45.

Evolution of saturation profile in an encroaching updip water drive. 188

Two-Phase Flows in Porous Media whereas in the other case its size grows and its motion accelerates with time. A variety of situations may occur depending on the medium properties, flow rates, etc. For a further discussion particularly with counterflow and partial transmission of water flux see (Morel-Seytoux, 1967). b. MOBILE INITIAL WATER OF GAS SATURATION.

In the earlier examples, it

was tacitly assumed, e.g., in the case of a water drive, that the initial water saturation was immobile. When 5 w i = S wr the saturation profile in block 2 is not affected until the water drive front has reached the boundary. If the initial water saturation is actually greater than S w r , the presence of the boundary immediately affects the saturation profile in block 2. Figure 46 illustrates the evolution of the saturation profile with time. As soon as water moves from the aquifer into block 1, a second front is formed and propagates in block 2 starting at the boundary. Its size is determined by the horizontal Evolution of Soturation Profile Boundary

i _ _ BJ0Ck__?·

Block l _

S

wl

13

ίιοπο«ιχ_~ I I 1-

Influence of Mobility of Iniliol Soturotion

BL,

BL 2 BL 2

FIG. 46. Influence of mobility of initial saturation on the saturation profiles. 189

H. J. MOREL-SEYTOUX

line construction OtO' and the velocity with which it moves in block 2 is proportional to the slope of the secant 02 O'. When the Buckley-Leverett front for block 1 reaches the boundary the front 02 O' has established a new initial saturation S0' in block 2. The front in block 2 in this instance is at first unstable and with time grows to the Buckley-Leverett value, BL 2 ', corre­ sponding to the initial saturation S0'. But this Buckley-Leverett front moves more rapidly than the front 02 O' and will catch up with it. When it does so and after a short period of time during which the new front is unstable, a new Buckley-Leverett front BL2 is established that corresponds to the initial saturation Swi. c. INTERACTION AND REFLECTION OF FRONTS. A few examples of fronts running into one another were already encountered in previous sections. In the case of an abrupt change in cross-sectional area with mobile initial water saturation (Fig. 46) the Buckley-Leverett front corresponding to the initial Interoction

of Front

Encroaching Water

—X

Fraction Flow Function

J » BL

B L I B t f o r · Cloth)

(After Oath)

FIG. 47. Effect of interaction of fronts on saturation profiles. 190

Two-Phase Flows in Porous Media saturation S0' in block 2 caught up with the front 02 O'. When two fronts run into one another, a resulting front is created immediately. Before colliding the two fronts had a common saturation. The size of the new front is the differ­ ence between the two saturations of the fronts that were not common to both fronts. In the case of Fig. 47, the fronts Ob and Oi collide. They had the saturation S0 in common. The resulting front at the time of the clash is front ib. This front is not necessarily stable. In previous sections unstable fronts were shown to grow with time. Here is a case where the front immediately decays into a Buckley-Leverett front ia and a continuous profile ab. The front ib is also mathematically acceptable, but so is the solution: front ia, continuous profile ab, which is the solution that requires a minimum of dis­ continuity compatible with singlevaluedness, and is therefore the solution. The case of reflection of fronts on a no-flow boundary is little different from the interaction of two fronts. Consider the downdip moving front oa of Fig. 48. As soon as it reaches the no-flow boundary the water saturation Cellar Oil : Reflection

at

Boundary

Q = -Q2 Water Injection

No-flow Boundary

Fractional Flow Curves

x

y

—^^-

*

^ · Minimum of Fw

FIG. 48. Reflection at a boundary. 191

H. J. MOREL-SEYTOUX

establishes itself to (1 — 5or), point n. The new front is the resultant combina­ tion of front oa and immobile front on, i.e., an. That front decays immediately into the Buckley-Leverett front am and the continuous profile mn. 2. Injection or Production under Unsteady-State Conditions In the central portion of the reservoir, water is injected into an updip moving stream of oil. Conditions are assumed such that the upstream flow rate is large enough to prevent downdip flow of water. Thus, at some time t after injection was initiated the saturation profile will be as indicated in Fig. 49. Up to that time the injection rate was AQW and the flow rate in block 2 was 02 = öi ~ Δ β , (168)

Water Injection

BL

FIG. 49.

BL

Influence of changing injection flow rate on saturation profiles. 192

Two-Phase Flows in Porous Media The water saturation at the boundary Sa is obtained as the solution of F2(Sa)=-AQJQ2

(169)

At time t the injection rate is lowered to the value |Ag w *| and this new injection rate is maintained thereon. A new flow rate Q2* prevails in block 2 Ö2* = Qi~ Δβ„*

(170)

and the water saturation at the boundary drops to the value Sb, solution of F2*(Sb)=-AQv/*/Q2*

(171)

The saturation profile at the time t* > t (see Fig. 49) consists of an imbibi­ tion profile (ocd) and a drainage tail which in this instance turns out to be a drainage front (ab). After a number of fluctuations in injection flow rate the saturation profile will picture an alternance of imbibition profiles and drainage profiles, mainly drainage fronts.

VII. Approximate Treatment of Three-Phase Flow

A. WATER-OIL-GAS SYSTEMS

The displacement of gas by oil is an efficient process because the viscosity ratio is very favorable. It is even more efficient when the general direction of flow is updip because the density difference is large. A typical oil-gas fractional flow curve is shown in Fig. 50. It is hard to distinguish between the Welge construction of the tangent to the curve from o and the secant oe. The displacement is practically pistonlike. It is therefore a good approximation to replace the front oa and the continuous profile ae by the front oe. In most instances the flow of oil is itself due to a water drive, natural or artificial. Thus water moves into a gas-oil zone and sets the oil and gas in motion. The oil in turn displaces the gas in a pistonlike manner and an oil bank is created between the water-oil front and the oil-gas front. The satura­ tion profile at a given time is illustrated in Fig. 51. 193

H. J. MOREL-SEYTOUX

| sgc

I

""' ~1°

s

*

0

Oil

swc

FIG. 50. Practically pistonlike displacement (Dougherty and Sheldon, 1964).

The less depleted the reservoir, i.e., the larger Soi the faster the oil-gas front propagates. It must be noted that the velocity of such a front is

Vt

/ 1 ~/.t \ v

\sgi - sj Φ

(172)

If the initial oil saturation is large, foi is not small compared to 1 and should not be neglected in Eq. (172). For small initial oil saturations, the front may move rather slowly. It may even move more slowly than the water-oil front, in which case the saturation profile of Fig. 51 does not apply [Dougherty and Sheldon, 1964]. In order not to have any three-phase flow region another approximation is necessary. It consists of choosing a solution for the water drive that uses more than the 194

Two-Phase Flows in Porous Media minimum amount of discontinuity necessary for single-valuedness. Material balance is satisfied but the solution is nevertheless approximate (Fig. 51) [Dougherty and Sheldon, 1964].

(a)

(b)

FIG. 51. Water drive-gas phase present (Dougherty and Sheldon, 1964).

B. STEAMFLOOD DISPLACEMENT PROCESS

As an example of applicability of the method, the displacement process in a laboratory steamflood is now briefly discussed. Detailed account of the procedure has been given by Pritchard (1966). On the basis of laboratory steamfloods (Willman et αί, 1961) it can be concluded (Figs. 52, 53, 54) that in the process, approximately: (1) the residual oil saturation is reduced; (2) a sharp temperature front exists in the core; (3) the steam front and the temperature front coincide at all times; and (4) very little oil flows after steam breakthrough. The following model of the steamflooding displacement process can then be accepted as a first approxi­ mation. In the steam-invaded zone the displacement is pistonlike. Behind the steam front the saturations are uniform and equal to (Sor)s for oil, S sr for steam and (5 wr ) s for water. The velocity at which the steam front propagates is denoted by Vf. Coming out of the steam front is a mixture of water and oil made of displaced oil, of displaced water, and of water condensed from steam. The 195

H. J. MOREL-SEYTOUX

2 4 6 8 10 Total Produced Fluids—Pore Volumes (Steam Considered as Equivalent Volume Condensate)

FIG. 52. Oil recoveries by steam injection, hot water flood, and cold water flood of crude oil (Willman et a/., 1961).

\ (Waterfloods) f—Water Breakthrough I (Steam Injection)

0

2

4

6

8

10

Total Produced Fluid-Pore Volumes (Steam Considered as Equivalent Volume of Condensate)

FIG. 53. Oil recoveries by steam injection, hot water flood, and cold water flood of nondistillable white oil (Willman et al., 1961).

rate of increase of volume occupied by steam is Kf^40Ssr. The injection rate of steam is Qs. The difference of these two rates gives the rate of condensation of steam. The flow rate of water condensed from steam coming out of the steam face is therefrom: (Qs — V{ Αφ8&τ)ρ5/ρν/. The rate of displacement of liquids, oil, and water, is equal to the rate of steam invasion, i.e., ν{Αφ8ΒΤ. Let Sw~ and SQ~ be the saturations existing at the steam front downstream face. Then the rate of displacement of water is v{A[SW- - (Swr)J = vtA4>[S„- - 1 + (Sor)s + Ssr] 196

(173)

Two-Phase Flows in Porous Media 400 \

\K

300 LX~

S 200

> ^

I 100

>]

b w

\

\ \ ,

\

OK

o^^

8 12 16 20 24 28 32 Distance from Core Inlet-inches Core 9 , 2 5 % Distillable 0il,AP=40psi

36

Legend 50 minutes 100 minutes 150 minutes 200 minutes 300minutes FIG. 54.

°

Temperature profiles during steam injection (Willman et al.9 1961).

Finally the total flow rate out of the steam face is Q = (Qs - v{ASsr =QA+

vfA4>S„(l - ^)

(174)

and the fraction of it that is water is (Q, - v, AtfrSJ H + νίΑφ[Ξν- + (Sor)s + Ssr - 1] Pw

/w =

Q (175)

Ö

or /w(s.-) = i

Αφ[\ - S w " ρ ^

+ ΑφΞ, r

\

(S„)M0

pJ

(176)

vt(0

With S sr , (Sor)s, and Vf(t) being known from experiments, the above equation takes the form /W(SW~) = Λ(ί) + B(t)SK (177) The value of Sw~, water saturation at the steam front downstream face, can be obtained as the abscissa of the intersection of the curve / w and the straight line A + BS. The values and variations of Sw~ versus Vf and there­ from versus t are then known. The displacement process ahead of the steam front is of the same nature as that described for unsteady injection conditions. The drive ahead of the steam front is one by a mixture of water and oil whose composition changes with time. 197

H. J. MOREL-SEYTOUX

In the laboratory floods the steam front did not propagate immediately. In the case v{ = 0 formula (176) yields / w = 1, which means that in the initial phase of the steamflood a standard cold water drive took place. In the second phase as the steam front moves along the core, the drive is of a mixed type and the saturation profile consists of an imbibition zone followed by a drainage front. A detailed account of the evolution of the saturation with time has been given by Pritchard (1966) in a report from which Figs. 55 and 56 are borrowed. Comparison of calculated values and experimental data for a sensitive quantity such as water-oil ratio is good and substantiates the validity of the model for the displacement process. To make recovery predictions in a reservoir the rate of advance of the steam front must be predicted. Based on the assumption that heat transfer occurs within the reservoir by convection only, that heat is lost to the over­ burden and underburden by conduction and that the temperature profile is a front, a heat balance (Marx and Langenheim, 1959) provides the answer vM = (C/Sn)Qsdrtt)ldt (178) for a linear case, where δη is the width, Qs is the steam injection rate, η(ί) is given by Eq. (179), and C is an appropriate constant [Marx and Langenheim, 1959] depending on reservoir thickness, fluid and rock heat capacities, etc., and η(ΐ) = ea2t erfc (at112) + 2at1/2n~1/2 - 1 (179) where a is a constant (Marx and Langenheim, 1959) depending on the

Location of the

o,l υ

*Ό

I

10

I

20

I

30

I

40

1

50

Distance, cm

1

60

Production Face 90.7 cm

1

70

1

80

"~*Ί

II I

1

90

1

100

FIG. 55. Saturation profiles from commencement of propagation of the steam front to breakthrough of the oil bank (Pritchard, 1966). 198

Two-Phase Flows in Porous Media 10.0

—

8.0 o σ (z O

6.0

a> cj σ»

E 4.0

o •o o

a. 2.0 1.0 0

20

40

60

80

100

120

Time , minutes FIG. 56. Water-oil rate versus time during steam flood (Pritchard, 1966).

reservoir and fluid characteristics. Knowing explicitly the velocity of the temperature front as a function of time the coefficients A and B in Eq. (177) are also known explicitly. Intersects of the straight line of time-dependent equation with the time-independent/ w curve yield pair of values t and 5W~ that define the function Sw~ = S w ~(0· It is a simple matter now to predict the time evolution of the saturation profile and to calculate the performances of the flood in a linear or radial geometry. V m . Two-Dimensional Two-Phase Flows In this case the problem of obtaining analytical solutions is so difficult that the numerical method of the finite-difference approximation has been almost exclusively used. The method has been used extensively in the petroleum industry but little in the field of hydrology. Apparently the only exception has been the work of Brutsaert (1970) who studied the well flow problem in an infinite aquifer from a two-phase point of view. ACKNOWLEDGMENTS

The author wants to thank the former students who worked their thesis or dissertation under his supervision for their very valuable research and results: Mr. Le Van Phuc, Republic of South Vietnam, Dr. R. L. Brustkern, Research Associate, Montana State University, and Mr. A. Noblanc, Engineer, Center for Oceanographic Research of Britany. The author had the pleasure to serve on the Ph.D. committee of Dr. D. McWhorter, Assistant Professor of Agricultural Engineering, Colorado State University and of Dr. W. Brutsaert, Associate Professor of Civil Engineering, New Mexico Institute of Mining and Technology. Exchange of ideas with all these excellent former students were very fruitful. 199

H. J. MOREL-SEYTOUX During 1970-1971 an interdisciplinary seminar brought together as good a team as one can gather anywhere, I believe, to discuss water movements in soil: Dr. A. Klute, Professor of Soil Science, Department of Agronomy, Dr. A. T. Corey, Professor of Agricultural Engineering, Dr. D. McWhorter, Assistant Professor of Agricultural Engineering, Dr. H. Gardner, and Dr. D. Heerman, Research Scientists, Agricultural Research Service and several students, notably Mr. P. Hamaker and Mr. R. Gilman from the Department of Agronomy and Mr. Sonu from the Department of Civil Engineering, Colorado State University. To all these colleagues I want to express my sincere thanks for a stimulating series of seminars. I want to extend my appreciation to Mr. J. Sonu, Graduate Research Assistant currently working on his Ph.D., for helping me diligently in shaping this chapter by preparing figures, performing various calculations, and proofreading the manuscript. I want to thank Mrs. M. Fox, of the Engineering Research Center for her excellent typing services and Miss H. Akari for the illustrations. I want to thank Chevron Oil Field Research Company and my former La Habra officemate, K. C. G. Pritchard, for permission to use material from unpublished research reports. The research on which much of this article is based was supported by the Department of the Interior, Office of Water Resources Research, as authorized under the Water Resources Research Act of 1964, and pursuant to Grant Agreement No. 14-01-0001-1886. The OWRR support is gratefully acknowledged. Finally I want to thank Dr. J. W. N. Fead, Head, Department of Civil Engineering, Colorado State University for his moral and financial support in this endeavor. REFERENCES

ADRIAN, D. D. (1964). " T h e Influence of Capillarity on the Flow of Air and Water in a Porous Medium,'* Tech. Rep. No. 38. Dep. of Civil Eng., Stanford Univ., Stanford, California. ADRIAN, D. D., and FRANZINI, J. B. (1966). Impedance to infiltration by pressure build-up ahead of the wetting front. / . Geophys. Res. 71, 5857-5863. AMES, W. F. (1964). "Nonlinear Partial Differential Equations in Engineering," 511 pp. Academic Press, New York. BEAR, J. (1972). "Dynamics of Fluids in Porous Media." Elsevier, Amsterdam. BRUSTKERN, R. L. (1970). "An Analytical Treatment of Two-Phase Flow during Infiltra­ tion," 105 pp. Ph.D. Thesis, Dep. of Civil Eng., Colorado State Univ., Fort Collins, Colorado. BRUSTKERN, R. L., and MOREL-SEYTOUX, H. J. (1970). Analytical treatment of two-phase infiltration. / . Hydraul. Div., Proc. Amer. Soc. Civil Eng. 96(Hy 12), 2535-2548. BRUSTKERN, R. L., and MOREL-SEYTOUX, H. J. (1972). Description of water and air move­ ment during infiltration. Water Resour. Res. (submitted). BRUTSAERT, W. (1970). "Immiscible Multiphase Flow in Ground Water Hydrology: A Computer Analysis of the Well Flow Problem," 119 pp. Ph.D. Thesis, Dep. of Civil Eng., Colorado State Univ., Fort Collins, Colorado. BUCKLEY, S. E., and LEVERETT, M. C. (1942). Mechanism of fluid displacements in sands. Trans. AIME 146, 107-116. CALHOUN, J. C. (1953). "Fundamentals of Reservoir Engineering." Univ. of Oklahoma Press, Norman, Oklahoma. CHILDS, E. C. (1969). "An Introduction to the Physical Basis of Soil Water Phenomenon." Wiley, New York. COREY, A. T. (1957). Measurement of water and air permeability in unsaturated soil. Soil Sei. Soc. Amer., Proc. 21(1), 7-10. 200

Two-Phase Flows in Porous Media CRAFT, B. C , and HAWKINS, M. F. (1964). "Applied Petroleum Reservoir Engineering." Prentice-Hall, Englewood Cliffs, New Jersey. DARCY, H. (1856). " L e s Fontaines Publiques de la Ville de Dijon." Dalmont, Paris. D E WIEST, R. J. M. (1967). " Geohydrology," 366 pp. Wiley, New York. D E WIEST, R. J. M., ed. (1969). " F l o w Through Porous Media," 530 pp. Academic Press, New York. DIXON, R. M., and PETERSON, A. E. (1971). Water infiltration control: A channel system concept. Soil Sei. Soc. Amer., Proc. 35, 968-973. DOUGHERTY, E. L., and SHELDON, J. W. (1964). The use of fluid interfaces to predict the behavior of oil recovery process. Soc. Petrol. Eng. J. 4(2), 171-182. ELRICK, D . E. (1961). Transient two-phase capillary flow in porous media. J. Phys. Fluids 4(5), 572-575. FITCH, R. A., and GRIFFIN, J. D . (1964). Experimental and calculated performance of miscible floods in stratified reservoirs. J. Petrol. Technol. 16(1), 1289-1298. FREE, J. R., and PALMER, V. J. (1940). Relationship of infiltration, air movement and pore size in graded silica sand. Soil Sei. Soc. Amer., Proc. 5, 390-398. GARNER, D . E., DONALDSON, J. K., and TAYLOR, G. S. (1969). Entrapped soil air in a field site. Soil Sei. Soc. Amer., Proc. 33, 634-635. GAUDET, J. P. (1972). "Effet de l'ecoulement de l'air sur une infiltration d'eau dans une colonne de sol verticale." Rapport de D.E.A., Universite Scientifique et Medicale de Grenoble. GOOLSBY, J. L., and ANDERSON, R. C. (1964). Pilot water flooding in a dolomite reservoir, the McElroy Field. / . Petrol. Technol. 16(12), 1345-1350. GREEN, D . W., DABIRI, H., WEINAUG, C. F . , and PRILL, R. (1970). Numerical modeling of

unsaturated ground water flow and comparison of the model to a field experiment. Water Resour. Res. 6(3), 862-874. HANDY, L. (1960). Determination of effective capillary pressures for porous media from imbibition data. Trans. AIME 219, 75-80. HILLEL, D . (1971). "Soil and Water. Physical Principles and Processes," 288 pp. Academic Press, New York. HOLM, L. W., CSASZAR, A. K., and BERNARD, G. G. (1965). Field test shows STP advantage in waterflood. Oil Gas J. 63, 72-76. HORTON, R. E. (1940). An approach toward a physical interpretation of infiltration capacity. Soil Sei. Soc. Amer., Proc. 5, 399-417. MCWHORTER, D . B. (1970). "Infiltration Affected by Flow of Air," 129 pp. Ph.D. Thesis, Dep. of Agr. Eng., Colorado State Univ., Fort Collins, Colorado. MCWHORTER, D . B. (1971). "Infiltration Affected by Flow of Air," 43 pp. Hydrol. Pap. No. 49. Colorado State Univ., Fort Collins, Colorado. MCWHORTER, D . B., and COREY, A. T. (1967). Similitude for flow of two fluids in porous media. Int. Hydrol. Symp., Fort Collins, Colo. pp. 136-140. MARTIN, J. C. (1960). Theory of one-dimensional two-phase flows of incompressible fluids. Prod. Mon. 24(9), 18-29. MARX, J. W., and LANGENHEIM, R. H. (1959). Reservoir heating by hot fluid injection. Trans. AIME 216, 312-315. MOREL-SEYTOUX, H. J. (1966). " F l o w of Immiscible Fluids in Porous Media," 433 pp. Tech. Memo. Chevron Res. Co., La Habra, California. MOREL-SEYTOUX, H. J. (1967). A study of quasi-linear noncapillary two-phase flow in porous media. Develop. Mech. 4, 1321-1335. MOREL-SEYTOUX, H. J. (1969). Introduction to flow of immiscible liquids in porous media. In " F l o w Through Porous M e d i a " (R. De Wiest, ed.), pp. 455-516. Academic Press, New York. 201

H. J. MOREL-SEYTOUX MOREL-SEYTOUX, H. J. (1971). " A Systematic Treatment of the Problem of Infiltration," OWRR Completion Rep. Ser. No. 23. Environ. Resour. Cent., Colorado State Univ., Fort Collins, Colorado. MOREL-SEYTOUX, H. J. (1972). A new analytical treatment for the infiltration problem. Proc. Joint Symp. Fundam. Transp. Phenomena Porous Media, 2nd, Guelph, Ont., pp. 257-265. Office of Continuing Education, University of Guelph, Canada. MOREL-SEYTOUX, H. J., and NOBLANC, A. (1971). Infiltration predictions by a moving strained coordinate method. Proc. Symp. Soil-Water Phys. Technol., Rehovot, Isr. MUSKAT, M. (1949). "Physical Principles of Oil Production." McGraw-Hill, New York. NOBLANC, A. (1970). " A Capillary Perturbation Study of Two-Phase Flow Infiltration." M.S. Thesis, Dep. of Civil Eng., Colorado State Univ., Fort Collins, Colorado. NOBLANC, A., and MOREL-SEYTOUX, H. J. (1972). Perturbation analysis of two-phase infiltration. ASCE / . Hydraul. Div., Proc. Amer. Soc. Civil Eng. 98 (HY9), 1527. PECK, A. J. (1965a). Moisture profile development and air compression during water uptake by bounded porous bodies: 2. Horizontal columns. Soil Sei. 99, 327-334. PECK, A. J. (1965b). Moisture profile development and air compression during water uptake by bounded porous bodies: 3. Vertical columns. Soil Sei. 100, 44-51. PHILIP, J. R. (1957). The theory of infiltration: 1. The infiltration equation and its solution. Soil Sei. J. 84, 163-178. PHILIP, J. R. (1969). Theory of infiltration. Advan. Hydrosci. 5, 215-305 PHUC, LE VAN (1969). "General One-Dimensional Model for Infiltration," 83 pp. M.S. Thesis, Dep. of Civil Eng., Colorado State Univ., Fort Collins, Colorado. PHUC, LE VAN, and MOREL-SEYTOUX, H. J. (1972). Effect of soil air movement and com­ pressibility on infiltration rates. Soil Sei. Soc. Amer., Proc. 36, 237-241. PIRSON, S. J. (1958). " O i l Reservoir Engineering," 735 pp. McGraw-Hill, New York. PRITCHARD, K. C. G. (1966). " A Comparison of Steam Flood Prediction Results with Experimental Data." Chevron Res. Co., La Habra, California. RICHARDS, L. A. (1931). Capillary conduction of liquids in porous media. Physics {New York) 1,318-333. SCHEIDEGGER, A. E. (1960). " T h e Physics of Flow Through Porous Media," 2nd Ed. Univ. of Toronto Press, Toronto. SMITH, C. R. (1966)." Mechanics of Secondary Oil Recovery," 154 pp. Reinhold, New York. SONU, J. (1972). Unpublished results, Ph.D. Thesis, dissertation, Dep. of Civil Eng., Colorado State Univ., Fort Collins, Colorado. SWARTZENDRUBER, D . (1969). The flow of water in unsaturated soils. In " F l o w Through Porous Media" (R. De Wiest, ed.), pp. 215-292. Academic Press, New York. VACHAUD, G., VAUCLIN, M., WAKIL, M., and KHANJI, D . (1971). Infiltration ä debit

constant dans une colonne verticale de sol stratifie. Influence de la pression de l'air. Congr. Int. Ass. Hydraul. Res., 14th 5, 439-448. VACHAUD, G., VAUCLIN, M., WAKIL, M., and KHANJI, D. (1972). Effects of air pressure

on water flow in an unsaturated, stratified, vertical column of sand. J. Water Resour. Res. 9 (1). To be published. WELGE, H. J. (1952). A simplified method for computing oil recovery by gas or water drive. Trans. AIME 195, 91-98. WILLMAN, B. T., WALLEROY, V. V., RUNBERG, G. W., CORNELIUS, A. J., and POWERS, L. W.

(1961). Laboratory studies of oil recovery by steam injection. J. Petrol. Technol. 13,681-689. WILSON, L. G., and LUTHIN, J. N . (1963). Effect of air-flow ahead of the wetting front on infiltration. Soil Sei. 96, 290-294. YOUNGS, E. G., and PECK, A. J. (1964). Moisture profile development and air compression during water uptake by bounded porous bodies: 1. Theoretical introduction. Soil Sei. 98,290-294. 202