vvEPA
United States
Environmental Protection
Agency
Robert S. Kerr Environmental
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/S-92/002 March 1992
ENVIRONMENTAL
RESEARCH BRIEF
Multiphase Chemical Transport in Porous Media
J.R Guarnaccia1, P.T. Irnhoff1, B.C. Missildine2, M. Oostrom2,
M.A. Celia1 •*, J.H. Dane2-*, P.R. Jaff61-*, and G.F. Finder3-*
Introduction
The contamination of groundwater by both lighter and heavier
than water nonaqueous phase liquids (NAPLs) is a well-
documented problem. Trichloroethylene is an example of a
commonly encountered heavier than water or dense NAPL
(DNAPL). Because DNAPLs are more dense than water, they
sink through the aqueous phase until reaching an impermeable
barrier, leaving a residual phase behind the advancing front.
The trapped residual in the saturated zone remains essentially
immobile, slowly dissolving into the flowing groundwater. This
residual DNAPL acts as a long-term source of groundwater
contamination.
This investigation focuses on the simulation of the emplacement
and subsequent dissolution of a DNAPL in a saturated
groundwatersystem. The final output of this project consists of
a model that includes two distinct flow and transport processes.
First is the emplacement of DNAPL in the subsurface. Here
two-phase flow dominates as the DNAPL displaces the water.
In addition, as the trailing edge of the slug of DNAPL migrates
into the system, it leaves behind an immobile residual held in
place by capillary forces. Contamination via this scenario
occurs over a relatively small portion of the area of interest and
encompasses relatively short time scales (on the order of
Department of Civil Engineering and Operations Research,
Princeton University, Princeton, NJ 08544
Department of Agronomy and Soils, Auburn University,
Auburn, AL 36849
Office of the Dean of Engineering and Mathematics,
University of Vermont, Burlington, VT 05405
Principal Investigators
days). The second flow and transport process involves the
dissolution of the immobile residual DNAPL through interfacial
mass exchange and transport of the dissolved contaminant in
the water phase. Contamination via this scenario can occur
over a relatively large areal extent; and because dissolution is
a slow process, the source for the contaminant can persist for
long periods of time (on the order of years).
Numerical simulations of this problem are computationally
intensive involving the iterative simultaneous solution of coupled
nonlinear partial differential equations over fine time and space
discretizations. In addition, validity of the model is dependent
on knowledge of a series of physically based relationships that
have to be determined experimentally. This communication
presents an overview of the tasks undertaken to develop an
efficient numerical simulatorto study the DNAPL emplacement/
dissolution problem. These tasks include: 1) the development
of the governing eq uations which describe the flow-dissolution-
transport process; 2) the experimental determination of the
empirical constitutive relationships between capillary pressure
and relative saturation for one representative DNAPL,
trichloroethylene (TCE); 3) the experimental determination of
an empirical modelforTCE dissolution; 4) an examination of the
numerical schemes necessary to model two-phase flow and
associated dissolved transport, and the subsequent
development of a numerical simulator from this study; and
finally 5) the verification of the numerical simulator using
analytical and experimental results.
Governing Equations/Constitutive Relations
The mass balance equations describing the simultaneous flow
of two semi-miscible phases, one wetting (w, e.g., water) and
{§£> Printed on Recycled Paper
-------
one nonwetting (n, e.g., DNAPL), under isothermal,
incompressible conditions are written as a set of four coupled
equations (following Finder and Abriola, 1986):
• two-phase (a) mass balance equations:
at
= w, n
two dissolved constituent (i, a) mass balance
equations:
V.[eSaDaVpia]=piaqa+pi(X
where the constituents (i, a) are a DNAPL species (o) dissolved
fn the wetting phase (o, w) and a water species (w) dissolved in
the nonwetting phase (w, n). The terms in equations (1) are:
Pp
mass concentration of species i in the a phase,
[M/L3]
p° « a phase density, [M/L3], in general a function of
composition
Stt « a phase saturation, defined such that Sw+Sn= 1
£ « the porosity of the porous medium
v° « the mass average a phase velocity [L/T] vector
defined as:
(vp° -pag), a = w.n
(2)
withk - the intrinsic permeability tensor [L2],
the relative permeability which is a function of Sw
IO
H°
p«
g
P.
kinematic viscosity of the a phase [FT/L2], in
general a function of fluid composition
pressure in the a phase [F/L2]
gravitational acceleration vector [L/T2]
the a phase hydrodynamic dispersion tensor
[L2/T]
external source [1/T], positive for influx and
negative for efflux
interface mass exchange term [M/L3T]
defined as:
(It? - p?),. cc = w,n
(3)
C*™ = the mass transfer rate coefficient [1/T] between
the wetting and nonwetting phases, in general a
function of fluid saturation, composition and
velocity; physical and chemical properties of the
species (i) and phases (w and n); and the porous
media properties
K" = the solubility limit [M/l,3] of species i in the a
phase, in general a function of fluid composition
p™ = p ™ + p %, [M/L3T], the total mass flux into or out of the
a phase where in this case p = - p
To close the system, two nonlinear material relationships must
be provided. These are the relationships between capillary
pressure, P (P = Pn- Pw), and saturations, SQ, and between
relative permeabilities, kra, and saturations, Sa. This;
development employs the van Gersuchten k - Sff - P . (K-S-P)
model described in Luckner et al. (1989) and written' here as:;
S QW~ "wr
w -;——
- S
n — — 77; —
l.U -
km(sw)=(snp[i-(i-snl)1H
(4a)
(4b)
(4c)
(4d)
2m (4e)
where Sw is the effective wetting phase saturation, Sw is the
residual wetting phase saturation, &rn is the effective nonwetting
phase saturation, Snr is the residual nonwetting phase
saturation, h = (Pc/pwg) is the capillary pressure head [L], rj
[dimensionless] andwa [1/L] are porous medium dependent
parameters to be fit to data, and m = (1 -1/7}). Measurement of
the fitting parameters, tj and a, is described below; appropriate
data are also given below. For the special case where one of
the fluid phases is present at a saturation at or below its resid ual
saturation, there is only one mobile fluid, and the capillary
pressure - saturation becomes meaningless. Treatment of this
-------
case, which has direct relevance to the problem of dissolving
residual NAPLs, is discussed in the Numerical Methods section
below. Treatment of hysteresis in the capillary pressure -
saturation relation is also discussed in that section.
An Improved Method lor the Determination of
Capillary Pressure - Saturation Curves Involving
Trichloroethylene, Water, and Air
An essential component in understanding and simulating
multiphase fluid flow is the accurate determination of the
hydraulic properties of the different fluids involved. It has been
standard procedure to use pressure cells (Lenhard and Parker,
1987; Turrin, 1988; Demond and Roberts, 1991) to determine
capillary pressure (Pc) - saturation (Sw) curves, where P = P" -
Pw = 20/r (o is the fnterfacial tension and r is the raSius of
curvature atthe interface of two fluids). Because the interiacial
tension between trichloroethylene (TCE) and air at 20°C is only
30 mN/m, and between TCE and water 38 mN/m, as compared
to 72 mN/m between water and air, the height of the pressure
cell may be critical during the displacement of one fluid by
another if Sw changes rapidly with small changes in P , as is the
case for coarser materials. The main objective of this part of the
study was therefore to explore an improved way to determine
Pc- S^ curves for TCE/air and TCE/water for a sand during
imbibition and drainage cycles (hysteresis loops). Two additional
factors of importance, for simulation and cleanup purposes, are
the PC values at which the nonwetting fluid starts to displace the
wetting fluid and the values for the residual saturation of TCE in
TCE/air and TCE/water systems. These values were determined
as well.
Materials and Methods
Aim long column experiment was designed to allow accurate
P - Sw information to be collected at any given point of interest.
The glass column (i.d. = 7.5 cm), with Teflon end caps, was filled
uniformly with Flintshot 2.8 Ottawa sand. The outlet at the
bottom of the column was connected to a supply (or drainage)
bottle by Teflon tubing. This bottle, partially filled with TCE, was
also used to adjust thefluid pressures in the column by lowering
or raising it.
Forthe TCE/air combination, the sand column was subjected to
the following cycles:
Saturation from the bottom, by slowly raising the bottle,
until TCE was ponded on the surface.
TCE displaced by air by stepwise lowering the bottle.
Air displaced by TCE by stepwise raising the bottle.
Upon reaching equilibrium after each step change (no
more flow), dual-energy gamma radiation measurements
were taken at the desired locations to determine the
volumetric content of TCE; Ow (Oostrom and Dane, 1990).
PC values were obtained from knowledge of the height of the TCE
level in the supply/drainage bottle. Corresponding S values
were calculated from Sw= 9w/e. A simplified schematic of the
experimental setup is presented in Figure 1, with the TCE level
CM
91.5 -
81.5 -
71.5 -
61.5 -
51.5 -
41.5 -
31.5 -
21.5 -
11.5 -
S
A
N
D
Q£ n
?H • U
85.7
— .
T
C
E
u
V
. E
R
F
L
0
W
Figure 1. Simplified schematic of experimental setup. The surface of the sand column was 94.0 cm above the reference level (bottom of
column). The TCE level in the supply/drainage bottle (indicated as overflow) was for this particular case at 85.7 cm.
-------
in the overflow tuba for this particular case at 85.7 cm above
the reference level (bottom of sand column). By matching the
corresponding Pc and
obtained.
S values,
w
- Sw data points were
Upon completion of the TCE/air experiments, a 2 cm layer of
water was maintained on the soil surface to displace the TCJE by
water, or vice versa, and P - Swcurves were determined for
TCEAvater in a similar manner as described for TCE/air. The
dual-energy gamma radiation system now determined both the
volumetric TCE content, 9 , and the volumetric water content,
V
The P - S data were fitted with the van Genuchten curve fitting
procedure (equations 4). The Pc entry value for the nonwetting
fluid was taken as 1/a, where a is a curve fitting parameter.
Pesults and Discussion
TCE saturation data (S ) during the displacement of TCE by air
(TCE draining), are shown in Figure 2. The solid line is the
curve as fitted by the van Genuchten procedure (a = 0.148, r\ =
7.244, S -0.059, and saturated Sw= 0.817). The data points
obtained"during the displacement of air by TCE (TCE imbibing)
and the fitted van Genuchten curve are shown in Figure 3 (a =
0.220, ;j - 6.657, S - 0.059, and saturated Sw - 0.742)>. To
appreciate the amowit of hysteresis, the fitted van Genuchten
functions are shown, without the data points, in Figure 4. It can
easily be seen in Figures 2, 3, and 4 that Sw changes from its
maximum to its minimum value, and vice versa, over a capillary
pressure difference of about 5.5 cm of equivalent water pressure.
The data also showthat ignoring hysteresis can have a profound
effect on S . Similar graphs were obtained during the dis-
placement of TCE by water (Figure 5: water imbibing, a -
0.216, TJ = 13.752, S = 0.0, and saturated Sw = 0.763) and for
the displacement of water by TCE (Figure 6: water draining, a
= 0.078, TJ = 10.457, S = 0.093, and saturated Sw = 0.770).
The fitted van Genuchten functions for the water/TCE system
are separately displayed in Figure 7. The change in Sw with
capillary pressure is again very rapid, especially for the (water)
imbibition curve. The amount of hysteresis in this case is even
more pronounced.than for the TCE/air system. ,
The 1/a value for the TCE/air system indicates that air will not
displace TCE unless its pressure is 6.8 cm higherthan the TCE
pressure. For the TCE/water system, the pressure in the TCE
must be at least 12.8 cm higher than in the water, before it will
displace it. ,
Residual TCE saturations, measured at the end of the TCE/air
and TCE/water experiments at a number of locations along the
column, were on average 0.085 and 0.050, respectively. As
was to be expected, the residual TCE-value is greater in the
TCE/air system than in the TCE/water system.
100 n
8-
O.Q
0,2 ;o.4
degree of
0,6 0,8
TCE saturation
1.0
Rgure 2. Capillary pressure (P) - saturation data during the displacement of TCE by air. The solid line represents the fitted van Genuchtan
curve (a = 0.148, TJ ='7.244, residual SOT = 0.059, saturated Sw = 0.817).
-------
100 q
8-
6-
2-
S
o
PH
lo -
0.0 0.2 0.4 0.6 0.8
degree of TCE saturation
1.0
Figure 3. Capillary pressure (P,) - saturation data during the displacement of air by TCE. The solid line represents the fitted van Genuchten
curve (a = 0.220, n =°6.657, residual S = 0.059, saturated S = 0.742).
v ' ' ' ' '
100 -n
0.2 0.4 0.6 0.8
degree of TCE saturation
1.0
Figure 4. Capillary pressure (Pc) - saturation curves during drainage and imbibition of TCE for a TCE/air system.
-------
10 -i
f-c
*cd
a
o
o
PH
6-
3-
2-
0.
>. 6>°
\
<
o
o
o
o
c
t
o
}
0
o
o
0 0.2 0.4 0.6 0.8 1.0
degree of water saturation
Rgura 5. Capillary pressure (P ) - saturation data during the displacement of TCE by water. The solid line represents the fitted van Genuchten
curve (a~0.216, tj = 13.752, residual8^= 6.0, saturatedSw = 0.763).
2-
fn
£ 10
a
o
o
Oc
8-
B-
7-
B-
B-
4-
3-
0.0 0.2 0.4 0.6 0.8 1.0
degree of water saturation
Rgura 6. Capillary pressure (P ) - saturation data during the displacement of water by TCE. The solid line represent!:! the fitted van Genuchten
curve (a = 0.078, tj - 10.457, residual Sm = 0.093, saturated Sw = 0.770).
-------
2-
t-l
.£ 10
cd
a
o
•4-
3-
2-
X.
drainage
imbibition
0.0 0.2 0.4 0.6 0.8
degree of water saturation
1.0
Figure 7. Capillary pressure (P ) - saturation curves during imbibition and drainage of water for a TCE/water system.
Once a capillary pressure (Pc) - volumetric TCE content (9 ), or
saturation (SJ, curve has been determined from measurements
taken at a point, as was done in this study (the collimated
gamma radiation beam was 6 mm in diameter), the P - 9
relationship can be calculated for a soil sample of a given
height. As an example, we assumed a pressure cell with a 6 cm
high sample, with cross sectional area A, initially saturated with
TCE. The reference level for the gravitational head and the
outlet of the pressure cell was taken at the midpoint of the
sample. The PC- 8w relationship displayed in Figure 2 was
chosen as the proper relationship for a point by point analysis.
Forgiven (assumed) air pressures applied to the pressure cell,
we then calculated the PC distribution in the soil samplerthe TCE
pressure (Pw) was computed at different vertical locations
assuming a linear variation in Pw across the sample.
Corresponding volumetric TCE values were taken from the
curve in Figure 2. The average 6 over the whole sample was
then calculated by applying the trapezoidal rule. These average
sample values are graphically displayed in Figure 8, together
with some of the point measured data. It should be noted that
the air entry value of the relationship determined by the gamma
radiation system is clearly defined (P = 5.4 cm of water). The
curve determined on a 6 cm high sample, however, starts to
drain at the top when the applied air pressure is equal to 1 cm
of water pressure (Pc =1 at z = 0 and therefore PC = 5.4 at z =
3). If the Pc- 6^ relationship had been determined on a 6 cm high
sample, the air entry value would have been determined to be
1 cm of water pressure. Due to the early drainage of the top of
the sample and the late drainage of the bottom part, the curve
shows a much more gradual change in TCE content with P than
the one determined at a point. As a result substantial differences
exist in 6w for the same PC .
Dissolution of Trichloroethylene at Residual
Saturation in Groundwater
In many cases it appears that the mass transfer between the
residual DNAPL and water is fast enough so that it can be
assumed water leaving a region of residual saturation has
dissolved concentrations of the DNAPL at the solubility limit.
This is the typical assumption used in models for the transport
of the dissolved organic species in the water phase. However,
for nonuniform distributions of residual saturation, very small
residual saturations, or high water velocities, the dissolution
process might be mass-transfer limited. This part of this study
focuses on the dynamics of the dissolution process in these
situations for a particular chemical, trichloroethylene (TCE).
Dissolution Model
Dissolution of TCE at residual saturation in groundwater will be
a function of many parameters including the water velocity, size
and shape of the resid ual blobs or ganglia, temperature, and the
effects of other chemical constituents and organic matter in the
system. Perhaps the most problematic variable is the structure
of the residual TCE ganglia.
TO circumvent the problem of describing the interfacial area
between the residual TCE and the water, we use a lumped-
parameter model. At a particular location vector X in a porous
-------
A point measurements
o 6-cm high sample
.050 .100 .150 .200 .250
volumetric TCE content
.300 .350
Rguro 8. Capillary pressure (P ) - volumetric TCE content (6w) during the displacement of TCE by air as measured at a single point (triangles) and as
calculated for a 6 cm high sample (circles) from the data obtained at a single point. ,
matrix, the mass transfer for TCE is assumed to be described
by:
described below. The governing equation for the transport of
dissolved TCE in the column is a simplified form of equation
pn
- = C
™
- KS1
(5)
where, as before, the superscripts denote phases [nonwetting
(n) and wetting (w)] and subscripts denote species [TCE (o) and
water (w)]. Implicit in this equation is a functional dependence
of C*0 on the interfacial contact area between the water and
TCE appropriate for some representative volume. We choose
to express this functional dependence by developing a relation
between C*" and 9 (TCE volumetric content, 9 - Sn x e), where
8 is now a representative measure of the inferfacial contact
area. Such a relation will likely be different for different pore
structures and thus different porous media. C™ will also vary
with other system parameters including flow velocity and the
presence of other chemical constituents and organic matter. In
the work presented here, we discuss the relationship between
C*™ and both 8 and the water velocity for a particular sand.
To study the usefulness of this lumped-parameter model, we
obtain local measures of E and 3Sn/3t in a vertical sand column
Impregnated with TCEat residual saturation. The determination
of e and 3S /3t is facilitated by gamma attenuation which is
e (l-
(6a)
_
at
where the hydrodynamic dispersion tensor is now a scalar
quantity, Dw, and the term on the far right replaces the interface
mass exchange term, "pj1. An order of magnitude analysis and
numerical studies suggest that for the experimental conditions
encountered in this work all terms in equation (6a) are small
exceptforthe advection and interface mass exchange terms. In
this case (6a) can be reduced to
asn
at
(6b)
Since e9S /dt can be measured quite accurately along the
length of the column by gamma attenuation, equation (6b) is
-------
integrated from the end of the column where p W(L, t) is
measured to any desired location x. The resulting expression
for P0w(x, t) when substituted into equation (5) yields
Cwn(x,t) =
pne
asn
8t
r ' l i
/•L
P 1 /}^1
ii i i ii i 1 r
Ho^'y ' e(l-Sn)vw c 3,
L •'x
adx' -K£
(7)
Thus, Cwn(x,t) is a function of parameters measured
experimentally: e(x), pow(L,t), and 3S /3t (x,t).
Experimental Method
Five different experiments were run at Darcy velocities ranging
from 0.1 to 1.4 m/day while maintaining a constant room
temperature of 21.5 °C ± 1.0 °C. The Darcy velocity remained
essentially constant throughout each experiment.
Concentrations of TCE in the column effluent were measured in
each of the experiments by extracting 200 u.l water samples into
500 p.l of methanol. At each sampling time, three such samples
were collected and then analyzed by injecting 1 uJ of the mixture
into a gas chromatograph equipped with a flame ionization
detector.
S andewere measured along the column by gamma attenuation.
The primary system components are a gamma radiation source,
450 mCi of 241 Am, and a Nal (Tl) scintillation detector, each
mounted securely on a movable platform. The sample cell is
fixed between the source and detector; and with the collimation
of the radiation beam, a 1 mm thick horizontal slice containing
48% of the cross-sectional column area can be scanned. The
remainder of the system is described by Ferrand, et al. (1986).
The system here is exactly the same as theirs except that the
7Cs source contained in their system was removed, allowing
for more accurate measurements.
Of primary importance is to achieve accurate measures of the
porosity and TCE saturation. Average values of both measures
are determined for each 1 mm horizontal slice of the column
which is irradiated. For the particular calibration and
measurement schemes developed for our work, we find the
following theoretical estimates forthe errors in these measures:
Parameter
Bound for
Systematic Error
±0.003
±0.008
Random Error (One
Estimated Standard
Deviation)
±0.002 to ±0.004
The superscript kdenotes the value of the parameter measured
in a particular horizontal slice k which has been irradiated. The
ranges given for the random error in Snk are indicative of the
errors expected over the range of counting times used with the
gamma radiation system in the five experiments. Note that
since e^ is measured once at the beginning of the experiment,
any error in its value will appear as a systematic error in
subsequent analyses.
A very uniform quartz sand with grain diameters ranging between
0.30 mm and 0.42 mm was used. This sand was washed and
incinerated to remove any fines and organic matter. The oven-
dry sand was packed tightly in an aluminum cell 8.25 cm in
diameter and 3 cm deep on top of a perforated stainless steel
plate covered by a Teflon mesh screen with 0.149 mm openings.
Typical porosities obtained were between 0.34 and 0.37. A 0.5
bar ceramic plate was placed on the sand surface and supported
from above so that the sand column was in compression. An CD-
ring provide a seal between the walls of the aluminum cell and
the ceramic, and a glass fiber filter disk was inserted between
the sand and ceramic plate to enhance the hydraulic contact
between the ceramic and the sand.
The column was gently flushed first with CO2 and then with
deaired deionized water from below, leaving the porous media
saturated with water. TCE was then carefully added from below
until a head of 62.5 cm (±1 cm) of water was established across
the ceramic plate. At this point, the TCE was typically at a
saturation of = 80% throughout the column. To displace the
mobile TCE, two pore volumes of water were pumped through
the column from above at the Darcy flow rate chosen for the
experiment. At the completion of this flush, a complete scan of
the column was made to measure the residual TCE saturation
(Sm) for each 1 mm horizontal slice of the porous media. Over
theentirefive experiments, valuesof S forthese 1 mmthick
sections ranged between 10-19%, wittfs^ typically between
11-14%. The dissolution experiment was then started by
resuming the flow of water and initiating both the gamma
radiation and eff I uent concentration measu rements. The gam ma
radiation system, mounted on a movable platform, automatically
scanned the column at 1 mm vertical intervals continuously,
pausing at each location to collect data. Water samples were
removed intermittently from the column effluent for analysis of
TCE concentrations. The experiment was stopped when there
was no detectable TCE in the column from the gamma radiation
measurements.
A comparison between the gamma attenuation measurements
and measurements of TCE concentrations in the effluent water
provides a mass balance check. For each experimentthe mass
of TCE leaving the column in the water phase was greater than
the mass measured with gamma attenuation. This error ranged
between 12-38% forthe five experiments. Evidence from a test
independent of these experiments indicates that during the
initial TCE invasion of the porous media, a significant mass of
TCE could leak around the O-ring sealing the ceramic plate
against the sand at the top of the column. This extra TCE
cannot be measured with gamma attenuation, and is the likely
cause of these mass balance errors. The computed C^are not
affected by these errors, since they are based on local measures
of e3Sn/3t and po". [Values of E3Sn/3t are measured directly with
gamma attenuation. p0w(x, t) are computed by integrating
(pnasn/3t)/vw(1 - Sn) from each x to the end of the column (x = L)
and adding this to the measured pow(L, t). Equation (7) shows
this computation.] Furtherdetailsof the experimental techniques
and a discussion of these errors can be found in Imhoff, et al.
(1992).
-------
Results
The results can best be illustrated by first examining one
experiment in detail. Here, we will look at the experiment where
a Darcy velocity of 0.53 m/day was used. A plot of the changing
TCE saturation along the column is shown in Figure 9. As clean
water enters from the top, TCE is dissolved in this region and
the TCE saturation decreases. Once the TCE saturation
reaches essentially zero at the very top of the column, a
dissolution front of length X* has formed. Although it is
difficult to observe in this figure, this dissolution front
moves downward through the column, remaining essentially
the same length. After about 200 pore volumes, essentially all
of the TCE entrapped in the column is dissolved away.
Figure 10 shows the changing TCE saturation versus time for
the location 3 mm from the bottom of the column (shown in
Figure 9 with a square around the data). Until the dissolution
front reached this location, there was no measurable TCE
transport due to mobilization of the residual ganglia. Also, there
is no increase in the slope of the dataf or small saturations which
might indicate mobilization of dissolving ganglia. A seventh
degree polynomial was fit through the Sn-t data and 3Sn/9t was
found by taking the derivative of this polynomial. The use of
such a high degree polynomial allowed for an accurate fit
to the curvature.exhibited by the clata at very high and low
S . Such curvature is not significant in Figure 9, but is
more pronounced for other column locations. Similar
analyses were performed on the data collected at these
other locations.
Figure 11 shows the computed Cm versus changing 6nfor all
column locations. The different data symbols indicate data
taken from different column locations, where some repetition of
the symbols was necessary to plot all data. The outliers come
largely from data locations near the very top and bottom of the
column where it was more difficult to obtain accurate measures
of 6 . In the region of the column where the dissolution front
forms ( x < X* ), for a given 9n the numerical value of C*"
decreases as one moves down iri the column. Once the
dissolution front is formed (x > X*), the C™ versus S relation
remains essentially invariant with column location. This same
general trend was observed in the other experiments, and it is
discussed in detail in Imhoff, et al. (1992). It is believed that the
most accurate data are from the region (x>X*), and these data
were used in all subsequent analyses.
0.20 n
Pore Volumes
Flushed
Distance From Top of Column (mm)
Rgure 9. Changing profile of TCE saturation, SB, during experiment. Darcy velocity = 0.53 m/day. See text for further explanation.
10
-------
0.160
0.105 -
n
0.050 -
-0.005
*
*
i i I
60 120
Time(hrs)
180
Figure 10. Changing TCE saturation, Sn, versus time for location 3mm from bottom of the column. Darcy velocity - 0.53 m/day. Error bars
represent ± one estimated standard deviation.
I
225
150 •
75-
0.00
MOVING FROM TOP OF
COLUMN DOWNWARD
0.02
LOCATIONS X > X*
0.04
0.06
Figure 11. Changing mass transfer rate coefficient, C1™, with TCE volumetric content, 6 , for all column locations. Darcy velocity = 0.53
m/day.
11
-------
The trend of decreasing Cm with decreasing 9n is primarily due
to the changing contact area between the water and TCE
phases. Thus, as 6 decreases the contact area decreases and
the mass transfer process slows.
The data from all experiments were put in dimensionless form
using the following parameters:
Sh =
Sc =
D0W
U.W
PWD0W
where 9 is the TCE volumetric content, Sh is the Sherwood
number "(the ratio of the mass transfer velocity to the diffusion
velocity), Sc is the Schmidt number (the ratio of the diff usivity of
momentum to the diffusivity of mass), and Re is the Reynolds
number (the ratio of inertial forces to viscous forces) (Cussler,
1984). The new terms in these expressions are dp = mean grain
diameter [L], and D w = molecular diffusivity of TCE in water
[L2/T].
Sh is plotted versus Re for three representative values of 6n in
Figure 12. The vertical bars indicate the range of data collected
for all x > X" at each Re. For a constant value of 9a , as Re
increases Sh increases in an approximately linear fashion up to
Re=0.01. Thus, as the interstitial water velocity increases, C™
increases although there appears to be some limit to this effect.
The increase in C™ with Re for Re & 0.01 is likely due to a
thinning of the diffusion zone between the surface of residual
TCE phase and the bulk water. This effect is illustrated in Figure
13, where as Re increases, the thickness of the stagnant film (or
diffusion zone) decreases, i.e., I becomes smaller in Figure 13.
A second trend is the increase in C1™ with 9n, which was noted
above for Darcy velocity = 0.53 m/daiy. Here, this same trend
is seen for all Re, indicating that the interfacial contact area
increases as 9 increases.
n
The data were fit to a power-law model of the same form as
Miller, et al. (1990) using least squares regression:
Sh
6, 9nB2ReB3Sc1/2
(8)
:3.18
. 0.54
:0.98
The coefficient of multiple determination (r2) = 0.94. The fitted
parameters are best estimates; confidence intervals could not
Sh
0.31
0.2-
0.1
0.0
en s 0.03
= °'02
6 = 0.01
n
0.00
0.01
0.02
Re
Rgure 12. Sherwood number, Sh, versus Reynolds number, Re, for three values of TCE volumetric content 9n. Vertical bars represent the range of data
collected for aH column locations.
12
-------
TCE? ,,;'<•
S j-fjf V ff
FILM
BULK AQUEOUS PHASE
pw
r
z=o z=/
Figure 13. Stagnant film model for mass transfer. The interfacial region is idealized as a hypothetical film which is unstirred. Mass transfer
involves diffusion across this film from pure TCE to the bulk aqueous phase. K * is the solubility limit of TCE in water p w is
the concentration of TCE in the bulk water. ° o
be determined for these parameters since an analysis of the
residuals from this model indicates the model errors are
correlated (Imhoff, et al., 1992).
The Schmidt number was included in this model for comparison
with the work of Miller, et al. (1990); as with this work, only 9
and Re were varied in their study. The best estimates for the
fitted parameters determined by Miller, et al. (1990) were: 3, =
12 ± 2, 8 = 0.60 ± 0.21, and B3 = 0.75 ± 0.08, where the range
specified for each parameter represents a 95% confidence
interval. There is good agreement for the value of 6,, but
appreciable differences for the other parameters. These
differences may be explained by several factors: Miller, et al.
(1990) used a glass bead media and a different technique for
establishing the initial residual TCE (Sn); measured dissolution
for discrete Sn (and not as Sn slowlyndecreased with time as
here); and used toluene as the dissolving NAPL rather than
TCE (different contact angles may result forthe different NAPLs
yielding different interfacial areas for mass transfer for the
sameSJ. In addition, the variation of C™" with location along the
sample column was observed in this study but remains
unexplained. Miller, et al. (1990) were not able to measure
variations of C1™ along their sample columns. These issues are
explored in detail in Imhoff, et al. (1992). However, additional
experimental work using both different porous media and NAPLs
with different contact angles is required to fully address these
issues. Thus, the mass transfer model developed in this study
is valid only for this porous media and TCE.
Importance of Nonequlllbrlum Dissolution
Given the added complexity in incorporating a mass transfer
model in a numerical simulation of NAPL transport, it is
reasonable to investigate the conditions for which this is
necessary. The experimental work described above provides
insight into this question. In these experiments, mass transfer
was sufficiently fast so that even after the formation of a
dissolution front (see Figure 9), aqueous phase TCE
concentrations reached the solubility limit over distances of
less than 3 cm for Darcy velocities between 0.1-1.4 m/day. S
ranged from = 0% to 10-19% over these dissolution,fronts"
Clearly, if conditions similar to these occur in nature, modeling
nonequilibrium dissolution forthe prediction of aqueous phase
TCE concentrations will be unnecessary. However, if instead
the objective is to capture the dynamics of the changing TCE
saturation, then use of the local equilibrium assumption will give
very different results. Local equilibrium between the water and
TCE phases would require the complete removal of TCE from
an element in a numerical simulator before dissolution could
begin in the adjacent downstream element. Thus, use of the
local equilibrium assumption would not predict the formation of
the dissolution front illustrated in Figure 9.
Even if interest is limited to the prediction of aqueous phase
TCE concentrations, modeling nonequilibrium dissolution could
be important for many real-world situations. Two cases for
which the situation would be different from that considered in
these experiments are: (1) TCE ganglia of similar size and
shape but distributed nonuniformly in the porous media, and (2)
TCE ganglia uniformly distributed throughout the media but
varying widely in size and shape. The numerical model described
below can be used to investigate the importance of modeling
nonequilibrium dissolution for case (1). Various distributions of
residual TCE could be created in a hypothetical aquifer modeled
with the two-d imensional simulator. By using the nonequilibrium
dissolution model in this simulator, downstream TCE
concentrations can be predicted. Comparing these with those
generated from the same simulator but which assumes local
13
-------
equilibrium between the TCE and water will provide information
on the importance of modeling nonequilibrium dissolution.
Case (2) is more difficult to analyze. In this study, the use of a
very uniform sand likely created TCE ganglia of similar size and
shape throughout the media. A natural aquifer will often
contain lenses of coarse and fine material. Wilson, et al. (1990),
using mlcromodels, show how a large NAPL ganglion can be
isolated in a coarse lens embedded in such media. Since the
Interfacial contact area per unit NAPL volume decreases with
Increasing ganglion size, mass transfer can be expected to also
decrease for the larger NAPL ganglia in this media. Although
the form of the dissolution model developed from our
experimental study should be correct for this case, the
parameters may be different.
Thus, while this experimental study provides important
information on the dissolution of TCE in the saturated zone,
more work is required to elucidate those situations where it is
Important to model nonequilibrium dissolution forthe prediction
of aqueous phase TCE concentrations.
Numerical Methods and Model Development
Coupled with the experimental investigations, numerical
methods were developed to solve the equations that describe
two-phase flow in porous media with associated dissolved
transport and mass exchange. The approach taken included
fundamental developments and analysis, as well as construction
of a practical two-dime nsional, two-phase simulator that includes
mass exchange and miscible transport in the aqueous phase.
Standard approximation methods, including finite difference
and finite element methods, were investigated and analyzed in
the context of multiphase flow simulation. From this came a set
of recommendations for use of these approximations. Matrix
solution techniques were also investigated, as were methods
for solution of the miscible transport equations. This work is
summarized below, and the numerical algorithm implemented
to solve the full immiscible/miscible problem is described.
Tobegin, acomplete analysis of the model multiphase equation,
Richards equation, was undertaken. After extensive
investigation of solutions using finite difference and finiteelement
methods, the following conclusions were drawn: (1) Mass
balance problems must be adequately controlled in these types
of equations. Thus the mass-conservative approximation
referred to as the Modified Picard method was recommended.
(2) Infinite element approximations, mass lumping is necessary
to eliminate oscillations that occur ahead of steep moisture
fronts. Thus a lumped finite element approximation, using a
Modified Picard linearization, was recommended. Details of
these results may be found in the paper of Celia et al. (1990a).
Based on these results, a formulation for the two-phase case
wasdeveloped. Detailsofthenumericalformulationareprovided
In Celia and Binning (1992). This algorithm was used to
examine the movement of air in unsaturated soils under
Infiltration. From solution of the two-phase equations, it was
observed that even though Richards equation is often a good
approximation for the liquid movement, the motion of the air
phaseis still an importantfactorforoontaminanttransport in the
air phase. That is, while the motion of the air does not affect
significantly the motion of the water, it can have important
implications for transport of volatile contaminants in the
unsaturated zone. Figure 14 illustrates a transient infiltration
front with associated pathlines for particles of air. Notice that
the air above the wetting front is flowing upward, toward the land
surface, while that ahead of the front is moving downward. As
the front progresses downward, air at progressively larger
depths is pushed upward. This bifurcated flow seerns to be
characteristic of the air phase in unsaturated soils.
One of the outstanding problems in multi-dimensional numerical
simulations for multiphase flow and transport systems is the
matrix solution step. Because the set of flow equations usually
leads to symmetric matrices, a variety of iterative solvers can be
applied with confidence. The pap«r of Bouloutas and Celia
(1991) describes an efficient conjugate gradient algorithm for
these equations. However, when solving the miscible transport
portion of the problem, the matrices are inherently unsymmetric;
and matrix solution becomes a more difficult issue. One
possible approach is to employ a merthod that symmetrizes the
resulting matrix. The Eulerian-Lagrangian Localized Adjoint
Method (ELLAM) of Celia et al. (1990b) is one possible choice.
This has the advantage of also providing enhanced accuracy
due to the Lagrangian approach to advection and also possesses
demonstrable mass conservation properties. This method has
beendeveloped extensively for one-dimensional systems (Celia
and Zisman, 1990), but the two-dimensional implementation is
not yet completed.
A second avenue for efficient snatrix solution involves
approximate factorization techniques. These include operator
factorization or alternating-direction (AD) methods and iterative
relaxation matrix factorization (for example block successive
relaxation). A major motivation for development of these
schemes is to reduce the solution of one large problem to a.
series of solutions of smaller indepandent subproblerns which
can be computed in parallel. One example of AD methods is the
AD collocation (ADC) method of Celia et al., (1987) and Celia!
and Pinder (1990). In addition to being highly parallelizable,
because the ADC method has collocation as its basis, it provides
the added benefit of continuous velocity fields from pressure
solutions on rectangular grids. This is convenient when solving
miscible transport equations coupled to the flow (phase)
equations. However, because this method employs operator
factorization, it cannot directly accommodate cross derivative,
tensorial terms inherent in the governing equations. Therefore,
a collocation based matrix factorization algorithm has been
developed for the two-phase flow and transport equations!
described above. Because this istha algorithm employed in the
computercode delivered to the EPA, a moredetailed description
is provided below.
' ' I
Collocation Discretization/Linearization
Equations (1 a) are rewritten in terms of the dependent variables^
Sw and Pw, by noting that the phase saturations sum to unity, Sy,
+ S =1, by substituting the multiphase extension of Darcy's law,
equation (2) into equation (1a), and by relating Sw to the
capillary pressure, PC(SW) = P" - Fw, when both phases are
mobile, e.g., V Pn = IV Pw + (dPc/dSw)V Sw] (see for example
Aziz and Settari, 1979, pages 133-134). Equations (1b) aro
written in terms of the dependent variable mass concentration,
Pi"-
14
-------
100.0
Water Content vs. Time
80.0 — .
60.0 —
40.0 —r ;
20.0 —
0.0
moisture content
+ 0.06
+0.09
+ 0.12
+ 0.15
+ 0.19
+ 0.22
+ 0.25
+ 0.28
+0.31
o.o
20.0
80.0
100.0
40.0 60.0
Time (min)
Figure 14. Infiltration of water into a column of porous material mat is ininany Tinea with air and water. Each vertical slice of the figure
represents the state of the column at a particular time. The gray shading indicates the moisture content. The lines show the path of
particles of air as the water infiltrates into the column.
The dependent variables, Sw, Pw and p.", are approximated in
space using a linear combination of bicubic Hermite polynomial
basis functions. For two-dimensional problems this results in
four degrees of freedom at each node for each variable (the
function, the two first order derivatives and cross derivative).
The other variables which are space dependent are
approximated using bilinear Lagrange basis functions. The
nonlinear terms, p", k^ and (dP/dS ) in equations (1 a) and S ,
vaand Da in equations (1b), are evaluated at the node and
interpolated into the element using bilinear Lagrange basis
functions. Because the mass exchange term, p;°, represents
the main coupling between the flow and transport equations (a
function of all the dependent variables), it is evaluated using
hermite interpolated dependent variables.
Equations (1) are discretized in time using an implicit backward
differencing scheme. The equations are linearized by the
typically robust and relatively simple to implement Picard iteration
technique, wherein the functions of the dependent variables are
dated at the known iteration level and successively updated
(see for example Huyakorn and Finder, 1983).
The discretized equations are reduced to a set of linear algebraic
equations by employing the orthogonal collocation method: a
method of weighted residuals wherein the weighting function is
the Dirac delta function. This is equivalent to driving the
residual to zero at specified points in the domain which are
denoted as collocation points. Thus no formal integrations are '5
required; and generation of the system matrix is computationally
analogoustothefin'rtedifference method. Orthogonal collocation
results when the equations are written at the four gauss
quadrature points in each element.
System Matrix Generation/Solution
Equations (1) define the nonlinear relationship between phase
flow and constituent transport. In general, the flow equations
(1a) define the distribution of the phases given phase
composition; while the transport equations (1b) define the
phase composition given phase distribution and velocity fields.
While rigorous evaluation of this system of fourequations would
require a simultaneous iterative numerical technique, if we note
that the flow and transport equations are only weakly coupled
through the phase density, viscosity and exchange terms, a
less computationally intensive way to solve the system is .to
decouple the equations. This so called decoupled approach
has been used in the modelling of saltwater intrusion problems
[Celia and Pinder, 1990] and has been applied to multiphase
flow and multicomponent transport problems by Reeves and
Abriola (1988). The decoupled approach can be summarized
as follows:
1.
Given the most recent solution to the transport
equations (1b), the phase flow equations (1a) are in
general solved simultaneously for Pw and S (see for
example Aziz and Settari, 1979, pages 133-1*34). The
transition between two-phase flow and one-phase
flow (DNAPL at residual) is incorporated by forciagk
-------
to go to zero faster than dP^dS^^ goes to infinity. In
those portions of the domain where DNAPL is absent
only the water phase balance equation (1 a) is solved
forPw.
2. Given the flow solution which defines the phase
distribution and velocity field, solve the uncoupled,
weekly nonlinear, transport equations (1 b) for pow and
Pw°-
3. Recalculate the flow equations with the new
compos'itton distribution. Update composition using
the newflow solution. Repeat until changes are small
in some appropriate norm.
This approach allows for smaller systems of equations to be
solved at any one time, allows concurrent execution of each
constituenttransport equation, and allows forthe incorporation
of additional constituents with little restructuring of the solution
algorithm.
The system matrices created by using the above discretizations
are tri-diagonal (this results from the collocation numbering
scheme; see for example Lapidus and Pinder, 1982, pages 337
and 338). During each iteration on the nonlinearities, spatial
decoupling is approximated by using a relaxation matrix splitting
scheme wherein the main diagonal blocks are treated implicitly
and the off diagonal blocks are treated explicitly. This results
In an approximate solution to the system of equations over each
nonlinear iteration. The effect of the relaxation scheme is to
reduce the solution to a series of independent subproblems or
sweeps which can be computed in parallel. For a more detailed
explanation of the parallel solution algorithm see Eppstein et al.
(1992) and Guarnaccia (1992).
Numerical Treatment of Dissolution
A final consideration in the solution algorithm is how dissolution
of residual organic is incorporated into the model. Two major
issues must be addressed. First is the emplacement of the
residual which requires treatment of hysteresis in Sw(Pc). The
hysteresis algorithm used in the simulator is based op the
hysterelic K-S-P models proposed by Kool and Parker (1987)
and Luckner et al. (1989). Specifically, the K-S-P model
presented in equations (4) is augmented to include hysteresis
by defining drainage and imbibition curves based on case-
spocific parameter sets. The result is that water imbibition
curves generally terminate at Sw < 1, where by definition k =
0. In these cases, the phase mass balance equation for the
NAPL reduces to a balance between a change in saturation and
the dissolution rate term. Details of the hysteresis algorithm
can be found in Guarnaccia and Pinder (1992).
The second major issue is the incorporation of the mass
exchange term. The model uses equation (3) with C™ defined
by equation (8). To force the model to assume local equilibrium
between the NAPL and water phases, the parameter S, in
equation (8) would be increased to a very high value.
Model Vertflcatton
Theflowportfonof the code was verified by comparing numerical
results with the exact integral solutions of McWhorter and
Sunada (1990) for two-phase flow. The integral solutions
provide comparison for various horizontal, immiscible phase
displacement scenarios in the presence of capillary pressure.
The K-S-P model utilized was the van Genuchten model
(equation 4). Figure 15 shows a close agreement between the
integral and the numeric solutions for a unidirectional
displacement of a wetting phase by a nonwetting phase. The
fluid and soil properties for this problem are: viscosity ratio =
2.0, S = 0.05, a = 0.098 cm"1,7j = 2.5, 8 = 0.1, k = 5.0 X 10'7
cm2; "The initial condition is S (x, 0) = 1.0; and boundary
conditions are: inflow boundary - S (0, t) = 0.6, outflow boundary
-Sw(100,t) = 1.0.
The mass exchange/transport portion of the code was verified
by numerically simulating the dissolution experiment described
above and illustrated in Figure 9. Tha experiment was designed
to study the dissolution of residual TCE in a uniform sand
column by flushing the system with clean water and tracking the
dissolution front as a function of time. The input parameters
used in the numerical model follow.
The initial conditions for residual TGE saturation are provided
in Figure 16 (t = 0 PV). The initial conditions on the dissolved
TCE species were set at the solubility limit. The mass transfer
rate coefficient was obtained from equation (7). A spatially
variable porosity (average 0.35) was used, based on gamma
attenuation measurements. Other physical parameters were
obtained from the literature at the temperature at which the
experiment was performed. The one-dimensional system was
modelled as a strip of elements with no flow boundaries on the
vertical sides and uniform conditions along horizontal
boundaries. The top horizontal boundary haddirichlet conditions
specified on both water pressure (homogeneous) and the
dissolved TCE species (time varying decay in concentration).
The bottom horizontal boundary had Neumann conditions
specified on both water pressure (such that the Darcy velocity
was 0.53 m/day) and on the dissolved TCE species
(homogeneous). The model length was 25 mm and the spatial
discretization used was 1 mm. The simulation time was 100
hours and the time discretization used was 100 seconds.
Figure 16 shows a plot of the saturation profiles of the
experimental and the numerical resultsfortimes in pore volumesi
(PV) of 50, 100 and 180 PV. The match in simulation results
was excellent from an engineering point of view. Note the close
agreement in terms of mass balsmce. The results predict
complete TCE removal in less than 100 hours or about 200 pom
volumes of water flushed. However, it is apparent that C™ Is
generally too lowtocapturethe precisse shape of the experimental
fronts. The discrepancy is probably due to two basic sources
of error: first, equation (7) was fit to data from five separate
experiments and we are attempting to reproduce only one of
these five; and second, it was fit to data in the lower portion of
the column.
Summary and Conclusions
An efficient collocation based numerical algorithm designed to
solve nonlinear multiphase flow and multicomponent transport
problems in porous media has been described. Specifically,
the target problem is the description of the emplacement and
subsequent removal through dissolution of contaminant NAPLs
in near-surface saturated groundwater environments. The
16
-------
0.6
40 60 80
distance from origin
100
Figure 15. Unidirectional displacement of a wetting phase by a nonwetting phase: comparison between the analytical solution of McWhorter
and Sunada (1990) and the numerically simulated solution.
0.20
6 8 10 12 14 1618 20 22 24
distance from top of column (mm)
Figure 16. Comparison between experimental (black symbols) send numerical results (white symbols) for a dissolution experiment where
residual TCE was removed from a column by flushing it with clean water. This figure uses the same data as that presented in
Figure 9.
17
-------
method is relatively simple to implement and exploits parallel
execution to dramatically reduce computer turnaround time of
the burdensome computations.
Capillary pressure - saturation curves were determined for a
Flintshot 2.8 Ottawa (medium) sand. The fluids involved were
TOE and air, and TCE and water. Saturation values were
determined at desired points along aim long column with a
dual-energy gammaradiatron system. Capillary pressure values
were determined at the same points from known fluid pressures
at hydraulic equilibrium. Matching the saturation and capillary
pressure values at a given point then resulted in P0-Sw data
points. To determine the extent of hysteresis in the Pc-S
relationships, the wetting fluids were subjected to drainage and
imbibition cycles. For both fluid systems the saturation changed
rapidly with capillary pressure changes ranging from 2.5 to 10
cm. This small range of Pc-values makes the use of the
standard procedure to determine Pc-Sw relationships, using a
pressure/tension apparatus and a porous medium sample of
often at least 5 cm in height, undesirable. Hysteresis was
shown to be substantial for both fluid systems and should be
accounted for in computer modeling. The pressures at which
tho nonwetting fluid starts to displace the wetting fluid was 6.8
cm of water for the TCE/air system and 12.8 cm for the TCE/
water system. These values are not only important from a
practical point of view, but also for simulation purposes. The
existence of well-defined displacement pressures suggests
that the use of the Brooks and Corey representation of Pe-Sw
curves may be favored overthe van Genuchten representation
as used in this study.
A new technique was developed to accurately measure and
observe the dissolution of a residual TCE phase in a uniform
sand. Aqueous phase concentration measurements from the
column effluent andgamma attenuation measurements suggest
that the dissolving TCE phase was not mobilized as individual
gangliadecreasedinsize. Masstransferwasshowntodecrease
with decreasing residual saturation, but increase with increasing
watervelocity, at least to a limit of Re=0.01. A power-law model
fit the data well; and although it is only valid for this particular
sand, it can be used by modelers and others to (1) investigate
the importance of modeling nonequilibrium dissolution and (2)
aid in the estimation of the cleanup-time associated with
removing TCE via dissolution.
References Cited
Aziz, K. and A. Settari, 1979, Petroleum Reservoir Simulation,
Applied Science Publishers, London.
Bouloutas, E.T. and M.A. Celia, 1991, "Fast Iterative Solversfor
Linear Systems Arising inthe Numerical Solution of Unsaturated
Flow Problems," to appear, Advances in Water Resources.
Celia, M.A., L.R. Ahuja, and G.F. Pinder, 1987, "Orthogonal
Collocation and Alternating-Direction Methods for Unsaturated
Flow," Advances in Water Resources, 10,178-187.
Celia, M.A. and G.F. Pinder, 1990, "Generalized Alternating-
Direction Collocation Methods for Parabolic Equations: 2.
Transport Equations with Application to Seawater Intrusion
Problems," Numerical Methods for PDE's, 6(3). 215-230,
Celia, M.A., E.T. Bouloutas, and R.L.. Zarba, 1990a, "A General
Mass-Conservative Numerical Solution for the Unsaturated
Flow Equation," Water Resources Research, 26, 1483-1496.
Celia, M.A., T.F. Russell, I. Herrera, and R.E. Ewing, 1990b,
"An Eulerian-Lagrangian Localized Adjoint Method for the
Advection-Diffusion Transport Equation," Advances in Water
Resources, 13, 187-206.
Celia, M.A. and P. Binning, 1992, "Two-Phase Unsaturated
Flow: One Dimensional Simulation and Air Velocities," Water
Resources Research, in revision.
Celia, M.A. and S. Zisman, 1990, "An Eulerian-Lagrangian
Localized Adjoint Method for Reactive Transport in
Groundwaters," Invited Paper, Proc. Eight Int. Conf. Comp.
Meth. Water Resources, Gambolati et al. (eds.), Springer, 383-
392.
Cussler, E.L, 1984, Diffusion: Mass Transfer in Fluid Systems,
Cambridge University Press, Cambridge, pages 227-228.
Demond, A.H., and P.V. Roberts, 1991, "Effects of llnterfacial
Forces on Two-Phase Capillary Pressure-Saturation
Relationships," Water Resources Research, 27(3), 423-437.
Eppstein, M.J., J.F. Guarnaccia and D.E. Dougherty, 1992,
"ParallelGroundwaterComputations Using PVM," Proceedings
oftheNinth International Conference on Computational Methods
in Water Resources, Denver, June (in press).
Ferrand, LA., P.C.D. Milly, and G.F. Pinder, 1986, "Dual-
Gamma Attenuation for the Determination of Porous Medium
Saturation with Respect to Three Fluids," Water Resources
Research, 22 (12), 1657-1663.
Guarnaccia, J.F., 1992, Ph.D. Dissertation, Princeton University,
Princeton, NJ, in preparation.
Guarnaccia, J.F. and G.F. Pinder, 1992, "A New Two-Phase
Flow and Transport Model with Inlerphase Mass Exchange,,"
Proceedings of the Ninth International Conference on
Computational Methods in Water Resources, Denver, June (in
press).
Huyakorn, P.S., and G.F. Pinder, 1983, Computational Methods
in Subsurface Flow, Academic Press, New York, pages 156-
158.
Imhoff, P.T., Jaffe, P.R., and Pinder, G.F., 1992, "An
Experimental Study of the Dissolution of Trichloroethylene in
Saturated Porous Media," Princeton University Water Resources
Program Report WR-92-1.
Kool, J.B. and J.C. Parker, 1987, "Development and Evaluation
of Closed-Form Expressions for Hysteretic Soil Hydraulic
Properties," Water Resources Research, 23 (1), 105-114.
Lapidus, L. and G. F. Pinder, 1982, Numerical Solution of
Partial Differential Equations in Science and Engineering, John
Wiley, New York.
18
-------
Lenhard, R.J. and J.C. Parker, 1987, "Measurement and
Prediction of Saturation-Pressure Relationships in Three-Phase
Porous Media Systems", Journal of Contaminant Hydrology, 1,
407-423.
Luckner, L, M.T. van Genuchten and D.R. Nielsen, 1989, "A
Consistent Set of Parametric Models forthe Two-Phase Flow of
Immiscible Fluids in the Subsurface," Water Resources
Research, 25(10), 2187-2193.
McWhorter and Sunada, 1990, "Exact Integral Solutions for
Two Phase Flow," Water Resources Research, 26(3), 399-
413.
Miller, C.T., M.M. Poirier-McNeill. and A.S. Mayer, 1990,
"Dissolution of Trapped Nonaqueous Phase Liquids: Mass
Transfer Characteristics," Water Resources Research, 26 (11),
2783-2796.
Oostrom, M., and J.H. Dane, 1990, "Calibration and Automation
of a Dual-Energy Gamma System for Applications in Soil
Science," Agronomy and Soils Departmental Series No. 145,
Alabama Agricultural Experiment Station, Auburn University.
Pinder, G.F. and LM. Abriola, 1986, "On the Simulation of
Nonaqueous Phase Organic Compounds in the Subsurface,"
Water Resources Research, 22(9), 1095-1195.
Reeves, H.W. and L.M. Abriola, 1988, "A Decoupled Approach
to Simulation of Flow and Transport of Non-Aqueous Organic
Phase Contaminants Through Porous Media," in Celia, M. A. et
al. (eds), Computational Methods in Water Resources, Vol 1
Modelling Surface and Subsurface Flows, New York, Elsevier,
147-152.
Turrin, R.P., 1988, "Pressure-Saturation Relations for Water,
TCE and Air in Sand," MS. Thesis, Princeton University .
Wilson, J.L., S.H. Conrad, W.R. Mason, W. Peplinski, and E.
Hagan, 1990, "Laboratory Investigation of Residual Liquid
Organicsfrom Spills, Leaks, and Disposal of Hazardous Wastes
in Groundwater," United States Environmental Protection
Agency Report EPA/600/6-90/004, Washington, D. C.
Disclaimer
The information in this document has been funded wholly or in
part by the United States Environmental Protection Agency
under Cooperative Agreement CR-814946 with Princeton
University. This Document has been subjected to the Agency's
peer and administrative review and has been approved for
publication as an EPA document.
Quality Assurance Statement
All research projects making conclusions or recommendations
based on environmentally related measurements are required
to participate in the Agency Quality Assurance (QA) program.
This project was conducted under an approved QA project
plan. The procedures specified in this plan were used without
exception. Information on the plan and documentation of the
QA activities are available from the respective principal
investigators (J.H. Dane for the capillary pressure-saturation
experiments, and P.R. Jaffeforthe solubilization experiments).
19
&U.S. GOVERNMENT PRINTING OFFICE: 1992 - SJ8-080/40202
-------
s-g
-i (0
T33'
3. o>
< co
w co
CO
co
CD
8
o
o
m
T)
<
.§
IIs-
-Q
!
o
s
8
------- |