Development of a Data

Evaluation / Decision

Support System for

Remediation of Subsurface

Contamination

United States

Environmental Protection

Agency

Office of Research and

Development

Washington DC 20460

EPA/600/R-01/044

July 2001
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

National Risk Management Research Laboratory

Office of Research and Development

U.S. Environmental Protection Agency

Cincinnati, OH 45268
ii

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 Technology. 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 Protection 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.
iii

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 environmental risks in the future.

The National Risk Management Research Laboratory (NRMRL) is the Agency’s

center for investigation of technological and management approaches for preventing and reducing risks from pollution that threatens human health and the environment. 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 regulations and strategies at the national, state, and community levels.

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. Also, this report contains

the results of the one-dimensional analysis as well as the development of the twodimensional state and estimation equations that are being used in ongoing research. 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 Development to assist the user community.

Stephen G. Schmelling, Acting Director

Subsurface Protection and Remediation Division

National Risk Management Research Laboratory
iv
v

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 concentrations 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 distribution coefficient) and parameters (hydraulic conductivity, dispersion coefficients,

specific storage) are known and certain (not random variables). The state equations 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 onedimensional 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.
vi
vii

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
viii

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 (σ+

y

) versus

distance for strategy 60a ......................................................................................... 27

Figure 13: Prediction error normalized by the posterior covariance (σ+

y

) versus

time for strategy 60a................................................................................................ 28

Figure 14: Prediction error normalized by the posterior covariance (σ+

y

) versus

distance for strategy 60b. ........................................................................................ 29

Figure 15: Prediction error normalized by the posterior covariance (σ+

y

) 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 (κ) in control volume ................... 35

Figure 19: Representation of non-equilibrium NAPL dissolution at the model scale. ................ 38
ix

Tables

Table 1: Parameter Values and Prior Statistics for Simulated Problem................................. 14

Table 2: Chemical Constants Used in the Simulated Problem .............................................. 15

Table 3: 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
x

SI Conversion Factors

English (US) Metric (SI)

Multiply Units by Factor to get Units

Area: 1 ft2 0.0929 m2

1 in2 6.452 cm2

Flow rate: 1 gal/min 6.31 x 10-5 m3/s

1 gal/min 0.0631 L/s

1 MGD 43.81 L/s

Length: 1 ft 0.3048 m

1 in 2.54 cm

Mass: 1 lb 453.59 g

1 lb 0.45359 kg

Volume: 1 ft3 28.316 L

1 ft3 0.028317 m3

1 gal 3.785 L

1 gal 0.003785 m3

Temperature: oF - 32 0.55556 oC

Concentration: 1 gr/ft3 2.2884 g/m3

1 gr/gal 0.0171 g/L

1 lb/ft3 16.03 g/L

Pressure: 1 lb/in2 0.07031 kg/cm2

1 lb/in2 6894.8 Newton/m2

Heating value: Btu/lb 2326 Joules/kg

Btu/scf 37260 Joules/scm
xi

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 comments. 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.
1

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 groundwater 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 byproducts. 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
2

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:

p

n

n V

V S = (1)

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 saturation usually 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 measurements 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
3

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:

∂

∂

= +

=

=

y

t yu yu

y x

y

A G

I

B

,, ,,,

,

,

α αυ

α

α

b gb g

b g bg

b g

t

0

(2)

In this equation y(x, t) is a vector of states, α(x) is a vector of time-invariant parameters, u(x, t) is a vector of known

