Ground water solute transport modeling

Material Information

Ground water solute transport modeling case study benchmarks and validation methodology
Long, David Scott
Publication Date:
Physical Description:
viii, 167 leaves : illustrations ; 29 cm


Subjects / Keywords:
Groundwater -- Pollution ( lcsh )
Groundwater flow -- Mathematical models ( lcsh )
Groundwater flow -- Mathematical models ( fast )
Groundwater -- Pollution ( fast )
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references (leaves 128-134).
General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Civil Engineering.
General Note:
Department of Civil Engineering
Statement of Responsibility:
by David Scott Long.

Record Information

Source Institution:
University of Colorado Denver
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
30656454 ( OCLC )
LD1190.E53 1993m .L66 ( lcc )

Full Text
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 by David Scott Long
Ail rights reserved.
MATHCAD is a registered trademark of Mathsoft, Incorporated
MS-DOS 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
Edward B. Nuhfer

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
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 steady-state flow field. There is good
general agreement between models in the second case, but not in the
first case with a grid-interior source. Each of the two numerical codes
studied have a tendency for specific shape bias, which is evident in all
sub-cases. Sub-cases 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' spatial-concentration volumes are comparable to

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 far-field, theoretical
difficulties with the short-duration or near-field spatial-concentration
volume are evident due to the non-Fickian 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.

Professional Technical Reviewers:,
Brian J. Long
Duane E. Long
Roger S. Stillwater
Manuscript Reviewers:
Hannah L. Kelminson
Marlene V. Long

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

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
5.1.3 Case Study/Benchmark 1 MOC Model Results and
5.1.4 Case Study/Benchmark 1 SUTRA Model Results and
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
6. Summary and Conclusions......................................119
6.1 Summary....................................................119
6.2 Conclusions................................................120
6.3 Conclusion and Interpretation Summary......................122
Appendix A. Ground Water Modeling Term Definitions..............135
Appendix B. Sample Output Data for Case 2A......................139

1. Introduction
Ground water solute contaminant transport modeling plays an
important role in investigations and predictionsincluding site
characterizations, risk-based exposure assessments, compliance
reporting, remedial engineering, planning trade-offs, 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 multi-disciplinary 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 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. Man-made 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 F-H 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,

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 clean-up 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 pre-RCRA 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.

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; El-Kadi, 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,
under-documented 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 1-1.
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).

