United States
Environmental Protection
Agency
Robert S. Kerr Environmental
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/S2-917042 Sept. 1991
v/EPA Project Summary
Modeling Multiphase Organic
Chemical Transport in Soils and
Ground Water
J.C. Parker, A.K. Katyal, J.J. Kaluarachchi, R.J. Lenhard, T.J. Johnson,
K. Jayaraman, K. Unlu and J.L. Zhu
Ground-water contamination due to
surface spills or subsurface leakage of
hydrocarbon fuels, organic solvents
and other immiscible organic liquids is
a widespread problem which poses a
serious threat to ground-water re-
sources. In order to model the move-
ment of such materials in the subsur-
face, it is necessary, in general, to con-
sider flow of water, nonaqueous phase
liquid (NAPL) and air, and transport of
individual chemical components, which
may move by convection and disper-
sion in each phase. A mathematical
model was developed for multiphase
flow and multlcomponent transport in
porous media with water, NAPL and air
or any subset of these phases. Numeri-
cal procedures for solving the system
of coupled flow equations, based on
various formulations of the governing
equations, were compared. Accurate
representation of three-phase perme-
ability-saturation-capillary pressure (k-
S-P) relations is crucial to model
multiphase fluid movement and accu-
rate models for interphase mass parti-
tioning are critical to describe species
transport. A detailed physically-based
model for hysteresis in three-phase k-
S-P relations was described. Simplified
models, which consider effects of
nonwettlng fluid entrapment, were
shown to provide a reasonable com-
promise between accuracy, on the one
hand, and efficiency and robustness,
on the other. Laboratory studies of light
and dense NAPLs in a 1 x1.5 meter
sand tank, involving measurements of
water and NAPL pressures and satura-
tions and component concentrations,
are described. These studies were used
to validate the mathematical model for
multiphase flow and transport.
This Project Summary was developed
by EPA's Robert S. Kerr Environmental
Research Laboratory, Ada, OK, to an-
nounce 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
Spills and subsurface leaks of immis-
cible organic liquids are a frequent cause
of ground-water contamination. This project
was undertaken to improve our under-
standing of contaminant migration arising
from such events and to develop improved
quantitative tools for their description. At-
tention was focused on three fronts in-
volving: (1) the development of efficient
and robust numerical methods for simulat-
ing simultaneous multiphase flow and
transport, (2) developing and numerically
implementing theoretical and empirical
constitutive models governing multiphase
flow and transport processes, and (3) per-
forming laboratory-scale experimental stud-
ies to validate the mathematical models
developed in conjunction with the first two
objectives.
Various formulations of the governing
equations for multiphase flow are derived
and implemented numerically and com-
pared. The most efficient, accurate and
robust formulation employed fluid pres-
sures as the primary variables with satu-
ration time derivatives treated as implicit
Printed on Recycled Paper
-------
functions of pressures, rather than by
chaining to obtain pressure time deriva-
tives. A new algorithm is presented and
discussed which enables substantial re-
ductions in computational effort by auto-
matic elimination and inclusion of elements
in the global solution, based on the de-
gree of transience in the local phase equa-
tions.
The theoretical foundation for modeling
coupled multiphase flow and multispecies
transport with equilibrium interphase mass
transfer is described which leads to a sys-
tem of flow equations for each bulk phase
and phase-summed species transport
equations for each species. A numerical
formulation for solving the system of equa-
tions is presented based on partial
decoupling of flow and transport equa-
tions. Results of several hypothetical nu-
merical simulations are presented to
verify the model and to demonstrate its
applicability to specific problems. The as-
sumption of local equilibrium-controlled
mass transfer is subsequently relaxed; and
a formulation is derived which enables the
form of the phase-summed equilibrium
model to be retained by introducing ap-
parent partition coefficients that depend
on the mass transfer rates and first-order
mass transfer coefficients, as well as the
true equilibrium partition coefficients. Nu-
merical simulations are presented to verify
the formulation and to demonstrate ef-
fects of nonequilibrium mass transfer. Re-
sults of laboratory experiments are pre-
sented that indicate interphase mass trans-
fer can be described accurately at the
laboratory scale by a first-order mass trans-
fer relation and that mass transfer coeffi-
cients may be estimated from empirical
relations previously reported in the litera-
ture.
Numerical and experimental studies
were undertaken to develop and test con-
stitutive models for permeability-saturation-
capillary pressure (/c-S-P) relations gov-
erning three-phase flow. A rigorous, physi-
cally-based hysteretic /c-S-P relation is de-
scribed, as well as various simplified mod-
els. Numerical comparisons of the models
are presented, and results of static and
dynamic laboratory column experiments
are compared with the models to evaluate
their accuracy. Two-dimensional labora-
tory experiments involving simulated spills
of light and dense organic liquids were
also conducted. The design of the experi-
mental tanks and procedures for measur-
ing fluid saturations and pressures and
component concentrations are described.
Experimental results are compared with
numerical simulations to provide valida-
tion of the coupled flow and transport
model.
Methodology and Results
Formulation of Multiphase
Flow Model
The flow of water, air and a nonaqueous
phase liquid (NAPL) in soils is described
by a system of highly nonlinear and
strongly coupled partial differential equa-
tions, which in general must be solved
numerically. Various solution methods may
be derived, based on different forms of
the governing equations and using vari-
ous numerical techniques. Five different
equation formulation methods are investi-
gated for solution of the oil (i.e., NAPL)
and water flow equations with gas at con-
stant pressure via an upstream-weighted
finite element method with Newton-
Raphson iteration for nonlinear terms.
Method 1: Primary variables are oil and
water pressure with saturation time de-
rivatives treated as implicit functions of
pressure. Method 2: Primary variables are
oil and water pressure with saturation time
derivatives chained to yield split oil and
water pressure time derivatives. Method
3: Primary variables are oil-water and air-
oil capillary pressures. Method 4: Primary
variables are oil pressure and water satu-
ration which is a "semi-diffusive" form of
the governing equations obtained by ex-
panding pressure derivative terms into
saturation gradients and pressure-satura-
tion derivative terms. Method 5: Semi-dif-
fusive formulation in terms of water pres-
sure and total liquid saturation.
The semi-diffusive forms (Methods 4 and
5) were found to be very efficient and
accurate under conditions in which all three
fluid phases are present (i.e., air, oil, wa-
ter). However, instability was encountered
in circumstances involving oil infiltration in
the presence of a water table, especially
when oil reached the water table. Method
3 exhibited moderate problems with mass
balance and nonconvergence for some
problems. A drawback to the capillary pres-
sure formulation is that stipulation of
boundary conditions is often difficult since
individual phase pressures cannot be con-
trolled independently on boundaries. The
pressure-pressure formulations are most
convenient in this regard. Method 2 exhib-
ited moderate convergence difficulties and
moderate to severe mass balance errors.
Method 1 exhibited consistently superior
mass balance and convergence behavior
and good efficiency relative to the other
methods.
Adaptive Solution Domain
Method
In most practical multiphase flow prob-
lems, large changes in fluid pressures and
saturations do not occur throughout the
spatial domain at a given time step. Com-
putational effort is thereby inefficiently
spent solving equations in areas where
little activity occurs, rather than concen-
trating effort in the more active zones.
However, since the location of active zones
changes with time (e.g., during oil infiltra-
tion or as phases reach a "field capacity"
following redistribution), schemes designed
to take advantage of these variations must
be capable of automatically adapting to
the current conditions. One method for
doing so involves adaptive grid refinement
in which the numerical mesh is automati-
cally refined in active zones and expanded
in inactive zones. Another general ap-
proach, which empbys a fixed grid, varies
the solution approach within the domain
to gain efficiency. A new Adaptive Solu-
tion Domain (ASD) finite element method
is described in which elements are cat-
egorized as either "active" or "inactive" at
a given iteration. If the element is classi-
fied as active, then it is included in the
global matrix assembly, otherwise it is ex-
cluded. The criteria for changing a node
from inactive to active is based on the
changes in water and oil pressure and
saturation at connected nodes from the
last converged time step to the current
iteration. The criteria for switching an ac-
tive element to an inactive status is based
on differences in the variables at the con-
nected nodes between the converged val-
ues at the previous and the current time
step.
Example problems involving multiphase
flow are presented to illustrate the ASD
method and to compare its performance
with a conventional full domain finite ele-
ment approach. The first example involves
infiltration and redistribution of oil in a
one-dimensional soil column, and the sec-
ond example involves a two-dimensional
domain. Substantial reductions in compu-
tational effort are obtained, which can
range from small improvements to orders
of magnitude, depending on the initial and
boundary conditions for the particular prob-
lem and the duration of the simulation.
For analyses of short-term infiltration be-
havior from small area sources, the im-
provements are generally greatest. Major
reductions can also be obtained for long-
term simulations for cases with steady or
near steady boundary conditions for the
oil phase.
-------
Constitutive Relations for
Three-Phase Flow
In order to model three-phase flow, func-
tional relationships between permeabilities,
saturations and capillary pressures (/c-S-
P) must be known. A model for three-
phase saturation-capillary pressure rela-
tions is derived assuming wettability in the
order: water, NAPL, air. To model hyster-
esis in three-phase saturation-capillary
pressure relations, the concept of appar-
ent saturation may be used. Apparent wa-
ter saturation is defined as the sum of
water and of air and oil trapped in the
water phase. Apparent liquid saturation is
defined as the sum of water and oil satu-
rations and of trapped air. Apparent water
saturation is assumed to be a function of
oil-water capillary pressure, and apparent
liquid saturation is assumed to be a func-
tion of air-oil capillary pressure. Three-
phase apparent saturation functions are
related to two-phase air-water capillary
pressure functions via a scaling proce-
dure that assumes interface curvatures
are related to the interfacial tension of the
fluid pairs. Hysteresis in apparent satura-
tion functions is described using an em-
pirical method based on a variant of van
Genuchten's parametric model. Trapped
fluid saturations are computed as a func-
tion of saturation history using an empiri-
cal relation reported in the petroleum en-
gineering literature. Relative permeability
relations are derived from the saturation-
pressure relations using the semi-theoreti-
cal pore structure model of Mualem.
Due to the computational complexity of
the complete hysteretic three-phase rela-
tions, various simplifications were investi-
gated. A simplified model was derived
which considers hysteresis in saturation-
capillary pressure relations and in relative
permeability-saturation relations due to
nonwetting fluid entrapment, but which dis-
regards hysteresis in apparent saturation-
capillary pressure relations. Numerical
simulations involving infiltration and redis-
tribution of hydrocarbons in unsaturated
soils with fluctuating water tables were
performed which indicated that the simpli-
fied model closely approximates the more
complex model.
Measurements of fluid saturations and
pressures as a function of saturation path
history were taken in two different porous
media to evaluate the hysteretic phenom-
ena. An experimental apparatus consist-
ing of alternating plexiglass sleeves con-
taining treated and untreated porous ce-
ramic rings was designed and constructed
to perform two- and three-phase measure-
ments. Model parameters were obtained
by best-fitting the model to experimental
water and total liquid saturation path his-
tories or by fitting to two-phase air-water
data only and estimating scaling coeffi-
cients from interfacial tension data. Con-
sistent model parameters were obtained
using both calibration methods, and close
agreement was observed between pre-
dicted and experimental water and total
liquid saturation path histories using the
model.
Column experiments were performed to
verify the hysteretic /c-S-P relation under
dynamic flow conditions involving mea-
surements of fluid pressures and satura-
tions during transient flow events. Experi-
ments were performed in two-phase air-
water systems involving water infiltration,
redistribution and fluctuating water table
conditions, and for three-phase flow in-
volving water and oil infiltration and redis-
tribution with a fluctuating water table. A
dual gamma attenuation device was used
to measure fluid saturations, and hydro-
phobic and hydrophilic tensiometers were
used to measure phase pressures. Using
independent calibrations of the hysteretic
k-S-P relation, numerical simulations of
the experiments were performed and com-
pared to the observed data. The detailed
hysteretic model was found to describe
the transient flow quite accurately. The
simplified model exhibited deviations from
the data near the infiltration boundary but
provided good predictions at the distances
further from the source.
Multiphase Transport with
Equilibrium Interphase Mass
Transfer
A model for coupled multiphase flow
and multicomponent transport is derived.
Flow and transport equations are compu-
tationally decoupled by time-lagging inter-
phase mass transfer terms and composi-
tion-dependent phase densities in the flow
equations. Phase-summed component
transport equations are derived for trans-
port by convection, diffusion and hydrody-
namic dispersion in water, NAPL and air
phases, assuming local equilibrium inter-
phase mass transfer. Mass transfer rates
are computed by back-substitution in the
individual phase transport equations after
solving the phase-summed equation at
each time step. Components are assumed
to exhibit linear equilibrium partitioning
among air, water and NAPL phases and
between water and solid phases. First-
order decay is also considered for each
component. Reactions between compo-
nents are not considered.
A numerical model is developed for
three-phase flow of air, water and NAPL
(with options for variable or constant air
pressure) and for transport of up to five
partitionable components in a vertical two-
dimensional domain. The solution proce-
dure involves: 1) solving flow equations
using current mass transfer rates and
phase densities, 2) solving phase-summed
component transport equations serially for
each component, 3) back-substitution in
phase transport equations to compute in-
terphase mass transfer rates, and 4) up-
dating phase densities based on current
phase composition. An upstream-weighted
finite element method is used to solve the
flow and transport equations. Several ex-
ample problems are performed to verify
the numerical model and to investigate
the behavior of light and dense NAPLs.
Transport with Nonequlllbrlum
Interphase Mass Transport
The phase-summed formulation of the
transport equations is generalized to ac-
count for nonequilibrium phase partition-
ing by introducing the concept of apparent
partition coefficients. Under transient field
conditions, at a given location in time and
space, actual phase concentration ratios
may differ from the true equilibrium ratios
due to nonequilibrium interphase mass
transfer. Mass transfer rates are assumed
to be described by first-order mass trans-
fer relations driven by the difference be-
tween the actual and equilibrium concen-
trations in a given phase relative to an
adjacent phase. Apparent partition coeffi-
cients, defined as ratios of actual
(nonequilibrium) phase concentrations are
shown to be a function of the equilibrium
partition coefficient as well as the current
concentrations and mass transfer rates.
The phase-summed transport equation,
which invokes apparent partition coeffi-
cients, thus becomes nonlinear due to the
dependence of apparent partition coeffi-
cients on concentrations and mass trans-
fer rates.
A numerical solution for multiphase flow
and multicomponent transport was devel-
oped for the case of nonequilibrium inter-
phase mass transfer. The solution proce-
dure is identical to that followed for the
equilibrium case, except that apparent par-
tition coefficients are employed in the
phase-summed transport equations, which
are solved iteratively with mass transfer
rates updated at each Picard iteration. A
numerical verification problem was run to
verify the nonequilibrium model's consis-
tency with the equilibrium model, with
which it became asymptotically identical
as the mass transfer coefficients were in-
creased.
-------
A series of laboratory column experi-
ments was conducted to investigate the
nonequilibrium interface mass transfer of
soluble constituents from a multicompo-
nent NAPL at residual saturation during
steady state water flow. Experiments were
performed for different aqueous phase
pore velocities, NAPL saturations, and
NAPL compositions. Breakthrough curves
obtained from the experiments were nu-
merically inverted using a nonlinear opti-
mization program to determine the oil-wa-
ter mass transfer rate coefficients. An em-
pirical correlation between mass transfer
coefficient, pore water velocity, and mean
grain size was evaluated. When the d10
grain diameter was used in the correla-
tion, the mean estimated rate coefficient
was very close to the mean of the mea-
sured values, with a standard deviation of
about 30%.
Intermediate Scale Laboratory
Investigations
Experiments were conducted to simu-
late flow and transport of a light
nonaqueous phase liquid (LNAPL) and a
dense nonaqueous phase liquid (DNAPL)
in an unconfined aquifer with two-dimen-
sional planar symmetry. A steel reinforced
tank, 1.01 m tall x 1.5 m long x 0.085 m
thick with amber transparent plastic sides
resistant to the organic chemicals was
used to perform these experiments. The
flume was instrumented with a dual gamma
attenuation device enabling readings from
66 locations through the narrow dimen-
sion of the cell. At 16 locations, water and
oil pressures were measured using hydro-
philic and hydrophobic ceramic cups, re-
spectively, instrumented with pressure
transducers and tied to a data acquisition
system. Twenty-four Teflon tube inserts
served as aqueous phase sampling ports
in the lower sections of the cell. Aqueous
phase samples were analyzed on a gas
chromatograph using a packed column.
The LNAPL experiment involved five
stages: (1) water drainage by lowering the
water table from an initially saturated con-
dition, (2) oil infiltration on a line source at
the soil surface, (3) fluid redistribution with
a water table gradient, (4) water infiltration
along the entire upper surface, and (5)
NAPL entrapment due to raising the water
table. The experiment was simulated nu-
merically. Soil properties were estimated
from the water drainage stage of the ex-
periment, and fluid and component prop-
erties were independently determined.
Good correspondence was evident be-
tween observed and simulated water satu-
rations at the end of Stage 1 and for
water and oil saturations at the end of
Stage 2. Predicted oil drainage from the
unsaturated zone, vertical oil penetration
at the water table and lateral spreading
during Stages 3 and 4 were somewhat
greater than observed, possibly reflecting
the simplified treatment of hysteresis. Pre-
dicted aqueous concentrations of benzene
and toluene at the end of Stage 3 also
showed greater vertical penetration than
observed due to the greater predicted oil
phase penetration. At the end of Stage 5,
observed and simulated oil saturations
exhibited small upward displacements, re-
flecting the fact that oil had already dis-
tributed to dose a residual state at which
its permeability was very low so that move-
ment was impeded. Aqueous sampling
ports at higher elevations exhibited high
concentrations which were corroborated
by the model.
An experiment involving DNAPL, com-
posed of a mixture of tetrachloroethylene,
soltrol, toluene and iodoheptane, was per-
formed in the same apparatus. The ex-
periment involved three stages: (1) water
drainage from an initially saturated state,
(2) NAPL infiltration on a strip source at
the soil surface, and (3) NAPL redistribu-
tion with a water table gradient. Water
saturations above the water table were
generally slightly underestimated by the
model. The model predicted an oil phase
plume largely confined to a deformed sphe-
roidal region above the water table, al-
though some oil is predicted to have sunk
through the aquifer and spread over the
tank bottom. Between the water table and
the bottom pool, oil saturations after drain-
age were predicted to be less than 1%.
Measurements could not be made near
the tank bottom, so confirmation of oil
spreading along the tank bottom could not
be confirmed. However, nonzero oil satu-
rations at measurement locations 15 cm
from the tank bottom suggest oil penetra-
tion below the water table. The predicted
and measured oil distributions in the un-
saturated zone agree with one interesting
exception. The observed data show a small
amount of oil apparently spread laterally
downgradient along the top of the capil-
lary fringe at oil saturations in the vicinity
of 1%. The simulation does not predict
such lateral migration along the water
table. Whether this is due to physical ef-
fects disregarded by the model (e.g., non-
zero spreading pressures at oil-water in-
terfaces) or to subtle heterogeneities due
to packing anomalies is uncertain.
Substantial concentrations along the
capillary fringe downgradient of the source
confirm the lateral migration phenomena
noted based on oil saturation measure-
ments. Sampling locations under the wa-
ter table and immediately beneath the
source indicate no dissolved phase con-
taminants, although gamma measure-
ments indicate DNAPL at or below these
depths. This suggests downward migra-
tion of DNAPL through the saturated zone
occurs through isolated fingers rather than
as a uniform front.
Conclusions and
Recommendations
Solution of the nonlinear coupled differ-
ential equations governing multiphase flow
and multicomponent transport is a
computationally arduous task. Numerical
results of this study have indicated that
the accuracy and robustness of multiphase
flow problems is quite sensitive to the
method of formulating the governing equa-
tions. In particular, it has been shown that
chaining of saturation time derivative terms
can lead to mass balance errors and other
difficulties that are greatly reduced by
means of an alternative formulation with
saturation time derivatives treated as im-
plicit functions of pressure. Marked reduc-
tions in computational effort were achieved
by an adaptive solution procedure that
takes advantage of the fact that flow may
be near steady state in large parts of the
physical domain. Implementation of
schemes to selectively eliminate individual
phase equations should prove even more
effective in reducing computational effort
as well as improving solution robustness.
Other facets affecting numerical solution
efficiency, accuracy and robustness, such
as matrix solution methodology, should
be pursued particularly to facilitate practi-
cal extensions to three-dimensional prob-
lems.
Solution of coupled multiphase flow and
transport problems introduces additional
numerical difficulties. The semi-decoupled
phase-summed approach employed in this
study is an efficient formulation, although
the efficiency undoubtedly comes at the
cost of accuracy and robustness to some
degree. The most critical aspect of the
decoupled approach is the back-calcula-
tion of interphase mass transfer rates,
since accumulated small errors can even-
tually destabilize the solution. Such prob-
lems are greatly diminished by suppress-
ing mass transfer calculations during peri-
ods of highly transient oil flow. Since com-
positional changes are generally small over
short time periods, mass balance errors
incurred by this approach are very small.
Future studies to investigate alternative
solution formulations and to develop algo-
rithms for automatic component mass bal-
ance controls on the solution should be
pursued (e.g., "re-injection" of mass bal-
-------
ance errors or iteration of transport and
flow with corrected phase transfer terms).
Accurate representation of k-S-P rela-
tions is crucial in order to predict
multiphase fluid movement in the subsur-
face, and accurate models for mass inter-
phase partitioning are critical to describe
species transport. A detailed model for
hysteresis in k-S-P relations has been de-
veloped and has proven to accurately de-
scribe static laboratory S-P data and tran-
sient flow response. Simplified models
which consider effects of nonwetting fluid
entrapment provide a reasonable compro-
mise between accuracy, on the one hand,
and efficiency and robustness, on the
other. Direct measurements of three-phase
relative permeabilities for non-monotonic
saturation histories is a formidable task,
but one that would be highly desirable to
undertake in the future.
Existing empirical relations for estimat-
ing column scale mass transfer rate con-
stants from basic physical properties of
the system yielded surprising accuracy.
However, under common field conditions,
it may be shown that mass transfer limita-
tions should be of minor importance. Field
scala heterogeneity, however, may lead
to preferred flow paths that produce ap-
parent nonequilibrium effects controlled by
diffusive mass transfer between 'last" and
"slow" zones.
Since explicit treatment of heterogene-
ity is generally impractical, future studies
should address the feasibility of describ-
ing field-scale behavior using effective
large-scale mass transfer relations. Like-
wise, effective k-S-P relations at the field
scale that implicitly consider effects of het-
erogeneity may well differ from laboratory-
scale relations. Substantial future efforts
will be needed to more fully understand
and model the behavior of large scale
heterogeneous systems.
GOVERNMENT PRINTING OFFICE: 1992 - 648-0*0/40222
-------
-------
-------
J.C. Parker, A.K. Katyal, J.J. Kaluarachchi, R.J. Lenhard, T.J. Johnson,
K. Jayaraman, K. Unlu andJ.L. Zhu, are with Virginia Polytechnic Institute and
State University, Blacksburg, Virginia 24061-0404.
J.S. Cho and LG. Swaby are the EPA Project Officers (see below).
The complete report, entitled "Modeling Multiphase Organic Chemical Transport in
Soils and Ground Water," (Order No. PB91- 231-514/AS; Cost: $31.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 Officers can be contacted at:
J.S. Cho
Robert S. Kerr Environmental Research Laboratory
U.S. Environmental Protection Agency
Ada, OK 74820
and
LG. Swaby
Office of Research and Development
U.S. Environmental Protection Agency (RD 675)
Washington, DC 20460
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/042
------- |