forcing functions, and ν(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. G is the model forcing

operator. The operators I and B 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 α 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(α, t) of states defined at the Ny nodes of the computational

grid, where α 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.

Figure 1. Conceptual representation of the distributed NAPL source problem.
4

The measurement equation relates field measurements to the states, the parameters, and a measurement error term

at discrete times and locations (n):

z nz nn n =+ = M by N , ; , ,... α ω g 1 2 (3)

where M n is a measurement operator which describes how the measurement z n is related to the states and parameters

and ωn

 is an additive random measurement error. The quantity Mn(y, α) is just the model’s prediction of the measurement

zn , which is obtained at location x n 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 α,ν, and ωn

 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 fit” 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:

J ,,   ,  , 

 

αυ α α

αα αα

υ υ

ω

α

υ

y z y Cz y

C

C

T

T

T

b g bg bg =− −

+− −

+

−

−

−

1

2

1

2

1

2

1

1

1

M M

(4)

The first term of this performance index is a standard weighted least-squares measure of the difference between a vector

composed of all Nz measurements (z) and the vector M(y,α) composed of the corresponding Nz discretized model

predictions. The weighting factor is the measurement error covariance matrix (Cω). 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 (υ), weighted by the model error covariance

(Cυ). 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 estimation (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

NAPL sources. The success of the method is evaluated by applying it to a synthetic problem. Section 2 provides a review

of previous research on NAPL processes and related characterization issues. Section 3 explains the formulation of a 1-

D NAPL 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.
5

Figure 2. Mass exchange among phases: aqueous, vapor, sorbed, and NAPL for chemical j.

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 al., 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 (foc).

water











































soil

soil

air



































ca







cn cv

air

cv







NAPL

sorbed cs

Single Species at Equilibrium

sat

a j a j

n j

sat

s j d a j

sat

v j H a j

c c

c

c K c

c K c

, ,

,

, ,

, ,

1

=

=

=

=

sat

a j j j a j

j

j

j

n j

sat

s j d j j a j

sat

v j H j j a j

c c

m

m

c

c K c

c K c

, ,

,

, ,

, ,

χ γ

χ

χ γ

χ γ

=

= =

=

=

∑

Multiple Species at Equilibrium
6

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:

s j d a j c K c , = , (5)

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:

sat

j j

N

j

j j

eff

a j c m m c j

γ



 



 = ∑

=1 , (6)

where mj is the moles of chemical j in the NAPL mixture, cj

sat is the pure phase solubility of chemical j in water, ceff is the

effective solubility, and mj 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.
7

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, Sr

 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 al. (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, regularization 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
8

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.
9

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 multidimensional 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:

∂

∂ =

∂

∂ =

∂

∂ =

c

m

s

x,t

t

x

x,t

t

x

x,t

t

x

b g b g

b g b g

b g b g

A

A

A

c s

m s

s s

csm

csm

csm

,, ,

,, ,

,, ,

(7)

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 As 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):

∂

∂ = − ∂∂

+

∂

∂ − ∂∂ = c x,t

t

c x,t

t

c x,t

x n

s

t

x 2

j jj

e a

j

c s S

csm bg bg bg ν D b g

2 1 A ,, , (8)

subject to the conditions:
10

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:

dm x ,t

dt c x ,t Q x

j

j s

j

eff

s

b g =− = bg b g =

Am csm ,, ,

123 , ,

(9)

where cj

eff (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

c x ,t f m ,t

c t

c L,t

,m ,x

c x m ,x

j

js j

j

j

j js

j js

b g d i

b g

b g

d i

b g d i

− =

=

=

U

V

|

W

|

=

= =

=

0

0 0

0

0 0

123

,

,

, ,

B

I

c

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 j on the dissolved concentrations

(cj

), the NAPL masses (mj), and the sorbed concentrations (sj) 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.

Figure 3. The one-dimensional NAPL point source problem.

j 1 c (x,t )

j

c (x,t )2

point source of

L

mass m (t) j

xs 0
11

Raoult’s Law, which equates the fugacity, the chemical potential at a given temperature and reference state, of two

phases at equilibrium:

1,2,3

( , )

( , ) ( , ) ( ) ( )

=

= =

∑

j

m x t

m x t

c x t c l c l j

j s

sat j s

j j j

sat

s j

eff

j χ γ γ (10)

The effective solubility of each constituent, is related to the pure phase aqueous solubility, cj

sat(l), the mole fraction of the

constituent in the NAPL mixture,χj

, and the activity coefficient of the constituent within the mixture, γj. Note that the mole

fraction depends non-linearly on the relative mass (in moles) of each constituent in the NAPL mixture, mj

. This

formulation assumes use of the pure phase aqueous solubility relative to the liquid form of the pure phase (as indicated

by the l ). 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:

(1 ) j s, j b e s = c ρ − n (11)

where ρs is the density of the soil itself (in gm/cm3), ne is the effective porosity (dimensionless) and ρb = ρs(1-ne). By

differentiating equilibrium sorption with respect to time and using (11), we produce:

t

c

n K

t

s a j

s e d j

j

∂

∂ = − ∂∂ ,

, ρ (1 ) (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:

s d j

e a

e

j K

n S

n

R ,

(1 ) 1 ρ − = +

(13)

The state equation for the sorbed phase may be re-written as:

∂

∂ = − ∂∂ = st

R nS

c

t

x j

j ea

j

d i1 A b g s s csm ,, , (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.

2

2 ( , ) ( , ) ( , )

x

c x t

R

D

x

c x t

R

v

t

c x t j

j

j

j

j

∂

∂

+

∂

∂ = − ∂∂

(15)

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:

d j om j om K = K ⋅ f , , (16)

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):
12

log K = −0.75log c (l) + 0.44 sat

om j

(17)

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, Nt time steps, and Nc constituents, the total number of scalar equations is {Nt * Nc * (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:

1,2,3

1,2,...

( , , ) ( , )

( , ) ( , ) ,

,

,

, ,

=

=

= + = +

j

n N

x

c x t z

c x t

z c x t

z

n j s n

n j n n s

n j n n

n j n j n n ω M n c, m ω

(18)

where ωn,j 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 locations). In this formulation, the measurement error will be

multiplicative as cn,j becomes small relative to zs (the half saturation constant), and it will be additive as cn,j 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:

n n j t nx n j c nj z , ( , ) = , + ω (19)

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:

Jx c c

C

xx xx

,  ,  

 

 

cm z C z

mm mm

C

s

T

T

m

s sxs s s

b g bg bg =− −

+− −

+− −

−

−

−

1

2

1

2

1

2

1

1

1

M M ω

(20)
13

The first part of this performance index is a standard weighted least-squares term that measures the difference between

measured dissolved concentrations z and 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 (moj) 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 α. This requires that the explicit concentration solution be linearized about the

best available estimate (obtained on iteration k):

    c c  c

= +

∂

∂α

α − = α bg b g α α 0 (21)

whereαà is the vector of estimates:











αk

1

2

3

s k

m

m

m

x

=

L

N

M

M

M

M

O

Q

P

P

P

P

The new α value ( 1 àαk + ) which minimizes (20) is obtained from the iterative Gauss-Newton algorithm (Tarantola, 1987):





 



 − ∂

∂ − +

∂

∂

















+

∂

∂

∂

∂ = + −

−

− −

+ à ( à ) ( à )

à

1

à

1

1

à

1

à

1 α α

α

α

α α α

α α

α

ω

α

α

α

ω

α

k

k

k

k

T

k k

T

k C C C z

c

c

c c c (22)

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:

ε

α δα α

α

c c( ) c( )

j

+ − = ∂∂ (23)

where i ij δα = εδ (δij = 1 if i=j and 0 otherwise) j=1,2,...Na, and ε is some very small value (relative to α). The sensitivity

derivatives expressed in (22) are matrices of dimension Nz by Nα (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 Nα+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
14

the Bayesian Cramer-Rau bound (or lower bound) on the estimate uncertainty when the prior probability density pα(α)

and posterior probability density pα|z(α|z) are Gaussian (McLaughlin and Townley, 1996):

[ ]

δα

α

δ

α α

α α

α

α

α

α

α

ω

α

α

∂

∂ ≈

∂

∂

∂

∂ ≈

















+

∂

∂

∂

∂ ≈

+

− − − +

c

c

c C c C

C c C c C

k k

T

y

k k

T

à à

1

à

1

à

1

(24)

The + on the left-hand-side covariances indicates that these are updated or posterior covariances. Note that these

relationships do not depend on the values of the measurements themselves, only on the estimated sensitivity derivatives

at measurement. The information (or inverse covariance) contained in the estimations (the posterior Cα) is equal to the

sum of the information in the measurements (the first term) and the prior information (the second term). As Cα decreases,

the confidence in the prior increases and the prior parameter information is weighted more heavily than the measurement

information. The posterior covariance of the states depends only on the prior parameter covariance, the measurement

error covariance, and the sensitivity derivatives ∂c ∂α . As Cα is more uncertain, the posterior Cy will also be more

uncertain.

3.4 Simulated Problem

The estimation algorithm (22) could be applied to any set of measurements and prior statistics. In this analysis we have

created a simulated set of measurements to evaluate the algorithm. Although it is also possible to apply the algorithm to

real field data, using synthetic data has a number of advantages: 1) with synthetic data the “true” values are known, thus

it is possible to assess the accuracy of the algorithm; 2) we can generate many synthetic problems to test the algorithm

under a variety of conditions; 3) with a synthetic data set we can test many different sampling strategies; and 4) getting

access to real field data is difficult and time consuming. The disadvantage of using synthetic data is that the model

generating the data exactly matches our model in the estimation algorithm. In reality, we would expect that the forward

model does not capture all of the complex processes existing in nature. However, as a development tool, tests on

synthetic data are still preferable to tests on real data; if the algorithm is not effective with the simulated data set, it is

unlikely to work at a real field site.

First, a set of “true” parameter values was generated from assumed statistical properties for all random variables. Each

variable was assumed to be normally distributed with given mean and variance. This particular problem has four

parameters: the initial log mass of each of three NAPL constituents and the location of the source. The true parameter

values and associated prior statistics are shown in Table 1. Note that the estimation algorithm works with the natural log

of the total initial normalized (unitless) moles [ln ( ) ref

oi

m

m ] of each constituent at the source, mref is taken to be unity. This

insures that the algorithm will produce a positive value for the mass. Table 1 also presents each parameter in more

commonly accepted mass units.

The prior mean and variance were chosen to produce a problem with significantly different amounts of each NAPL

constituent at the source. Note that the prior variance is fairly high (i.e., two standard deviations is 40 kg). The

constituents of the NAPL mixture are assumed to be known. We used a set of chemical properties consistent with

benzene (C6H6), toluene (C6H5CH3), and o-xylene(1,2-(CH3)

2

C6H4). These constituents are LNAPLs that are commonly

studied and often present in mixture with each other. They are also useful for this problem because their pure phase

aqueous solubilities vary by orders of magnitude.

Table 1. Parameter Values and Prior Statistics for Simulated Problem

Prior mean Prior variance Simulated true value

Benzene Benzene Benzene Benzene Benzene initial initial initial initial initial mass ln ( ) ref

o

m

m 1 7.3 (120 kg) 8.5 (400 kg2

) 8.3 (312 kg)

Toluene Toluene Toluene Toluene Toluene initial initial initial initial initial mass ln ( ) ref

o

m

m 2 6.8 (80 kg) 8.4 (400 kg2

) 7.3 (132 (kg)

o-Xylene -Xylene -Xylene -Xylene -Xylene initial initial initial initial initial mass ln ( ) ref

o

m

m 3 5.9 (40 kg) 8.2 (400 kg2

) 5.3 (23 kg)

Location Location Location Location Location of source (xs

) 40 meters 1,600 meters2 69 meters
15

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.

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 2. Chemical Constants Used in the Simulated Problem

Benzene Toluene o-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 l w

sat b g at 25∞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)

Table 3. Physical Constants Used in the Simulated Problem

Ground-water velocity (v) 0.1 m/day

Porosity (ne

) 0.3

Dispersion coefficient coefficient coefficient coefficient (D) 1m2

/day

Bulk density density density density density of soil (ρ) 2.65 gm/cm3

Fraction Fraction Fraction Fraction Fraction organic organic organic organic organic matter in soil (fom) 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 (∆x) 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

over time 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 chromatographic 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.
16

Figure 4. True model outputs for simulated NAPL point source problem.

0 5 10 15 20 25

0

0.80.60.40.2

1

time (years)

normalized concentration

Source Concentration

benzenetolueneo−xylene

0 50 100 150 200 250 300

0

0.80.60.40.2

1

Solute Concentration Profile at 5 Years

distance (meters)

normalized concentrationsource benzene toluene o−xylene

0 100 200 300 400

0

0.80.60.40.2

1

Solute Concentration Profile at 15 Years

distance (meters)

normalized concentrationsource benzene toluene o−xylene

0 100 200 300 400 500

0

0.80.60.40.2

1

Solute Concentration Profile at 25 Years

distance (meters)

normalized concentrationsource

benzenetolueneo−xylene
17

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, o-xylene). 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 (ω) was assumed to be normally distributed with a mean of 0.0 and a variance (σ ω

2 ) 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

o-xylene 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
18

0 200 400 600

5

25201510

distance (meters)

time (years)

6

0

0 200 400 600

5

25201510

distance (meters)

0 200 400 600

5

25201510

distance (meters)

0 200 400 600

5

25201510

distance (meters)

0 200 400 600

5

25201510

a

time (years)

5

0 200 400 600

5

25201510

b

0 200 400 600

5

25201510

c

0 200 400 600

5

25201510

d

0 200 400 600

5

25201510

0 200 400 600

5

25201510

time (years) 10

0 200 400 600

5

25201510

0 200 400 600

5

25201510

0 200 400 600

5

25201510

time (years)

3

0

0 200 400 600

5

25201510

0 200 400 600

5

25201510

0 200 400 600

5

25201510

Figure 5. Map of measurement strategies for the estimation algorithm. Each star indicates a concentration measurement taken for each chemical constituent. The strategy number indicates the number of measurements, varying from 5 to 60.
19

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.
20

Figure 6. (continued)
21

spacing better characterizes the concentration profile from the simulated NAPL source. In both strategies, measurements 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 (ρ1) 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(c1

,c2), to the product of the variances of sample (Jenkins and Watts, 1968). In Table 4 these coefficients are given

by:

∑

−

= +

+ ∂

∂

∂

∂

− =

1

1 1

1

1, 1

1 Nt

t t t

t

T

i

t

i

t

i

c C c

N σ σ

α α

ρ

α

 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:

4

( ) ( ) ( ) ( ) 1 1 2 2 3 3 s s m m m m m m x x

A J

σ σ σ σ σ σ σ σ + + + + + + + = (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.
22

0 50 100 150 200 250 300

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

distance (meters)

Predicted Model and Measurements at 11 years

benzene 

toluene 

o−xylene

0 50 100 150 200 250 300

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

distance (meters)

Predicted Model and Measurements at 19 years

benzene 

toluene 

o−xylene

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.
23

0 5 10 15 20 25

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

time (years)

Predicted Model and Measurements at 200 Meters

benzene 

toluene 

o−xylene

0 5 10 15 20 25

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

time (years)

Predicted Model and Measurements at 250 Meters

benzene 

toluene 

o−xylene

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.
24

0 50 100 150 200 250 300 350 400

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

distance (meters)

Predicted Model and Measurements at 11 years

benzene 

toluene 

o−xylene

0 50 100 150 200 250 300 350 400

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

distance (meters)

Predicted Model and Measurements at 19 years

benzene 

toluene 

o−xylene

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.
25

0 5 10 15 20 25

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

time (years)

Predicted Model and Measurements at 350 Meters

benzene 

toluene 

o−xylene

0 5 10 15 20 25

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

normalized concentration

time (years)

Predicted Model and Measurements at 400 Meters

benzene 

toluene 

o−xylene

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.
26

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.
27

0 50 100 150 200 250 300 350 400 450 500

−3

−2

−1

0

1

2

3

distance (meters)

normalized prediction residuals

Normalized Prediction Errors at 19 Years

benzene 

toluene 

o−xylene

0 50 100 150 200 250 300 350 400 450 500

−3

−2

−1

0

1

2

3

distance (meters)

normalized prediction residuals

Normalized Prediction Errors at 11 Years

benzene 

toluene 

o−xylene

Figure 12. Prediction error normalized by the posterior covariance (σ+

y

) versus distance for strategy 60a.
28

0 5 10 15 20 25

−3

−2

−1

0

1

2

3

time (years)

normalized prediction residuals

Normalized Prediction Errors at 250 Meters

benzene 

toluene 

o−xylene

0 5 10 15 20 25

−3

−2

−1

0

1

2

3

time (years)

normalized prediction residuals

Normalized Prediction Errors at 200 Meters

benzene 

toluene 

o−xylene

Figure 13. Prediction error normalized by the posterior covariance (σ+

y

) versus time for strategy 60a.
29

0 50 100 150 200 250 300 350 400 450 500

−3

−2

−1

0

1

2

3

distance (meters)

normalized measurement residuals

Normalized Prediction Errors at 19 Years

benzene 

toluene 

o−xylene

0 50 100 150 200 250 300 350 400 450 500

−3

−2

−1

0

1

2

3

distance (meters)

normalized measurement residuals

Normalized Prediction Errors at 11 Years

benzene 

toluene 

o−xylene

Figure 14. Prediction error normalized by the posterior covariance (σ+

y

) versus distance for strategy 60b.
30

0 5 10 15 20 25

−3

−2

−1

0

1

2

3

time (years)

normalized measurement residuals

Normalized Prediction Errors at 400 Meters

benzene 

toluene 

o−xylene

0 5 10 15 20 25

−3

−2

−1

0

1

2

3

time (years)

normalized measurement residuals

Normalized Prediction Errors at 350 Meters

benzene 

toluene 

o−xylene

Figure 15. Prediction error normalized by the posterior covariance (σ+

y

) versus time for strategy 60b.
31

Table 4. Lag One Correlation Coefficients for Normalized Prediction Errors

11 years 19 years 200 (350) meters 250 (400) meters

Strategy 60a

Benzene 0.88 0.89 0.99 0.99

Toluene -0.51 -0.26 0.98 0.88

o-xylene -0.64 -0.50 0.98 0.90

Strategy 60b

Benzene 0.88 0.90 0.99 0.99

Toluene -0.51 -0.26 0.70 0.60

o-xylene -0.64 -0.50 0.76 0.70

Figure 16. Comparison of algorithm performance for varying measurement strategies.

0 50 100 150 200 250 300 350

0

0.2

0.4

0.6

0.8

1

cost of measurement strategy (thousands of dollars), JC

average of normalized posterior st. deviation, JA

Performance of Measurement Strategies

 5a

 5b

 5c

 5d

10a

10b

10c

10d

30a

30b

30c

30d

60a

60b 60c

60d

INFERIOR STRATEGIES

PARETO OPTIMAL STRATEGIES
32

The total cost of each sampling strategy depends on the number of sampling locations (Nzx) and the number of times a

sample is taken (Nzt):

C w zx sNzt J = a N + a (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 Marginal Marginal Marginal Marginal cost of uncertainty uncertainty uncertainty uncertainty uncertainty reduction reduction reduction reduction reduction / percent percent percent percent percent reduction reduction reduction reduction reduction

5d ó

10d $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.
33

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.
34

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 (θn) or a fraction of the total

void space (Sn).

n

M

S

V

M

n n

n

T

n n

n

ρ

ρ θ

=

=

(28)

At a larger scale (i.e., multiple pores, Figure 17b) we can repeat this calculation of θn 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

Section 5

Multi-Dimensional Problem Formulation

Figure 17. Definition of NAPL saturation as a macroscale property.

     a) Residual Blob  

b) Multiple Pores

c) Continuum

