1

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

Analytic Element Modeling of Ground-Water Flow and

High Performance Computing

National Risk Management

Research Laboratory

Ada, OK 74820

United States

Environmental Protection

Agency

Research and Development

ENVIRONMENTAL

RESEARCH BRIEF

EPA/600/S-00/001 May 2000

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 threedimensional ground-water flow 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 a confined 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 new 3D 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

approximate, but analytic solution to the regional ground-water

flow problem. The analytic functions are called “analytic elements,”

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.) are discretized

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 locally the hydrologic features 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.

Prior to 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. The three-dimensional

solution to flow near a well in an unconfined aquifer is 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
2

is highly suitable for parallel processing as the evaluation of the

head or the 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-Dimensional Unconfined

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

truly three-dimensional near the 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

for the 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 differentiable,

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 χ 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 nor the location of the seepage face is known in advance.

In Figure 2, we compared our solution to other solutions 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 Figure 2) compares
3

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.

Hall's surface

Dupuit-Forchheimer 

 surface

<

AEM, BIEM, FEM

 solutions

seepage

 face

R = 195 cm

χ

φw = 30.4 cm

rw = 12.2 cm

 = 0

φR = 122 cm

Figure 2. Comparison with the sand box experiment of Hall (1955). Current analytic solution is denoted by the AEM solution.
4

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.

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 groundwater 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 a two-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 twoand 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-
5

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 twoand three-dimensional solutions. Under conditions of unconfined

flow, three-dimensional solutions may be embedded in a DupuitForchheimer model by defining a 3D zone near one or more

partially penetrating wells. Heads obtained from the DupuitForchheimer 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 for a Partially

Penetrating Well

The analytic element solution for a partially penetrating well,

developed as part of this research project, is superior to 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 or fluxes at points in the aquifer. We

tested the performance increase for each of these 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 or workstation) and an MPP machine.

The code requires a GUI for efficient operation. An example GUI

is the public domain code WhAEM (U.S. EPA, 2000), which has

been developed for capture zone delineation purposes under

EPA’s wellhead protection program. WhAEM is a Windows 95/

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 for the implementation of advances in the analytic element

method. ModAEM does not currently include any threedimensional 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

the well. For the axisymmetric case of flow toward a single well we

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 linear overall 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 threedimensional solutions, for instance, for a shallow creek or lake,

may be expected to be similarly complicated. Kelson [1998]

proposes a technique for approximate three-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 project was conducted under an approved quality assurance

program plan and the procedures specified therein were used.
6

Information on the plan and documentation are available from the

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 CR823637 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., and S.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 Groundwater Modeling

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 Groundwater Flow near Wells, 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 WhAEM2000: source water

assessment for a glacial outwash aquifer, Vincennes, Indiana,

EPA/600/R-00/022.
7

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