GROUND WATER SOLUTE TRANSPORT MODELING
CASE STUDY BENCHMARKS
AND
VALIDATION METHODOLOGY
by
David Scott Long, PE
B.S.C.E., University of Utah at Salt Lake City, 1977
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Civil Engineering
1993
1993 by David Scott Long
Ail rights reserved.
MATHCAD is a registered trademark of Mathsoft, Incorporated
MSDOS is a registered trademark of Microsoft Corporation
This thesis for Master of Science
degree by
David Scott Long
has been approved for the
Department of
Civil Engineering
by
Edward B. Nuhfer
Date
Long, David Scott (M.S., Civil Engineering)
Ground Water Solute Transport Modeling pase Study Benchmarks
and Validation Methodology
Thesis directed by Professor Jonathan T. H. Wu
ABSTRACT
The primary objective of this thesis is to perform benchmark
modeling of ground water solute transport from a rectangular continuous
source. Benchmarking is performed in two case studies by comparing
MOC and SUTRA numerical solution computer programs, and TACT, a
coded analytical solution. Also presented, is the general theory and
practice methodology for ground water contaminant transport code
selection, validation, and use. The first case study is a typical generic
site, and the second case study focuses on an actual site, using data
derived from literature and field data, at Liberty Aircraft Plant, New York.
Both case studies model the transient contaminant transport in an areal
homogenous aquifer with a steadystate flow field. There is good
general agreement between models in the second case, but not in the
first case with a gridinterior source. Each of the two numerical codes
studied have a tendency for specific shape bias, which is evident in all
subcases. Subcases vary run options, grid size, dispersivity, source
volume, and source type. In numerical scenarios, SUTRA is typically
most accurate near sources/sinks, continuing out to the advective center
region, while MOC is most accurate at the leading edge and outer
contours. However, the analytical solution (TACT) is the standard for
uniform flow and aquifer properties, except in the proximity of the source.
The numerical codes' spatialconcentration volumes are comparable to
IV
each other and to the analytical solution in the second case. Source
volume input option types are found to influence the results in a
magnitude similar to possible aquifer parameter variation, due to
apparent local divergent flow mass balance difficulties around the grid
interior source cell control area. While the analytical solution for a
continuous rectangular source is exact in the farfield, theoretical
difficulties with the shortduration or nearfield spatialconcentration
volume are evident due to the nonFickian dispersion/diffusion shape
distribution in a divergent Lagrangian velocity field at the source.
This abstract accurately represents the content of the candidate's thesis.
I recommend its publication.
Signed
v
ACKNOWLEDGMENTS
Professional Technical Reviewers:,
Brian J. Long
Duane E. Long
Roger S. Stillwater
Manuscript Reviewers:
Hannah L. Kelminson
Marlene V. Long
vi
CONTENTS
Chapter EflOfi
1. Introduction...................................................1
1.1 The Role of Ground Water Modeling...............................2
1.2 Problem Statement The Ground Water Model
Credibility Question........................................4
1.3 Objectives....................................................11
1.4 Approach......................................................11
1.4.1 Approach for Solute Transport Code Selection.................12
2. Ground Water Solute Contaminant Transport Modeling
Background and State of the Art.............................13
2.1 Ground Water Flow Governing Equations and Model Type
Selection................................................... 15
2.2 Solute Contaminant Transport Governing Equations..............19
2.3 Review of the Site Modeling Process...........................22
2.4 Establishing Ground Water Model Quality Assurance
and Control.................................................24
2.5 Solute Contaminant Transport Model Evaluations................25
2.6 Industry Model Usage and Needs................................27
3.0 Solute Transport Analytical Solution..........................32
4.0 Solute Transport Numerical Solution Modeling Techniques.......39
4.1 SUTRA Model Description.......................................40
4.2 MOC Model Description.........................................46
4.3 Practical Considerations of Numerical Modeling................51
4.3.1 Model Grid Scale and Grid Types..............................51
4.3.2 Boundary Conditions, Sources, and Sinks......................52
4.3.3 Transient Flow Numerical Modeling............................56
4.3.4 Aquifer Parameter Estimation.................................57
5.0 Benchmarking Case Studies.....................................59
VII
5.1 Case Study/Benchmark 1.....................................60
5.1.1 Problem Definition for Case Study/Benchmark 1.............60
5.1.2 Case Study/Benchmark 1 Analytical Model Results and
Discussion.................'.............................66
5.1.3 Case Study/Benchmark 1 MOC Model Results and
Discussion...............................................66
5.1.4 Case Study/Benchmark 1 SUTRA Model Results and
Discussion...............................................68
5.1.5 Case Study/Benchmark 1 Comparison and Conclusions.........78
5.2 Case Study/Benchmark 2......................................85
5.2.1 Problem Definition for Case Study/Benchmark 2.............91
5.2.2 Case Study/Benchmark 2 Results, Comparison, and
Conclusions..............................................96
6. Summary and Conclusions......................................119
6.1 Summary....................................................119
6.2 Conclusions................................................120
6.3 Conclusion and Interpretation Summary......................122
References......................................................128
Appendix A. Ground Water Modeling Term Definitions..............135
Appendix B. Sample Output Data for Case 2A......................139
viii
1. Introduction
Ground water solute contaminant transport modeling plays an
important role in investigations and predictionsincluding site
characterizations, riskbased exposure assessments, compliance
reporting, remedial engineering, planning tradeoffs, and management
evaluations. Ground water modeling provides a means for organizing
data and conceptualizing the hydrogeologic system dynamics, and for
using that information to make predictions and to evaluate hypotheses.
Successful modeling requires selecting a "validated" math model,
applying a correct model type, interpreting results, verifying the
application, and calibrating results with field data from a site
characterization study. A multidisciplinary effort may be required to
properly apply all of the technical principles to a conceptual model.
Model validation composed of verification, calibration, quality control,
and documentation is required to assure that all model analysis is
plausible and defensible in a court of law. This is true of the conceptual
model, the math model computer program, and the applied site model.
Model validation requires historical confidence of a validation/verification
methodology to select the appropriate math model code, and it requires
detailed judgment in the site calibration. (Essential terms are defined in
Appendix A.)
A general lack of confidence persists concerning the accuracy of
ground water modeling results because of appreciable differences in
some comparative model application studies as well as some difficulty in
actual field representations. Model theoretical/conceptual problems and
computer program coding errors continue to be reported. As a result,
there exists a need for verification, benchmarking, and validation of the
modeling computer code.
1
1.1 The Role of Ground Water Modeling
Ground water is an important resource supplying half of this
nation's drinking water (Wentz, 1989). Once contaminated, it is difficult to
restore. Manmade chemicals have been detected in 30% of the nation's
water supplies. Modeling can be useful in site characterization,
exposure assessment, and remediation processes (EPA, 1990,1992), by
predicting rate, direction, and concentration distribution. Contaminant
transport modeling use thus far has met varying degrees of success. At
least general information on important processes and trends can usually
be obtained.
Landfills or open dumps were the final destination of most solid
waste and treated hazardous waste before the Resource Conservation
and Recovery Act (RCRA). The EPA RCRA regulations, Title 40 Code of
Federal Regulations Section 264 Subparts FH and N, establish
permitting and the minimum technical requirements for currently
operating landfills (CFR, 1990; EPA, 1990). This section also sets forth
the operational and monitoring requirements for leachate detection.
Ground water modeling may support remedial actions and designs at
older and existing RCRA landfills. These are typically older landfills that
have been closed or have had lateral extension portions closed at the
same general site acreage. The RCRA corrective action process
requires that a RCRA permitted facility assessment (RFA) be performed to
gather information on the site units and releases (EPA, 1991,1985). A
decision is made on whether to proceed with further investigation. Next,
a RCRA facility investigation (RFI) is made from site inspections and
sampling. Fate, pathway, and transport assessments are then made.
(This is a policy similar to CERCLA's RI/FS process.) A thorough
hydrologic examination of the site is required to support the feasibility of
various remedies and the selected proposed remedy for a site.
If the landfill site was closed before the Resource Conservation
and Recovery Act (RCRA), it would not have RCRA landfill design,
2
/
operation, and monitoring requirements. A site predating RCRA would
then legally fall under the Comprehensive Environmental Recovery,
Compensation and Liability Act (CERCLA). However, under CERCLA no
action is taken unless the site is bad enough to be placed on the National
Priority List (NPL). CERCLA provides the EPA with authority to evaluate
and rate sites and determine the need for a response. The enactment
plans are provided in the National Contingency Plan (NCP) and research
and documentation is provided by the EPA Office of Research and
Development (ORD) (EPA, 1985). The NCP Section F (as required
under Section 104 of CERCLA) establishes requirements for a Remedial
Investigation and Feasibility Study (RI/FS) to be carried out at all
federally funded cleanup sites. In this evaluation, sufficient information
must be provided to allow management and engineering decision
makers to select potentially applicable technologies. Procedures for
evaluating releases, remedies, and use of funds are found in the
CERCLA Section 105 Hazardous Substance Response Trust Fund
(Superfund). The NPL has about 1200 sites, with 300 of these as former
Municipal Solid Waste Landfills (MSWLS), that received hazardous
waste before RCRA restrictions.
A risk assessment analysis of EPA survey data for 19,000
hazardous waste sites (derived from EPA, 1984a; 1989a) which
assumed average preRCRA landfill conditions, concluded that the risk
from old abandoned sites is substantial at approximately 17 x 10+3
persons per lifetime at risk nationally (Long, 1992). This assessment
reveals that more attention should be applied to those older sites that
have not made the NPL. The risk from new regulated landfills (EPA,
1990,1986) are very low and are acceptable to protect human health
and the environment. The calculated risk is 1.25 x 10'6 persons per
lifetime at risk nationally (average basis, no confidence interval) for new
facilities permitted and designed under RCRA since November 8,1985;
and this value is considerably less than for abandoned sites.
3
1.2 Problem Statement The Ground Water Model
Credibility Question
Full confidence in ground water modeling results is apparently
lacking (Van der Heijde, 1990). Application problems and differences
have been reported in some comparative studies (Hamilton et al., 1985;
McLauglin, 1987; Freyberg, 1988; Grondin, 1990; Reddi, 1990; Leppert
et al.,,1993). Problems in model conceptualization and coding have also
been found (Goode & Konikow, 1989; ElKadi, 1988; Zomorodi, 1990).
Field testing is an important part of establishing ground water model
reliability. Yet, less than a third of the saturated flow model codes have
been tested and less than 15% of the saturated transport model codes
have been field tested (EPA, 1991). Some of the credibility question can
be attributed to improper conceptual assumptions, misapplied models,
underdocumented model codes, or misunderstanding of anticipated
accuracy from site models with scant data; but are there other problems?
The general scope of the theoretical problem may be better
understood by comparison to the historical acceptance of finite element
models in the structural engineering field. This is outlined in Table 11.
There are similarities in the basic formulated equations. The largest
difference appears to be uncertainty in the hydraulic conductivity which
can typically vary by two orders of magnitude compared to a modulus of
elasticity which will typically vary by less than 10%. Another aspect of
property variability is the prodigious spatial heterogeneity of many
hydrologic and geologic parameters compared to the relatively
homogeneous nature of most other materials. Enormous actual
variations are possible between conceptualized boundary conditions
and initial conditions that adversely affect application of the governing
equations. Site model conceptualization uncertainties can lead to errors
as great as those associated with unknown aquifer characteristics
(Leppert et al., 1993).
4
Table 11 Ground Water Models Compared to
Structural Models Large potential variations in ground
water parameters create model result uncertainty_____________
Ground Water Models Structural Models
Credibility unproved Credibility accepted
Q = KA i q = Ki orv = Ki/n P = K U stress = Ex strain
K Factor variation: 1
Heterogeneous material  Geologic variation  Hydrologic variation  Temp. & Density variation * Homogeneous material
Boundary conditions Many possibilities Boundary conditions Few possibilities
Q = Volumetric flow rate K = Hydraulic conductivity A = Crosssectional area i = Gradient q = Specific discharge velocity v = Velocity L = Load K = Stiffness E = Modulus of Elasticity U = Displacement
5
Another part of the credibility problem is the lack of ground water
model code maturity. Programming code errors were found recently in
some of the major numerical codes that have been used for years. This
problem is not unique to the ground water modeling field. This should
not be surprising because code problems were found in the major
structural finite element programs years after release, although mature
code errors typically have a smaller impact on most results. Errors were
recently found in the analytical equations for solute transport in a few
published sources that might be used to check the numerical codes.
Some errors are to be expected ("To err is human"), so the importance of
independent validation/verification cannot be overstated. A diligent
engineer must therefore do at least sanity checks for the problem at
hand. Source code validation and site model verification and calibration
processes should be done for proper modeling. In general, verification,
benchmarking, and validation of both numerical and analytical ground
water transport codes is still in its infancy.
There is apparently still some debate over the validity of defining a
heterogeneous system in a deterministic manner. Lack of confidence in
some areas of solute transport theory also persists (Anderson, 1984;
Beljin, 1988). It is recognized in the field that more work is required in
these problem areas:
Dispersion, scaling, and application of Fick's Law
Laboratory results in comparison to field values
Some chemical & biological processes and rates
Immiscible transport
Fractured media transport
Unsaturated media transport
Stochastic analysis
Phased transport
Fixation and ion exchange with transport
Verification, validation, and benchmarking of ground water
contaminant model codes
6
Complete and computer usable programs for general analytical
solutions require further work. An example is an analytical model code
for solute transport with nonlinear equilibrium sorption and ion exchange
(Konikow & Bredehoeft, 1978). Some purely analytical codes are
available for ground water transport in the literature, but they may need to
be set up and verified as public domain codes, easily obtainable for
convenient general professional use. However, proprietary and
commercial analytic codes can be purchased.
Discrepancies between theory and laboratory observations have
been reported (EPA, 1989b). There is also some debate over
representing dispersion in the form of Fick's Second Law of diffusion in
the nearfield spatial and temporal domain. Another reported
phenomena is that dispersivity increases with distance near the source,
as part of the "scale effect". For large times and distances, dispersivity
values have been used in successful calibration contour fits (Anderson,
1984). Large dispersivities in the field are the result of aquifer
heterogeneity creating macroscopic dispersion.
Although most of the public domain codes are available from the
United States Geological Survey (USGS) and EPA, these are not always
in usable form as received. Few of the government codes investigated
ran as delivered. They usually require modification, debugging, and
compiling. Sometimes this is due to differences in FORTRAN compilers
and computer operating systems. Usually job control setup with
input/output file designations must be added. (Compilation of one
prominent code yielded hundreds of errors and warnings.) For this
reason, it is advisable to acquire a compiled version through a software
agency, whenever feasible.
Documentation of the mathematical theory for the primary public
domain codes is usually good. However, the documentation for the user
input and run processes is generally poor. Errors have been found in
many of the input requirement descriptions of the codes used in this
thesis; and this casts an overall credibility question on the quality of the
7
code. Trial and error runs are usually required to determine the proper
input format. Test runs and completely documented benchmarking in
existing programs is frequently terse or nonexistent.
General modeling problems collected in the Geraghty & Miller
Survey of modelers (Rumbaugh & Ruskauff, 1993) are listed as follows in
order of the number of respondents:
Lack of good quality data
Modeling is of poor quality or too complex
Lack of trained/experienced/knowledgeable modeling staff
Models are difficult to calibrate or verify
Modeling results are often trusted blindly
Modeling is too expensive and too timeconsuming
Boundary conditions are difficult to understand and apply
Inappropriate code applied to problem
Models contain too many assumptions
Nonmodelers often misunderstand the results
Models exhibit a high degree of uncertainty that is difficult to
quantify
Model results are difficult to communicate to others
Problems are often poorly conceptualized
Heterogeneous systems are difficult to model
Models are often oversimplified
Modeling results are often misinterpreted
Misuse of codes
Models are difficult to validate
Lack of standardization for terms, codes, or procedures
Model limitations are often poorly communicated
Enduser client expectations are often too high
Modeling objectives are poorly defined
Nonuniqueness of the calibrated model
Lack of documentation of the model application
Lack of proper peer review for QA/QC
The industry now is interested in independent reviews of a
model's fundamental theoretical basis, tests, and intercomparison with
other modeling programs. Usually an orderofmagnitude result is
sufficient for professional judgments regarding concentration or travel
time, especially with data and degradation process uncertainty (ElKadi,
8
1988); but greater accuracy may be necessary for legal compliance and
decisions. Gone are the days when validation of a model code by a
developer is completely accepted without independent proof of model
code validation, especially in expensive remediation decisions and legal
defenses. See Figure 1 1.
The main types of model testing involved in validating of model
codes are:
Verifying the governing equation's theory, application,
algorithms and program coding
Benchmarking specifically defined problems by comparison
with: ,
 Other numerical solution models
 An analytical solution
 Field data case studies and history matching
Sensitivity studies and test cases
9
Figure 11 The Computer Code Model Validation
Process
D. Long S3
10
1.3 Objectives
The primary objective of this thesis is to perform benchmark
modeling of ground water solute transport from a rectangular continuous
source in two case studies. Benchmarking is performed by comparing
MOC and SUTRA numerical solution computer programs, and TACT, a
coded analytical solution. The secondary purpose, which supports the
primary purpose, is to research the general theory and practice
methodology for ground water contaminant transport code selection,
validation, and use.
The expected contributions of this study are to:
Perform benchmarking of two selected premier numerical program
packages, comparing each with the other and with a coded
analytical solution on two different case sites, modeling transient
solute transport in a homogenous quasi threedimensional areal
steadystate flow field with a continuous rectangular source.
Provide a general overview of current ground water modeling
solute transport governing equations and code attributes for code
model selection and modeling practice; and general principles of
model validation methodology.
1.4 Approach
The general theory and practice methodology for ground water
contaminant transport code selection, use, and validation is investigated.
Research was conducted into the theoretical equations relative to flow
and transport model types and stateoftheart practice to lay groundwork
for modeling in this thesis. Two benchmark case studies are modeled
with premier public domain numerical computer program packages,
SUTRA & MOC, selected from the criteria in Section 1.4.1. These codes
11
are compared to each other and to an analytical solution. Both cases are
represented as a homogeneous steadystate flow field with a continuous
rectangular source. The effect of changing different sensitive parameters
and run options on a gridinterior source are presented in Case 1. Case
2 is an actual case history, with data from the literature, compared by
parameter calibration. Subcases vary run options, grid size, dispersivity,
source volume, and source type. Comparison plots are made of the
spatialconcentration volumes, the outer contours, the centerline
concentration profiles, and sections through the concentration volumes.
Result differences are interpreted relative to the governing equations.
1.4.1 Approach for Solute Transport Code Selection
The solute contaminant transport model code selection for this
study was made based on the following criteria:
1. The computer model code must be available as public
domain.
2. The computer code is available for use on a MSDOS based
486DX.
3. The executable computer model code must be of a nominal
cost (not over $300).
4. The code has gained acceptance in the profession and in the
scientific community.
5. Some general code validation, verification, testing, or
benchmarking has been published in the literature.
12
2. Ground Water Solute Contaminant Transport Modeling
Background and State of the Art
There are several types of modeling methods. The theoretical
basis for solving flow and transport models may be summarized as
follows:
I. Differential Equation Numerical modeling
A. Finite Element (FE)
B. Finite Difference (FD)
C. Integrated Finite Difference (IFD)
D. Boundary (Integral) Element
E. Analytical Element (limited application)
II. Differential Equation Analytical Exact Solution modeling
III. Differential Equation Stochastic and Probabilistic modeling
IV. Compartmental Mass Phase Balance modeling
(limited application)
V. Method of Characteristics (MOC) and moving particle tracking
modeling (transport only)
VI. LaplaceTransformGalerkin (LTG) modeling
(linear transport only)
The differential equations of the first three classes are composed
primarily of the Darcy, water mass balance, and momentum equations.
Finite element and finite difference numerical solutions offer the greatest
versatility in solving complex problems. Finite Difference (FD) and
Integrated Finite Difference (IFD) methods solve partial differential
equations by the standard method of finite differences. Finite Element
methods solve the differential equations by integration and then forming
first order differential equations by weighted residuals or variational
means. The analytical element is a new, limited application method that
requires flow field input. Analytical solutions are equations with
13
deterministic, closedform solutions of the governing differential
equations.
Probabilistic and stochastic methods are considered
indeterministic. These include socalled geostatistical methods such as
Kriging. Kriging smoothes out random noise and estimates
heterogeneous spatially distributed parameters. Stochastic partial
differential equations are described by one or more dependent or
independent random variables (De Marsily, 1986). Stochastic methods
use a random function. One of the best stochastic methods is a spectral
analysis by way of Fourier transforms. The Monte Carlo method is
popular and statistically evaluates an ensemble of solutions from varying
parameters at each location, but can be computationally intense. The
probabilistic and stochastic methods can be used in combination with the
preceding methods (types I and II). Compartmental models may use
either differential or algebraic equations formulated for a specific
subsystem or subprocess (Calabrese & Kostecki, 1993). These are
generally oriented toward the chemical reaction systems and are weak in
the macroscopic physical formulation. An example of a compartmental
model is the use of a mass phase equilibrium to find soil matrix
concentration.
Particle tracking and MOC are popular transport methods that
diminish numerical dispersion. There is also a new approach, the
LaplaceTransformGalerkin (LTG) method that might be considered a
separate class. The LTG method handles the time domain with Laplace
transforms so time stepping and numerical oscillations are eliminated.
The LTG approach only applies to linear transport with steadystate
velocity fields (Allen et al., 1993).
14
2.1 Ground Water Flow Governing Equations and Model Type
Selection
Ground water models can be spatially categorized as follows:
Twodimensional areal
Twodimensional profile
Quasi threedimensional
Threedimensional
The timedomain run states are classified as follows:
Transient
Steadystate
Ground water flow conceptual models are divided into two main
systems, 1. the aquifer viewpoint and 2. the flow system viewpoint
(Anderson & Woessner, 1992):
1. Aquifer viewpoint
Twodimensional models
Both confined and unconfined
Stores and transmits water to wells
Uses leakage term for confining beds
Transmissivity in the z direction undefined
Governing equation form:
(2.1)
L = Kz
hih
m
15
+W = Outflow (W = Inflow) [ L/T] .
Kz = Vertical hydraulic conductivity [L/T]
hs = Source head [ L ]
m = Aquifer thickness [ L ]
S = Storage coefficient or storativity [ ]
h = head at x, y [ L ]
t = time [ T ]
Tx = Longitudinal transmissivity [ L2 / T ]
Ty = Transverse transmissivity [ L2 / T ]
L = Leakage term
The unconfined aquifer equation is defined as the Boussinesq
equation:
R= Inflow [L/T]
+R = Outflow [L/T]
Sy = Specific Yield [ ]
m = Thickness [ L ]
h = head at x, y [ L ]
t = time [ T ]
Kx = Longitudinal hydraulic conductivity [ L / T ]
Ky = Transverse hydraulic conductivity [ L / T ]
2. Flow System viewpoint
Threedimensional and twodimensional profile
Threedimensional distribution of:
 Head
 Hydraulic conductivity
 Storage coefficients
Governing equation form:
(2.2)
(2.3)
16
+W = Outflow (W = Inflow) [ L/T]
Ss = Specific storage [ 1 / L ]
h = head at x, y [ L ]
t = time [ T ]
Kx = Longitudinal hydraulic conductivity [ L / T ]
Ky = Transverse hydraulic conductivity [ L / T ]
Kz = Vertical hydraulic conductivity [ L / T ]
m = Thickness [ L ]
The socalled "Dupuit assumptions" are applied to the unconfined
aquifer in the aquifer viewpoint (Davis & DeWiest, 1966):
1. Flow lines are horizontal and equipotential lines are vertical,
2. Horizontal hydraulic gradient is equal to the slope of the
potentiometric surface to the aquifer base.
In the flow system viewpoint, the equipotential lines bend and pass
through all layers, and the heads are calculated in all layers.
Twodimensional models are used for aquifers confined,
confinedleaky, unconfined, and mixed. Head distribution and the
aquifer storage coefficient is required. In quasi threedimensional
models and twodimensional models of confinedleaky aquifers, the
leakage is represented as a function of a leakance term (Anderson &
Woessner, 1992) where:
Leakance = Kz/m
The confining layers are not explicitly represented and it is assumed that
there is no confining bed storage release. Pumping tests will provide
transmissivities and storage coefficients for confined aquifers, or specific
yield for unconfined aquifers for use in the Boussinesq equation 2.2.
Anisotropy of transmissivities can usually be input in terms of x and y and
aligned on model axes. If there is less than two orders of magnitude in
17
the difference between the aquifer and confining beds, a transient nature
in the water flux, or a deep aquifer compared to the twodimensional flow,
a full threedimensional model may be preferred.
Full threedimensional and profile models are used for unconfined
aquifers to simulate important vertical head gradients in the flow system
viewpoint. Some finite element codes are able to adjust the water table
boundary translational nodes during solution.
Profile models are a special twodimensional case where vertical
flow is important. These models are set along a flow line and assume
that all flow is in the plane of the profile and that there is no outofplane
flow. The governing equation 2.4 for a profile model is similar to the
above aquifer viewpoint equation 2.1 with some exceptions. The model
codes solve the previous aquifer viewpoint equation; so the user must set
up a conceptual transformation. The z parameters are substituted in for y
parameters and transmissivities are calculated by multiplying K by the
unit thickness of the profile. Recharge and storage coefficients must be
adjusted by unit thickness/delta y. The specific storage coefficient is
usually near zero, except at the water table nodes, where specific yield is
used. The following governing equation then applies to a profile model:
Ss = Specific storage [ 1 / L ]
+W = Outflow (W = Inflow) [ L / T ]
h = head at x, y [ L ]
t = time [ T ]
Kx = Longitudinal hydraulic conductivity [ L / T ]
Kz = Vertical hydraulic conductivity [ L / T ]
m = Thickness [ L ]
(2.4)
18
2.2 Solute Contaminant Transport Governing Equations
In general, mass transport with reactions for an element must meet
the law of conservation of mass:
{change in mass storage with time} =
+ {mass inflow} {mass outflow}
 {mass decay rate} + {mass production rate}
For uniform twodimensional transport the governing equation
utilizing dispersion/diffusion from Fick's Second Law and advection from
Darcy's Law in terms of mass balance* may be expressed as follows:
^ a2c ^ a2c ac ac
Dxy + Dyry Vx 
dx2 dy dx dt
(2.5)
With linear sorption and decay, the equation may be expressed as (Bear,
1979):
a2c a2c
Vxttx + VxCIy
Ridx2
RfSy2
ac ac
 VxA.C =
Rfdx dt
(2.6)
Adding in the source and sink term, the equation is (based on Goode &
Konikow, 1989; Goode, 1990):
ac a2c
dt VxCU Ridx2
i rcciw
Rimne
4VitXy
a2c
Rfdy2
dC
v*A.C
Ridx
(2.7)
19
C = Solution concentration [ M / L3]
C = Source solution concentration [ M / L3]
vx = Linear longitudinal velocity [ L / T ]
W = Outflow (W = Source) [ L / T ]
Rf = Linear sorption retardation factor [ ]
ne= Effective porosity [ ]
m = Aquifer thickness [ L ]
ax = Longitudinal dispersivity [ L ]
ay = Transverse dispersivity [ L ]
X = Decay rate constant [ 1 / T ]
The velocity vector field is supplied from the flow portion of the model.
The relationship of the flow and transport operations is shown in Figure
21. If fluid density changes during the simulation time, then the flow and
transport equations must be solved in conjunction with one another.
20
Figure 21 SystemLevel Flow Chart for Ground
Water Flow and Transport
21
2.3 Review of the Site Modeling Process
The first step in the modeling process is determining if there is a
need for mathematical modeling and what purpose a site model will
serve. See Figure 22. Once this has been defined, then the right math
model code for the job can be selected. The computer code model
should be verified, validated, and benchmarked, if not already done.
Next, a conceptual model must be developed from the site hydrologic
and geologic elements from available field data. This includes
identification of the lithology, stratigraphy, and structural features. Site
model setup involves selecting aquifer parameters, grid sizing/grid
placement, time steps, boundary conditions, and initial conditions. The
flow system can be defined from both physical hydrologic data and
geochemical data (Reily et al., 1987). Identification of geologic facies is
useful in defining hydrostratigraphic character. Once field data becomes
available, it should be compared to the model results. Calibration is the
triaianderror or automated process of adjusting a site model's input
parameters to provide results that match field measurement points. This
would include a mass balance check. The calibrated model still has a
degree of uncertainty because of the combination of unknown spatial
and temporal distribution of parameter values. A sensitivity analysis
should be performed to determine the amount of uncertainty in the
calibrated model by varying the parameters (Anderson & Woessner,
1992). The model is then used to simulate the future response of the
system. Predictive uncertainty exists in simulating magnitude, timing,
and distribution of responses to stresses. Next, is a model audit which is
a peer review and independent check. Site model verification refers to
checking the model code input and the entire preceding process. (The
term "verification has also been used to only refer to the step of checking
model reproducibility with different stresses.) Postaudits are done later
with new field data to see if the initial model predictions were correct.
22
Figure 22 The Site Modeling Process
O.Long S3
23
2.4 Establishing Ground Water Model Quality Assurance and
Control
The Quality Assurance (QA) plan should include a review of the
model selection, its conceptualized application, the calibration, a
sensitivity analysis, and the documentation. Implementation of Quality
Control (QC) requires the establishment of formal procedures for
checking the theoretical soundness and application of a model. Cost
and time may limit the degree to which QA/QC is implemented in context
with management modeling objectives. A Quality Assurance plan for the
development and application of models should include the following
points (based on Van der Heijde, 1987; 1990):
Objectives and quality level for
 Validity
 Verification
 Uncertainty
 Accuracy
 Completeness
Comparability
Documentation
 Versions
 Input
 Run history
 Problem definition
 Methodology
 Assumptions
 Sensitivity analysis
 Protocols for parameter estimation
Operational Standards and Procedures
 Personnel assignments
 Formal documentation
 Express limitations of model
 Validation/Verification of modeling codes
 Benchmarking of modeling codes
24
Internal and External Auditing
 Peer review
 Evaluation of calibration goals
 Verification check of input/output
 Independent model
 Predictions
2.5 Solute Contaminant Transport Model Evaluations
A project started in 1984 called the Hydrologic Code Inter
Comparison (HYDROCOIN), assessed ground water flow modeling
strategies and code comparisons (Beljin, 1988). In the project, no
significant solution differences were observed between finite difference
and finite element methods. Emphasis was placed on knowledgeable
and experienced users with good parameter data. Accuracy was found
to depend on selection of the proper model with its model theory and
correct definition of aquifer properties, boundary conditions, and initial
conditions. Misuse of models and lack of good data was discussed as
the subject of some controversy in the past. Lack of accurate field data
has made it difficult to correlate real case results to site models. The
HYDROCOIN project reported that full model validation may be elusive.
Comparisons of code results were ranked in order of general
represented accuracy:
1. Head and Pressure
2. Velocity
3. Flow Direction
4. Time Traveled
The EPA funded a study comparing GEOFLOW and MOC to
analytical models. This study as well as other studies concluded that
MOC was numerically accurate and most frequently matched transport
data.
25
A study done for the EPA evaluated the four following
saturated/unsaturated models:
SUTRA
PRZM
HYPE
FEMWATER/FEMWASTE
The study made comparisons between the model's predicted results and
limited field data. The evaluation concluded that SUTRA was the best for
most saturated/unsaturated flow problems. It should be noted that PRZM
has documentation especially suited to chemical concentration flow in
the unsaturated root zone. There are other codes that are designed for
dealing specifically with unsaturated transport.
A comparison was performed on four tested unsaturated flow and
transport models in a transient twodimensional homogeneous soil
section problem:
FEMWATER
UNSAT2D
SUTRA
VS2DT
From this comparison it was concluded that VS2DT was the most efficient
and accurate (Carovillano, 1993), although the other model results were
close.
There are other recent rankings of miscible fate models, however,
evaluations of these models have not gained much interest due to the
human factors of subjectivity and tradition. As with structural models, it is
unlikely that the large task of a conclusive comprehensive evaluation
would be undertaken in a totally objective manner which addresses the
merits of each code. As a concluding comment, even though some
computer codes are ranked lower in a comparison rating, that does not
always mean they are significantly less useful in providing desired
26
information one code does not perform all functions or do the best in
all scenarios.
2.6 Industry Model Usage and Needs
Geraghty & Miller (Rumbaugh & Ruskauff, 1993) conducted a
survey of ground water modelers in 1992 to determine the types of
software being used and the benefits and problems. A total of 876
surveys were received, with 67% of the respondents working as
consultants, 6% in industry, 10% in academia, and the rest in
government. Most of the respondents, 81%, were in model application;
and 8% were in management, 8% were reviewers, and only 3% were
oriented toward testing and software development.
Modelers identified the following top ten benefits of ground water
modeling in order of importance as;
1. Ability to predict future system response
2. Modeling requires development of a conceptual model
3. Evaluation/design of remedial alternatives
4. Graphical presentations
5. Ability to analyze "what if" scenarios
6. Analysis of system sensitivity/uncertainty
7. Models are good planning tools
8. Modeling forces the integration of all available data
9. Modeling identifies data requirements and deficiencies
10. Provides solutions to complex problems
The survey identified 194 commercial and public domain codes
and 180 proprietary software codes. Sixtyone percent of the codes were
used by less than five modelers. Out of the 194 commercial/public
domain codes, eight were identified by over 100 surveys. The top ten
software packages were identified as follows:
27
Software Function Public
MODFLOW Flow Model yes
MOC Transport Model yes
ModelCad Preprocessor no
SURFER Plotter no
PLASM Flow Model yes
SUTRA Transport Model yes
Random Walk Transport Model yes
AQTESOLV Well Test Curve Match no
PREMOD Preprocessor no
MODPATH Transport Model yes
Noteworthy codes used by some local consultants and experts
found in recent literature are listed as follows:
Software___ Function____________________________ Public
MODFLOW Flow Model, 3D FDM yes
PLASM Flow Model, 2D FDM yes
AQUIFEM1 Flow Model, 2D FEM tri elem yes
MODPATH Transport Model, Particle Track, 3D FDM yes
FLOWPATH Transport Model, Particle Track, 2D FDM yes
MOC Transport Model, 2D FDM yes
MOCDENSE MOC Transport Density Model, 2D FDM yes
BIOPLUMEII MOC Transport Biodegrad, 2D FDM yes
SEFTRAN Transport Model, 2D FEM no
RNDWALK Transport Model, 2D FDM yes
SUTRA Trans DensityTemp, 3D FEMIFDM yes
MT3D MOC Transport, Chem React, 3D FDM no
VS2DT Unsaturated Flow, 2D, IFD yes
The preeminent ground water modeling codes related to
contaminant transport are listed in Figure 23. This list is only intended to
be representative of some good codes, not entirely inclusive.
28
Figure 23 Some Premier Ground Water Modeling
Codes Related to Solute Contaminant Transport
Flow Codes:
MODFLOW 3D FDM PLASM 2D FDM AQUIFEM1 2D FEM AQUIFEMN 3D FEM
Transport Codes:
M0C2D FDM SUTRA 3D FEMIFDM Density & Heat MODPATH 3D FDM Particle Track FLOWPATH 2D FDM Particle Track RNDWALK 2D FDM MT3D 3D FDM MODFLOW link, Cham React HST3D 3D FEMFDM Density & Heat SEFTRAN 2D FEM BIOPLUMEII 2D FDM FEMWASTE 2D FEM MOCDENSE 2D FDM Variable Density
Unsaturated Tran snort Codes:
VS2DT2D IFDM
Special Purpose Codes:
HELP Landfill Seepage AQTESOLV Well Test Curve Matching PRZM2 Pesticide Root Zone ModelCad Preprocessor SURFERPlotter PREMOD Preprocessor TOUGH Multiphase Fracture Flow
Geochemical Codes:
MINTEQA2/PRODEFA2 WATEQ4F PHREEQE
Parameter Estimation Codes:
MODFLOWP GEOEASE MODMAN
V^Sng'BT
29
MODFLOW is a popular threedimensional finite difference flow
program with a "strongly implicit procedure" or "slice successive over
relaxation solution" technique. PLASM is a twodimensional finite
difference model using an interactive alternating direction implicit
procedure. AQUIFEM1 is a twodimensional finite element linear
triangle model using Crout's direct solution technique. AQUIFEMN is a
threedimensional version of AQUIFEM1. MODPATH is a three
dimensional semianalytic particle tracking code that is only for steady
state flow and uses heads input from MODFLOW. FLOWPATH is a two
dimensional finite difference code for steadystate problems with particle
tracking. MOC is a mature twodimensional finite difference flow and
transport code solved by the Method of Characteristics (MOC).
BIOPLUME and MOCDENSE are extensions to MOC for biodegradation
and variable density. Another mature hybrid FE/IFD flow and transport
code, SUTRA, is able to handle both variable density and temperature.
One major area under development is unsaturated flow. VS2DT is
a recently released program that has been recommended. It is a static
steadystate model that can handle nonlinearity. SUTRA can also be
modified and recompiled for unsaturated flow, but has limitations
common to coupled models. The HELP model is a useful tool for landfill
design (EPA, 1984b).
A new commercial solute transport code worth mentioning that has
a growing number of users is MT3D. This code is a 3D solute transport
model with chemical reaction capacity that uses an enhanced version of
MODFLOW to calculate flow. It is sold as a commercial program for
$450. MT3D has three solution options:
1. Method of Characteristics
2. Modified Method of Characteristics
3. Hybrid of the two above methods
30
The method of characteristics eliminates the need for numerically
calculating dispersion by finite difference. MT3D also has the standard
finite difference mode. It processes output files for plotting with programs
like SURFER. SURFER is a commercial contour mapping program that
uses Kriging.
Another new program is MODFLOWP. MODFLOWP is a three
dimensional inverse code that runs MODFLOW for statistical parameter
estimation and calibration. GEOEASE is a public domain code for
estimating variograms and performs Kriging for parameter estimation.
MODMAN is a parameter optimizer for MODFLOW.
Geochemical modeling codes compute the equilibrium of dissolved
species from point to point dealing with saturation, oxidation, reduction,
dissolution, precipitation, and ion exchange. Some prominent codes are
PHREEQE, MINTEQA2/PRODEFA2, and WATEQ4F. There are other
limited specialty geochemical codes for calculating chemical equilibrium
along a flow line.
In the Geraghty & Miller survey, ground water modelers expressed
their software needs, listed from the greatest to the least:
1. Postprocessors
2. Preprocessors
3. Quick, efficient, userfriendly models
4. Transport models
5. Unsaturated flow & transport models
6. Database interfaces for models
7. Inverse models
8. Multiphase models
9. Gas flow volatilization models
10. Fracture flow & transport models
11. Transport models with geochemical reactions
12. Expert systems for model design
13. Finiteelement flow and transport models
31
3.0 Solute Transport Analytical Solution
Ground water solute transport analytical solutions are available
that provide temporal and spatial solute transport in homogeneous and
isotropic systems. Classical solutions are primarily based on the
assumption that dispersion can be described by similarity to molecular
diffusion. Molecular diffusion is the mixing caused by random molecular
motion due to thermal kinetic energy. Fick's First Law states that
chemical mass flux is proportional to the concentration gradient. In a
simple aqueous system it is expressed as:
J = Dd{ VC} (3.1)
J = Chemical mass flux [ moles / L2 T ] or [ M / L? T ]
Dd= Fluid diffusion coefficient [ L2 / T]
C = Concentration [ moles / L3 ] or [ M / L3 ]
V = Gradient of scalar with respect to distance [ 1 / L ]
Transport is in the direction of decreasing concentration. It is assumed
that mechanical dispersion behaves like diffusion because the observed
macroscopic phenomena is similar. The coefficient of hydrodynamic
dispersion is then the sum of these two effects. Fick's law, when defined
in terms of mass transport (Fick's Second Law) and merged with
advective mass transport in conjunction with Darcy's Law, can be
converted to a form describing advective/dispersive flow in a porous
media. The concentration timeanddistance shape is based on a normal
distribution. The governing differential equation for dispersion by
conservation of mass assumes:
32
Dispersion behaves like diffusion
Homogeneous and isotropic media
No phase transfer
Constant average velocity throughout flow field
Dispersion/advection equation is then expressed as:
D = Hydrodynamic dispersion coefficient [ L2 / T ]
C = Solute concentration [ M / L3]
t = Time [ T ]
x = Longitudinal coordinate [ L ]
v = Longitudinal velocity [ L / T ]
V = Gradient of scalar with respect to distance [ 1 / L ]
In onedimensional form, rearranging for a semiinfinite medium and
source at x= 0, the dispersion equation 3.2 becomes (Freeze & Cherry,
1979):
ac _a2c ac
aT=Da?vaT (33]
A closed form "exact" solution for a onedimensional advection
dispersion flow source was derived by Ogata & Banks (1961):
(3.2)
The errorfunction is defined from:
33
The experimental BreakThroughCurve (BTC) results, at a distance from
the source, were matched with experimental data (Ogata & Banks, 1961).
The first term of the Ogata & Banks equation is commonly used
alone:
The second term is small for most cases and has been dropped by many
authors. However, it was discovered that some error may be introduced
near the source when ignoring the second term in certain value
combinations with large dispersivity. Dropping the second term removes
the true asymmetric shape of the BTC.
Analytical solutions have been published for two and three
dimensional transport in a uniform flow field (Wilson & Miller, 1978a,
1978b); Hunt, 1978; Bear & Verruijt, 1987; and Domenico & Schwartz,
1990). Analytical solutions are also available for solute transport with
linear sorption and radioactive decay. They can also be used in Monte
Carlo techniques. Hunt (1978) derived a continuous point source
equation, and Wilson & Miller (1978a, 1978b) derived a continuous line
source equation. The Ogata & Bank's (1961) equation 3.4, Wilson &
Miller's (1978a, 1978b) equation, Hunt's (1978) equation; and
contributions by Gershon & Nir (1969), Cleary & Ungs (1978), and Crank
(1979), were formulated into the Domenico & Robbins' (1985) three
dimensional advectiondispersion rectangular constant source equation:
(3.6)
(3.7)
34
x = Longitudinal distance traveled
y = Transverse distance traveled
z = Depth distance traveled
v = Longitudinal velocity
Y = Source width
Z = Source depth
Co = Initial concentration
ax = Longitudinal dispersivity
ay = Transverse dispersivity
z = Vertical dispersivity
erf = Error function
erfc = Complementary error function
exp = Natural exponential
A check was made of this equation by integrating the simpler Hunt point
source equation over an area of a finite source. Domenico & Robbins
(1985) found that there was a good match, with a small exception. Near
the source, within one source width range distance, the concentration
could differ; and this was by less than ten percent for the example case
tested. The point source equation and the rectangular source equation
must converge in the farfield. The Domenico and Robbins' equation
allowed adaptation to a grid for coding. The Y term was separated into y
+ Y/2 and y Y/2, and the Z term separated into z + Z/2 and z Z/2 terms,
where Y and Z are the source width and depth. When sources are at the
water table, Z replaces Z/2. The coordinate geometry is shown in Figure
31. As the transverse and vertical dispersivities approach zero, the
OgataBanks equation is returned. When these dispersivities or x are
zero, the last four terms for Z and Y are undefined for equation 3.7. The
first term is multiplied by 2 for each of the pairs of Z or Y terms dropped. If
z spreading is prevented as in a confined aquifer, the last pair of Z error
function terms are dropped and the first term is multiplied by 2. The
length that z spreading occurs in the x direction is (Domenico & Robbins,
1985):
35
(3.8)
x' = (H Z)2
Oz
When distances are greater than x', x in the z term denominator of
equation 3.7 is replaced with x' from equation 3.8.
Domenico incorporates the decay arguments (Bear, 1979) into the
above equation 3.7, resulting in the following equation:
x = Longitudinal distance traveled
y = Transverse distance traveled
z = Depth distance traveled
v = Longitudinal velocity
Y = Source width
Z = Source depth
Co = Initial concentration
X = Decay rate constant
ax = Longitudinal dispersivity
ay = Transverse dispersivity
az = Vertical dispersivity
erf = Error function
erfc = Complementary error function
exp = Natural exponential
C(x,y,z,t) =
(3.9)
36
Other forms of this equation are presented in Bear (1979), De Marsily
(1986), and in Goode & Konikow (1989). This equation 3.9, including the
second term in equation 3.4, plus retardation and decay, was coded
using FORTRAN into the program TACT ThreeDimensional Analytical
Solution Contaminant Solute Transport Code for Saturated Matrix with
Dispersion and Retardation (Long, 1993).
A special case of the equation 3.7 along the centerline plane of
symmetry at the water table, which may be referred to as the Centerline
ThreeDimensional Dispersion Equation, is as follows:
Equation 3.10 can be used for hand Calculation checks of equation 3.7.
This equation 3.10 was coded in TACT as a check against the grid
equation. In addition, the onedimensional equation 3.4 was coded in
TACT and was compared against onedimensional BTC's.
The analytical equation error functions were checked
independently against published values from equation 3.5. These
results are graphed in Figure 32. The TACT code threedimensional
results were independently verified against the governing equations
programmed into MATHCAD, which calculated identical results.
(3.10)
37
Function
Figure 31 Analytic Solution Transport Coordinate
Geometry
O. Long S3
Figure 32 Error Function & Complementary Error
Function Check Plot
38
4.0 Solute Transport Numerical Solution Modeling
Techniques
Transport numerical modeling techniques typically employed are
finite difference, finite element, or combinations of these. These methods
use time and space discretization to solve the differential equations.
Finite element models define the variation of head within an
element by means of a shape (interpolation) function. The finite element
method has an advantage in some cases with better element simulation
of fault zones, moving water tables, point sources, and fitting of irregular
shaped aquifer units and boundaries. Finite difference models compute
the head at the element center node representing an average of the
surrounding element area. The finite difference method might be
considered a special case of the finite element method.
MOC is a finite difference model and SUTRA is a hybrid model
using both finite element and finite difference methods. MOC and
SUTRA were selected for use in this thesis because they were prominent
in meeting the approach criteria laid out in Section 1.4.1. Supporting
reasons for selecting these codes can be seen from the particular
attributes discussed in the following sections.
39
4.1 SUTRA Model Description
SUTRA (Voss, 1984) is the acronym for the Saturated/Unsaturated
Transport computer program that simulates fluid movement and transport
of solutes or thermal energy in a saturated environment. This model
employs a quasi twodimensional hybrid integrated finite difference (IFD)
method with a finite element (FE) method using weighted residuals.
SUTRA is capable of calculating concentration and thermal transport or
variation with time in a subsurface system. It was one of the first
programs to provide timedependent transport solutions. SUTRAPLOT
provides for plotting the output.
The FE method is employed for the fluid and solute mass balance
equations. All remaining nonfluid terms are calculated using IFD
methods from the mesh. An advantage of the FE method is that irregular
quadrilateral elements with irregular thicknesses can be used (there may
be some transient storage error).
SUTRA is designed primarily for twodimensional subsurface flow.
This is usually sufficient for most problems because the media thickness
is much less than the other two dimensions. In SUTRA, the three
dimensional aspect is introduced by adding thickness to a two
dimensional model, which functions as a pseudo thirddimension. The
spatial coordinate system may be specified as Cartesian or Cylindrical.
Solute transport is simulated through numerical solution of solute
mass balance where solute concentration is allowed to affect fluid
density. Sorption of solutes may be represented by Linear, Freundlich,
or Langmuir isotherms. In a like manner, thermal energy transport is
simulated using the numerical solution of an energy balance equation.
Introduction of solute mass volume and decay is also allowed. Sources,
sinks, and boundary conditions of the fluid flow can be specified as
constant or variable with time.
While the title SUTRA implies unsaturated as well as saturated
transport, unsaturated analysis is not a part of the delivered program in
40
executable form. Unsaturated runs require that the program code be
modified and recompiled, and is not intended as a primary application of
the code. The model dimensional and time discretization must be very
fine to satisfactorily simulate unsaturated flow. There are other
unsaturated zone model codes that can be used, such as EPACML, a
composite semianalytical model, or VS2DT, an integrated finite
difference model code.
The standard means of hydrodynamic dispersion for isotropic
media assumes direction independence for longitudinal and transverse
dispersivity. SUTRA also allows for velocitydependent dispersion for
anisotropic aquifers.
The finite element mesh is composed of quadrilateral elements.
Each of the four comers has a node representing the midpoint of a z
edge (see Figure 41). In SUTRA, constant parameter values are
assigned to each mesh node, either elementwise, nodewise or cell
wise. Cells are considered to be centered on a node. Cell boundaries
are half way between element sides (Figure 42). The spatial
distributions of parameters are assigned as follows:
 K elementwise
 h nodewise
 S cellwise
41
f:'9ure *
ent ($*'fjmt
SS' ,9^6 e'ein, Me
Sh
9na
42
Figure 42 Definition of TwoDimensional Quadrilateral
Elements, Cells, and Nodes (Voss, 1984)
43
The SUTRA governing equation is based on compressible fluid
flow, and its approximate simplified solute form is as follows (based on
Voss, 1984; Goode, 1990):
ac_ 1 a pmneDij dxj ac
dt pRunne dxi " Vi Ridxi
+(cc)w
pRimne
(4.1)
pile = (pile)o + poSop[p po]
C = Solution concentration [ M / L3 ]
C' = Source/sink solution concentration [ M / L3 ]
Cs = Sorbed solution concentration [ M / L3 ]
vj = Linear longitudinal velocity [ L / T ]
W = Outflow flux (W = Inflow) [ L/T]
Rf = Linear sorption retardation factor [ ]
ne = Effective porosity [ ]
m = Aquifer thickness [ L ]
X = Decay rate constant [ 1 / T ]
Oijmn = Dispersivity [ L ]
Dy = Dispersion coefficient [ L ]
ij = Dimensional subscripts 1,2 [ ]
m,n = Directions of flow subscripts [ ]
v = Velocity vector magnitude [ L / T ]
p' = Source/sink solute density [ M / L3]
p = Solute density [ M / L3]
Po = Initial solute density [ M / L3]
Sop = Specific pressure storativity [ 1 / L ]
h = Head [ L ]
t = time [ T ]
Density changes to water caused by Total Dissolved Solids (TDS)
or temperature can affect flow results. In the past, solute density was
considered to not be a problem in most aquifer water, except for salt
water and some high concentration landfill leachates. However, recent
44
field studies have indicated that densityinduced flow can be a significant
factor in some cases, even for tracers at low concentrations, and
especially for buoyancy induced flow (Istok & Humphrey, 1993).
Concentration density should be considered when [(Ap/p) / (dh/dx)] >
0.33. The magnitude of this effect can be seen in the following example:
Salt water has a TDS of 36,000 mg/L for a density of 1.025 gm/cm3. With
a slope of 0.025, the force in the horizontal direction is as great as the
force in the vertical direction. The resulting flow force and direction (at
45, not lateral) is a vector sum of these two forces. Usually the higher
vertical conductivities will restrict this vertical component flow to the
aquifer. Also, seasonal variations in water temperature can change
viscosity considerably, affecting hydraulic conductivity.
45
4.2 MOC Model Description
MOC (Konikow & Bredehoeft, 1978) is a twodimensional finite
difference flow model with solute transport in a homogeneous saturated
aquifer. MOC uses a blockcentered (node in the center) finite difference
grid for flux and transport calculations. With constant head, it can .
calculate the spatial distribution of concentration with time as caused by
advection and hydrodynamic dispersion, and chemical reactions. MOC
has two versions that allow solving finite difference flow with the (1)
Alternating Direction Implicit (ADI) procedure, or (2) Strongly Implicit
Procedure (SIP). MOC uses the method of characteristics to solve the
solute transport equation with particle tracking. Transport is simulated by
particle movement which is uncoupled from the flow portion of the
program to prevent numerical dispersion. This is a complicated process
of using a "material derivative" from a set of ordinary differential
equations that are equivalent to the transport equation. This allows each
point to migrate at a distance proportional to the time increment length
and to the velocity at the point location. MOC's cell definitions for particle
velocity are shown in Figure 43, and MOC's algorithm is shown in Figure
44.
MOC contains the following chemical reactions in a single
advectiondispersion solute transport equation:
Firstorder, irreversible ratereaction (radioactive decay)
Reversible equilibrium reactions
 Linear sorption
 Freundlich sorption
 Langmuir sorption
 Monovalent ion exchange
 Divalent ion exchange
 Monovalentdivalent ion exchange
 Divalentmonovalent ion exchange
46
Figure 43 Definition of MOC Particle Movement in
Relation to Grid & Node (Konikow & Bredehoeft, 1978)
EXPLANATION
Initial location al particle
O New location ol particls
Flow line and direction of (low
Computed path of particle
Part of hypothetical finite
difference grid showing relation of
flow field to movement of points.
X
*',/1 1.71
ji.i /,/ W vM.
1
Vi,/.I J^/*1 T */1. /!
EXPLANATION
Node of finitedifference grid
o Location of particle p
 X or V component of velocity
Area of influence for interpolating velocity
in X direction at particle p
Area ol influence far interpolating velocity
in f direction at particle p
Part of hypothetical finitedifference grid
showing areas over which bilinear Interpolation Is
used to compute the velocity at a point. Note that
each area of Influence Is equal to onehalf of the
area of a cell.
EXPLANATION
' Node of finitedifference (pld
Location of particle
Halation between possible Ini
tial locations of points and Indlcas of ad
jacent nodes.
47
Figure 44 Main Level MOC Algorithm
(Konikow & Bredehoeft, 1978)
48
Sorption is explicitly incorporated into the particle tracking using a
linearized velocity retardation for each time step after accounting for
advection and decay (Goode & Konikow, 1989).
MOC uses the "flowequationremoved" form of the transport
equation and is considered to be "fluidmassconservative" (Voss 1984).
This is stated as follows (based on Konikow & Grove, 1984; Goode &
Konikow, 1989; Goode, 1990):
3C
dt
_1_____3_
Rfirnie dxi
imiaDii
ac
Jaiy
ac
Ridxi
_3h 3m
S + Wnc
dt dt
CW
Rfmne
XC
(4.2)
where:
Dij (Xijmn
VmVn
mn = (mne)o + S[h ho]
Tij ah
nine dxj
49
C = Solution concentration [ M / L3]
C" = Source/sink solution concentration [ M / L3 ]
vj = Linear longitudinal velocity [ L / T ]
W = Outflow flux (W = Inflow flux) [ L / T ]
Rf = Linear sorption retardation factor [ ]
ne = Effective porosity [ ]
m = Aquifer thickness [ L ]
X = Decay rate constant [ 1 / T ]
ottimn = Dispersivity [ L}
Dy = Dispersion coefficient [ L2 / T ]
ij = Dimensional subscripts 1,2
m,n = Directions of flow subscripts
v = Velocity vector magnitude [ L / T ]
pb = Bulk density [ M / L3 ]
S = Storage coefficient [ ]
h = Head [ L ]
ho = Initial head [L ]
t = time [ T ]
The last term is modified for the condition of nonlinear sorption. The
MOC saturated volume is adjusted in the most recent version for transient
flow to account for changes in thickness or porosity that in turn affect the
flow head and velocity calculations. For a flux boundary condition, solute
travel time is affected by mass storage errors, but very little by
transmissivity that is dependent upon head (Goode, 1990). This same
equation (4.2) can be used in Lagrangian models.
50
4.3 Practical Considerations of Numerical Modeling
4.3.1 Model Grid Scale and Grid Types
The mass balance equation that includes solute dispersion is
difficult to solve numerically. It can suffer from numerical dispersion
errors due to discretization. Backwards finite difference schemes,
considered mass conservative, have truncation errors on the order of (Ax
grid)2. These introduce numerical dispersion unless the mesh size is
very small (Axgrid 2 x Dispersivity). Centered finite difference
schemes, considered not mass conservative, have a truncation error on
the order of (Axgrid)3. These are subject to upstream and downstream
oscillations unless the grid size is small (Axgrid < 2 x Dispersivity) (De
Marsily, 1986; Pisani & Tosi, 1993). However, this size is not always
practical, so ways of dealing with these problems have been found in
some code formulations. In many situations, the characteristic grid length
should be no longer than ten times the characteristic dispersivity, with
four times preferred (Anderson & Woessner, 1992). The time step should
be selected so that the time step times the velocity is less than the
characteristic grid length (Axgrid). The element aspect ratio should be
no greater than 1.5.
A blockcentered finite difference element grid has flux portions
along the boundaries on the edge of the blocks and specified head at the
center of the element. Meshcentered finite difference elements and
finite elements have both flow and specified heads at the nodes located
at element grid intersections.
51
4.3.2 Boundary Conditions, Sources, and Sinks
Boundary condition assignment is the most delicate portion of
modeling. The flow pattern of steadystate simulations is primarily
determined by boundary conditions. When transient simulation stresses
reach the boundary, equation solutions may be affected adversely.
There are three main types of boundaries:
Physical boundaries
 Impermeable layer
 Surface water
Distant boundaries
Hydraulic boundaries/artificial boundaries
 Divides
 Streamlines
Boundaries are selected so as to shrink the grid domain. These
boundaries can be:
Specified head
 Dirichlet conditions
Specified flow flux including noflow
 Neumann conditions
 Function of flux derivative
Headdependent flow boundaries
 Cauchy of mixed conditions
 Flux calculated from head
Care must be taken to define conceptually realistic boundaries and ones
that the models can handle numerically.
A specified head boundary requires a large supply of water to
maintain the head. The modeled site may draw or discharge water from
the constant head cell, in a steadystate situation, the heads are set at
the water table elevation. An example of a site model is shown in Figure
45.
Specified flow boundaries are simulated in the model by using
injection or extraction wells. Use of flux boundaries exclusively should
52
Figure 45 Example of Model Boundary Conditions
(Konikow, 1977)
0 1 2 MILES
1 11I1
0 12 KILOMETERS
EXPLANATION
Zerotransmissivity call
Constanthead cell
A Disposalpond call.
I Irrlgationracharga eall
L CanaMaakage cell
 Boundary or areas in which alluvium is absent or unsatureted
53
be avoided because these are derivatives of head and may allow a non
unique solution (Anderson & Woessner, 1992). At least some head
values are required for steadystate problems. Transient solutions have
initial heads, so flux boundaries could be justified. Noflow flux
boundaries are placed at an impermeable bedrock, a ground water
divide, a fault zone, or a streamline. Some models automatically assume
noflow boundaries but other models must have them defined with zero
conductivity.
Headdependent conditions are used to simulate leakage to and
from a river or lake and can occur not only at a boundary, but also in the
grid interior as sources and sinks. A "losing reach" of stream recharges
the aquifer and "gaining reaches" reduce aquifer flow volume.
Hydraulic boundaries can be used to produce a steadystate flow
field for calibration purposes (Anderson & Woessner, 1992). Distant
boundaries are located far enough out so as to not affect the stressed
area of transient problems (as seen in Ward, 1987).
An impermeable boundary may be defined as a layer where the
hydraulic conductivity increases by a twoorderofmagnitude. This
contrast in conductivity causes refraction of flow lines or parallel line flow
in the higher conductivity layer. The flow in the lower conductivity layer is
primarily perpendicular to this. In this case, ignoring the confining layer
horizontal flow will result in less than 5% error in head (Anderson &
Woessner, 1992). In profile seepage models, the lateral boundaries
should be at least three times the depth of the flow system (3(KX/KZ)1/2 in
anisotropic media).
A modeling difficulty is simulating recharge to the aquifer. One
way is to use an unsaturated zone model to find infiltrated water reaching
the water table. Infiltrated water will eventually reach the ground water
unless evapotranspiration is significant. Unsaturated zone models are
not always used due to lack of required field parameters and complexity.
Timed output from an unsaturated onedimensional model could be used
as input to a saturated model. Coupled models could account for
54
movement of the water table with recharge. Most fieldscale soil profile
models with both unsaturated and saturated regions may not be
compatible (Frind & Verge, 1978). This is due to the smaller spatial and
temporal discretization required for unsaturated flow. An apparently
successful coupled model is presented by Fedors et al. (1993).
The mass balance (water budget) that must be represented in the
flow model is summarized as follows:
Sources inflows & recharge
 precipitation
 overland flow
 surface water bodies
 flow from streams
 some underflow
Sinks outflows & discharge
 flow to streams
 evapotranspi ration
 pumping
 spring flow
 some underflow
These sources and sinks may be placed at the boundaries or throughout
the grid. Flux points are represented by pumping or injection well nodes.
In twodimensional areal models and quasi threedimensional models,
an injection node represents the full depth of the aquifer. Head gradient
distribution inside a finite difference well cell are not accurate because
the water addition/extraction is made to the entire cell area rather than to
a point. The Thiem equation can be used to estimate head at the well for
steadystate and even transient simulations, because the effects of
storage are quickly overcome. Recharge is also a source in two
dimensional models. Water table recharge may constitute a boundary
condition rather than only a source term in some profile and three
dimensional models that do not use the Dupuit assumptions (Anderson &
Woessner, 1992). Interior sources and sinks may not be strictly
considered boundary conditions. Specified heads may also be placed
within the grid interior to represent lakes and rivers.
55
4.3.3 Transient Flow Numerical Modeling
Transient models are flow simulations which are time dependent.
Transient cases have special problems because many program model
packages have assumed the boundary conditions do not change
significantly from initial conditions. Transient conditions may require
heads to vary at hydraulic boundaries. Most programs do not allow the
heads to change at these boundary points and can cause error at the
interior of the gridwork. A larger model grid of the flow region containing
a submodel with a denser mesh may be required.
Transient models require initial head conditions and then
calculation of new heads with each time step. The transient models may
be broken into four main types as follows:
Dynamic average steadystate initial conditions
 Head varies spatially and temporally
 Flow in = flow out
 Boundary conditions can cause influence problems
Steadystate initial head conditions
 Head constant
 Assumes no flow
 Good for datum superposition in drawdown simulations
 Influence of stress must not reach boundary
Dynamic cycling initial conditions
 For cyclic water table fluctuations
Arbitrary head distribution
 To match field measured heads after some time
The critical time step is Ate = S(Axgrid)2/4T (Anderson &
Woessner, 1992). Convergence to the analytical solution can be
achieved in approximately six time steps without using small time steps.
When the storage coefficients of some models are set to zero, a steady
state solution can be obtained in one time step. Transient leakage is not
represented in most flow codes.
56
4.3.4 Aquifer Parameter Estimation
The hydraulic conductivity can be determined from pumping tests
and slug tests of the stratigraphic layers. Lab test results are typically
several orders of magnitude smaller due to absence of larger scale
features with locally higher transmission and heterogeneous macro
dispersion that increases the larger scale hydraulic conductivity.
Horizontal anisotropy of hydraulic conductivity (Kx/Ky) may be estimated
from field measurements. Horizontal anisotropy of hydraulic conductivity
may vary by two orders of magnitude and vertical anisotropy (Kx/Kz) by
three orders of magnitude. Vertical anisotropy can be estimated from
pumping tests, but is usually unknown and is estimated and adjusted
during calibration. The model axes should be aligned along these
conductivity axes if possible. The equivalent hydraulic conductivity for
flow in the direction of a layered section is a weighted arithmetic mean,
and for flow at right angles to layering is a weighted harmonic mean. A
geometric mean is used for random spatial heterogeneity, which typically
occurs as a lognormal distribution. Porosity has a normal distribution.
Evapotranspiration can be estimated by lysimeters or vegetation types.
Storage coefficients and transmissivity can be determined from
well pump tests. The specific storage can be estimated from coefficients
of vertical compressibility (Domenico & Schwartz, 1990), and is usually
for a confined aquifer. In water table aquifers or unconfined conditions, it
is commonly assumed that storativity or storage coefficient is equivalent
to specific yield and this is approximately equal to the effective porosity.
In a steadystate condition, the storage coefficient is assumed to be zero.
Since field data is usually sparse, statistical, probabilistic, or
interpolation methods can be used for estimating the spatial distribution
of parameters. However, there are some objections to these methods.
One is that these methods are often used as a substitute for more realistic
geologic interpretations. Another is the lack of deterministic results.
Questions about possible mathematical difficulties have also been raised
57
for socalled "inverse" methods of calibration (GhohestaniBojd, 1990;
Van der Heijde, 1990). Kriging is a commonly used statistical
interpolation method with variograms that preserves the field
measurement data points, unlike other Stochastic methods. While these
methods may be useful tools, studies have shown that different types of
Kriging or Stochastics can lead to conflicting results due to different
parameter biasing (McKenna & Poeter, 1993). For example, results of
two Stochastic methods, the parametric turning bands method and the
nonparametric sequential indicator simulation method, gave divergent
spatial distribution results, although the arrival times were similar
(Ruskauff & Field, 1993). The potential for misuse of calibration methods
is great (Reddi, 1990); and there is no strong proof yet that these
procedures are better than "trial and error." The potential for non
statistical interpolations was not investigated. The literature indicates
that good model calibrations to field data do not always lead to good
model predictions (Freyberg, 1988). This is because the inverse problem
is nonunique and derived conductivities are sensitive to observed head
noise and errors.
58
5.0 Benchmarking Case Studies
The purpose of the two case studies is to evaluate the accuracies
of MOC and SUTRA numerical solution FE/IFD models with an analytical
solution model, TACT, in solving an areal transient (timedependent)
solute transport problem in a steadystate saturated flow field with
dispersion, sorption, and decay from a continuous rectangular source.
The flow field for the two cases are defined as steadystate two
dimensional or quasi threedimensional. This can be done because the
ratio of the field length to aquifer thickness is very large. Ogata & Banks'
and Domenico & Robbins' analytical equations were coded into TACT as
described in Section 3. The results in each case were plotted for
comparison with MOC and SUTRA.
The program code versions are:
TACT
 ThreeDimensional Analytical Solution Contaminant Solute
Transport Code for Saturated Matrix with Dispersion and
Retardation
 Version 1.1
 October 1993
 FORTRAN (M.S. 4.1)
SUTRA
 Saturated/Unsaturated Transport
 Version V06902D
 June 1990
 FORTRAN
MOC
 TwoDimensional Solute Transport Method of Characteristics
 Version IGWMCFOS 07 PC Version 3.0
 November 1989
 FORTRAN (M.S. 4.1)
59
5.1 Case Study/Benchmark 1
The Case 1 scenario is a generic model test case contrived in a
manner to exemplify the relative differences in the code computational
attributes. The aquifer values and scales are realistic and similar to a
SUTRA test case.
5.1.1 Problem Definition for Case Study/Benchmark 1
The problem analyzed was a continuous rectangular trench
source 1,000 feet wide at the water table with a concentration of 1,000
ppm for 16 years in 4 year increments.
The site model grid plan is shown in Figure 51 for SUTRA and
Figure 52 for MOC. The boundary conditions were 200 feet of constant
head along the upper entrance boundary and 50 feet of head along the
exit boundary. The MOC heads were adjusted for the shorter total grid
length to give approximately the same gradient. A constant gridinterior
source concentration was applied at an equivalent element width of
1,000 feet at node 264, in SUTRA as shown in Figure 51, and in MOC as
shown in Figure 52. The input data is summarized in Tables 51 and
52. Each subcase varies one or more of the parameters relative to the
other subcases. SubCase 1E could be considered baseline, because
the numerical and analytical source volume are equal.
The following parameter values were used in this case with
variations for five subcases:
Hydraulic conductivity: K = 5.0 x 104 feet/second
Porosity: n = 0.20
Retardation constant: Rf = 2.00
60
Hydraulic gradient:
i = AH
L
= 20050 =0.0075
20,000
Darcy's velocity:
q= Ki = 5.0x1 O'4 (0.0075) = 3.750x106 feet/second
Darcy's volumetric glow rate:
Q = A K i = 50 x 1,000. x 3.750 x 10'6 = 1.8750 x 10'1 feet3/second
Water velocity:
vw = = 3.750 x 106 = 1.875 x 10'5 feet/second
n 0.20
Contaminant velocity:
v = vffi = 1.875 x 105 = 9.375 x 106 feet/second
R, 2.0
Dispersion coefficients: a* = 200 feet
a, =50 feet
Time:
t = 16 yr x 365.25 day/yr x 24 hr/day x 60 min/hr x 60 sec/min
= 5.0492 x 10 seconds
Concentration: Co = 1,000 ppm = 0.001 pound/pound water
The relationship between the retardation factor, Rf, and the sorption
isotherm distribution coefficient, Kd, is as follows:
Rf. 1+ l^]p.Kd
. n J
or
Kd =
n
Hp*.
(Rf1)
(5.1)
61
centerline concentration is too low near the source. The other high
volumetric flow subcases have high centerline concentrations. For a
large source, MOC's centerline concentrations for decay (SubCase 1D)
from the advective center to the leading edge are known to be high
because they exceed the analytical solution onedimensional results in
this region. See Figure 513. MOC's centerline concentration at some
distance from the source to the advective center (advective front) is also
high, but the degree is less certain.
67
5.1.4 Case Study/Benchmark 1 SUTRA Model Results and
Discussion
Several different execution option combinations are available for
SUTRA. The different SUTRA execution options that were tried are as
listed in Table 53.
Table 53 SUTRA Execution Options
SUTRA Run Option Types Constant Concen tration Source Flow Damping Steadystate Solute Transport
Type A X
Type B X X
Type C X X X
Type D X X
Type E X* X
Type F X X X X
D. Long 93
* Flow concentration increased to achieve approximate source
concentration when added to flow volume.
It was found that a constant source concentration may be required
when the source flow is low relative to the flow capacity per cell in a
generally advective field. Damping is required in this case to prevent
sidelobe and upstream concentration oscillations. Steadystate
transport is good for long time durations beyond transient effects. It was
determined from calibration testing that Type C is the best execution
option for this case.
The SUTRA threedimensional and longitudinalvertical slice
profile solution plots for SubCase 1A and 1E are shown in Figures 514
through 517. The SUTRA model results for centerline concentrations
comparing all subcases are presented in Figure 518. In SubCase 1E
SUTRA is showing a tendency to extend its leading edge.
68
Figure 53 Case 1A Analytic ThreeDimensional Plot
E
Q.
a
o
(0
c
0)
a
c
o
O
Figure 54 SubCase 1A Analytic ThreeDimensional
Longitudinal Slice Profile Plot Longitudinalvertical slices
from the preceding threedimensional plot.
Range Distance 1000's Feet
69
Concentration ppm Concentration ppm
Figure 55 Case 1 Analytic Centerline Profile
Concentration Subcase Comparison Plot
Figure 56 Case 1 Analytic Comparison of One
Dimensional and ThreeDimensional Profiles
70
Figure 57 SubCase 1A MOC ThreeDimensional Plot
c
o
+*
Se
C Q.
8
C
o
o
Figure 58 SubCase 1A MOC ThreeDimensional
Longitudinal Slice Profile Plot Longitudinalvertical slices
from the preceding threedimensional plot.
Range Distance 1000's Feet
71
Figure 59 SubCase 1E MOC ThreeDimensional Plot
Range
Distance
1000's Feet
o
<9
i_
C
0)
o
c
o
E
Q.
a
Figure 510 SubCase 1E MOC ThreeDimensional
Longitudinal Slice Profile Plot Longitudinalvertical slices
from the preceding threedimensional plot.
Range Distance 1000's Feet
72
Figure 511 Case 1 MOC Centerline Profile
Concentration SubCase Comparison Plot
73
Concentration ppm Concentration ppm
Figure 512 Case 1 MOC Source Volume Sensitivity
Centerline Profile
Figure 513 Case 1C & 1D MOC Comparison to
Analytic OneDimensional Centerline Profiles
1000
900
800
700
600
500
400
300
200
100
0
Case 1C MOC
* e Case 1D MOC Case 1C Analytic 1 D
i,\ N
\
\ i o Case 1D Analytic 1 D .
\ <
KV LJ
L\ 1
Lj
*i H C i 8 i O
5 6 7 8 9 10 11 12 13 14 15
Range Distance 1000's Feet
74
Concentration ppm
Figure 514 Case 1A SUTRA ThreeDimensional Plot
c
o
C o.
0) &
o
c
o
o
Figure 515 SubCase 1A SUTRA ThreeDimensional
Longitudinal Slice Profile Plot Longitudinalvertical slices
from the preceding threedimensional plot.
Range Distance 1000's Feet
75
Figure 516 SubCase 1E SUTRA ThreeDimensional
Plot
Figure 517 SubCase 1E SUTRA ThreeDimensional
Longitudinal Slice Profile Plot Longitudinalvertical slices
from the preceding threedimensional plot.
Range Distance 1000's Feet
76
Figure 518 Case 1 SUTRA Centerline Profile
Concentration SubCase Comparison Plot
Range Distance 1000's Feet
77
5.1.5 Case Study/Benchmark 1 Comparison and Conclusions
A comparison of centerline concentration profiles for SubCases
1A, 1B, and 1C are plotted in Figure 519. Comparisons of centerline
profiles for SubCases 1D and 1E are plotted in Figure 520. In Sub
Case 1A, SUTRA has good agreement with the analytic solution while
MOC has a higher centerline concentration. SubCase 1 B's source
volume is half of SubCase 1A. SubCase 1C, with higher dispersivity
(and high source volume), has the most pronounced concentration
deviation of the numerical solutions from the analytical solution (see
Figure 519, SubCase 1C). Generally a finer grid relative to dispersivity
is purported to mathematically give better numerical results. SubCase
1D dispersion spreading length increases only by the square root of the
dispersivity. SubCase 1E numerical models have a source volumetric
flux equal to the analytic solution average cell flow.
For SubCases 1A 1D, SUTRA's centerline concentration results
at the advective center tend to fall behind MOC's and closer to the
analytical solution's advective center. SUTRA's leading edge
longitudinal and transverse concentration is extended, and this
occurrence is observed in all subcases and execution options. Some
dispersion extension is to be expected due to artifacts of damping, but
sensitivity runs indicate this is not the only cause. A cutoff could be
selected at some low percentage (<10%). MOC's centerline
concentration in SubCases 1A 1D tends to elevate the concentration
between the source and the advective center (advective front) and its
advective center is ahead of the analytical solution's advective center. It
is surprising in subCase 1A, that although the shapes of SUTRA and
MOC are different, the total spatialconcentration volume integration
(SubCase 1 A) differ by only 5%. However, their volumes deviate by
26.3% in SubCase 1E. They are both significantly larger than the
analytical solution spatialconcentration volume for SubCase 1A, as
would be expected with higher source volumes. In SubCase 1E, MOC's
78
spatialconcentration volume is 8.1% higher, but SUTRA is 36.6% larger
than the analytical solution spatialconcentration volume. Here in Sub
Case 1E, SUTRA's leading edge protraction problem is more evident.
SUTRA's manual warns of a longitudinal dispersivity increase of up to
onehalf of the grid length when damping is used, especially in highly
advective velocity fields. The concentration contour plots for SubCase
1E are shown in Figure 521. Both of the numerical codes' outer
contours are conservative here.
SUTRA and MOC execution options and parameters were varied
for sensitivity assessments. It was discovered for this case that source
volume significantly affects the results. The reason for this is speculated
to be that the injection node requires enough source volume to generate
volumetric flow and velocity which is initially dependent upon flow field
directions and gradients into the adjoining two to four elements. This is
why the higher volumetric flow rate is used in some subcases. The
analytical solution assumes the volumetric flow rate in all cells are equal
to the source flow rate (divided by its number of cells) and has no other
mass accounting for source volume relative to the aquifer. The analytical
solution deals directly with threedimensional advection and dispersion
movement from this control area. A control volume investigation of both
the theoretical numerical and analytical model equations for this effect
would be helpful in resolving ground water model differences. A three
dimensional or a profile model in the vicinity of the source may be useful
in describing this phenomena. Studies of model dimensional effects
have been performed by Burnett and Frind (1987), and were found to be
significant. However, these source difficulties do not invalidate the two
dimensional assumptions for the farfield where the length to aquifer
thickness is large.
While the analytical solution should be considered the most
accurate in the farfield, it should not be construed as complete "truth"
near the source. In development of the dispersion term in the analytical
advectiondispersion equation, the coefficient of mechanical dispersion
79
is combined with the coefficient of molecular diffusion. In most cases the
ground water velocity is high enough that the mechanical dispersion is
the dominant effect. A Fickian model requires a Gaussian distribution,
but many test observations reveal a true skewed timeconcentration
curve near the source. Contaminant distribution is proportional to the
mean travel distance rather than to the square root of travel distance,
which is currently used (EPA, 1989b). Thus, the actual spread could be
greater than predicted by the analytical solution. In addition, it can be
observed that upstream flow is possible due to diffusion in slowvelocity
fields or due to vertical flow gradient dispersion from an elevated source.
In a revealing discussion by De Marsily (1986), he states that dispersivity
should not be defined as a function of the medium or time, but as
something else, such as a function of the variance in a Lagrangian
velocity field. There is no compelling reason why a Fickian model with
Gaussian distribution should be expected nearfield. This is because
dispersion mixing is confined to two dimensions far from the source, but
is in threedimensional expansion near the source. Therefore, the
Eulerian velocity field is only valid for farfield and later transport times
where asymptotic behavior takes place. In fact, Ogata and Banks (1961),
in their analytical solution to the mass transport differential equation,
describe the actual asymmetrical nature of the breakthrough curve near
the source. Gershon and Nir (1969) discuss errors due to boundary
conditions. A new analytical dispersion equation is needed for nearfield
and early transport times. Analytical equations differentiating nearfield
and farfield concentrations have been proposed for short duration slug
releases (Baetsle, 1969; and Charbeneau et al., 1992). It is suggested
that a similar adaptation be developed for the continuous source
analytical solution used in this simulation. This nonFickian effect
impacts benchmarking with numerical solutions in this case because the
full analytical spatialconcentration volume is affected. Future attempts
are expected that revise the advectiondispersion equation with a non
80
Fickian term. This should not, however, detract from the analytical
equation's current utility in generally accurate farfield calculations.
In conclusion, choices in MOC and SUTRA run options in this case
can have a significant effect on the results; in many cases the magnitude
of change is as great as change caused by input aquifer parameter
ranges. Both MOC and SUTRA models with source flows greater than
the field average, are conservative in the sense that total spatial
concentration volumes are greater than the spatialconcentration volume
from the analytical solution and the leading edge concentrations extend
beyond the analytical solution's. When the numerical code source
volume flows are equal to the analytical solution average flow, the
advective centers are lower; but the outer length contours are generally
equal to, or extended further than, the analytical solution. The analytical
solution equation appears to be underestimating the spatial
concentration volume near the source in this case, due to nearfield non
Fickian behavior.
81
Concentration ppm Concentration ppm Concentration ppm
Figure 519 SubCases 1A, 1B, 1C Centerline
Concentration Profile Comparison
Case 1A Centerline Concentration Profile Comparison
Range Distance 1000's Feet
Case 1B Centerline Concentration Profile Comparison
Range Distance 1000's Feet
Case 1C Centerline Concentration Profile Comparison
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
Range Distance 1000's Feet
82
Concentration ppm Concentration ppm
Figure 520 SubCases 1D, 1E Centerline
Concentration Profile Comparison
Case 1D Centerline Concentration Profile Comparison
1000 T i 1 i1 ! M i j ! r 1
, i ! ! i j riM
! 4 i : ( 1 l  i i 1 SUTRA r ANALYTIC .
1 ! 1/ ] U i %, i  l\l\ \L _i 1 1 i
900 // i \! b T\! ! n
\ v "\ 
LlL \ i\ \ K \l v 1 i 
onn J ML w j
I  LL i i 1
1UU n ^ 1 I \ _L jl i
1 2 3 4 S 6 7 B 9 10 11 12 13 14 15 16 17 18 19 20 21 22
Range Distance 1000's Feet
Case 1E Centerline Concentration Profile Comparison
Range Distance 1000's Feet
83
MOC
SUTRA
00
Uio4'
LTIO4'
I.1M04*
1104*
0000
sao>
f
simuciET
M0CC1E
Analytical (TACT)
on
If
c 3
as
<0
= 6
f.8
=?
cd m
CO w
Â§ 3
m
CO'S
O
a SL
2 
CD
< 52
o C
c H
3 5
CD >
Concentration ppm
Distance 1000 Feet / Unit
O
o
5.2 Case Study/Benchmark 2
A field case of hexavalent chromium contaminant transport from
Liberty Aircraft Plant in Long Island, New York has been documented by
Perlmutter and Lieber (1970). The hexavalent chromium is waste from
anodizing aluminum by electrolytic action in chromic acid. Plating
wastes were deposited into pits in contact with the saturated zone (on
seasonal average). Hexavalent chromium was first discovered in a well
right next to the disposal pits in 1942. The plume was modeled for the
year 1949, just before contact with the Massapequa Creek caused non
uniform flow. In 1949, the plume extended about 850 feet wide, and
3,850 feet end to end or 3,675 feet from the edge of the source length
wise. The maximum concentration reading was 40 mg/L. Estimates of
the plume position from field measurements in a plan view are shown in
Figure 522, and in crosssections on Figure 523. The field geologic
section is shown in Figure 524. The soil was Upper Pleistocene glacial
deposits of permeable fine sands with lenses of clay, and silt. This media
exhibited considerable heterogeneity, and local permeability
discontinuities affecting dispersion and advection by branching and
agglomeration can be seen in the cross sections; however, the larger
scale gross state is represented as homogeneous.
Monitoring wells were installed and samples collected at five foot
depth intervals. The upper unconfined saturated thickness ranges from
60 to 140 feet. Chlorides, nitrates, and sulfates were also found in higher
concentrations than uncontaminated water on the island. This is a
relatively good case for testing ground water transport models except that
exact source concentrations are not known. Records indicate that
approximately 52 pounds of chromium in 200,000 to 300,000 gallons of
water per day were disposed, on an average in the earliest years, but
concentrations may have been higher and the actual amount was
probably not steadystate. The aquifer was found to be primarily
85
86
Figure 523 Case 2 CrossSection Plume Location
Estimates from Field Measurements Liberty Aircraft Plant
in Long Island, New York; Slices D, E, & F (Perlmutter and Lieber,
EXPLANATION
Contaminated water
26
r
Teat well and field number
5
Line of equal concentration of contaminant, 194962,
in milligrams per liter
Dashed ttAfn approximate
0
SOO
i
1000 FEET
87
Figure 524 Case 2 Field Geologic Section Liberty
Aircraft Plant in Long Island, New York; Slice A (Perlmutter and
Lieber, 1970).
SOUTH
*4
Teat well end field number
Munfer at bottom of well if head at ocrt+n drpth. *r
fetl above auan tea Ur*l
Approximate direction of groundwater flow
O f>
Gravel
Sand
Sand, allt, and
aandj elaj
I
Clay
Contaminated water
88
medium to coarse sand with 35% porosity and a hydraulic conductivity of
1.0 to 3.0 x 10'3 ft/sec. The vertical conductivity was estimated to be five
to ten times larger.
In 1946, 0.05 mg/L was the recommended limit for hexavalent
chromium. A median value for chromium recently in normal ground
water at another North American location is 0.00043 mg/L (Matthess,
1982). Accuracy of laboratory testing at that time did not allow qualifying
chromium below 0.03 mg/L. Hexavalent chromium is toxic while trivalent
chromium may have little toxicity in small doses. Oxides of chromium are
neurotoxins. The chronic and subchronic reference doses for hexavalent
chromium are 5. x 10'3 mg/Kg/Day and 2. x 102 mg/Kg/Day respectively.
The chronic and subchronic reference dose numbers for trivalent
chromium are 1. and 10. mg/Kg/Day. These values are based upon rat
studies (ATSDR, 1990). If an extract from waste contains more than 5.0
mg/liter chromium independent of valence, it is considered a hazardous
waste (CFR Title 40,1990). The current SDWA level is at 0.1 mg/L
(Federal Register 56 (20), 1991). The chromium soil levels above 400
mg/liter, regardless of valence, are classified as a hazardous waste
(Federal Register 55 (145), 1990). This is interesting, since naturally
occurring chromium can sometimes exceed this level (Matthess, 1982).
The most common naturally occurring chromium is as chromite (FeCr204)
ore (Nerbergall et al., 1972). However, pollution can usually be
distinguished by chromium being present in the form of yellow chromate
(CrC>42) at pH ~ 6, and orangered dichromate (C^Ot2) at pH ~ 2 to 6
(Chambers et al., 1991). Regulatory concern is for the possibility of
microbial oxidation of trivalent chromium back to the much more toxic
hexavalent form (Chambers et al., 1991, Hanson et al., 1993).
Hexavalent compounds do not precipitate in the pH range from 1.0 to 9.0,
so stay very mobile until sorption occurs below pH 5 to 6 (Chambers et
al., 1991).
Sorption of chromium was studied by Ku et al. (1978).
Polynuclear complexes form in the presence of more than one metal and
89
pH determines speciation. Hydroxyl complexes dominate at high pH's
and the mobile ion at lower pH's. Ku removed free oxide crystalline
forms using DCB (dithionitecitratebicarbonate). Amorphous forms of
chromium were stripped by an ammonium oxalate. Both filtrates were
acidified with HNO3 to mobilize the metals. Extracts were then analyzed
by direct aspiration on an atomicabsorption spectrophotometer. The
maximum sorption capacity was found to be 19 milligrams per kilogram
of soil. Ku postulated that the iron oxides and the amorphous iron
oxyhydroxides provide a matrix upon which heavy metals may be
attached sorbed, coprecipitated, or occluded. The chromium will also
leach from the soil materials after terminating the source depending upon
Eh, pH, relative ion saturation phases, anaerobic conditions, and
fluctuating water table. Based on this sorption information, the
retardation coefficient was calculated to be a maximum of 2.28 for this
thesis study. Allen et al. (1993) indicate general agreement with the
sorption process and the amount of sorption originally described by Ku.
However, using a single constant retardation coefficient is a simplification
of a much more complicated process. The initial phase may be high
decay and adsorption due to the presence of oxygen and later decreases
in a flow progression to a suboxic regime. This could be thought of as a
twocomponent transport with chained decay reactions. As the oxidation
progresses from the greatest free energy through the redox zonations,
transport of redoxsensitive chromium ions are affected. Hexavalent
chromium is reduced at low pH's to trivalent chromium by Fe(ll) on the
surfaces of minerals like biotite and magnetite (Davis et al., 1993).
Hexavalent chromium may also be reduced by organic acids or bacteria
in aerobic and anaerobic environments (Nyer, 1985). Dilute aqueous
forms of hexavalent chromium are in chromate (CrQ*2') and bichromate
(HCrCV) anions. Trivalent chromium forms various hydroxide complexes
depending on pH (Domenico & Schwartz, 1990); Cry= Cr3* + Cr(OH)2+ +
Cr(OH)2+ + Cr(OH)3 as a mole balance. (Pourbaix (1973) has
developed a chromium EnpH redox potential chart.) The sorption of
90
chromate and bichromate ions on the surfaces of oxide minerals like
hydrous iron and aluminum oxides is a function of pH range with different
ion complexing and dissociation. The sorption of hexavalent chromium
is also highly dependent upon competitive sorption by ions such as
sulfate, bicarbonate, and silicate.
5.2.1 Problem Definition for Case Study/Benchmark 2
The site model grid plan is shown in Figure 525 for SUTRA and
Figure 526 for MOC. Data for model input is composed from already
cited references plus others (Pinder, 1973; Cleary, 1978; Sawyer, 1978,
Wagner et al., 1982; Cherry et al., 1984; Wood et al., 1984; Reily et al.,
1987; Beljin, 1988, and Yong et al., 1992). The input data is presented in
Tables 54, 55, and 56. In SUTRA, constant source concentration and
no damping is used. In MOC, the SIP method is used.
91
Figure 525 Case 2 Site Model Grid for SUTRA Liberty
Aircraft Plant in Long Island, New York.
fc x
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
5? p F F TO sSS F 1 m a 3 F F
P In i i in
A*
60
77
94
111
128
145
162
170
196
213
230
9A7
264
261
298
w;
332
Wi AX sXA VA XV AV HE 22S XX! .w
(Node Numbers Shown) 357
D.Long 93
92
Figure 526 Case 2 Site Model Grid for MOC Liberty
Aircraft Plant in Long Island, New York.
D.Long 93
Constant Head
& Leakance
Source Concentration
Constant Head
& Leakance
Periphery = Zero
Transmissivity &
Zero Head
Constant Head
& Leakance
(Node Numbers Shown)
93
Table 54 Case Study/Benchmark 2 Model Input Liberty
Aircraft Plant in Long Island, New York.___________________
Parameter Type Parameter Value Symbol
Grid Element Dimension A & C: 225 ft x 225 ft B: 200 ft x 200 ft w or Axgrid
Source Width A: 225 ft B & C: 200 ft Y
Source Concentration 73.00 ppm Co
Source Volume 0.200 ft3/sec
Source Time Duration 7.671 yrs 2.42078 x 108 t =7.671x365.25x24 x60x60
Porosity 0.35 n
Aquifer Thickness 110ft m
Hydraulic Conductivity 2.0 x 10*3 ft/sec K
Retardation Factor A & B: 2.28 C: 2.00 Rf
Sorption Distribution Coefficient A &B: 4.170x10*3 C: 3.257 x 103 K,, [n/(1 n)p J x {R,1}
Soil Bulk Density Soil Solid Density 107.4 lb/ft3 165.3 lb/ft3 Bd Bs
Decay Half Life Large tl/2
Dispersion Coefficient Longitudinal 70 ft a,
Dispersion Coefficient Transverse 0.17x70 = 14.0 ft ay
Recharge Besides Constant Head Nodes 0
Average Hydraulic Gradient 0.00455 i = delta H/delta L = 60.042.6/3825
Average Darcy's Velocity 9.10 x 10'6 ft/sec q = K i
Average Darcy's Volumetric Flow Rate per Cell A & C: 0.225 ft3/sec B : 0.200 ft3/sec Q = m w Ki
Average Water Velocity 2.60 x 10*5 ft/sec vw=q/n
Contaminant Velocity 3.99 x 10'6 ft/sec v = vw / Rf
D. Long 93
94
Table 55 Case Study/Benchmark 2 Model SubCase Input
Values Liberty Aircraft Plant in Long Island, New York.___
Case 2A 2B 2C
Source ft3/sec 0.225 0.200 0.225
Grid Size ft 225 200 225
Retardation Factor 2.28 2.28 2.00
D. Long 93
Table 56 Case Study/Benchmark 2 Model SubCase
Calibration and Sensitivity Trial Values Liberty Aircraft Plant in
Long Island. New York._________________________________
Parameter Type Values
Element Grid Dimension ft 165, 175, 200, 225
Source Width 165, 175, 200, 225
Retardation Factor 1.74, 2.0, 2.07, 2.28
Source Volume ft3/sec 0.20, 0.225, 0.36, 0.46
Soil Bulk Density 132.2, 107.4
Dispersivities:
Longitudinal 50, 70
Transverse 7.5,14.0
Gradient 0.00837, 0.007096, 0.00680, 0.00455
Source Type Definition Injection, Boundary
D. Long 93
95
5.2.2 Case Study/Benchmark 2 Results, Comparison, and
Conclusions
The source in the numerical models can be specified as a part of
the boundary or offset as an injection well. The effect of input source type
choice is plotted in Figures 527 and 528. It can be seen that the profile
results are significantly affected. The gridboundary source was selected
for use in Case 2. The total source volumetric flow rate in the source cell
is equal to the average analytical control area flow in all subcases (see
Tables 55). In MOC, the source volume does not need to be specified
when it is part of the boundary. The parameters were varied for
calibration and sensitivity trials between probable values in Table 56.
As a result of sensitivity trials in this study, it was found that the choice of
model source type can influence the outcome by the same magnitude as
the possible variability in aquifer values. At one time, MOC had some
inaccuracy with the flow field and mass solute balance near radial flow
sinks (ElKadi, 1988) as would be evident with an offset source. This
same situation is observed and discussed in Case 1. A similar problem
is occurring at an injection source. A transient storage improvement was
made to MOC (Goode, 1990), but it is not known whether all storage
difficulties were solved. It is suspected that MOC mass volumes are not
accurately handled in a local radially divergent velocity vector field. (The
effect of increased grid density around these areas was not investigated
in this study.)
The centerline concentrations of all three solutions for SubCases
2A, 2B, and 2C are plotted in Figures 529 and 530. The same general
trend is observed for the numerical solutions in SubCases 2A, 2B, and
2C as in some subcases of Case 1. As before, MOC is found to be
somewhat high in the advective center region, but MOC's outer contour
comes close to the analytical's. SUTRA's centerline matches reasonably
well except that the advective center region is low and the leading edge
is protracted and never terminates. The SUTRA SubCase 2C solution
96
has a similar centerline profile until the leading edge region, where it
also has difficulty by extending too far. The length of the analytical plume
was found to decrease with shrinking source width.
The threedimensional plots and longitudinalvertical three
dimensional slice profile plots for the three models of SubCases 2A and
2B are presented in Figures 531,532, 534, and 535. A crosssection
of SubCase 2A is plotted in Figure 533. The contour plots for the Sub
Cases 2A and 2B are plotted in Figures 536 through 541. In general for
these subcases, all of the codes predict the total spatialconcentration
volume reasonably well. A comparison of the results is shown in Table
57. Since the actual spatialconcentration volume is misshapen due to
aquifer heterogeneity and unsteady source, the numerical code volumes
are compared to the analytical solution volumes. MOC's SubCase 2A
total spatialconcentration volume is only 3.6% less than the analytic
solution model. MOC's SubCase 2A spatialconcentration volume
shape is surprisingly close to the analytical solution. Its shape varies the
most from the analytical shape by a higher reading along the centerline.
SUTRA's SubCase 2A spatialconcentration volume is 4.1% lower in
comparison to the analytic solution. In SubCase 2B, MOC's spatial
concentration volume is less than the analytical by only 5.5% and
SUTRA is lessor by 7.7%. The analytical solution length is the shortest
and widest of the three models (see Table 58). MOC is closer than
SUTRA to the analytical outer contours and volumes in all Case 1 sub
cases. In all subcases, MOC's width and length is also closest to the
field data. In contrast, SUTRA has a narrow width for all subcases, and
the longest length of the three models. While the concentration volume
shape details of all three models are a little different, the total volumes
match well.
Since no damping is used in SUTRA for Case 2, numerical
dispersion is assumed to be the sole cause of the low concentration
longitudinal extension. An attempt to remove SUTRA's leading edge
numerical dispersion with a superposition subtraction was not
97