saturation defined at every point in continuum

     C                   



 



Ca,j

   Mn

δ 







 











 







 

VT

VT
35

The mass balance for this control volume is expressed mathematically (adapted from Abriola and Pinder 1985):

( )( )( ) ∑=

+ ∇ • − ∇ • ∇ = = ∂

∂ Ni

i

i i i i i i I I

t 1

ρκθ κηκ , ρκθ κ vκηκ , θ κDκ , ρκ ηκ , κ , β ,κ , (29)

where:

ρκ volume averaged density of phase κ in the control volume

θκ volume fraction of phase κ in the control volume

ηκ,i mass fraction of species i in phase κ in the control volume

vκ velocity of phase κ

Dκ,i dispersion coefficient of species i in phase κ

Iκ,i total inter-phase exchange of species i to phase κ in the control volume

Iβ,κ,i inter-phase exchange of species i from phase β to phase κ in the control volume

Figure 18. Mass balance of each species (i) in each phase (κ) in control volume.

rate coefficients, μ (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 (κ = 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 Nj

chemical constituents [j]). For each species (j) in each phase (κ), the total mass (Mκ,i) 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 (Fκ,i) and a dispersive flux (θκ,i) 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 Iβ,κ,i (i.e., exchange of mass from phase β to phase κ).

F F z, ,i z, ,i # # +

F F x, ,i x, ,i # # +

Fz, ,i

Fy, ,i F F y, ,i y, ,i # # +

Fx, ,i

x, ,i

z, ,i z, ,i # # +

x, ,i x, ,i # # + z, ,i

y, ,i y, ,i y, ,i # # +
36

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:

η

θ

κ

κ

κ

κ κ

κ

κ

κ

,

,

,

i

i

i

i=1

N

=1

N

i=1

N

i=1

N

i

i

∑

∑

∑

∑

=

=

=

=

1

1

0

I I

I

(30)

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+Nj

) 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 (ηv,i = ηk,v = 0). We make the further