Table 1-1 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 = Cross-sectional area i = Gradient q = Specific discharge velocity v = Velocity L = Load K = Stiffness E = Modulus of Elasticity U = Displacement

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

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 near-field 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 set-up 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

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 non-existent.
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 time-consuming
Boundary conditions are difficult to understand and apply
Inappropriate code applied to problem
Models contain too many assumptions
Non-modelers often misunderstand the results
Models exhibit a high degree of uncertainty that is difficult to
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
End-user client expectations are often too high
Modeling objectives are poorly defined
Non-uniqueness 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 inter-comparison with
other modeling programs. Usually an order-of-magnitude result is
sufficient for professional judgments regarding concentration or travel
time, especially with data and degradation process uncertainty (El-Kadi,

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

Figure 1-1 The Computer Code Model Validation
D. Long S3

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 three-dimensional areal
steady-state 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 state-of-the-art 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

are compared to each other and to an analytical solution. Both cases are
represented as a homogeneous steady-state flow field with a continuous
rectangular source. The effect of changing different sensitive parameters
and run options on a grid-interior source are presented in Case 1. Case
2 is an actual case history, with data from the literature, compared by
parameter calibration. Sub-cases vary run options, grid size, dispersivity,
source volume, and source type. Comparison plots are made of the
spatial-concentration 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
2. The computer code is available for use on a MS-DOS based
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.

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
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. Laplace-Transform-Galerkin (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

deterministic, closed-form solutions of the governing differential
Probabilistic and stochastic methods are considered
indeterministic. These include so-called 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 sub-process (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
Particle tracking and MOC are popular transport methods that
diminish numerical dispersion. There is also a new approach, the
Laplace-Transform-Galerkin (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 steady-state
velocity fields (Allen et al., 1993).

2.1 Ground Water Flow Governing Equations and Model Type
Ground water models can be spatially categorized as follows:
Two-dimensional areal
Two-dimensional profile
Quasi three-dimensional
The time-domain run states are classified as follows:
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
Two-dimensional 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:
L = Kz

+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
-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
Three-dimensional and two-dimensional profile
Three-dimensional distribution of:
- Head
- Hydraulic conductivity
- Storage coefficients
Governing equation form:

+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 so-called "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.
Two-dimensional models are used for aquifers confined,
confined-leaky, unconfined, and mixed. Head distribution and the
aquifer storage coefficient is required. In quasi three-dimensional
models and two-dimensional models of confined-leaky 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

the difference between the aquifer and confining beds, a transient nature
in the water flux, or a deep aquifer compared to the two-dimensional flow,
a full three-dimensional model may be preferred.
Full three-dimensional 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 two-dimensional 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 out-of-plane
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.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 two-dimensional 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 + Dy-r-y Vx -
dx2 dy dx dt
With linear sorption and decay, the equation may be expressed as (Bear,
a2c a2c
Vxttx + VxCIy-
ac ac
- Vx-A.C =
Rfdx dt
Adding in the source and sink term, the equation is (based on Goode &
Konikow, 1989; Goode, 1990):
ac a2c
dt VxCU Ridx2
i rc-ciw

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
2-1. If fluid density changes during the simulation time, then the flow and
transport equations must be solved in conjunction with one another.

Figure 2-1 System-Level Flow Chart for Ground
Water Flow and Transport

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 2-2. 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 set-up 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
triai-and-error 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.) Post-audits are done later
with new field data to see if the initial model predictions were correct.

Figure 2-2 The Site Modeling Process
O.Long S3

2.4 Establishing Ground Water Model Quality Assurance and
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
- 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

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

A study done for the EPA evaluated the four following
saturated/unsaturated models:
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 two-dimensional homogeneous soil
section problem:
From this comparison it was concluded that VS2DT was the most efficient
and accurate (Carovillano, 1993), although the other model results were
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

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. Sixty-one 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:

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
AQUIFEM-1 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 Density-Temp, 3-D FEM-IFDM 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 2-3. This list is only intended to
be representative of some good codes, not entirely inclusive.

Figure 2-3 Some Premier Ground Water Modeling
Codes Related to Solute Contaminant Transport

Flow Codes:

Transport Codes:

Unsaturated Tran snort Codes:

Special Purpose Codes:
HELP Landfill Seepage AQTESOLV Well Test Curve Matching PRZM2 Pesticide Root Zone ModelCad Preprocessor SURFER-Plotter PREMOD Preprocessor TOUGH Multiphase Fracture Flow

Geochemical Codes:

Parameter Estimation Codes:


MODFLOW is a popular three-dimensional finite difference flow
program with a "strongly implicit procedure" or "slice successive over-
relaxation solution" technique. PLASM is a two-dimensional finite
difference model using an interactive alternating direction implicit
procedure. AQUIFEM-1 is a two-dimensional finite element linear
triangle model using Crout's direct solution technique. AQUIFEM-N is a
three-dimensional version of AQUIFEM-1. MODPATH is a three
dimensional semi-analytic 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 steady-state problems with particle
tracking. MOC is a mature two-dimensional 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
steady-state 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 3-D 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

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, user-friendly 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. Finite-element flow and transport models

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 time-and-distance shape is based on a normal
distribution. The governing differential equation for dispersion by
conservation of mass assumes:

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 one-dimensional form, rearranging for a semi-infinite medium and
source at x= 0, the dispersion equation 3.2 becomes (Freeze & Cherry,
ac _a2c ac
aT=Da?-vaT (3-3]
A closed form "exact" solution for a one-dimensional advection-
dispersion flow source was derived by Ogata & Banks (1961):
The error-function is defined from:

The experimental Break-Through-Curve (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
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 advection-dispersion rectangular constant source 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
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 far-field. 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
3-1. As the transverse and vertical dispersivities approach zero, the
Ogata-Banks 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,

x' = (H Z)2
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) =

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 Three-Dimensional 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
Three-Dimensional 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 one-dimensional equation 3.4 was coded in
TACT and was compared against one-dimensional BTC's.
The analytical equation error functions were checked
independently against published values from equation 3.5. These
results are graphed in Figure 3-2. The TACT code three-dimensional
results were independently verified against the governing equations
programmed into MATHCAD, which calculated identical results.

Figure 3-1 Analytic Solution Transport Coordinate
O. Long S3
Figure 3-2 Error Function & Complementary Error
Function Check Plot

4.0 Solute Transport Numerical Solution Modeling
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.

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 two-dimensional 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 time-dependent transport solutions. SUTRA-PLOT
provides for plotting the output.
The FE method is employed for the fluid and solute mass balance
equations. All remaining non-fluid 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 two-dimensional 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 third-dimension. 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

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 semi-analytical 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 velocity-dependent dispersion for
anisotropic aquifers.
The finite element mesh is composed of quadrilateral elements.
Each of the four comers has a node representing the mid-point of a z-
edge (see Figure 4-1). In SUTRA, constant parameter values are
assigned to each mesh node, either element-wise, node-wise or cell-
wise. Cells are considered to be centered on a node. Cell boundaries
are half way between element sides (Figure 4-2). The spatial
distributions of parameters are assigned as follows:
- K element-wise
- h node-wise
- S cell-wise

f:'9ure *
ent ($*'fjmt
SS' ,9^6 e'ein, Me

Figure 4-2 Definition of Two-Dimensional Quadrilateral
Elements, Cells, and Nodes (Voss, 1984)

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

field studies have indicated that density-induced 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.

4.2 MOC Model Description
MOC (Konikow & Bredehoeft, 1978) is a two-dimensional finite
difference flow model with solute transport in a homogeneous saturated
aquifer. MOC uses a block-centered (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 4-3, and MOC's algorithm is shown in Figure
MOC contains the following chemical reactions in a single
advection-dispersion solute transport equation:
First-order, irreversible rate-reaction (radioactive decay)
Reversible equilibrium reactions
- Linear sorption
- Freundlich sorption
- Langmuir sorption
- Monovalent ion exchange
- Divalent ion exchange
- Monovalent-divalent ion exchange
- Divalent-monovalent ion exchange

Figure 4-3 Definition of MOC Particle Movement in
Relation to Grid & Node (Konikow & Bredehoeft, 1978)
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.
*',/-1 1.7-1
j-i.i /,/ W- vM.
Vi,/.I J^/*1 T */1. /!

Node of finite-difference 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 finite-difference 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 one-half of the
area of a cell.
' Node of finite-difference (pld
Location of particle
Halation between possible Ini-
tial locations of points and Indlcas of ad-
jacent nodes.

Figure 4-4 Main Level MOC Algorithm
(Konikow & Bredehoeft, 1978)

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 "flow-equation-removed" form of the transport
equation and is considered to be "fluid-mass-conservative" (Voss 1984).
This is stated as follows (based on Konikow & Grove, 1984; Goode &
Konikow, 1989; Goode, 1990):
Rfirnie dxi
_3h 3m
S + W-nc
dt dt
Dij (Xijmn
mn = (mne)o + S[h ho]
Tij ah
nine dxj

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.

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 (Ax-grid 2 x Dispersivity). Centered finite difference
schemes, considered not mass conservative, have a truncation error on
the order of (Ax-grid)3. These are subject to upstream and downstream
oscillations unless the grid size is small (Ax-grid < 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 (Ax-grid). The element aspect ratio should be
no greater than 1.5.
A block-centered 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. Mesh-centered finite difference elements and
finite elements have both flow and specified heads at the nodes located
at element grid intersections.

4.3.2 Boundary Conditions, Sources, and Sinks
Boundary condition assignment is the most delicate portion of
modeling. The flow pattern of steady-state 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 no-flow
- Neumann conditions
- Function of flux derivative
Head-dependent 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 steady-state situation, the heads are set at
the water table elevation. An example of a site model is shown in Figure
Specified flow boundaries are simulated in the model by using
injection or extraction wells. Use of flux boundaries exclusively should

Figure 4-5 Example of Model Boundary Conditions
(Konikow, 1977)
0 1 2 MILES
1 ----1--1I------1
Zero-transmissivity call
Constant-head cell
A Disposal-pond call.
I Irrlgation-racharga eall
L CanaMaakage cell
----- Boundary or areas in which alluvium is absent or unsatureted

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 steady-state problems. Transient solutions have
initial heads, so flux boundaries could be justified. No-flow flux
boundaries are placed at an impermeable bedrock, a ground water
divide, a fault zone, or a streamline. Some models automatically assume
no-flow boundaries but other models must have them defined with zero
Head-dependent 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 steady-state 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 two-order-of-magnitude. 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 one-dimensional model could be used
as input to a saturated model. Coupled models could account for

movement of the water table with recharge. Most field-scale 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 two-dimensional areal models and quasi three-dimensional 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
steady-state 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.

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 sub-model 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 steady-state initial conditions
- Head varies spatially and temporally
- Flow in = flow out
- Boundary conditions can cause influence problems
Steady-state initial head conditions
- Head constant
- Assumes no flow
- Good for datum superposition in draw-down 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(Ax-grid)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.

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 log-normal 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 steady-state 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

for so-called "inverse" methods of calibration (Ghohestani-Bojd, 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
non-parametric 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 non-unique and derived conductivities are sensitive to observed head
noise and errors.

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 (time-dependent)
solute transport problem in a steady-state saturated flow field with
dispersion, sorption, and decay from a continuous rectangular source.
The flow field for the two cases are defined as steady-state two-
dimensional or quasi three-dimensional. 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:
- Three-Dimensional Analytical Solution Contaminant Solute
Transport Code for Saturated Matrix with Dispersion and
- Version 1.1
- October 1993
- FORTRAN (M.S. 4.1)
- Saturated/Unsaturated Transport
- Version V06902D
- June 1990
- Two-Dimensional Solute Transport Method of Characteristics
- Version IGWMC-FOS 07 PC Version 3.0
- November 1989
- FORTRAN (M.S. 4.1)

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 5-1 for SUTRA and
Figure 5-2 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 grid-interior
source concentration was applied at an equivalent element width of
1,000 feet at node 264, in SUTRA as shown in Figure 5-1, and in MOC as
shown in Figure 5-2. The input data is summarized in Tables 5-1 and
5-2. Each sub-case varies one or more of the parameters relative to the
other sub-cases. Sub-Case 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 sub-cases:
Hydraulic conductivity: K = 5.0 x 10-4 feet/second
Porosity: n = 0.20
Retardation constant: Rf = 2.00

Hydraulic gradient:
i = AH
= 200-50 =0.0075
Darcy's velocity:
q= Ki = 5.0x1 O'4 (0.0075) = 3.750x10-6 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 10-6 = 1.875 x 10'5 feet/second
n 0.20
Contaminant velocity:
v = vffi = 1.875 x 10-5 = 9.375 x 106 feet/second
R, 2.0
Dispersion coefficients: a* = 200 feet
a, =50 feet
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
Kd =

centerline concentration is too low near the source. The other high
volumetric flow sub-cases have high centerline concentrations. For a
large source, MOC's centerline concentrations for decay (Sub-Case 1D)
from the advective center to the leading edge are known to be high
because they exceed the analytical solution one-dimensional results in
this region. See Figure 5-13. 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.

5.1.4 Case Study/Benchmark 1 SUTRA Model Results and
Several different execution option combinations are available for
SUTRA. The different SUTRA execution options that were tried are as
listed in Table 5-3.
Table 5-3 SUTRA Execution Options
SUTRA Run Option Types Constant Concen- tration Source Flow Damping Steady-state 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
side-lobe and upstream concentration oscillations. Steady-state
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 three-dimensional and longitudinal-vertical slice
profile solution plots for Sub-Case 1A and 1E are shown in Figures 5-14
through 5-17. The SUTRA model results for centerline concentrations
comparing all sub-cases are presented in Figure 5-18. In Sub-Case 1E
SUTRA is showing a tendency to extend its leading edge.

Figure 5-3 Case 1A Analytic Three-Dimensional Plot
Figure 5-4 Sub-Case 1A Analytic Three-Dimensional
Longitudinal Slice Profile Plot Longitudinal-vertical slices
from the preceding three-dimensional plot.
Range Distance 1000's Feet

Concentration ppm Concentration ppm
Figure 5-5 Case 1 Analytic Centerline Profile
Concentration Sub-case Comparison Plot
Figure 5-6 Case 1 Analytic Comparison of One
Dimensional and Three-Dimensional Profiles

Figure 5-7 Sub-Case 1A MOC Three-Dimensional Plot
C Q.
Figure 5-8 Sub-Case 1A MOC Three-Dimensional
Longitudinal Slice Profile Plot Longitudinal-vertical slices
from the preceding three-dimensional plot.
Range Distance 1000's Feet

Figure 5-9 Sub-Case 1E MOC Three-Dimensional Plot
1000's Feet
Figure 5-10 Sub-Case 1E MOC Three-Dimensional
Longitudinal Slice Profile Plot Longitudinal-vertical slices
from the preceding three-dimensional plot.
Range Distance 1000's Feet

Figure 5-11 Case 1 MOC Centerline Profile
Concentration Sub-Case Comparison Plot

Concentration ppm Concentration ppm
Figure 5-12 Case 1 MOC Source Volume Sensitivity
Centerline Profile
Figure 5-13 Case 1C & 1D MOC Comparison to
Analytic One-Dimensional Centerline Profiles
Case 1C MOC

* e Case 1D MOC Case 1C Analytic 1 -D
i,\ N
\ i- o Case 1D Analytic 1 -D .
\ <
L\ 1
-*i H C i 8 i O
5 6 7 8 9 10 11 12 13 14 15
Range Distance 1000's Feet

Concentration ppm
Figure 5-14 Case 1A SUTRA Three-Dimensional Plot
C o.
0) &
Figure 5-15 Sub-Case 1A SUTRA Three-Dimensional
Longitudinal Slice Profile Plot Longitudinal-vertical slices
from the preceding three-dimensional plot.
Range Distance 1000's Feet

Figure 5-16 Sub-Case 1E SUTRA Three-Dimensional
Figure 5-17 Sub-Case 1E SUTRA Three-Dimensional
Longitudinal Slice Profile Plot Longitudinal-vertical slices
from the preceding three-dimensional plot.
Range Distance 1000's Feet

Figure 5-18 Case 1 SUTRA Centerline Profile
Concentration Sub-Case Comparison Plot
Range Distance 1000's Feet

5.1.5 Case Study/Benchmark 1 Comparison and Conclusions
A comparison of centerline concentration profiles for Sub-Cases
1A, 1B, and 1C are plotted in Figure 5-19. Comparisons of centerline
profiles for Sub-Cases 1D and 1E are plotted in Figure 5-20. In Sub-
Case 1A, SUTRA has good agreement with the analytic solution while
MOC has a higher centerline concentration. Sub-Case 1 B's source
volume is half of Sub-Case 1A. Sub-Case 1C, with higher dispersivity
(and high source volume), has the most pronounced concentration
deviation of the numerical solutions from the analytical solution (see
Figure 5-19, Sub-Case 1C). Generally a finer grid relative to dispersivity
is purported to mathematically give better numerical results. Sub-Case
1D dispersion spreading length increases only by the square root of the
dispersivity. Sub-Case 1E numerical models have a source volumetric
flux equal to the analytic solution average cell flow.
For Sub-Cases 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 sub-cases 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 cut-off could be
selected at some low percentage (<10%). MOC's centerline
concentration in Sub-Cases 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 sub-Case 1A, that although the shapes of SUTRA and
MOC are different, the total spatial-concentration volume integration
(Sub-Case 1 A) differ by only 5%. However, their volumes deviate by
26.3% in Sub-Case 1E. They are both significantly larger than the
analytical solution spatial-concentration volume for Sub-Case 1A, as
would be expected with higher source volumes. In Sub-Case 1E, MOC's

spatial-concentration volume is 8.1% higher, but SUTRA is 36.6% larger
than the analytical solution spatial-concentration 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
one-half of the grid length when damping is used, especially in highly
advective velocity fields. The concentration contour plots for Sub-Case
1E are shown in Figure 5-21. 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 sub-cases. 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 three-dimensional 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 far-field where the length to aquifer
thickness is large.
While the analytical solution should be considered the most
accurate in the far-field, it should not be construed as complete "truth"
near the source. In development of the dispersion term in the analytical
advection-dispersion equation, the coefficient of mechanical dispersion

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 time-concentration
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 slow-velocity
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 near-field. This is because
dispersion mixing is confined to two dimensions far from the source, but
is in three-dimensional expansion near the source. Therefore, the
Eulerian velocity field is only valid for far-field 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 near-field
and early transport times. Analytical equations differentiating near-field
and far-field 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 non-Fickian effect
impacts benchmarking with numerical solutions in this case because the
full analytical spatial-concentration volume is affected. Future attempts
are expected that revise the advection-dispersion equation with a non-

Fickian term. This should not, however, detract from the analytical
equation's current utility in generally accurate far-field 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 spatial-concentration 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 near-field non-
Fickian behavior.

Concentration ppm Concentration ppm Concentration ppm
Figure 5-19 Sub-Cases 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

Concentration ppm Concentration ppm
Figure 5-20 Sub-Cases 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

Analytical (TACT)
c 3
= 6
=? cd m
CO w
§ 3
a SL
2 -
< 52
o C
c H
3 5
CD >
Concentration ppm
Distance 1000 Feet / Unit

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 5-22, and in cross-sections on Figure 5-23. The field geologic
section is shown in Figure 5-24. 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 steady-state. The aquifer was found to be primarily


Figure 5-23 Case 2 Cross-Section Plume Location
Estimates from Field Measurements Liberty Aircraft Plant
in Long Island, New York; Slices D, E, & F (Perlmutter and Lieber,
Contaminated water
Teat well and field number
Line of equal concentration of contaminant, 1949-62,
in milligrams per liter
Dashed ttAfn approximate
1000 FEET

Figure 5-24 Case 2 Field Geologic Section Liberty
Aircraft Plant in Long Island, New York; Slice A (Perlmutter and
Lieber, 1970).
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 ground-water flow
O f>
Sand, allt, and
aandj elaj
Contaminated water

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 orange-red 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

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 (dithionite-citrate-bicarbonate). 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 atomic-absorption 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
two-component transport with chained decay reactions. As the oxidation
progresses from the greatest free energy through the redox zonations,
transport of redox-sensitive 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 En-pH redox potential chart.) The sorption of

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 5-25 for SUTRA and
Figure 5-26 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 5-4, 5-5, and 5-6. In SUTRA, constant source concentration and
no damping is used. In MOC, the SIP method is used.

Figure 5-25 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

Wi AX sXA VA XV AV HE 22S XX! .w

(Node Numbers Shown) 357
D.Long 93

Figure 5-26 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)

Table 5-4 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 Ax-grid
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.0-42.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

Table 5-5 Case Study/Benchmark 2 Model Sub-Case 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 5-6 Case Study/Benchmark 2 Model Sub-Case
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
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

5.2.2 Case Study/Benchmark 2 Results, Comparison, and
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 5-27 and 5-28. It can be seen that the profile
results are significantly affected. The grid-boundary 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 sub-cases (see
Tables 5-5). 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 5-6.
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 (El-Kadi, 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 Sub-Cases
2A, 2B, and 2C are plotted in Figures 5-29 and 5-30. The same general
trend is observed for the numerical solutions in Sub-Cases 2A, 2B, and
2C as in some sub-cases 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 Sub-Case 2C solution

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 three-dimensional plots and longitudinal-vertical three-
dimensional slice profile plots for the three models of Sub-Cases 2A and
2B are presented in Figures 5-31,5-32, 5-34, and 5-35. A cross-section
of Sub-Case 2A is plotted in Figure 5-33. The contour plots for the Sub-
Cases 2A and 2B are plotted in Figures 5-36 through 5-41. In general for
these sub-cases, all of the codes predict the total spatial-concentration
volume reasonably well. A comparison of the results is shown in Table
5-7. Since the actual spatial-concentration volume is misshapen due to
aquifer heterogeneity and unsteady source, the numerical code volumes
are compared to the analytical solution volumes. MOC's Sub-Case 2A
total spatial-concentration volume is only 3.6% less than the analytic
solution model. MOC's Sub-Case 2A spatial-concentration 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 Sub-Case 2A spatial-concentration volume is 4.1% lower in
comparison to the analytic solution. In Sub-Case 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 5-8). MOC is closer than
SUTRA to the analytical outer contours and volumes in all Case 1 sub-
cases. In all sub-cases, MOC's width and length is also closest to the
field data. In contrast, SUTRA has a narrow width for all sub-cases, 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