EPA/600/2-91/015
May 1991
APPROXIMATE MULTIPHASE FLOW MODELING BY CHARACTERISTIC
METHODS
by
James W. Weaver
Processes and Systems Research Division
Robert S. Kerr Environmental Research Laboratory
United States Environmental Protection Agency
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
UNITED STATES ENVIRONMENTAL PROTECTION AGENCY
ADA, OKLAHOMA 74820
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completi
1. REPORT NO.
EPA/600/2-91/015
2.
PB91-190959
4. TITLE AND SUBTITLE
APPROXIMATE MULTIPHASE FLOW MODELING BY CHARACTERISTIC
METHODS
5. REPORT DATE
May 1991
&. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
James W. Weaver
8. PERFORMING ORGANIZATION REPORT NO
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Robert S. Kerr Environmental Research Laboratory
U.S. Environmental Protection Agency
P.O. Box 1198
Ma, Oklahoma 78420
10. PROGRAM ELEMENT NO.
ABWD1A
11 CONTRACT/GRANT NO
In-House
12. SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Laboratory
U.S. Environmental Protection Agency
P.O. Box 1198
Ada, Oklahoma 78420
13. TYPE OF REPORT AND PERIOD COVERED
Project Report/Summary
14. SPONSORING AGENCY CODE
EPA/600/15
15. SUPPLEMENTARY NOTES
Project Officer: James W. Weaver
FTS: 743-2420
16. ABSTRACT
The flow of petroleum hydrocarbons, organic solvents and other
liquids that are immiscible with water presents the nation with
some of the most difficult subsurface remediation problems. One
aspect of contaminant transport associated releases of such liquids;
is the transport as a water-immiscible liquid phase. In thj^T
document approximate models of immiscible flow are presented for
two- and three-phase flow. The approximations are constructed by
representing the flow by hyperbolic equations which have method of
characteristics solutions. This approximation has the additional
benefit of being based on the fundamental wave behavior of the
flow, which is revealed by the solutions of the models. An
important result is that for three-phase flow, two flow regimes
exist. The first is characterized by the displacement of one of
the liquids into a bank which moves ahead of the other liquid. The
-second is characterized by almost complete bypassing of a liquid by
the other. The occurrence of the flow regimes is dependent on the
organic liquid properties, soil type and the initial amounts of the
fluids present.
7.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS
COSATi field,Group
8. DISTRiaUTION STATEMENT
RELEASE TO THE PUBLIC
19 SECURITY CLASS iTIns Report)
UNCLASSIFIED
21 NO OF
128
20 SECURITY CLASS (Tins page'
UNCLASSIFIED
22 PRICE
EPA Form 2220-1 (R.». 4-77) PREVIOUS EDITION is OBSOLETE
-------
DISCLAIMER
The information in this document has been funded wholly or in part by the United States
Environmental Protection Agency. It has been subjected to the Agency's peer and administrative
review, and ft has been approved for publication as an EPA document Mention of trade names or
commercial products does not constitute endorsement or recommendation for use.
All research projects making conclusions or recommendations based on environmentally
related measurements and funded by the United States Environmental Protection Agency are required
to participate in the Agency Quality Assurance Program. This project did not involve environmentally
related measurements and did not involve a Quality Assurance Plan.
-------
FOREWORD
EPA is charged by Congress to protect the Nation's land, air and water systems. Under a
mandate of national environmental laws focused on air and water quality, solid waste management
and the control of toxic substances, pesticides, noise and radiation, the Agency strives to formulate
and implement actions which lead to a compatible balance between human activities and the ability of
natural systems to support and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the Agency's center of expertise for
investigation of the soil and subsurface environment. Personnel at the Laboratory are responsible for
management of research programs to: (a) determine the fate, transport and transformation rates of
pollutants in the soil, the unsaturated and the saturated zones of the subsurface environment; (b)
define the processes to be used in characterizing the soil and subsurface environment as a receptor of
pollutants; (c) develop techniques for predicting the effect of pollutants on ground water, soil, and
indigenous organisms; and (d) define and demonstrate the applicability and limitations of using
natural processes, Indigenous to the soil and subsurface environment, for the protection of this
resource.
This report describes a mathematical technique for simulating the flow of water, organic liquids
and air hi the unsaturated zone. A suite of assumptions are made in order to reduce the governing
equations to simple forms. From the solutions of the latter, the fundamental behavior of such systems
is revealed. Critical insight into the flow of fluids in the subsurface is gained from this work.
Clinton W. Hall
Director
Robert S. Kerr Environmental Research Laboratory
in
-------
ABSTRACT
The flow of petroleum hydrocarbons, organic solvents and other liquids that are immiscible
wfth water presents the nation with some of the most difficult subsurface remediation problems. One
aspect of contaminant transport associated with releases of such liquids is transport as a water-
Immiscible liquid phase. Conventional finite element and finite difference models of multiphase flow
may have extreme data and computational requirements. Sites with insufficient data or modeling for
screening purposes may warrant using a simplified approach. In this document approximate models
of immiscible flow are presented for two- and three-phase flow. The approximations are constructed
by representing the flow by hyperbolic equations which have method of characteristics solutions. This
approximation has the additional benefit of being based on the fundamental wave behavior of the flow,
which is revealed by the solutions of the models. An important result is that for three-phase flow, two
flow regimes exist The first is characterized by the displacement of one of the liquids into a bank
which moves ahead of another liquid. The second is characterized by almost complete bypassing of a
liquid by the other. The occurrence of the flow regimes is dependent on the organic liquid properties,
soil type and the initial amounts of the fluids present. Two-phase flows consisting of pulse applications
of water result in overlapping simple waves and contact discontinuities. These models form the basis
for future extension of the three-phase model to include pulse boundary conditions.
This report covers a period from October 1988 to September 1990, and work was completed
as of September 1990.
IV
-------
TABLE OF CONTENTS
Foreword Hi
Abstract iv
Rgures vi
Tables vii
1. Introduction 1
2. Conclusions and Recommendations 3
3. Research Scope and Background Literature 6
Review of Related Literature 7
Sharp Front Approximation for Infiltration 7
Method of Characteristics Solutions for Multiphase Flow 7
4. Physical Principles Governing Row in Porous Media 9
5. Model Formulation 18
Classical Method of Characteristics Solutions for Systems of Hyperbolic Equations 23
Generalized Method of Characteristics Solutions for Discontinuities 26
Riemann Problem Definition and Solution 28
On the Occurrence of Continuous and Discontinuous Waves 33
Summary of Approximate Governing Equations 41
Numerical Solution Methods 43
Construction of the Saturation Composition Space for Three-Phase Flow 44
Qualitative Nature of Three-Phase Flow 47
Construction of Saturation Profiles 48
Two-Phase Injections 50
Example 1 Oleic Liquid Bank Formation 50
Example 2 Oleic Liquid Bypassing 56
Example 3 Oleic Liquid Bank Formation, Second Type 56
Example 4 FlowofTCE 66
Example 5 Horizontal Flow 66
Example 6 Oleic Liquid/Air Injection 71
Parameter Variability 73
Single-Phase Water Injection 79
6. Generalization of Riemann Problem Solutions 82
Model Equations for the Infiltration of Water Subject to Air Phase Resistance 82
Generalized Total Flux 83
Overlapping Characteristic Patterns 85
7. Discussion of Results 92
References 95
Appendices 98
1 Derivation of the Three-Phase Fractional Flow Equations 98
2 Kinematic Model Solution 105
3 Horizontal Flow 109
4 Derivation of the Two-Phase Fractional Flow Equations 111
5 Input Data Sets 114
-------
FIGURES
1 Typical Capillary Pressure Curve 12
2 Typical Two-Phase Relative Permeability Curve 14
3 Sharp Front Approximation to True Spreading Front 22
4 General Boundary Transition Effects in Two-Wave Systems 34
5 Multiple Boundary Transitions in Two-Wave Systems 35
6 Abrupt Transition in Two-Wave Systems 36
7 Monotonically Increasing Wave Generating Function 37
8 Monotonically Decreasing Wave Generating Function 39
9 Composite Wave Generating Function 40
10 Example 1 X Eigenvalues 42
11 Example 1 X Eigenvalues 42
12 Example 1 Saturation Composition Space 45
13 Example 1 Saturation Routes 49
14 Example 1 Base Characteristic Plane 51
15 Example 1 Saturation Profile at 24 Hours 53
16 Example 1 Saturation Profile at 48 Hours 54
17 Example 2 Bypassing Profile 57
18 Example 2 Base Characteristic Plane 58
19 Examples Bank Profile 59
20 Example 3 Saturation Routes 60
21 Profile Transition from Bypassing to Bank Forming 61
22 Example 4 TCE Saturation Composition Space 67
23 Example 5 Horizontal Flow Saturation Composition Space 68
24 Example 5 Comparison of Horizontal and Vertical Bank Profiles at 1 Day 69
25 Example 5 Comparison of Horizontal and Vertical Bank Profiles at 2 Days 70
26 Example 6 Oleic Fluid Injection 72
27 Saturation Composition Space for Hydraulic Conductivity of 3.35e-3 cnVs 74
28 Effect of Saturated Hydraulic Conductivity on Bank Profiles 75
29 Effect of Saturated Hydraulic Conductivity on Bypassing Profiles 76
30 Saturation Composition Space for Brooks and Corey's Lambda - 0.30 77
31 Effect of Brooks and Corey's Lambda on Bank Profiles 78
32 Effect of Brooks and Corey's Lambda on Bypassing Profiles 80
33 Water-Air Infiltration Base Characteristic Plane 87
34 Water-Air Infiltration Profile at 0.2 Days 88
35 Water-Air Infiltration Profile at 0.3 Days 90
36 Water-Air Infiltration Profile After Overtaking 91
VI
-------
TABLES
1 Mass Balance Errors for Examples 1 to 3 62
2 Shock Inequalities Satisfied by Examples 1 to 3 63
VII
-------
SECTION 1
INTRODUCTION
The focus of the research described in this report Is the transport of oily liquids within the
subsurface. Subsurface contamination by oily wastes is recognized as a ubiquitous problem which
presents some of the most challenging remediation difficulties. The three-phase flow of oily liquids in
the unsaturated zone comprises one aspect of this problem. These fluids act as carriers for, or are
themselves, hazardous organic chemicals. Motion of the organic liquid phase can enhance transport
of hydrophobic organics, by orders of magnitude over that occurring when water is the only mobile
liquid.
Transport in the unsaturated zone is characterized by transient water flow caused by discrete
rainfall events of varying Intensity, duration, and inter-storm arrival time. Varying water fluxes play a
central role in multiphase subsurface flow, as the rainfalls provide a driving force for the hydrologic
Inputs to the system. A focus of the research discussed in this report is the effect of transient flows in
the water, organic, and air phases.
Since the governing equations for this problem are a system of coupled non-linear partial
differential equations, numerical methods must be used for their solution. Existing complex models of
multiphase flow can be CPU time intensive and their proper application is data intensive. Numerical
solution techniques sometimes introduce difficulties, since they are not approximation-error free.
During the initial phase of site investigation, or for alternatives screening, a simplified approach is
warranted, because data from typical RCRA and CERCLA sites Is usually inadequate to justify a
complex model. The approach that is taken in this report is to represent the flow by a system of
hyperbolic equations, which describe the wave behavior of the flow syslem. These equations can be
solved by the method of characteristics. As will be shown, this technique results in models which have
semi-analytic solutions. The latter are important because their analytic character allows rapid
computation of many alternative scenarios. The models also have intrinsic value for checking
-------
numerical simulators. Most important, however, Is the Insight they provide into the fundamental
character of the flow.
As stated above, the phases modeled in this report are water, a water-immiscible organic
liquid and air. The term oleic liquid is used to Indicate the organic liquid. The commonly used
acronym NAPL (Non Aqueous £hase Liquid) is avoided, because some nonaqueous liquids are
completely misdble wtth water (i.e. alcohols) and as such are not considered here. The word oteic Is
used because tt means "oily,' and carries the connotation of immisdbility wtth water. The terminology
and the models presented herein are not restricted to petroleum based oleic liquids, but also include
organic solvents and other liquids.
The remainder of the report is organized as follows. Section 2 contains a summary of the
conclusions and recommendations. Section 3 describes the background and motivation for the work.
Related models are briefly discussed. Section 4 reviews the multiphase, porous media, flow theory
upon which the models are based. In section 5 the classical and generalized method of
characteristics solutions of hyperbolic equations are presented. The solution of a "Riemann problem"
of mathematical physics is presented. The solution methodology is applied to one-dimensional, three-
phase flow in soils. The fundamental character of the flow Is shown through a series of Illustrative
examples. Since Riemann problems have restrictive boundary and Initial conditions, the
generalization to pulse boundary conditions is discussed in section 6. Application of this theory is
made to the infiltration of water subject to air phase resistance. A discussion of the results is
presented In section 7.
-------
SECTION 2
CONCLUSIONS AND RECOMMENDATIONS
The application of characteristic methods to one-dimensional flow problems reveals the
fundamental behavior of three-phase systems. Although the full governing equations are parabolic,
the approximate hyperbolic equations which are solved in this report describe a fundamental portion of
the flow phenomena. The following conclusions are drawn from the work.
1) For three-phase systems wtth constant Injection conditions and uniform initial conditions,
flow occurs in two regimes. When the injection consists mostly of water, mobile oleic liquids
(organic liquids or so-called NAPLs) are displaced into "banks" or are bypassed. The banks are
regions in the profile where the oleic liquid fills a higher fraction of the pore space. The banks
move ahead of, and are driven by, the incoming water. The bypassing profile is characterized
by the injected water moving past the oleic liquid, without causing it to be displaced deeper into
the profile. The occurrence of the regimes is determined by the oleic liquid properties, media
properties, and the initial amount of all fluids present. Numerical solution of the model
quantifies these phenomena.
2) That the character of the results depends on the initial condition suggests that for cases
where the oleic liquid has drained to a low saturation in the upper part of the soil profile,
incoming rainfall will likely bypass the oleic liquid. This conclusion holds even for oleic liquids
with relatively high mobilities in the subsurface. Conversely, water injections wtth relatively high
initial oleic liquid saturations are dominated by the bank forming regime.
3) For systems wtth injections of oleic liquid, the dominant flow regime is water bypassing.
This conclusion confirms the intuition that the oleic liquid, which occupies a mid-range of the
pore space due to Its wetting behavior, does not displace water from the small pores.
Preferential wetting causes the water to be too tightly held to the media to be easily displaced
by a non-wetting liquid.
-------
4) In all cases, the behavior of the oleic liquid Is dependent on Its transport properties (density
and viscosity) relative to water. As the oleic liquid becomes more mobile, water injections favor
the bank formation regime, because the oleic liquid is of high enough mobility to match the
speed of the incoming water. The character of the results, however, always reflects preferential
wetting as indicated in conclusion number three.
The following recommendations are derived from this work.
1) The three-phase solutions should be extended to simulate conditions with more general
boundary and initial conditions. The methodology presented for two-phase flow simulation In
section 6 appears to be adaptable for this purpose.
2) Laboratory validation of the existence of sharp fronts would be valuable for demonstrating
the usefulness of the model. Absolutely sharp fronts are not necessary, for the model,
because the sharp fronts can be shown theoretically to move at the same speed as the true
fronts. Experimental verification of the bypassing and bank forming regimes is needed as this
work proposes their existence from purely theoretical considerations. A non-destructive
imaging technique, such as dual gamma ray attenuation, is necessary for tracking the time
dependent fluid saturations. Also, experimental work would indicate the importance of
hysteresis in the constitutive relations governing the flow.
3) The laboratory work should include investigation of the saturation conditions at the ground
surface. Somewhat arbitrary conditions are used in the work discussed in this report. A related
fesue is that the injections used here force flow to be unidirectional. Experimental work
indicates countercurrent flow is common, so the model should also be adapted for this
condition.
4) Useful models for field problems should be constructed from the immiscible transport
model developed herein and the extension proposed in Kern 1. Such models must include
interphase partitioning phenomena. Previous work on kinematic models demonstrates that
such phenomena are amenable for inclusion into this type of model.
-------
5) Mufti-dimensional extension of the models should be considered, since one-dimensional
flow is a restrictive condition which doesn't apply to heterogeneous field sites. The one-
dimensional models can be useful, however, for spills occurring over large areas.
-------
SECTION 3
RESEARCH SCOPE AND BACKGROUND LITERATURE
The specific motivation for this work is that previous attempts to model multiphase flow of
contaminants by simplified methods are limited to plug flow or to kinematic conditions. Plug flow
models require the assumption that the fluids occupy a fixed portion of the pore space. Kinematic
models allow more realistic fluid distributions, but are limited to gravity driven flows. As a result, an
Intense rainfall, high viscosity oleic phase, and/or tow permeability soil may cause the model to break
down (Charbeneau et at., 1989). In the present work, the general approximate framework of the
kinematic approach Is retained and extensions to overcome the model's limitations are sought The
models that are presented allow some dynamic phenomena to be included along with the kinematic or
gravity phenomena.
In this report two types of solutions are presented. First is the solution to a classical problem of
mathematical physics, called a Riemann problem. Mathematically, a Riemann problem is defined as a
hyperbolic system of equations with constant injection and uniform initial conditions. Although this is
of limited direct applicability to field problems, the solution of the Riemann problem is important for
several reasons: Mathematically, there has been a great deal of work done on these solutions for
systems other than those presented here. In this report, the application of this work to multiphase flow
is explored. Second, the solutions of the Riemann problem reveal some of the fundamental character
of the flow. Results will be presented which Illustrate several typical flow regimes for three-phase
flows. Third, one-dimensional solutions of the Riemann problem form the basis for multi-dimensional
front tracking models of petroleum reservoirs (e.g., Qlimm et al., 1983) which may be adaptable for
subsurface contaminant transport
The second type of solution presented In the report is for systems which can be represented by
injection conditions which consist of discrete pulses. The model system for this application is the
infiltration of water, subject to air phase resistance. The solution of this problem consists of an
extension and generalization of the Riemann problem solution techniques. The solution which is
6
-------
presented is the first step in developing a model for a surface spill of an otete contaminant, followed or
preceded by a sequence of rainfall events.
REVIEW OF RELATED LITERATURE
Models which are related to method of characteristics solutions presented In this report are
reviewed below. These include sharp front approximations and method of characteristics solutions for
one and two phases.
Sharp Front Approximation for Infiltration
Green and Ampt (1911) used a sharp front at the leading edge of infiltrating water, and the
assumption of constant moisture content above the front to develop an analytical solution for constant
head infiltration. Mull (1970) developed a sharp front model of oleic liquid profiles to simulate non-
uniform redistribution of the oleic liquid after the release ends. Each profile in Mull's model consisted
of oil filling a fixed amount of the pore space from the surface to the front. A series of such uniform
profiles constituted the model's result. Reible et al. (1990) used a sharp front at the leading and
trailing edges of an infiltrating oleic liquid to develop a similar simple model. The water saturation was
assumed to be residual so that the oleic liquid displaces air only. Parameters from the model were
fitted to laboratory column data.
Method of Characteristics Solutions
Smith (1983) and Charbeneau (1984) developed similar models for approximate solution of
Richards' (1931) equation. These models assumed that the flow Is kinematic (or equivalent^ driven
by gravity), so the governing equation can be approximated by a hyperbolic equation which allows
solution by the method of characteristics. Weaver and Charbeneau (1989) extended this methodology
-------
for Infiltration of an oleic liquid with a constant water saturation. The last three models have an
advantage over the previous models, because they do not require a uniform moisture or oleic liquid
profile. All of the models mentioned above share the characteristic that they are essentially single
phase flow models. The presence of the other phase(s) is Included implicitly, but their flow is not
t
modeled. Charbeneau et al. (1989) developed a kinematic model which Included the simultaneous
unsteady flow of water and oteic liquid; the presence of air in the unsaturated zone was included, but
Us flow not explicitly considered. The scenario modeled was a surface release of an oleic liquid,
followed by a series of discrete rainfall events. The model, being limited to gravity flow, breaks down
when fluxes exceed the gravity flux. An advantage of the approach, and also that of Weaver and
Charbeneau (1989), Is that the model executes rapidly on small computers.
The simplest solution of a Riemann problem for multiphase flow is the two-phase Buckley-
Leverett (1942) solution, which consists of a simple wave and a contact discontinuity. The work
presented here is essentially a generalization of the Buckley-Leverett model to three-phase flow
including the effects of gravity. Helfferich (1981) presented the theory of "coherence" for multiphase,
multicomponent porous media flows. For systems of hyperbolic equations, the method of coherence
is equivalent to the method of characteristics solution of the Riemann problem. Helfferich (1986)
believed that the method of coherence contains the solutions of systems of hyperbolic equations as a
special case, but applies in general to other systems also. Pope et al. (1978) presented solutions with
straight paths across the saturation composition space, which is an essential difference wtth the work
presented here. Dougherty and Sheldon (1964) presented a similar Riemann problem solution, but
restricted their work to only one of the dominant flow regimes presented in the present report, and to
horizontal flow.
8
-------
SECTION 4
PHYSICAL PRINCIPLES GOVERNING FLOW IN POROUS MEDIA
The equations thai govern multiphase flow fa porous media are coupled, nonlinear, partial
differential equations (PDEs). There is a mass conservation equation for each fluid, auxiliary
relationships which incorporate various physical phenomena into the PDEs. By nature, these relations
are empirical, since the phenomena are highly complex and thus are not suited to exact mathematical
expression. Taken together these equations form a mathematical model of the flow (e.g., Peaceman,
1977).
A typical mass conservation equation may be written in one dimension as
dJ
im _m
dt dz " pm
3
where m [M/L ] is the mass per bulk volume of porous media, J [M/T] is the mass flux, and p
3
[M/TL ] is the source strength per unit volume. Throughout the rest of this work, the liquids will be
assumed Incompressible and the medium non-deformabte. The resulting continuity (or phase
conservation) equation is
as.
i source
where q. [L /L \\ is the volume flux of fluid i per unit area, TI [L /L j is the porosity, p* [1/T] Is the i
3 3,
strength of fluid i divided by the fluid density of i, and S [L /Lj is the fraction of the pore space filled by
fluid i. S. te called the degree of saturation of fluid i, or simply the saturation. In a multiphase flow
system,
-------
IS,- 1
where i represents each fluid present. For water, an oleic liquid and air, equation 3 becomes S + S +
Wr w
S - 1. In order to model the fluids wfth the phase conservation equation (2), the composition of each
a
phase is assumed to have no effect upon any fluid's transport properties. Fluids for which partitioning
phenomena can significantly alter the saturation during the period of simulation are not suited to this
approach.
The flux term of equation 2 can be related to fundamental properties of the system through
use of Darcy type equations. Specifically,
kkri(S)
where z is directed positive downward, k [L ] is the intrinsic permeability of the media, k Is the relative
ri
permeability of the medium to fluid i, n. [M/LT] is the dynamic viscosity, p [M/L ] the density, g [L/T]
te the acceleration due to gravity, and P [M/L'T] is the fluid pressure. The quantity K (S) [L/T] is
i 0!
called the effective hydraulic conductivity of the media to fluid i and is defined by
k k (S) p g
Ke|(S) - " ' - KS| kri(S) 5
where Kg| [L/T] is the fully saturated hydraulic conductivity to fluid i. The last equation underscores
that the intrinsic permeability, density, and viscosity play an important role in determining the
conductivity, and thus the fluxes, in multiphase flow systems. The relative permeability Is the ratio of
conductivity at any saturation to the fully saturated conductivity and thus is dimensionless. In these
10
-------
last two equations no subscript appears on the saturation, because k (S) may or may not depend on
the amounts of the other fluids present.
Equation 4 is an extended form of Darcy's law. It incorporates the physical phenomena
occurring in a porous medium containing immiscible fluids into a form that can be used in the
continuity equations. Thus, knowledge of the geometry of each individual pore and of the fluid
Interactions within the matrix is not needed. As outlined below, this information Is contained in an
Integrated sense in the relative permeability and capillary pressure functions. The dependence of
these functions on saturations of two or more fluids results in coupling of the conservation equations.
In multiphase systems the fluid pressures are related by the capillary pressure, P . This
c
fundamental relationship is defined for each pair of fluids present by
P - P - P
c rrw w
where the pressure, P , Is the pressure in the wetting phase - the fluid most strongly attracted to the
^N
solid. P is the pressure in the nonwetting fluid. For simple geometries (i.e., single circular tubes),
flW
there are analytic expressions for the capillary pressure (e.g., Bear, 1972, pg. 446). Soils, however,
are composed of a variety of irregular, tortuous pores. The relationship expressed by equation 6 still
holds, but must be expanded to represent the variation in P which can occur over a representative
C
elementary volume. In contrast to media composed of uniform geometries, certain sized pores are
filled by certain fluids, leading to varying contributions to the overall pressure difference.
A typical capillary pressure function for a strongly wetted system Is shown hi figure 1. This
curve shows that the capillary pressure is infinite when the wetting phase is at the residual saturation.
In three-phase systems, there are three possible capillary pressure curves, one for each pair of fluids.
Only two of the capillary pressures are independent, however, as the third can be obtained from the
other two.
11
-------
Capillary
Pressure
Drainage
Imbibition
0.0 1.0
Water Saturation
Figure 1 Typical two-phase capillary pressure curve showing hysteresis.
12
-------
Numerous characteristics of both the capillary pressure and relative permeability functions can
be related to the underlying physical phenomena (Weaver, 1984, pp. 57-60, 62-76). Interpretation of
directly measured capillary curves reveals Information concerning these phenomena. Notably,
capillary pressure (and relative permeability) are subject to hysteresis. That is, the capillary pressure
depends on the displacement direction. As indicated by the arrows tn figure 1, It te higher during
drainage (displacement of a wetting fluid by an nonwetting fluid) than during imbibition (the opposite
displacement).
The relative permeability function, k (S) te the fraction of the fully saturated hydraulic
conductivity existing at any saturation. Numerically, II Is a number from 0 to 1. K quantifies the notion
that the presence of other fluids reduces the hydraulic conductivity of the medium to a given fluid.
Physically, this occurs due to reduction in the numbers of pores available to the fluid and due to the
increased tortuosity of those remaining. These effects are illustrated by the typical two-phase relative
permeability curves in figure 2. The permeabilities drop off rapidly In a characteristic manner, when
the saturation is reduced. At the trapped or residual saturations, S the relative permeability is rero,
implying that all of the fluid exists as a discontinuous phase. Residual wetting fluid Is held in the
smallest pores, while residual nonwetting fluids are found in larger pores as isoiated droplets.
Due to the considerable difficulties of measurement, three-phase relative permeability is rarely
measured (Stone, 1970). Alternatively, various models of relative permeability are used. One way to
develop these Is to begin with conceptual models of the porous medium, presumed distributions of
fluids wtthin the medium, and a solution of laminar flow through the medium (Bear, 1972, pg. 463).
The Burdine equations form one such model (Burdine, 1953, Wyllle and Gardner, 1958); for the
water, otelc liquid, and air phases, the relative permeability equations are:
r w
s
w
1
s
wr
S
wr
*-
f " '
dS
. dP2 /
0 c
7a
13
-------
1.0
Relative
Permeability
0.0
0.0
Saturation
1.0
Figure 2 Typical two-phase water/oil relative permeability curve showing hysteresis.
14
-------
ro
's -
0
1 -
S +S
s
or
s
nr
-------
where S te the effective water saturation, P is the minimum capillary pressure required to force a non-
e ce
wetting fluid into a fully saturated medium, and X is a dimensionless parameter that is Indicative of the
pore size distribution. Using equation 8 in equation 7 results in the following model of drainage relative
permeability for water, the oleic phase and air:
fS -S
, w wr
krw-ll -S J
wr
12 +
where e - ^ '
fs -s I2ffs + s - s le'2fs -s E'2
. _p_ pi o w wr w wr
kro-l 1-S J II 1- S J ' I 1-S
or wr wr
fs -s ]2r fs +s - s e'2
L. a ar o w wr
kra"l 1-S J I1'I 1- S
ar wr
Despite its apparent analytic nature, certain parameters must always be measured to fit the model to a
specific situation. For equations 9, 10, and 11 the parameters are the residual water, oleic liquid, and
air saturations, and the Brooks and Corey exponent. The values of S , X, and P are obtainable
wr ce
through a procedure outlined by Brooks and Corey (1964). They are essentially the result of a
capillary pressure curve measurement. An extensive tabulation of these parameters for many different
soil types has been developed, making this an attractive model (Brakensiek et al., 1981, and Rawis et
al., 1983).
The wetting fluid relative permeability (equation 9), is the same as would result from a two-
phase integration of the Burdine equations. This is because in a strongly wetted system, the strongest
wetting phase (say water) is most strongly attracted to the solid. Thus, water is assumed always to
be found in the smallest pores and grain contacts, independent of the non-wetting fluid(s)
16
-------
composition. The most nonwetting fluid (air) is repelled by the surfaces so it is found in the largest
pores. The attraction of the oleic liquid to the surfaces is Intermediate between that of water and air,
the amounts present of which determine the size of the pores occupied by the oleic liquid. Visual
studies confirm this distribution for some systems (e.g., Schwille, 1988).
17
-------
SECTION 5
MODEL FORMULATION
The three-phase flow model is based on the assumption that the medium is uniform and the
flow is one dimensional. Each phase has a residual or trapped saturation which remains constant.
The relative permeabilities are described by equations 9 through 11. All phase transport properties are
assumed to remain constant. The model describes immisdble transport without interphase
partitioning phenomena, or any sources/sinks in the domain of Interest. Assuming incompressible
flow in non-deforming uniform media, the phase conservation equations (2) for the water, oleic liquid,
and air are
dS dq
+ _ - 0 12
dt dz
dS dq
o o
TI . + « 0 13
dt dz
and
14
dt dz
Mass is conserved in all three phases, since the source terms are set to zero. Summing gives,
S+S+S+ q+q+
dll a ° WJ dzl a °
18
-------
By defining the total flux, a as, q- q +q + q, and recalling that the three fluids fill the pore space,
I.e., S + S + S « 1, equation 15 reduces to
o O W
dqt
__ -0 16
dz
Implying that for incompressible, multiphase flow the total flux is independent of depth. Equation 16
makes no statement concerning the time dependency of a. For the remainder of section 5, a is taken
as a constant, while time dependent total fluxes are discussed in section 6. Normalized fractional
flows, f., can be defined as
f.-q/q.
where the subscript refers to fluid i. Note that the fractional flow represents the fraction of the total flux
which is contributed by the flux of fluid i. Obviously,
f +f +f
a o w
A consequence of depth-independent total flux and filled pore space is that any one of the three
conservation equations for the fluids may be eliminated from the system of equations. For example,
the air equation (14) can be written as,
d(1-S -S ) d(1-f -f )
* W O W O -
r\ +a -0 17
dt dz
which is immediately reduced to an identity by noting the phase conservation equations for water and
oteic liquid. Thus the resulting model consists of a system of only two coupled mass conservation
equations in two unknowns. In the following discussion, the water and oleic liquid equations are used
19
-------
hi developing the solution. At times, use of the water/air and oleic liquid/air formulations is convenient
and is noted in the report.
Expanding the spatial derivatives in terms of the saturations, S and S , in the water and oleic
Hquid conservation equations (12 and 13) gives
as
w
at
+ s
df dS df dS
w w wo
dS dz
I W
dS dz
o
-0
18
and
ds
dt
df ds df dS
00 O W
as
dS dz
w
-0
19
where the fractional flow derivatives can be expressed tn terms of fundamental multiphase transport
parameters, as discussed below. Writing the equations in matrix form gives
20
dt
dz
where the saturation composition vector, S, is given by
and the matrix, A(S),
20
-------
q
A(S)-_I
TI
r df
w
aS
w
af
Q
L dS
w
df ,
W
as
o
df
Q
as J
0
Specific values of Swill be denoted by S(S , S ,S ), where numerical values replaces , S , andS (S
W 0 3 W O a a
is included in this notation only for convenience, tt is recognized that S is not properly a part of the
a
vector S). The model equation 20 Is a two by two system of quasi-linear partial differential equations.
Semi-analytical solutions of the model equations are sought, as they are useful In revealing the
behavior of the system. Consequently, attention Is restricted to uniform porous media.
The fractional flows and their derivatives can be related to fundamental transport properties by
the use of Darcy's law equations for the phases. Fractional flow functions in terms of water and oleic
Hquid parameters are presented In appendix 1. For the system of equations presented here, there are
two independent expressions for the fractional flow of the oleic liquid,
f . -f Jk ,k ,f
o1 oV ro rw
f -< o(k -k -f )
o2 o2y ro ra wr
where k Is the fraction of the saturated conductivity to fluid I. The relationships are linear and easily
solved for f (f «f »f ) and f , given S. The partial derivatives of the fractional flows are likewise
ov o o1 o2; w" r
given by systems of two linear equations In two unknowns (appendix 1). All coefficients appearing In
equation 20 are either constants or are determined by S, which Is the only unknown appearing in that
equation.
As indicated in appendix 1, the model neglects the gradients of the capillary pressure by
omitting them from the fractional flow functions. This assumption is critical for eliminating the parabolic
character of the phase conservation equations, so that the method of characteristics can be applied to
this problem. One of the main effects of the capillary gradient Is to smooth sharp fronts. Figure 3
21
-------
Saturation
Depth
""- F
Rgure 3 Sharp front approximation to a true spreading front. The sharp front is selected by
mass conservation to match the speed of the true front between points 1 and 2.
22
-------
shows the relationship between a true smooth front and the sharp front approximation used tn this
solution. Because of the way that the speed of the sharp front Is determined, the mean displacement
speed of the true smooth front matches that of the sharp front (Charbeneau, 1984). Row visualization
experiments conducted at the Robert S. Kerr Environmental Research Laboratory and elsewhere
(Reible et al.t 1990) suggest that infiltrating fronts indeed remain sharp for a variety of situations.
CLASSICAL METHOD OF CHARACTERISTICS SOLUTIONS FOR SYSTEMS OF
HYPERBOLIC EQUATIONS
The following section discusses solution of systems of quasilinear equations. The primary
references on this material are Jeffrey and Taniutl (1964), Lax (1957, 1973), RozdestvensWi and
Janenko (1980), and Smoller (1983). For readers whose main interest is the results of the model, the
next three sections may be skimmed. The definition of the Riemann problem must be noted, however,
because this is the type of problem solved in this section. Application of the theory to multiphase flow
problems begins to be discussed In the section entitled "On the Occurrence of Continuous and
Discontinuous Waves." A summary of the equations used for each part of the solution is given in the
next section ("Summary of Approximate Governing Equations").
For an arbitrary quasilinear system of equations, a characteristic form can be derived as follows.
*^
The k toft eigenvector, 1 of the matrix A(S) is defined by
K
lk(A(S)-Xkl)-0
Al_
where the scalar, X , is by definition the k eigenvalue of A(S), and I is the identity matrix. For
systems of two equations in two unknowns, the eigenvalues of matrix A(S) are denoted by X and X ,
and defined as
23
-------
21
and
v;
22
Note that X , X , or X refer to eigenvalues, while X refers to Brooks and Corey's exponent (equation 8).
1 ^ K
The components of the teft eigenvectors are determined by solving
and arbitrarily fixing the value of 1 (1) or 1 (2). Multiplying equation 20 by i , and noting the definition
K K K
of the left eigenvector, gives
dS dS
The substantial derivative of S, DS, is the sum of the temporal and spatial derivatives,
dt dt dz
Thus the classical method of characteristics form for equation 20 is simply seen to be
lkDS-0
23
24
-------
along paths, which are called characteristics, defined by
where k - 1 , 2 for the two by two system. There may be some systems for which equations 23 and 24
have analytical solutions, resulting in dosed-form relations for the independent variables. For three-
phase porous media flow, however, an analytical solution te not known and numerical solutions of
these equations must be found. Application of the method of characteristics, however, reduces the
system of n coupled partial differential equations to a system of 2n coupled ordinary differential
equations. In general, the latter are easier to solve numerically than the former. The relationship
between the kinematic model solution for three-phase flow presented by Charbeneau et al. (1989) and
Weaver (1988) and the classical method of characteristics solution Is discussed in appendix 2.
For a two by two system, equation 23 becomes
k.1,2 25
dS1 V2>
in each k direction specified by equation 24. Since there are two eigenvalues of A(S), all four ordinary
differential equations represented by equations 24 and 25 must be solved simultaneously to determine
S(z,t). Such solutions are defined as classical solutions of equation 20 by the method of
characteristics.
Two important caveats must be noted. First, a necessary condition for the existence of classical
solutions determined from (equations 23 and 24) is that the system remains hyperbolic throughout the
solution domain. The condition for hyperbolicity is that the eigenvalues remain real and distinct
throughout the solution domain. For the multiphase flow problem considered here, this will be
demonstrated below. Second, the classical solutions are guaranteed to exist only in the 'small," i.e., in
a small neighborhood of the initial data. Proper initial data does not prevent discontinuities from
forming within the solution domain. At these discontinuities, the partial derivatives appearing in
equation 20 do not exist, so neither do the classical solutions. In order to overcome this difficulty,
25
-------
solutions which satisfy integral forms of the mass conservation equations are used. In the present
work, these are used where the sharp fronts replace the true fronts. The classical and generalized
solutions are patched together to construct the complete solution, which consists of regions with
smooth saturation variation separated by sharp fronts.
GENERALIZED METHOD OF CHARACTERISTICS SOLUTIONS FOR
DISCONTINUITIES
When the saturation derivatives, dS/dt and dS/9z, fail to exist wtthin the solution domain, the
classical method of characteristics solution presented In the previous section also fails to exist. The
differential mass conservation equations 12 to 14 are no longer meaningful and must be supplanted
with generalized (also called integral or weak) solutions. By Integrating the mass conservation
equations (e.g., Charbeneau, 1984} an integral solution that conserves mass is found. For example,
a discontinuity in the water phase obeys
s
dt
w2 w1
Sw2" Sw1
26
where s is the speed of the discontinuity, z is the location of the sharp front, and the subscripts 1 and
2 refer to points on either side of the discontinuity (Figure 3). Analogous equations exist for
discontinuities in the oleic and air phases.
Discontinuous solutions are notorious for being non-unique and abundant examples of this
behavior exist (e.g., Jeffrey and Taniuti, 1964, pp 119-121). Generalized entropy conditions (shock
Inequalities) are used to select physically realistic discontinuous solutions. Two types of discontinuous
solutions are important for this report, k-shocks and contact discontinuities. By obeying the shock
inequalities, a discontinuous solution is assured to be physically realistic and therefore the proper
discontinuous solution. The shock inequalities define k-shocks, and are given by (Smoller, 1983)
26
-------
VS2>
28
where the subscripts 1 and 2 on S refer to points behind and ahead of the discontinuity, respectively.
For a two by two system wfth X < X , either
29
and
s - 0 33
27
-------
where the bracket, < , >, is defined as the inner product of two vectors. Linear degeneracy in a
characteristic field implies that the system, although nonlinear, behaves like a linear system where
shock speeds always match characteristic speeds. For the multiphase problem VX , has not been
K
calculated explicitly. Instead of showing that the characteristic fields are linearly degenerate by
satisfying equation 33, the examples below demonstrate that the characteristic speed of the smooth
wave matches, to wtthin numerical error, the shock speed. In fact, the contact discontinuities
presented below are 'almost" 1-shocks, since they satisfy equation (29) and violate equation (30) only
because s - X (S ).
RIEMANN PROBLEM DEFINITION AND SOLUTION
An important class of problems is defined by equation 20 and the boundary and initial data,
S(O.t)
S(z,0)
where S is the boundary saturation composition and S is the initial saturation composition. These are
called Riemann problems, and are the focus of this section. For the Riemann problem, the
characteristic form given by equations 23 and 24 simplifies as discussed below.
Momentarily restricting the discussion to small neighborhoods of the Initial data for hyperbolic
systems, the following definition is made. A function w Is defined as a k-Riemann Invariant of the
system, if ft is a smooth function, w(S ,S ), that satisfies
W O
"
4h ll,
where r is the k right eigenvector of A(S). The k right eigenvector is defined by
K
28
-------
(A(S) - Xl)r - 0 35
The components of the right eigenvector r are determined by solving the homogeneous
f\
equations.
By assigning r (1) - 1, and using the first of the above equations, r (2) is found to equal -a /(a - X ).
K K 1 fc 1 I K
Using the second equation to determine r (2) can be shown to give the same result.
K
Details of the following results for hyperbolic systems are presented by Smoller, and form the
basis for construction of the solutions. When S is a C continuous solution of 20 in a domain D, and all
k-Riemann invariants are constant, then S is defined as a X -simple wave. When the X -simple wave
K K
depends only on a parameter £ (z-z )/(t-t ), the wave is defined as a X -centered simple wave. The
wave is centered at the point (z ,t ). Simple waves are important for the Riemann problem, because
they are one one of the main features of the solutions.
For equation 20, a Riemann invariant (equation 34) can be written as
w
which leads to
29
-------
&L rk(1)
36
When the Riemann invariants are constant, i.e., within a simple wave, the implicit function theorem
shows that
dw dw dS
Q
dS +dS dS "
w o w
Now equation 36 can be rewritten in terms of the unknown values of the variables, S and S , as
\nr O
dS
°- . 37
dS r.(1)
w kv '
Equation 37 applies throughout the simple wave, notably across the characteristics. For a given
eigendirection, k, equation 37 is a single non-linear ordinary differential equation. The non-linearity
results from the right eigenvectors being functions of the saturations. The direction in saturation
composition space given by this equation varies for a given eigendirection, because the eigenvalue
and eigenvector vary with the saturation composition, S.
A similarity solution of equation 20 is possible because the problem has no characteristic time (t)
or depth scales (z) (Charbeneau, 1988). The two independent variables z and t collapse into a single
Independent variable £ - zA. Equation 20 is transformed by the relations
dt d£ dt "
as ds &
3z " d£ dz
30
-------
resulting In
(A(S)-$I)S -0 38
where the initial data is transformed as S(o) = S , and S(») - S . When £ Is not an eigenvalue of A(S),
S must be equal to zero. Thus the saturation change vector equals zero and solutions are admitted
which consist of regions in z-t space with no change in saturation composition. Such regions of
constant state are called plateaus.
As a consequence of the definition of the right eigenvector, H £ is an eigenvalue of A(S), then S
must be a right eigenvector of A(S). This relationship is noted by comparing the definition of the right
eigenvectors (equation 35) with equation 38. Further, for a given eigenvalue, X , the eigenvectors
K
satisfy a bi-orthogonallty relation, (e.g., Young and Gregory, 1988)
V,-o "
So, for the Riemann problem where r - S ,
k s
l.r. = t. S
k k k
The classical method of characteristics form (equation 23), which states that I DS - 0 along dz/dt -
K
X , Is satisfied when DS - 0, since i * 0, and as shown In the last equation 1 S * 0. For the
k k k ^
Riemann problem
dz
DS - 0 along - X 39
dt k
31
-------
The characteristic form for the Riemann problem (equation 39) shows that the saturation change
vector DS is zero along the characteristics, implying that saturation compositions remain fixed along
each characteristic. Since K - V(S), the characteristics must be straight lines. From a computational
K K
viewpoint straight characteristics are convenient features, since equation 24 has a simple analytical
solution and numerical integration is not needed.
Equation 37 provides a powerful and necessary tool for characterizing the possible saturation
compositions. By picking an eigendirection, k, the saturation composition across the wave is
determined by solving equation 37. A map of ail possible states, S, which satisfy equation 37 given a
starting saturation composition can be drawn. Such S maps are discussed below in the section
entitled, "Construction of the Saturation Composition Space for Three-Phase Row." Since during the
construction of the S map, all waves are assumed to be continuous, modification of the routes is
required before using the map to construct solution profiles containing shocks.
The variation of saturation along a shock can be determined by equating the jump equations.
For example, the shock speeds are the same for discontinuities in water and oleic liquid saturations,
so
dt
f
w
S
I w
-f*
w
-S*
w J
-i
T\
f -f*
0 O
S -S*
I o oJ
where (S* ,S*) is a fixed point in the saturation space, and f* and f* the values of f and f evaluated at
wo o o w o
(Sw,S^). As a result, the water and oleic liquid saturations across the discontinuity are related by,
f -r
S - S* + :*% (S - S*)
o o f -r v w w'
w w
40
which is used to adjust the routes across the saturation space when the solution becomes
discontinuous.
32
-------
ON THE OCCURRENCE OF CONTINUOUS AND DISCONTINUOUS WAVES
Since there are two model equations, there are two eigenvalues of the matrix A(S). Thus two
waves are created for each change in boundary condition. Following Hetfferich (1986), figure 4 shows
a smooth transition from constant state 1 to constant state 2. Since each value of saturation
composition on the boundary is associated with two characteristic directions, two waves are created
by the boundary transition. One wave corresponds to X and the other corresponds to X . Note that
this te not a Riemann problem and the general method of characteristics solution (equations 23 and 24)
must be integrated throughout the triangular region ABC to find the solution shown in this figure. The
waves shown are X-simple waves which define a continuous solution of equation 20. When more
transitions are included, more regions with overlapping characteristics appear (figure 5). In a Riemann
problem the boundary transition between the states is reduced to a single point, so there is immediate
resolution Into two distinct waves. Since they originate from a common point, the waves are X -
centered simple waves. This condition is illustrated in figure 6. All of these figures assume that X and
X decrease smoothly from state 1 to state 2 so that the characteristics form typical fan-shaped simple
wave patterns.
Either one or both of the waves may be replaced wholly or in part by a discontinuity which
corresponds to k-shocks or contact discontinuities respectively. Since the discontinuities can be
created instantly, there may be no continuous portion of the wave, as in the case of the k-shock.
Rgures 7 to 9 illustrate, in a schematic fashion, functions which generate simple waves and
discontinuities for both smooth and abrupt boundary transitions. In figure 7a a monotonically
Increasing X(S) function is shown. If there is a gradual transition from a low (S ), to a high (S ),
saturation composition, a simple wave is generated along the boundary, as each eigenvalue increases
with Increasing S. Since characteristic speeds are equal to the eigenvalues, each successive
characteristic has a greater slope than Its predecessor, causing the characteristics to separate (figure
7b). When the transition zone shrinks to a point (i.e., an abrupt transition of a Riemann problem) the
smooth wave still exists, as shown In figure 7c. Rgure 7c shows a centered simple wave as the
characteristics originate from a common point. Solutions for soil moisture (Charbeneau, 1984) and
33
-------
Constant State
Multiple
Simple
Waves
After Helffrich (1986)
Figure 4 General boundary transition effects In two-wave systems. A gradual transition
between constant state(1) and constant state(2) creates two overlapping simple
waves (triangle ABC).
34
-------
Figure 5 Multiple boundary transitions in two-wave systems. Multiple transitions create regions
with overlapping simple waves which originate from different boundary transitions
(DEFG).
35
-------
Starting
Curve
Constant State
Figure 6 Abrupt transition in two-wave systems. When the boundary transitions shrink to a
point there is immediate resolution into simple waves.
36
-------
A
S1 S2
t
Figure 7 Monotonically increasing wave generating function. Figure 7a shows eigenvalues
which are monotonically increasing with S. Figure 7b shows a gradual boundary
transition (in time) from constant state, S
to constant state, S which results in a
simple wave. Rgure 7c shows a similar abrupt transition.
37
-------
oleic liquid transport with constant water saturation (Weaver and Charbeneau, 1989) use this type of
function, and show this type of result.
In the case of monotonically decreasing X(S) (figure 8a), characteristics cross at some point A
Inside the z-t space (figure 8b). At point A, a discontinuity forms, since the characteristics each carry a
distinct saturation composition value and there cannot exist any point with multiple saturation
compositions. If such points existed, they would have, simultaneously, multiple water, oleic liquid and
air saturations. This situation is clearly meaningless. A discontinuity replaces the multiple-valued
saturation composition and is analogous to wave breaking. When the transition is abrupt (figure 8c),
the variation Is immediately resolved into a discontinuity. This situation corresponds to a complete
displacement of one composition by another, a condition which is solely determined by the X(S)
function. Plug flows are thus admitted as possible solutions of the hyperbolic system.
Discontinuities in the kinematic models at the leading edge of the infiltrating fluid are of the type
presented in figure 8c. Note, however, that a monotonically decreasing function is followed when the
transition on figure 7a is from S to S . Thus one X(S) function can generate both continuous and
discontinuous waves depending on the direction of the transition.
When the X(S) function has the form shown in figure 9a, the resulting wave displays both
continuous and discontinuous behavior. In figure 9b there is a portion of the wave which remains
continuous, corresponding to the Increasing portion of the X(S) function. When X(S) becomes a
decreasing function, the wave breaks as in figure 9b. For the abrupt transition (figure 9c), the
discontinuity originates on the boundary. Since the speed of the discontinuity equals the characteristic
velocity on one side, it is a contact discontinuity by definition.
The precise saturation composition where the wave breaks is determined by the particular
problem being solved. In some cases (presented below), no discontinuity forms because the Initial
condition is reached before the wave breaks. In other cases the Inttial condition helps determine the
saturation composition at the edge of the breaking wave.
38
-------
A
Figure 8 MonotonicaJly decreasing wave generating function. Figure 8a shows eigenvalues
which are monotonically decreasing with increasing S. Figure 8b shows a gradual
boundary transition (in time) from constant state, S . to constant state S which results
in overlapping characteristics (point A). Figure 8c shows a similar abrupt transition
which immediate resolution into a discontinuity.
39
-------
Figures Composite wave generating function. Figure 9a shows eigenvalues which first
Increase, then decrease with S. Figure 9b shows a gradual boundary transition (In
time) from constant state S1 to constant S^. The wave contains both a continuous and
a discontinuous portion. Figure 9c shows a similar abrupt transition.
40
-------
Rgure 10 shows the X eigenvalues as a function of water saturation for example 1, which is
described below. Injection at high water saturation follows a path with initially Increasing eigenvalues.
At some point this wave breaks, causing a contact discontinuity to form. Wave breaking Is associated
wtth the decreasing eigenvalues. By comparison wtth figure 9a, the solutions for the multiphase
Riemann problem are expected to consist of shocks and continuous waves which terminate In contact
discontinuities. Figure 11 shows the X eigenvalues, which behave similarly to the X eigenvalues.
The X eigenvalues, however, decrease very steeply in the vicinity of the maximum water saturation
(approximately 0.85). Most injections at high water saturation experience decreasing X eigenvalues
as the initial condition is approached, resulting in discontinuous paths.
SUMMARY OF APPROXIMATE GOVERNING EQUATIONS
For a Riemann problem, equation 20 along with initial data of the form
S(0,t)
S(z,0)
Is reduced to the following forms. For continuous saturation composition variation, solution of equation
37
r (2)
w
gives the saturation compositions across simple waves. The continuous wave speeds are given by
equation 24
41
-------
20.
15--
Eigenvalue 10--
5--
0-
0 0 0'.2 Qi4 0.6 0.8
Water Saturation
1 0
Figure 10 Example 1 X eigenvalues plotted as a function of water saturation.
20.
15--
Eigenvalue 10--
5--
0-
0.0 0.2 0.4 0.6 0.8 1 0
Water Saturation
Rgure 11 Example 1 X eigenvalues plotted as a function of water saturation.
42
-------
1
"x
dt
Since the characteristics within simple waves are straight lines, equation 24 has a simple analytical
solution. When saturation composition variations are discontinuous, solution of the non-linear
algebraic equation 40
wfrf'v^
w w
gives the saturation compositions across the discontinuities. The speeds of the discontinuities are
given by equation 26
dt
w2' w1
S -S
w2 w1
NUMERICAL SOLUTION METHODS
Most of the model equations are reduced to ordinary differential equations by the application of
the method of characteristics. Equations which do not have analytical solutions are solved by a Runge-
Kutta-Fehlberg method (RKF, Fehlberg, 1969). Fehlberg's methods have the advantage that they
contain automatic step size control based on a specified truncation error tolerance. For this work a
third order method is used for the solution and an embedded fourth order method is used for the step
size control (called RKF3(4)). Equation 37 is solved with the RKF3(4) solver for the saturation
composition variation across the simpie waves. Since the characteristic speeds are constant for
simple waves and the speeds of the discontinuities are also constant, equations 24 and 26 represent
straight lines and need not be solved numerically.
43
-------
Discontinuous paths across the saturation space diverge from continuous paths when the latter
are not straight (Helfferich 1970). Since equation 40 is a single non-linear algebraic equation, it is
solved by a method which combines bisection, inverse quadratic interpolation, and the secant method
(Press et al., 1986). The routine automatically chooses the most appropriate technique. The specific
usage of equation 40 is discussed below. Note however, that for breaking waves, the origin (i.e., the
point (S*, S*)), of a solution of equation (40) depends on the point at which the wave breaks. The
solutions are thus problem dependent so that a general solution cannot be mapped in advance. This
situation is in direct contrast to the continuous case where the saturation paths can be mapped for any
initial and injection conditions.
CONSTRUCTION OF THE SATURATION COMPOSITION SPACE FOR THREE-
PHASE FLOW
The input data for all the examples appears in appendix 5. In the first example, the oleic liquid's
mobility is approximately one-third that of the water, since Its density is less than that of water and its
viscosity is higher. The media is a silt loam soil, with average parameters from Brakensiek et al.
(1981). Since one potential application of the model is to the unsaturated zone, parameters
representing actual soils are used. Soils, in general, are expected to have wider pore size
distributions than aquifer materials, so the results presented here may differ somewhat from those
determined for aquifer materials. Because of the width of the pore size distribution, the relative
permeabilities drop rapidly as the saturation drops, due to the relatively high viscous dissipation in the
smaller pores. Significant flow of the fluids thus occurs when saturations are high. These conditions,
in fact, represent severe conditions for the mathematical model.
Figure 12 shows the solution of equation 37 for a number of initial saturation composition
vectors. The results are plotted on a triangular domain where each point represents a particular value
of S - S(Sw, SQ, Sfl). This representation of the solution is not unique to the Riemann problem, as any
three-phase flow solution could be plotted this way. The saturations are read off the graph by noting
44
-------
Riemann Problem Saturation Composition Space
Example 1
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0518
Matrix Properties
Lambda - 0.2099
Fluid Properties
Sa -= 100%
W-Density «=
O-Density
A-Density ~
W-Viscosity
O-Viscosity
A-Viscosity
,0000
.7000
,0012
1.0019
2.0039
0.0170
Sw - 100%
= U%
So - 100%
Figure 12 Example 1 saturation composition space showing all paths (variations in saturation
composition) which occur as solutions of the Riemann problem. The triangle vertices
are points of 100% saturation, while the edges are fines of zero saturation.
45
-------
that each side represents zero saturation of a fluid and the vertex opposite that side represents full
saturation. For example, the top vertex represents completely air-filled pore space, S - 1, and the
8
bottom of the triangle represents the two-phase water/oleic liquid system where no air is present in the
pore space, S - 0. The point at the middle of the triangle represents equal saturation of ad fluids,
8(0.3333, 0.3333, 0.3333).
The inner triangle bounds the region where all three fluid saturations are above residual. Air is
assumed to have a trapped (called residual) saturation during the injection processes considered here.
The trapped air saturation is assumed to be the air saturation which results in water relative
permeability of one-half, a value which comes from observation (Bouwer, 1966). In between the
boundaries, each path drawn represents part of a classical Riemann problem solution, as they are
determined from solution of equation 37. At each point within the saturation space, there are two
possible directions of saturation change, one associated with each eigenvalue. So the inner triangle is
filled with a grid representing all possible saturation composition variations which are solutions of
equation 37.
The calculations required for developing the entire saturation composition space typically
require less than 20 minutes of CPU time on a VAX 11/785 minicomputer. Solutions for specific
problems like those discussed below require much less CPU time as they do not require construction
of the entire saturation composition space.
The solutions for the X -waves originate at points wtth oleic liquid saturations just above
residual. Beginning at residual saturation is not possible, because the three-phase equations are
singular at SQ « SQr and the two-phase water-air equations have no component (see equation 37) which
points toward increasing oleic liquid saturation. The X -waves are determined via solution of the
water/oleic liquid equation 37. This equation is used instead of oleic liquid/air or water/air formulations,
because the paths end at the Sw - S^. edge. Thus the end point of the interval of integration is known.
The solution for the X -waves uses equation 37 and its water/air analogue to complete the paths. The
solutions begin at points with SQ« 0.15 and solved toward SB - S^ wtth the water/air equations, then
46
-------
solved toward S - S wtth the water/oieic liquid equation. The partial derivatives of the fractional flow
for the water/air system are found In a fashion similar to that presented in appendix 1. Solving across
the domain wtth different sets of equations (water/air, water/oteic liquid, oleic liquid/air) provides a
check of the solution. Here, all three sets produced identical paths.
QUALITATIVE NATURE OF THREE-PHASE FLOW
Before discussing specific Riemann problem solutions which can be constructed from figure 12,
the qualitative nature of three-phase flow in this domain will be discussed. These results are
representative of results obtained with other oleic liquids, soils and total fluxes. Paths resulting from
solution of equation 37 with the small eigenvalue, X , begin on the residual oleic liquid side of the
triangle and end on the residual water side. These paths are all bowed toward the residual air side of
the triangle, indicating that water/air and oleic liquid/air injections beginning on the left and right sides
of the triangle, respectively, tend to cause displacement of air from the pore space. The eigenvalues
associated with an injection at S - 0.849, S « 0.0501 = S , and S - 0.1009 are shown in figure 10.
' w o or a "
Since the eigenvalues decrease in magnitude after the water saturation decreases below
approximately 0.70, a X -wave wtth this path becomes discontinuous. Breaking of the wave results in
a contact discontinuity, which determines the exact saturations at the wave front.
Paths originating along the base of the inner triangle are associated with the larger eigenvalue,
X and follow three distinct patterns. First, those originating near the left vertex (S - 1) are nearly
fc W
parallel to the left edge of the triangle. In waves associated wtth such paths there is no appreciable
change In oleic liquid saturation even though all three phases flow. Although the oleic liquid Is flowing
and its saturation is above residual, there is unsteady flow in only the water and air phase along these
paths. Second, paths originating near the right vertex (S - 1) display similar behavior wtth regard to
the water phase, in that only the oleic liquid and air saturations vary. The third type of path is
Intermediate between the extremes and shows variation in all three saturations. These latter paths
47
-------
cover only a small region of the space. This behavior agrees with intuition as there is only a relatively
small region where all three relative permeabilities are significant (e.g., Leverett and Lewis, 1941).
CONSTRUCTION OF SATURATION PROFILES
Solutions are constructed from figure 12 by determining an injection condition-plateau-initjal
condition route on the diagram. This terminology follows Helfferich (1981) as portions of paths are
used to construct routes for each specific problem. Construction of the complete set of paths crossing
the saturation space is not necessary for a specific problem; only the route Is needed. The wave
associated with X (the smaller eigenvalue) is followed from the Injection condition toward the plateau.
This wave is called a X -wave or slow wave, and is followed first, because with smaller eigenvalues
this wave must be traversed before the faster wave (the X -wave). At the plateau, the solution
switches to the X -wave to complete the route to the initial condition. The route determines the specific
saturation compositions which exist in the solution. The occurrence in space and time of the
saturations along the route is determined by equation 24 on the continuous portion of the wave. When
discontinuities form, the route is adjusted as discussed below. The adjustment alters the saturation
composition at the plateau and thus the wave speeds. The speeds of the discontinuous portion of the
waves (shocks) are determined by equation 26. The continuous and discontinuous paths for this
example are shown in figure 13. When the routes are straight lines, the continuous and discontinuous
routes are identical, and only one X route exists. In this example the X routes are almost coincident.
48
-------
Riemann Problem Saturation Composition Space
Example 1 Saturation Routes
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0518
Matrix Properties
Lambda -= 0.2099
Fluid Properties
W-Density -=
O-Density =
A-Density -=
W-Viscosity
O-Viscosity
A-Viscosity
1.0000
0.7000
0.0012
- 1.0019
- 2.0039
- 0.0170
100%
njection
Condition
Sw - 100%
S =0%
a
So
100%
Figure 13 Example 1 saturation routes showing continuous solution of equation 37 and the
correction due to equation 40, which are nearly identical.
49
-------
TWO-PHASE INJECTIONS
Example 1 Oleic Liquid Bank Formation
In example 1, water and air are Injected into the soil at a saturation composition of S - S
(0.8490, 0.0501, 0.1009). Air te injected for Illustrative purposes hi these examples; later, injections
without air are discussed. The oleic liquid saturation at injection is slightly above residual (SQr -
0.0500), because the solution along the side of the inner triangle is degenerate. That is, the side is a
two-phase region, with no path which enters the three-phase region. Also, the three-phase equations
are singular along the edge (i.e. k - OatS »S in equation 53) and cannot be used. In this example
the total flux equals 2.00 m/d. The water, oleic liquid and air fluxes are determined from the fractional
flow equations (52 and 53), given the injection saturation composition. Selection of the injection
saturation composition is an important aspect of the problem solution. Allowing air to be above
residual serves two purposes. First, water may not fill the pore space during injection (rainfall). This is
particularly the case when air actually flows out of the domain (McWhorter, 1971). The second reason
for this type of injection is that it causes the X and X waves to separate from each other completely
so that there is no overlap. This condition is similar to that occurring in ion-exchange chromatography
solutions (Helfferich, 1970) which were studied in preparation for application to three-phase flow.
A depth-time plot of the solution, which is called the base characteristic plane, is shown in figure
14. Shown on the plot are the X - and X -waves. In this case, the X -wave has a continuous portion
which is illustrated by the fan-shaped characteristic pattern. This wave is a X -centered simple wave
which terminates in a contact discontinuity. The X -wave is discontinuous and is demonstrated below
to be a 2-shock. The plateau emerges from the origin and expands with time because of the
difference in speed between the contact discontinuity and the 2-shock. Ahead of the 2-shock is the
initial condition; at any time the location of the 2-shock corresponds to the maximum influence of the
injection.
50
-------
-24-
Depth
in Meters
-34-
-44-
-5-
Injection Condition
Initial
Condition
0.0
0.5 1.0 1.5
Time in Days
Xiwave
2.0
Contact Discontinuity
2-Shock
Characteristics
Figure 14 Example 1 base characteristic plane (depth-time plot of the solution). The solution
consists of a centered simple wave, contact discontinuity and 2-shock.
51
-------
Rgure 15 shows a depth-saturation profile for the solution at 24 hours after the beginning of the
Injection. Depth-saturation profiles compliment the base characteristic plane as they show the fluid
saturations at a given time in the solution. The water and total liquid saturation are plotted directly
against the depth. The oleic liquid saturation is the difference between the total liquid and water
saturations. Likewise, the air saturation is determined by subtracting the total liquid saturation from
one. Recall that at the boundary (depth - 0) the oleic liquid is very near its residual saturation, and
water and a small amount of air are injected. The general nature of this solution is that the injection at
the surface causes the oleic liquid to be displaced into a bank moving ahead of the water front (CEFQ
on figure 15). At this time, the water front is located at a depth of 0.73 meters (CD). The water
saturation at the front remains constant at 0.7821. The oleic liquid displacement is not complete, as
some of it is left behind the water front at saturations above residual (ABDE, the distance between A
and B is the oleic liquid injection condition which is dose to residual). The water saturation decreases
from the injection (A) to the water front due to the oleic liquid that is bypassed (AD). The smooth
variation in saturation above the water front corresponds to the X -centered simple wave on figure 14.
Most of the oleic liquid is displaced, however, into a bank moving ahead of the water front. In
this example the oleic liquid bank is bounded above by the contact discontinuity and below by the 2-
shock (figure 14). The bank corresponds to the plateau In the saturation space where S - 8(0.1000,
0.8034, 0.0966). The water saturation at the plateau equals the initial water saturation (S - 0.1000),
w
because the X -wave portion of the route is a straight line parallel to the water axis. The bank front
(FG) moves faster than the water front (CD), and so the bank expands with time. (Compare the
solution at 24 hours with that for 48 hours, as shown on figures 15 and 16.) The bank forms because
the flux of the oleic liquid is high enough so that it can move ahead of the incoming water. The leading
edge of the oleic liquid bank (FG) is discontinuous, because the X eigenvalues decrease from the
plateau to the initial condition. The details of how the discontinuous route is established follow.
When the X1 - wave breaks, equation 40 is used to determine the route from the water front to
the plateau. The saturation at the wave front, where the wave breaks, is the known (starred)
saturation composition in equation 40, and the discontinuous plateau saturation is unknown. The
saturation composition at the front is determined in the following way. By Iteration, a characteristic
52
-------
Bank-1 Data Set Profile at 1 Day
Depth
(m)
Oj
-1-
-3-
-4-
- 1
Injection Condition
1 ' > i i i i i
/\ R
/:,i
c
Plateau
F G
Initial Condition
1 1 1 1 i 1 1 1
0.0 0.1 0.2 0.3 0.4 0.5 0-6 0.7 0.8 0.9 1.0
Saturation
Water
Total Liquid
Figure 15 Example 1 saturation profile at 24 hours showing fluid saturations versus depth of
penetration. Tne oleic liquid Is being displaced into a bank (CEFQ) by the Incoming
water and air.
53
-------
Bank-1 Data Set Profile at 2 days
Injection Condition
V-
-1-
Depth -
(m)
-3-
-4-
1
i i 1 i 1 i r 1 / i
f \
1
Plateau
i
Initial Condition
1 1 1 1 i 1 1 1 -
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Saturation
Water
Total Liquid
Floure16
54
-------
speed is found that matches the speed of the discontinuity (Jump condition, equation 26) from the front
composition to the plateau (Initially estimated by the Intersection of the continuous waves). The
discontinuous route Is found by solving equation 40 with the trial frontaJ saturations. An outer iteration
is used to determine the value of the frontal saturations and plateau saturations which satisfy both the
Jump condition (equation 26) and equation 40 simultaneously. Typically two or three Iterations are
required to find the correct compositions.
For most of the saturation space, the paths from a plateau to an initial condition experience
decreasing eigenvalues. So the route from plateau to the initial condition is discontinuous and must
also be determined by solution of equation 40. In this case the starred saturation composition Is the
initial composition, and once again the plateau saturation is unknown. Solution of equation 40 for both
front to plateau and plateau to initial condition gives a new composition route across the saturation
space. The point where the two lines meet is the new plateau, and both equations are satisfied
simultaneously. In practice the plateaus diverge slightly from the continuous wave plateaus. The
slight divergence is enough to significantly affect the mass balance of the solution. Figure 13 shows
the continuous and the discontinuous routes, which are almost identical for this example.
The location in S-space of the intersection between the two waves is critical in determining the
character of the solution, if the intersection occurs after the X -wave breaks, a discontinuity forms in
the X -wave before the plateau is reached. The X -wave from the plateau to the Initial condition is a
breaking wave because the eigenvalues decrease in the direction of the initial condition. This
sequence describes the bank profile discussed above. When the plateau is reached before the X -
wave breaks, however, a second type of profile is produced. Here there Is a continuous wave solution
from the injection condition all the way to the plateau, and again a discontinuity from the plateau to the
Initial condition. The discontinuity in the oteic phase saturation is very slight in this case, and the
solution is characterized by bypassing of the oteic liquid by the water.
55
-------
Example 2 Qleic Liquid Bypassing
Rgure 17 shows the bypassing profile that is produced when the Initial condition is 8(0.5000,
0.1000, 0.4000), and the injection condition is the same as for example 1. In essence, the effective
conductivity erf oteic liquid present initially is Insufficient to match the speed of the incoming water. All
of the oteic Bquid is bypassed by the water. The base characteristic plane for exampte 2 is shown In
figure 18. The X -centered simple wave merges with the plateau before a discontinuity can form, so
there is no contact discontinuity and no oieic Bquid bank In this example. As in exampte 1, the X2-wave
is a 2-shock, the position of which determines the maximal influence of the Injection. There is a slight
discontinuity in oteic liquid saturation (0.1019 above vs. 0.1000 below) at this front (A). This small
Jump moves at the same speed as the water front.
Example 3 Oleic Liquid Bank Formation. Second Type
When the initial condition is nearer the center of the saturation space, the X path is curved, thus
all three fluid saturations in the bank differ from the initial condition. Figure 19 shows a resulting profile
with this initial condition. For such a displacement process the difference between the continuous and
discontinuous routes through the saturation space is significant (figure 20, contrast to figure 13).
The three types of solutions discussed above (bypassing, bank forming with the plateau water
saturation equal to the initial water saturation, and bank forming with different plateau and initial water
saturations) are related to each other by the locations of the X routes In the saturation composition
space (figure 21). At low initial air saturation, as used in the above examples, the character of the
solution changes gradually from bypassing to bank forming as the initial oteic liquid saturation
Increases along the line labeled A-A' on figure 21. When the air saturation is higher, as the olete Bqukj
saturation increases along line B-B', the zone where the X paths merge is relatively abrupt, and so is
the transition between the types of solutions, initial conditions tying near this zone represent problems
56
-------
Bypassing Data Set Profile at 1 Day
0.0-
-0.5--
-1.0- -
-1.5--
sr -**
-2.5--
-3.0--
-3.5--
-4.0-
Initial
Condition
Injection Condition
H 1
Plateau
0.0 0.2 0 4 0.6
Saturation
0.8
Water
Total Liquid
1 0
Figure 17 Example 2 bypassing profile showing ftuW saturation versus depth. The Incomina
water and air bypasses the otete HquW.
57
-------
Injection Condition
-14-
-24-
Depth
in Meters
-5-
-34-
-44-
0.0
Initial
Condition
015 110 1.5
Time in Days
2-Shock
Characteristics
2.0
Figure 18 Example 2 base characteristic plane (depth-time plot of the solution). The solution
consists of a centered simple wave, and 2-shock.
58
-------
Depth
(m)
-2--
-3--
-4--
-5-
Bank-2 Data Set Profile at 1 Day
Injection Condition
1 1 1
Plateau
Initial
Condition
o.o 0.2 0:4 o;e
Saturation
0.8
1 0
Water
Total Liquid
19 Example 3 bank profile with differing Initial and plateau water saturations.
59
-------
Riemann Problem Saturation Composition Space
Example 3 Saturation Routes
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0518
Sa - 100%
Matrix Properties
Lambda - 0.2099
Fluid Properties
W-Density - 1.0000
O-Density » 0.7000
A-Density » 0.0012
W-Viscosity -= 1.0019
O-Viscosity - 2.0000
A-Viscosity - 0.0170
Injection
Condition
Sw
100%
Sa = 0%
So
100%
Rgure 20 Example 3 saturation routes showing the continuous solution of equation 37 and Its
tfiscorrtinuous correction with equation 40.
60
-------
Riemann Problem Saturation Composition Space
Example 1
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sax - 0.0518
Sa - 100%
Matrix Properties
Lambda - 0.2099
Fluid Properties
W-Density - 1.0000
O-Density - 0.7000
A-Density - 0.0012
W-Viscosity - 1.0019
O-Viscosity 2.0039
A-Viscosity -
Sw - 100%
Sa =
So
100%
Figure 21 Profile transition from bypassing to bank forming. Une AA' crosses a region with a
smooth variation from bank to bypassing profiles. Une BB' crosses a region with an
abrupt variation from bank to bypassing proffles.
61
-------
where the method is essentially unstable, because the paths heading down toward the bottom of the
triangle diverge in many directions. A small change in the initial data causes a large change in the
solution.
Low mass conservation error in the examples demonstrates the existence of the solution to the
model equation 20. Since the governing equations are derived from mass conservation equations,
mass conservation is required in the solution. In other numerical methods, higher mass balance errors
are sometimes tolerated. Here, however, correct mass balance provides some assurance that the
solution has the correct plateaus, X -simple waves, shocks and contact discontinuities. These features
l\
must be programmed into the model, instead of arising purely from the solution. This difficulty is the
reason why so much emphasis is placed on the qualitative nature of the solutions in the previous
sections. Incorrect choices of solution features always lead to "solutions" which do not conserve
mass. So very low mass balance errors are required for these Riemann problem solutions. Table 1
shows the mass balance errors for examples 1 to 3. All of the errors are less than 0.2% and most are
one or two orders of magnitude lower. The highest mass balance error occurs for example 3, where
the plateau water saturation is different from the initial water saturation. In some regard this case is
the most difficult computationally.
Table 1 Mass Balance Errors for Examples 1 to 3
Example
1
2
3
Table 2 shows the eigenvalues and shock speeds for the major examples (1 to 3) presented in
this report. Examples 1 and 3 each have two shocks in their solutions: a contact discontinuity, and a 2-
shock. The contact discontinuity is "almost" a 1-shock. Example 2 has only one shock, which is a 2-
shock. The results presented in this table demonstrate that the discontinuities in the solutions obey the
appropriate shock inequalities (equations 28 to 32).
62
Water Phase
% error
-0.0019
-0.0056
-0.0700
Oleic Liquid
% error
-0.0358
0.0652
-0.1570
-------
Table 2 Shock Inequalities Satisfied by Examples 1 to 3
Example 1
a Incoming water front
Shock speed, s - 0.7014 m/d
X^)- 0.701 3
XjfSg)- 0.1453e-11
X2(S1)-21.94
X(S) - 14.24
Shock Inequalities satisfied: X (S ) < s < *«(SJ
Shock type: "almost" 1-shock, contact discontinuity
b. Oleic liquid bank front
Shock speed, s 1.5921 m/d
X1 (8^-0.14536-11
X^S). 0.27706-12
X (S )« 14.24
X2(S ) - 0.6706e-2
Shock Inequalities satisfied: ^0^2^ < 6
X/C \ ^ tt 1 /C '
.^w^J ^ o ^ o\ 4 i
Shock type: 2-shock
63
-------
Table 2 (Continued)
Example 2
Oleic liquid bypassing
Shock speed, s « 2.6684 m/d
X1 (8^.0.96486-1
X1(S2)-0.6072e-5
X^S^ - 38.31
X2(S2).0.1727e-3
Shock inequalities satisfied: X (S_) < s
Shock type: 2-shock
Example 3
a. Incoming water front
Shock speed, s « 2.8900 m/d
X^S^-2.890
X^Sg). 1.021
X^S^-29.58
l /c \ _ *ao *ap
o\ o' OC.wfc
Shock inequalities satisfied: X (S ) < s < X (S )
Shock type: "almost" 1-shock, contact discontinuity
64
-------
Table 2 (Continued)
Example 3 (continued)
b. Oteic liquid bank front
Shock speed, s - 3.8951 m/d
XJSJ-1.021
-32.32
- 0.6853e-1
Shock Inequalities satisfied: *o(S J < s
Shock type: 2-shock
65
-------
Example 4 Flow of TCE
When the mobility of the oteic Iqukl Is greater than water, the character of the results changes
sightly. For the same sofl as examples 1 to 3, but wtth TCE as the otete Iquld, figure 22 shows the
saturation composition space. Here the mobility of the TCE is almost three times that of the water. Wtth
higher oteic phase mobility, there is a greater tendency for otete Iquld banks to form as the diagram has
fewer bypassing paths (paths parallel to the S - S edge). More Initial conditions result In banks than
previously (examples 1 to 3), and only the lowest otete liquid Initial conditions result In bypassing profiles.
The oteic liquid, although now more mobile than water, still occupies the middle range of the pore size
distribution. Much of the qualitative character of the oteic Iquld flow Is due to this fact, largely
Independent of the oteic Kquld properties.
Example 5 Horizontal Flow
When the flow Is horizontal, some terms drop out of the fractional flow equations (appendix 3).
Given the same soil type, oteic Squid, and total flux as example 1, the loss of the gravity terms causes a
slight change in the saturation composition space, resulting in slightly changed saturation profiles.
(Compare figures 12 and 23). Comparison of the gravity and non-gravity forms of the fractional flow
equations (52 and 53 vs. 60 and 61) shows that the effect of gravity depends on the density differences,
oteic liquid relative permeability, and air relative permeability. For example, the saturation composition
space (figure 23) for horizontal flow te almost identical to that for vertical flow at low air saturation. Profiles
for both vertical flow (example 1) and horizontal flow are shown In figures 24 and 25. The only difference
in the scenarios is the absence of gravity in the horizontal flow case. At one day after the beginning of
the injection, the water and oteic liquid fronts are deeper when gravity is included. Obviously, gravity
adds a driving force to the vertical system that ts absent from the horizontal. Figure 25 shows the profiles
at the end of 2 days. Only very slight differences in saturation are seen In the continuous portions of the
profiles as the air saturations are low.
66
-------
Riemann Problem Saturation Composition Space
TCE Saturation Composition Space
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0518
Matrix Properties
Lambda - 0.2099
Fluid Properties
W-Density - 1.0000
O-Density » 1.4600
A-Density - 0.0012
W-Viscosity - 1.0019
O-Viscosity - 0.5699
A-Viscosity - 0.0170
Sa - 100%
Sw - 100%
Sa = 02
So - 100%
Figure 22 Example 4 TCE saturation composttkxi space showing, by comparison to figure 12, a
compressed zone of bypassing profiles (ABC).
67
-------
Riemann Problem Saturation Composition Space
Example 5 Horizontal Flow
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0518
Sa - 100%
Matrix Properties
Lambda - 0.2099
Fluid Properties
W-Density »
O-Density «
A-Density
W-Viscosity
O-Viscosity
A-Viscosity
1.0000
0.7000
0.0012
- 1.0019
- 2.0039
- 0.0170
Sw - 100%
Sa = 0%
So - 100%
Figure 23 Example 5 horizontal flow saturation composition space.
68
-------
Bank-1 Data Set Profile at 1 Day
-1-
Depth ,
(m) ~2'
-3-
-4-
' / \
A. A;/ j
' f
AA AA' «
t
B B';
BB BB '
i i J- .- i
0.0 0.2 0.4 0.6
Saturation
0.8
Water (Vertical)
Total Liquid (Veritcal)
1.0
Water (Horizontal)
Total Liquid (Horizontal)
Figure 24 Example 5 comparison of horizontal (A-A1 and B-B') and vertical (AA-AA1 and BB-BB')
bank profiles at 1 day.
69
-------
Bank-1 Data Set Profile at 2 Days
w-
-1-
Depth ,
(m)
-3-
-4-
J
/ \
*L _A:| r
1 :
AA AA1 {
i .
i
Bn' '
i
I i
IBB BE'
i
i
i i i i
0.0 0,2 0.4 0.6
Saturation
0.8
Water (Vertical)
Total Liquid (Veritcal)
1 0
Water (Horizontal)
Total Liquid (Horizontal)
Rgure 25 Example 5 comparison of horizontal (A-A1 and B-B') and vertical (AA-AA' and BB-BB'l
bank profiles at 2 days.
70
-------
Example 6 Oleic Liquid/Air Injection
The previous examples have focused on Injections wtth the oleic liquid at residual. When the
Injection consists of oleic liquid and air, wtth water at residual, the solution is constructed exactly as
described above, only wtth different initial and injection conditions. Returning to figure 12, injections from
the right edge (S - S ) can be seen to favor bypassing of the water by the oieic liquid over water bank
formation. The explanation for these phenomena follows. The eigenvalues (figure 10) Increase slowly as
the water saturation increases, reaching a maximum at S 0.85. Thus the contact discontinuity occurs
only after the oleic liquid saturation decreases to a tow value. As before, If the Initial condition lies on a X
path which intersects the X path before the discontinuity occurs, then the solution Is characterized by
bypassing. In this case, it is water that is bypassed, and the large number of paths parallel to tt>e S - S
W WrT
edge results in mostiy bypassing profiles. The water bank profiles result mostly from X paths which are
parallel to the S - S edge. Only a relatively few X paths of this type exist. The latter paths are
associated with very tow Initial oleic fluid saturations.
Injection at S(0.10, 0.80, 0.10) with an initial condition of 8(0.700, 0.075, 0.225) results In
displacement of the water into a water bank (figure 26). The oleic liquid saturation drops very rapidly in
the upper part of the slow wave (EH). The oleic liquid saturation at the front is tower than the water
saturation for a corresponding water Injection. This behavior Indicates the tendency for the water to
remain in place In the small pores. Only when water is in the large pores does the oleic liquid succeed in
displacing It.
71
-------
Oil Displacing Water
E Injection Condition F
0 0 0:2 014 0.6
Saturation
0.8
Water
Total Liquid
u.u-
-0.5-
-1.0-
-1.5-
Depth -2.0-
-2.5-
-3.0-
-3.5-
4.0.
: A :
Ate r
B G
.
Plateau
D C
Initial Condition
..... _L 1 1 1 _
1 0
Figure 26 Example 6 oteic fluW Injection showing displacement of water into a bank (ABCD) by
the incoming oteic liquid (EFQH) and air.
72
-------
PARAMETER VARIABILITY
Subsurface modeling is subject to uncertainty, as parameters may vary or be difficult to measure
precisely. In this section the effects of variation of saturated hydraulic conductivity and Brooks and
Corey's X are discussed.
By Increasing the saturated hydraulic conductivity by one standard deviation for the average silt
loam soil (Brakensiek et al., 1981), the saturation composition space is seen to shift (compare figures 27
and 12). The X paths show tess of a dip toward the S - S edge than figure 12. The previous solution of
I 8 ol
example 1 at two days fe compared to the solution with higher hydraulic conductivity In figure 28. As
expected, the water and oleic liquid fronts are deeper with Increased hydraulic conductivity. The
continuous saturation profiles reflect the variation In the saturation composition space. Slightly higher
water and oleic liquid saturations are encountered above the water front, since less mobile air Is
apparently retained in the profile as the hydraulic conductivity increases. In figure 29, the bypassing
profile of example 2 is compared with Its counterpart with higher hydraulic conductivity. Very slight
Increases In water saturation occur which result In a significantly higher front speed. This phenomena
implies a certain instability in the result as slight perturbations In saturation lead to large changes in front
position.
Increasing Brooks and Corey's X by one standard deviation for this soil type gives the saturation
composition space shown in figure 30. (Recall that a "X* wtth no subscript refers to a parameter of the
Brooks and Corey Model (equation), while a "X* with a subscript refers to an eigenvalue.) With a higher X
the pore size distribution Is more uniform, although In this case the distribution Is stilt relatively wide. The
apparent effect of the increase is to Increase the number of X paths wtth constant oleic liquid saturation.
More initial conditions tead to bypassing conditions than with tower Xs. Figure 31 shows the effect of the
change in X on the bank profile result from example 1. With higher X the water and oleic liquid fronts are
deeper into the profile, and the continuous water saturation profile shows a greater drop in saturation from
the surface to the water front, for example. The generally increased speeds are due to the fact that
given a saturation above residual, the relative permeability increases wtth X (equations 9 to 11). This
73
-------
Riemann Problem Saturation Composition Space
Hydraulic Conductivity 3.35e-3
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0518
100%
Matrix Properties
Lambda -= 0.2099
Fluid Properties
W-Density =
O-Density =
A-Density «=
W-Viscos ity
O-Viscosity
A-Viscosity
.0000
.7000
.0012
1.0019
2.0039
0.0170
Sw - 100%
Sa = 0%
So - 100%
Figure 27 Saturation composition space for hydraulic conductivity of 3.35e-3 cm/s
74
-------
Bank-1 Data Set Profile at 2 days
Depth
(m)
-i.
-2-
-3-
-4-
-5-
j
1 1 1 TTT r
A. A'/,' i;
1
AA AA1 I !
B B'
BB BB '
i '
i
i
_i .._. _i 1 1
0 0 0.2 0.4 0.6
Saturation
0.8
1.0
Water (he = 4.56e-4)
Total Liquid (he = 4.56e-4)
Water (he = 3.35e-3)
Total Liquid (he = 3.35e-3)
Figure 28 Effect of saturated hydraulic conductivity on bank profiles, showing Increase In depth
of water and oteic liquid fronts with increase In hydraulic conductivity (AA-AA1 vs. A-A'
and BB-BB1 vs. B-B').
75
-------
Bypassing Data Set Profile at 1 Day
0
-1--
Depth _.
(m)
-3--
-4
B
BB
0.0 0.2 Q'.4 0-6
Saturation
B'
BB',
0.8
1.0
Water (he = 4.56e-4)
Total Liquid (he = 4.56e-4)
Water (he = 3.35e-3)
Total Liquid (he = 3.35e-3)
Flaure29
76
-------
Riemann Problem Saturation Composition Space
Lambda 0.30
Residual Saturations
Swr - 0.0370
Sor - 0.0500
Sar - 0.0666
Matrix Properties
Lambda -= 0.3000
Fluid Properties
100%
W-Density -
O-Density -
A-Density -=
W-Viscosity
1.0000
0.7000
0.0012
- 1.0019
O-Viscosity - 2.0039
A-Viscosity
Sw
100%
sa - 0%
So
100%
Figure 30 Saturation composition space for Brooks and Corey's lambda - 0.30 showing a shift
away from the S «S axis (compare with figure 12).
o or
77
-------
Bank-1 Data Set Profile at 2 days
Depth
(m)
-1-
-2-
-3-
-4-
-
-f-
y
/ 1
A A' .'I 5
1 !>
AA AA1' ii
B B'
BB BB1
ill
0.0 0.2 0.4 0.6
Saturation
0.8
1.0
Water (lambda = 0.21)
Total Liquid (lambda = 0.21)
Water (lambda = 0.30)
Total Liquid (lambda = 0.30)
Figure 31 Effect of Brooks and Core/s lambda on bank profiles showing that Increased X leads
to deeper water (AA-AA1 vs. A-A') and otelc liquid fronts (BB-BB' vs B-R-\
(BB-BB1 vs. B-B1).
78
-------
effect is evident on the bypassing profiles shown in figure 32. The profiles show only a very slight change
In saturation as X increases, but with higher X much higher front speeds result.
SINGLE-PHASE WATER INJECTION
Injections at S - S(1 - S -S ,S , S J correspond to injections of water only into systems
or flu or GU
containing oleic liquid and air. This S corresponds to a single-phase flow condition. Thus relative to the
three-phase flow problem, this injection is doubly degenerate. The siow (X ) route leaving this point lies
along the bottom side of the inner triangle. This route must be followed before a fast (X ) route is
followed. The Implication of this fact is that until the plateau is reached the solution follows a degenerate
(two-phase) route. Any injection of water with both air and oleic liquid at residual results in saturation
profiles with no mobile air from the injection point to the water front.
The routes from the plateau to the initial condition take the solution into three-phase saturation
space. The X eigenvalues, however, are identically zero at S -S . Thus both the x -waves and X-
tL fi of 1 fc
waves begin with characteristic speeds Identical to zero, and the waves must overlap. (In the previous
examples, X < X at the plateau, so the waves separate).
Two additional problems are encountered in this transition to three-phase flow. First the X paths
have a minute bend toward the S - S axis, which not visible on any of the saturation composition
W WrT
triangles shown. The shape of this bend is similar to the visible bend In all the X paths. Numerical
Integration of equation 37 is very difficult across the bend.
The second problem may be related to either finite precision arithmetic or the loss of hyperbolidty
of the governing equations in part of the region of the S - S side of the triangle. Holden (1990)
G &T
discusses the possible loss of hyperbolidty for multiphase flow systems. The existence of real and
79
-------
Bypassing Data Set Profile at 1 Day
Depth
(m)
o
-1-
-2-
-3-
-4-
-5-
i
i
1
i i
i i
B B1
BB BB^j
i 1 1
t
f
I
I
j
1
i
f
t
t
t
\
<
i
i
i
*
0.0
0.2
Saturation
0.8
Water (lambda = 0.21)
Total Liquid (lambda =
Water (lambda = 0.30)
Total Liquid (lambda =
1.0
0.21)
0.30)
Floure
80
-------
distinct X1 and X^ eigenvalues, and thus hyperbolic character of the governing equations throughout the
rest of the domain is demonstrated by the solution of equation 37. ft the real eigenvalues did not exist, the
paths could not be drawn as shown in figures 12, 22, 27, and 30. Extension of the paths to the triangle
sides revealed truncation errors in a double precision version of the code. Calculation using 128 bit
arithmetic (VAX FORTRAN quad precision) allows the paths to be determined as shown. Further
increases in precision are not possible with RSKERL computing resources. Thus failure of the equations
to remain hyperbolic within the solution domain cannot be ruled out.
Since no assumptions are made about the type of the partial differential equations when deriving
the jump equation (Charbeneau, 1984), it holds regardless of the hyperbolictty of equation (20). The X
portion of the route is discontinuous, despite the initial increase of eigenvalues from x « 0 along the route
(figure 11), which otherwise would lead to a continuous wave solution. The continuous X -wave does not
exist because the x - and X -waves overlap. The X.-wave must have a continuous portion, since the
displacement of the oleic liquid by water cannot be complete. The solution for injection at S(1-S - S
S , S ) must be analogous to the two-phase Buckley-Leverett solution which requires dele liquid
bypassing.
The only difficulty remaining in profile construction relates to the location in S of the plateau.
Since discontinuities exist in the solution, the discontinuous route differs from the continuous route when
the continuous route is curved. Where the X routes are straight at their intersection with the S * S side,
^ 8 Eu
solutions can be reliably obtained by the methodology discussed above.
81
-------
SECTION 6
GENERALIZATION OF RIEMANN PROBLEM SOLUTIONS
Although solutions of Riemann problems provide much insight into fundamental fluid flow
phenomena, the restrictive nature of their boundary and initial conditions limits their usefulness, hi this
section, the limitations on the boundary conditions are removed for the two-phase flow of air and water
during infiltration. All of the principles involved also apply to three-phase flow problems. After
presenting the model equations, two major issues are discussed: first is the determination of the total
flux and second is the interaction of overlapping characteristic patterns.
MODEL EQUATIONS FOR THE INFILTRATION OF WATER SUBJECT TO AIR
PHASE RESISTANCE
The problem considered is the infiltration of water, subject to resistance due to the air phase
flow. The air is assumed to be displaced downward by the infiltrating water. Mass conservation for
water,
dS dq
_jw. =0 41
at dz
the incompressible flow condition,
and the saturation condition,
82
-------
give a single first-order hyperbolic equation to describe the system. The same procedure Is followed
as was applied to the three-phase flow problem presented previously In section 5. The fractional flow
equation and its first derivatives are given in appendix 4. After substituting In the fractional flow
function, the approximate governing equation is
as q at as
_ w. +2] _ w. _ W..O 44
at TI as az
w
The continuous wave solution of the single hyperbolic equation 44 is
dS dr 3S
w. + w. * 0 45
at dt az
along paths defined by
46
eft TI as
w
The characteristic lines defined by equation 46 are straight lines onty tf the total flux, a , does not vary
with time. For this problem a is allowed to vary with time, so curved characteristics must be tracked.
GENERALIZED TOTAL FLUX
One reason why Riemann problem solutions become attractive is that there often can be
found a single point tn time and space where the total flux is known. Usually this is an injection
condition, as was the case in the previous section and in the classic BucWey-Leverett solution, where
83
-------
constant flux boundary conditions are applied. When the injection is not constant, then difficulty is
encountered in determining the total flux. One time when the total flux cannot be considered constant
te when a constant flux condition abruptly ends. After this occurs there is zero flux of the injected fluid
at the boundary, but the total flux dearly is not zero throughout the domain. A second case occurs
when there is a constant rate rainfall at the boundary. The amount of water that can enter the profile is
limited by the infiltration capacity of the soil, which is a function of the saturated hydraulic conductivity,
relative permeability, antecedent water saturation, cumulative infiltration, and other factors (Bear,
1972). Thus the water flux entering the profile varies with time, even when the precipitation rate Is
constant. An expression for the water flux can be used to determine the total flux, assuming that at
the surface the water flux equals the total flux. Although this is recognized as being only an
assumption, it is expected that in most cases little error is introduced into the water flow.
The total flux is determined by summing and then integrating Darc/s law equations for water
and air. This procedure is taken from Morel-Seytoux (1973). Summing the Darcy's law equations
gives the total flux, for downward flow
kk
H
w
-dP
+ P Q
dz Kwy
kk
ra
47
By adding an identity and rearranging, the following is obtained
w
-dP
kk
ra
dP
dz
48
which is further manipulated to give
w
dP
_
dz
k
w ra
k + p. k
a rw 'w ra
'dP
49
This results in an expression which can be integrated for the total flux, a,
84
-------
f2(kra'
where
v -V-
1 ra' rw k(u k + u. k )
a rw W ra
i k
w
k + jj. k
w ra a r>v
Each inl&gra! ap-pearing in equation 50 can be evaluated numerically at any time necessary in the
soJution, since they depend only on the saturation profile. The total flux so found applies
instantaneousJy to the who^le profile, because of equation 16. The integrals are evaluated numericalK
by a seven-point Gaussian quadrature routine (Abramowttz and Stegun, 19-651
OVERLAPPING CHARACTERISTIC PATTERNS
Each change in boundary condition of a quasi-Pi near system generates a set of waves. Since
there is one equation solved for this problem, each boundary condition change creates one wave. The
derivative of the fractional flow function determines the smoothness of the solution, as tn a sense, this
derivative is the eigenvalue of the system. For this problem the fractional flow derivative has a shape
sJmftar to that sho-wn on figure 10. Recaing the discussion in section 5, each wave consists of a
continuous portion which terminates in a contact discoiTtinurtv,
85
-------
Rgure 33 shows the base characteristic plane for a problem which consists of a finite duration
rainfall of constant precipitation rate which both begins and ends with abrupt boundary transitions.
The beginning of the rainfall creates the first wave. This wave is centered at the beginning of the
rainfall. There is a transition from the high boundary saturation to the relatively tow initial saturation.
Characteristics originate from this point. The characteristic associated with the highest possible water
saturation has zero velocity, and thus remains at the ground surface. Each successive characteristic
associated with lower saturations has a faster velocity than the previous characteristic. At some point,
however, the wave "breaks' as the characteristics cross. After that point, decreasing the saturation
decreases the speed of the characteristic. As before, without a discontinuity, the solution would
suggest that one point (z,t) is associated with two values of water saturation. Clearly this te not
physically realistic and the classical method of characteristics is supplemented by a generalized
solution for the contact discontinuity.
Figure 34 shows the saturation profile at 0.2 days after the beginning of the rainfall. Initially,
the profile contained a uniform water saturation of 0.51. By 0.2 days, the water front has progressed
60 cm into the soil. The water saturation decreases slightly from its maximum of 0.9328 at the surface
to 0.9208 at the contact discontinuity (front). Physically, the immobility of the maximum water
saturation suggests that water cannot completely displace the mobile air. If the water entered the
profile at its maximum saturation, all mobile air would be displaced. As It is, some of the mobile air
above residual remains behind the water front. Flow Into the profile can only take place when a small
amount of mobile air remains behind the water front. The speed of the characteristics is directly
proportional to the magnitude of the total flux (equation 46), so when the total flux Is constant, as
during the initial portion of the rainfall, the characteristics are straight lines. When the total flux is
decreased, due to decreasing infiltration capacity, the characteristics begin to curve toward the time
axis. After this time, equation (34) must be solved numerically for the position of each characteristic.
The additional numerical integration adds significantly to the computational burden of the solution.
If the water is supplied at such a rate that the water saturation at the surface is tess than the
maximum (S^ « 1. - S^J, then the characteristic associated with the boundary has a nonzero speed.
Thus the characteristic moves into the subsurface. A plateau region emerges behind this
86
-------
Base Characteristic Plane
Silt Loam Soil
6 Hour Duration Rainfall
0.0
depth
(meters)
-1 . 0--
-1 .4--
Rgure 33 Water-air Jnffftrstion bass characteristic plane showing depth of fronts and characters
wtth time. Fronts are created when the rafrrfafl begins at 0.0 days and when tt ends at
0.25 davs.
-------
DEPTH
(meters)
TIME =0.2 DAY
Injection
w * w-
-0.2-
-0.4-
0 fi
\J w*
-0.8-
-1.0-
-1.2-
-1.4-
C
Initial
Condition
i 1
A B
-
-
-
-
-
i t
0.0 0.2 0.4 OJ6 0^8
WATER SATURATION (%)
1.0
Figure 34 Water-air infiltration profile at 0.2 Days showing rainwater (ABC) and slight decrease
In water saturation from the injection condition to the front (CB).
88
-------
characteristic. The water saturation in the plateau must remain constant at the boundary value,
because there is only enough water supplied to sustain this value.
When the rainfall ends, there is a second abrupt change in the boundary condition. The
change results in the formation of a second wave. This wave is governed by the same fractional flow
derivative as the first wave. Note that the transition from tow saturation to high saturation follows the
same general shape of the function as did the transition from high to low saturation. Thus the
boundary saturation, which In this case is taken as residual, te again associated with an immobile
characteristic. Increasing the water saturation gradually increases the speed of each successive
characteristic, so that the smooth, fan-shaped characteristic pattern begins to emerge. Once again a
saturation is reached where the characteristic speed begins to decrease. Since the saturation cannot
be multiple valued, a discontinuity must also appear in the second wave. Note that a discontinuity
must also appear between the first and second waves, because characteristics from the second wave
would touch characteristics from the first wave which are associated with different water saturations.
The profile at 0.3 days (figure 35) shows the second wave originating at the ground surface
and terminating in a shock which separates the two waves at 12 cm Into the profile. The water
saturation at the surface is residual and the saturation increases rapidly to 0.8520 just behind the front.
The front is a contact discontinuity since its speed matches the characteristic speed just above the
front. Ahead of the second wave lies the first which shows the same slight decrease in saturation as it
did before the rainfall ended.
As indicated by figure 36, the second wave's discontinuity rapidly overtakes the original (first
wave) front. In reality, the second wave's front will be strongly smoothed by both gravity and capillary
phenomena. The sharpness of the front should be viewed as an artifact of the model, but is essential
for the conservation of mass in the solution.
89
-------
TIME =0.3 DAY
DEPTH
(meters)
\J . \J-
0.2-
0.4-
0.6-
0.8-
1.0-
1.2-
1.4-
1 1 1 «-o
E \
D\
C
A
Initial
Condition
1 1
B
1 1
0.0 0.2 0.4 0.6 0.8
WATER SATURATION (%)
1 0
Figure 35 Water-air infiltration profile at 0.3 Days showing the second shock (DC) created when
the injection drops to zero (E). ' "waiea wnen
90
-------
TIME =5.0 DAYS
0.0-
Injection Condition
DEPTH
(meters)
0.2--
-0.4--
-0.6--
-0.8--
-1.0--
-1.2--
-1.4--
D
Initial
Condition
00 0.2 0.4 0.6 0.8
WATER SATURATION (%)
1.0
Figure 36 Water-air infiltration after the original front has been overtaken by the second (AD).
91
-------
SECTION 7
DISCUSSION OF RESULTS
The application of characteristic methods to one-dimensional flow problems reveals the
fundamental behavior of three-phase systems. Although the full governing equations are parabolic, the
approximate hyperbolic equations which are solved hi this report describe a fundamental portion of the
flow phenomena. The efficiency that is achieved by simplifying the governing equations allows results to
be developed which apply to any pair of injection and initial conditions for a given oleic phase and soil
type (e.g., figure 12). The semi-analytic nature of the solution is primarily responsible for the generation of
these "universal" results. The following conclusions are drawn from the work:
For systems with constant injection conditions and uniform initial conditions (Riemann problems),
flow occurs in two regimes. When the injection consists mostly of water, mobile oleic liquids (organic
Hquids or so-called NAPLs) are displaced into 'banks' or are bypassed. The banks move ahead of, and
are driven by, the incoming water. The bypassing profile is characterized by the water moving past the
oleic liquid, without causing H to be displaced into the profile. The occurrence of the regimes is
determined by the oleic liquid properties, media properties, and the initial amount of all fluids present.
Although the type of boundary condition used in these solutions is not what is expected in the
field, mathematically, It causes the hyperbolic problem to be a Riemann problem, for which there exists a
large body of precise mathematical results. The flow regimes which exist in these solutions are believed
to be fundamental for those expected from the solution of the hyperbolic system with general initial and
boundary conditions.
The dependency of the results on the initial condition suggests that for cases where the oleic
liquid has drained to a low saturation in the upper part of the soil profile, incoming rainfall will likely bypass
the oleic liquid. Near the surface, this condition is likely to occur when a rainfall follows the end of an oleic
liquid release by some hours. The oleic liquid has enough time to begin draining from the surface, so the
incoming water experiences oleic liquid saturations which are low. Thus the bypassing regime is likely to
be favored. Conversely, water injections at relatively high oleic liquid saturations are dominated by the
92
-------
bank formation regime. Following the previous discussion, when the rainfall begins infiltrating there is
likely to be bypassing near the surface. As the water front moves deeper into the profile, however, a
saturation composition may be reached where the oleic liquid forms a bank. Thus the event may be
characterized by both regimes. Verification of this behavior awaits the development of the genera) model.
These composite flow regimes point out the necessity of full understanding of the three-phase Riemann
problem before attempting to solve general problems.
For systems with injections of mostly oleic liquid, the dominant flow regime is water bypassing.
This conclusion confirms the intuition that the oleic liquid, which occupies a mid-range of the pore space
due to its wetting behavior, does not displace water from the small pores. Preferential wetting causes the
water to be too tightly held to the media to be easily displaced by a non-wetting liquid. Thus spills of oleic
liquid to soil profiles, which contain water at the so-called field capacity, are likely to bypass the pre-
existing water.
In all cases, the behavior of the oleic liquid is dependent on its properties (density and viscosity)
relative to those of water. As the oleic liquid becomes more mobile, water injections favor the bank
formation regime somewhat, because the oleic liquid Is of high enough mobility to match the speed of the
Incoming water. The character of the results, however, always reflects preferential wetting as indicated in
the conclusion stated above. Thus in the example presented, the bypassing regime was only slightly
larger for the higher mobility TCE case than the lower mobility oil, because of preferential wetting.
The conclusions slated lead to the following recommendations:
The three-phase solutions should be extended to simulate conditions with more general boundary
and initial conditions. The methodology presented for two-phase flow simulation appears to be adaptable
for this purpose. Models could then be created for discrete pulses of water which represent rainfalls
preceding and following oleic liquid releases. For each liquid entering the profile, air phase resistance
could be included by following the procedures developed in section 6. Alternatively, the kinematic
formulation could be used when the flow is gravity driven. The next step fe to develop useful models for
field problems which Include interphase partitioning phenomena. Previous work on kinematic models
demonstrates that such phenomena are amenable for inclusion into this type of model. Multi-dimensional
extension of the models should be considered, since one-dimensional flow is a restrictive condition, which
93
-------
doesn't apply to heterogeneous field sites. The one-dimensional models can be useful, however, for
spills occurring over large areas or for screening applications.
Laboratory validation of the existence of sharp fronts would be valuable for demonstrating the
usefulness of the model. However, absolutely sharp fronts are not necessary for the model, because the
sharp fronts can be shown theoretically to move at the same speed as the true fronts. The relative
permeability equations used in the models were derived for drainage conditions. Experimental work
would indicate the importance of hysteresis in the constitutive relations governing the flow and the
necessity of including imbibition relative permeabilities.
The laboratory work should include investigation of the saturation conditions at the ground
surface. Somewhat arbitrary conditions are used in the work discussed in this report. A related issue is
that the Injections used here force flow to be unidirectional. Experimental work indicates countercurrent
flow is common, so the model should also be adapted for this condition.
94
-------
REFERENCES
Abramowttz, M., and I. A. Stegun, Handbook of Mathematical Functions. Dover, New York, 1965.
Be31". ->-. Dynamics of Fluids in porous Media American Elsevier, New York, 1972.
Bouwer, H., Rapid field measurements of air entry value and hydraulic conductivity of soil as
significant parameter in flow system analysis, Water Resources Research. 2, 729-738, 1966.
Brakensiek, D. L, R. L Engleman, and W. J. Rawls, Variation within texture dasses of soil water
parameters, Transactions of the American Society of Agricultural Engineers. 335-339, 1981.
Brooks, R. H. and A. T. Corey, Hydraulic Properties of Porous Media. Colorado State University
Hydrology Paper No. 3. Ft. Collins, Colorado, 1964.
Buckley, S. E. and M. C. Leverett, Mechanism of fluid displacement in sands, Transactions American
Institute of Mining Engineers. 146, 107-116, 1942.
Burdine, N. T., Relative permeability calculations from pore size distribution data, Transactions
American Institute of Mining Engineers. 198, 71-78, 1953.
Charbeneau, R. J., Mutticomponent exchange and subsurface solute transport: Characteristics,
coherence, and the Riemann problem, Water Resources Research. 24(1), 57-64, 1988.
Charbeneau, R. J., Kinematic models for soil moisture and solute transport, Water Resources
Research. 20, 699-706, 1984.
Charbeneau, R. J., J. W. Weaver, and V. J. Smith, Kinematic Modeling of Multiphase Solute Transport
In the Vadose Zone. United States Environmental Protection Agency,Report, EPA/600/2-89/035,
June 1989.
Dougherty, E. L and J. W. Sheldon, The use of fluid interfaces to predict the behavior of oil recovery
process, Society of Petroleum Engineers Journal. 4(2), 171-182, 1964.
Fehlberg, E., Low-order Classical Runoe-Kutta Formulas With Stepslze Control and Their Application
to Some Heat Transfer Problems. NASATR R-315,1969.
Qlimm, J., B. Undquist, O. McBryan, and L Padmanabhan, A front tracking reservoir simulator, five-
spot validation studies and the water coning problem, in The Mathematics of Reservoir
Simulation. R.E. Ewing, editor, Society of Industrial and Applied Mathematics. Philadelphia,
Pennsylvania, 107-136,1983.
Green, W. H. and G. A. Ampt, Studies on soil physics, Journal of Agricultural Science. 4, 1-24,1911.
Helfferich, F. G., Theory of mutticomponent, multiphase displacement in porous media, Society of
Petroleum Engineers Journal, 21, 51-62, 1981.
95
-------
Helfferich, F. G., Multic»mponent wave propagation: Attainment of coherence from arbitrary starting
conditions, Chemical Engineering Communications. 44, 275-285, 1986.
Hetfferich, F. G., and G. Klein, Mutticomponent Chromatooraphv. Theory of Interference. Marcel
Dekker, New York, 1970.
Holden, L, On the strict hyperbolicity of the Buckley-Leverett equations for three-phase flow In a
porous medium, SIAM Journal of Applied Mathematics. 50(3), 667-682,1990.
Jeffrey, A., and T. Taniuti, Non-Linear Wave Propagation. Academic Press, New York, 1964.
Lax, P. D., Hyperbolic systems of conservation laws II, in Communications on Pure and Applied
Mathematics. Interscience Publishers, New York, Vol. 10,1957.
Lax, P. D., Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves.
Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania, 1973.
Leverett, M. C., and W. B. Lewis, Steady flow of gas-oil-water mixtures through unconsolidated sands,
Transactions. American Institute of Mining Engineers. 142,107-152,1941.
Kreyszig, E., Advanced Engineering Mathematics. 3rd edition, Wiley, New York, 1972.
McWhorter, D. B., Infiltration Affected bv Flow of Air. Colorado State Hydrology Paper No. 49, Ft.
Collins, Colorado, 1971.
Morel-Seytoux, H. J., Two-Phase Flows in Porous Media, In Advances In Water Resources. Vol. 9, V.
T. Chow, editor, Academic Press, New York, 119-202,1973.
Mull, R., The migration of oil-products in the subsoil with regard to ground-water pollution by oil,
Proceedings of the Fifth International Conference on Advances In Water Pollution Research. San
Francisco and Hawaii, S. A. Jenkins, editor, Pergammon, Oxford, Vol. 2. HA 7(A)/1-8, 1970.
Peaceman, D. W., Fundamentals of Numerical Reservoir Simulation. Elsevier Scientific, Amsterdam,
1977.
Pope, G. A., L. W. Lake, and F. G. Helfferich, Cation exchange In chemical flooding: Part 1 - Basic
theory without dispersion, Society of Petroleum Engineers Journal. 18(6), 418-434,1978.
Press, W. H, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterting, Numerical Redoes The Art of
Scientific Computing. Cambridge University Press, 1986.
Rawls, W. J., D. L Brakensiek, and N. Miller, Green-Ampt infiltration parameters from soils data,
ASCE Journal of Hydraulic Engineering. 109, 62-70, 1983.
Richards, L A., Capillary conduction of liquids through porous mediums, Phvsies. 1, 318-333,1931.
Rieble, D. D., T. H. Illangasekare, D. V. Doshi, and M. E. Malhiet, Infiltration of immiscible
contaminants in the unsaturated zone, Ground Water. 28(5), 685-692,1990.
96
-------
Rozdestvenskii, B. L and N. N. Janenko, Systems of Quasllinear Equations and Their Applications to
Gas Dynamics. American Mathematical Society, Providence, Rhode Island, 1980.
Schwille, F.. Dense Chlorinated Solvents in Porous and Fractured Media. Lewis, Chelsea, Michigan,
1988.
Smtth, R. E., Approximate soil water movement by kinematic characteristics, Soil Science Society of
America Journal. 47, 3-8, 1983.
Smoller, J., Shock Waves and Reaction-Diffusion Equations. Springer-Verlag, New York, 1983.
Stone, H. L., Probability model for estimating three-phase relative permeability and residual oil data,
Journal of petroleum Technology. 20(2), 214-218, 1970.
Weaver, J. W., A Critical Review of Immiscible Qroundwater Pollutant Transport. Masters Thesis, The
University of Texas at Austin, 1984.
Weaver, J. W., Kinematic Modeling of Murtiphase Subsurface Transport. Doctoral Dissertation, The
University of Texas at Austin, August, 1988.
Weaver, J. W., and R. J. Charbeneau, A kinematic model of oil drainage in soils, submitted for
publication in Water Resources Research. 1989.
Wyllie, M. R. J. and G. H. F. Gardner, The generalized Kozeny-Carman equation: tts application to
problems of multiphase flow in porous media, World Oil. 146, 210-227, 1958.
Young, D. M., and R. T. Gregory, A Survey of Numerical Mathematics. Volume II, Dover, New York,
1988.
97
-------
APPENDIX 1
DERIVATION OF THE THREE-PHASE FRACTIONAL FLOW EQUATIONS
The three-phase fractional flow functions are derived as follows for the water-oleic liquid
system. Two independent expressions for the fractional flow exist and are derived as follows.
Beginning with the flux equation for the oleic liquid,
-k k
ro
dP
o
+ p g sina
I dz °
51
where a is the angle between the flow direction and the horizontal. When flow is downward sina is
equal to -1, and equation 51 is the same as equation 4. Using equation 51 and noting the capillary
pressure relationships
P -P -P
o w cwo
dP dP dP
o cwo w
dz
dz dz
the flux equation for water can be used to eliminate the derivative of the oleic liquid pressure In
equation 51
-k k
ro
dP q n
cwo ^w w
- (p.-
dz
k k
rw
Since this model neglects the contribution of the capillary pressure gradient to the flux the first
relationship for the oleic liquid fractional flow is
98
-------
-k k
o1 * o
ro
gsina + _^J!f
w
52
rw o
The second of the relationships is derived from the air phase flux equation
-k k
ra
9P
a
Idz
+ P g sina
0
and capillary pressure relationships,
P - P - P
a o coa
dP 3P dP
a o coa
dz dz
Using the oleic liquid flux equation to eliminate the air phase pressure derivative yields,
-k k
ra
a
dP q n
coa ^o o
- (P
dz k k
ro
which, when the capillary gradient is neglected, is conveniently rewritten as
k k
1 -
o2~ o
M-
ro a
53
The two functions can be written in compact notation as
99
-------
f .Ak + Bf k /k
o1 ro w ro rw
54
and
o2
<1-fw+Ckra)kro
(k + Dk )
v ro ra'
55
wtth
w
K -K
SW |1 SO
sina
w
K -K
so sa n
a )
sina
Equations 52 and 53 comprise a system of two linear equations for the two unknowns, f (f «f «f )
andf .
w
The partial derivatives of the fractional flows, which are required for the mass conservation
equations 18 and 19, are likewise given by a system of two linear equations In two unknowns. The
derivatives with respect to water are,
.
o1
9f . dk
o1 r o
dS dk dS
w row
df . dk
o1 r w
dk dS
rw w
df . df
o1 w
df dS
w w
56
100
-------
dfo2dkro d'o2dkra \2 d< w
+ + 57
row raw w w
While the fractional flow derivatives with respect to the oil saturation are,
df . df . dk df . df
ol 01 ro o1 w
58
dS dk dS df dS
o ro o wo
df . df . dk df . dk df _ df
o2 o2 ro o2 ra o2 w
59
dS dk dS dk dS df dS
o ro o ra o wo
The partial derivatives of the fractional flows with respect to the independent variables appearing In
equations 56 through 59 are (from equations 54 and 55)
df . -Bf k
o1 w ro
dk k 2
rw rw
df . f
o1 w
A + B~;
-.1 k
dk r w
ro
df . Bk
o1 ro
df " k
w rw
101
-------
9f Ck
o2 ro
D(1 - f + Ck )k
w ra ro
dk (k +Dk ) (k + Dk )'
ra ro ra' * ro ra'
df _ (1 - f + Ck ) (1 - f + Ck )k
o2 v w ra w ra ro
3k (k + Dk ) (k + Dk )
ro v ro r a' v ro r a7
df
o2
-k
r o
df (k + Dk )
w v ro ra'
To complete the evaluation of the partial derivatives, the relative permeability (fraction of saturated
conductivity) functions are needed. For this report, the functions are based on a Brooks and Corey
Integration of the Burdine equations and have the form,
rw
S S
w wr
1 -S
wr
ro
S - S
o or
1 -S
or
S +S - S
o w wr
1 - S
wr
e-2
S - S
w wr
1 -S
wr
e-2
k -
ra
S - S
a ar
1 -S
ar
1 -
S +S - S
o w wr
1 - S
wr
e-2
For simplicity, only drainage functions are considered here. The required partial derivatives of these
functions are given by,
102
-------
dk
r w
dS
w
1 S
wr.
S - S
w wr
1 -S
wr
e-1
dk
ro
dS
E - 2
1 -S
wr
s s
o or
1 -S
or )
S +S - S
o w wr
1 - S
wr
e-3
(1 -S )
or'
S +S - S
o w wr
1 - S
wr
£-2
S S
w wr
1 -S
wr
e-2
dk
r o
w
£ - 2
1 -S
1 wrJ
S - S
o or
1 S
I or J
S +S - S
o w wr
1 - S
1 wr J
E-3
S - S
w wr
1 -S
1 wr J
E-3
dk
ra
1 - S - S - S
o w ar
1 -S
I ar J
2
1 - S
ar
tL
r 2- E
1 -S
1 wr
1 -S -S -S
o w wr
1 -S
1 ar J
rs +s - s
o w wr
1 - S
i wr J
1 -
E- a
S +S - S
o w wr
1 - S
1 wr J
E-2
dk dk
ra r a
dS dS
o w
Given a saturation composition, the evaluation of equations 56 through 59 is straightforward, since
they also comprise a two by two system. Once all of the fractional flow derivatives are known for a
given saturation composition, the eigenvalues and eigenvectors of matrix A(S) appearing in equation
20 can be determined and the method of characteristics solutions can be found.
103
-------
At times, water/air and oleic fluid/air formulations are convenient for solving the approximate
governing equations by the method of characteristics. The equations for these systems can be found
by switching appropriate subscripts in equations 52 and 53. The partial derivatives of the resulting
equations are found in a manner similar to that presented above.
104
-------
APPENDIX 2
KINEMATIC MODEL SOLUTION
The solution of the model equation 20 for kinematic conditions has been presented by
Charbeneau et al. (1989), and Weaver (1988). The solutions obtained In the previous work are
developed from ad hoc considerations. In this appendix, the classical method of characteristics
technique is used to derive the kinematic model solution. The two solutions are shown to be identical.
M all fluids flow at fluxes less than the effective conductivities given a saturation composition, then the
flow is driven by gravity only. No pressure gradient term appears In the Darcy's law equations and the
approximate governing equations have the form
as
dt
as
IT- -o
dZ
where the matrix, A*(S), is
A*(S)
and the components are
dK
a
11 dS
w
dK
"21 dS
w
105
-------
a -
22 dS
o
The zero appears in matrix A*, because the water relative permeability in a strongly water wet system
Is a function of water saturation only. There is no derivative of water relative permeability with respect
to oil saturation. The two eigenvalues of matrix A* are equal to the diagonal entries (e.g., Kreyszig,
1972), so A-1 - a11 and \2 « a^.
The left eigenvectors of the matrix are defined as
One set of vectors satisfying this relationship is
',-[1,0]
Equation 23, i k DS = 0, can be used to determine the classical method of characteristics solution of
the kinematic model.
For the first eigenvalue, X , i - i and equation 23 becomes
.,(.)
Idt dz J Ut dz
which reduces to
106
-------
dS ]
a _JK -0
Idt az )
The last result can be rewritten in terms of the total derivative,
0 68
Dt
along characteristics defined by dz/dt = X - (1 h\) dK /dS . Note that the solution for X is the same as
I iWr inr I
the solution for the water phase when no oteic fluid is present.
For the second eigenvalue, X , t - 1 and equation 23 becomes
2 K 2
fas as 1 ._.fas as } .
L_Jtt + a^_jK+I (2) IB + a^_fl - 0
lat az ) lat az J
«a_ *«**_?!+t*«**jy.o
22"811 dZ
which is equal to
DS -a_ fas dS
Dt a22~an dt dz
along characteristics defined by dz/dt = X2= (1/n) aKro/aSo- The last equation Is simplified by adding an
Identity which, upon rearrangement, gives
Dt a-an at az az
107
-------
The sum of the first two terms inside the parentheses is identically zero, since these terms are the
mass conservation equation for water,
DS 95
2 - -a^ w. 69
Dt dz
The solution given by equations 68 and 69 is the same as the solution reported previously by
Charbeneau et al. (1989).
108
-------
APPENDIX 3
HORIZONTAL FLOW
When the flow system is horizontal, the coefficients A and C appearing in equations 54 and 55
are equal to zero. The fractional flows become
f 4 - B f k / k 60
o1 w ro rw
and
ro
f - 61
(k + Dk )
ro ra
The required partial derivatives of the horizontal fractional flow equations are, (from equations 60 and
61)
df , -Bf k
o1 w ro
dk k 2
rw rw
62
df . f
o1 w
BT~ 63
dk r w
ro
df . Bk
o1 ro
. 64
df k
w rw
109
-------
df DM - f )k
o2 w' ro
m . 65
dk ( k + Dk )
ra ro ra
dfo2
dkro
df _ -k
o2 r o
66
67
df (k + Dk )
w ro ra
The effect of gravity depends on the magnitude of coefficients A and C, which in turn depend on the
density difference between the fluids and on the oleic fluid and/or air relative permeability. The latter
dependency indicates that the soil properties play a significant role in determining the importance of
the gravity term. When the oleic fluid and air relative permeabilities are low, the effect of gravity is
minimized. Both expressions for the fractional flow, equations 60 and 61, depend on the gravity effect,
as do some of the partial derivatives (equations 63, 65, and 66).
110
-------
APPENDIX 4
DERIVATION OF THE TWO-PHASE FRACTIONAL FLOW EQUATIONS
The two-phase fractional flow functions are derived as follows for the water-air system.
Beginning wtth the flux equation for water,
-kk
rw
w
w
Uz
sina
70
and noting the capillary pressure relationships
P -P -P
a w cwa
and
dP dP dP
a cwa w.
dZ * dZ dZ
the flux equation for air can be used to eliminate the derivative of the water liquid pressure in equation
70
V
-k k
rw
^w
dP q \i
cwa ^a i
. dz kkn
- " (P^, + PQ)9 sina
W a
3
Since this model neglects the contribution of the capillary pressure gradient to the flux, the relationship
for the water fractional flow is
111
-------
k k
_Dtf_
f * ra w
w
kk
(P -
k._ K.. q.H "a
1 +
ra w
which can be written as
Ak k + Bk
rw rp r w
w* k + Bk
ra rw
where A «= (K - K >i /M^/cr, and B = ft /n . The partial derivative of fractional flow with respect to
the water saturation is needed for equation 46, and is given by
df df dk df dk
w. w. r w ffi r a
dS " dk dS + dk dS
w r w w raw
The necessary partial derivatives are
df Ak + B B(Ak k +Bk )
w. ra v rw ra r_w_
dk * k + Bk ., _, .2
r w ra rw (k + Bk )
v ra iw
df
Ak
ra
(Ak k +Bk )
rw ra rw
dk " k + Bk " , _, ,2
r a ra rw k + Bk j
v ra PAT
The relative permeabilities for the two-phase system are
k -
rw
S - S
w wr
1 -S
wr
112
-------
ra
s -s ]
, w wr
'- 1-s
wrJ
2
1 -
S - S I
w wr
1 -S
1 wr J
e-2'
The required partial derivatives of these functions are given by
dk
r w
dS
w
c
1 - S
wrJ
S - S
w wr
1 -S
wr
e-1
dk
dS
w
S 1
wr
S -S
w wr
1 -S
wr.
1 -
S - S
w wr
1 -S
wr
e-2
1 -
S -S
w wr
1 -S
wr.
2- e
1 -S
wr
1 -
S - S
w wr
1 -S
wr
e-21
113
-------
APPENDIX 5
INPUT DATA SETS
Units for the following input data sets are as follows:
Saturated hydraulic conductivity
Viscosity
Density
Total Flux
Beginning and ending time
Angle with the horizontal
The remaining input parameters are dimensionless:
Brooks and Corey's lambda (X)
Porosity
Residual saturations (all fluids)
Saturations (all fluids)
The "calculated parameters" have the following dimensions:
centimeters/second
centipoise
grams/cubic centimeter
meters/day
days
degrees
Brooks and Corey's epsilon (e)
Residual air saturaion
Saturated Conductivities (ail fluids)
dimensionless
dimensionless
meters/day
114
-------
Example 1 Oleic Liquid Bank Formation
MULTIPHASE COHERENCE CALCULATIONS
A**********************-*,
RIEMANK PROBLEM DATA SET SILT LOAM FROM BRAKENSIEK
BROOKS AND COREYS LAMBDA 0.2100
SATURATED HYDRAULIC CONDUCTIVITY 0.4560E-03
POROSITY 0.4840
RESIDUAL WATER SATURATION 0.3700E-01
RESIDUAL OIL SATURATION 0.5000E-01
WATER VISCOSITY
OIL VISCOSITY
AIR VISCOSITY
1.002
2.004
0.1800E-01
WATER DENSITY
OIL DENSITY
AIR DENSITY
1.000
0.7000
0.1205E-02
TOTAL FLUX
ANGLE WITH HORIZONTAL
BEGINNING TIME
ENDING TIME
2.000
-90.00
O.OOOOE+00
24.00
INITIAL CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.1000
0.5000
0.4000
INJECTION CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
CALCULATED PARAMETERS
BROOKS AND COREYS EPSILON
RESIDUAL AIR SATURATION
SATURATED WATER CONDUCTIVITY
SATURATED OIL CONDUCTIVITY
SATURATED AIR CONDUCTIVITY
0.8490
0.5010E-01
0.1009
12.52
0.5185E-01
0.3940
0.1379
0.2643E-01
115
-------
Example 2 Oleic Liquid Bypassing
MULTIPHASE COHERENCE CALCULATIONS
***********************************'
RIEMANN PROBLEM DATA SET SILT LOAM FROM BRAKENSIEK
BROOKS AND COREYS LAMBDA 0.2100
SATURATED HYDRAULIC CONDUCTIVITY 0.4560E-03
POROSITY 0.4840
RESIDUAL WATER SATURATION 0.3700E-01
RESIDUAL OIL SATURATION 0.5000E-01
WATER VISCOSITY
OIL VISCOSITY
AIR VISCOSITY
1.002
2.000
0.1800E-01
WATER DENSITY
OIL DENSITY
AIR DENSITY
1.000
0.7000
0.1205E-02
TOTAL FLUX
ANGLE WITH HORIZONTAL
BEGINNING TIME
ENDING TIME
5.000
-90.00
O.OOOOE+00
24.00
INITIAL CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.4000
0.7500E-01
0.5250
INJECTION CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.8490
0.5010E-01
0.1009
CALCULATED PARAMETERS
BROOKS AND COREYS EPSILON
RESIDUAL AIR SATURATION
SATURATED WATER CONDUCTIVITY
SATURATED OIL CONDUCTIVITY
SATURATED AIR CONDUCTIVITY
12.52
0.5185E-01
0.3940
0.1382
0.2643E-01
116
-------
Example 3 Oleic Liquid Bank Formation, Second Type
MULTIPHASE COHERENCE CALCULATIONS
RIEMANN PROBLEM DATA SET SILT LOAM FROM BRAKENSIEK
BROOKS AND COREYS LAMBDA 0.2100
SATURATED HYDRAULIC CONDUCTIVITY 0.4560E-03
POROSITY 0.4840
RESIDUAL WATER SATURATION 0.3700E-01
RESIDUAL OIL SATURATION 0.5000E-01
WATER VISCOSITY
OIL VISCOSITY
AIR VISCOSITY
1.002
2.004
0.1800E-01
WATER DENSITY
OIL DENSITY
AIR DENSITY
1.000
0.7000
0.1205E-02
TOTAL FLUX
ANGLE WITH HORIZONTAL
BEGINNING TIME
ENDING TIME
2.000
-90.00
O.OOOOE+00
24.00
INITIAL CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.6000
0.2000
0.2000
INJECTION CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.8490
0.5010E-01
0.1009
CALCULATED PARAMETERS
BROOKS AND COREYS EPSILON
RESIDUAL AIR SATURATION
SATURATED WATER CONDUCTIVITY
SATURATED OIL CONDUCTIVITY
SATURATED AIR CONDUCTIVITY
12.52
0.5185E-01
0.3940
0.1379
0.2643E-01
117
-------
Example 4 Flow of TCE
MULTIPHASE COHERENCE CALCULATIONS
RIEMANN PROBLEM DATA SET SILT LOAM FROM BRAKENSIEK
BROOKS AND COREYS LAMBDA 0.2100
SATURATED HYDRAULIC CONDUCTIVITY 0.4560E-03
POROSITY 0.4840
RESIDUAL WATER SATURATION 0.3700E-01
RESIDUAL OIL SATURATION 0.5000E-01
WATER VISCOSITY
OIL VISCOSITY
AIR VISCOSITY
1.002
0.5700
0.1800E-01
WATER DENSITY
OIL DENSITY
AIR DENSITY
1.000
1.460
0.1205E-02
TOTAL FLUX
BEGINNING TIME
ENDING TIME
2.000
O.OOOOE+00
24.00
INITIAL CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.1000
0.5000
0.4000
INJECTION CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
CALCULATED PARAMETERS
BROOKS AND COREYS EPSILON
RESIDUAL AIR SATURATION
SATURATED WATER CONDUCTIVITY
SATURATED OIL CONDUCTIVITY
SATURATED AIR CONDUCTIVITY
0.8490
0.5010E-01
0.1009
12.52
0.5185E-01
0.3940
1.011
0.2643E-01
118
-------
Example 5 Horizontal Flow
MULTIPHASE COHERENCE CALCULATIONS
**********!
RIEMANN PROBLEM DATA SET SILT LOAM FROM BRAKENSIEK
BROOKS AND COREYS LAMBDA 0.2100
SATURATED HYDRAULIC CONDUCTIVITY 0.4560E-03
POROSITY 0.4840
RESIDUAL WATER SATURATION 0.3700E-01
RESIDUAL OIL SATURATION 0.5000E-01
WATER VISCOSITY
OIL VISCOSITY
AIR VISCOSITY
1.002
2.004
0.1800E-01
WATER DENSITY
OIL DENSITY
AIR DENSITY
1.000
0.7000
0.1205E-02
TOTAL FLUX
ANGLE WITH HORIZONTAL
BEGINNING TIME
ENDING TIME
2.000
O.OOOOE+00
O.OOOOE+00
24.00
INITIAL CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.1000
0.5000
0.4000
INJECTION CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.8490
0.5010E-01
0.1009
CALCULATED PARAMETERS
BROOKS AND COREYS EPSILON
RESIDUAL AIR SATURATION
SATURATED WATER CONDUCTIVITY
SATURATED OIL CONDUCTIVITY
SATURATED AIR CONDUCTIVITY
12.52
0.5185E-01
0.3940
0.1379
0.2643E-01
119
-------
Example 6 Oleic Liquid/Air Injection
MULTIPHASE COHERENCE CALCULATIONS
t*********************t
RIEMANN PROBLEM DATA SET SILT LOAM FROM BRAKENSIEK
BROOKS AND COREYS LAMBDA 0.2100
SATURATED HYDRAULIC CONDUCTIVITY 0.4560E-03
POROSITY 0.4840
RESIDUAL WATER SATURATION 0.3700E-01
RESIDUAL OIL SATURATION 0.5000E-01
WATER VISCOSITY
OIL VISCOSITY
AIR VISCOSITY
1.002
2.004
0.1800E-01
WATER DENSITY
OIL DENSITY
AIR DENSITY
1.000
0.7000
0.1205E-02
TOTAL FLUX
ANGLE WITH HORIZONTAL
BEGINNING TIME
ENDING TIME
2.000
-90.00
O.OOOOE+00
24.00
INITIAL CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.7000
0.7500E-01
0.2250
INJECTION CONDITION
WATER SATURATION
OIL SATURATION
AIR SATURATION
0.1000
0.8000
0.1000
CALCULATED PARAMETERS
BROOKS AND COREYS EPSILON
RESIDUAL AIR SATURATION
SATURATED WATER CONDUCTIVITY
SATURATED OIL CONDUCTIVITY
SATURATED AIR CONDUCTIVITY
12.52
0.5185E-01
0.3940
0.1379
0.2643E-01
120
------- |