simplifications that the soil (i = s) and water (i = a) species do not change phases (Ik,a = Ik,s = 0), and only exist in the solid

and aqueous phases, respectively (ηa,s = ηn,s = ηs,a = ηn,a = 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*Nj

) mass fractions and 3 volume fractions:

mass fractions

ηs,s soil in the solid phase

ηa,a water in the aqueous phase

ηs,j chemical j in the solid phase

ηa,j chemical j in the aqueous phase

ηn,j chemical j in the NAPL phase

for j = 1,2,...,Nj

 chemicals

volume fractions

θa aqueous phase

θs solid phase

θn 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 vs = 0 , Ds,i = 0. Thus the mass of soil in the control volume

doesn’t change over time. The soil species balance becomes:

( , ) = 0

∂

∂ s s s s t ρ θ η (31)
37

The water species is also assumed to be present only in a single phase. The species balance is:

( , ) + ∇ • ( , ) − ∇ • ( , ∇( , )) = 0

∂

∂ a a a a a a a a a a a a a a a t ρ θ η ρ θ v η ρ θ D η (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 , Dn,i = 0. The only movement of chemicals occurs in the aqueous phase. For each

chemical j the species balance is:

( ( )) 0

( ) ( )

, ,

, , , ,

− ∇ • ∇ =

+ + + ∇ • ∂

∂

a a a j a j

a a a j n n n j s s s j a a a a j

D

t

ρ θ η

ρ θ η ρ θ η ρ θ η ρ θ v η (33)

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:

( ) 0

1

+ , , = ∂

∂ ∑=

N j

j

s s s a j I

t ρ θ (34)

Again, an immobile medium is assumed. In the aqueous phase, water and chemical constituents are present:

( ) ( ) 0

1

, ,

1

+ ∇ • − , , − = ∂

∂ ∑ ∑ = =

j N j

j

n a j

N

j

a a a a a s a j I I

t ρ θ ρ θ v (35)

In the NAPL phase, only the chemical constituents are present. Note that this phase is immobile:

( ) 0

1

− , , = ∂

∂ ∑=

N j

j

n n n a j I

t ρ θ (36)

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 (Is,n,j = In,s,j = 0). We are left

chemical mass only transfers between the aqueous-solid (Ia,s,j and Is,a,j) phases and the aqueous-NAPL (In,a,j and Ia,n,j)

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 (ca,j). Assuming a dilute aqueous solution (i.e., the contaminant

contributes a negligible amount to the density of the aqueous phase), ca,j may be expressed as a function of the mass

fraction:

a

a j

a j

c

ρ

η ,

, = (37)

where ρa

 is the density of water.
38

The transfer of chemical mass (j) from the aqueous to the solid phase is then described as:

s a j

a j

a s j d j s e

j I

t

c

I K n

t

s

, ,

,

, , , (1 ) = − ∂

∂ = = − ∂∂

ρ (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., ρs ≈ ρs,s and θs ≈ θs,s). 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).

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 19

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:

