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

-------