vvEPA
United States
Environmental Protection
Agency
National Risk Management
Research Laboratory
Ada, OK 74820
Research and Development
EPA/600/S-00/001
May 2000
ENVIRONMENTAL
RESEARCH BRIEF
Analytic Element Modeling of Ground-Water Flow and
High Performance Computing
H.M. Haitjema1, V.A. Kelson2, and K.H. Luther3
Abstract
Several advances in the analytic element method have been
made to enhance its performance and facilitate three-
dimensional ground-waterflow modeling in a regional aquifer
setting. First, a new public domain modular code (ModAEM)
has been developed for modeling ground-water flow on a
regional scale by use of the AEM method. The code can be
compiled as Fortran 90 for use on a single processor machine,
such as a PC, or it can be compiled with a high performance
Fortran (HPF) compiler for use on a massively parallel
processor machine (MPP). The performance increase from
the MPP machine enhances the applicability of the code to
large regional aquifers including local three-dimensional
flow effects. Second, new algorithms were developed, coded,
and tested to expand an earlier solution to a partially
penetrating well in aconfined aquifer to cases of unconfined
flow. The new three-dimensional solution to a partially
penetrating well includes the computation of the phreatic
surface and, when appropriate, a seepage face. The solution
also includes local aquifer stratification. The newSD solution
is suitable for embedding in a Dupuit-Forchheimer model
based on the analytic element method. However, the 3D
solution is mathematically involved and computationally
demanding; it has not been included in ModAEM.
Introduction
The analytic element method (AEM) has been shown to be
uniquely suited for modeling large regional aquifers. The method
employs superposition of analytic functions to create an
' School of Public and Environmental Affairs, Indiana University,
Bloomington, IN 47405.
2 Now at Wittman Hydro Planning Associates, Inc., Ellettsville, IN
47429.
3 Now at Department of Mathematics and Computer Science,
Valparaiso University, Valparaiso, IN 46383.
approximate, but analytic solution to the regional ground-water
flowproblem. The analyticfunctions are called "analyticelements,"
each of which represents some flow feature in the aquifer. For
instance, a line sink is used to model a section of a stream. Since
only the hydrologic features (streams, lakes, etc.)arediscretized
in the analytic element method and not the flow domain in
between these features (no finite difference grid or finite element
mesh), the AEM is relatively insensitive to scale. Regionally, a
rather coarse representation of streams and lakes may be
employed, while locallythe hydrologicfeatures can be represented
in more detail in order to obtain the required model precision. In
fact, three-dimensional solutions can be embedded in the regional
Dupuit-Forchheimer model (horizontal flow model) to further
enhance the realism of the model in the area of interest. This
flexibility in scale and model resolution makes the analytic
element method efficient for modeling large regional aquifers.
Priorto this research, local three-dimensional flow in the analytic
element method was mostly limited to flow near a partially
penetrating well in a confined aquifer. As part of the current
research project, the solution to a partially penetrating well has
been expanded to cases of unconfined flow. Thethree-dimensional
solution to flow near a well in an unconfined aquiferis significantly
more complex than the earlier confined flow solution due to the
presence of a phreatic surface and, in many cases, a seepage
face along the upper portions of the well bore. In addition,
embedding the 3D unconfined solution is not a simple matter of
superposition, as it was for the case of confined flow.
Large regional models with much local detail require many
analytic elements which all contribute to the flow. This leads to
many calculations to arrive at the head or flow vector at a point in
the aquifer. The computational effort is particularly large when
some of the analytic elements represent three-dimensional flow,
like the one for a partially penetrating well. To facilitate large and
complex models that remain computationally fast, we tested the
suitability of the AEM for parallel processing. In principle, the AEM
-------
is highly suitable for parallel processing as the evaluation of the
head orthe flow vector at a point involves very many independent
calculations that can be executed in parallel. As part of this
research, we developed a new public domain code which can be
compiled for a single processor machine (e.g., PC) or an MPP
machine. The new code, ModAEM, is modular in its design to
facilitate future expansions by others.
The current research is reported in detail in two different Ph.D.
theses. The development and testing of the three-dimensional
solution to a partially penetrating well in an unconfined aquifer is
described in "Analytic Solutions to Three-DimensionalUnconfined
Groundwater Flow near Wells" (Luther, 1998). The development
and testing of the code ModAEM and the testing of the AEM on
an MPP machine is found in "Practical Advances in Groundwater
Modeling with Analytic Elements" (Kelson, 1998). The latter
thesis also reports on two alternative procedures to locally
improve the realism of the regional analytic element method. The
first procedure is the extraction of a local U.S. Geological Survey
MODFLOW model out of a regional analytic element model. The
second procedure is to model locally three-dimensional flow in an
approximate fashion using only two-dimensional functions. The
idea is to subdivide the aquifer into layers and model Dupuit-
Forchheimer flow in each layer, while modeling vertical flow
based on leakage between the layers. Both theses are published
by the Graduate School of Indiana University. An electronic copy
can be downloaded from the resources library of the OpenAEM
webpage(http://www.indiana.edu/~aem/). An additional reference
for the 3D well research is Luther and Haitjema [1999].
Three-Dimensional Flow near a Partially Penetrating
Well
Most aquifers of interest for water resource purposes are relatively
thin compared to their horizontal (lateral) extent. Under these
circumstances, resistance to vertical flow may be ignored. This
is known as the Dupuit-Forchheimer assumption. Locally,
however, for instance near a partially penetrating well, resistance
to vertical flow is important and a fully three-dimensional solution
to the flow problem may be needed. An analytic solution to flow
near a partially penetrating well in a confined aquifer has been
implemented before in the code GFLOW [Haitjema, 1995]. This
solution uses a vertical line sink at the center of the well bore to
model the flow into the well. The line sink has a variable sink
density distribution to account for the inflow variation along the
well bore. In fact, the sink density near the ends of the line sink
(and thus well bore) is infinite. The no-flow conditions along the
(horizontal) aquifer top and bottom are realized by use of a
truncated series of image line sinks, in combination with two
semi-infinte line sinks [Haitjema and Kraemer, 1988].
The complete three-dimensional solution is simply superimposed
onto the many two-dimensional solutions that generate the
regional Dupuit-Forchheimer flow field. The resulting solution is
trulythree-dimensional nearthe well and is a Dupuit-Forchheimer
approximation (2D) everywhere else.
As part of the current research project, a three-dimensional
solution for flow near a partially penetrating well in an unconfined
aquifer was developed and tested. This solution is substantially
more involved than the earlier one for confined flow. In the first
place, the unconfined solution involves a phreatic surface near
the well whose shape and location is not a priori known. In many
cases there is also a seepage face along the upper part of the well
bore, the extent (indeed existence) of which is not known in
advance. Finally, the three-dimensional solution cannot be simply
superimposed onto the two-dimensional solutions, as was done
forthe case of confined flow. To further enhance the realism of the
solution, we also implemented aquifer stratification in the vicinity
of the well.
Representation of the Phreatic Surface
We set out to create an analytic solution for flow near a partially
penetrating well that is everywhere regular and differentiate,
including at the boundaries of the flow domain. To accomplish
this, the generating functions for flow, which often contain
singularities, were placed outside the actual flow domain. Thus,
instead of representing the phreatic surface by a facet surface of
single or double layers (sink or doublet panels), as is done in the
classical boundary element method, we introduced a continuous
scalar function %(x,y,z), which is zero along the phreatic surface.
This auxiliary function x may be viewed as the three-dimensional
version of the Zhukovski function used in complex potential
theory to represent a phreatic surface in the vertical plane. In
words, the function % is the hydraulic conductivity times the
pressure head, and can be expressed in terms of the specific
discharge potential (hydraulic conductivity times the head) minus
the hydraulic conductivity times the elevation of the point where
the potential is evaluated. The phreatic surface is then defined by
the surface where the function % is equal to zero, since that is
where the pore pressure (pressure head) is zero. We placed point
sinks and sink rings of different radii above the phreatic surface
to enforce the condition that the flow normal to the phreatic
surface is zero. The sink densities of the line sink along the well
axis and the auxiliary sink distributions above the phreatic surface
are determined in such a manner as to satisfy the following
conditions. Along the well bore, equipotential conditions are
imposed, and along the %=0 surface, no flow conditions or a
specified recharge rate is enforced. These conditions are met
nearly exactly at a number of collocation points on these surfaces,
while they are approximate in between. Due to the non-linear
nature of the equations that must be solved to arrive at the sink
density parameters, an iterative solution procedure is employed
by which the location of the phreatic surface is adjusted during
iterations.
In Figure 1, streamlines toward three wells are depicted in a 3D
diagram. All streamlines start from a vertical array of points,
located in the vertical sheet to the right in the figure. The
horizontal grid is the aquifer bottom, while the thin shaded zone
above the wells is a side view of the (slightly) curved phreatic
surface. The vertical arrays of points above the wells are point
sinks outside the flow domain to enforce no flow conditions along
the phreatic surface. Some streamlines are shown above the
phreatic surface, which enter these point sinks while emanating
from points just above the phreatic surface. The latter streamlines
are, of course, outside the flow domain and shown only to
illustrate the no-flow conditions along the phreatic surface.
Representation of the Seepage Face
If the head in the well occurs below the top of the well bore, a
seepage face may develop along the upper portion of the well
bore. The pressure head along the seepage face is zero, but the
flow normal to the seepage face is not zero. Consequently, only
the condition %= 0 is enforced along the seepage face. The line
sink along the well bore opposite the seepage face is used to
enforce this condition in a number of collocation points along the
seepage face. The solution procedure for the sink density
parameters is iterative since neither the location of the phreatic
surface northe location of the seepage face is known in advance.
In Figure2, we compared our solution to othersolutions presented
in the literature. The problem is a sand box model of flow toward
a well in an unconfined aquifer as published by Hall in 1955. Our
analytic element solution (denoted by AEM in Figure2) compares
-------
, , 1 ....._; i.._.j.—• :__^. I t—L ., , i 4-MW4X
^£=3~f-::±::-i- T-^---V-^r-^^--\----v-^-V^---\^Av1V^:
Figure 1. Streamlines toward three partially penetrating wells in an unconfined aquifer. Dots above the phreatic surface are point sinks to enforce
no-flow conditions along the phreatic surface.
AEM, BIEM, FEM
solutions
Dupuit-Forchheimer
surface
Figure
2. Comparison with the sand box experiment of Hall (1955). Current analytic solution is denoted by the AEM solution.
-------
well with a published boundary integral equation model (BIEM)
[Liggett and Liu, 1983] and a published finite element model
(FEM) [Neuman and Witherspoon, 1970]. All three numerical
solutions, however, exhibit a phreatic surface that is lower than
observed in the experiment. It is noted that published numerical
solutions that include unsaturated flow do match more closely
with Hall's experiment [Hall,1955]. The horizontal lines inside the
well bore, opposite the seepage face, and the horizontal lines
above the phreatic surface, represent sink rings used for the AEM
solution. Sink rings were used inside the well rather than line
sinks because of the relatively large well diameter.
Representation of Aquifer Stratification
Particularly useful are solutions to flow near wells in a stratified
aquifer. Such solutions exhibit the preferential paths of ground-
water flow in the more permeable strata. Conceptually, the
aquifer is subdivided into a number of horizontal layers of varying
thickness and hydraulic conductivity: the aquifer strata.
Consequently, the discharge potential (hydraulic conductivity
times the head) jumps across the interfaces between the layers.
The jump in the potential across an interface is modeled with
double layers, which are dipole distributions over the interface
with their axis normal to the interface (also referred to as doublet
distributions). We used doublet discs with a radially varying
doublet density to account for the radial variation in the potential
jump.
In Figure 3, some streamlines are depicted in a cross-section
over a well in a stratified aquifer. There are three strata with a
conductivity of 10, 50, and 10 feet per day, respectively.
The solutions for stratified aquifers are not as robust as those for
homogeneous aquifers. Since we use radially varying doublet
densities we can only model radially symmetric flow, hence flow
toward a single well. Also, contrasts in the hydraulic conductivity
in excess of a factor of 10 yield inaccurate solutions for the
velocities near the interfaces. Finally, when the phreatic surface
intersects the interface, some local inaccuracies in the flow field
near the intersection may occur.
Embedding 3D Solutions in a Dupuit-Forchheimer
Model
Under confined flow conditions, the three-dimensional solution to
a well reduces to atwo-dimensional (Dupuit-Forchheimer) solution
at some distance from the well (3 to 4 times the aquifer thickness).
Near the well, where the flow is three-dimensional in nature, the
effect of remote flow features (e.g., streams) is two-dimensional
due to the horizontal aquifer boundaries. Consequently, the two-
and three-dimensional solutions in confined aquifers can simply
be superimposed. In so doing, three-dimensional flow is taken
into account where it occurs, near partially penetrating wells,
while all other flow features are represented by two-dimensional
functions. This makes for very efficient modeling of three-
= ft/d k = 50 ft/d
Figure 3. Streamlines toward a partially penetrating well in a stratified unconfined aquifer. The two interfaces between different hydraulic
conductivities are modeled with doublet discs with a radially varying doublet density.
-------
dimensional flow in regional aquifers. Under unconfined flow
conditions, however, the presence of a three-dimensional phreatic
surface near the well prevents the simple superposition of two-
and three-dimensional solutions. Under conditions of unconfined
flow, three-dimensional solutions may be embedded in a Dupuit-
Forchheimer model by defining a 3D zone near one or more
partially penetrating wells. Heads obtained from the Dupuit-
Forchheimer solution are defined on the boundary of that zone,
which contains the partially penetrating well(s), as fully penetrating
with the same pumping rate(s). The solution procedure is iterative
in case the pumping rate of a well is not known in advance; for
instance, when the head at a well is specified instead of the
pumping rate. For a well with a known (specified) pumping rate,
no iterations are needed to couple the 2D and 3D solutions.
Efficiency of the 3D Solution fora Partially
Penetrating Well
The analytic element solution for a partially penetrating well,
developed as part of this research project, issuperiorto previous
solutions in that the discharge potential (thus the piezometric
head) is everywhere differentiable, including at the well bore and
at the phreatic surface. However, the solution requires many
analytic functions (line sinks, point sinks, sink rings, doublet
discs, etc.) to enforce the boundary conditions at collocation
points along the well bore, phreatic surface, and interfaces
between aquifer strata. In fact, to obtain accurate solutions to
non-axisymmetric flow, a least square method was used to
approximate the conditions at many more control points than the
number of degrees of freedom (e.g., sink density parameters)
available at the auxiliary functions. To maximize performance,
some experimentation with the number, nature, and location of
the auxiliary functions was required. In summary, the solutions
are involved, computationally intensive, and not always
straightforward in their implementation.
Analytic Elements on a Massively Parallel Processor
Machine
Analytic element models perform three computationally intensive
tasks: (1) generating a set of equations to solve for the strength
parameters of the analytic elements, (2) solving these equations,
and (3) calculating heads and orfluxes at points in the aquifer. We
tested the performance increase foreachofthese tasks separately
on Indiana University's Intel-Paragon MPP machine with 64
processors. Tasks (1) and (3) consist of a large number of
independent calculations: the computation of the influence of
each analytic element at a number of points in the aquifer. All of
these computations can be carried out in parallel with little or no
communication between processors (no overhead).
Consequently, we found that the computational speed scaled up
linearly with the number of processors employed. Task (2), the
solution of a set of equations with a full coefficient matrix, can also
be speeded up using parallel processors, but not without adding
overhead. The computational speed up of the matrix solution
process, therefore, appeared to be sub-linear. Yet the overall
speed up of analytic element models on parallel processors is
nearly linear, since the matrix solution procedure accounts for
only a small portion of the overall computational effort during
analytic element modeling.
A Modular Parallel Analytic Element Code
We developed a new public domain analytic element code,
ModAEM, with the following characteristics: (1) a modular design
similar to MODFLOW, (2) containing a limited set of basic but
versatile analytic elements, (3) designed to run on both a single
processor machine (e.g., PC orworkstation) and an MPP machine.
The code requires a GUI for efficient operation. An example GUI
is the public domain code WftAEM (U.S. EPA, 2000), which has
been developed for capture zone delineation purposes under
EPA's wellhead protection program. W/?AEM is a Windows 957
98/NT program, which uses ModAEM as its computational engine.
ModAEM has been written in Fortran 90 (F90) and has been
expanded with directives for High Performance Fortran (HPF).
When compiled using a standard F90 compiler it will create a
single processor code, while when compiled with HPF it will run
on an MPP machine. ModAEM is intended as a publicly available
basis forthe implementation of advances in the analytic element
method. ModAEM does not currently include any three-
dimensional analytic elements, such as the partially penetrating
well discussed above.
Summary and Conclusions
We developed and tested a new three-dimensional solution for
flow toward a partially penetrating well in an unconfined aquifer.
The solution includes a phreatic surface and seepage face near
thewell. Forthe axisymmetric case of flowtoward asinglewellwe
also included aquifer stratification. The analytic solution employs
singularity distributions outside the flow domain. Consequently,
the resulting discharge potential is everywhere differentiable,
including at the domain boundaries (e.g., well bore, seepage
face, and phreatic surface). The solution, however, is complex
and computationally intensive.
To increase the computational efficiency of the AEM, we tested
the AEM performance on an MPP machine. The majority of the
computations scaled up linearly with the number of processors
used due to the independent nature of most AEM computations.
The nearly linearoverall performance increase is unusual. Analytic
element models appear particularly suitable for MPP architectures.
We developed a public domain modular analytic element code,
ModAEM, which will run on both a serial machine (e.g., PC) and
an MPP machine. The code uses a limited number of versatile
analytic elements to perform most groundwater modeling tasks.
ModAEM is written in F90 and contains directives for HPF
compilers to adapt it to MPP machines.
We found that the three-dimensional solution to the partially
penetrating well in an unconfined aquifer is involved and
computationally intensive. Development of additional three-
dimensional solutions, for instance, for a shallow creek or lake,
may be expected to be similarly complicated. Kelson [1998]
proposes a technique forapproximatethree-dimensional solutions
in analytic element models, which is conceptually similar to the
procedure used in MODFLOW. The aquifer is subdivided into a
number of horizontal layers with resistance to vertical flow
represented as resistance to flow between the layers ("quasi
three-dimensional flow"). Each layer supports Dupuit-Forchheimer
flow, hence only two-dimensional solutions to flow are being
used. The approach has been demonstrated for relatively high
vertical resistances, such as occur when the vertical hydraulic
conductivity is less than one-tenth of the horizontal conductivity.
However, it is uncertain if accurate solutions can be obtained
under more homogeneous conditions.
Acknowledgments
Dr. Stephen Kraemer served as the U.S. EPA Project Officer. We
thank the following reviewers for their valuable comments:
Dr. David Steward of Kansas State University; Mr. Brent Swartz
of Lockheed Martin Technical Services, Inc.; Dr. Charlie Fitts of
the University of Southern Maine; and Dr. Mark Bakker of the
University of Nebraska-Omaha.
Quality Assurance Statement
This projectwas conducted underan approved quality assurance
program plan and the procedures specified therein were used.
-------
Information on the plan and documentation are available fromthe
Principal Investigator.
Notice
The U.S. Environmental Protection Agency through its Office of
Research and Development partially funded and collaborated in
the research described here under Cooperative Agreement CR-
823637 with Indiana University. The work has been subjected to
the Agency's peer and administrative review and it has been
approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement
or recommendation for use.
References
Haitjema, H.M., 1995. Analytic Element Modeling of Groundwater
Flow. Academic Press.
Haitjema, H.M., andS.R. Kraemer, 1988. A new analytic element
function for modeling partially penetrating wells, Water
Resources Research, 24(5):683-690.
Hall, H.P., 1955. An investigation of steady flow toward a gravity
well, La Houille Blanche, 10, 8-35.
Kelson, V.A., 1998. Practical Advances in GroundwaterModeling
with Analytic Elements, Ph.D. dissertation, Indiana University,
Bloomington (H. Haitjema, advisor).
Liggett, J.A., and P.L.-F. Liu, 1983. The Boundary Integral
Equation Method for Porous Media Flow. George Allen and
Unwin Ltd., London, U.K.
Luther, K.H., 1998. Analytic Solutions to Three-Dimensional
Unconfined GroundwaterFlownearWells, Ph.D. dissertation,
Indiana University, Bloomington (H. Haitjema, advisor).
Luther, K., and H.M. Haitjema, 1999. An analytic element
solution to unconfined flow near partially penetrating wells,
Journal of Hydrology, 226:197-203.
Neuman, S.P., and P.A. Witherspoon, 1970. Finite element
method of analyzing steady seepage with a free surface,
Water Resources Research, 6(3):889-897.
U.S. EPA, 2000. Working with W/7AEM2000: source water
assessmentforaglacial outwash aquifer, Vincennes, Indiana,
EPA/600/R-00/022.
-------
United States
Environmental Protection Agency
Center for Environmental Research Information
Cincinnati, OH 45268
Official Business
Penalty for Private Use $300
EPA/600/S-00/001
BULK RATE
POSTAGE & FEES PAID
EPA
PERMIT No. G-35
------- |