j

N

j

j j

sat

j

eff

a j

j

c c m m γ















 = ∑=1

, (39)

Figure 19. Representation of non-equilibrium NAPL dissolution at the model scale.
39

This is the same form as presented in Section 3. This form of competitive dissolution directly relates NAPL saturation (Sj)

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 (σx

) associated with diffusion is σx = (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:

a j a n a j

eff

a n j a j I c c I , , , , , , = −μ( − )θ = − (40)

We consider μ 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:

[ ]b

x t n μ(x,t) = S ( , ) (41)

where b is the estimated parameter describing the power law dependence of μ 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:
40

a s

s s n I

t

M

t

M ,

, = ∂

∂ = ∂

∂ (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, Ia,s .

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:

t

M

t

a

a a a ∂

∂ = −∇ • = ∂

∂ (θ ) (θ v ) (43)

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:

( ) ( ) 0 ( )

, , ,

, , , + ∇ • − ∇ • ∇ = ∂

∂

+

∂

∂

+

∂

∂

a a a j a a j a j

n j s j a j a v c D c

t

c

t

M

t

M

θ θ

θ (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):

t

c

c c n S R

v c D c

t

c

a j

a j a e a j

eff

a j

a a a j a a j a j

a j a

∂

∂ = − − −

+ ∇ • − ∇ • ∇

∂

∂

,

, ,

, , ,

,

( ) ( 1)

( ) ( ) ( )

μ θ

θ θ

θ

(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:

( ) ( ) 1 ( ) ( )

, , , , ,

,

a j

eff

a j

j

a e a j a j

j a e

a j

j

a j a c c

R

S n D c

R S n

c

R

v

t

c = − ∇ + ∇ • ∇ + − ∂

∂ μ (46)

This equation is subject to the homogeneous boundary and initial conditions:

( ,0) 0

( , ) 0          on the domain boundary  

,

,

=

Ω = Ω

x a j

a j

c

c t

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 (Da,j), are all

functions of the water saturation (Sa), and the macro mass transfer coefficient (μ) 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):

a s

N

j

a j

d j s s e

s I

t

c

K n

t

M j

,

1

,

, , (1 ) = ∂

∂ = − ∂

∂ ∑=

ρ (47)
41

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 ρa to be treated as a constant. After applying the

third assumption:

∑ ∑ = =

+ ∇ • ∂

∂

+ <<

j N j

j

n a j a a a a a

N

j

s a j v

t

I I

1

, ,

1

, , ρ (θ ) ρ (θ ) (48)

the water phase balance becomes identical to equation (43), the ground-water flow equation. In terms of saturations, the

flow equation is:

( ) ( ) a e a

a e S n v

t

S n = −∇ • ∂

∂ (49)

The flow equation is completed by substituting in the generalized form of Darcy’s Law. This yields the non-linear

Richards equation.

[ ( ) ] ( ) Kk S h

t

S n r a

a e = ∇ • ∇

∂

∂ (50)

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 Sa 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 Sa

 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):

d

aSe f n aS r k ( ) ≈ [ ( )] (51)

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:

∇( ) Kk (S )∇h = 0 r a (52)

The boundary conditions applied to (52) are:

( , ) 0

( , )

2

1

• Ω =

Ω =

n q t

h t hc

G G

In other words, these are simple no-flow (Ω2) or constant head (Ω1) 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:

( ) ( ) 0

1

+ , − , = ∂

∂ ∑=

N j

j

a j a

eff

n Sn ca j c S

t ρ μ (53)

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:

e n n j

e j j n

j

j

j

n S mw

n S mw

m

m

ρ

ρ =

∑ (54)
42

After some simplification, the NAPL dissolution equation becomes:

( ) ( ) 0

1

,

1

, =















 −













 +

∂

∂ ∑ ∑ = =

j N j

j

a j

N

j

sat

a j

j

j j j

n n

n

n n a c c

mw

S

S

mw S S

t

ρ γ

ρ

ρ μ (55)

This equation requires an initial condition of the fractional NAPL saturation:

S (x,t 0) S (x) j = = jo

and is completed with the volume averaged quantities:

∑

∑

∑

=

=

=

=

=

=

j

j

j

N

j

n j

n

N

j

j j

n

n

N

j

j j

n

S S

S

S mw

mw

S

S

1

1

1

ρ

ρ

(56)

where ρj

 is the density of pure phase chemical j and mwj 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 Sj

 rather than Sn:

0 ,

, =

















+ − ∂

∂

j

sat a j

a j

j

j

n n

n

a

j c

c

mw

S

S

mw S

t

S

ρ ρ

μ (57)

Numerical Solution of State Equations

Along with boundary and initial conditions, equations (46), (52), and (57) are the final state equations for our multidimensional 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 (Sj

, Sa, ca,j) is discretized at grid nodes and represented with the

trial basis functions. Other parameters (D, v, Rj, μ) are considered constant over elements. The spatially discretized (but

continuous in time) state equations are presented in Table 6.
43

The finite element discretization produces a set of ordinary differential equations in the discretized states Sj(t), Sa(t), and

ca,j(t) with the coefficient matrices A, B, F, G, H, K, and L. All of the coefficient matrices are dependent on the

discretization (i.e., the basis functions and finite element grid). In addition, all of the coefficient matrices except A and B

are dependent on the states. Table 6 includes the discrete steady-state ground-water flow equation that solves for the

flow field as a function of water saturation. Since we have not made the hydraulic head a state in our formulation, the

ground-water flow equation does not produce a solution for a state, but is nonetheless necessary for calculating the

velocity field for the other state equations.

As mentioned above, we use a fully-implicit finite difference discretization in time for the state equations. We use the

fully-implicit formulation to maximize stability of the solution. Due to the non-linearities in the equations (i.e., the

parameters depend on the states), each equation must be solved iteratively. On each iteration, the set of discrete

equations for each state is solved with a bi-conjugant gradient algorithm that uses a partial Cholesky decomposition preconditioner. Other, and probably more efficient, solvers could potentially be used for this problem, however we are

limited to solvers that can address non-symmetric matrices. We chose the bi-conjugant gradient method primarily for

ease of implementation. For the iterative part of the algorithm, we utilize a Picard method, both for ease of execution and

efficiency. The Picard method, as opposed to gradient-based search algorithms (e.g., Newton-Raphson), does not

require evaluation of derivatives. The derivatives for this problem would have to be calculated numerically and could be

quite costly. The Picard method likely requires more iterations than a gradient-based method, but for this problem could

prove to be more stable.

Table 6. State Equations for Multi-dimensional NAPL Source Problem

Final State Equations

( ) ( ) 1 ( ) ( )

, , , , ,

,

a j

eff

a j

j

a e a j a j

j a e

a j

j

a j a c c

R

S n D c

R S n

c

R

v

t

c = − ∇ + ∇ • ∇ + − ∂

∂ μ      for each j

c , (x,0) = 0 a j















 = − − ∂

∂

j

sat a j

a j

j

j

n n

n

a

j c

c

mw

S

S

mw S

t

S

ρ ρ

μ ,

,      for each j

S (x,0) S (x) j = jo

∇( ) Kk (S )∇h = 0 r a

S (x,0) S (x) a = ao

Spatially Discretized State Equations

= B [ ] F + H − ( ) ( )

,

, 1 t

dt

d t

a j

a j ❃

❃      for each j

[ ] ( ) ( ) ( )

,

1 t t

dt

d t

j a j

j ✳ ❃

✳ = A G + K −      for each j

L❈= P

with    H(Sj, Sa, b), F(Sj, Sa, b), G(Sa, Sj), K(Sj, b), L(Sa)

A and B constants
44

In summary, the steps of the solution to the state equations are:

1. Initialize the model with a velocity field, initial Sj

, initial Sa, and initial ca,j=0.

2. Calculate the effective solubility, mwn

, ρn, Sn, and μ 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 z p j = ca p j (x p ,t p ) +ω p j          p =1, 2,...N (58)

where p is the measurement index. Again the measurement error ωp,j is considered a normally distributed random

variable with zero mean and variance Cω.

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:

( , τ) (τ) τ

1 ( )

( , τ) (τ) τ

1 ( )

] à[

1 ] à[

1

] à[ 1 ] à[

] à[ 1 ] à[

] , [ à 1 ] ( )

