United States
Environmental Protection
Agency
National Risk Management
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/SR-99/110
January 2000
4>EPA Project Summary
Analytic Element Modeling of
Coastal Aquifers
Mark Bakker, Stephen R. Kraemer, Willem J. de Lange, and Otto D.L. Strack
Four topics were studied concerning
the modeling of groundwater flow in
coastal aquifers with analytic elements:
(1) practical experience was obtained
by constructing a groundwater model
of the shallow aquifers below the
Delmarva Peninsula USA using the
commercial program MVAEM; (2) a
significant increase in performance was
obtained by implementing the theory
for variable density flow in a computer
program that ran on a supercomputer
using vectorization; (3) a new
representation for the density variation
was developed that can simulate the
change from brackish to fresh water
more accurately; and (4) it was shown
that for a specific example of a bell-
shaped transition zone a Dupuit model
gives accurate results unless the bell-
shape is too narrow compared to the
thickness of the aquifer.
This Project Summary was developed
by EPA's National Risk Management
Research Laboratory's Subsurface
Protection and Remediation Division,
Ada, OK, to announce key findings of
the research project that is fully
documented in a separate report of the
same title (see Project Report ordering
information at back).
Introduction
Salt water intrusion is a potential threat
to drinking water supplies in the coastal
areas of the USA. Cities along the coast
are increasing their pumping of
groundwater, which may result in an
increase of chlorides in the well water due
to the upconing of brackish groundwater.
Even a small concentration of chlorides
(> 0.2 grams/liter) will violate the drinking
water standard. To compare, sea water
contains about 18 grams of chlorides per
liter. The potential upconing of brackish
groundwater may be studied by the
simulation of groundwater flow with
numerical models that are designed
specifically for the modeling of flow in
coastal aquifers.
Groundwater flow in coastal aquifers is
affected by the difference in density
between fresh and salt water. The fresh
water is separated from the salt water
through a brackish transition zone, in which
the salinity (and thus the density) of the
water varies from that of salt water to that
of fresh water; salt water is approximately
2.5% heavier than fresh water. The
modeling of the flow generated by variations
in density, the variable density flow, was
the subject of this project. Specifically, it
was investigated whether groundwater flow
in coastal aquifers can be modeled under
the Dupuit approximation in combination
with analytic elements.
The Dupuit Approximation
The construction of numerical models of
three-dimensional variable density flow is
limited by the availability of density data
and by the speed of digital computers. The
Dupuit approximation may be adopted to
reduce the computation time and for the
usual reason of simplicity. The resistance
to flow in the vertical direction is neglected
in Dupuit models, and vertical flow is
governed by continuity of flow only; this
leads to a hydrostatic pressure distribution
in the vertical direction. The Dupuit
-------
approximation is reasonable for flow in
aquifers of large horizontal extend, also
referred to as aquifers with shallow flow or
regional flow. It is possible to construct
regional groundwater flow models of
coastal aquifers by adopting the Dupuit
approximation. The Dupuit theory for
variable density flow may be combined
with any method to model Dupuit flow in a
single density model. The flow field is
modeled with analytic elements for this
project.
The Dupuit theory for variable density
flow assumes that the density distribution
in the aquifer is known at some time. In
practice, the density is known only at a
number of isolated points. We represent
the density distribution with a three-
dimensional function that interpolates
between the points of known density. We
use the multiquadric radial basis
interpolatorforthis purpose. It is noted that
the multiquadric interpolator is a form of
Kriging; the interpolator is identical to
Kriging with a linear variogram if the shape
factor in the interpolator is set equal to
zero.
The density distribution (and thus the
flow field) changes over time as the salt
moves with the groundwater flow. The flow
is incompressible in our model and the
flow field represents the flow at a given
time. The evolution of the density
distribution in the aquifer may be simulated
by numerical integration through time.
During each time step, the velocity field is
fixed and the salt is moved with the
groundwater flow. At the end of a time step
the velocity field corresponding to the new
density distribution is computed and the
process is repeated. This procedure is
also known as successive steady-state
solutions. Second order processes that
affect the salinity distribution, such as
diffusion, are neglected. It may be expected
that this is reasonable for relatively short
times, on the order of 25 years, as is of
interest for most engineering problems.
The Dupuit theory for variable density
has been implemented in the proprietary
computer program MVAEM prior to this
study. MVAEM has been used successfully
to simulate piezometric heads in multi-
aquifer systems in the coastal aquifers of
The Netherlands. The computation time
involved in the construction of these large
regional models was significant. Multi-layer
models that include thousands of points
where the salinity is specified take on the
order of hours to compute a solution on a
high end PC (in 1996). A large portion of
the computation time is used to compute
the effect of the variation in density on the
flow. The computation times ofthetransient
simulations involve a repeated computation
of the solution at different times plus a
large number of evaluations of the velocity
in the aquifer and is of the order of days or
more, as compared to the hours it takes to
obtain one steady-state solution. These
large computation times diminish the
practicality of the modeling approach.
The use of a Dupuit model to simulate
variable density flow raises a number of
additional issues. Transient simulations
require an accurate representation of the
velocity field and the numerical integration
requires an analysis of stability and
convergence. Furthermore, it must be
investigated how accurate a Dupuit model
is to simulate the change of a salt
distribution over time as compared to a
model in which the vertical resistance to
flow is not neglected. Dupuit models are
generally used for regional (shallow) flow
and may become inaccurate when the
shallowflowassumption is not appropriate,
for example near a partially penetrating
well where the flow field changes rapidly
over a distance of several times the aquifer
thickness. And finally, the assessment of
the upconing of saltwater belowa pumping
well needs analysis of the physical stability
of the transition zone itself. A selection of
the issues raised above has been
investigated in this project.
Objectives
The objective of this report is to
investigate the performance of the Dupuit
theory for variable density flow combined
with analytic elements to model
groundwater flow in coastal aquifers. The
four specific areas of interest are (1) the
application of a Dupuit analytic element
model to simulate flow in a coastal aquifer
of the USA, (2) the increase in performance
that may be obtained by using a
supercomputer, (3) the representation of
the density distribution, and (4) the error
introduced by adopting the Dupuit
approximation.
Results and Conclusions
Topic 1. The analytic element
modeling of groundwater flow in
the first confined aquifer beneath
the Delmarva Peninsula.
An analytic element model was
constructed to simulate ground-water flow
in the shallow aquifers beneath the
Delmarva Peninsula (DELaware,
MARyland, VirginiA). The modeling of the
Delmarva peninsula was greatly hampered
by the unavailability of measurements of
the salinity of the groundwater. Only six
measurements were obtained (from
available records) of brackish orsaltwater.
All other measurements indicated fresh
water. The six measurements were
insufficientto construct a three-dimensional
picture of the salinity distribution below the
Delmarva Peninsula. To overcome this
problem, publications of salinity in the
Chesapeake Bay and the Atlantic Ocean
were used to construct a conceptual model
of the salinity distribution below the
peninsula. This distribution results in a
horizontal variation from fresh water on the
peninsula to salt water in the bay and
ocean. The variation inthevertical direction
is assumed to be negligible. This is known
to be unrealistic, but there are not enough
data to propose otherwise.
The resulting model of the fresh water
head was compared to nine measurements
of the head on the eastern part of the
peninsula (the counties of Sussex, DE,
Wicomico, MD, and Worchester, MD). The
comparison showed that the simulation of
fresh water heads in that area are
reasonable, but general conclusions for
the entire model cannot be drawn.
Additional data (head, discharge and
chloride measurements) are necessary to
calibrate the model and to improve the
conceptual model of the density distribution.
The modeling study demonstrated the
many challenges in building groundwater
models of flow systems in coastal aquifers.
The biggest constraints in the Chesapeake
Bay area were the availability of density
data and water table measurements. To
assess the full capabilities of the approach,
it should be applied to an area where more
data are available, perhaps an area outside
the United States.
Topic 2. The reduction of
computation time by the use of a
supercomputer.
The Dupuit approximation for variable
density flowwas implemented in an existing
analytic element model (SLWL) that was
written in FORTRAN. The variable density
routines were written to run on the
CrayC916, a vector machine. The Dupuit
formulation for variable density flow, as
well as the analytic element formulation for
groundwater flow, are suited ideally for
implementation on supercomputers. Both
formulations result in large sums of
complicated functions that are easily
vectorizable. The performance of the
vectorized code was an order of magnitude
better than the un-vectorized code. The
vectorized code performed at a speed of
over 300 billion floating point operations
per second (300 Gflops).
It is noted that high performance
computing seems to move away fromvector
machines to massively parallel machines.
This will have no adverse effect on the
-------
implementation of analytic element codes
on supercomputers. The large sums of
complicated functions can just as easily be
modified to run on massively parallel
machines as on vector machines.
Topic 3. An accurate
representation of the density
distribution.
Forthe investigation ofthefirsttwo topics,
the density distribution was represented
by the three-dimensional multiquadric
interpolator function. The multiquadric
interpolator includes a shape factor A that
controlsthe smoothness ofthe interpolator.
In practice, A is often chosen close to zero.
The main reason forthis choice is to improve
control ofthe interpolator, especially at the
transition from brackish water of variable
density to fresh or salt water of constant
density. It was demonstrated that this
choice will result in reasonable simulations
of the density (and thus the fresh water
head, as was the objective of the modeling
study ofthe Delmarva Peninsula), but that
the resulting velocity vectorwas unrealistic
near the points where the density is
specified. This makes it difficultto simulate
the change of the salinity distribution
through time; during such a simulation the
salt moves with the groundwater and the
velocities have to be accurate.
A new function to represent the density
was introduced to solve this problem. The
representation consists of a number of
surfaces of constant density; the elevations
ofthe surfaces are approximated with two-
dimensional multiquadric interpolators and
the density varies linearly between them. It
was shown that the smoothness parameter
A in the two-dimensional multiquadric
interpolator must be of the order of the
average distance between control points
to obtain reasonable results forthe velocity
field. It is noted that this was not practical
for the original three-dimensional
interpolator.
Topic 4. Implications of adopting
the Dupuit approximation for
variable density flow.
The new density distribution was used to
compare a Dupuit solution with an exact
solution, a solution in which the vertical
resistance to flow is not neglected. The
problem chosen for comparison consisted
of a bell-shaped transition zone. An exact
solution was derived and compared to the
Dupuit solution. The comparison showed
that forthis specific case the Dupuit solution
overestimates the specific discharge vector
and that the flow field becomes inaccurate
when the width ofthe bell-shaped transition
zone is less than two times the thickness
ofthe aquifer.
The newdensity distribution, as proposed
under "Topic 3," represents the transition
from variable density zones (brackish
water) to constant density zones (fresh or
salt water) accurately. The newdistribution
does introduce complications that remain
to be resolved, however, such as the
intersection ofthe transition zone with the
base ortop ofthe aquifer. The new density
distribution is highly useful for the
comparison of Dupuit solutions with exact
solutions, because it is relatively easy to
obtain exact solutions for flow in the vertical
plane corresponding to the new density
distribution. A comparison between exact
and Dupuit solutions may be used to
determine the full range of applicability of
the Dupuit approximation for variable
density flow.
-------
Mark Bakker was with the University of Minnesota during this project but is currently with
the University of Nebraska, Omaha, NE68182.
Willem J. de Lange is with the National Institute of Inland Water Management RIZA,
Lelystad, The Netherlands.
Otto D.L Stack is with the University of Minnesota, Minneapolis, MN 55455.
Stephen R. Kraemer was the EPA Project Officer while at the EPA National Risk
Management Research Laboratory in Ada, OK 74820, and is currently at the EPA
National Exposure Research Laboratory in Athens, GA 30605.
The complete report, entitled "Analytic Element Modeling of Coastal Aquifers" (Order No.
PB2000-101970; Cost: $29.50, subject to change) will be available from:
National Technical Information Service
5285 Port Royal Road
Springfield, VA 22161
Telephone: 703-487-4650
An electronic version of this Project Summary and the Full Report will be available for
download on the Internet at (http://www.epa.gov/ada/kerrcenter.html)
United States
Environmental Protection Agency
Technology Transfer and Support Division (CERI)
Cincinnati, OH 45268
BULK RATE
POSTAGE & FEES PAID
EPA
PERMIT No. G-35
Official Business
Penalty for Private Use $300
EPA/600/SR-99/110
------- |