United States
Environmental Protection
Agency
Robert S. Kerr Environmental
Research Laboratory
Ada OK 74820
Research and Development
EPA/600/S2-91/015 Sept. 1991
EPA Project Summary
Approximate Multiphase Flow
Modeling by Characteristic
Methods
James W. Weaver
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 asso-
ciated with releases of such liquids Is
transport as a water-lmmlsclble liquid
phase. Conventional finite element and
finite difference models of multiphase
flow may have extreme data and com-
putational requirements. Sites with In-
sufficient data or modeling for screening
purposes may warrant using a simpli-
fied approach. In this document, ap-
proximate models of Immiscible flow
are presented for two- and three-phase
flow. The approximations are con-
structed by representing the flow by
hyperbolic equations which have
method of characteristics solutions.
This approximation has the additional
benefit of being based on the funda-
mental wave behavior of the flow, which
is revealed by the solutions of the mod-
els. An Important result Is that for three-
phase flow, two flow regimes exist. The
first Is characterized by the displace-
ment 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 re-
gimes Is dependent on the organic liq-
uid properties, soil type, and the Initial
amounts of the fluids present. Two-
phase flows consisting of pulse appli-
cations of water result in overlapping
simple waves and contact discon-
tinuities. These models form the basis
for future extension of the three-phase
model to include pulse boundary condi-
tions.
This Project Summary was devel-
oped by EPA's Robert S. Kerr Environ-
mental Research Laboratory, Ada, OK,
to announce key findings of the research
project that Is fully documented In a
separate report of the same title (see
Project Report ordering Information at
back).
Introduction
The focus of the research described in
the full report is the transport of oily liquids
within the subsurface. Subsurface contami-
nation by oily wastes is recognized as a
ubiquitous problem which presents some
of the most challenging remediation diffi-
culties. 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 hy-
drophobic 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 vary-
ing intensity, duration, and inter-storm ar-
rival time. Varying water fluxes play a
central role in multiphase subsurface flow,
as the rainfalls provide a driving force for
the hydro log ic inputs to the system. A fo-
cus of the research is the effect of tran-
sient flows in the water, organic, and air
phases.
Printed on Recycled Paper
-------
Since the governing equations for this
problem are a system of coupled non-
linear partial differential equations, numeri-
cal 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 approxima-
tion-error free. During the initial phase of
site investigation, or for alternatives screen-
ing, 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 is to represent the flow by a
system of hyperbolic equations, which de-
scribe the wave behavior of the flow sys-
tem. 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 simu-
lators. Most important, however, is the in-
sight they provide into the fundamental
character of the flow.
As stated above, the phases modeled
are water, a water-immiscible organic liq-
uid, and air. The term oleic liquid is used to
indicate the organic liquid, usually such a
liquid is called a NAPL (Non-Aqueous
Phase Liquid). The word oleic is used be-
cause it means "oily" and carries the con-
notation of immiscibility with water. The
terminology and the models presented
herein are not restricted to petroleum based
oleic liquids but also include organic sol-
vents and other liquids.
The contents of the report are as fol-
lows: Section 2 contains a summary of the
conclusions and recommendations. Sec-
tion 3 describes the background and moti-
vation 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 equa-
tions is presented. The solution of a "Ri-
emann 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 illus-
trative examples. Since Riemann problems
have restrictive boundary and initial condi-
tions, the generalization to pulse boundary
conditions is discussed in Section 6. Appli-
cation of this theory is made to the infiltra-
tion of water subject to air phase resistance.
A discussion of the results is presented in
Section 7. These sections are summarized
in the foltowing pages.
Research Scope
The specific motivation for this work is
that previous attempts to model multiphase
flow of contaminants by simplified meth-
ods 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 mod-
els allow more realistic fluid distributions
but are limited to gravity driven flows. As a
result, an intense rainfall, high viscosity
oleic phase, and/or low permeability soil
may cause the model to break down. 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.
Two types of solutions are presented.
First is the solution to a classical problem
of mathematical physics, called a Riemann
problem. Mathematically, a Riemann prob-
lem is defined as a hyperbolic system of
equations with constant injection and uni-
form initial conditions. Although this is of
limited direct applicability to field problems,
the solution of the Riemann problem is
important for several reasons: Mathemati-
cally, there has been a great deal of work
done on these solutions for systems other
than those presented here. In the 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
are presented which illustrate several typi-
cal flow regimes for three-phase flows.
Third, one-dimensional solutions of the Ri-
ernann problem form the basis for multi-
dimensional front tracking models of
petroleum reservoirs which may be adapt-
able for subsurface contaminant transport.
The second type of solution is for sys-
tems which can be represented by injec-
tion conditions which consist of discrete
pulses. The model system for this applica-
tion is the infiltration of water, subject to air
phase resistance. The solution of this prob-
lem consists of an extension and generali-
zation of the Riemann problem solution
techniques. The solution which is presented
is the first step in developing a model for a
surface spill of an oleic contaminant, fol-
lowed or preceded by a sequence of rain-
fall events.
Model Formulation
The three-phase flow model is based
on the assumption that the medium is uni-
form and the flow is one-dimensional. Each
phase has a residual or trapped saturation
which remains constant. The relative
permeabilities are described by equations
that are based on the Burdine, and Brooks
and Corey theories. All phase transport
properties are assumed to remain con-
stant. The model describes immiscible
transport without interphase partitioning
phenomena, or any sources/sinks in the
domain of interest.
As indicated in Appendix 1 of the full
report, the model neglects the gradients of
the capillary pressure. One of the main
effects of the capillary gradient is to smooth
sharp fronts. 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. Flow visualization experiments con-
ducted at the Robert S. Kerr Environmen-
tal Research Laboratory and elsewhere
suggest that infiltrating fronts indeed re-
main sharp for a variety of situations. This
assumption is critical for eliminating the
parabolic character of the phase conser-
vation equations, so that the method of
characteristics can be applied to this prob-
lem. Application of the method of charac-
teristics 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.
Two important caveats must be noted.
First, a necessary condition for the exist-
ence of classical solutions is that the sys-
tem remains hyperbolic throughout the
solution domain. The condition for
hyperbolicrty is that the eigenvalues re-
main real and distinct throughout the solu-
tion domain. For the multiphase flow
problem considered here, this is demon-
strated in Section 5. 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 discon-
tinuities, the partial derivatives appearing
in the conservation equations do not exist,
so neither do the classical solutions. In
order to overcome this difficulty, 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 com-
plete solution, which consists of regions
with smooth saturation variation separated
by sharp fronts.
Discontinuous solutions are notorious
for being non-unique, and abundant ex-
-------
amptes of this behavior exist. Generalized
entropy conditions (shock inequalities) are
used to select physically realistic discon-
tinuous solutions. Two types of discontinu-
ous solutions are important for the report,
k-shocks and contact discontinuities. By
obeying the shock inequalities, a discon-
tinuous solution is assured to be physically
realistic and therefore the proper discon-
tinuous solution. The shocks appearing in
the solutions presented are shown to be
proper k-shocks. The second type of shock
is a contact discontinuity, for which the
shock speed matches the characteristic
speed on one side of the shock.
Numerical Solution Methods
Most of the model equations are re-
duced to ordinary differential equations by
the application of the method of character-
istics. Equations which do not have ana-
lytical solutions are solved by a Runge-
Kutta-Fehlberg method. Fehlberg's meth-
ods 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.
Discontinuous paths across the saturation
space diverge from continuous paths when
the latter are not straight. The equation
governing this behavior is a single non-
linear algebraic equation; it is solved by a
method which combines bisection, inverse
quadratic interpolation, and the secant
method. The routine automatically chooses
the most appropriate technique.
Construction of Saturation
Profiles
Solutions are constructed in two basic
steps: First, an injection condition-plateau-
initial condition route on the saturation com-
position space diagram (Figure 1) is
determined. As explained in the report the
solutions consist of two waves which may
be continuous, continuous and terminate
in a contact discontinuity, or discontinuous
(k-shock). The wave associated with XI
(rr\e 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 V
wave). At the plateau, the solution switches
to the Xj-wave to complete the route to the
initial condition. The route determines the
specific saturation compositions which ex-
ist in the solution. When discontinuities
form, the route is adjusted as needed. The
adjustment alters the saturation composi-
tion at the plateau and thus the wave
Riemann Problem Saturation Composition Space
Example 1 Saturation Routes
Residual Saturations
$„ - 0.0370
Sw » 0.0500
S.,« 0.05/5
Matrix Properties
Lambda - 0.2099
Fluid Properties
W-Density « 1.0000
O-Density M 0.7000
A-Density » 0.00/2
W-Viscosity* 1.0019
O-Viscosity « 2.0039
A-Viscosity*0.0170
100%
Sw m 100%
Sa'0%
S0 = /00%
Figure 1. Example 1 saturation routes showing the continuous solution and the correction
due to the discontinuities. In this case, the routes are nearly identical and not
distinguishable on the figure.
speeds. When the routes are straight lines,
the continuous and discontinuous routes
are identical, and only one Xj route exists.
In this example, the X, routes are almost
coincident. The second step is to deter-
mine the occurrence in space and time of
the saturations along the route.
Example 1 Olelc Liquid Bank
Formation
In Example 1, water and air are in-
jected into the soil with a total flux of 2.0 m/
d, and a saturation composition of S - S
(Sw, S0, S.) - S (0.8490, 0.0501, 0.1009),
where Sw, S0, and S. are the percent of the
pore space (saturation) filled by water, oleic
liquid and air, respectively. The oleic liquid
saturation at injection is slightly above re-
sidual (S,,, - 0.0500), because the solution
along the side of the inner triangle is de-
generate. That is, the side is a two-phase
region, with no path which enters the three-
phase region. Also, the three-phase equa-
tions are singular along the edge and
cannot be used.
A depth-time plot of the solution, which
is called the base characteristic plane, is
shown in Figure 2. Shown on the plot are
the X, - and Xa-wayes. 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 discon-
tinuity. The Vwave is discontinuous and
is 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 condi-
tion; the location of the 2-shock corre-
sponds to the maximum influence of the
injection.
Figure 3 shows a depth-saturation pro-
file for the solution at 24 hours after the
beginning of the injection. Depth-satura-
tion profiles complement the base charac-
teristic 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. Like-
wise, 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 dis-
-------
Injection Condftion
-1 ••
-2 ••
Depth
(m)
-3 ••
-4 ••
0.0
Initial
Condition
\1 Wave
0.5 1.0 1.5
Time in Days
Contact Discontinuity
2-Shock
Characteristics
2.0
Figure 2. Example 1 base characteristic plane (depth-time plot of the solution). The solution
consists of a centered simple wave, contact discontinuity, and 2-shock.
placed into a bank moving ahead of the
water front (CEFG on Figure 3). 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 com-
plete, 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 close 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,-cen-
tered simple wave on Figure 2.
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 2). The bank corresponds to the
plateau in the saturation space where S -
S(0.1000,0.8034, 0.0966). The water satu-
ration at the plateau equals the initial water
saturation (Sw - 0.1000), because the V
wave portion of the route is a straight line
parallel to the water axis. The bank front
(FG on Figure 3) moves faster than the
water front (CD), and so the bank expands
with time. 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,, eig-
envalues decrease from the plateau to the
initial condition.
Example 2 Oleic Liquid Bypassing
Figure 4 shows the bypassing profile
that is produced when the initial condition
is 8(0.5000, 0.1000, 0.4000), and the in-
jection condition is the same as for Ex-
ample 1. In essence, the effective
conductivity of oleic liquid present initially
is insufficient to match the speed of the
incoming water. All of the oleic liquid Is
bypassed by the water. The base charac-
teristic plane for Example 2 is shown in
Figure 5. The X,-centered simple wave
merges with the plateau before a disconti-
nuity can form, so there is no contact
discontinuity and no oleic liquid bank in
this example. As in Example 1, the Xj-
wave is a 2-shock, the position of which
determines the maximal influence of the
injection. There is a slight discontinuity in
oleic 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.
More examples are presented in the
report which illustrate the effects of oleic
liquid properties, media properties (hydrau-
lic conductivity and pore size distribution),
gravity and the injection of the oleic liquid.
In addition to developing the solution for
one injection-initial condition pair, the equa-
tions can be solved for the entire suite of
possible conditions. The report illustrates
several full complete saturation composi-
tion spaces, showing all possible varia-
tions in saturations which are solutions of
the classical Riemann problem. From these,
general conclusions can be drawn con-
cerning the displacement processes.
Generalization of Riemann
Problem Solutions
Although solutions of Riemann prob-
lems provide much insight into fundamen-
tal fluid flow phenomena, the restrictive
nature of their boundary and initial condi-
tions limits their usefulness. In Section 6 of
the report, 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. Two major is-
sues are discussed: first is the determina-
tion of the total flux and second is the
interaction of overlapping characteristic pat-
terns.
One reason why Riemann problem so-
lutions become attractive is that there of-
ten can be found a single point in time and
space where the total flux is known. Usu-
ally this is an injection condition, as was
the case for the examples presented and
in the classic Buckley-Leverett solution,
where constant flux boundary conditions
are applied. When the injection is not con-
stant, then difficulty is encountered in de-
termining the total flux; i.e., the total flux
cannot be considered constant when a
constant flux condition abruptly ends. After
this occurs, there is zero flux of the in-
jected fluid at the boundary; but the total
flux clearly is not zero throughout the do-
main. 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 satu-
rated hydraulic conductivity, relative per-
meability, antecedent water saturation,
cumulative infiltration, and other factors.
Thus, the water flux entering the profile
varies with time, even when the preciprta-
-------
Injection Condition
-1 --
-2--
-3--
E
Plateau
Initial
Condition
0.0
0.5
Saturation
Water
Total Liquid
1.0
Figure 3.
Example 1 saturation profile at 24 hours showing fluid saturation versus depth of
penetration. The oleic liquid is being displaced into a bank (CEFG) by the incoming
water and air.
Injection Condition
A 7 Wave
0.5 1.0 1.5
Time in Days
2.0
2-Shock
Chat acteri sties
Figure 4.
Example 2 base characteristic plane (depth-time plot of the solution). The solution
consists of a centered simple wave, and 2-shock.
tion 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 as-
sumption, it is expected that in most cases
little error is introduced into the water flow.
Integrating the Darcy's law equations re-
sults in an expression for the total flux
which can be evaluated at any time in the
solution. Because of the time variation in
the total flux, the characteristics which pre-
viously were straight lines now curve. Thus,
much additional complexity is introduced
into the numerical implementation of the
solution.
Discussion of Results
The application of characteristic meth-
ods to one-dimensional flow problems re-
veals the fundamental behavior of
three-phase systems. Although the full gov-
erning equations are parabolic, the ap-
proximate hyperbolic equations which are
solved in this report describe a fundamen-
tal 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. The semi-ana-
lytic nature of the solution is primarily re-
sponsible for the generation of these
"universal" results. The following conclu-
sions are drawn from the work.
For systems with constant injection con-
ditions and uniform initial conditions (Ri-
emann problems), flow occurs in two
regimes. When the injection consists mostly
of water, mobile oleic liquids (organic liq-
uids 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 it 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 hy-
perbolic problem to be a Riemann prob-
lem, for which there exists a large body of
precise mathematical results. The flow re-
gimes which exist in these solutions are
believed to be fundamental for those ex-
pected 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 tow
-------
saturation in the upper part of tha soil
profile, incoming rainfall will likely bypass
the oleic liquid. Near the surface, this con-
dition is likely to occur when a rainfall
follows the end of an oleic liquid release by
several hours. The oleic liquid has enough
time to begin draining from the surface, so
the incoming water experiences oleic liq-
uid saturations which are low. Thus the
bypassing regime is likely to be favored.
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 dis-
place water from the small pores. Prefer-
ential wetting causes the water to be too
tightly held to the media to be easily dis-
placed 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 (den-
sity and viscosity) relative to those of wa-
ter. 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, al-
ways reflects preferential wetting as indi-
cated 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.
Injection Condition
-1 --
Depth
-2 --
.3 --
Plateau
Initial
Condition
0.0
0.5
Saturation
Water
Total Liquid
1.0
Figure 5 Example 2 bypassing profile showing fluid saturation versus depth. The incoming
water and air. bypasses the oleic liquid.
.S. GOVERNMENT PRINTING OFFICE: 1992 - 648-OHO/40224
-------
-------
James W. Weaver (also the Project Officer) is with Robert S. Kerr Environmental
Research Laboratory, Ada, OK 74820.
The complete report, entitled "Approximate Multiphase Flow Modeling by Characteristic
Methods" (Order No. PB91-190959/AS; Cost: $23.00, subject to change) will be
available only from:
National Technical Information Service
5285 Port Royal Road
Springfield, VA 22161
Telephone: 703-487-4650
The EPA Project Officer can be contacted at:
Robert S. Kerr Environmental Research Laboratory
U. S. Environmental Protection Agency
Ada, OK 74820
United States
Environmental Protection
Agency
Center for Environmental Research
Information
Cincinnati, OH 45268
BULK RATE
POSTAGE & FEES PAID
EPA PERMIT NO. G-35
Official Business
Penalty for Private Use $300
EPA/600/S2-91/015
------- |