United States
Environmental Protection
Agency
National Risk Management
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/SR-97/102
October 1997
4>EPA Project Summary
NAPL: Simulator Documentation
Joseph Guarnaccia, George Pinder, and Mikhail Fishman
A mathematical and numerical model
is developed to simulate the transport
and fate of NAPLs (Non-Aqueous Phase
Liquids) in near-surface granular soils.
The resulting three-dimensional, three
phase simulator is called NAPL. The
simulator accommodates three mobile
phases: water, NAPL and gas, as well as
water- and gas-phase transport of NAPL
contaminants. The numerical solution
algorithm is based on a Hermite
collocation finite element discretization.
Particular attention has been paid to the
development of a sub-model that
describes three-phase hysteretic
permeability-saturation-pressure (k-S-P)
relationships, and that considers the
potential entrapment of any fluid when it
is displaced. In addition rate-limited
dissolution and volatilization mass
transfer models have been included. The
overall model has been tested for self-
consistency using mass balance and
temporal and spatial convergence
analysis. The hysteretic k-S-P and mass
exchange models have been tested
against experimental results. Several
example data sets are provided,
including a setup of the artificial aquifer
experiments being conducted at the
EPA's Subsurface Protection and
Remediation Division of the National Risk
Management Research Laboratory in
Ada, OK.
This Project Summary was developed
by EPA's National Risk Management
Research Laboratory's Subsurface
Protection and Remediation Division,
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 physical problem which is addressed
is the contamination of a pristine porous
medium as the result of releases of organic
liquids, commonly referred to as Non-
Aqueous Phase Liquids (NAPLs), in near-
surface heterogeneous granular soils. The
organic liquids can be either lighter than
water or heavierthan water. The soil domain
is idealized as consisting of three
interrelated zones: a vadose zone which is
in contact with the atmosphere, a capillary
fringezone, and a water-table aquiferzone.
A particular problem of interest may include
allthreezonesorasubsetthereof. Granular
soils are stable (non-deforming) and
relatively chemically inert (the soil particles
do not interact with the soil fluids).
Therefore, the soil is idealized as containing
a high percentage of quartz particles and
only a minor percentage of clay particles
and organic matter.
A conceptual illustration of surface-
release-generated NAPL migration in the
vadose, capillary fringe and aquifer zones
is provided in Figure 1. There are three
fundamental mechanisms for NAPL
migration. First, the NAPL infiltrates into
the soil and migrates both vertically and
laterally underthe influence of gravitational
and capillary forces. The distribution of the
NAPL liquid is a function of fluid properties
(density, viscosity, interfacial tension,
wetting potential and variable chemical
composition), soil properties (grain size
distribution, mineral content, moisture
content, porosity, hydraulic conductivity and
spatial heterogeneity), and system forcing
history. If the source is periodic in nature,
then during drying periods, not all the NAPL
will drain from the pore space, leaving
-------
03
0)
to
£
n
precipitation
NAPL source area
soil
moisture
profile
t
water content
vadose zone
capillary fringe zone
.
aquifer zone
t
residual NAPL
dissolved NAPL
vaporized NAPL
mobile NAPL
water table
groundwater flow
Figure 1. Definition illustration of NAPL contamination in near-surface soils due to an intermittent release.
behind an immobile residual, held in place
by capillary forces. If the NAPL is more
densethan water, it will migrate through the
capillary fringe and continue its vertical
migration until either the mobility becomes
zero (all the NAPL liquid is at the immobile
residual state) orthe NAPL front encounters
an impenetrable geologic horizon.
The second contaminant transport
mechanism is dissolution and consequent
advection in the down ward-flowing water-
phase, with precipitation providing the water
source in the vadose zone. In the case of
a DNAPL, flowing ground water picks up
dissolved NAPL constituents.
The third transport mechanism is
transport as a vapor NAPL constituent in
the soil gas, where the increased gas-
phase density induces downward
movement. Partitioning between the gas-
and water-phase contaminants further
enhances the migratory potential of the
NAPL constituents.
DESCRIPTION OF MODEL
Mathematical Statement of the
Model NAPL
The model NAPL is designed to solve the
following system of governing equations
along with initial and boundary conditions.
Mass Balance Equations
Following the development of Pinder and
Abriola (1986), the mass balance law for
each fluid constituent, an ordered pair(i,a)
representing a species iin a fluid phase is:
dt
M + eS
P"
(1)
where the five constituents (i,a) relevant to
this simulator are identified as: (w,I/I/), a
water species in the water phase;(n, I/I/), a
NAPL species in the water phase; (n,N), a
NAPL species in the NAPL phase; (n,G), a
NAPL species in the gas phase; and (g,G),
a gas species in the gas phase.
Other symbols occurring in equation 1 are
used to represent the following:
e is the porosity of the porous medium.
S^ is the saturation of the a-phase.
r* is the mass concentration of species in
1 the a-phase [M/L3].
na is the mass average velocity of a-phase,
a vector [L/T\.
Da is the dispersion coefficient for the a-
phase, a symmetric second-ordertensor
[L2/n
Q* is the point source(+) or (-) a-phase
mass [1/7].
k* is the decay coefficient for species in the
a-phase [1/7].
p is the source or sink of mass for a
species in the a-phase [M/L3T] due to
interphase mass exchange (i.e.,
dissolution, volatilization and
adsorption).
The exchange of mass for each constituent
in equation 1 is defined by:
Ps=(
where
(2)
represents dissolution mass transfer
of NAPL species from the NAPL phase
to the water phase;
-------
Ey, representsvolatilization masstransfer
of the NAPL species from the water
phase to the gas phase;
E^ representsvolatilization masstransfer
of the NAPL species from the NAPL
phase to the gas phase;
En, represents adsorption mass transfer
of the NAPL species from the water
phase to the soil.
Asixth mass balance equation is required
to describe the NAPL species mass which
is adsorbed onto the soil. This equation is
written as:
(3)
where
ps is the density of the soil [MIL3] and
K)S is the mass fraction of the adsorbed
NAPL on the solid [dimensionless].
The balance of equation 3 is replaced by
the following linear equilibrium relationship:
where
Kd is a distribution coefficient [L3/M\.
To ensure global mass conservation, the
following definitions and constraints on fluid
volume, density and mass exchange are
employed:
1. The a-phase saturations must sum to
one:
Sw+SN+SG=l (5)
2. The a-phase mass density [M/L3] is the
sum of the species mass
concentrations in the a-phase:
(6)
l=w,n,g
3. The sum of mass fluxes of all species
into the a-phase must equal the total
mass change in the a-phase:
(7)
4. The total mass change over all phases
must be zero:
(8)
5. The sum of the reacting mass must be
equal to the sum of the produced mass:
=0 a =
(9)
A set of fluid phase mass balance equations
can be generated by summing the balance
equation 1 for each species within the
phase, and by incorporating equations 2,6
and 7. The three resulting fluid-phase
mass balance equations are:
Water-phase:
n
= P
NAPL-phase:
(10)
(11)
Gas-phase:
i=w,n,g
a?
(12)
With this development, the physical problem
can be cast into a mathematical
representation consisting of five mass
balance equations. Of the balance
equations written, the following five are
used in the simulator:
1 to 3) The three fluid-phase balance
equations, equations 10,11 and 12. These
equations define the temporal and spatial
distribution and the flow properties of the
water-, NAPL-and gas-phases throughout
the domain.
4 and 5) The two NAPL species balance
equations, equation 1 with (i,a) = (n, I/I/) and
(n,G). These equations define the temporal
and spatial distribution of the NAPL species
as they are transported within and between
their respective phases.
PHASE ADVECTION
Fluid flow is defined by the parameter Vs
in equations 1 and 10 through 12. The
phase velocity is written in terms of the
multi phase extension of Darcy's law.
ry K K__,
oc = W,N,G
(13)
where
p« is the a-phase pressure
ya =pag is the specific weight of the
phase [MIT2],
g is the acceleration due to gravity [LIT1],
k is the intrinsic permeability [L2],
considered a scalar herein, and
km is the relative permeability.
The a-phase relative permeability is a
scaling factor, 0 < km < 1 , which accounts
for the case where the porous medium is
not fully saturated with the a-phase. This
parameter is in general a function of the
aphase saturation. Given the assumption
that the phase wetting-order is constant,
and follows, from mostto least, water-NAPL-
gas, the following functional dependence is
assumed:
(14)
where the relationships listed reduce to
their proper two-phase forms when
appropriate. The form of equation 13 for
each fluid phase of interest is detailed here
as:
kk
kk
eS
-(V(PW + P )-Y"VZ)
(15)
-------
MASS TRANSFER
Four types of mass transfer processes
are important in defining the physics of the
fate of NAPLs in near surface granular
soils:
(a) dissolution masstransferof pure phase
NAPL to the water phase;
(b) evaporation mass transfer of pure
phase NAPL to the gas phase;
(c) evaporation mass transfer of NAPL
species in the water phase to the gas
phase;
(d) adsorption mass transfer of NAPL
species in the water phase to the soil
phase.
Adsorption of NAPL species in the gas
phase directly to the soil phase is neglected.
Liquid-Liquid Mass Transfer
When the organic phase is at an immobile
residual state, saturation is no longer
considered a function of capillary pressure
since capillary pressure becomes
undefined. Consider the NAPL phase
balance equation 11 for the case of an
immobile residual with constant phase
density, constant porosity, and no external
sources or sinks. For these conditions
equation 11 reduces to:
as,
-= -E: -E
(16)
U I
This equation states that change in NAPL
saturation is due to mass transfer processes.
The dissolution model defining the mass
exchange term, Ew, is assumed to be a
first-order kinetic-type reaction of the form:
where
Cw[1/7] is the rate coefficient which
regulates the rate at which equilibrium is
reached, and
p^ [M/L3] is the equilibrium concentration
of the NAPL species in the water phase
(solubility limit^
In the simulator p^ is assumed to be a
measurable constant value.
To determine the parametric form of Cw,
the work of Imhoffetal. (1992) is employed.
They conducted column experiments
designed to study dissolution kinetics of
residual trichloroethylene (TCE) in a uniform
sand by flushing the system with clean
water and tracking the dissolution front as a
function of time. Using a lumped parameter
model, they derived the following power-
law relationship for CW:
where
p>0.5
and
P3«io
are dimensionless fitting parameters.
The parameter (3^ [1/7] is the rate
coefficient, and it is fit to available
experimental- or field-scale data.
Thevolatilization model definingthe mass
exchange term EG is assumed to follow a
similar model as for dissolution, i. e.:
where
CG [1/7] is the rate coefficient which
regulates the rate at which equilibrium is
reached, and
p~MG [M/U] is the constant equilibrium vapor
concentration of the NAPL species in
the gas phase (vapor solubility limit).
The rate coefficient CG is assumed to have
the form:
C? = pr(eSN)p' (20)
where
(32 is the same as forthe dissolution model,
and
(3^ [1/7] is fit to available data.
Consider now the volatilization of a
dissolved NAPL species in the water phase
to the gas phase. Assuming that the water
phase is at residual saturation in the vadose
zone, and thatthere are no external sources
or sinks of mass, then equation 10 can be
written as:
where
the exchange term Ewis defined in
,-,G
equation 15 and r!y governs the
volatilization mass transfer of a dissolved
NAPL species in the water phase to the gas
phase:
4
(22)
where
H is the dimensionless Henry's law
coefficient which is defined at equilibrium
conditions as follows:
(23)
and
C^ [1/7] is the mass transfer rate
coefficient which is assumed to be
defined by the power law:
CGw=pGW(eSw)p2 (24)
where
the fitting parameter (32 is assumed to be
the same as for the liquid-liquid mass
transfer models, and pGW [1/7] is fit to the
available data.
Liquid-Solid Mass Transfer
Finally, mass exchange due to
adsorption, £/ is assumed to be defined
by a linear equilibrium model:
GOJj=Kdpjf (25)
where
Kd is the distribution coefficient [L3/M]
defined as a function of the organic
carbon content of the soil and the relative
hydrophobicity of the dissolved NAPL
species:
(26)
Kd=focKoc
where
foc is the mass fraction of organic carbon
and
Koc is the organic carbon partition
coefficient. Combining equations 3 and
25 yields the following definition for E^:
'3pr '
(27)
where
pb =[l e]ps is the bulk density of the
soil.
-------
MODEL TESTING AND
EXAMPLE PROBLEMS
A series of pseudo-one-dimensional
example problems were presented in order
to evaluate convergence and mass balance,
and to give the user an indication of
appropriate discretization for a given set of
input data. In addition to addressing self-
consistency, four example problems were
designed to simulate specific physical
experiments:
1. a three-phase LNAPL spill and
redistribution experiment (Van Geel
and Sykes, 1995a, 1995b);
2. a three-phase DNAPL spill and
redistribution experiment conducted at
the EPA's Subsurface Protection and
Remediation Division of the National
Risk Management Research
Laboratory in Ada, Oklahoma;
3. an experimental investigation of the
dissolution of residual DNAPL in a
saturated sand (Imhoff et al, 1992);
4. an experimental investigation of DNAPL
vaportransport in an unsaturated sand
(Lenhard et al., 1995).
The first example simulates water
drainage in a one-dimensional soil column,
1.2 meters long and initially saturated with
water. The boundary conditions are: at the
top, open to the atmosphere, and at the
bottom, specified water head (fine sand =
10 cm, coarse sand = 60 cm). The columns
are then allowed to drain underthe influence
of gravity until quasi-static conditions
prevail. Figure 2 presents the results for
the simulation in the first example.
The second example is summarized as
follows:
1. Fortime = 0 to 100 s, DNAPL is rejected
at a constant volume rate of 0.03 cm3/s
2. for time = 100 s to the end of the
simulation, P°= atmospheric.
Simulation results for two different
discretizations are presented in Figure 3. A
time step of order 2 seconds was required
to obtain a solution for this problem, as
larger time steps caused convergence
problems. Figure 3 shows the total liquid
saturation solution at initial conditions (T =
0) and at T = 200 s(100 safterthe DNAPL
source was removed). One can see that
while both discretizations capture the sharp
DNAPLfront, the x = 2.5 cm solution exhibits
oscillations behind the front. Figure 3
presents saturation results at time = 5000
s, after the DNAPL has migrated to near
static, residual state. It is apparent that for
these model parameters, a grid spacing of
approximately 1 cm is required. Figures
3.2c and 3.2d present mass balance results
for this simulation. In general the model
performs well with respect to mass balance
except when boundary forcing is changed,
and several time steps are needed to
accommodate the discontinuity imposed.
Mass Transfer Model
Here we investigate the convergence
attributes of the kinetic mass transfer model.
Consider the following one-dimensional
water flow and contaminant transport
problem. The initial conditions are set such
that the domain is saturated, and there is a
zone of residual DNAPL, SN= 0. 1 5, uniformly
distributed from x=25 cm to the bottom.
The boundary conditions are set such that
there is a constant influx of clean water at
the top at a rate of 0.008 cm/s, and an
equivalent efflux of contaminated water at
the bottom. Relevant mass transfer and
transport parameters are:
|32=0.5,
and
w i
aL =lcm,
and
pf = 0.001gm/cm3.
The results for different values of the
exchange rate coefficient, fflN , are
presented in Figures 4a and 4b. As shown
in the Figure, a distinct dissolution front is
created, the shape of which is a function of
the size of fflN . High values effectively
approximate the equilibrium partitioning
approximation and produce a sharp front,
while low values produce a broad front.
From a numerics standpoint, the dissolution
front should be resolved over several
elements to minimize oscillations in the
solution which can cause erroneous NAPL
saturations upstream of the source area.
Spatial convergence is illustrated in
Figures 4c and 4d. For a constant rate
coefficient, PjVN=24/d, the model
exhibits convergence as the mesh is refined.
BIBLIOGRAPHY
Imhoff, P.T., P.R. Jaffe, and G.F. Pinder,
An experimental study of the dissolution
of trichloroethylene in saturated porous
media, Princeton University Water
Resources Program Report: WR-92-1,
1992.
Lenhard, R.J., M. Oostrom, C.S. Simmons,
and M.D. White, Investigation of density-
dependent gas advection of
trichloroethylene: Experiment and a
model validation exercise, J.
Contaminant Hydrology, 19, 47-67,
1995.
Pinder, G.F., and L.M. Abriola, On the
simulation of nonaqueous phase organic
compounds in the subsurface, Water
Resour. Res., 22 (9), 109s-N9s, 1986.
Van Geel, P.J., and J.F. Sykes, Laboratory
and model simulations of a LNAPL spill
in a variable-saturated sand medium: 1.
Laboratory experiment and image
analysis techniques, Journal of
Contaminant Hydrology, 17(1), 1-26,
1995a.
Van Geel, P.J. and J.F. Sykes, Laboratory
and model simulations of a LNAPL spill
in a variable-saturated sand medium: 2.
Comparison of laboratory and model
results, Journal of Contaminant
Hydrology, 17(1), 27-53, 1995b.
-------
120
100
80
I 60
13 40
20
0
0
Experimental S-P relationship
0.2
n = 6.49
a = 0.0203/cm
Swr = 0.17
K = 8.64m/d
0.4 0.6
water saturation
0.8
(a)
Computed steady-state moisture profile
(b)
dx = 20 cm
dx = 10 cm
dx = 5 cm
0.4 0.6 0.8
water saturation
120
110
100
90
S 80
6
13 70
60
Experimental S-P relationship
0.2
n= 11.434
a = 0.0849/cm
Swr = 0.08
K = 170 m/d
0.4 0.6
water saturation
(c)
120
110
100
90
80
70
60
50
Computed steady-state moisture profile
(d)
dx = 5 cm
dx = 2 cm
WT
0 0.2 0.4 0.6 0.:
water saturation
Figure 2. Analysis of appropriate grid spacing to compute capillary rise for different soil-types. Parts (a) and (b) are for a relatively fine
sand, and parts (c) and (d) are for a relatively coarse sand.
-------
120
110
100
90
80
70
60
50
WT
T = 0, dx = 2.5cm
T = 0, dx = 1cm
T = 200,dx = 2.5cm
T = 200,dx= 1cm
0 0.2 0.4 0.6 0.8
total liquid saturation
120
100
80
t
60
40
0 0.2
time = 5000
(b)
Sw, dx = 2.5 cm
Sw, dx = 1cm
Sn, dx = 2.5 cm
Sn, dx = 1cm
0.4 0.6
saturation
2
1.5
1
1
0.5
0
1(
(c)
EARLY TIME
dx= 1cm
dt=2s
r
1= 100
)0 200 300 400
time
2r
1.5
cd
X>
0.5
(d)
LATE TIME
dx= 1cm
dt=2s
1000 2000 3000 4000 5000
time
Figure 3. Results of a one-dimensional, three-phase, DNAPL injection and redistribution simulation, highlighting spatial convergence
and mass balance.
-------
t
120
100
80
60
40
20
0
(a)
dx = 2.5cm
ex = 100/d
ex = 24/d
ex = 10/d
0.2 0.4 0.6 0.8 1
NAPL saturation/O.I 5
t
120
100
80
60
40
20
0
(b)
dx = 2.5cm
ex = 100/d
ex = 24/d
ex = 10/d
0.2 0.4 0.6 0.8
Dissolved concentration / 0.001
(c)
ex = 24/d
dx = 2.5cm
dx = 5cm
dx= 12.5cm
(d)
0.2 0.4 0.6 O.S
NAPL saturation/O.I 5
ex = 24/d
dx = 2.5cm
dx = 5cm
dx = 12.5cm
0.2 0.4 0.6 0.8
Dissolved concentration / 0.001
Figure 4. Computational analysis of the dissolution model. Parts (a) and(b) illustrate the effect that the rate constant has on the solution.
As the dissolution front sharpens, oscillation appears indicating that a finer grid spacing is required. Parts (c) and (d) illustrate
spatial convergence for ex=24/d. For the parameters chosen a grid spacing of approximately 5 cm is appropriate.
-------
Joseph Guarnaccia and George Finder are with Research Center for Groundwater
Remediation Design, University of Vermont, Burlington, VT 05401, and Mikhail
Fishman is with Robert S. Kerr Environmental Research Center, Ada, OK 74820.
Thomas Short is the EPA Project Officer (see below).
The complete report, entitled "NAPL: Simulator Documentation," ( Order No. PB95-X;
Cost: $X. 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:
U. S. Environmental Protection Agency
National Risk Management Research Laboratory
Subsurface Protection and Remediation Division
P.O. Box 1198
Ada, OK 74820-1198
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/SR-97/102
------- |