1

, [ à

1

dtd n t

n

t C T

n

dtd c t

c

t C T

c

b

T

j N

j

ao jo

jo S

T

jo jo

ao ao

ao S

T

ao ao

p

a j c pj z p C T Nz

p

p

a j c pj z

j N

j

J

υ υ

υ υ





 









 





− +

− +

− − + −

∑

=

′ − ′ −

′ + ′ − ′

′ − ′ −

′ + ′ − ′

− − ∑

=

∑ −

=

=

∫ ∫

∫ ∫

υ

υ

ω

b b C b b

S S C S S

S S C S S

(59)
45

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:

ref

a

a

a

ref

j

j

j

S

S S

S

S

S

ln

ln

′ =

′ =

(60)

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:

(61)

∫ ∑ −

=

+ ′ + ′

∫ + + +

∂

∂

∑

=

+

∫

′ +

′ + ′ +

∂

∂ ′

∑

=

+

− +

− +

− − + −

∑

=

′ − ′ −

′ + ′ − ′

′ − ′ −

′ + ′ − ′

− − ∑

=

∑ −

=

=













































































 









 





∫ ∫

∫ ∫

dt

j N

j

t j t a t T

a

t dt c t t a j t

t

t a j t t T

c j

j N

j

dt

t n t a j t a t

t a t t j t

t j

t T

n j

j N

j

dtd n t

n

t C T

n

dtd c t

c

t C T

c

b

T

j N

j

ao jo

jo S

T

jo jo

ao ao

ao S

T

ao ao

p

a j c pj z p C T N z

p

p

a j c pj z

j N

j

J

1

1

( )] à ( )] exp[ à ( ) exp[ à 2

( ) à ( ) à ( ) , ( ) à ( ) , à

( ) à ( ) , à

1

2

( ) à ( ) , à ( ) exp[ à ( )]à

( )] à ( ) exp[ à ( )] à exp[

( ) à

( ) , à

1

2

( , τ) (τ) τ

1 ( )

( , τ) (τ) τ

1 ( )

] à[ 1 ] à[

1

] à[ 1 ] à[

] à[ 1 ] à[

] , [ à 1 ] ( )

1

, [ à

1

S S

F c H

c

B

K S c

S G S

S

A

b b C b b

S S C S S

S S C S S

λ

λ

λ

υ

υ

ω

υ

υ

υ υ

υ υ
46

The hats (^) 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 (δ) 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 (δJ) and set each partial variation equal to zero. The first variation of J is represented symbolically as:

a

a

c j

c j

n j

n j

c

c

n

n

N

j t T j

t T j

a t T

a t T

N

j

jo

jo

ao

ao

N

j

j

j

a

a

N

j

a j

a j

J J J J J b

b

J S

S

J

S

S

J S

S

J S

S

J S

S

J S

S

J

c

c

J J

j

j j j

δλ

λ

δλ

λ

δλ

λ

δυ

υ

δυ

υ

δ δ

δ δ δ δ δ δ δ

∂

∂

+

∂

∂

+

∂

∂

+

∂

∂

+

∂

∂

+

∂

∂ ′ +

∂ ′

∂

+

′ ∂ ′

∂ ′ +

∂ ′

∂ ′ +

∂ ′

∂ ′ +

∂ ′

∂ ′ +

∂ ′

∂

+

∂

∂ =

∑

∑ ∑ ∑

= =

=

=

= = = =

,

,

,

1 ,

1 1 1

,

, (62)

The set of coupled non-linear equations resulting from setting δJ = 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 EulerLagrange 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).
47

Parameter Update Equations

































∂

∂

+

∂

∂

+

∂

∂ = + ′ − ′ ∫ ∫ ∑ ∑

T T

j

c j

T

c j

T

T

a j

j

n j

T

T

a j

T

b a j dt

b

H

b

F dt c

b

K b b C S S c

0 0

, , , , , à à à à

à à à

)] à à à [exp( à λ λ λ

































∂ ′

∂

+ ′ − ′ 















∂ ′

∂

+

∂ ′

∂ ′ = ′ − ′ ∫ ∑ ∫ ∑

T

j

n j

ao

T

T

a j

T

c j

ao

T

j

c j

ao

T

T

ao ao S a j dt

S

G dt S S

S

H

S

F S S C c ao

0

,

0

, , , à à

)] à à [exp( à à à à à à λ λ λ

∑

∫ ∑ ∫ ∑

− =

















































∂ ′

∂

+

∂ ′

∂

+ ′ − ′ 















∂ ′

∂

∂ ′

∂ ′ = ′ −

′

′

j

n j

T

S

T

j

n j

jo

T

jo

T

a j

T

a j

T

j

c j

jo

T

T

a j

jo

T

jo jo S

C A t

dt

S

G

S

K dt S S c

S

F

c

S

H S S C

jo

jo

( 0) à

à à à

)] à à à [exp( à à

à

à à

,

0

, ,

0

, ,

λ

λ λ

Model Error Equations

= − ∑

j

n n j n C , à υà υ λ

= − ∑

j

c c j c C , à υà υ λ

Adjoint (Backward) Equations











 = − − ′ − ′ − ∂

∂

∑− −

a j n j c j

m

m

c j a j

c j B C z c S S K F

t , , , ,

, 1 1 à à à à ] à à [ à ] exp[ à

à

λ λ

λ

ω      for each j

, = 0 t=T

λc j















 ′ − ′ + ′ 













 − ∂ ′

∂ − +

∂ ′

∂ = ∂

∂ − ] à exp[ à ] à à exp[ à à à à à

à

,

, 1

j

T

a j a

j j

n j

n j K S S S

S

K G

S

G A

t

λ λ

λ     for each j

, = 0 t=T

λn j

∑

∑

















+ +

∂ ′

∂

+

∂ ′

∂ − ′ ′ − ′

































∂ ′

∂

+

∂ ′

∂ = − ′

j

n j

T T

a j

T T

a j

T

a

T

a

T

a j

T

a

j

c j

T

a

T

a

T

a j

T

a a

c G c K

S

K

S

G S S S

S

H

S

F S c

, , ,

, ,

à à à à à

à à

] à à ] exp[ à exp[

à à à

exp[ à ] à à

λ

λ λ

= 0 a t=T

λ

Table 7. Set of Euler-Lagrange Equations for Multi-dimensional NAPL Source Problem
48

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 (Nz

 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.
49

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

Section 6

Conclusions
50

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 twodimensional 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.
51

Appendix A: Table of Variables ariables ariables ariables and Operators Operators Operators Operators

Variables Units

as cost per chemical sample analysis [$]

aw cost per well installation [$]

cκ,i concentration of i in phase κ [M/L3

]

cj

eff effective solubility for j in water [M/L3

]

cj

sat pure phase saturated solubility for j in water [M/L3

]

Cy state covariance function matrix -variesCω measurement error covariance matrix -variesCυ model error covariance function matrix -variesCα prior parameters covariance matrix -variesDκ,i dispersion coefficient of i in phase κ [L2

/T]

f

om fraction organic matter [M/M]

F advective flux [M/T/L3

]

i index for species (n=chemical, a=water, v=air, s=soil) -noneI

β,κ,i interphase exchange of i from phase β to phase κ [M/T/L3

]

I

κ,i total interphase exchange of i into phase κ [M/T/L3

]

j index for NAPL chemical -noneJ performance index -noneJA performance measure for accuracy -noneJC performance measure for cost [$]

k estimation iteration index -nonekrw relative permeability -noneK hydraulic conductivity [L/T]

Kd,j sorption coefficient for chemical j [(M/L3

)/(M/M)]

Kom organic matter partitioning coefficient [(M/L3

)/ (M/L3)]

L end-point of one-dimensional problem domain [L]

mj moles of j in NAPL phase [moles]

Mκ,j mass in VT of j in phase κ [M]

n measurement index -nonene effective porosity [L3

/L3]

Nα number of parameters -noneNκ number of phases -noneNj number of chemical constituents -noneNi number of species (air, water, soil, chemicals) -noneNe number of equations -noneNt number of time steps -noneNx number of spatial steps -noneNz number of measurements (Nz = Nzz + Nzt) -noneNzx number of measurement locations -noneNzt number of measurement times -noneQ flux of water through the one-dimensional source [L3

/T]

R retardation factor -nones sorbed concentration (mg chemical/mg organic matter) [M/M]

Sa water saturation [L3

/L3]

Sn,j NAPL fractional saturation of chemical j [L3

/L3]

Sr residual NAPL saturation [L3/L3]
52

Ss specific storage [1/L]

t time dimension (days) [T]

u model forcing function -variesvκ velocity of phase κ [L/T]

Vn volume of NAPL in control volume [L3

/L3]

Vp volume of pore spaces in control volume [L3

/L3]

Vv volume of voids in control volume [L3

/L3]

VT total volume [L3

]

x spatial dimension [L]

xs location of point source in one-dimensional problem [L]

y model state vector -variesz measurement vector -variesα model parameter vector -variesβ phase index (NAPL=n, vapor=v, aqueous=a, solid=s) -noneχ mole fraction of constituent in NAPL mixture [mol/mol]

γ activity coefficient of constituent in NAPL mixture -noneθ dispersive flux [M/T/L3

]

κ phase index (NAPL=n, vapor=v, aqueous=a, solid=s) -noneηκ,i mass fraction of i in phase κ [Μ/Μ]

μ macro mass transfer rate coefficient [1/T]

θκ volume fraction of phase κ [L3

/L3]

ρκ volume averaged density of phase κ [Μ/ L3]

ρ1 lag one correlation coefficient -variesσ standard deviation -variesυ process/model error vector -variesω measurement error vector -variesOperators

A state equation operator

Ac state equation for chemical mass in the dissolved phase

As state equation for chemical mass in the sorbed phase

An state equation for chemical mass in the NAPL phase

b RHS operator of forcing terms from solved state equation

B boundary conditions for state equation

D LHS operator from implicit form of solved state equation

F forward operator

G model error operator

I initial conditions for state equation

M 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.
53

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 Water 30(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 twodimensional 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.

References
54

James, A.I., Graham, W.D., Hatfield, K., Rao, P.S.C., and Annable, M.D. 1997. “Optimal estimation of residual nonaqueous 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):1201-1211.

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. 1986. “A comparison of several methods for the solution of the inverse problem in two-dimensional steadystate groundwater flow modeling.” Water Resour. Res. 22(5):1543-1570.

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):1551-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.

32(5):1131-1161.

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):1006-1021.

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-1195.

Powers, S., Abriola, L., Dunkin, J., and Weber, W. 1994. “Phenomenological models for transient NAPL-water masstransfer 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(10):2691-2705
55

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.

Res. 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.