&EPA
United States
Environmental Protection
Agency
Office of Research and
Development
Washington DC 20460
EPA/600/R-01/044
July 2001
Development of a Data
Evaluation/Decision
Support System for
Remediation of Subsurface
Contamination
-------
EPA/600/R-01/044
July 2001
Development of a Data Evaluation/Decision
Support System for Remediation of Subsurface
Contamination
Freddi Eisenberg
Dennis B. Mclaughlin
Parsons Laboratory for Water Resources
and Environmental Engineering
Massachusetts Institute of Technology
Cambridge, MA 02139
Cooperative Agreement Number CR-821894
Project Officer
Mary E. Gonsoulin
National Risk Management Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Ada, Oklahoma 74820
U.S. Environmental Protection Agency
Region 5, Library (PL-12J)
71 West Jackson Boulevard, 12th f to*
Chicago. IL 60604-3590
National Risk Management Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, OH 45268
Printed with vegetable-based ink on
paper that contains a minimum of
50% post-consumer fiber content
processed chlorine free.
-------
Notice
The U. S. Environmental Protection Agency through its Office of Research and
Development partially funded and collaborated in the research described here
under Cooperative Agreement No. CR-821894 to Massachusetts Institute of Tech-
nology. It has been subjected to the Agency's peer and administrative review and
has been approved for publication as an EPA document. Mention of trade names
or commercial products does not constitute endorsement or recommendation for
use.
All research projects making conclusions or recommendations based on
environmentally related measurements and funded by the Environmental Protec-
tion Agency are required to participate in the Agency Quality Assurance Program.
This project was conducted under an approved Quality Assurance Project Plan.
The procedures specified in this plan were used without exception. Information on
the plan and documentation of the quality assurance activities and results are
available from the Principal Investigator.
-------
Foreword
The U.S. Environmental Protection Agency is charged by Congress with
protecting the Nation's land, air, and water resources. Under a mandate of national
environmental laws, the Agency strives to formulate and implement actions leading
to a compatible balance between human activities and the ability of natural systems
to support and nurture life. To meet this mandate, EPA's research program is
providing data and technical support for solving environmental problems today and
building a science knowledge base necessary to manage our ecological resources
wisely, understand how pollutants affect our health, and prevent or reduce environ-
mental risks in the future.
The National Risk Management Research Laboratory (NRMRL) is the Agency's
center for investigation of technological and management approaches for prevent-
ing and reducing risks from pollution that threatens human health and the environ-
ment. The focus of the Laboratory's research program is on methods and their
cost-effectiveness for prevention and control of pollution to air, land, water, and
subsurface resources; protection of water quality in public water systems; remediation
of contaminated sites, sediments and ground water; prevention and control of
indoor air pollution; and restoration of ecosystems. NRMRL collaborates with both
public and private sector partners to foster technologies that reduce the cost of
compliance and to anticipate emerging problems. NRMRL's research provides
solutions to environmental problems by: developing and promoting technologies
that protect and improve the environment; advancing scientific and engineering
information to support regulatory and policy decisions; and providing the technical
support and information transfer to ensure implementation of environmental regula-
tions and strategies at the national, state, and community levels.
This work addresses the feasibility of using dissolved concentration measure-
ments to estimate the spatial distribution of each component of an immobile (or
residual) NAPL mixture in a saturated field-scale system. Also, this report contains
the results of the one-dimensional analysis as well as the development of the two-
dimensional state and estimation equations that are being used in ongoing re-
search. Results from the two-dimensional analysis could be extended to a general
three-dimensional problem and the application of the algorithm to field data. This
report is published and made available by EPA's Office of Research and Develop-
ment to assist the user community.
Stephen G. Schmelling, Acting Director
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
-------
Abstract
Subsurface contamination frequently originates from spatially distributed sources
of multi-component nonaqueous phase liquids (NAPLs). Such chemicals are
typically persistent sources of ground-water contamination that are difficult to
characterize. This work addresses the feasibility of using dissolved concentration
measurements to estimate the spatial distribution of each component of an
immobile (or residual) NAPL mixture in a saturated field-scale system.
We pose the characterization process as a state estimation problem (Mclaughlin,
1995; Bennet, 1992). The estimated states are the dissolved contaminant concen-
trations and the NAPL saturations. These states vary with time and space. The
time-dependence of solute concentrations originating from competitive dissolution
is, in fact, an important source of information for the estimation process. The
estimated parameters are the mass transfer rate coefficients that govern the
transfer of constituents between the NAPL and dissolved phases. In the present
analysis, we assume that all other states (hydraulic head, organic matter distribu-
tion coefficient) and parameters (hydraulic conductivity, dispersion coefficients,
specific storage) are known and certain (not random variables). The state equa-
tions used in the estimation process account for competitive dissolution from a
multi-component source, ground-water flow in a heterogeneous medium, and
dissolved phase contaminant transport. We assume mass exchange of the NAPL
constituents among the dissolved, sorbed, and NAPL phases. Dissolution at the
local scale is assumed to be at chemical equilibrium. However, we represent this
process at the field scale with a first order kinetic model. A linear measurement
equation relates dissolved concentration measurements to the estimated states
and parameters. This equation includes an additive Gaussian error term.
The estimation procedure is based on a Bayesian generalized least-squares
approach. The estimated states and parameters are considered random functions
of space and/or time with prior (or expected) distributions. The estimation algorithm
minimizes a performance index that quantifies the effects of measurement errors,
model errors, and uncertainties in prior information. The state equations and other
constraints are adjoined to the performance index.
In order to investigate the possibility of using this approach at a field site, we
applied the estimation to simulated problems. The first of these is a simple one-
dimensional point source with three NAPL components. One advantage of this
problem is that it is computationally cheap enough to enable repeated runs. Thus it
lends itself to Monte Carlo analysis of varying measurement strategies. The second
problem is a two-dimensional spatially distributed source. This problem requires a
significantly different formulation from the one-dimensional analysis. This report
contains the results of the one-dimensional analysis as well as our development of
the two-dimensional state and estimation equations that are being used in ongoing
research. Results from the two-dimensional analysis could be easily extended to a
general three-dimensional problem and the application of the algorithm to field
data.
This report was submitted in fulfillment of cooperative agreement No. CR-821894
by Massachusetts Institute of Technology under the sponsorship of the United
States Environmental Protection Agency. This report covers the period from
10/01/93 to 06/16/98, and work was completed as of 06/17/2000.
-------
Contents
1. Introduction 1
2. Literature Review 5
2.1 Research Relevant to NAPLs 5
2.2 Research Relevant to Site Characterization 7
3. One-Dimensional Problem Formulation 9
3.1 State Equations 9
3.2 Measurement Model 12
3.3 Estimation Equations 12
3.4 Simulated Problem 14
4. One-Dimensional Results 17
4.1 Measurement Strategies 17
4.2 Parameter and Model Fit 17
4.3 Evaluation of Uncertainty 21
4.4 Evaluation of Measurement Strategies 21
5. Multi-Dimensional Problem Formulation 34
5.1 State Equations 35
5.2 Measurement Model 44
5.3 Estimation Equations 44
6. Conclusions 49
Appendix A: Table of Variables and Operators 51
References 53
VII
-------
Figures
Figure 1: Conceptual representation of the distributed NAPL source problem 3
Figure 2: Mass exchange among phases: aqueous, vapor, sorbed, and NAPL 5
Figure 3: The one-dimensional NAPL point source problem 10
Figure 4: True model outputs for simulated NAPL point source problem 16
Figure 5: Map of measurement strategies for the estimation algorithm 18
Figure 6: Fit of estimated parameters to true values 19
Figure 7: Measurement residuals for strategy 60a versus distance 22
Figure 8: Measurement residuals for strategy 60a versus time 23
Figure 9: Measurement residuals for strategy 60b versus distance 24
Figure 10: Measurement residuals for strategy 60b versus time 25
Figure 11: Concentration profiles for each chemical with superimposed
measurement strategies 26
Figure 12: Prediction error normalized by the posterior covariance (cr*y) versus
distance for strategy 60a 27
Figure 13: Prediction error normalized by the posterior covariance (o+y) versus
time for strategy 60a 28
Figure 14: Prediction error normalized by the posterior covariance (cr^) versus
distance for strategy 60b 29
Figure 15: Prediction error normalized by the posterior covariance (cr^) versus
time for strategy 60b 30
Figure 16: Comparison of algorithm performance for varying measurement strategies 31
Figure 17: Definition of NAPL saturation as a macroscale property 34
Figure 18: Mass balance of each species (i) in each phase (K) in control volume 35
Figure 19: Representation of non-equilibrium NAPL dissolution at the model scale 38
VIII
-------
Tables
Table 1: Parameter Values and Prior Statistics for Simulated Problem 14
Table 2: Chemical Constants Used in the Simulated Problem 15
Table3: Physical Constants Used in the Simulated Problem 15
Table 4: Lag One Correlation Coefficients for Normalized Prediction Errors 31
Table 5: Marginal Cost of Uncertainty Reduction for Pareto Optimal Sampling Strategies. .. 32
Table 6: State Equations for Multi-dimensional NAPL Source Problem 43
Table 7: Set of Euler-Lagrange Equations for Multi-dimensional NAPL Source Problems .... 47
Appendix A: Table of Variables and Operators 51
IX
-------
SI Conversion Factors
Multiply
Area:
Flow rate:
Length:
Mass:
Volume:
Temperature:
Concentration:
Pressure:
Heating value:
English (US)
Units by
1 ft2
1 in2
1 gal/min
1 gal/min
1 MGD
1 ft
1 in
1 Ib
1 Ib
1ft3
1 ft3
1 gal
1gal
°F-32
1 gr/ft3
1 gr/gal
1 Ib/ft3
1 Ib/in2
1 Ib/in2
Btu/lb
Btu/scf
Factor to get
0.0929
6.452
6.31 x10'5
0.0631
43.81
0.3048
2.54
453.59
0.45359
28.316
0.028317
3.785
0.003785
0.55556
2.2884
0.0171
16.03
0.07031
6894.8
2326
37260
Metric (SI)
Units
m2
cm2
m3/s
Us
Us
m
cm
g
kg
L
m3
L
m3
°C
g/m3
9/L
9/L
kg/cm2
Newton/m2
Joules/kg
Joules/scm
-------
Acknowledgments
This project was initiated by Dr. M. E. Gonsoulin and Dr. C. Enfield, U.S.
Environmental Protection Agency. We thank Dr. S. G. Schmelling of the U.S.
Environmental Protection Agency, Dr. R. Michler of the University of North Texas,
and Dr. P. B. Bedient of the Rice University, for their review and valuable com-
ments. Also, we thank Ms. S. J. Elliott of the U.S. Environmental Protection Agency
and Ms. M. Williams of Computer Sciences Corporation for providing graphical
support and developing the final format of the report.
XI
-------
Section 1
Introduction
Nonaqueous phase liquids (NAPLs) are a class of environmental contaminants that exist as a separate phase in the
subsurface due to their low solubility in water. Pools or even pore scale droplets of NAPL in the subsurface behave as
long-term sources of dissolved contamination to ground water. These chemicals are among the most common ground-
water contaminants in North America as a result of decades of production, use, and subsequent release into the
environment (Cohen and Mercer, 1993). Furthermore, NAPL contaminants are frequently considered toxic at dissolved
concentrations far lower than their aqueous solubility (Pankow and Cherry, 1996). Our current understanding of NAPL
behavior in the subsurface relies primarily on complex multi-phase physical models that are most easily applied at the
pore scale and in glass-bead laboratory experiments. The heterogeneities encountered in the field combined with
difficulties in directly observing NAPLs in the field make accurate determination of NAPL location very difficult (Bedient,
1994). Ambiguities in the NAPL source location complicate remediation efforts. To effectively remove, treat, or contain
the source of subsurface contamination, spatial characterization of NAPL distributions is essential.
NAPLs are usually divided into two categories: D(dense)NAPLs and L(light)NAPLs. The density or lightness of the NAPL
is evaluated relative to water. This distinction is made because of the difference in their behavior in the subsurface.
LNAPLs tend to "float" on the water table and thus exist primarily in the unsaturated zone while DNAPLs sink to the
bottom of the aquifer and thus exist primarily in the saturated zone. LNAPLs are usually petroleum products or by-
products. Typical LNAPL chemical constituents are benzene, toluene, ethyl benzene, and xylene. It has been estimated
that over 75,000 underground storage tanks annually release 11 million gallons of gasoline to the subsurface (Hall and
Johnson, 1992), and this figure doesn't include spills or disposal of by-products.
The primary types of DNAPLs found in the environment are halogenated solvents, coal tar and creosote, PCB oils, and
complex DNAPL mixtures. Of these, the most widespread contamination is associated with halogenated solvents
(Cohen and Mercer, 1993), especially chlorinated solvents. These chemicals, most notably trichloroethene, trichloroethane,
and carbon tetrachloride, are used in chemical production, degreasing, dry cleaning, and many other manufacturing
processes. In 1990, twenty-nine billion pounds of chlorinated solvents were produced in the United States (U.S.
International Trade Commission, 1991). It is estimated that these DNAPL constituents contaminate literally thousands of
sites in the United States (NRC, 1990).
Other types of DNAPLs are also important sources of contamination in some areas. Coal tar is a by-product of coal
gassification and steel industry coking operations. In 1990, 158 million gallons of crude coal tar were produced in the
United States (U.S. International Trade Commission, 1991). Creosotes, made up of coal tar distillates, are primarily used
in wood-treating operations. Although current use of creosotes in wood preserving is decreasing, in the first half of this
century creosote was the only available preservative. Because waste management practices at the time involved direct
release into the environment, most wood-treating sites have DNAPL contamination (Cohen and Mercer, 1993). Although
PCB use is now severely restricted, in the past PCBs were commonly used in transformer oil production. Of the 1.25
billion pounds of PCBs produced from 1929 to 1977, 440 million pounds were either placed in land disposal sites or
directly released to the environment (Lavigne, 1990). Of the top 20 contaminants detected in any medium at hazardous
waste sites, 10 are considered DNAPL constituents (U.S. EPA, 1991). These statistics help to illustrate the significance
of DNAPLs as environmental contaminants in the United States.
Compounding the detrimental effect of widespread NAPL contamination are the low toxicity levels associated with most
NAPL constituents. Many NAPL chemicals are considered harmful to humans at part per billion levels. Thus small
amounts of NAPL phase contamination may result in large amounts of dissolved contamination. For example, the
aqueous solubility of pure phase trichloroethylene (TCE) is 1100 mg/L. To produce a dissolved phase contaminant
plume of TCE of 25 thousand cubic meters at 0.01 mg/L requires only 5.3 liters of TCE (Cohen and Mercer, 1993). Note
that this is still twice the Maximum Contaminant Level (MCL) set for TCE by the U.S. EPA. Actual plumes tend to be
much bigger than this. For example, at the Massachusetts Military Reservation on Cape Cod, a plume of TCE and
trichloroethane (TCA) from the sewage treatment area comprises 40 million cubic meters of ground water. The estimated
volume of NAPL that contributed to this contamination was only 1.5 cubic meters, or approximately seven 55-gallon
-------
drums (Mackay and Cherry, 1989). In general, the MCLs for NAPL chemicals are very low, making even very small
amounts of NAPL cause for concern.
This research addresses NAPL characterization by focusing on a subset of the problem: the case of a NAPL present in
the saturated zone at residual saturation. Residual saturation is defined as the pore volume fraction of NAPL that is
held in place by capillary forces (i.e., it is not flowing). The pore volume fraction of NAPL, the NAPL saturation, is defined
as:
Sn=^ (D
p
Sn is the NAPL saturation, Vn is the volume occupied by the NAPL in a control volume of the porous medium, and Vp is
the volume of the pore spaces in the control volume. Thus Sn is a unitless measure that ranges from 0 (no NAPL present)
to 1 (pores completely filled with NAPL). Residual saturationnusually ranges from less than 0.1 to 0.5 (Mercer and Cohen
1990). The value of residual saturation depends on a number of factors descriptive of both the subsurface media and the
chemical properties of the NAPL: 1) pore size distribution, 2) wettability of the NAPL, 3) viscosity of NAPL, 4) interfacial
tension, 5) density of the NAPL, and 6) ground-water flow gradients. A NAPL contaminant (e.g., from a spill or leaky tank)
will percolate through the subsurface leaving behind residual volumes in the pore spaces. The contaminant will flow
preferentially through higher conductivity areas in a heterogeneous medium, resulting in heterogeneous distribution of
NAPL saturation.
NAPL constituents are distributed throughout the subsurface not only via movement of the NAPL itself, but also through
mass exchange among different phases. A NAPL source in the unsaturated zone potentially exchanges mass among
four phases: the aqueous or dissolved phase, the vapor phase, the solid phase, and the NAPL phase. The physical
processes affecting NAPL distribution in the unsaturated zone are further complicated by the movement of the water
table. In the saturated zone, the problem is simplified because the vapor phase may be neglected. NAPL constituents
found in the saturated zone may include free-flowing and residual DNAPLs as well as LNAPLs held at residual. From a
modeling perspective, characterization of a NAPL source in the saturated zone is generally an easier problem to
address.
Characterization of NAPL in the field is difficult and expensive under the best circumstances and nearly impossible under
the worst. Even measuring the presence of NAPL in a sample is not simple, requiring careful sample preservation and
analytic methods. One of the complications is differentiating between contaminant that is sorbed to soil and contaminant
that exists in the NAPL phase. In some cases direct observation of NAPL may be possible in monitoring wells; for
example, coal tar may be easily identifiable as black gooey blobs in a sample of contaminated water or soil. However,
sinking a monitoring well directly over a NAPL blob requires quite a bit of luck, and many NAPL constituents (e.g.,
solvents) are not visually noticeable in any case. Most characterization of NAPL contamination must rely on other types
of measurements; for example, dissolved concentrations or concentrations on soils (Feenstra et al., 1991). Even in this
case, it is not clear what dissolved concentration levels are required to indicate the presence of NAPLs. Often very low
aqueous concentrations of constituents are detected at known DNAPL sites. A general rule of thumb is to suspect NAPL
contamination when constituents are present at one percent of aqueous solubility (U.S. EPA, 1992).
Compounding the difficulties of source characterization, field measurements of dissolved concentrations and head are
sparse in practice—perhaps on the order of 50 to 100 measurements for a given site. Extensive sampling campaigns,
such as the USGS Cape Cod tracer test (100,000 concentration measurements and 1,400 head measurements), are
rare, primarily because they are costly and time consuming (LeBlanc, 1991; Hess, 1992). Taken alone, field measure-
ments provide an incomplete picture of dissolved contamination and an ambiguous characterization of potential NAPL
sources. In this research project we characterize the spatial distribution of residual NAPL saturation by using physical
models of subsurface processes to interpret relevant field measurements, including measurements of dissolved and
sorbed contaminant concentrations, hydraulic conductivity, and piezometric head.
Figure 1 conceptually represents the NAPL source characterization problem studied here. The two darker blobs on the
left-hand side of the picture are NAPL phase sources composed of three contaminants (i.e., red, blue, and green). Both
the total amount and chemical composition of residual NAPL vary spatially. The variability in NAPL composition affects
the contaminant dissolution rate. In general, the contaminant present as the largest fraction of NAPL mass will dissolve
the fastest. The clouds emanating from the NAPL sources in Figure 1 are dissolved phase plumes. At a given
measurement location (the well in the figure), the dissolved concentration of each constituent will change over time as the
source composition changes. Thus the time history of concentration measurements from the well contains information
about the contaminant mass present in the source. These measurements are related to the NAPL source by the physical
processes of fate and transport in the subsurface.
The task of using models to interpret field measurements may be viewed as a state estimation problem (McLaughlin,
1995). The environmental variables used to characterize dissolution, ground-water flow, and contaminant transport in the
-------
/
Tim«
Figure 1. Conceptual representation of the distributed NAPL source problem.
vicinity of a NAPL source are called state variables. The temporal and spatial evolution of these variables is described
by a set of partial differential equations called state equations. These equations are typically based on conservation
laws and associated constitutive relationships that describe the physical system of our problem. Further details for the
specific problem considered in this report are provided in Chapters 3 and 5. The spatially distributed state equations for
this problem have the following general form:
at
y(x, /)=!(«)
(2)
8(y,a)=0
In this equation y(x, f) is a vector of states, ct(x) is a vector of time-invariant parameters, u(x, t) is a vector of known
forcing functions, and v(x, t) is a vector of unknown model errors. The model errors are assumed to be zero in
deterministic modeling applications but may be treated as random functions in stochastic modeling applications. A is a
(usually) non-linear spatial operator that describes the relationships among the variables. Gj is the model forcing
operator. The operators land 6 are the initial and boundary conditions, respectively, for the state equation.
In this investigation, the states y are NAPL residual saturations and solute concentrations. The parameters a include the
ground-water flow velocities, mass transfer phase partitioning coefficients, and initial conditions for all of the states. The
state equations are the set of equations for dissolved phase contaminant transport, competitive dissolution, and sorption
as well as physical constraints on some variables. The initial and boundary condition operators are specified
concentration conditions.
The state equations must generally be solved numerically, after discretizing all states and parameters over a
computational grid that covers the region of interest (specific solution methods are discussed further in Chapters 3 and
5). The resulting spatially discretized solution is a vector y(a, t) of states defined at the A/y nodes of the computational
grid, where a is now interpreted as a spatially discretized parameter vector. For convenience we assume that the
parameters and states are discretized over the same computational grid.
-------
The measurement equation relates field measurements to the states, the parameters, and a measurement error term
at discrete times and locations (n):
zn = Mn(y,a) + a>n ; n = l,2,...Nz (3)
where M'„ is a measurement operator which describes how the measurement zn is related to the states and parameters
and con is an additive random measurement error. The quantity AYn(y, a) is just the model's prediction of the measurement
zn, which is obtained at location xn and time tn. The measurement operator may be as simple as a linear function of the
states (additive error) or something much more complicated. In this investigation the measurement equation relates field
measurements of dissolved concentrations to the true values.
In our approach to inverse estimation the uncertain variables a, v, and con are presumed to be random fields with specified
prior statistics (mean and covariance). The goal of the estimation procedure is to find the values of these variables that
give the "best fif to available measurements, while satisfying the dynamical constraints imposed by the state equations.
The overall quality of a given set of estimates and states is measured by the following discrete generalized least-squares
performance index:
1 r /\ — iT --,_i r ^ — 1
T[«-«] Ca'[a-a]
z (4)
The first term of this performance index is a standard weighted least-squares measure of the difference between a vector
composed of all N^ measurements (z) and the vector M(y,a) composed of the corresponding N discretized model
predictions. The weighting factor is the measurement error covariance matrix (CJ. The second term forces the estimates
to stay close to specified prior discretized parameter values. This term is weighted by the prior parameter covariance
(C ). The third term measures the magnitude of the discretized model error (u), weighted by the model error covariance
(Cj. The "hats" over the states and parameters indicate that they are estimates of the true values. These are called
least-squares estimates because they minimize the generalized performance index of Equation (4) (Mclaughlin and
Townley, 1996). The least-squares estimation approach described here is discussed in more detail in Mclaughlin (1995)
and in the later sections of this report.
It is important to note that the model error and prior terms of the performance index "regularize" the estimation problem
by making it less sensitive to anomalous measurements. Least-squares optimization problems usually have more than
one local solution, making minimization of the performance index numerically difficult. By restricting the domain of the
minimization algorithm, regularization alleviates some of the ill-posedness inherent to non-linear least-squares estima-
tion (Sun, 1997). Regularization is frequently required to obtain stable, physically reasonable estimates from real field
data (Mclaughlin and Townley, 1996).
The remainder of this report describes the development of a method for estimating the spatial distribution of residual
MAPI sources. The success of the method is evaluated by applying it to a synthetic problem. Section 2 provides a review
of previous research on MAPI processes and related characterization issues. Section 3 explains the formulation of a 1-
D MAPI point source estimation problem. This problem was designed to assess the feasibility of characterizing a NAPL
source from solute measurements before attempting a more complicated multi-dimensional problem. It also provides
some useful insights on the design of NAPL monitoring systems. Section 4 presents the results of the 1-D analysis and
evaluates the success of the estimation algorithm. Section 5 formulates for the multi-dimensional distributed source
estimation problem. The two-dimensional formulation turns out to be significantly different from that of the 1-D problem.
Section 6 discusses conclusions from this research and identifies promising topics for further work.
-------
Section 2
Literature Review
Up until the last decade, most investigation of immiscible fluids in the subsurface was conducted from the perspective of
petroleum exploration and extraction. In fact, the term NAPL was coined in 1981 (Pankow and Cherry, 1996). The first
significant research on non-dissolved phase contamination in the subsurface is usually attributed to Fredrick Schwille
(1988) who conducted laboratory experiments to observe how DNAPLs flow into porous media. Much of the current
hydrologic study of NAPLs concerns the small-scale behavior of immiscible fluids in water or models of multiphase flow
(Abriola and Pinder, 1985; Pinder and Abriola, 1986; Faust et at., 1989; El-Kadi, 1992). In this section, we focus on NAPL
research concerning the distribution and fate of residual NAPL in the saturated zone. Special attention is paid to
modeling residual NAPL dissolution (Powers et al., 1994; Mayer and Miller, 1996; Sleep and Skyes, 1989). Little if any of
the current literature addresses the problem of systematically characterizing NAPL sources from measurements of
dissolved concentrations. The second part of this chapter reviews the relevant parameter estimation literature. The
overall goal of this section is to provide a framework for the variational state estimation method used in this research.
2.1 Research Relevant to NAPLs
Physical Processes of Residual NAPL
Residual NAPL mass may undergo three phase transformations in the subsurface: dissolution (NAPL phase to dissolved
phase), sorption (NAPL phase to solid/soil phase), and volatilization (NAPL phase to air phase). Figure 2 is a conceptual
representation of these processes at a pore scale. This research will only address dissolution and sorption, those
processes likely to occur in the saturated zone.
Many researchers have explored sorption mechanisms and rates for organic chemical sorption onto subsurface media
(soils). When the organic content of the soil is high enough, sorption onto the organic matter present on soil particles is
usually the dominant sorption mechanism. (Schwartzenbach et al., 1993). Thus sorption is often assumed directly
related to the organic content of the soil (f ).
Single Species at Equilibrium
c =
L
Multiple Species at Equilibrium
c = y r csat
°a,y Aj'jaJ
Figure 2. Mass exchange among phases: aqueous, vapor, sorbed, and NAPL for chemical j.
-------
Most researchers present sorption as either an equilibrium or a first order kinetic process. If the time scales of sorption
are much smaller than the time scales of transport, an equilibrium sorption model is reasonable. Chemical equilibrium
suggests a simple linear relationship between the dissolved contaminant concentration and the sorbed concentration:
Where ca is the concentration of the chemical in water, cs is the sorbed concentration, and Kd is the partitioning coefficient
(see Appendix A for definitions of mathematical symbols). The partitioning coefficient varies based on the hydrophobicity
of a given chemical. In ground-water transport of contaminants, sorption is usually incorporated into the transport
equation as a retardation factor on contaminant movement.
Dissolution may also be represented as an equilibrium process, in this case between the aqueous and NAPL phases. In
the case of a single component NAPL, the dissolved phase concentration is given by the pure phase saturated solubility
of the chemical. Solubilities may be determined experimentally or estimated based on the structure of the chemical
(Schwartzenback et al., 1993). For most chemicals of interest, solubility values may be found in the literature. DNAPL
constituents generally vary in solubility from 104 to 10'2 mg of chemical /liter of water (Cohen and Mercer, 1993). Solubility
may be affected by temperature, salinity, the presence of cosolvents, and the presence of dissolved organic matter such
as humic and fluvic acids. In a NAPL mixture (i.e., more than one chemical present), the equilibrium dissolved
concentration is called the effective solubility. NAPLs form a mixture dissolved competitively based on their pure phase
aqueous solubilities, their activity coefficient in the NAPL mixture, and their mole fraction in the NAPL mixture. This
relationship is generally based on Raoult's Law:
(6)
where m( is the moles of chemical j in the NAPL mixture, c/*' is the pure phase solubility of chemical j in water, celf is the
effective solubility, and m} is the activity coefficient of j in the mixture. If the chemical components of a NAPL mixture are
structurally similar, the activity coefficient is near unity (Banerjee, 1984). Since the effective solubility of NAPL mixture
components is related to their mole fraction in the mixture, the effective solubility of a chemical will change over time as
mass is dissolved and the mole fraction in the NAPL mixture changes. Generally, the mixture left behind becomes
increasingly less soluble (Mercer and Cohen, 1990).
Modeling NAPL Dissolution
Early models of NAPL exchange with the aqueous phase tended to assume equilibrium dissolution. In lab column
experiments, NAPL is seen at aqueous solubility for flow velocities of 10-100 cm/day (Anderson et al., 1992). However,
concentrations equal to effective solubilities are rarely found in the field, even when NAPLs are known to be present
(Pankow and Cherry, 1996). Mackay and Cherry(1985) suggests that field concentrations are generally found at less
than 10% of effective solubility. This discrepancy may be caused by heterogeneity in the NAPL distribution (leading to
dilution), heterogeneity in the porous medium, mixing of stratified water in sampling wells, and/or mass transfer
limitations that make dissolution a non-equilibrium process (Feenstra and Cherry, 1988; Powers et al., 1991). Brusseau
(1991) concluded that the equilibrium assumption is appropriate for subsurface pore velocities of less than 0.2 cm/hour,
but that a non-equilibrium model is more appropriate for velocities greater than 1 cm/hour. The results of Burr (1994)
indicated that on a field scale, the selection of an equilibrium (versus non-equilibrium) sorption model is not significant in
the face of aquifer heterogeneity. Some research suggests rate limited or non-equilibrium dissolution is an appropriate
model when NAPL is present in large pools, at low saturation values, and at high ground-water flow velocities (Powers
et al., 1992). Research by Imhoff et al. (1993) concludes that the narrow width of dissolution fronts observed in the
laboratory suggests that the equilibrium relationship is achieved over small spatial scales (although this is affected by
heterogeneity of the medium).
Non-equilibrium models of NAPL dissolution feature mass transfer rate coefficients that may be developed in various
ways. A number of studies of residual NAPL dissolution indicate that dissolution rates may be estimated as a function of
interfacial area and the Darcy flux of the ground water (Powers et al., 1992; Miller et al., 1990; Parker et al., 1991). In
other studies, the transfer rate coefficients are functions of the NAPL saturation itself (Miller et al., 1998). Thus as the
amount of NAPL decreases, dissolution slows. This tailing process is often seen in the field and is partially responsible
for longevity of NAPL sources. In all cases, the rate coefficients are fit to empirical (usually column) data. Thus they apply
over small scales of homogeneous media.
-------
Identification of Residual NAPL
Residual DNAPL is difficult to measure and characterize in the subsurface. The presence of NAPL at a field site may be
assessed by inspection if the volume content in the medium is sufficiently high (e.g., streaking is visible or free product
is found in wells). However, residual NAPL in soil is usually well below this level. If a soil sample is taken within a NAPL
zone, the NAPL saturation may be estimated in the laboratory (Feenstra et al., 1991). The total amount of contaminant
is measured (i.e., sorbed and NAPL phase). Then the amount of NAPL phase present is back-calculated by estimating
the maximum amount expected to be sorbed onto soil and comparing this value with the measured value. While this
method may provide an assessment of the presence of residual NAPL, 1) it depends upon accurate estimation of
sorption coefficients and 2) it depends on having a soil sample from a NAPL zone. Since DNAPLs are generally deep in
the subsurface, it is difficult (and uncommon) to have such samples.
It is also difficult to predict NAPL saturation based on porous media characteristics. Wilson (1990) found that residual
saturation cannot be predicted from soil texture because heterogeneity within the media and minor difference in texture
(e.g., lenses) make a large difference in saturation. Further, residual saturation values in heterogeneous media tend to
be larger than in homogeneous media. Powers et al. (1992) found that graded sand held more NAPL in residual than
uniform sand with the same mean grain size.
A limited number of experimental releases of NAPL have provided some insight into the distribution of residual NAPL in
the subsurface. Kueper et al. (1993) studied the spatial distribution of residual NAPL after a release of PCE in the Border
aquifer. They found residual NAPL saturation ranging from one to thirty-eight percent of available pore space. Spatial
variability was very dependent upon soil grain size and the maximum nonwetting saturation along the drainage path (i.e.,
where continuous mobile NAPL was present). Their results also show that the ultimate distribution of residual NAPL is
dependent on the rate of release of NAPL at the surface. Even for a homogeneous medium, S should be expected to
vary significantly. These results are similar to the range of reported values for NAPL retention of one to fifty percent (for
widely varying media) by Mercer and Cohen (1990). Therefore, even a general knowledge of the initial NAPL source is
insufficient to characterize the resulting residual distribution.
2.2 Research Relevant to Site Characterization
The use of inverse/state estimation methods for subsurface characterization has been reviewed by a number of authors
(Yeh, 1986; Kuiper, 1986; Carrera, 1987; Ginn and Cushman, 1990; Sun, 1994; McLaughlin and Townley, 1996). In
particular, Zimmerman et at. (1998) recently compared the effectiveness of seven different inverse methods at solving
the same test problems. It is not within the scope of this section to describe all inverse/state estimation methods. Instead
we provide a context for our current research. Neuman (1973) describes two general categories of inverse methods:
direct methods and indirect methods. The indirect methods systematically fit the output values of a model to observed
values. The goal is to minimize the residual between the observed values and model predictions by adjusting parameters
in the forward model. Direct methods, by contrast, assume that states (the results of the forward model) are known at all
locations on the model grid. Then the solution of parameters may be posed as a linear problem. In ground-water
problems, however, states are only known at a few locations, if at all. Thus most researchers focus on indirect methods.
Indirect methods are inherently non-linear and suffer problems of instability and non-uniqueness. Generally, regulariza-
tion terms are added to the minimization to make the method successful, as discussed in Section 1.
In their review paper on inverse methods, McLaughlin and Townley (1996) characterize inverse methods by four criteria:
1) the way spatial variability is described, 2) the forward equation used to relate measurements to parameters, 3) the
performance criterion chosen, and 4) the solution technique used. They show how a number of inverse algorithms can
be described in a maximum a posteriori framework. These include linear methods (Dagan, 1985; Hoeksema and
Kitanidis, 1984; Sun and Yeh, 1992) and non-linear methods (e.g., maximum a posteriori methods, non-linear least
squares, maximum likelihood methods, pilot point methods, and extended Kalman filtering).
Reid and McLaughlin (1994) used a generalized least-squares approach to estimate a spatially distributed log
transmissivity function from measurements of piezometric head. They approximate the transmissivity function by a
discrete expansion in a series of basis functions called "representers." The representers are related to the cross
covariances between measurements and log transmissivities. Reid (1996) expanded on her previous work by estimating
three-dimensional hydraulic conductivity as well as head, and concentration using measurements of all three variables.
She again parameterized the hydraulic conductivity field with representer expansions. The algorithm was successfully
applied to data from a field site to test the effectiveness of the approach. Half of the measurements from the site were
used in the algorithm, and the other half were withheld to test the accuracy of the estimator's prediction.
The work of Sun (1997) furthered the above research by developing an estimation algorithm that includes sorption of the
chemical contaminant onto soil. His technique estimates hydraulic conductivity, dissolved concentration, piezometric
head, and chemical/soil partitioning coefficient fields. The predictions are conditioned on measurements of head,
hydraulic conductivity, dissolved concentration, and fraction of organic matter in soil. Sun's method follows Reid's in its
-------
basis function expansion within an iterative Gauss-Newton estimator. However, Sun advances the sophistication of the
approach by using an adjoint state approach to calculated the sensitivity coefficients and cross-covariances. Like Reid,
Sun applies his method to a field site by using half of the measurements to condition the algorithm and the other half to
test the accuracy of the estimation.
Valstar (1998) takes a similar approach to Reid and Sun. He developed an algorithm to estimate piezometric head and
dissolved concentration fields from measurements of both states. He uses a variational approach to state estimation,
similar to the approach used here. However, his problem assumes a conservative contaminant source with a known
initial condition and an unsteady flow field. Valstar also includes model error in his performance index to account for
uncertain forcing in the ground-water flow equation. He solves his estimation equations using representers for all states
and parameters and adjoint variables.
As far as we are aware, very few attempts have been made to characterize NAPL contamination using inverse/state
estimation methods. Jin et al. (1995) used a non-linear least squares technique to estimate residual NAPL saturations in
a partitioning tracer column experiment. They estimated the average residual saturation over the entire column by fitting
predictions and observed measurements. Mackay and Wilson (1995) conducted a similar study, again with partitioning
tracers in a column experiment.
James et al. (1997) developed a method to estimate residual NAPL distribution using partitioning tracer data. The
method is based on a stochastic derivation of covariances and cross-covariances among tracer temporal moments,
residual NAPL saturation, pore water velocity, and hydraulic conductivity fields. They perform both an unconditional and
a conditional analysis. The conditional analysis incorporates measurements with a Bayesian estimator technique. This
algorithm was tested on a three-dimensional synthetic data set (a 5x4x2 meter domain) modeled after the partitioning
tracer test conducted at Hill Air Force Base (Annable, 1997). Using temporal moments (zeroth and first moments) for 72
simulated breakthrough curves, they created both a spatial prediction of residual NAPL and an estimate of the
uncertainty of the prediction. Overall, the method reduced uncertainty in the estimated NAPL distribution by 25%. This
method is similar, at least in its overall objective, to the research described in this report. However, it requires data from
a partitioning tracer test, which are not available at most field sites. Nonetheless, the results of James et al. (1997) make
a useful comparison to the results of our approach.
-------
Section 3
One-Dimensional Problem Formulation
The first part of this research is concerned with the feasibility of estimating NAPL source locations and compositions from
dissolved-phase measurements. The feasibility question is important since the processes involved are non-linear and
the sensitivity of the estimation to the measurements is unknown. Estimation of residual saturation from dissolved
concentration measurements can succeed only if Sn is sufficiently sensitive to the measurements. Otherwise the
estimation problem will be ill-posed (i.e., the solution may be non-unique or very sensitive to small measurement errors).
The success of the estimation is also affected by the degree of non-linearity of the forward model. Our estimation method
is based on a series of local linearizations of the forward model in the neighborhood of the current parameter estimates.
If the forward model is very non-linear, even a local linearization may produce a poor parameter estimate.
In order to examine the feasibility of our estimation objective we constructed a simplified one-dimensional (1 -D) problem.
This 1 -D problem is also useful for exploring the relative effectiveness of varying sampling strategies. While a multi-
dimensional or distributed source problem would be too computationally expensive to run for multiple scenarios, the 1-D
problem is very amenable to exploratory analysis.
In our 1-D formulation, the residual NAPL mixture is considered to be a point source of dissolved phase contamination
at an uncertain location. We treat the source as a time-varying concentration boundary condition for a ground-water
transport model. The source concentrations are calculated from the dissolution of a NAPL mass. The unknown
parameters are the initial masses of each constituent, or chemical, in the NAPL mixture and the location of the source.
Measurements of dissolved concentration are used to estimate these parameters. The remainder of this section
describes the mathematical formulation of the state equations, the measurement equation, and the estimation algorithm
as well as the simulated problem, which represents the "true" set of states estimated in Section 4.
3.1 State Equations
The state equations for this 1 -D problem describe the movement and transfer of contaminant mass among the dissolved,
sorbed, and NAPL phases. These equations may be written in terms of the notation introduced in Section 1:
— = Am(c,s,m,x,)
(7)
ds(x,t)
\ / =
at
The three states are the dissolved concentration (c), the NAPL phase chemical in moles (m), and the sorbed phase
chemical (s). These states are related to each other through the operators, Ac, Am, and A which are dependent upon all
of the states as well as the source location, xs. Note that the model error is assumed to be zero in this 1-D example.
The dissolved phase portion of the model is described by the standard 1 -D advection/dispersion equation including linear
sorption (Bear, 1979):
dc{(x,t) dc^t] _d2Ci(x,t) 1 ds.
~T1 ^^7=4(c>s>m>*.
(8)
subject to the conditions:
-------
c,(0,f) =
on the domain bounded by x=0 and x=L (note that the concentration boundary conditions are imposed at the uncertain
source location, xs, as well as at x=0 and x=L). The ground-water flow velocity (v), the dispersion coefficient (D), and the
effective porosity (ne) appearing in this state equation are constant and known. Our system is considered fully saturation,
thus Sa, the water saturation equals one, except at the source location. The subscript yon the dissolved concentrations
(c), the NAPL masses (rrij), and the sorbed concentrations (S|) refers to a particular chemical in the NAPL mixture (three
different chemicals are included in this analysis). The given initial concentration is zero at all locations. Thus t=0
represents a time just prior to emplacement of the NAPL source. Exchange between the dissolved phase and the sorbed
phase enters into the state equation directly, while the exchange between the dissolved and NAPL phases is accounted
for in the first boundary condition, which is imposed at location xs. The other two boundary conditions are standard
specified concentration conditions on both edges of the domain.
Figure 3 depicts the dissolved phase transport in this 1 -D point source problem. The problem domain is from 0 to L in the
x dimension. From the location of the point source, xs, concentration profiles develop throughout the domain, becoming
more dispersed over time.
I
I X
1 point source of
mass m (t)
Figure 3. The one-dimensional NAPL point source problem.
If we assume local chemical equilibrium between the aqueous and NAPL phase concentrations of each constituent, the
state equation for the NAPL phase takes on the form:
(9)
= 1,2,3
where c^" (moles of chemical/liter of water) is the equilibrium aqueous concentration with the NAPL phase, also called the
effective solubility. This is the dissolved concentration that enters into the boundary condition for Equation 8.
Q (volume/time) is the flux of water through the source. The effective solubility for each constituent is described by
10
-------
Raoult's Law, which equates the fugacity, the chemical potential at a given temperature and reference state, of two
phases at equilibrium: / f\
.ft _.. ._. . "• ; \Jt., I I
7 = 1,2,3
The effective solubility of each constituent, is related to the pure phase aqueous solubility, csat(/), the mole fraction of the
constituent in the NAPL mixture,^, and the activity coefficient of the constituent within the mixture, y. Note that the mole
fraction depends non-linearly on the relative mass (in moles) of each constituent in the NAPL mixture, m. This
formulation assumes use of the pure phase aqueous solubility relative to the liquid form of the pure phase (as indicated
by the / ). This is necessary since we have used the liquid reference state in the definition of fugacity. Thus at
temperatures under which the chemical would be a gas or solid, the sub-heated or super-cooled liquid aqueous solubility
should be used.
When the equilibrium relationship between the NAPL and aqueous phase is used and the model error is zero, the NAPL
mass at any time is known as long as the initial mass is known. In other words, the initial NAPL masses at the point
source, mo1,mo2, and mo3 completely describe the state m.
If we further assume that the aqueous and sorbed phases are in equilibrium, the sorbed phase state equation is
described by the linear sorption model (5), which relates sorbed concentration to dissolved concentration. The sorbed
concentration, in turn, is related to the state s by:
where ps is the density of the soil itself (in gm/cm3), ne is the effective porosity (dimensionless) and pb = ps(1-ne). By
differentiating equilibrium sorption with respect to time and using (11), we produce:
^sj dcaj
dt ^s" e' d'J dt (12)
Kd ([mg chemical/mg soil]/[mg chemical/liter water]) is the aqueous/sorbed phase partitioning coefficient. Equilibrium
sorption is traditionally incorporated into equation (8) as a retardation factor (R) that modifies the velocity and dispersion
terms of the equation, where R is defined as:
j*d'J (13)
The state equation for the sorbed phase may be re-written as:
dc,
(14)
When equation (14) is substituted into equation (8), the result is the standard form of the 1-D advection-dispersion
equation with linear sorption.
a/ /Rj a* /Rj a*2
The variation in the sorption rates for the different constituents depends entirely on their hydrophobicity as represented
by the partitioning coefficient. Kd is related to the pure phase aqueous solubility of each chemical. We determine the Kd
based on the partitioning of chemical j to the organic matter in the soil:
K - K • f
d-,} om,j J om /* o\
Kom is the organic matter partitioning coefficient (mg chemical/mg organic matter) and fom is the fraction of organic matter
present in the soil. We calculate the Kom using the linear free energy relationship between the aqueous solubility and the
Kom (Schwartzenback et al., 1993):
11
-------
log Kom = -0.75 log cf (/) + 0.44 (1 7)
The variables in the calculation of the partitioning coefficient are all assumed to be known and constant.
Equations 8 through 17 comprise the set of state equations used in our 1-D example. Due to the phase equilibrium
assumptions, the only unknown states in the model are the dissolved concentrations and the only unknown parameters
are the three initial NAPL masses and the source location. The total number of scalar equations to be solved depends on
the spatial and temporal discretization of the problem. For example, the transport equation (8) must be solved for each
node point of the spatial grid and for each time step and each constituent. The ordinary differential equation (ODE)
describing mass transfer from the NAPL source (9) must be solved for each time step and each constituent. If there are
Nx spatial steps, N, time steps, and Nc constituents, the total number of scalar equations is {N, * N * (Nx + 1 )} = Ne. These
Ne equations are solved simultaneously at the discretized times and locations. Solutions are calculated at equally spaced
time steps but varying spatial steps in order to increase the resolution near the source location. The advection/dispersion
equation is solved using a fully-implicit Galerkin finite element solution. The mass balance ODE for the NAPL phase is
solved with a fourth-order Runge Kutta technique. This solution insures that the total mass released from the point
source does not exceed the initial source mass for any constituent. As mentioned above, the state equation solution
generates a set of concentration profiles, in space and time, for each constituent in the NAPL mixture. The amount of
mass present in the NAPL at any given time can also be tracked, but this information is not necessary for the estimation,
since only the initial NAPL mass is considered a parameter in the estimation algorithm.
3.2 Measurement Model
Equipment limitations, variability in sampling techniques, small-scale spatial variability, and human error all contribute to
the deviation of a measurement from the "true" value predicted by the state equation solution of the preceding section.
Therefore we assume that the vector of field measurements used for estimation (z) are corrupted by measurement noise:
j = U,3
where con| is normally distributed with a mean of zero and a variance C . The subscript n is the measurement index, and
Nz is the total number of measurement instances (times and (ocations).ln this formulation, the measurement error will be
multiplicative as cnj becomes small relative to zs(the half saturation constant), and it will be additive as cnj becomes large
relative to zs. This insures that measurements do not take on negative values, since concentration is a positive value in
reality. For the 1-D problem, we simplify this equation by assuming zs is equal to zero. Thus the measurement model
reduces to an additive error term on the forward model:
3.3 Estimation Equations
The goal of the estimation algorithm is to determine the values of the initial NAPL saturations and (in the 1-D case) the
location of the NAPL point source. The best estimate of the spatially discretized states is the one that minimizes the
discrete version of the generalized least squares criterion given earlier:
+-[m-m]TCm1[m-m]
2L J mL J (20)
.
12
-------
The first part of this performance index is a standard weighted least-squares term that measures the difference between
measured dissolved concentrations zand the corresponding model predictions M (c). The second and third terms force
the estimates to stay close to the prior estimates of the initial masses (m^) and the source location (xs). The hats on the
masses, source location, and dissolved concentrations indicate that they are estimates. This formulation assumes that
all constituents are measured at the same times and locations and neglects model error.
Equation (20) may be minimized by setting the variation (or total differentiation) of J equal to zero and solving
explicitly for the unknown variation in a. This requires that the explicit concentration solution be linearized about the
best available estimate (obtained on iteration k):
where a is the vector of estimates:
,(a-6c)
or
(21)
«t =
m,
m.
The new a value (A ) which minimizes (20) is obtained from the iterative Gauss-Newton algorithm (Tarantola, 1987):
"t+i
ok
+ c:
ok
Ok
(a*-a)
(22)
Ok
This method uses the first derivative of the measurement predictions with respect to estimated parameters, as well as an
approximation of the second derivative. When no model error is included in the problem and an explicit expression for the
state equation is substituted into the measurement operator, this method is identical to the variational method (described
in Section 5). Applications of Gauss-Newton solutions to subsurface problems may be seen in Reid (1996) and Sun
(1997).
With each iteration k of the algorithm, new values of the states and the sensitivity derivatives are calculated. The
derivatives required in the linearization are calculated numerically with a finite difference approach:
dc _ c(a + da) - c(a)
9a e
(23)
where<5a,. = eStj (8|( = 1 if i=j and 0 otherwise) j=1,2,...Na, and e is some very small value (relative to a). The sensitivity
derivatives expressed in (22) are matrices of dimension Nz by Na (i.e., number of measurements by number of
parameters). Since the derivatives are calculated from differences, each set of derivatives requires that the state
equations be once for each parameter. Thus each iteration of the estimation requires Na+1 solutions. Equation 22 is
solved iteratively until the values of the parameters converge (i.e., the parameter estimates stop changing significantly).
In this problem we use a convergence tolerance of 0.01 normalized log moles on the initial mass parameters and one grid
node for the location of the source.
State estimation theory also produces expressions for the updated, or posterior, covariances for the parameters and
states (Schweppe, 1973). These equations evaluate the uncertainty associated with the final estimates and represent
13
-------
the Bayesian Cramer-Rau bound (or lower bound) on the estimate uncertainty when the prior probability density pa(
-------
For the purposes of our problem, the NAPL source could be envisioned as residual LNAPL in a region of a moving water
table (i.e., the NAPL is trapped in the unsaturated zone near the water table and then the water table moves above it).
The chemical properties necessary to the estimation are presented in Table 2.
Table 2. Chemical Constants Used in the Simuiated Problem
Benzene Toluene 0-Xylene
Molecular weight 78.1 (g/mol) 92.1 (g/mol) 106.2 (g/mol)
Density 0.88 (g/cm3) 0.87 (g/cm3) 0.88 (g/cm3)
Aqueous solubility 0.023 (mol chem/L water) 0.0056 (mol chem/L water) 0.0017 (mol chem/L water)
(C;'(/)at25°C)
Kom(calculated) 46.8 (mol chem/mol org) 134.1 (mol chem/mol org) 323.0 (mol chem/mol org)
Retardation (calculated) 5.13 (unitless) 12.8 (unitless) 29.5 (unitless)
The forward model also requires some physical constants in order to calculate the states (dissolved concentrations).
Again, these values are assumed to be known and constant and were chosen based on common aquifer properties.
They are presented in Table 3.
Table 3. Physical Constants Used in the Simulated Problem
Ground-water velocity (v) 0.1 m/day
Porosity (ne) 0.3
Dispersion coefficient (D) IrrrVday
Bulk density of soil (p) 2.65 gm/cm3
Fraction organic matter in soil (fm) 0.01
The problem discretization is designed to produce a Courant number (v*dt/dx) of 1 and a Peclet number (v*dx/D) of 1
near the source. This yields a spatial discretization (Ax) of 10 meters (near the source) and a temporal discretization of
100 days. The model was run over a domain of 220 spatial steps (Nx) and 100 time steps (Nt). The simulated "true"
results are presented in Figure 4. This figure shows the concentration of each chemical at the source (normalized by its
aqueous solubility) as well as the concentration profile over the entire domain for three different times. The effects of the
different solubilities are readily seen in this figure. The source concentration shows the results of competitive dissolution,
and the domination of the most soluble constituents. In the concentration profiles, the peak of each constituent separates
overtime due to different retardation coefficients; the more soluble (less hydrophobic) constituents are retarded less than
the more hydrophobic (less soluble) constituents. Thus the varying solubilities of the constituents have a chromato-
graphic effect, separating the concentration profiles of the constituents from each other.
Given the parameters of this simulated problem, the total initial NAPL mass at the source is 458 kg. The total volume of
this mass, calculated from the relative proportions of each constituent, is 0.54 cubic meters. With the given discretization
and porosity, the total volume of voids in a grid cell near the source is 6 cubic meters. If the NAPL were distributed evenly
throughout the cell, the saturation would be 0.09, on the low end for residual NAPL in the saturated zone. Clearly the
saturation value is dependent entirely on the assumed NAPL distribution (e.g., if the source were distributed over a bulk
volume of two cubic meters the saturation would be 0.9). These calculations are merely intended to provide a perspective
on the magnitude of the sample problem. The next section discusses the application of the estimation algorithm to this
sample problem.
15
-------
00
CD
o o o o
UOUB4U90UOO P9Z!|BLUJOU
CO
P9ZJIEIUJOU
2
Q.
c
'o
Q.
1
3
E
"
Q.
3
O
"0
0)
§
O)
16
-------
Section 4
One-Dimensional Results
This section discusses the results of the estimation algorithm described in Section 3. In general, the algorithm accurately
predicted the location and mass of the NAPL point source. We used the algorithm both to explore the feasibility of a
larger NAPL estimation problem and to examine the efficacy of different measurement strategies. Although the 1-D
problem does not represent the real complexities of a multi-dimensional problem, it may be viewed as a rough solution
to source estimation along a streamline. Based on the success of this algorithm, we have formulated a multi-dimensional
distributed source problem (presented in Section 5.)
4.1 Measurement Strategies
The results of the estimation algorithm are entirely dependent upon the chosen measurement strategy; i.e., where and
when samples are taken. At a field site there is likely to be only one measurement design, but with a simulated problem,
we can test the algorithm for different strategies. The advantage of the 1 -D problem is that it is small enough to enable
repeated runs of the estimation algorithm, providing a way to determine the accuracy and cost of estimates obtained with
different sampling configurations. To explore this capability, we have created 16 different measurement strategies.
There are four strategies each for sets of 5 (Strategies 5a-d), 10 (Strategies 10a-d), 30 (Strategies 30a-d), and 60
(Strategies 60a-d) measurements. Each measurement consists of a dissolved concentration sample for each of the
three contaminants (i.e., benzene, toluene, oxylene). The sampling strategies were not chosen in a rigorous optimal
fashion, but rather to illustrate the tradeoff between accuracy and cost for a range of different sampling alternatives. The
strategies were varied to explore the importance of the total number of samples, the sampling locations and the sampling
times. The sixteen strategies in this analysis are shown in Figure 5.
For each strategy, measurements were derived from the measurement presented in Section 3. The additive random
measurement error (co) was assumed to be normally distributed with a mean of 0.0 and a variance (a^) of 0.05 in units
of normalized concentration. In cases where the measurement error was negative and greater than the true value at that
location, the measurement was given a value of zero, thus preventing physically impossible negative concentrations.
The measurement units are normalized (non-dimensionalized) aqueous concentrations. The concentrations are
normalized by the pure phase aqueous solubility of each chemical. Since the solubilities of benzene, toluene, and
oxylene differ by orders of magnitude, non-normalized concentrations could not be displayed on the same plot.
The estimates obtained for each measurement strategy are derived from a particular set (or realization) of the random
problem inputs (initial masses, source location, and measurement errors). The observed performance of the estimator
(as measured by differences between "true" and estimated values) depends, in part, on the particular realization
selected. An extended Monte Carlo (multi-replicate) analysis would be necessary to rigorously test the effects of different
measurement strategies. Nevertheless, the single replicate results given here provide a feeling for the effect of sampling
configuration on estimation accuracy.
4.2 Parameter and Model Fit
Figure 6 shows the error between the true and estimated parameter values (normalized by the true values) for each
measurement strategy. In general, all the strategies do a good job of estimating the initial NAPL mass and source
location. The worst estimates have an error of half the true value, and the best are nearly perfect fits. Of the four
parameters being estimated, the location of the source is usually the least accurate. The mass of the benzene (the most
soluble chemical) is usually estimated better than any other parameters. In these results, it is clear that the number of
measurements is not the most important factor in the success of a sampling strategy. For example, in this particular
realization of measurements, Strategy 5a (with only 5 measurements) outperforms many of the other strategies. In
contrast, Strategy 10a has twice as many measurements, but performs more poorly than 5a. This effect is the result of
measurement spacing; i.e., the 5a measurements are spaced at 50-meter increments while the 10a measurements are
spaced at 10-meter increments. Thus 5a covers a distance of 250 meters while 10a covers 100 meters. The larger
17
-------
o
m
o
o
CD
O
0
o
o
CM
O
i i
o
o
CD
O
O
O
0
CM
O
<•*
O
O
CD
O
O
O
0
CM
O
O
0
CD -^
CD
0 0)
^^ f—
•* -B
o
o c.
o co
w to
TJ
O
nomom inomoin inomoin inomoin
M CM i- i- CM CM i- ••- CMCMi-T- CM CM T- i-
.^-^
O
O
CD
O
O
O
O
CM
O
^jj^jc
8
CD
O
O
O
O
CM
O
O
8
0
o
0
o
CM
O
-
*•*
O
S
o o5
O p
CD
o
o c
O CO
w to
'•o
o
nomom inomoin inomoin inomoin
N CM •«- •>- CM CM -i- i- CM CM T- T- CM CM i- -i-
1
O
O
CD
O
O
O
0
CM
O
1
O
O
CD
O
0
o
o
CM
O
I
CD
O
O
O
O
CM
O
+ +
o
S
0 ft
ol
o eg
^ o5
T3
0
nomom inoinoin inomoin in o m o m
N CM T- T- CM CM *- ••- CM CM •«- i- CM CM T- i-
*
*
*
*
*
O
O
CO
O
o
0
o
CM
0
i
O
o
CD
8
8
CM
O
t
§
o
o
0
o
CM
0
H
o
o
CO —~
CO
8 |
CD
O
o c
O CO
CM ft
a
.nomoin moinoin inomoin inomoin
N CM i- i- CM CM T- 1- CM CM -i- T- CM CM •>- i-
(sjeaA) awn (sjeaA) aiuii (sjeaA) aujji (sjeaA) ami.;
O O O
in T- co CD
Figure 5. Map of measurement strategies for the estimation algorithm. Each star indicates a concentration mea-
surement taken for each chemical constituent. The strategy number indicates the number of measure-
ments, varying from 5 to 60.
18
-------
o.s
S 0.4
i 0.3
ro.2
' 0.1
Parameter Estimates for Strategy 5a
rr
0.5
50.4
!0.3
'0.2
'0.1
Parameter Estimates for Strategy 5b
benzene toluene o-xylene location
benzene toluene o-xylene location
0.5
; 0.4
.30.3
ro.2
' 0.1
Parameter Estimates for Strategy 5c
Parameter Estimates for Strategy 5d
0.5
;0.4
J0.3
ro.2
10.1
benzene toluene o-xylene location
Parameter Estimates for Strategy 10a
o1—
0.5
;0.4
i 0.3
ro.2
'0.1
benzene toluene o-xylene location
Parameter Estimates for Strategy 10b
_
benzene toluene o-xylene location
benzene toluene o-xylene location
Parameter Estimates for Strategy 10c
:0.4
.§0.3
'0.1
benzene toluene o-xylene location
0.5
:0.4
.3 0.3
a
Parameter Estimates for Strategy 10d
0.2
0.1
benzene toluene o-xylene location
Figure 6. Fit of estimated parameters to true values. The absolute value of the difference between the estimate and
the true value for each is normalized by the true value of the parameters.
19
-------
Parameter Estimates for Strategy 30a
Parameter Estimates for Strategy 30b
0.5
0.4
0.3
0.1
n
i 1 . i 1
0.5
i
=50.4
§0.3
|0.2
0.1
n
-
-
I I
benzene toluene o-xylene location
benzene toluene o-xylene location
Parameter Estimates for Strategy 30c
Parameter Estimates for Strategy 30d
0.5
0.4
0.3
0.2
0.1
n
|^-'~ .(
•';-'Sf~^f.
- -- V
--V-1"-
£-&' .-'„':-
•-£j?f-~-^s '
^0.5
=3-0.4
I °'3
lo.2
-§
0.1
n
'
,
benzene toluene o-xylene location
benzene toluene o-xylene location
0.5
;0.4
!0.3
Parameter Estimates for Strategy 60a
0.2
10.1
benzene toluene o-xylene location
0.5
; 0.4
iO.3
.'0.2
Parameter Estimates for Strategy 60b
benzene toluene o-xylene location
Parameter Estimates for Strategy 60c
r
Jfl.3
|0.2
B 0.1
n
-----
0.5
;0.4
J0.3
ro.2
'0.1
Parameter Estimates for Strategy 60d
benzene toluene o-xylene location
benzene toluene o-xylene location
Figure 6. (continued)
20
-------
spacing better characterizes the concentration profile from the simulated NAPL source. In both strategies, measure-
ments are taken at a single time (19 years). Similar comparisons may be made among all sixteen strategies from the
results shown in Figure 6.
Figures 7 and 8 present the fits of the model predictions to the measurements for Strategy 60a, and Figures 9 and 10
present the same information for Strategy 60b. Each of these strategies samples two wells at many times, and many
locations at two times, for a total of 60 measurements. The figures present four concentration profiles, two each over
distance and time, for each strategy. The predicted model is displayed with a solid line (for each chemical), while the true
model is represented by a dotted line. Measurement locations and times are marked with stars. The predicted model
matches the true model more closely for Strategy 60a than for Strategy 60b. The difference in the results of the two
strategies is primarily due to the location of the wells; the samples in Strategy 60a do a better job of capturing the
contaminant profile than do those of 60b. The wells in Strategy 60b do not sample the concentration profiles for toluene
and o-xylene at all, and never intercept the peak for the benzene. In contrast, the wells in Strategy 60a sample across
the peak concentration of all three constituents. This assessment is seen clearly in Figure 11, which superimposes the
measurement strategies (white stars) on the true concentration profiles for each chemical. The better placement of
sampling wells in strategy 60a is also reflected in the convergence of the estimation algorithm. Under strategy 60a the
estimator converged within five iterations while strategy 60b had not converged within tolerance after 15 iterations.
4.3 Evaluation of Uncertainty
Figures 14-15 present the prediction error (model prediction minus true value) for the Strategy 60a and 60b profiles. The
error is normalized by the posterior concentration standard deviation. If these prediction errors were normally distributed,
we would expect approximately 95% of them to fall between +2 and -2 (two standard deviations of a unit normal
distribution). For the most part, this is what we observe. Moreover, the prediction errors obtained for our problem are
very persistent over both time and space. This pattern reflects the effect of the three initial NAPL masses on dissolved
concentration at all subsequent times and locations. Note that where little contamination is present (either in the
predicted or true models) the prediction residuals do exhibit a more scattered pattern. This is most apparent in Figure 19
at early times. As a quantification of the persistence, we can evaluate the approximate lag one concentration correlation
coefficients (p^ for the four distance and time profiles shown in Table 4. The measure of how correlated a residual is with
its nearest or lag one (in space or time) neighbor is generally calculated as the ratio of the square of the covariance,
cov(c,,c2), to the product of the variances of sample (Jenkins and Watts, 1968). In Table 4 these coefficients are given
by:
3c
1
i = 1,2,3 (25)
Note that we have used the derived posterior in the numerator of this expression. The derived correlation coefficients
confirm what is apparent in the figures: there is a high positive correlation among the normalized residuals in the distance
profiles, and the correlation is slightly stronger in the Strategy 60a residuals than in the Strategy 60b residuals.
In large part the success of the estimation algorithm depends on how well the linearization of the forward model
approximates the true model. For this particular problem, the forward model is not only non-linear, but is also not
monotonic. While this characteristic will not necessarily hamper the effectiveness of the algorithm in estimating states
and parameters, it does affect the accuracy of the calculated posterior covariances. Sensitivity analyses of forward
model to the estimated parameters suggest that the local linearization may not necessarily be valid over the range of
parameter values covered by the posterior covariance. Since the posterior covariances are based on this linearization,
the calculated covariances may not be a good representation of the uncertainty in the estimates.
4.4 Evaluation of Measurement Strategies
In Figure 16 we compare the tradeoff between estimation accuracy and sampling cost for different measurement
strategies. The accuracy performance measure used in this figure is the average normalized posterior standard
deviation of the parameters:
*> } (26)
A value of one indicates almost no improvement over the prior, and a value of 0.5 represents an average of 50 percent
improvement on the prior. Note that JA is independent of the measurements themselves, since it depends only on the
sensitivity derivatives at measurement locations and times.
21
-------
Predicted Model and Measurements at 11 years
benzene
toluene
o-xylene
50
100 150 200
distance (meters)
250
300
Predicted Model and Measurements at 19 years
50
100 150 200
distance (meters)
250
300
Figure 7. Measurement residuals for strategy 60a versus distance. The dotted lines show the "true" model value
while the solid lines are the predicted model based on the estimates.
22
-------
Predicted Model and Measurements at 200 Meters
10 15
time (years)
20
25
Predicted Model and Measurements at 250 Meters
10 15
time (years)
20
25
Figure 8. Measurement residuals for strategy 60a versus time. The dotted lines show the "true" model values while
the solid lines are the predicted model based on the estimates.
23
-------
Predicted Model and Measurements at 11 years
benzene
toluene
o-xylene
50
100
150 200
distance (meters)
350
400
Predicted Model and Measurements at 19 years
benzene
toluene
o-xylene
100
150 200 250
distance (meters)
300
350
400
Figure 9. Measurement residuals for strategy 60b versus distance. The dotted lines show the "true" model values
while the solid lines are the predicted model based on the estimates.
24
-------
Predicted Model and Measurements at 350 Meters
1
0.9
0.8
I 0.7
£
v 0.6
o
8o.5h
•o
0>
.M 0.4 h
tO
f o.3i-
c
0.2-
0.1 -
0
0
benzene
toluene
o-xylene
10 15
time (years)
20
Predicted Model and Measurements at 400 Meters
1
0.9
0.8
I 0.7
&
-f-»
§0.6
o
§0.5
T3
CD
.^ 0.4
CO
| 0.3
c
0.2 -
0.1 -
0L
0
benzene
toluene
o-xylene
10 15
time (years)
Figure 10. Measurement residuals for strategy 60b versus time. The dotted lines show the "true" model values while
the solid lines are the predicted model based on the estimates.
25
-------
?
(SJ313U1) 3DUB;Slp (SJ3i3Ul) 33UBISJP (SJ9)3U|) 33UB4.SJP
euezueq euenjo^ eua|Ax-o
Figure 11. Concentration profiles for each chemical with superimposed measurement strategies. The white stars
show the measurement locations and times for strategies 60a and 60b.
26
-------
Normalized Prediction Errors at 11 Years
o
2
CO
CO
1 1
g
o
T3 Oi
o> ^
Q.
T3 *
CD .
N '
i "1
o
c
-2
! ! ! ! ! ! ! ^ f
*
**** * ***,
**** *
^^" _ ^. -jjt— __ ^t- 5(fe_J^-j( a^. ^ Jf" j^_ -^r^*' "^ ^-"^ --^ -3 ^-^— z£ ^ St-"^-^ 4- *— $ -W- jfe- _.
* "* -#-^1 ^iK'^'ii/^Jt-*.**'*-*
f*****^* ^*******
r* * ****
******
-
1 1 1 1 1 1 1
benzene
toluene
o-xylene
1 E H!h* * M~i 8"^
* 5
-
i i
0 50 100 150 200 250 300 350 400 450 50
distance (meters)
Normalized Prediction Errors at 19 Years
3
•esiduals
-•• IN;
— i
O
o
^ 0
£ *
Q. 1
1
CO -1
0
c
-2
o
1 1 1 J 1 1 1
#**£*********#
* * W-
' ^ -A- ^ 3t(
* *** *****
* * ** * *****
o -?fe_ _»jf it *-* -4^- ^fe-i*-A ^u x -^ ^^r-* -^ ^- y * ^-.^-a^ . ^ *_^_ ^ __ 9fc-^_
* ^^^ * :-^ * ^
' * It * * ^
* ************
-
*
-
1 1 1 1 I 1 1
i i
benzene
toluene
o-xylene
*
£ jfe 3t-3£-*£ *$• jf ^ ^_i t
*^ * il
((*********"
-
-
0 50 100 150 200 250 300 350 400 450 500
distance (meters)
Figure 12. Prediction error normalized by the posterior covariance (a+) versus distance for strategy 60a.
27
-------
Normalized Prediction Errors at 200 Meters
10
15
time (years)
20
CO
1 1
2
o
T3
(D
N
O
c
*
*
Normalized Prediction Errors at 250 Meters
benzene
toluene
o-xylene
10 15
time (years)
20
25
Figure 13. Prediction error normalized by the posterior covariance (cry) versus time for strategy 60a.
28
-------
Normalized Prediction Errors at 11 Years
CO
CO
0)
£
1 °
0)
T3
CD
•— —1
"5
o
-2
* *
*
o I ,_
0
*
**
* *
*
benzene
toluene
o-xylene
*********
50 100 150
200 250 300
distance (meters)
*
350 400 450 500
Normalized Prediction Errors at 19 Years
"* *
UJ
«
13
***
£ 1 -
c
CD
CD
3
CO
CD
CO
O
C
L_
k~*—******-~-#v
* * * *
**
*
*
*
benzene
toluene
o-xylene
* *
*• -*•
.**:
*,**$*
,,
**
**"*
_**
-2 -
-3
0 50 100 150 200 250 300 350 400 450 500
distance (meters)
Figure 14. Prediction error normalized by the posterior covariance (o~+) versus distance for strategy 60b.
29
-------
Normalized Prediction Errors at 350 Meters
_cn
to
£ 1
c
CD
-------
Table 4. Lag One Correlation Coefficients for Normalized Prediction Errors
11 years 19 years 200 (350) meters
Strategy 60a
Benzene
Toluene
o-xylene
Strategy 60b
Benzene
Toluene
o-xylene
0.88
-0.51
-0.64
0.88
-0.51
-0.64
0.89
-0.26
-0.50
0.90
-0.26
-0.50
0.99
0.98
0.98
0.99
0.70
0.76
250 (400) meters
0.99
0.88
0.90
0.99
0.60
0.70
Performance of Measurement Strategies
0.8
<0
I
TJ
CO
O
0-6
d>
N
§0.4
CD
2
CO
0.2
10c
5c
5d
5b
5a
10d
10b
10a
30d
*
30b
INFERIOR STRATEGIES
60d
30c
60b 60c
* *
PARETO OPTIMAL STRATEGIES
30a
-*
60a
—*
i
50 100 150 200 250
cost of measurement strategy (thousands of dollars), J
300
350
Figure 16. Comparison of algorithm performance for varying measurement strategies.
31
-------
The total cost of each sampling strategy depends on the number of sampling locations (N^) and the number of times a
sample is taken (Nrt):
Jc=awNa+atNtt <27>
where as = $100 per sample chemical-suite analysis and aw = $10,000 per each new well installation required.
Assessing the sampling strategies with these two measures allows identification of so-called Pareto optimal strategies.
In this context, Pareto optimal sampling strategies are startegies which cannot be improved upon without sacrificing
either accuracy or cost. Applying this criterion, we see that the Pareto optimal strategies are 5d, 10d, 30a, and 60a. All
other strategies are inferior to these — either incurring higher cost without any associated reduction in uncertainty or
more uncertainty without any associated reduction in cost. The worst strategy in terms of performance is 10b, with no
improvement on the prior. Table 5 shows the marginal cost of reduction in uncertainty, relative to Strategy 5d, for the four
optimal strategies.
Table 5. Marginal Cost of Uncertainty Reduction (relative to Strategy 5d) for Pareto Optimal Sampling Strategies
Strategy Marginal cost of uncertainty reduction / percent reduction
5d —
lOd $6,490
30a $22,190
60a $49,530
As mentioned earlier, the 16 sampling strategies were chosen to illustrate a range of different sampling options. What all
four optimal strategies have in common is that they sample at multiple times—from 5d which samples one well five times
over the course of 2 years, to 60a which samples two wells 10 times each. Results obtained for the set of strategies
displayed in Figure 16 reveal several insights about sampling design:
• Location
Perhaps the most obvious, and important, factor in determining success in estimation is whether the
sampling is conducted near the contamination (i.e., whether the measurements contain any information
about the source). This was seen to be the primary difference between Strategies 60a and 60b, but may
also be seen in the dominance of 5d over 5c, 10a over 10b, and 30a over 30b. These pairs of strategies only
differ in fortuitous well placement. The wells of 60b, 10b, and 5c miss most of the contaminant concentration
profiles.
• Characteristic Scale
Another important factor is whether the sampling configuration captures enough of the concentration profile
to characterize the source. Since the estimation algorithm can use information about how the NAPL
chemicals interact with each other, measurements that show a significant difference in the relative
concentration of each chemical produce better estimates. Therefore the spatial coverage of the strategy is
very important. A good example of this is the dominance of Strategy 5a over 10a. Strategy 5a, with half as
many measurements, covers roughly twice as much distance as Strategy 10a and produces a much higher
reduction in uncertainty. This factor is also seen in the dominance of 5a over 5b and 5b over 30a.
• Multiple Times
Another way of getting information about the changes in relative chemical concentration is sampling over
multiple times. In general, it is less expensive to sample an existing well than to install a new one, so
sampling at multiple times (rather than multiple locations) will always be preferable in terms of cost. To be
useful, the sampling times must be far enough apart to capture the interaction among the different NAPL
constituents. All of the Pareto optimal strategies sample at multiple times. Another example is Strategy 10d,
which improves over Strategy 5a by adding a second sampling time, even though the spatial coverage of
sampling is smaller.
32
-------
It should be noted once more that these results are based on a single realization of the NAPL source problem. A different
source (i.e., with different initial masses in a different location) could yield different relative performance results for the
sampling configurations examined in this section. To thoroughly explore the characteristics of these sampling strategies,
the estimation algorithm should be evaluated for a number of different realizations. This was not done in the present
study because the computational cost associated with running a large ensemble of optimization problems for many
different sampling strategies is quite large. However, preliminary experiments with a small number of different
realizations confirm the general sampling design principles outlined above.
33
-------
Section 5
Multi-Dimensional Problem Formulation
The distributed source problem takes a decidedly different form than the 1 -D point source problem of Section 3. Instead
of a fictional non-dimensional point source, we are concerned with source masses that exist in a specified volume.
Furthermore, sources may potentially exist at all locations in the problem domain (i.e., a distributed source). This invites
expression of the source as a continuous property. Figure 17 illustrates the development of this idea. At the pore scale
(Figure 17a), residual NAPL mixtures of mass (Mn) occupy a given fraction of the total volume (9n) or a fraction of the total
void space (Sn). _ Mnpn
0=-
s=-
VT
(28)
At a larger scale (i.e., multiple pores, Figure 17b) we can repeat this calculation of 9 or Sn with a new residual NAPL
mass and new total volume. It is clear that the volumes Figure 17a and b have very different residual NAPL saturations.
If we average the NAPL content over too large a volume, we lose descriptive detail (e.g., overestimate NAPL content in
areas where NAPL is not present and underestimate it where NAPL is present). Therefore, to go from the microscopic
or pore scale to the continuum or macroscopic scale (e.g., on the order of centimeters), we utilize the concept of a
representative elementary volume (REV) (Bear, 1979). We then say that the volume fraction of NAPL (Sn) is a
continuous property, with values at any mathematical point in the problem domain (Figure 17c). Using this continuum
approach, we derive state equations for the system by imposing mass balance constraints on a hypothetical control
volume in Section 5.1.
The second half of this section describes our development of the estimation algorithm for the multi-dimensional
distributed source problem. Whereas in the 1 -D problem we estimated four parameters (i.e., the initial point source mass
and location), in the multi-dimensional problem we are estimating a number of states that vary over space and time. The
state variables for this problem are the residual NAPL saturation for each chemical, the dissolved concentration of each
chemical in the ground water, and the water saturation. Additional unknown parameters are the macro-mass transfer
b) Multiple Pores
c) Continuum
saturation defined at every point m continuum
Figure 17. Definition of NAPL saturation as a macroscale property.
34
-------
rate coefficients, |i (described in section 5.1). As in the 1-D problem, the states and parameter estimation uses
measurements of dissolved chemical concentration in combination with the physical/chemical model described by the
state equations.
5.1 State Equations
The general description of multiphase fluid behavior may be derived from a statement of mass balance among phases
and species. We consider a control volume (Figure 18) in which there are four potential phases (K = vapor [v], solid [s],
aqueous [a], and NAPL [n]) and some number of species within those phases (i = air [v], water [a], soil [s], and N.
chemical constituents 0]). For each species (j) in each phase (K), the total mass (MKl) present in the control volume (VT)
at any given time is related to the fluxes into and out of the control volume as well as the exchange of mass among
phases within the control volume. We consider an advective flux (FKj) and a dispersive flux (dK) in each direction, for
each species in each phase. Each species may be present in any of the phases. The exchange of mass from one phase
to another is represented by lpici (i.e., exchange of mass from phase p to phase K).
Figure 18. Mass balance of each species (i) in each phase (K) in control volume.
The mass balance for this control volume is expressed mathematically (adapted from Abriola and Pinder 1985):
where:
3.K.I
(p
(pKeK
K .
. - v
K .
(29)
volume averaged density of phase K in the control volume
volume fraction of phase K in the control volume
mass fraction of species i in phase K in the control volume
velocity of phase K
dispersion coefficient of species i in phase K
total inter-phase exchange of species i to phase K in the control volume
inter-phase exchange of species i from phase p to phase K in the control volume
35
-------
This statement assumes no external sources or sinks of mass (e.g., no pumping, no chemical decay). Equation (29) is
completed with the following constraints:
=1
K=/
i-l
1=1
These constraints insure that mass is conserved over species and phases. The first requires that in a given phase, the
mass fractions of each species present sum to the total mass of the phase. The second equation constrains the volume
fractions such that the sum of the volume of all phases is not greater than the total volume. The third and fourth
constraints apply to the inter-phase exchange rates. The total change of mass of a given phase must be equal to the
contributions (i.e., change in mass) of each species to that phase (third constraint). Finally, for a given species, mass is
conserved over all phases or the total change in mass is zero (fourth constraint).
When equation (29) is summed over each phase and each species and the constraints (30) are applied, the result is a
set of coupled partial differential equations describing all of the mass fractions and volume fractions in the system. In the
current formulation there are 4*(3+N.) mass fractions (e.g., with 3 chemical constituents there would be 24 mass
fractions). However, our initial description of the problem allows some simplifications. Since we are considering a NAPL
source in the saturated zone only, we neglect the vapor phase and air species (T)V.= TIM= 0). We make the further
simplifications that the soil (i = s) and water (i = a) species do not change phases (lKa = lks = 0), and only exist in the solid
and aqueous phases, respectively (t\as= iln8 = t|sa= rina= 0). Practically, this means we are neglecting water in the
sorbed phase and water in the NAPL phase.' We are now left with only 2 + (3*N() mass fractions and 3 volume fractions:
mass fractions
T)8S soil in the solid phase
r|aa water in the aqueous phase
T|SJ chemical j in the solid phase
TJJJ chemical j in the aqueous phase
T|nJ chemical j in the NAPL phase
for j = 1,2,...,^ chemicals
volume fractions
6a aqueous phase
6S solid phase
9n NAPL phase
Species Balances
Conservation of mass in the control volume leads to a continuity equation for each of the (2 + j) species. Each of the
species balances takes advantage of the last constraint in (30): the inter-phase exchange of a species over all phases
is equal to zero. As stated earlier, we assume that soil is only present in the solid phase. We also assume an immobile
medium, i.e., there is no movement of the solid phase and vt = 0 , D$l = 0. Thus the mass of soil in the control volume
doesn't change over time. The soil species balance becomes:
36
-------
The water species is also assumed to be present only in a single phase. The species balance is:
= 0 (32)
The chemical species, on the other hand, are present in all three phases. Since our problem statement describes a
residual or immobile NAPL, vn = 0 , Dni = 0. The only movement of chemicals occurs in the aqueous phase. For each
chemical j the species balance is:
Phase Balances
Whereas the species balances assure that all of the mass of a particular species is conserved, the phase balances
assure that the total mass of all species in a phase is equal to the mass lost or gained from that phase. Each of these
balances takes advantage of the fact that the sum of all dispersive fluxes in a phase must equal zero. In addition, each
phase balance invokes the first constraint that the sum of all the species mass fractions within a phase is equal to one.
In the solid phase, soil and chemical species are present:
— (p8) + ?I .=0 (34)
-} \rsvs/ £j*s,a,j "
Again, an immobile medium is assumed. In the aqueous phase, water and chemical constituents are present:
,y =0 (35)
In the NAPL phase, only the chemical constituents are present. Note that this phase is immobile:
=o (36)
n,atj
7=1
Equations (31) through (36) require additional information in order to be solved. This information is how chemical mass
moves among the three phases, since water and soil mass do not change phase. While expressing the species and
phase balances, we use the aqueous phase as a reference state and assume that all mass transfer occurs through this
phase. Therefore there is no direct mass transfer between the solid and NAPL phases (lsnj= I f = 0). We are left
chemical mass only transfers between the aqueous-solid (las) and lsaj) phases and the aqueous-NAPL (lna) and lanj)
phases.
Interphase Exchange
As in the 1-D formulation in Section 3, sorption is considered to be a linear equilibrium process. We believe this is a
reasonable assumption since 1) both soil and water exist everywhere throughout the model domain and 2) the time
scales of sorption (on the order of seconds) are fast relative to both advective transport time scales and the time steps
of the numerical model (on the order of tens of days). The equilibrium solid/water partitioning relationship is expressed
in terms of the aqueous concentration of chemical j (caj). Assuming a dilute aqueous solution (i.e., the contaminant
contributes a negligible amount to the density of the aqueous phase), caj may be expressed as a function of the mass
fraction:
n =-±L (37)
^ Pa
where p is the density of water.
37
-------
The transfer of chemical mass (j) from the aqueous to the solid phase is then described as:
(38)
The partitioning coefficient (Kd) is in units of (mass of j on soil per mass of soil) / (mass of j in water per volume of water).
This mathematical statement assumes that the contribution of the chemical to the mass and volume of the solid phase
is very small compared to the contribution of the soil (i.e., ps = pss and 9s = 6SS). A more detailed discussion of the solid
water partitioning coefficients may be found in Section 3.
The exchange of mass between the NAPL and aqueous phases is a rather different process. Whether or not we assume
a chemical equilibrium between these two phases, it is unlikely that dissolved chemical concentrations would reach
chemical equilibrium over a large (i.e., measurement) scale. The low concentrations relative to equilibrium are attributed
to irregular distribution (see Figure 19 below), heterogeneous flow patterns, dilution, and/or non-equilibrium mass
transfer. Recently, a number of studies have supported the idea that most of the seemingly non-equilibrium dissolution
of the NAPL phase may be explained by heterogeneity of the NAPL distribution and in the porous media (Sorens et al.,
1998; Unger et al., 1998). To capture this effect, we consider the dissolution/precipitation process at two scales: the
local scale (e.g., centimeters) and the model scale (e.g., tens of meters).
Area of
equilibrium
NAPL Source
O-
Grid node in model
c,j not at equilibrium
with NAPL phase
-O
Figure 19. Representation of non-equilibrium NAPL dissolution at the model scale.
Figure 19 illustrates our conceptualization of dissolved phase concentrations at the numerical model grid scale. In a
given model grid cell, areas nearest to NAPL sources will be at chemical equilibrium (effective solubility), however areas
with no NAPL present will have dissolved concentrations less than equilibrium values. The overall concentration in the
cell will not necessarily reach equilibrium. At the beginning of this section, Figure 17 described a continuum of NAPL
saturation (i.e., Sn defined at every point in space). Figure 19 relates to the discretization of this property for modeling
purposes. As will be discussed in Section 5.1.5, we must solve the state equations numerically on a discrete grid.
Furthermore, in order to model processes at a field scale (hundreds of meters to kilometers), the model discretization
must be large compared to scales of chemical equilibrium. If variability in NAPL distribution like that shown in Figure 1 9
is averaged over a model cell, the model will underestimate the NAPL saturation or even predict no NAPL present at all.
Our model formulation addresses this variability by representing mass transfer at the model scale as a kinetic process
but representing dissolution at the local scale as chemical equilibrium.
In sum, our formulation applies different models of dissolution at two scales. At the local scale, we assume chemical
equilibrium as modeled by Raoult's Law:
(39)
38
-------
This is the same form as presented in Section 3. This form of competitive dissolution directly relates NAPL saturation (Sp
to the dissolved concentrations. The equilibrium assumption holds as long as the time scales of phase exchange are
long enough to cover distances of interest. At typical ground-water velocities chemical equilibrium between the dissolved
and NAPL phases would be expected to persist at the centimeter scale.
At the model scale, however, we no longer assume chemical equilibrium. A number of processes other than chemical
kinetics may contribute to the non-equilibrium concentrations seen at field sites:
• The time scales needed for diffusion through an entire grid cell are large. The size of the model grid cells
alone means that it takes much longer to come to chemical equilibrium throughout an entire cell. For a
typical organic chemical, the aqueous molecular diffusion coefficient is 10'5 cm2/s. The coefficient roughly
provides a time scale for chemical equilibrium since the spread (ax) associated with diffusion is ax = (2Dt)1/2.
At this rate, molecules travel approximately 1 cm in each direction from the source in a day. An area the size
of a model grid cell (e.g., 1 meter) would take 6 days to reach equilibrium in the absence of advective
transport.
• Heterogeneity in the porous medium and subsequent heterogeneity in the velocity field lead to mixing and
dilution of saturated contaminant concentrations.
• In some areas of the porous medium, mass transfer from the NAPL phase may be diffusion limited (e.g., in
dead-end pores or other areas where advection is not the contaminant transport mechanism) rather than
advective transport limited.
• Heterogeneity in the distribution of NAPL phase throughout a model cell will also affect the extent of dilution.
In addition, the rate of dissolution of the NAPL is related to the interfacial (or contact) area between the NAPL
and aqueous phases. A more evenly distributed source would be expected to dissolve more quickly than a
highly heterogeneous source.
Most researchers currently use a first-order equation to model aqueous-NAPL phase mass transfer. The question, then,
is how to develop a mass transfer rate coefficient (Miller et al., 1998). To date, most experimental work towards this end
involves small-scale (a few centimeters) column experiments of a single NAPL component in a homogeneous medium.
Based on the uncertainty in current approaches to modeling NAPL dissolution, especially as applied to large scales, we
have chosen to use a simple model:
Ia,n,j=-»(ce?j-caJ)ea=-In>aJ (40)
We consider jo. a macroscopic mass transfer coefficient (in units of 1/time). A good review of empirical models for mass
transfer rates may be found in Miller et al. (1998). The mass transfer rates are developed as functions of Reynolds
number (velocity, viscosity), Sherwood number (diffusion, chemical rate coefficient), Schmidt number (diffusion,
viscosity), varying grain size measurements, NAPL saturation, uniformity indexes, contact area between NAPL and
water, and various assumptions about the shape of NAPL blobs. These relationships are far too complicated for our
purposes in this analysis. We assume that m is a simple function of NAPL saturation:
(*> »P (41)
where b is the estimated parameter describing the power law dependence of \i on Sn. This mass transfer coefficient only
explicitly varies with the NAPL saturation, and allows the model to quantify the effects of dilution and heterogeneity on
equilibrium dissolution. Guiguer and Frind (1994) use a similar form of the mass transfer coefficient as a generic
correlation to NAPL saturation. A recent study by Unger et al. (1998) found that the Sn based dissolution model of
Guiguer and Frind yielded comparable results to a more complicated model that included a Reynolds number and grain
size correlation (Mayer and Miller, 1996). Given these results and the lack of established dissolution models in the
literature, we believe that the simple Sn-based model is the best choice for this research.
These mass balances and constitutive relationships complete the set of state equations. Now we proceed to make a
number of additional assumptions that further simplify our mathematical description of NAPL phase exchange and
transport.
Simplification of Equations
The first species balance equation (31) tracks soil mass, which is not particularly the focus of this analysis. Since we are
assuming that the soil does not move or exist in phases other than the solid phase, this equation becomes trivial; any
change in mass of the solid phase is due to the change in chemical mass sorbed. Rearranging the equation yields:
39
-------
(42)
Note that the problem formulation assumes all sorbed chemical mass has come from the aqueous phase. The rate of
change of sorbed mass is equal to the interphase exchange expression, las .
Equation (32) may be simplified by applying the dilute solution assumption: the mass of chemical j in the water is very
small compared to the mass of water itself. The fraction of water in the aqueous phase is approximately equal to one.
When this term is constant, the diffusive flux term of the equation drops out. By further assuming an incompressible fluid,
the density of the aqueous phase is approximately constant. Then the water species equation becomes a sum of fluxes
(water mass balance) over the control volume, recognizable as the groundwater flow equation:
In the discussion below, equation (35), the aqueous phase balance, will reduce to the same expression.
Rearranging the chemical species balance (33) for each chemical constituent leads to the transport equation for each
chemical constituent in the aqueous phase. Again we assume a dilute solution, which allows the mass balance of
equation (33) to be expressed in terms of concentrations:
+ M + V.vc-V.e = 0 (44)
The interphase exchange relations (38) and (40) by definition describe the change in mass in the sorbed phase and the
NAPL phase. Substituting in these expressions and applying the standard definition of a retardation factor (this
expression is given in Section 3):
^^ + V • (O.V.C.J ) - V • (S0DaJ Vc., )
fc (45)
The final form of our transport equation comes from applying the chain rule to the first and second terms of the above
equation, using the relation from equation (43), and rewriting in terms of water saturations rather than volume fractions:
--Vv-DVc+e -'
This equation is subject to the homogeneous boundary and initial conditions:
caj(£l,t) = Q on the domain boundary Q
cfl>,(x,0) = 0
Note that the effective solubility of chemical j in the transport equation is a function of the residual NAPL saturation. Other
parameters in equation (46) are also a function of states: velocity (va), retardation (Rj), and dispersion (Daj), are all
functions of the water saturation (Sa), and the macro mass transfer coefficient (u.) is a function of the NAPL saturation
(Sn). The solid phase balance (34) reduces to a trivial expression of the interphase exchange of chemical into the sorbed
phase (the same expression that the soil species balance produced):
v si \ ^°"'J r (47)
^p"*1-".)—'- ( '
40
-------
The next phase balance (35) produces the standard ground-water flow equation after applying the assumptions: a) the
chemical constituent in the aqueous phase is a dilute solution, b) water is incompressible, and c) the interphase
exchange terms are much smaller than the other terms in the expression and may be neglected. The first two
assumptions are the same as those used in equation (43), and allow pa to be treated as a constant. After applying the
third assumption:
(48)
the water phase balance becomes identical to equation (43), the ground-water flow equation. In terms of saturations, the
flow equation is:
The flow equation is completed by substituting in the generalized form of Darcy's Law. This yields the non-linear
Richards equation.
(50)
at
where h is hydraulic head, K is the saturated hydraulic conductivity, and kr is the relative permeability. In a system of
constant boundary conditions, no (or little) recharge, and a nearly incompressible medium, the ground-water flow field
will only change as the water saturation changes. As NAPL dissolves, the total volume occupied by NAPL will decrease
and the volume occupied by the water will subsequently increase. However, the change in Sa over time is likely to be
small since 1) the residual NAPL volumes are small compared to S and 2) NAPL dissolution occurs slowly — on the order
of tens to hundreds of years. Thus over a model time step (e.g., hundreds of days), Sa will be approximately constant.
For this reason, we use a steady-state ground-water flow equation in our model. Meanwhile, small changes in S can
result in large changes in the relative permeability. The empirical relationship between water saturation and permeability
is usually represented as a power law (e.g., Van-Genuchten, Brooks-Corey):
kr(Sa)~[f(neSa)]d (5D
where f and d are fitting parameters. Then the steady-state ground-water flow equation with a variable hydraulic
conductivity that depends on the NAPL dissolution model through the dependence on Sa:
V(Kkr(Sa)Vh)=Q (52)
The boundary conditions applied to (52) are:
In other words, these are simple no-flow (Q2) or constant head (Q,) boundaries. More complicated boundary conditions
could be applied to the model at a marginal computational cost. Although there are technically no initial conditions
necessary to solve (52), it does require a value for Sa to initialize the relative permeability.
The final phase balance (36) describes the NAPL dissolution process. Writing the equation in terms of saturations and
applying the dissolution interphase exchange model, we produce:
From this point, we further simplify the equation by combining it with Raoult's law and coverting moles of NAPL into
fractional saturations using the identity:
m ^neSiPjmw» (54)
41
-------
After some simplification, the NAPL dissolution equation becomes:
PA£
-iu
= 0
(55)
This equation requires an initial condition of the fractional NAPL saturation:
and is completed with the volume averaged quantities:
(56)
where p( is the density of pure phase chemical j and mw) is its molecular weight.
Since we are interested in the dissolution of each individual component from the NAPL mixture, we can rewrite (55) with
the time derivative expressed in terms of S. rather than Sn:
(57)
PA mW, "'J
Numerical Solution of State Equations
Along with boundary and initial conditions, equations (46), (52), and (57) are the final state equations for our multi-
dimensional distributed source NAPL problem. Together they produce a system of (2j+1) coupled partial differential
equations summarized in Table 6. Note that the dissolution and transport equations both depend on the dissolved
concentration and the residual NAPL saturation. Additionally, all of the state equations contain parameters that also
depend on the states.
We solve the set of equations using a Galerkin finite element method (over linear triangular elements) combined with a
fully-implicit finite difference temporal discretization. It would be easy to discuss at some length the relative merits of
different numerical methods (and there are many such books). For this analysis, suffice it to say that the finite element
approach offers flexibility for coupled problems and minimizes problems of numerical dispersion. The Galerkin technique
(as opposed to the co-location or subdomain methods) is the most commonly used in ground-water modeling (Celia and
Gray, 1992). This method uses the basis functions of the trial space (in our case linear triangular basis functions) as the
weighting functions in the residual minimization and is equivalent to a variational principle, when one exists for the given
problem. As applied to our state equations, each state (S^ Sa, ca) is discretized at grid nodes and represented with the
trial basis functions. Other parameters (D, v, R, jo.) are considered constant over elements. The spatially discretized (but
continuous in time) state equations are presented in Table 6.
42
-------
Table 6. State Equations for Multi-dimensional NAPL Source Problem
Final State Equations
V « «• "y
as,
foreachj
Pj
(
-------
In summary, the steps of the solution to the state equations are:
1. Initialize the model with a velocity field, initial S., initial S , and initial c =0.
j a a,j
2. Calculate the effective solubility, mwn, pn, Sn, and |i for the current NAPL saturations.
3. Solve the dissolution model for NAPL saturations at new time = t +1.
4. Solve the ground-water flow equation using the water saturations for time = t +1 and update the velocity field
for time = t + 1.
5. Return to step 2 and recalculate parameters dependent on NAPL saturation. Iterate until convergence of the
NAPL saturations for time = t + 1.
6. Input NAPL saturations, water saturation, and effective solubility (for time = t +1) from the dissolution model
into the transport model.
7. Solve the transport equation for the new dissolved concentration (time = t + 1) of each chemical.
8. Return to step 2 using concentration value and saturations for time = t + 1. Iterate until NAPL saturations
and dissolved concentrations for time = t + 1 converge.
9. Advance one time step. Saturations and concentrations from step 8 become values for time = t. Return to
step 2.
10. Continue algorithm (steps 2-9) until final time step.
5.2 Measurement Model
The measurement model for the multi-dimensional problem is the same as the one used in the 1 -D problem. In this initial
approach to the problem, we posit the linear measurement model of dissolved concentrations (state variables) only:
z -c (x t } + (O o = 12N (58)
^ p,j t-a,pjVAp>' p' *"pj r ^ff-f—^z
where p is the measurement index. Again the measurement error oop is considered a normally distributed random
variable with zero mean and variance CM.
5.3 Estimation Equations
The performance index for the multi-dimensional distributed source problem may be expressed in the following spatially
discrete continuous time form:
Ni N
T TC T/ T P "P -\T fr,p^—\r P *P -,
./= I Xk/-
-------
Although we have constrained the state estimates to be near some prior values, we need to take additional care to insure
that the estimator will not produce physically impossible values for the states. Both the water and NAPL saturations must
be positive values ranging between zero and one. In fact, since these states are volume of void fractions, at any given
location they must sum to one. We address this requirement by transforming the saturations to ensure that they are
positive and adding a constraint to the performance index to ensure that physical volume limitations are respected. Both
the water and the NAPL saturations are transformed by introducing the new variables:
s?
(60)
S'= In-
iref
which are log saturations. Sref is a normalization factor equal to 1. The performance and the model state equations are
now written in terms of the new transformed states.
We use a Lagrange multiplier approach to augment the performance index so that the various physical constraints on the
problem can be incorporated into the least-squares minimization. We adjoin (2j + 1) constraints to the performance
index. For each chemical j, we want to constrain the estimate of the dissolved concentration and the NAPL saturation by
the corresponding state equations. These nonholonomic constraints are identical to the spatially discrete form of the
state equations shown in Table 6, with the addition of a model error term. The final constraint is an algebraic expression
that relates the water saturation to the NAPL saturations. Combined with the transformation on the state variables
discussed above, this constraint forces the saturation estimates to sum to one. The performance index then becomes:
AT,
NJ
I
7=1 />=!
_z
caj]
1
-lr P
j=\
; [S'ao-S'o]
jo
Jo
C 1(M)(T)dfa/T
vc J
f f wJ
J J
f f 1)^(0 C"1 (/,t)L,(T)<*«/
J J L v n J
(61)
Nj
21
7=1
dt
K(/)exp[S'0(0]ca?y-(0 + i)w(0
exp[S'0(0]
NJ
I exp[S'y-(/)]
7=1
-1
•dt
45
-------
The hats (A) over each variable represent the estimate of that variable.
When J is at its minimum value, the state and parameter estimates will be the best fit to the measurements, model, and
prior statistics. We accomplish this minimization with a variational approach. This approach allows us to elegantly
minimize a least-squares type cost function while conforming to a number of constraints. The variational method is
somewhat analogous to the differential method, exploring the effect of a displacement on a function or functional. The
variational operator (8) represents an infinitesimal perturbation of the function, while the derivative (d) represents an
infinitesimal change in the function with respect to an independent variable. These differences seem subtle, but result in
important properties of the variation: 1) at the stationary point (minimum or maximum) each partial variation of the
function must equal zero and 2) if the boundaries of the function are prescribed, the variation with respect to the
boundaries must also equal zero. (Daley, 1991). To find the stationary value of our performance index, we take the first
variation (SJ) and set each partial variation equal to zero. The first variation of J is represented symbolically as:
^ X^, -L.^JfiVj-^ ^ XV , d-/ XV ,
OC„ • ~r ~ r~ CXj T" / ~~" Q3 ; T ~ :CAJ __ "1"
+
t=T
(62)
^A oJ ™ _,/ oJ ,.. oJ ~ oJ ~ oJ 01 oJ f. oJ ™.
t=T
The set of coupled non-linear equations resulting from setting 5J = 0 are called the Euler-Lagrange equations. They
include: 1) a set of update equations for the parameters, 2) a set of state equations, 3) a set of adjoint equations, and 4)
a set of model error equations. Sometimes only the adjoint equations are called the Euler-Lagrange equations, but we
will use this term to refer to the entire group of equations. Each of the partial variations in (62) corresponds to an Euler-
Lagrange equation:
• The variations with respect to the parameters (including the state initial conditions) lead to parameter update
equations;
• The variations with respect to the adjoint variables lead to the state equations;
• The variations with respect to the states lead to the adjoint equations;
• The variations with respect to model error lead to model error equations;
• The variations with respect to the state terminal values lead to the terminal conditions on the adjoint
variables.
We evaluate these equations at the current estimate of all variables. After some manipulation, the final Euler-Lagrange
equations for this problem are given in Table 7.
The parameter update equations depend on estimates of the states, the adjoint variables, and the sensitivity derivatives
(evaluated at current estimates). The model error equations depend on the adjoint variables. The adjoint equations are
subject to a terminal conditional and are thus called the backward equations. The state equations (subject to initial
conditions) are called the forward equations. The adjoint equations depend on the adjoint variables, the states, and the
sensitivity derivatives. It is important to note that while the state equations (Table 7) are non-linear in the state variables,
the adjoint equations are linear in the adjoint variables. Due to the coupled and non-linear nature of the Euler-Lagrange
equations, solving them is not simple. The equations must either be de-coupled or solved iteratively (or both). There are
many approaches to solving two-point boundary value problems (like the Euler-Lagrange equations); here we limit the
discussion to three methods.
Probably the simplest method for solving the Euler-Lagrange equations is the direct approach. In this method, the adjoint
equations are solved using an estimate of the states, and then the state equations are solved using the adjoint variable
estimates. The adjoint solution sweeps backwards (starting from the terminal condition) and the state solution sweeps
forward from the initial condition. Each iteration begins with the estimates from the previous iteration and this process
continues until the convergence criterion is met.
Another solution approach is to use a Kalman filter formulation to transform the two-point boundary value problem into an
initial value problem. This approach is based on an assumed linear relationship between the states and the adjoint
variables, which is used to derive a Kalman gain. The jumps calculated from the Kalman filter are used to solve the
adjoint equations, and then the adjoint variables are used to solve the state equations. The development of the Kalman
filter can be quite computationally expensive since it requires propagation of the state covariance matrix through time
(Bennet, 1992).
46
-------
Table 7. Set of Euler-Lagrange Equations for Multi-dimensional NAPL Source Problem
Parameter Update Equations
0 J
0 J
Sb
?' =S' -C f
>flo "° s'°°\l\
-^ -C /
JO
*
Model Error Equations
Adjoint (Backward) Equations
6t
; - S]
t=T
= 0
f ro/Tr
Afl=-exp[5J
«r
= 0
T a//rx
for each j
foreachj
a...
47
-------
The representer approach de-couples the Euler-Lagrange equations with a series of linear expansions of the state
estimates, parameter estimates, and adjoint variable estimates—called representer expansions. The representer
solution is, in effect, a solution to the linearized problem. In general, this method is more stable than the direct approach
(Bennet, 1992) but requires more computational effort. Where each iteration of the direct approach requires one solution
of the state equations and one solution of the adjoint equations, the representer solution must solve for the representers
for each measurement (N times). An advantage of this method is that uncertainty evaluation of the estimates (posterior
covariances) are readily derived from the representer expansions. This approach is used by Bennet (1992) and Valstar
(1998).
In this research, we used direct approach to solving the Euler-Lagrange equations. Each iteration consists of two model
runs, one forward run and one adjoint run. The general sequence of the procedure follows:
1. Initialize the estimator.
(a) Set the initial model error to zero.
(b) Set the initial parameter values and initial conditions on the states to the prior values.
(c) Solve the state equations (see Section 5.1.5) to produce initial estimates for the states.
2. Using the state estimates, solve the adjoint equations to produce estimates of the adjoint variables.
3. Using the adjoint variable estimates, calculate the model error variables.
4. Solve the update equations for the initial conditions and parameters using the model error estimates and the
adjoint variable estimates.
5. Solve the state equations based on the values from the update equations.
6. Repeat steps 2 through 6.
7. Estimator is done when the estimates of the initial conditions and parameters converge.
This completes the development of the multi-dimensional NAPL estimation problem. The algorithm described in this
section may be conditioned with computer-simulated data or actual field data of varying source configurations. The next
logical steps of this research involve formulating the prior statistics for the states and testing this algorithm on different
sets of data. This future work is discussed further in Section 6.
48
-------
Section 6
Conclusions
This research explores the possibility of using readily available measurements of dissolved contaminant concentrations
to characterize a NAPL source. We approach this task as a state estimation problem, viewing both the dissolved
concentrations and NAPL saturations as states that vary over time and space. As an inverse/state estimation method,
our research takes the following form:
• Our method is based on a stochastic description of subsurface variability. States and parameters are
considered to be random variables characterized by means and covariances.
• For the multi-dimensional problem, we have selected a set of state equations that model NAPL dissolution
as a first order kinetic process. This produces a mass flux source term for the ground-water transport
equation.
• We use a generalized least-squares performance index that fits predicted state values to measurements
while including regularization terms derived from prior statistics for uncertain model parameters and model
errors.
• In the 1-D analysis, we minimize the performance index with an iterative Gauss-Newton algorithm that
produces estimates of the states and parameters.
• In the multi-dimensional analysis, we use a variational approach to minimize the least-squares performance
index.
Although results of the multi-dimensional distributed source problem are not presented in this report, we feel that our
formulation shows promise based on successful applications of similar methods (Reid, 1996; Sun, 1997; Valstar, 1998).
The 1 -D point source problem presented in Sections 3 and 4 demonstrates the feasibility of our hypothesis. Additionally,
this formulation is computationally cheap enough to allow numerous estimation runs and thus exploration of varying data
gathering strategies.
Our approach is based on a number of assumptions and limitations that have been discussed throughout this report.
Many of these simplifications were made to ease the computational burden of the estimator. Others were necessitated
by a lack of precedent or time to develop new methodologies. Almost all of these limitations represent areas for further
. research. The primary limitations are:
• Simplification of physical processes: In many places we make assumptions that allow us to represent
physical processes in more convenient terms. For example, we assume a steady-state description of
ground-water flow that neglects recharge. Additionally, our model of NAPL dissolution is not specifically
related to porous media properties, nor does it account for diffusion limited zones (e.g., dead-end pores)
within the media. As is often the case with issues of model complexity, it is not clear how crucial these
simplifications are to an accurate description of physical reality. Further analysis that compares our models
to more sophisticated (and more computationally intensive) models would provide a great deal of insight into
the nature of these assumptions.
• Omission of potentially significant processes: The models we have constructed do not address a
number of physical processes. For example, processes of chemical decay and biodegradation are not
included in the ground-water transport mode. Inclusion of these processes may or may not improve our
representation of reality, but it would most certainly increase the computational burden on the estimator.
• Form of prior information/covariance structure: The random variables in this analysis (states and
parameters) are all described with assumed prior statistics. In some cases (e.g., NAPL saturation) we have
no data on which to base the prior. Furthermore, our methodology assumes a Gaussian structure to the
random variables. There is no particular reason to believe that this is a correct assumption.
• Application to synthetic data: Testing the estimation algorithm on synthetic data has both advantages and
disadvantages. In a synthetic problem, unlike a field problem, we know the answer that the estimator is
trying to find. This allows us to assess the accuracy of the estimator. On the other hand, there is no
49
-------
guarantee that the synthetic data, which after all is generated by the same simplified models we use in the
estimator, bear any relation to real field data. One way to address this concern is to generate synthetic data
with a more complicated forward model than the one used in the estimator. Another method is to test the
estimator with an appropriate set of real field data.
• Simple incorporation of model error: Model error enters our performance index as an additive term. A
more sophisticated representation of process error might relate it to the model states in such a way that it
represents a physical process (e.g., unknown recharge). This could lead to model error terms that enter
the state equations in a non-additive fashion.
Over the near future, we hope to complete the multi-dimensional problem outlined in this report. This work involves
developing and implementing the computer code for the estimator. We plan to test the estimator on a synthetic two-
dimensional problem. Once the estimator is complete we will be able to further test the sensitivity of the algorithm to
different problems. The eventual goal of this work is to expand the estimator into three dimensions and apply the
estimator to field measurements from an actual contaminated site. We believe this method will provide much insight into
characterization of subsurface NAPL contamination.
50
-------
Appendix A: Table of Variables and Operators
,
eff
f
om
F
i
I.
Variables
as cost per chemical sample analysis
aw cost per well installation
cti concentration of i in phase K
effective solubility for j in water
pure phase saturated solubility for j in water
state covariance function matrix
measurement error covariance matrix
model error covariance function matrix
prior parameters covariance matrix
dispersion coefficient of i in phase K
fraction organic matter
advective flux
index for species (n=chemical, a=water, v=air, s^soil)
interphase exchange of i from phase (3 to phase K
total interphase exchange of i into phase K
index for NAPL chemical
performance index
performance measure for accuracy
performance measure for cost
estimation iteration index
relative permeability
hydraulic conductivity
sorption coefficient for chemical j
organic matter partitioning coefficient
end-point of one-dimensional problem domain
moles of j in NAPL phase
mass in VT of j in phase K
measurement index
ne effective porosity
Na number of parameters
NK number of phases
N number of chemical constituents
N. number of species (air, water, soil, chemicals)
Ne number of equations
Nt number of time steps
Nx number of spatial steps
N number of measurements (N = N + N )
Z v Z ZZ ZK
N^ number of measurement locations
Nrt number of measurement times
Q flux of water through the one-dimensional source
R retardation factor
s sorbed concentration (mg chemical/mg organic matter)
Sa water saturation
Sn NAPL fractional saturation of chemical j
S residual NAPL saturation
K
Kd
K
,
M .
KJ
Units
[$]
[$]
[M/L3]
[M/L3]
[M/L3]
-varies-
-varies-
-varies-
-varies-
[LVT]
[M/M]
[M/T/L3]
-none-
[M/T/L3]
[M/T/L3]
-none-
-none-
-none-
[$]
-none-
-none-
[L/T]
[(M/L3)/(M/M)]
[(M/L3)/ (M/L3)]
[L]
[moles]
[M]
-none-
[LVL3]
-none-
-none-
-none-
-none-
-none-
-none-
-none-
-none-
-none-
-none-
-none-
[M/M]
[LVL3]
[LVL3]
[LVL3]
[LVT]
51
-------
[l/L]
[T]
-varies-
[L/T]
[V/V]
[U/V]
[V/V]
[L3]
[L]
-varies-
-varies-
-varies-
-none-
[mol/mol]
-none-
[M/T/L3]
-none-
[M/M]
[1/T]
[LVL3]
[L]
Ss specific storage
t time dimension (days)
u model forcing function
VK velocity of phase K
V volume of NAPL in control volume
n
V volume of pore spaces in control volume
V volume of voids in control volume
v
VT total volume
x spatial dimension
xs location of point source in one-dimensional problem
y model state vector
z measurement vector
a model parameter vector
P phase index (NAPL=n, vapor^v, aqueous=a, solid=s)
X mole fraction of constituent in NAPL mixture
Y activity coefficient of constituent in NAPL mixture
$ dispersive flux
K phase index (NAPL=n, vapor=v, aqueous=a, solid=s)
T|KI mass fraction of i in phase K
u, macro mass transfer rate coefficient
6K volume fraction of phase K
pK volume averaged density of phase K
Pj lag one correlation coefficient
0 standard deviation
i) process/model error vector
ro measurement error vector
Operators
/( state equation operator
^ state equation for chemical mass in the dissolved phase
^/ state equation for chemical mass in the sorbed phase
j4n state equation for chemical mass in the NAPL phase
/ RHS operator of forcing terms from solved state equation
£* boundary conditions for state equation
V LHS operator from implicit form of solved state equation
y forward operator
0 model error operator
V initial conditions for state equation
1X measurement operator
UNITS NOTATION: $ = dollars, M= mass, L - length, T - time, mol = moles. The units notation "varies'
indicates a general variable that may represent specific variables of differing units.
-varies-
-varies-
-varies-
-varies-
-------
References
Abriola, L.M. and Pinder, G.F. 1985. "A multiphase approach to the modeling of porous media contamination by organic
compounds 1. Equation development." Water Resour. Res. 23(1): 11-18.
Anderson, M., Johnson, R., and Pankow, J. 1992. "Dissolution of Dense Chlorinated Solvents into Ground Water: 1.
Dissolution for a Well-Defined Residual Source." Ground Water 30(2):250-256.
Annable, M.D. 1997. "Use of partitioning tracers for measuring residual NAPL: Results from a field scale test." J.
Environ. Eng. In press.
Banerjee, S. 1984. "Solubility of Organic Mixtures in Water." Environ. Sci. Technol. 18(1984):587-591.
Bear, J. 1979. Hydraulics of Groundwater. McGraw-Hill (New York).
Bedient, P.B. 1994. Ground Water Contamination. Prentice-Hall (Englewood Cliffs, New Jersey).
Bennet, P.B. 1992. Inverse Methods in Physical Oceanography. Cambridge University Press (New York).
Brusseau, M.L. 1991. "Rate Limited Sorption and Nonequilibrium Transport of Organic Chemicals in Low Organic
Carbon Aquifer Materials." Water Resour. Res. 27(6):1137-1145.
Burr, D.T. 1994. "Nonreactive and reactive solute transport in three-dimensional heterogeneous porous media: mean
displacement, plume spreading, and uncertainty." Water Resour. Res. 30(3):791-815.
Carrera, J. 1987. "State of the art inverse problem applied to flow and solute transport problems." In Groundwater Flow
and Quality Modeling. NATO ASI Series, p. 549-583.
Celia, M.A. and Gray, W.G. 1992. "Numerical Methods for Differential Equations: Fundamental Concepts for Scientific
and Engineering Applications." Prentice Hall, 436 pp.
Cohen, R.M. and Mercer, J.W. 1993. DNAPL Site Evaluation. C.K. Smoley Ed. (CRC Press, Florida).
Dagan, G. 1985. "Stochastic modeling of groundwater flow by unconditional and conditional probabilities: The inverse
problem." Water Resour. Res. 21(1):65-72.
Daley, R. 1991. Atmospheric Data Analysis. Cambridge University Press Ed. (Cambridge, England).
El-Kadi, A.I. 1992. "Applicability of sharp-interface models for NAPL transport: 1. Infiltration." Ground Water30(6):849.
Faust, C., Guswa, J., and Mercer, J. 1989. "Simulation of three-dimensional flow of immiscible fluids within and below
the unsaturated zone." Water Resour. Res. 25(12):2449-2464.
Feenstra, S., Mackay, D., and Cherry, J. 1991. "A method for assessing residual NAPL based on organic chemical
concentrations in soil samples." Groundwater Monitor. Rev. Spring 1991:128-136.
Feenstra.S. and Cherry, J. 1988. Subsurface contamination by dense nonaqueous phase liquid (DNAPL) chemicals.
Presented at International Groundwater Symposium, International Association of Hydrogeologists. Halifax, Nova
Scotia, May 1 -4, 1988.
Ginn, T.R. and Cushman, J.H. 1990. "Inverse methods for subsurface flow: a critical review of stochastic techniques."
Stochastic Hydrol. Hydraul. 4(1990): 1-26.
Guiguer, N. and Frind, E.O. 1994. "Dissolution and mass transfer processes for residual organics in the saturated
groundwater zone." Transport and Reactive Processes in Aquifers. Balkema (Rotterdam) 475-480.
Hall, C.W. and Johnson, J.A. 1992. "Limiting factors in ground water remediation." J. Haz. Materials. 32(2):215-223.
Hess, K.M., Wolf, S.H. and Celia, M.A. 1992. "Large-scale natural gradient tracer test in sand and gravel, Cape Cod,
Massachusetts: 3. hydraulic conductivity variability and calculated macrodispersivities." Water Resour. Res.
28(8):2011-2028.
Hoeksema, R.S. and Kitanidis, P.K. 1984. "An application of the geostatistical approach to the inverse problem in two-
dimensional groundwater modeling." Water Resour. Res. 20(7): 1003-1020.
Imhoff, P., Jaffe, P., and Pinder, G. 1993, "An experimental study of complete dissolution of a nonaqueous phase liquid
in saturated porous media." Water Resour. Res. 30(2): 307-320.
53
-------
James, A.I., Graham, W.D., Hatfield, K., Rao, P.S.C., and Annable, M.D. 1997. "Optimal estimation of residual non-
aqueous phase liquid saturations using partitioning tracer concentration data." Water Resour. Res. 33(12):2621-
2636.
Jenkins, G., and Watts, D. 1968. "Spectral analysis and its applications." Holen-Day, Oakland, CA.
Jin, M., Delshad, M., Dwarakanath, V., McKinney, D.C., Pope, G.A., Sepehrnoori, K., and Tilburg, C.E. 1995.
"Partitioning tracer test for detection, estimation, and remediation performance assessment of subsurface nonaqueous
phase liquids." Water Resour. Res. 31 (5):1 201 -1 21 1 .
Kueper, B., Redman, D., Starr, R., Reitsma, S., and Mah, M. 1993. "A Field Experiment to Study the Behavior of
Tetrachloroethylene Below the Water Table: Spatial Distribution of Residual and Pooled DNAPL." Ground Water
31(5):756-766.
Kuiper, L.K. 1 986. "A comparison of several methods for the solution of the inverse problem in two-dimensional steady-
state groundwater flow modeling." Water Resour. Res. 22(5):1 543-1 570.
Lavigne, D. 1990. "Accurate, on-site analysis of PCBs in soil- A low cost approach." Proceedings of HMCRI's 11th
Annual National Conference and Exhibition Superfund '90. Hazardous Materials Control Research Institute,
Washington, D.C. p. 273-276.
LeBlanc, D.R. 1991. "Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts: 2.
Analysis of spatial moments for a nonreactive tracer." Water Resour. Res. 27(5):895-910.
Mackay, D.M. and Cherry, J. 1985. 'Transport of organic contaminants in groundwater." Environ. Sci. Technol.
19(5):384-392.
Mackay, D.M. and Cherry, J. 1989. "Groundwater contamination: pump and treat remediation." Environ. Sci. Technol.
23(6):630-636.
Mackay, D.M. and Wilson, R.D. 1995. "Direct determination of residual nonaqueous phase liquids in the saturated zone
using SF6 as a partitioning tracer." Environ. Sci. Technol. 29:1255-1258.
Mayer, A. and Miller, C. 1996. The influence of mass transfer characteristics and porous media heterogeneity on
nonaqueous phase dissolution." Water Resour. Res. 32(6):1 551 -1567.
Mclaughlin, D. 1995. "Recent developments in hydrologic data assimilation." Review of Geophysics, Supplement. U.S.
National Report to International Union of Geodesy and Geophysics 1991-1994. 977-984.
Mclaughlin, D. and Townley, L.R. 1996. "A reassessment of the groundwater inverse problem." Water Resour. Res.
Mercer, J.W. and Cohen, R.M. 1990. "A review of immiscible fluids in the subsurface: Properties, models,
characterization, and remediation." J. Contam. Hydrol. 6 (1990): 107- 163.
Miller, C.T., Christakos, G., Imhoff, P.T., McBride, J.F., Pedit, J.A., and Trangenstein, J.A. 1998. "Multiphase flow and
transport modeling in heterogeneous porous media: challenges and approaches." Advances in Water Resources.
21(2):77-120.
Miller, C., McNeill, M., and Mayer, A. 1990. "Dissolution of Trapped Nonaqueous Phase Liquids: Mass Transfer
Characteristics." Water Resour. Res. 26(11 ):2783-2796.
Neuman, S.P. 1973. "Calibration of distributed parameter groundwater flow models viels as a multi-objective decision
process under uncertainty." Water Resour. Res. 9(4):1 006-1 021.
National Research Council (NRC) 1990. Ground Water Models: Scientific and regulatory applications. Water Science
and Technology Board. National Academy Press. Washington, D.C.
Pankow, P. and Cherry, J. 1996. Dense Chlorinated Solvents and other DNAPLs in Groundwater: History, Behavior,
and Remediation. Waterloo Press (Ontario, Canada).
Parker, J.C., Katyal, A.K., and Kaluarachchi, J.J. 1991. Modeling multiphase organic chemical transport in soils and
groundwater. EPA/600/2-91/042. U.S. Environmental Protection Agency (Washington, D.C.)
Pinder, G.F. and Abriola, L.M. 1986. "On the simulation of nonaqueous phase organic compounds in the subsurface."
Water Resour. Res. 22(9): 1095-1 195.
Powers, S., Abriola, L., Dunkin, J., and Weber, W. 1994. "Phenomenological models for transient NAPL-water mass-
transfer processes." J. Contam. Hydrol. 16:1-33
Powers, S., Abriola, L., and Weber, W. 1992. "An experimental investigation of nonaqueous phase liquid dissolution in
saturated subsurface systems: steady state mass transfer rates." Water Resour. Res. 28(1 0):269 1-2705
54
-------
Powers, S., Loureiro, C., Abriola, L, and Weber, W. 1991. Theoretical study of the significance of nonequilibrium
dissolution of nonaqueous phase liquids in subsurface systems." Water Resour. Res. 27(4):463-477
Reid, L. 1996. A functional inverse approach for three-dimensional characterization of subsurface contamination. Ph.D.
thesis, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology.
Reid, D.B. and Mclaughlin, L.B. 1994. "Estimating continuous aquifer properties from field measurements: the inverse
problem for groundwater flow and transport." In Computational Methods in Water Resources, X. p. 777-784.
Dordrect. Kluwer Academic.
Schwartzenbach, R.P., Gschwend, P.M., and Innboden, D.M. 1993. Environmental Organic Chemistry. John Wiley &
Sons (New York).
Schweppe, F.C. 1973. Uncertain Dynamic Systems. Prentice Hall (New Jersey).
Schwille, F. 1988. Dense Chlorinated Solvents in Porous and Fractured Media. Lewis Publishers (Chelsea, Michigan).
Sleep, B. and Skyes, J. 1989. "Modeling the transport of volatile organics in variably saturated media." Water Resour.
fles.25(1):81-92.
Sorens, T., Sabatini, H., and Harwell, J. 1998. "Effects of flow bypassing and nonuniform NAPL distribution on the mass
transfer characteristics of NAPL dissolution." Water Resour. Res. 34(7):1657-1673.
Sun, C. 1997. A stochastic approach for characterizing soil and groundwater contamination at heterogeneous field
sites. Ph.D. thesis, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology.
Sun, N-Z, 1994. Inverse Problems in Groundwater Modeling. Kluwer, Dordrecht, Netherlands.
Sun, N-Z and Yeh, W.W-G 1992. "A stochastic inverse solution for transient groundwater flow: parameter identification
and reliability analysis." Water Resour. Res. 28(12):3269-3280.
Tarantola, A. 1987. Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation. Elsevier,
Netherlands.
Unger, A., Forsyth, P., and Sudicky, E. (1998) "Influence of alternative dissolution models and subsurface heterogeneity
on DNAPL disappearance times." J. Contam. Hydrol. 30(1998):217-242.
U.S. Environmental Protection Agency, 1991. Dense Nonaqeous Phase Liquids. EPA Groundwater Issue Paper, EPA/
540/4-91/002. U.S. Environmental Protection Agency, Washington, D.C.
U.S. Environmental Protection Agency, 1992. Estimating Potential for DNAPL Occurrence at Superfund Sites. EPA
Quick Reference Fact Sheet, EPA Publication 9355.4-07FS. U.S. Environmental Protection Agency, Washington,
D.C.
U.S. International Trade Commission, 1991. Synthetic organic chemicals, United States production and sales, 1990.
USITC Publication 2470, Washington, D.C.
Valstar, Johan. 1998. Stochastic Characterization of Subsurface Contamination at Industrial Sites. Ph.D. thesis,
Department of Civil Engineering, Technical University of Delft. The Netherlands.
Wilson, J.L., Conrad, S., Mason, W.R., Peplinski, W., and Hagan, E. 1990. "Laboratory investigation of residual liquid
organics." EPA/600/6-90/004, Robert S. Kerr Environmental Research Laboratory, Ada, Oklahoma.
Yeh, W.W.G.1986. "Review of parameter identification procedures in groundwater hydrology: The inverse problem."
Water Resour. Res. 22(2):95-108.
Zimmerman, D.A., deMarsily, G., Gotway, C.A., Marietta, M.G., Axness, C.L., Beauheim, R.L., Bras, R.L., Carrera, J.,
Dagan, G., Davies, P.B., Gallegos, D.P., Galli, A., Gomez-Hernandez, J., Grindrod, P., Gutjahr, A.L., Kitanidis, P.K.,
Lavenue, A.M., McLaughlin, D., Neuman, S.P., RamRoa, B.S., Ravenne, C., and Rubin, Y. 1998. "A comparison of
seven geostatistically based inverse approaches to estimate transmissivities for modeling advective transport by
groundwater flow." Water Resour. Res. 34(6): 1373-1414.
55
-------
U.S. Environmental Protection Agency
Region 5, Library (PL-12J)
77 West Jackson Boulevard, 12th Floor
Chicago, IL 60604-3590
-------
United States
Fnvironmental Protection Agency/ORD
National Risk Management
Research Laboratory
Cincinnati, OH 45268
Please make all necessary changes on the below label,
detach or copy, and return to the address in the upper
left-hand corner.
If you do not wish to receive these reports CHECK HERE[] ;
detach, or copy this cover, and return to the address in the
upper left-hand corner
PRESORTED STANDARD
POSTAGES FEES PAID
EPA
PERMIT No. G-35
Official Business
Penalty for Private Use
$300
EPA/600/R-01/044
------- |