EPA/600/B-18/089

June 2018

Working with WhAEM

Demonstration of Capture Zone Delineation for a
City Wellfield in a Valley Fill Glacial Outwash
Aquifer for Wellhead Protection

By

Stephen R. Kraemer, Ph.D.

U.S. Environmental Protection Agency
Office of Research and Development
Los Angeles, CA 90017
Henk M. Haitjema, Ph.D.
Haitjema Consulting, Inc.
Bloomington, IN 47401

Office of Research and Development

U.S. Environmental Protection Agency
Washington, DC 20460
ii

Notice
This document has been subject to the Agency’s peer and administrative review, and it has been approved for
publication as an EPA Report. Model verification test files are described in the WhAEM Quality Assurance
Project Plan.
The use of geographic maps and data in the example case study is for demonstration purposes and not
part of an approved source water protection plan.
Mention of trade names of commercial products does not constitute endorsement or recommendation for
use.
The information in this document has been funded (wholly) or (in part) by the United States Environmen-
tal Protection Agency under contracts (P.O. 8L-0425-NTSA, P.O. 0L-2294-NASA, P.O. EP-17-W-000076)
to Haitjema Consulting and under contract (68W01032) to Computer Sciences Corporation.
This document was typeset using WinEdt of Aleksander Simonic (www.winedt.com) and MikTEX(including
LaTEX2) software of Christian Shenk (www.miktex.org). TEX is a trademark of the American Mathematical
Society.
 iii

Abstract
The purpose of this document is to demonstrate through an example study the use of the ground water
geohydrology computer program WhAEMTM for Windows, which was developed by the U.S. Environmental
Protection Agency (EPA). The Wellhead Analytic Element Model (WhAEM) is an open source ground-
water flow model designed to facilitate capture zone delineation and protection area mapping in support
of the states’ and tribes’ Wellhead Protection Programs (WHPP) and Source Water Assessment Planning
(SWAP) for public water well supplies in the United States. Program operation and modeling practice is
covered in a series of progressively more complex representations of a well field tapping a glacial outwash
aquifer. WhAEM provides an interactive computer environment for design of protection areas based on
simple capture zone delineations (e.g., radius methods, well-in-uniform-flow solutions) and geohydrologic
modeling methods. Protection areas are designed and overlain upon U.S. Geological Survey Digital Line
Graph (DLGs), Digital Raster Graphic (DRGs) or other electronic base maps. Geohydrologic modeling for
steady pumping wells, including the influence of hydrological boundaries, such as rivers, recharge, no-flow
boundaries, and inhomogeneity zones, is accomplished using the analytic element method. Reverse gradient
tracelines of known residence time emanating from the pumping center are used to delineate the capture
zones. WhAEM has import and export utilities for DXF files and Shapefiles. WhAEM has on-line help
and tutorials. Install scripts and base maps are available for download from the EPA Center for Exposure
Assessment Modeling web site (www.epa.gov/exposure-assessment-models/whaem).
iv

Foreword
The EPA’s National Exposure Research Laboratory (NERL) has a valued reputation for supporting the
Agency’s mission of protecting human health and the environment with multidisciplinary expertise that
brings cutting-edge research and technology to address critical exposure questions and to develop approa-
ches for reducing harmful exposures. The NERL Systems Exposure Division (SED) integrates exposure data,
tools, methods and models to assess cumulative exposures to humans and ecosystems. This assessment uses
systems-based approaches to organize current knowledge, simulate potential futures and synthesize results.
SED translates information into decision-making frameworks that integrate human and ecological, multime-
dia, socio-economic and cumulative risk considerations to improve human well-being and the sustainability
of the built and natural environment. The holistic view on exposure science and its linkage to sustainable
solutions emphasizes cross-divisional and laboratory integration.
To support this mission, SED is comprised of dedicated scientific and administrative professionals who
deliver high-quality, timely, effective and understandable products and applications that meet our customers
needs, with a commitment to exceed their expectations.
The development and release of the updated Wellhead Analytic Element Model (WhAEM) is in support
of the state’s and tribe’s Source Water Assessment and Wellhead Protection Programs as required by the
Safe Drinking Water Act. WhAEM is a general purpose ground water hydrology desktop modeling tool based
on the innovative computational solution technique known as the analytic element method. WhAEM assists
in the delineation and mapping of the area contributing source water to pumping wells, such as public
water supply wells. The development of WhAEM has been supported by the EPA Office of Research and
Development and the EPA Office of Ground Water and Drinking Water since 1994.

Jay Garland, Ph.D.

Director
Systems Exposure Division
Cincinnati, Ohio
Contents

1 BACKGROUND 1
1.1 Source Water Protection Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Modeling Philosophy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 The Emergence of the WHPA Model in the 1990s . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 The Evolution of WhAEM to the Present . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.5 The Base Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.6 Exercise: Exploring a USGS quad map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 PROTECTION ZONES I – DISTANCE CRITERIA 15

2.1 Setbacks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.1 Exercise: Glacial outwash wellfield . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 PROTECTION ZONES II — RESIDENCE TIME CRITERIA 19

3.1 Calculated Fixed Radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.1 Exercise: Calculated fixed radius for city wellfield . . . . . . . . . . . . . . . . . . . . . 21
3.2 Well in Uniform Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.1 Exercise: Well in Uniform Flow for city wellfield . . . . . . . . . . . . . . . . . . . . . 24

4 SIMPLE WHPAs 33
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.2 Back-of-the-Envelope WHPAs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.2.1 Exercise: Simple WHPAs and no ambient flow . . . . . . . . . . . . . . . . . . . . . . 35
4.2.2 Exercise: Simple WHPAs with ambient flow . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2.3 Exercise: Simple WHPAs with stream boundary . . . . . . . . . . . . . . . . . . . . . 38
4.2.4 Exercise: Simple WHPAs with multiple wells . . . . . . . . . . . . . . . . . . . . . . . 41

5 PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA 43

5.1 Introduction to the Analytic Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2 Rivers and Line-sinks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.1 Exercise: Study area and Line-sinks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.3 Geologic Contacts and No-Flow Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.3.1 Exercise: No flow boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.4 Aquifer Inhomogeneities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.4.1 Exercise: Inhomogeneities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.5 Riverbed Resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

v
vi CONTENTS

5.5.1 Exercise: Resistance Linesinks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

6 SUMMARY 65
6.1 Good Modeling Practice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.1.1 Rules-of-Thumb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.1.2 Choosing the Right Tool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Bibliography 71

A Conversion Factors 73

B WhAEM Base Maps and the Internet 75

B.1 Vector Basemaps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
B.1.1 Getting Binary Base Maps (BBM) from the EPA Web Server . . . . . . . . . . . . . . 75
B.1.2 USGS Digitial Line Graphs (DLG) on EarthExplorer . . . . . . . . . . . . . . . . . . . 75
B.1.3 USEPA BASINS maps (Shapefile) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
B.2 Raster Basemaps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
B.2.1 USGS digital raster graphic (DRG) maps on National Map Viewer . . . . . . . . . . . 77
B.2.2 GoogleEarth maps (tif) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
B.3 USGS National Water Information System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
CONTENTS vii

Acknowledgements
Thanks to all participants involved in the evolution of the code, including the participants in the EPA
Regional Training Courses: (R3 Philadelphia, June, 1999, August 2001, December 2003, June 2010), (R4
Atlanta, March 2004),(R5 Chicago, May, 1999), (R8 Denver, February 2001), (R9 San Francisco, April 2002,
September 2005, December 2007), (R10 Seattle, January, 2000), (HQ, March, 2000). Thanks to our initial
external peer reviewers, Dr. Steve Silliman of the University of Notre Dame, and Dr. Todd Rasmussen of
the University of Georgia. Thanks to Dr. Jim Weaver of EPA-ORD for his early reviews. Dr. Robert Swank
provided much appreciated and thorough review of the 2007 update. Dr. Susan Dakin (www.sci-writer.com)
provided editorial review of Chapter 1. This current update benefitted from the peer review comments of
Douglas Rambo, Delaware Department of Natural Resources and Environmental Control, Rick Arnold, EPA
Region 8, Denver, CO, David Rebhuhn, Oak Ridge Associated Universities, Corvallis, OR, and Dr. Mark
Bakker, Technical University Delft, The Netherlands. Thanks for support from the EPA Office of Ground
Water and Drinking Water, Washington DC.
viii CONTENTS

This page is blank.

Chapter 1

BACKGROUND

The objective of this report is to introduce and demonstrate the the use of the geohydrology computer model
the Wellhead Analytic Element Model for Windows (WhAEM) for delineating capture zones for pumping
wells. The presentation covers both basic model operations and modeling practice. We use a capture zone
delineation case study to illustrate a stepwise and progressive modeling strategy. The first chapter provides
the introduction.

1.1 Source Water Protection Programs

Ground water is a primary source of public drinking water in many parts of the country. Ground water is
a renewable resource if properly managed, and its quality can be excellent, as a result of natural cleansing
by biologically active soil cover and aquifer media. In the past decades, it has come to the public attention
that ground water wells can be threatened by over-drafting and pollution. Source water protection is based
on the idea of: (1) protecting public water wells from contamination; and (2) working to maintain the
quality of the resource. Source water protection means taking positive steps to manage potential sources of
contaminants and contingency planning for the future by determining alternative sources of drinking water.
The Safe Drinking Water Act requires the United States Environmental Protection Agency (EPA) to develop
a number of programs that involve source water protection for implementation at the state level, such as the
Wellhead Protection Program (WHPP) and the Source Water Assessment Program (SWAP) [USEPA,2017].
The main steps of the source water protection process involve assessment of the area contributing water
to the well or wellfield, a survey of potential contaminant sources within this area, and evaluation of the
susceptibility of the well to these contaminants. Susceptibility evaluation involves the evaluation of conta-
minant release and the likelihood that contaminants will be transported through the soil and aquifer to the
well screen. This report will focus on the first step — protection zone delineation.
Designation of a wellhead protection area is thus a commitment by the community to source area mana-
gement. Delineation of the wellhead protection area is often a compromise between scientific and technical
understanding of geohydrology and contaminant transport, on one hand, and practical implementation for
public safety on the other.
The EPA Office of Ground Water and Drinking Water (OGWDW) established guidance on the criteria
and methods for delineating protection areas [USEPA,1993], [USEPA,1994]. The criteria include (1) distance,
(2) drawdown, (3) residence time, (4) flow boundaries, and (5) assimilative capacity. Methods for applying
the criteria range from the simple to the complex, including setback mapping, the calculated radius method,

1
2 CHAPTER 1. BACKGROUND

Table 1.1: Wellhead protection area delineation criteria and methods.

METHODS arbitrary calc. well in geo- geo- transport/
fixed fixed uni- hydrolo- hydrologic transfor-
radius radius form gic modeling mation
CRITERIA flow mapping modeling
field
distance x x
draw down x x x x
residence time x x x x
geohydrologic x x x
boundaries
assimilative capability x

and geohydrologic modeling (see Table 1.1). This case study document reviews all of the criteria and their
methods using WhAEM with emphasis on the residence time criterion and methods. The evaluation of
assimilative capability is an active area of research but not discussed in detail in this document.
The EPA OGWDW has created the online mapping tool called the Drinking Water Mapping Application
to Protect Source Waters (DWMAPS), which in concert with other state and local mapping tools, can be used
to update their source water assessments and protection plans (http://www.epa.gov/sourcewaterprotection/dwmaps).

1.2 Modeling Philosophy

The philosophy presented below is for deterministic mathematical modeling; a more complete discussion
can be found in Konikow and Bredehoeft [1992]. The reader is encouraged to explore alternative modeling
approaches based on stochastic theories [Vassolo et al., 1998] , [Feyen et al., 2001]. A discussion of the merits
of both approaches is presented in Molz [2015].
Modeling of flow and transport through geohydrologic systems is a complex undertaking. Our presen-
tation of good modeling practice revolves around stepwise problem solving. The mathematical model is an
abstraction of the real world, and all abstractions are simplifications. The challenge is to make the abstraction
represent the essence of the natural system as required to address the problem at hand. Representation of
the advective flow system is understood to be prerequisite to modeling transport and transformations.
Formally speaking, mathematical modeling of ground-water flow is concerned with finding solutions to
a governing differential equation subject to a set of boundary conditions. This procedure can be simple
or complex, depending on the ground-water flow processes included in the differential equation and the
complexity of the boundary conditions. That combination of processes and boundary conditions is referred
to as the “conceptual model.” An example of a very simple conceptual model is a single well in a homogeneous,
confined aquifer near a straight, infinitely long equipotential line (line of constant head). Such a conceptual
model is a severe abstraction of a well near a stream or lake boundary. The actual meandering surface water
boundary is replaced by an infinitely long straight line, the heads underneath are assumed to be equal to
the surface water elevation, aquifer properties are assumed uniform and homogeneous throughout the flow
domain, and the ground-water flow is assumed to be strictly horizontal. The virtue of this conceptual model
is the ease with which a solution can be obtained (Thiem’s solution for a well with an image recharge well
across from the equipotential line). Of course, the price paid is the lack of realism of this conceptual model,
1.2. MODELING PHILOSOPHY 3

compared with the complex reality of a partially penetrating well near a meandering stream with a silty
stream bottom in a heterogeneous aquifer of varying thickness.
A realistic conceptual model is often thought of as a representation of reality that includes most, if not
all, of the geologic and hydrologic complexities of an aquifer or aquifer system. The more realistic (read
complex) a conceptual model is, the more involved the modeling becomes, mostly for the following reasons:

1. Field data requirements can escalate when the modeling effort strives for maximum realism.

2. Computer codes needed to solve complex conceptual models are difficult to operate.

3. Non-uniqueness invades the modeling process — many different sets of parameters and boundary
conditions can represent the same observed water levels and fluxes.

While incorporation of more and more “realism” in the conceptual model might suggest an increasingly
accurate representation of the ground-water flow regime, the reality can be quite different. Increased data
requirements and the associated data uncertainty, as well as the difficulty of interpreting the results of
complex (multi-parameter) models, can defeat the theoretical improvement in model realism. This is true
not only for ground-water flow modeling. Albert Einstein is believed to have offered the following advice when
talking about models of the universe: “Things should be made as simple as possible . . . but not simpler.”
In other words, good modeling practice requires a balance between the realism of the conceptual model and
the practical constraint of parameterizing the resulting ground-water flow model.
When deciding on an appropriate conceptual model, it is important to take into account not only the
complexity of the geohydrology, but also the objective of the modeling exercise. In the context of wellhead
protection, the purpose of the modeling is to delineate a time-of-travel capture zone for one or more drinking
water production wells. For the purpose of wellhead protection, the delineation of the time-of-travel capture
zone is based on the assumption of steady-state ground-water flow and average ground-water travel times.
A time-of-travel capture zone, defined in this manner, is arguably not the most relevant entity when it
comes to protecting drinking water wells from potential toxic substances. The movement of contaminants
toward a drinking-water well is a complex contaminant transport problem where such processes as dispersion,
adsorption, and (bio)chemical reactions strongly influence the degree to which a well is actually exposed to
contaminants in the aquifer. Incorporation of all these aspects in the definition of a wellhead protection
zone would require sophisticated transport modeling for a whole series of different contaminants. The field
data requirements for such an endeavor would be staggering, and the technology to reliably describe these
contaminant transport processes is not completely developed. In any case, the definition of multiple wellhead
protection areas, one for each potential contaminant, is impractical. Consequently, many states base their
wellhead protection areas on time-of-travel capture zones that result from considering average ground water
travel times and steady-state (average) flow conditions.
A goal of the ground-water flow modeling is to find a proper balance between a simple conceptual model
and a realistic one. The most economical way to judge whether a conceptual model is not too simple is
to gradually increase its complexity and investigate the effect on the purpose of the modeling. A good
example is the classical dilemma of modeling two-dimensional flow or three-dimensional flow. Most modelers
feel that using a horizontal flow model (based on the Dupuit assumption) is inadequate to represent the
three-dimensional world. Intuition, however, offers poor guidance here. Only a comparison between a
three-dimensional solution and two-dimensional one can answer this question. This is true for all aspects
of conceptual model design. The need to include such factors as aquifer heterogeneities, varying recharge
rates, and resistance to flow in or out of surface water features can be assessed by comparing ground-water
4 CHAPTER 1. BACKGROUND

flow solutions with and without these features. In all cases it is important to look at the impact on the
time-of-travel capture zone, since that is the endpoint of the modeling. Once progressively more complex
conceptual models do not change the calculated time-of-travel capture zone in a meaningful way, the modeling
is considered complete. It may sometimes be preferable to trade data acquisition and model complexity for
a larger, thus more protective, wellhead protection zone. Such a conservative time-of-travel capture zone can
be obtained by making assumptions that are known to lead to a larger capture zone than the zone that may
result from a more detailed conceptual model. Also, modeling experience will lead to general rules of thumb
in formulation of the conceptual model. These topics are visited in the discussion in the final chapter.

1.3 The Emergence of the WHPA Model in the 1990s

Nearly all ground-water modeling projects are constrained by budget and time. This is particularly true for
time-of-travel capture zone delineation in the context of wellhead protection. Most states were confronted
with the task of developing wellhead protection programs for thousands of public drinking water supply
systems in only a few years. Conducting an extensive ground-water modeling campaign for each individual
drinking water well (or well field) was out the question in view of both the time involved and the cost.
The EPA recognized this reality from the start and proposed a series of simplified capture zone delineation
methods to facilitate a timely implementation of the States’ wellhead protection programs:

1. arbitrary fixed radius method

2. calculated fixed radius method

3. simplified variable shapes method

4. well in a uniform flow method

5. geologic mapping

6. ground-water flow modeling based on a simple conceptual model

7. ground-water flow modeling based on a sophisticated conceptual model

While only the last, all-out modeling exercise (no. 7) will completely satisfy the ground-water hydrologist,
the more approximate approaches are offered in order to get the job done. The rationale is that it is better to
define a less-than-perfect wellhead protection zone than none at all. The task of the ground-water hydrologist
is to select the best possible approach within the time and budget constraints of the wellhead protection
program.
The first two methods lead to circular wellhead protection zones and often over-protect some areas while
not protecting other areas near the well.
The well in a uniform flow method takes ambient ground-water flow into account (approximately) and
offers a significantly better approximation of the actual time-of-travel capture zone. This approach was
implemented in a computer code, the Wellhead Protection Area (WHPA) model, made public by the EPA
[Blandford and Huyakorn, 1993]. In addition, WHPA allowed for wells in flow fields generated by enclosed
(rectangular) boundaries, such as provided by combinations of canals (constant head) or geologic contacts (no
flow). WHPA used particle tracking to draw bundles of streamlines toward the well, including the boundary
streamlines that define the capture zone for the well. The model delineated a time-of-travel capture zone by
1.4. THE EVOLUTION OF WHAEM TO THE PRESENT 5

drawing streamlines with a constant ground-water residence time, from starting point to the well. WHPA
was menu driven, had a Windows-like graphical user interface under DOS, and required minimal computer
resources.

1.4 The Evolution of WhAEM to the Present

The WhAEM (pronounced wăm) project offers an evolutionary but significant step. In addition to facilitating
the simplified techniques (e.g., the fixed radius, simplified shapes, and uniform flow methods), the modeling
environment allows incorporation of the principal hydrologic boundaries (surface waters) in the regional
ground-water flow domain and calculates ambient ground-water flow patterns based on aquifer recharge due
to precipitation in the presence of these surface waters. Consequently, the ambient flow does not have to be
specified a priori and will not be limited to a uniform flow field. Also, surface water features and geologic
boundaries are not limited to infinitely long straight lines as in WHPA, but can be of arbitrary shape, further
improving the realism of the conceptual model. On the other hand, WhAEM currently lacks such features as
continuously varying aquifer bottom elevations, three-dimensional elements, multi-aquifer flow, and transient
flow, which are available in professional ground-water models. While these omissions limit the applicability of
WhAEM they also make the code easier to operate and thus easier to learn. Furthermore, WhAEM requires
less input data. Time-of-travel capture zone delineation with WhAEM thus is considerably cheaper than when
using a full-feature ground-water flow model, such as, the U.S. Geological Survey MODFLOW [McDonald
and Harbaugh, 1988]. Obviously, the price is a less sophisticated (potentially less accurate) definition of the
time-of-travel capture zone. The full implementation of WhAEM capabilities offers more realism than the
first five approaches in the EPA list, but this also requires more time and higher cost. In some cases, this extra
time and expense may be justified; in some cases, one of the simpler methods may be adequate; and in some
cases, a more powerful ground-water model than WhAEM must be employed. It is the responsibility of the
model user to determine whether WhAEM is an adequate tool for the delineation of a time-of-travel capture
zone for a particular public water supply system. That decision, as already explained, should be based on
both technical merit of alternative approaches and the practical requirements for timely implementation of
the wellhead protection program.
The release of the Capture Zone Analytic Element Model (CZAEM) for modeling ground-water flow,
[Strack et al., 1994], introduced a new analytical solution technique to the wellhead protection community
– the analytic element method. This method is based on the principle of superposition of many closed-
form analytic functions, each representing a hydrological feature, such as point-sinks for wells, line-sinks for
rivers, area elements for zones of effective recharge, and line doublets (double layers) for geologic contacts
[Strack, 1989], [Haitjema, 1995],[Strack,2017]. The two-dimensional analytic element models invoke the
Dupuit assumption [Dupuit, 1863], meaning that resistance to vertical flow is negligible and that the head
is constant with depth (zero vertical head gradients). The Dupuit assumption has been shown to be very
powerful for modeling flow in shallow aquifers. The Dupuit assumption is reasonable for capture zone
delineation when the capture zone is an order of magnitude wider than the saturated thickness of the
aquifer. The significance of the Dupuit assumption is that it allows the approximation of a three-dimensional
flow field with a two-dimensional one. CZAEM included analytic elements for wells, line-sinks, uniform
flow, and circular recharge zones, and user interaction was through a DOS command line. CZAEM had
sophisticated routines for drawing time-of-travel-based capture zone envelopes and time and source-water
sub-zones [Bakker and Strack, 1996].
In the first release of the modeling system WhAEM , a graphical user interface (GUI) and geographic ana-
6 CHAPTER 1. BACKGROUND

lytic element preprocessor, GAEP [Kelson et al., 1993] were linked to the solution engine CZAEM [Haitjema
et. al, 1995], [Haitjema et. al, 1994]. This version of WhAEM ran under Windows 3.1. The development
of GAEP opened the door for on-screen conceptual model building, where the model elements (e.g., wells,
line-sinks) are superimposed on top of binary base maps showing hydrography and roads. It is interesting
to note that the conceptual modeling approach pioneered by GAEP is now being implemented in advanced
modeling packages, such as the Groundwater Modeling System [Aquaveo, 2017]. WhAEM for Windows 95
was released in 1997. The release of WhAEM version 1 (April 2000) included an improved Windows GUI,
the use of U.S. Geological Survey (USGS) digital line graphs (DLGs) as the base map, and a new solu-
tion engine, ModAEM [Kelson, 1998], which included arbitrarily shaped no-flow boundaries. The release
of WhAEM version 1 supported delineation by fixed radius, calculated radius, well in uniform flow, and
steady flow in confined/unconfined aquifers with constant aquifer properties and recharge. WhAEM version
2 incorporated the GFLOW1 solution engine [Haitjema, 1995] and included new elements: inhomogeneities
and resistance line-sinks.
The current release, WhAEM version 3, has additional support for raster base maps, such as USGS
digital raster graphics (DRG) and shapefiles (SHP). Also, WhAEM has simple WHPA tools to include the
functionality of the EPA WHPA code. WhAEM runs under Windows 7 and 10. WhAEM is available for
download from the EPA Center for Exposure Assessment web site (https://www.epa.gov/exposure-assessment-
models/whaem).
Developments in analytic elements may be followed at the community web site (www.analyticelements.org).
A review of the early history and contribution of the analytic element research program is discussed in
[Kraemer, 2007].

1.5 The Base Map

Ground-water models are applied to regular Cartesian coordinate systems, even though the surface of the
earth is curved. Cartographers have solved that problem by projecting the surface of a spheroid, where
position is measured in latitude and longitude, to a rectangular coordinate system on a flat piece of paper
or computer screen. In addition, the origin of the ground-water model coordinate system is defined in a
standard way so that maps and spatial data from different sources can be overlapped.
WhAEM can use any Cartesian coordinate system, such as Universal Transverse Mercator (UTM), State
Plane, or local site coordinates. The WhAEM graphical user interface (GUI) provides both pre- and post-
processing in the base map coordinates, while the computational engine operates in local rectangular coor-
dinates for reasons of numerical accuracy; however, this is hidden from the user.
Standard base maps may include US Geological Survey (USGS) electronic vector and raster maps. The
digital line graph (DLG) series includes lines for hydrography (e.g., rivers, lakes, etc), roads, hypsography (to-
pography), and other line-features. DLG maps are available at convenient scales for ground-water modeling,
including the intermediate scale (1:100,000) and the large scale (1:24,000) and are available for download
(see Appendix B.2).
The digital raster graphics (DRG) series are scanned bitmaps of the paper topographic maps, and are
also available at intermediate and large scales; DRG maps are also available for download (See Appendix
B.2.1).
The USGS maps use the UTM projection [USGS, 1997]. In this grid system, the globe is divided into
60 north-south zones, each covering a strip 6 degrees wide in longitude. The zones are numbered from 1 to
60; zones 10 to 19 cover the conterminous U.S. (see Figure 1.1). In each zone, coordinates are measured
1.6. EXERCISE: EXPLORING A USGS QUAD MAP 7

Figure 1.1: UTM zones of the conterminous U.S.

north and east in meters. The northing values are measured continuously from zero at the Equator, in a
northerly direction. A central meridian through the middle of each 6-degree zone is assigned an easting value
of 500,000 meters. Grid values to the west of this central meridian are less than 500,000 meters; to the east
they are more than 500,000 meters. Base maps do not join smoothly across UTM zones.
For performance reasons, a compressed binary form of the DLG called binary base map (BBM) has been
developed. WhAEM reads BBM files from the disk each time it redraws a map on the screen. The EPA Center
for Exposure Assessment (CEAM) operates a web server for download of BBMs (www.epa.gov/exposure-
assessment-models/whaem2000-bbm-files-us) . The CEAM web server has a graphical index, and it can be
accessed with a browser to select specific maps for a project. This procedure is described in Appendix B.
Another important projection is the State Plane Coordinate System (SPCS), wherein each state has its
own zone(s). The number of zones in a state is determined by the area the state covers and ranges from one
for a small state such as Rhode Island to as many as five. The projection used for each state is variable;
some States use a transverse Mercator projection, while some States use a Lambert conformal projection.
The WhAEM map viewer is compatible with georeferenced vector maps in the BBM or shapefile formats,
and georeferenced raster maps in the geotiff format.

1.6 Exercise: Exploring a USGS quad map

We will demonstrate source water area delineation for a city well field in a glacial outwash valley of a large
river in midwestern USA. This exercise will be used to: (1) familiarize the reader with the basic operations
of WhAEM ; and (2) illustrate some basic ground-water modeling concepts.
The floodplain of the river is characterized by valley fill from glacial outwash, and consists of coarse sand
and fine gravel on top of underlying bedrock. The area is assumed to receive about 40 inches per year of
precipitation. The areal recharge to the aquifer is estimated to be 12.1 inches per year (in/yr), the average
hydraulic conductivity of the local aquifer to be 350 feet per day (ft/d), and the effective porosity to be 0.2
8 CHAPTER 1. BACKGROUND

Figure 1.2: WhAEM splash screen.

[Shedlock, 1980]. The wellfield pumped about 2.75 million gallons per day(mgd) in 1977, and 3.4 mgd in
1999, and about 3.8 mgd in 2005. You will learn more details about the wellfield as you progress through
the case study.
In this exercise, we will place the wellfield in its geographic context. You will locate the wellfield on the
digital U.S. Geological Survey 7.5 minute DRG topographic map (vincennes 24k.tif). The DRG is available
at the CEAM web data download area (www.epa.gov/exposure-assessment-models/whaem2000)).
We also assume user familiarity with the Microsoft Windows operating system, and basic mouse operation.

Step 1 Start WhAEM and familiarize yourself with the graphical user interface (GUI).

From the Start button, go to All Programs, start WhAEM . The “About WhAEM” splash screen
(See Figure 1.2) will appear briefly, and credits are given.
The workspace for WhAEM has standard windows features: a Menu bar on top; smart icons bar under-
neath the menu bar; and status bar on the bottom of the window (Figure 1.3).
Step 2 Create a new project database.
From the Menu Bar, drop down the Project options. Select New Database; then select your project
folder (e.g., where you stored your BBM’s, such as My Documents folder, WhAEM Projects folder, Vincennes
folder; and type
“vinbase” into the dialog box for the File name (See Figure 1.4). Click Open.
The “New Database Wizard” will prompt you for a project description and the units associated with
the input data. Enter ”Glacial outwash wellhead case study” as a project description (See Figure 1.5).
WhAEM requires consistent units, either meters and days or feet and days. The “Distance Units in the
Basemap Files” will be dictated by your source files. Select meters for the “Distance Units in Basemaps
Files” by clicking on the radio button. The provided basemaps are in UTM coordinates so you must select
meters. You are free to select the “Units for Computations”; the computer program will enforce consistent
units for future data entry. Select feet and days for this exercise (See Figure 1.5). After you have made
1.6. EXERCISE: EXPLORING A USGS QUAD MAP 9

Figure 1.3: WhAEM user interface.

Figure 1.4: Create new database.

10 CHAPTER 1. BACKGROUND

Figure 1.5: New database wizard.

these selections, click on Create Database.

Repeat: Your base map units are actually not a choice, but are dictated by the units of your electronic source. The
coordinates displayed on-screen will match your map units. The units for computation are your choice. For instance,
if you select feet and days for your units for computations, the hydraulic conductivity must be entered in feet per
day and well pumping rate must be entered in cubic feet per day. The GUI will convert the base map units to the
computational units when forwarding data to the solver, so that it works in consistent units.
The next window facilitates the definition of base maps (in binary format). The window “Currently Used
Base maps” will be empty at this point. Click on Add Map to display the available *.bbm files. Change
to the folder where the basemaps are stored. Click to select the file vincennes 24k.tif. Then click Open.
That file will be displayed in the box “Currently Used Base Maps”. Enter in the Base Filename (DOS) box
”vinbase”. (See Figure 1.6). Click OK and wait a few seconds. The map will be loaded and plotted on the
project workspace (See Figure 1.7).
From the smart icons bar, test the zoom in, zoom out, zoom to extents, zoom to window icons, scrolling,
and re-centering (See Figure 1.8).

Zoom In to the city, and find the city wellfield (See Figure 1.9). Note the UTM coordinates of the
pumping center (approximate center of the wellfield) reported in the Status Message bar in the lower right
corner of the Workspace (approximately 452650 m, 4280665 m). Pan-up the river to the north and locate
where the 400 foot topographic line crosses the river.
You can add vector maps to the basemap window. Go to Project , Project Settings, Add Map, and
shift-click to select all of the vi*.bbm in the examples folder.
You can toggle the maps on and off using View , Base Map, Raster Graphics or Vector Graphics.
Continue to the next chapter or exit the program by clicking Project and Exit.
1.6. EXERCISE: EXPLORING A USGS QUAD MAP 11

Figure 1.6: Project settings.

12 CHAPTER 1. BACKGROUND

Figure 1.7: The WhAEM base map view of USGS DRG.

Figure 1.8: The View smart icons.

1.6. EXERCISE: EXPLORING A USGS QUAD MAP 13

Figure 1.9: Zoom-in view of city wellfield.

.
14 CHAPTER 1. BACKGROUND

This page is blank.

Chapter 2

PROTECTION ZONES I – DISTANCE

CRITERIA

2.1 Setbacks
One of the basic criteria for defining a protection zone surrounding a pumping well or wellfield is distance.
The method of designating a setback based on a fixed radius is commonly used in state source water pro-
grams [Merkle et al.,1996]. For example, at one time, Indiana had a setback distance of 200 feet between
potential sources and public drinking water wells. The fixed radius is the first line of defense against surface
contaminants that could reach the wellhead and travel down an improperly sealed borehole into the ground
water near the well screen. A fixed radius wellhead protection zone is also used as a setback for potential
sources of microbial pathogens. The assumption is that pathogens entering ground water outside the setback
area will become inactive before entering the well. The actual survival time of pathogens in ground water,
however, is still the topic of ongoing research. However, a setback is not likely protective from a gasoline
spill and plume containing the additive methyl tertiary butyl ether. MTBE resists degradation in many
subsurface environments.

2.1.1 Exercise: Glacial outwash wellfield

Step 1 In this step you will locate a well in the center of the city wellfield, and draw the setback protection
area.

Open (or restore to the view) the project database vinbase, and make a duplicate database “setback1”:
(Project , Make Duplicate Database, setback1). WhAEM stores project information as it is entered,
as is common in database programs. This is in contrast to word processing programs, where you are required
to enter information, and then save (or autosave) to the file. Therefore, it is recommended that you make a
new project database anytime you start a new exercise.
Zoom to the city wellfield. Position the cursor near the pumping center (center of the wellfield) and read
off the UTM coordinates. Note that 452650 m, 4280665 m is about the pumping center.
To create a well, from the Menu select Element , then select New and click on Well. Move the cursor
to the approximate location of one of the five wells, and click the left mouse button. A dialog box will open
as in Figure 2.1. Type in a name for a label, e.g. ‘‘Well 1’’.

15
16 CHAPTER 2. PROTECTION ZONES I – DISTANCE CRITERIA

Figure 2.1: Well dialog box.

Figure 2.2: The Other tab.

Click on the Other tab, and fill-in a setback radius of 200 feet. See Figure 2.2.
Click OK to add the well to the database. Zoom into the well to see the dotted setback circle. Turn off
the setback radius for this pumping center. Now create the four other wells with their 200 f t setbacks.
Draw the wellhead protection area (Edit , Draw wellhead protection area). Name the protection
area “200 f t setback”. With the left mouse button, click a series of short segments along the outside envelope
of the dashed setback circles. The closing of the wellhead protection area is accomplished with a right-mouse
button click, and the wellhead protection area will show as a blue line. See Figure 2.3.
The purpose of drawing this area by hand is to separate the delineation of the capture zone (based on
hydraulics/hydrology) from the definition of the wellhead protection area (based on other considerations,
such as parcel or zoning boundaries). The draw utility allows the wellhead manager to fine tune the shape of
the wellhead protection area to correspond with on-the-ground realities. And as in this case, when multiple
capture zones are generated using various conceptual models, a “wellhead protection area” can be drawn
2.1. SETBACKS 17

Figure 2.3: 200 foot setback envelope.

that envelops the various capture zones for maximum protection. More on this later.
Continue to the next chapter or exit the program by clicking Project and Exit.
18 CHAPTER 2. PROTECTION ZONES I – DISTANCE CRITERIA

This page is blank.

Chapter 3

PROTECTION ZONES II —
RESIDENCE TIME CRITERIA

The justification of the residence time criteria as protective of subsurface source water is based on the
assumption that: (1) non-conservative contaminants (open to sorption or transformation processes) will
have an opportunity to be assimilated after a given mean time in the subsurface; or (2) that detection of
conservative contaminants (no sorption or transformation) entering the protection area will give enough lead
time for the drinking water community to develop a new water supply or take immediate remedial action.
In this chapter we will examine two simple methods involving the residence time criteria as applied to our
case study: (1) calculated fixed radius; and (2) well in uniform flow.

3.1 Calculated Fixed Radius

The fixed radius is calculated based on a simple, two dimensional, static water balance analysis, assuming
negligible ambient flow in the aquifer compared to the pumping rate of the well. If we assume radial flow
toward a well in an aquifer with a constant saturated thickness H(m), the cylindrical boundary of radius
R(m) is delineated by an isochrone of residence time t(days). This means that any water particle that enters
the cylinder or is present in the cylinder will travel no longer than t days before being pumped-up by the well
(See Figure 3.1). If the pumping rate of the well is Q(m3 /day), the areal recharge to the water table rate
due to precipitation is N (m/day), and the aquifer porosity is n(-), a water balance for the period t yields:

1
N πR2 t + nπR2 H ≈ Qt (3.1)
2

The first term represents the inflow due to aquifer recharge. Since the area ranges from 0 to πR2 during
the time span t only half of the area is included (the average area during time span t). The second term
represents the amount of water contained inside the cylindrical aquifer, and the term of the right hand side
is the total amount of water removed by the well. The radius R can, therefore, be expressed as:
s
Qt
R≈ 1 (3.2)
2 N πt + nπH

This approximate equation is adequate up to a radius R that is about 80% of the radius Rdivide of the total
capture zone envelope associated with the water divide [Haitjema, 2011]. The latter can be obtained from a

19
20 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

Figure 3.1: Water balance for radial flow to a well in a domain bounded by an isochrone of residence time t.

simple water balance

2
tπRdivide N = Qt (3.3)
so that s
Q
Rdivide = (3.4)
πN
Use of this equation is called the “recharge method” [USEPA,1993].
If t or N are relatively small the recharge term in equations 3.1 and 3.2 may be ignored:
s
Qt
R≈ (3.5)
nπH

Use of this equation is called the “volumetric method” [USEPA,1993]. Equation 3.5 therefore represents an
overestimation of the radius R for the isochrone, which is conservative for the purpose of wellhead protection.
The advantage of this simplification is that no determination of the (local) recharge is required: a parameter
that is often difficult to obtain.
The equations 3.2 through 3.5 only apply if the Dupuit assumption applies, that is for shallow, horizontal
flow. If the capture zone or isochrone radius is smaller than say twice the saturated aquifer thickness,
and the well is partially penetrating, then three-dimensional effects might become important and these
approximations can lead to an underestimation of the capture zone isochrone radius R [Haitjema,2006].
Such a non-conservative result is avoided by replacing the saturated thickness H in Equation 3.5 by the
length of the well screen.
The radially symmetric isochrone in Figure 3.1 will only occur in the absence of ambient flow or if the
well dominates the flow inside the isochrone. For an isochrone with t as long as 1 or 2 years and a high
capacity well, the assumption of radial flow is often reasonable. However, a low capacity well in a strong
ambient flow field will exhibit an elongated isochrone that extends in the upgradient direction well beyond
the circular representation of the isochrone computed by using Equation 3.2. Under these conditions, the
direction and magnitude of the ambient flow must be determined and isochrones should be constructed using
the solution for a well in a uniform flow field.
3.2. WELL IN UNIFORM FLOW 21

The isochrone calculations, whether with or without ambient flow, are based on average aquifer properties.
In reality, an aquifer is often stratified exhibiting different layers with different hydraulic conductivities
and porosities. An effective constant saturated thickness H might occur in a confined aquifer of constant
thickness. The saturated thickness of an unconfined aquifer will not be constant, even assuming a horizontal
base, due to drawdown of the water table. A conservative (protective) capture zone will be obtained by using
the smallest saturated thickness as H in Equation 3.5.

3.1.1 Exercise: Calculated fixed radius for city wellfield

Step 1 Use WhAEM to plot the Calculated Fixed Radius protection areas based on the recharge method
and the volumetric method.

We will continue to build on the previous exercise. Open the setback1.whm project. You may get
a warning about no specified head condition set yet. Acknowledge, and carry on. From the menu, select
Project then Make Duplicate Database and type vincfr1 in the File Name box. Click open.
We cannot ignore the effects of the mutual interferences of the multiple pumping wells. So for practical
purposes we will define a single pumping center point with a pumping rate the same as the sum of the
multiple pumping wells. We will assume the pumping center point of the wellfield is located at x = 452650 m,
y = 4280665 m, based on our previous exercise.
For the recharge method, use equation 3.4. Assume the well field pumps Q = 370, 000 f t3 /d , the areal
recharge to the water table N = 0.00276 f t/d (12.1 in/yr), and calculate R. If you do not have a portable
calculator handy, we suggest you use the Microsoft Windows calculator (Start > Calculator). Recall π = 3.14.
To draw the new radius on the screen, double click on the well, and edit the setback box and type in the new
radius. You do not need to fill out the other information for the well at this time. After clicking OK, drop
down the View menu and click on Refresh. You may have to zoom out to see the circle (setback zone).
Draw the wellhead protection area (recall the exercise in the previous chapter).
For the volumetric method, use equation 3.5. Assume Q = 370, 000 f t3 /d, aquifer saturated thickness
H = 67.5 f t, the aquifer porosity n = 0.2, and the pumping period is five years (1826.25 d). Draw the
wellhead protection area.
It may be useful to edit the vertices that make up the wellhead protection area envelope. Select a vertex
by right-clicking on a node; it will turn red. The options to move, delete, or refine can be accessed through
the Edit menu, or the smart icons.
Save these protection areas for future overlay. This is a two step process. First, export to DXF: Tools >
Export > DXF > Wellhead protection areas, accept the default layer name, and name the saved dxf
file. You can then bring into the basemap through Project > Project Settings > Add map.
Compare your results to Figure 3.2.
It is possible to toggle basemaps on and off through the commands View > Base Maps > Raster
graphics (also for Vector graphic basemaps).

3.2 Well in Uniform Flow

Ambient flow, resulting from far field aquifer recharge due to precipitation and ground water exchanges with
streams and lakes, is approximated by a uniform flow field (straight streamlines) in the immediate vicinity of
the well. The capture zone for the well in a uniform flow field will no longer be circular and centered about
the well, but will be a somewhat elongated oval-shaped domain in the direction of the uniform flow.
22 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

Figure 3.2: Calculated fixed radii for city wellfield based on recharge (R=6532 feet) and volumetric methods
(R=3992 feet).
3.2. WELL IN UNIFORM FLOW 23

To determine this capture zone (e.g. by use of WhAEM ), it is necessary to determine:

• direction of the ambient flow.

• the hydraulic gradient.

• the aquifer transmissivity.

• the pumping rate of the well.

• the desired maximum residence time.

Given three water levels in three wells (average, static), it is possible to estimate the direction of the
ambient flow and the hydraulic gradient [Heath,1983](See Figure 3.3):

a. Identify the well that has the intermediate water level.

b. Calculate the position between the well having the highest head and the well having the lowest head
at which the head is the same as that in the intermediate well.

c. Draw a straight line between the intermediate well and the point identified in step (b). This line
represents a segment of the water level contour along which the total head is the same as that in the
intermediate well.

d. Draw a line perpendicular to the water level contour and through the well with the lowest (or highest)
head. This indicates the direction of ground water movement in an isotropic aquifer.

e. Divide the difference between the head of the well in step (d.) and that of the contour by the distance
between the well and the contour. This gives the hydraulic gradient.

Alternatively, the reader is referred to EPA’s on-line three-point gradient calculator created by Dr. Jim
Weaver at (www3.epa.gov/ceampubl/learn2model/part-two/onsite/gradient3ns.html).
A 3-point method based on the squared heads (water levels) has been shown to reduce bias in the estimate
of the magnitude and direction of the uniform flow vector in heterogeneous unconfined aquifers, provided
the monitoring wells have a significant separation distance, and that there are no sources or sinks in between
[Cole and Silliman, 1996].
The scientific literature, the US Geological Survey, and state Geological Surveys are good resources for
studies of the regional aquifer system and the regional hydraulic gradient.
Aquifer transmissivity is measured by a pumping test. This is an expensive procedure, and transmissivity
is known to be scale dependent. Perhaps transmissivity measurements are available from locations not too
far from the well. Estimates of hydraulic conductivity based on geology (e.g., rock type, unconsolidated
media type, grain size) are available in text books. The USGS Ground Water Atlas is a good resource for
regional information about the aquifer system (water.usgs.gov/ogw/aquifer/atlas.html).
The well-in-uniform-flow solution is difficult to parameterize in practice. It is difficult to anticipate the
average hydraulic gradient and direction of flow given limited synoptic observations of a transient phenome-
non. It is also difficult to separate out the local effects of the well, unless it is out of production for a period
of time. For these reasons, caution is advised in the use of the well-in-uniform flow solution. The modeler is
encouraged to respond to the uncertainty by generating a number of reasonable solutions, and encapsulating
the union of solutions into a protection area. In other words, use the “wellhead area” drawing tool to draw
an envelope around the set of calculated capture zone realizations.
24 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

Figure 3.3: Procedure A for determination of an equipotential contour and direction of ground water flow in
homogeneous,isotropic aquifer (after [Heath,1983]).

3.2.1 Exercise: Well in Uniform Flow for city wellfield

Step 1 Determine the magnitude and direction of a uniform flow field near the city pumping center point.

First we will give this exercise a new name. Open the database from the previous exercise vincfr1. From
the Menu, Project , Make Duplicate Database, and enter vinuni1.
As a first estimate, we will attempt to characterize the uniform flow field by using the information
contained in the USGS DRG topographic map. The “lakes and ponds” in this area are actually abandoned
quarries, and the removed overburden gives a window into the elevation of the water table. We assume these
pits act like large piezometers, and field observations confirm that their water levels do in fact respond to
recharge events.
A note of caution is in order regarding extracting water elevations from USGS topographic maps. These
elevations reflect a snapshot in time (check the date of publication of the map), and may not reflect conditions
related to your study. The USGS map was photo-revised in 1989. The elevation values are accurate to a
value about one half the distance between contours; in this case 2.5 feet for areas represented with 5 foot
contour lines, and 5 feet for 10 foot contour lines [USGS, 1999]. The spatial accuracy of the topographic
lines are accurate approximately plus or minus 40 feet [USGS, 1999]. These uncertainties are acceptable for
this level of analysis.
Use the Heath method or the OnSite Calculator discussed in the previous section with the following three
points of known head: (1) Mirror Lake (407 f t); (2) Otter Pond (403 f t); and (3) the 400 f t topographic line
crossing of the Wabash River.

Step 2 Use the uniform flow element in WhAEM to fit the estimated uniform flow.
3.2. WELL IN UNIFORM FLOW 25

Figure 3.4: WhAEM view of test points used to calculate the uniform flow vector.
26 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

Figure 3.5: Aquifer settings dialog box.

Zoom to approximately the view shown in Figure 3.4, which includes Mirror Lake and Otter Pond.
From the Menu, select Element , then Uniform Flow. Acknowledge the warning. Click the cursor over
Mirror Lake. Enter in the Reference Head 407 f t. Enter the gradient (for our case about 0.001). Enter the
orientation (for our case about 150 degrees, which can be estimated based on the gradient calculator, the 3
point method, or trial and error).

Step 3 Enter the aquifer and contouring properties and test points.

From the Menu, select Model , then Settings, and then select the Aquifer tab. Set the base elevation
to 330 f t, the thickness to 100 f t, the hydraulic conductivity to 350 f t/d, and leave the porosity as 0.2.
(The geology of the study area is described in [Shedlock, 1980].) See Figure 3.5.
Now set the Contouring properties. Select the Contouring tab. Check the Show Contours box and
select the radial button for coarse (resolution). Enter the minimum contour elevation of 390 f t, the
maximum elevation of 430 f t, and a contour interval of 1 f t. See Figure 3.6.
Next we set some options for the solution process. Select the Solver tab. Make sure that the options:
View error log file after solve and View message log file after solve are both checked. While our uniform flow
model has only linear equations (see comment on the dialog box), you can set the number of iterations to 2
in order to achieve some “iterative refinement” of the solution. See Figure 3.7.
Click on OK.
To test our conceptual solution, we can introduce so-called “test points” at locations of known head, e.g.
Otter Pond (403f t and utm 452325, 4277311) and the 400 f t topographic intersection (utm 454815, 4283085)
with the Wabash River, by selecting Model , then add test points. Locate each point on the screen with
a click of the left mouse button, and enter the elevations. The test points will store the difference between
the model head and the observed head. Note: the head differences at the test points are visually displayed
after the “solve” with a scaled triangle, or can be queried by clicking on the test point (look at the status
bar at the bottom of the screen when clicking on a test point). The solve process is described below.
3.2. WELL IN UNIFORM FLOW 27

Figure 3.6: Contouring settings dialog box.

Figure 3.7: Solver settings dialog box.

28 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

Figure 3.8: Solve icons.

Figure 3.9: Hydraulic head contours of uniform flow field with test points.

Step 4 Solve the model and plot the contour map.

From the menu select Model , then Solve (Alternatively, use the F9 mapped function key, or the smart
icon, See Figure 3.8). A black box appears with progress indicators for the solution process. When the
solver is done, the base map reappears with an overlay of piezometric contours (blue lines) that are labeled
with their piezometric head elevations. On top of the base map the two log files, specified on the Solver
tab, are opened in Notepad. The file Error.log should have only an asterix on the start of the first line and
be empty otherwise. This file is used by the Solver to report problems (errors or warnings) with the input
data. Close the file after reviewing its contents. The file Message.log contains Solver activity messages and
a report of maximum errors in boundary conditions. In general, errors should be less than 1%. Close the
file after reviewing its contents.
The hydraulic head contour map of this solution is shown in Figure 3.9.
Note: since our model heads will be close to the test point heads, the triangles designating the head difference will
3.2. WELL IN UNIFORM FLOW 29

Figure 3.10: Well dialog box.

be too small to see. Change the direction or magnitude of the uniform flow field and you will see the triangles grow in
size. You can modify the test point appearance by selecting Test Points on the View menu.

Step 5 Enter the pumping rate and radius of the pumping center point.

Locate the city wellfield pumping center point either by double clicking on the existing well or creating a
new one (Element > New > Well). Enter in the well pumping rate of 370, 000 f t3 /d. Enter a radius of
4 f t to represent our pumping center. The city wellfield actually consists of a local cluster of at least 5 wells
that we will represent with the single pumping center. The entered value for the radius in WhAEM does
not effect the flow solution, but will provide the starting point for tracing particles during capture zone
delineation. Too small a radius might effect the accuracy of the particle tracing. See Figure 3.10.
Enter OK. You may get a box asking if you want to remove the current solutions results from the view
and invalidate the solution. Respond in the affirmative.

Step 6 Plot a capture zone based on 5 years of residence time.

On the well properties box, click on the “Other” tab. Then check the box Trace particles from Well
in the Other tab, and set the Trace Number of particles to 20. The number of pathlines will control the
resolution of the capture zone. For exploration the user may choose a low number, say 8 pathlines. For final
maps, the user may choose a larger number, say 30 pathlines. The well is assumed fully penetrating so we
have an option to set the starting elevation of the pathlines. Set the Starting elevation to 335 f e (that
is 5 f t above the aquifer base). Click on OK. Go back into the Model , Settings, and select the Tracing
tab. Check the box Compute particle paths. Accept the Default Step size, and change the Maximum
travel time to 1826.25 d(5yr). In the “Contouring” tab, click off the show contours box. Click OK and
click the calculator button on the tool bar (re-solve). You can add the tic marks by going to View , select
Pathlines, then Show Time Tics, and check Every 1 years. Draw the wellhead protection area around
the trace lines. The graphics should look similar to Figure 3.11. You are encouraged to try a number of
equally valid uniform flow field solutions, and vary the Model setting parameters to get insight into the
sensitivity of the solution to the parameterizations.
30 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

Figure 3.11: Capture zone for 5 years residence time based on tracelines for well in uniform flow field. The
previous CFR solutions, volumetric and recharge, are also shown for scale.
3.2. WELL IN UNIFORM FLOW 31

Figure 3.12: The 5 yr residence time tracelines are shown for the capture zone for a well in uniform flow
(with the gradient set zero). The volumetric CFR and recharge CFR are shown for comparison.

Save the protection areas for future overlay. First, export to shapefile: Tools , Export, Shapefiles,
WHP Areas, and name the file(s).
The solution is becoming more realistic. However, we still do not account for the presence of the Wabash
River and other hydrological features in the area! We shall do that in the next chapter.

Step 7 Compare solutions.

It is interesting to do a verification of the well-in-uniform flow solution in comparison to the CFR

volumetric method. Click on the well, enter the 5 year CFR volumetric method result previously calculated
(R = 3992 f t). Change the uniform flow gradient to a very small number, say 0.00001, and plot the 5 yr
capture zone. Notice that the uniform flow and the CFR capture zones become essentially equivalent, as
expected. See Figure 3.12.
32 CHAPTER 3. PROTECTION ZONES II — RESIDENCE TIME CRITERIA

This page is blank.

Chapter 4

SIMPLE WHPAs

4.1 Introduction
The material presented in this section is based on the work of Ceric [2000], Haitjema [2002], and Ceric and
Haitjema [2005] sponsored by the Indiana Department of Environmental Management.
A “Simple WHPA” is a wellhead protection area approximation based on an analytic solution for the cap-
ture zone of a pumping well. The USEPA recognizes several simple approximations for capture zones, among
them the arbitrary fixed radius capture zone and the calculated fixed radius capture zone. Ceric introduced
an improved calculated fixed radius approximation based on a dimensionless travel time parameter. While
the Simple WHPA can be calculated and graphed by hand [Haitjema, 2006], both the arbitrary fixed radius
and the improved calculated radius method are included in the Simple WHPAs option on the Tools menu
of WhAEM or as smart icons on the toolbar. In addition, some of the functionalities from the USEPA
computer program WHPA version 2.2 are reproduced under this option. Specifically, the Simple WHPA op-
tion might require a well-in-uniform-flow field solution. The user can optionally include a stream boundary
(infinitely long straight line) near the well. Finally, as in WHPA 2.2, capture zones can be delineated using
multiple wells in a uniform flow field near a (infinitely long straight) stream boundary. In all cases, the
stream boundary is modeled using image recharge wells (as is done in WHPA 2.2). The no-flow boundary
option in WHPA 2.2 is not reproduced as part of the Simple WHPAs feature, since it leads to internally
inconsistent solutions (this was also a problem in WHPA 2.2).
When selecting Simple WHPAs, three options are offered:

1. Arbitrary Radius. Draws a circle with a specified radius around selected well(s).

2. Calculated WHPA. Draws centric or eccentric radius around well, or draws streamlines for a well in a
uniform flow field, depending on the dimensionless time of travel parameter T̃ , to be discussed below.
In the case of a well in a uniform flow field a stream boundary can be specified.

3. Multiple Wells. To be used whenever calculated capture zones are to be generated for more than one
well. Allows for nearby wells (well interference).

The Simple WHPAs options will be displayed if the Show Simple WHPAs is checked on the View menu.
The Simple WHPAs can be permanently deleted using Edit , Delete All, Simple WHPAs.

33
34 CHAPTER 4. SIMPLE WHPAS

Figure 4.1: The magnitude of the discharge vector of the uniform flow field is calculated knowing the hydraulic
conductivity (k), the hydraulic gradient (i), and the saturated thickness (H). The simplified case is for a
well pumping at rate Q/ The maximum width of the capture zone envelope is w. The cross-sectional area A
is w times H. The hydraulic head is φ.

4.2 Back-of-the-Envelope WHPAs

How does the Simple WHPAs tool work? The calculations depend on several parameters, including the
magnitude and direction of the ambient flow near the well or well field, which is challenging to characterize.
The magnitude of the uniform flow is denoted by Qo [L2 /T ], and can be estimated from the hydraulic
gradient i [-] and the aquifer transmissivity kH (hydraulic conductivity k times saturated aquifer thickness
H) [L2 /T ]. See Figure 4.1. The magnitude of the uniform flow rate is calculated as

Qo = kiH (4.1)

The flow Qo is the total amount of water in the aquifer integrated over the saturated thickness, per unit
width of the aquifer.
The ambient flow is challenging to estimate, and the reader is referred to the previous discussion of
Section 3.2.
The shape and size of a simplified time-of-travel capture zone can be related to a dimensionless travel
time parameter, T̃ , defined as
T
T̃ = (4.2)
To
where T is the time-of-travel [T] and To [T] is a reference time defined as:
nHQ
To = (4.3)
2πQo 2
4.2. BACK-OF-THE-ENVELOPE WHPAS 35

where n is the aquifer porosity [-], and Q [L3 /T] is the pumping rate of the well.
When T̃ ≤ 0.1, the calculated fixed radius (CFR) centered on the well, including a safety factor for a
non-zero ambient flow field, is given by s
QT
R = 1.1543 (4.4)
πHn
When 0.1 < T̃ ≤ 1, the CFR is given by

R = Ls [1.161 + ln(0.39 + T̃ )] (4.5)

where Ls is the distance from the well to the stagnation point down gradient from the well given by

Ls = Q/(2πQo ) (4.6)

and where the eccentricity δ is given by

δ = Ls [0.00278 + 0.652T̃ ] (4.7)

When T̃ > 1, a uniform flow envelope, the so-called boat-shaped capture zone, can be defined as

x = y/tan(y/Ls ) (4.8)

where y is bounded by
−Q/2Qo < y < +Q/2Qo (4.9)
and clipped at the up-gradient distance Lu given by

Lu = Ls [T̃ + ln(e + T̃ )] (4.10)

and where e = 2.718.

The reader is referred to Ceric [2000], Ceric and Haitjema [2005] for details on the above analysis. This
calculation sequence and results are illustrated Figure 4.2.
We will explore the Simple WHPAs tool based on our knowledge of the calculated fixed radius and
well-in-uniform-flow for an example case study using WhAEM .

4.2.1 Exercise: Simple WHPAs and no ambient flow

Step 1 Start WhAEM . Create a new database called simple cfr.whm. Set the computational units to
f t and d, and add the base map vincennes 24k.tif (located in your My Documents, WhAEM projects,
Vincennes folder).
Step 2 Add the wellfield. Zoom in and add the single pumping center well (utm-x=452650, utm-
y=4280665, discharge=370000 f t3 /d, radius 4 f t).
Step 3 Calculate the Simple WHPA with no ambient flow. Click on Tools and then Simple WHPAs
then Calculated WHPAs or use the smart icon on the toolbar. Fill in the dialog as shown in Figure 4.3.
Continue, note that the T̃ value is less than 0.1, and Draw the Simple WHPA based on a calculated fixed
radius.
The 5-yr Simple WHPA is slightly bigger than the traditional 5-yr volumetric-calculated fixed radius of
Equation 3.5 due to the safety factor in Equation 4.4, as shown in Figure 4.4. In order to create the map
with multiple CFR solutions you can turn on and off previously exported DXF; to add a wellhead protection
area to the base map (Project > Project Settings > Add Map); to select the wellhead protection area
for display (View > Base Maps > Vector Graphics).
36 CHAPTER 4. SIMPLE WHPAS

Figure 4.2: Simplified delineation techniques for a well pumping at rate Q, in an ambient flow field Qo ,
with aquifer saturated thickness H and porosity n, and time of travel capture zones of time T , after [Hai-
tjema,2002].

Figure 4.3: Dialog box for calculated WHPA.

4.2. BACK-OF-THE-ENVELOPE WHPAS 37

Figure 4.4: 5-yr Simple WHPA with no ambient flow (red circle) compared to calculated fixed radius solutions
using volumetric equation and recharge equation.
38 CHAPTER 4. SIMPLE WHPAS

Figure 4.5: Dialog box calculated WHPA including ambient flow.

4.2.2 Exercise: Simple WHPAs with ambient flow

Step 1 Create a new database called simple uf.whm. Continue working in ft and days for computational
units, and add the base map vincennes 24k.tif. This can also be done by Project > Make a Duplicate
Database called simple uf.whm and Edit > Delete All > Simple WHPAs > Calculated WHPA.
Step 2 Set the model settings, such as aquifer properties (aquifer base at 330 f t, thickness at 1000 f t,
hydraulic conductivity at 350 f t/d, and porosity at 0.2), contouring (390 to 400 step 1 f t), and tracing
(1825 d).
Step 3 Add the wellfield as before (single pumping center well: utm − x = 452650, utm − y = 4280665,
discharge=370000 f t3 /d, radius 4 f t). Add other properties, such as 20 tracking particles starting from
elevation 335 f t. The choice of elevation is a bit arbitrary at this point; just need to make sure below the
water table elevation and above the aquifer base elevation. Choose the color of the tracelines such as red.
Step 4 Calculate the Simple WHPA with ambient flow. Fill in the dialog box shown in Figure 4.5 where
ambient flow is calculated using Equation 4.1. Fill in the next dialog box with a direction of flow of 150
degrees, and as shown in Figure 4.6. Note the T̃ value is greater than 1. You may need to modify the window
and zoom out and re-solve in order to view the full capture zone envelope: Edit > Delete All > Simple
WHPAs > Calculated WHPA; zoom-out; then repeat exercise.
Click the no response to the issue of a stream boundary, see Figure 4.7. and Draw.
The slight difference in the solutions for the Simple WHPA and the well-in-uniform-flow function is that
the former is solved under the assumption of constant transmissivity, while the latter allows for unconfined
flow conditions. See Figure 4.8.

4.2.3 Exercise: Simple WHPAs with stream boundary

The Simple WHPA tool can also be used to create capture zones for a well field including the influence of
a stream boundary. The reader should follow the instructions in the dialog boxes to create something like
4.2. BACK-OF-THE-ENVELOPE WHPAS 39

Figure 4.6: Dialog box simple well in uniform flow.

Figure 4.7: Dialog box simple stream boundary.

40 CHAPTER 4. SIMPLE WHPAS

Figure 4.8: Simple WHPA including ambient flow tracelines compared to well in uniform flow envelope.
4.2. BACK-OF-THE-ENVELOPE WHPAS 41

Figure 4.9: Simple WHPA including the influence of a stream boundary in comparison to capture zone for
a well-in-uniform-flow field.

Figure 4.9. The elevation of the Wabash River is around 400 f t near the wellfield. Notice that the 5-yr
capture zone has shrunk in size since the well is now getting water from the river.

4.2.4 Exercise: Simple WHPAs with multiple wells

The Simple WHPA tool can even be used to create capture zones for multiple wells within an ambient flow
field. The reader should follow the instructions in the dialog boxes (see Figure 4.10) to create something like
Figure 4.11. Notice that the trace lines from each well can be assigned a different color.
42 CHAPTER 4. SIMPLE WHPAS

Figure 4.10: Dialog box Simple WHPAs with multiple wells.

Figure 4.11: Simple WHPAs including multiple wells in ambient flow.

Chapter 5

PROTECTION ZONES III — FLOW

BOUNDARIES CRITERIA

In the previous analyses we de-emphasized the presence of hydrological boundaries. Instead, we lumped
their effects near the well into a uniform flow field. Or we represented rivers as infinite linear features. Our
geology might also have been oversimplified. In this chapter we will examine the influence of hydrological
and geological boundaries on the capture zone. To do so we will use the full features of WhAEM . We
will progressively refine the complexity of our hydrogeological conceptual model using the glacial outwash
example. But before doing so we will explain briefly how the analytic element methods work. The reader
is encouraged to explore in more depth the topic of system and boundary conceptualization outside of this
tutorial [Reilly, 2001]. The sufficient representation of geohydrologic complexity in capture zone delineation
for the purpose of wellhead protection is summarized in the following chapter.

5.1 Introduction to the Analytic Element Method

WhAEM uses the analytical element method to solve ground-water flow subject to hydrological boundary
conditions. Analytic elements are mathematical functions that represent geohydrologic features in a ground-
water flow model. For instance, a meandering creek or river might be represented by a string of straight line
elements that collectively approximate the geometry of the stream. In WhAEM each of these line elements is
associated with a “line-sink function” that represents mathematically the ground-water inflow or outflow for
that stream section. The line sink can be imagined as an infinite number of point sinks, or wells, integrated
along a finite line. The sink density of the line-sink equals the ground-water inflow/outflow into/out of the
stream section. Each stream section is likely to have a different ground-water inflow rate and thus a different
sink density. Sink densities are not known in advance. What are considered known are the heads underneath
the stream sections. For instance, in WhAEM it is assumed that the head underneath the center of each
line-sink is equal to the average water level in the stream section it represents. WhAEM has as many known
heads at line-sink centers as it has unknown sink densities for the line-sinks. The unknown sink densities can
then be calculated by solving a system of (linear) algebraic equations, see [Strack,1989] and [Haitjema,1995].
Other examples of analytic elements in WhAEM are wells (point sinks), line doublets or “double layers”
to represent aquifer inhomogeneities, and area elements to model area recharge due to precipitation (using
a negative sink density to create infiltration). Once all sink densities and doublet densities are calculated —
the solution procedure in WhAEM — the head and average ground-water flow rate can be calculated at any

43
44 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

point by adding the influence of all the line-sinks, point-sinks, line doublets, and areal elements. This is the
principle of superposition.
Note: the ground-water velocity is calculated by adding “derivative functions” for the analytic elements. These
“functions” have been obtained analytically and are programmed into WhAEM . For details on the mathematical
“functions” and on the mathematical procedures to enable the superposition of solutions for confined and unconfined
flow, the reader is referred to [Strack,1989] and [Haitjema,1995].
Analytic elements are defined in WhAEM by clicking the Element menu, and selecting New.

5.2 Rivers and Line-sinks

Line-sinks are used to represent streams (creeks, ditches, and rivers) and lake and coastal boundaries. A line-
sink in WhAEM has a constant sink density (constant groundwater extraction rate along the line element).
The sink density or “strength” of a line-sink that represents a stream section or a section of a lake boundary
is in general not known. In WhAEM this is handled by specifying a piezometric head at the center of the
line-sink, that is selected equal to the average water level in the stream section or lake that is represented by
the line-sink. For every unknown line-sink strength there exists one specified head. This leads to a system of
equations that is solved for the unknown strength parameters of the line-sinks. Every change in data in the
WhAEM model, therefore, requires a new solution procedure to recalculate the strength parameters. The
line-sink strengths represent the ground-water extraction rates of the line-sinks, in volume per time per unit
length of the line-sink (stream or lake boundary).
The use of line-sinks along lake boundaries implies that in WhAEM a lake receives or loses ground water
along its boundary only. This appears to be a good approximation of reality when the size of the lake is
large compared to the saturated thickness of the aquifer. For rivers and streams that have a width that is
large compared to the aquifer thickness, the realism of the line-sink representation might be improved by
positioning line-sinks along both sides of the river.
The number of line-sinks required to represent a stream or lake boundary depends on two factors: (1)
the need to follow the detailed stream geometry, and (2) the need to allow for variations in the ground-water
discharge or recharge along the stream or lake boundary. Representing stream geometry is more important
in the near field (immediate area of interest) than in the far field. Similarly, the need to use many small
line-sinks to follow discharge or recharge variations along the stream is most pressing in the near field,
particularly near a high capacity well drawing water from the stream. Consequently, surface water features
in the area of interest (near field) should be represented with more and smaller line-sinks than surface waters
in the far field, remote from the area of interest.
Line-sinks are defined in WhAEM by selecting the Element menu, and selecting New and then Line-
sinks.

5.2.1 Exercise: Study area and Line-sinks

Step 1 Locate the pumping center.

Create a new project database called hydro1 with basemap units meters and computational units f t and d.
Add the regional series of USGS 7.5 minute raster basemaps
Vincennes 24K utm16n27.tif,Bicknell 24K utm16n27.tif, Decker 24K utm16n27.tif,
EMtCarmel 24K utm16n27.tif, Frichton 24K utm16n27.tif, Heathsville 24K utm16n27.tif,
Iona 24K utm16n27.tif, Lawrenceville 24K utm16n27.tif, MonroeCity 24K utm16n27.tif,
5.2. RIVERS AND LINE-SINKS 45

Oakton 24K utm16n27.tif, Patoka 24K utm16n27.tif, Plainville 24K utm16n27.tif,

Russellville 24K utm16n27.tif, SandyHook 24K utm16n27.tif, StFrancisville 24K utm16n27.tif,
Washington 24K utm16n27.tif, Wheatland 24K utm16n27.tif to the project, and view.
Position the pumping center as described earlier, and create the discharge specified well (Qw = 370, 000 f t3 /d,
radius = 4 f t, utm-x= 452650 m, utm-y= 4280665 m ).
From the Menu, select Model , and reset your Settings, including Aquifer and Contouring parameters
as in previous exercises (base elevation 330 f t, hydraulic conductivity 350 f t/d, porosity 0.2, thickness 500 f t).
Step 2 Annotate the basemap with elevations.

It is useful to annotate your basemap with estimated surface water elevations at those streams that you
intend to include in the model (represent by line-sinks).
Zoom In to the location of the 400 f t contour north of the city. From the Menu, select Edit, and then
Add text. Move the text cursor to the proper location, say the center of the river and along the 400 f t
contour and click the left mouse button. Type the contour elevation 400 f t in the dialog box, check the
Hydrography Label option, and click OK. The label appears in the map and the text cursor reappears.
You need to repeat this procedure to add all the rest of the hydrography elevations, i.e. elevations where
topography contour lines cross the perennial surface water features.
It is often sufficient to add hydrography labels near the map boundaries and near confluences of streams.
When creating line-sink strings you will only be prompted for a head at the start of the string and at the end
of the string. It is good practice to start a linesink string at the higher elevation and work down. This will
ensure more meaningful data reporting in the graphical user interface (GUI). The WhAEM program linearly
interpolates between the starting head and ending head for the heads at vertices of the linesink stream.
We have created a text file with elevations to assist in the completion of your basemap. The label file is
a comma delimited ASCII text file (filename.TX) in which each line has the following format:

utm − x, utm − y, h, label (5.1)

where x and y are the (basemap) coordinates of the lower left corner of the text, where h = 1 for hydrography
labels (water levels in surface waters) and h = 0 for any other text label. The label is the water level or text
to be printed on the map. The text labels will be added to the current project database file. To import the
file containing text labels, select the Tools menu, then Import, and then Text Label File. For this case
study, select the file Hydro Example.tx. In order to view the hydrography labels on the basemap, you may
need to toggle them on (View > show hydrography labels).
See Figure 5.1.
It is an advanced procedure to change the size of the text labels on the basemap. From Tools select
WhAEM Database Viewer. Acknowledge the warning message. Then View > Base Tables > Text
Labels. Edit the cells and set the desired Font Size. The labels in 5.1 are 18 point. You can also change
the color of the labels.
Step 3 Represent the River with a string of Line-sinks.

At this point, it is recommended that you set your window to frame the region including the 410 f t
elevation and 400 f t elevation of the Wabash River. Then, from the menu select Element and then New,
and then Linesink, or alternatively, click on the Line-sink smart icon (See Figure 5.2). You will be prompted
for the name of the line-sink string about to be entered. Let’s start with Wabash River, hence type WabashR01.
Treat the linesinks as “far field” and the click the “along stream centerline” radio button. Next fill in the
Starting Head of 410 f t and Ending Head of 400 f t (See Figure 5.3). Click OK. Move the cursor to
46 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.1: Regional USGS DRG basemaps with hydro elevation labels.

Figure 5.2: Elements smart icons: wells, linesinks, inhomogeneities.

5.2. RIVERS AND LINE-SINKS 47

Figure 5.3: Linesink properties dialog box.

center of the Wabash at the 410 f t text label north of the city and click the left mouse button. Now move
the cursor down stream, clicking the left mouse button to enter linesink vertices. After clicking on the 400 f t
location, click the right mouse button to end the linesink entry. Click OK.
Note: As mentioned, it is recommended to enter line-sink strings along streams going down-stream (following the
direction of the stream flow). This is particularly important when querying line-sinks that are part of stream networks,
as it will assure that the stream flow reported actually occurs at the vertex that is highlighted.
That completes your first line-sink string. At the center of each line-sink in the string, a head will be
specified by WhAEM based on a linear interpolation of the starting head and ending head along the line-sink
string.
Continue representing the Wabash River with linesinks. Represent the south/east side of the river (click
the “along surface water boundary” radial) with a linesink string spanning the 400 f t and 397.5 f t elevations
(estimated half way between the 400 f t and 395 f t stream elevations.
We also recommend you represent the north/west side of the Wabash River with a string of line-sinks,
particularly near the wellfield. Start again at the of 400 f t contour, and continue to the estimated location
of 397.5 f t.
Return to Wabash River centerline linesinks from 397.5 f t to 395 f t and between 395 f t and 390 f t, and
between 390 f t and 380 f t.
You may change the location or properties of line-sinks as follows. Click on the line-sink vertex in the
string you want to edit. The well or line-sink string will become highlighted (red). Click on the Edit menu
and select Properties, Move, Delete, or Refine. For line-sinks you will be asked whether you want to
delete only the selected vertex or the entire string. When selecting Refine, a vertex is added midway in the
line-sink on either side of the vertex that was selected (highlighted in red).
Let’s refine the line-sink distribution of the Wabash River near the wellfield. This is the area of induced
recharge from the river to the well, and the solution will be improved with refinement. Click on a vertex of
the line-sink string. Next click on Refine on the Edit menu. You may also click on the refine string smart
icon (the third to the right of the printer icon). See Figure 5.4. Notice the two added vertices (dots) on either
48 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.4: Smart icons for editing vertices: delete, move, refine.

side of the vertex you clicked on. Repeat this for one or two other vertices of the Wabash River line-sink
strings near the well and on both sides of the river.

Step 4 Surround the wellfield with hydrologic (line-sink) boundaries.

Create line-sinks representation of the major perennial streams. You should attempt to locate head-
specified boundaries (line-sinks) in all directions projected out from the well-field, putting more detail (i.e.
using more line-sinks) in the near field and less detail in the far field. In our steady state model, we include
the major perennial streams (Wabash River, White River) that flow all year round.
Note that perennial streams show up as solid blue lines on USGS topographic maps, and dashed if intermittent or
ephemeral.
Try to use short length line-sinks in the near field and longer length line-sinks in the far field. Try to
keep the total amount of line-sinks below 200 for efficient model operation. You can check this by clicking
Model and Model size.
Your layout of elements should look something like Figure 5.5.
Step 5 Add areal recharge element.

Now let’s add areal recharge based on the USGS suggested value of 12.1 in/yr (0.00276 f t/d) [Shedlock,
1980]. To do this, from the Menu, select Element > New > Inhomogeneity. Type in the recharge rate of
0.00276 f t/d. See Figure 5.6. Click OK. Use the mouse and click to define the boundary vertices. Surround
the line-sink elements with an inhomogeneity polygon. We will further the discussion of the inhomogeneity
elements later in this chapter.
Your model space should look something like Figure 5.7.

Step 6 Solve and plot the 5-year capture zone for the well field.

In Model , Settings, select the Contouring tab and check the Show Contours box of detailed
resolution with minimum contour of 350 f t and maximum contour of 500 f t and interval 1 f t. Select the
Tracing tab and check the Compute Particle Paths. Set the maximum travel time to 1825 d and use
default step size. Click OK. Double click on the well and select the “other” tab. Check the box Trace
Particles From Well. Select Number of Particles 20 and starting elevation 335 f t and choose pathline
color red. Click OK and then Solve or Run. Zoom the window down to the capture zone and re-solve.
Draw the wellhead protection area. View > Pathlines > Time-of-Travel tics > every one year. The
result should look similar to Figure 5.8.
Note: if no pathlines appear, ensure they are turned on under the View menu and that your aquifer parameters
are correct; if heads are drawn down to the aquifer bottom no pathlines will appear!

Step 7 Check the solution for integrity.

From the Model menu, select View Results Table, click on View and on Line-sinks. You are looking
at the vertices of the line-sink strings. Look further to the right in the table to ensure that the “Specified
5.2. RIVERS AND LINE-SINKS 49

Figure 5.5: Hydrography labels and line-sinks representing the Wabash River and White River.
50 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.6: Inhomogeneity element properties box.

Head” and “Modeled Head” for all line-sinks are nearly the same (less than one percent error). This ensures
a solution that meets the boundary conditions, at least at the line-sink centers.
It is sometimes possible to compare the cumulative “discharge” of a line-sink string to the observed base
flow of a stream. That observation could be based on the difference between base flow separation estimates
at two USGS gauging stations, or a base flow separation at a USGS gauge on a head waters reach. This
procedure would provide an additional test of the reasonableness of the model water balance (which is insured
in the analytic element method) in comparison to the field, provided the catchment area between the two
gaging stations is part of the nearfield in the model domain.
It is also useful to compare the model-predicted heads to observed heads at monitoring wells. Fortunately,
the USGS has reported a series of measured water levels in the study area for the period January 23-25,
1978 [Shedlock, 1980]. These water levels have been stored in a “Test Point” file. To import this file, select
Tools and then Import the file vin-mw.tp. Solve to generate an updated result. The view will show the
test point locations as triangles, with the size and direction related to the difference between modeled head
and observed head (See Figure 5.9). Click on a triangle to report its error on the bottom status bar. Notice
we have some large errors in the model to the south. A summary of the test point statistics can be shown
by Model > View Calibration Statistics. The conceptual model will be refined in the next section.
It is good practice to keep separate test point files to distinguish between differences in water level
observations, such as differences in dates of observations, differences in well screen elevations, differences in
quality of the measurements, etc. There is an Edit > Delete All > Test Points sequence to assist bringing
different sets of test points to the view.
Calculated heads are expected to differ from observed heads for many reasons, most importantly because
the model aquifer is merely an abstraction of the real aquifer system. A good model will show test point
5.2. RIVERS AND LINE-SINKS 51

Figure 5.7: Hydrography labels, line-sinks, and inhomogeneity (recharge) element (blue rectangle).
52 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.8: 5 year capture zone including the influence of hydrologic boundaries.

Figure 5.9: Modeled head contours compared to field observations at monitoring wells.
5.3. GEOLOGIC CONTACTS AND NO-FLOW ELEMENTS 53

triangles that point both up and down, whereby the differences between model head and observed head are
not too large on average. A good model will also show a rather random distribution of up and down arrows,
instead of having all up arrows clustered in one area and all down arrows in another.
In an ideal, homogeneous, single layer aquifer, some general rules apply for heads that are too high or too
low. Heads that are consistently too low close to the well field suggest too low a transmissivity in the model.
Heads that are consistently too high in areas of ground- water mounding might also indicate too low a model
transmissivity or too high an areal recharge rate (or both). When the aquifer is significantly stratified or
consists of several vertically stacked interconnected aquifers, discrepancies between observed and calculated
heads can also be explained based on the relative depth of the piezometers (wells) in the field. In most
cases, shallow piezometers in a stratified aquifer tend to show higher heads than predicted in a homogeneous
aquifer, while deep piezometers tend to show lower heads. Only close to discharge areas, such as streams or
lakes, can this trend be reversed.
Modeled heads that are too high or too low can also be the result of improper representation of surface
waters. For instance, the headwaters of a creek can artificially raise the ground-water table in WhAEM by
infiltrating large amounts of water (that these headwaters don’t have to give). Such stream sections and
their line-sinks should be removed from the WhAEM model, particularly if they show up as ephemeral or
intermittent steams on USGS topo maps.

5.3 Geologic Contacts and No-Flow Elements

The geology of our study area deserves a closer look at this point as we refine the conceptual model. The
area is underlain by unconsolidated sediments as shown in Figure 5.10. These sediments were deposited by
glacial processes in the Wabash River valley during the Pleistocene Epoch and reworked by wind and water.
The outwash pinches-out at the margins of the river valley, and generally rests directly on the underlying
bedrock, as shown in Figure 5.11. A couple of rock “islands” outcrop just south of the wellfield.
Ideally, we should represent the low permeable rock and highly permeable outwash explicitly in our
model. As a first step, we will represent in WhAEM open strings of line elements to create a no-flow
boundary positioned along the boundary of the outwash or channel deposits to approximate the effect of an
absolute reduction in permeability. These no-flow boundaries should be extended into the far field.
WhAEM conceptual models with and without these no-flow boundaries can be viewed as bounding cases
of the real situation, at least for the property of heterogeneous hydraulic conductivity distributions. If
these two extreme models do not significantly affect the time of travel capture zones, no further (more
sophisticated) modeling seems necessary. On the other hand, if the two extreme models lead to drastically
different time of travel capture zones, then further refinement of the model is justified.

5.3.1 Exercise: No flow boundaries

Step 1 Create the no-flow boundary at the topographic contact between outwash and rock. Solve.

First, create a duplicate database as noflow1. You may find it helpful to add the USGS DLG hypsography
(topography) maps (VI1hp06.bbm, VI1hp07.bbm) to the basemap which already includes the USGS DRG
(vincennes 24k.tif). Locate the topographic feature Bunker Hill about 2km south of the wellfield. From the
menu, select Element > New > and Horizontal Barrier. Name the feature “Bunker Hill” and click OK.
54 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.10: Surficial geology of the study region. [Shedlock,1980][Gray et al.,1970].

5.3. GEOLOGIC CONTACTS AND NO-FLOW ELEMENTS 55

Figure 5.11: Generalized geologic cross section of the Wabash River Valley near the study region [Shed-
lock,1980][Gray et al.,1970].
56 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.12: Superposition of no-flow elements along the geologic contact inferred from the topographic
contours for Bunker Hill.

Use the mouse to create vertices along the northern boundary of the outcrop (left click); use the right click
to end.

Step 2 Solve and draw the 5 year capture zone.

Double click on the wellfield, and release around 20 particles. Check the Model > Settings > and
Solve. In View , set the Pathlines to Show time tics to every 1 year. Draw the wellhead protection
area (See Figure 5.12).
Apparently, the presence of a no-flow boundary has a significant effect on the flow patterns in the outwash,
particularly on the capture zone for our wellfield (compare Figure 5.8 and Figure 5.12). If in reality there
might be a reduction of water available from the till and rock formations to the east. It makes sense to
explore solutions for the problem where the till and rock formations are given a finite hydraulic conductivity,
as we shall see in the next section.

5.4 Aquifer Inhomogeneities

Domains with differing properties are enclosed within polygons defining the inhomogeneity using the WhAEM in-
terface. These domains can have different hydraulic conductivities, aquifer bottom elevations, effective poro-
sities, and/or areal recharge due to precipitation. WhAEM uses mathematical features called line-doublets
to simulate the polygon segments. Domains can be nested, but should not overlap. There are special consi-
derations when domains share segments or are abutted. Also, you should avoid sharp corners when designing
the domains. A line-sink should be either inside or outside the inhomogeneity domain. Extra vertices are
needed near high capacity wells. More details are available with the online help (Element > New > Inho-
mogeneity, then F1). See Figure 5.13. It may be useful to edit the vertices that make up the inhomogeneity
5.4. AQUIFER INHOMOGENEITIES 57

Figure 5.13: Drawing inhomogeneities using the mouse.

envelope. Select a vertex by right-clicking on a node; it will turn red. The options to move, delete, or refine
can be accessed through the Edit menu, or the smart icons.

5.4.1 Exercise: Inhomogeneities

Step 1 Create inhomogeneities using WhAEM to represent the outwash.

We will build an inhomogeneity domain to represent the outwash of the Wabash River with a hydraulic
conductivity of 350 f t/d, in contrast to the rock uplands area with a hydraulic conductivity of 10 f t/d.
The first thing to do is to make a duplicate database called inhomo1. Now go to Model , then Settings,
and enter a new hydraulic conductivity as the background value of 10 f t/d.
We no longer want the simple recharge distribution, so select the rectangular inhomogeneity polygon
with clicking the left mouse button close to a vertex/corner, and then hit the delete key to erase the entire
string.
Also delete the no flow strings created in the previous exercise.
If not already included, it will be helpful to add the USGS DLG of hypsography for the study area to the
set of vector maps in the basemap. Go to Project , then Project Settings, and add maps VI1hp06.bbm,
VI1hp07.bbm.
To build the outwash inhomogeneity, select Element > New > Inhomogeneity. You will be presen-
ted with the opportunity to parameterize the inhomogeneity (See Figure 5.14). Label the inhomogeneity
“outwash”. Enter the new hydraulic conductivity for the outwash domain k = 350 f t/d. Add the recharge
N = 0.00276 f t/d. Click OK. Use the hypsography basemaps (vi1hpf06, vi1hpf07) to guide your drawing of
the polygon, following the sharp transition from outwash to rock (See Figure 5.15).

Step 2 Refine the linesinks to overlap the inhomogeneities.

58 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.14: Dialog window for entering inhomogeneity properties.

The vertices of linesinks and inhomogeneities should overlap where these elements cross (See online
help). Therefore, refine the linesink vertices such that the linesinks are completely inside and outside the
inhomogeneities. Also note that inhomogeneities should not have sharp corners and the refine and move
tools can be used to smooth the corners. See Figure 5.15.

Step 3 Solve the model and plot the 5 year capture zones.

You are encouraged to check for improvement of the model solution by comparing the model predicted
heads to the test points ( Tools > Import > Test Point File > vin mw.tp). The capture zone solution is
shown in Figure 5.16. It is an additional exercise to explore the influence of the background rock hydraulic
conductivity on the local 5-yr capture zone. Once the background rock hydraulic conductivity is less than
two orders of magnitude less than outwash hydraulic conductivity it performs similar to a no-flow boundary.
The capture zone is relatively insensitive to the exact value of the rock hydraulic conductivity.

Advanced Use: When only the recharge differs, inhomogeneity domains can have arbitrarily large line-doublet sides, a
domain can overlap any other feature, including other inhomogeneity domains and line-sinks, and can have any shape,
including sharp corners. When only porosity changes/jumps, but not hydraulic conductivity or base elevation, the
inhomogeneity domains cannot overlap other inhomogeneity domains, but can take any shape, including sharp corners
and large sides, and can intersect line-sinks, and the domains may share sides (See online help description).
5.4. AQUIFER INHOMOGENEITIES 59

Figure 5.15: Layout of linesinks, wellfield, and aquifer inhomogeneity domain. Red circles show intersections
of linesinks with inhomogeneity.

Figure 5.16: 5 year capture zone solution including the influence of aquifer inhomogeneities.
60 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.17: Conceptual model of stream/lake resistance.

5.5 Riverbed Resistance

How do you know if riverbed resistance is important in your field study? You have resistance if there is a
significant difference between water levels in the stream and heads below the stream. What causes resistance?
Resistance is caused by vertical flow near a stream/lake, silt and muck on the stream/lake bottom, a clay or
till layer underneath the stream/lake, or an aquitard between aquifers (multi-aquifers).
A conceptual model of stream/lake resistance is shown in Figure 5.17. Resistance c [T ] is equal to the
thickness of the layer (δ) divided by the hydraulic conductivity (kc ) of the layer.
δ
c= (5.2)
kc
The vertical specific discharge qz [L/T ] is equal to the head in the aquifer (φa ) minus the head in the stream
(φs ) divided by the resistance (c).
φa − φs
qz = (5.3)
c
The discharge per length of the stream segment σ[L3 /T /L] is equal to the specific discharge times the width.

σ = wqz (5.4)

What values might one expect of the resistance parameter? Major rivers with sandy bottoms exhibit
resistance between 0 and 1 d. Small streams with silty bottoms have resistance in the range between 1 and
10 d. Stream or lakes in till formations (overlying a deeper aquifer) exhibit resistance up to 100 d. Resistance
greater than 100 d means almost no interaction with the aquifer. Resistance at a specific site is dependent
on both the hydraulic conductivity of the riverbed and the thickness of the bed layer.
The theoretical basis of the implementation of resistance, depth, and “effective width” in the ground
water solver GFLOW1 of the WhAEM framework is described in the document Dealing with resistance
to flow into surface waters.pdf Haitjema [2005]. The user will appreciate that the methodology has been
automated to calculate the effective width parameter associated with an effective leakage zone underneath
the surface water body, given a resistance, layer depth, and surface water width. This effective leakage zone
depends on the bottom resistance and the aquifer transmissivity, and may be calculated. You must do this
calculation by hand if you check the radio button “Unknown” for the line-sink location. If you check the
radio button “Along surface water boundary” or “Along stream center line”, the actual surface water width
5.5. RIVERBED RESISTANCE 61

Figure 5.18: Dialog box for entry of resistance linesinks.

must be entered as the Width parameter (the parameter B in the referenced document). The Solver will
automatically calculate the correct effective Width parameter, which will be reported back to the GUI as
the width of the line-sink in the Model Results table (on the Model menu).
Note: There is another situation where the width parameter on the Linesink String Properties dialog may be set
to the actual stream or lake width, and a resistance entered to represent resistance to three-dimensional flow.

5.5.1 Exercise: Resistance Linesinks

Step 1 Calculate the resistance of the bottom materials in the Wabash River.

The USGS report [Shedlock, 1980] includes a discussion of the leakage from the Wabash River to the
outwash aquifer in the vicinity of the city pumping center. They report a leakage layer of average hydraulic
conductivity 0.03 ft/d, and of local values of 0.07 and 0.7 f t/d, and of average thickness of 1 f t. Let’s work
with the value of kc = 0.07 f t/d and δ = 1 f t, and calculate the resistance using equation 5.2.

Step 2 Estimate the width of the River near the pumping center.

Estimate the width of the Wabash River near the pumping center location using the raster base map
vincennes 24k.tif. You may elect to use a small footprint geoviewer software package, such as the USGS
dlgv32pro, now known as Global Mapper evaluation version (www.globalmapper.com).

Step 3 Enter resistance line-sinks to represent the River near-field.

Let’s build off of the previous database. Open the project database inhomo1.whm. Create a duplicate
database resist1. Click on the linesink string representing the Wabash River south. Unactivate the “Treat
as farfield” box. Enter the resistance, depth, and width parameters as calculated above (See Figure 5.18).
Do the same for the Wabash River north.
62 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.19: 5-year capture zones including resistance of Wabash River (c = 14.3 d, δ = 1 f t, width B
=615 f t, 36 particles released from the well
.

Step 4 Solve the model and plot the 5 year capture zone.

(See Figure 5.19).

5.6 Discussion
In the face of uncertainty, it is prudent to present more than one piezometric contour map and more than
one time-of-travel capture zone. The impact of the most sensitive parameters on the capture zone should be
illustrated by presenting the different capture zones and clearly explaining the different parameter choices
that lead to them. It is much more convincing to show that, for instance, a 50 f t/d difference in hydraulic
conductivity does not noticeably change the capture zone than argue at length why the particular value
chosen for the model is the right one. You can annotate the graphics (piezometric contour maps or capture
zones) directly on the screen by use of the Add Text option on the Edit menu. The maps are printed
using the Printer setup and Print options on the Project menu. You can also export the graphics using
Tools to DXF or Shapefile for further post processing.
Haitjema investigated the influence of complexity of the conceptual model on the capture zones for the
study area [Haitjema,1995]. He used the analytic element code GFLOW [Haitjema,1995] to include the
influence of the low hydraulic conductivity of the rock outcrops. The finite difference code MODFLOW
[McDonald and Harbaugh, 1988] was used to represent localized three-dimensional flow, anisotropy, and
resistance between the river and the aquifer. In carrying-out this modeling, progressively more complex
conceptual models were tested. Results indicated that three-dimensional effects did not impact the capture
zones, but different values for the resistance of the Wabash River bottom did. The reader is directed to
Chapter 6.1 in [Haitjema, 1995] for a full discussion. In Figure 5.20, “composite capture zones” are depicted,
based on the many capture zones generated. The capture zones were produced in GFLOW using only 2D
5.6. DISCUSSION 63

flow, but included resistance to flow from or into the Wabash River.
Based on the exercises in this guide for the city wellfield, it may be shown that the 5 year capture zone
solution is sensitive to the existence of the rock outcrops but insensitive to the exact value of the hydraulic
conductivity of the rock outcrops — a no-flow boundary (assuming the rock is impermeable) may capture the
essence. However, we did learn that the 5 year capture zone is sensitive to the resistance of the Wabash River
channel — the capture zone extended all the way to the north side of the river when resistance was given
a conservative (high end) value. Thus, more investment in characterizing the riverbed resistance through
field observations and modeling simulations might be justified. Alternatively, one can assume the worst and
accept the capture zone to the north side of the river as realistic.
Remember, the case study for the city wellfield was for illustration purposes, and is not to be misunder-
stood as the approved State source water assessment.
Now some other complexities. We have not yet introduced the concept of horizontal anisotropy in the
hydraulic conductivity distribution. Horizontal anisotropy is known to influence the capture zone [Cole and
Silliman, 1997]. There is the temptation to represent a highly fractured rock aquifer as an equivalent porous
medium with anisotropic hydraulic conductivity, at least at the regional scale; this should be done with
caution in fractured carbonate rock aquifers, [Kraemer, 1990],[Podgorney and Ritzi,1997]. WhAEM should
not be used in aquifers with a highly anisotropic hydraulic conductivity distribution.
We also have not introduced the concept of dispersion to this point. The capture zones presented in this
document represent average water residence times in the subsurface. These residence times more closely repre-
sent the breakthrough time of the mean concentration of a conservative tracer at the well. If you wanted to re-
present the first arrival of the tracer in your delineations, and thus capture the initial breakthrough of a poten-
tial contaminant, you would need to incorporate macroscale dispersion (www3.epa.gov/ceampubl/learn2model/part-
two/onsite/conc.html). Characterizing dispersion in capture zone models is a topic of active research [Maas,
1999].
We have limited this guidance to manual methods of model calibration. The reader is referred to Hill [
1998] for a discussion of automated methods for parameter estimation and model calibration. Automated
non-linear parameter estimation is implemented in UCODE (igwmc.mines.edu/freeware/ucode/) and PEST
(www.pesthomepage.org/).
There are other professional analytic element modeling packages, e.g., AnAqSim, GFLOW, and MLAEM.
The TTim model is an open source project for multi-aquifer transient flow systems using analytic elements
((github.com/mbakker7/ttim). The analytic element modeling community maintains a web page with links to
the above analytic element ground water models (www.analyticelements.org).
64 CHAPTER 5. PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA

Figure 5.20: Composite time of travel capture zones with 2-,5-, and 10-year isochrones for the city well field
using the full-feature analytic element model GFLOW (after [Haitjema,1995]).
Chapter 6

SUMMARY

6.1 Good Modeling Practice

Fortunately, we do not always have to repeat the procedure of hypothesis testing with increasingly complex
conceptual models for each individual modeling project. After a while some general conclusions can be drawn
by comparing simple and complex conceptual models. In many cases the validity of certain assumptions has
already been investigated by others and published in the literature. It is not practical to present here a
complete list of valid and invalid conceptual models, covering all possible geohydrological circumstances. It
is the responsibility of the modeler to refer to the ground-water literature and make these choices on a case
by case basis (See Haitjema [2006]). We do, however, provide a few general rules-of-thumb that may be of
help in making some initial decisions regarding conceptual models.

6.1.1 Rules-of-Thumb
1. River Proximity. The residence-time-based protection zone, or capture zone, is not likely influenced
p
by surface water when the distance of the well to the surface water feature is greater than 4 QT /(πHn),
where Q is the pumping rate of the well, H is the aquifer saturated thickness, n is the aquifer porosity,
and T is the residence time associated with the well capture.

2. Three-Dimensional Flow. Three-dimensional flow modeling may be necessary when one of the
following circumstances apply:
p
(a) The well or well field is within a distance of 2 kh /kv H from a hydrological boundary, where kh
and kv are the horizontal and vertical hydraulic conductivity, respectively, and where H is the
average saturated aquifer thickness.
p
(b) The width of the capture zone in a two-dimensional model is less than 2 kh /kv H, where kh and
kv are the horizontal and vertical hydraulic conductivity, respectively, and where H is the average
saturated aquifer thickness.

Note that in the event that the (2D) capture zone width is about equal to the aquifer thickness, the
three-dimensional capture zone is likely to be a slender, three-dimensional stream tube occuring at
some depth below the water table. A conservative (protective) time-of-travel capture zone is obtained
for this case by using a wellhead protection zone obtained from a horizontal flow model assuming a
confined aquifer of the same thickness as the length of the well screen.

65
66 CHAPTER 6. SUMMARY

3. Multi-Aquifer Flow. The subsurface is usually stratified (layered), with each layer having its own
hydraulic conductivity. If the hydraulic conductivities of the various formations are all within the same
order of magnitude, the layers form one stratified aquifer, see item 4 in this list. If, on the other hand,
one or more layers have a hydraulic conductivity that is more than an order of magnitude lower than
the more permeable layers such layers are referred to as aquitards. Aquitards divide the subsurface into
more than one aquifer, forming multi-aquifer systems. In case of a multi-aquifer system, a multi-aquifer
model (quasi three-dimensional model) may be necessary if all of the following circumstances apply:

(a) The well occurs in a multi-aquifer zone and is not screened in all aquifers.
(b) The length of the capture zone, based on flow in the aquifer(s) accessed by the well, is in the
√
order of (or smaller) than 4 T c, where T (f t2 /d) is the transmissivity of the aquifer and c (d)
is the largest resistance of one of the adjacent aquitards. Resistance (c) is equal to the aquitard
thickness (b) divided by its hydraulic conductivity (k) or c = b/k.
√
(c) A hydrological boundary (stream or lake) is present within a distance of 4 T c from the well,
where T (f t2 /d) is the transmissivity of the aquifer and c (d) is the largest resistance of one of
the adjacent aquitards.

Multi-aquifer flow often greatly complicates the delineation process due to uncertainties in aquifer
interaction. In many cases, aquifers do not interact through a homogeneous leaky aquitard, but through
discrete openings in clay layers, whose presence and location may or may not be known. If a well is
screened in only one aquifer, the largest time-of-travel capture zone will occur if interaction with
adjacent aquifers is ignored. If the directions of flow in the various aquifers are expected to be not too
different, the single aquifer capture zone may offer a conservative time-of-travel capture zone.

4. Aquifer Stratification. An aquifer is considered stratified if the hydraulic conductivities of the

various geological formations (layers) are different but within the same order of magnitude. Aqui-
fer stratification will not affect the capture zone width, but may significantly affect the capture zone
length. The groundwater flow velocities in the various aquifer layers are proportional to their hydraulic
conductivities and inversely proportional to their porosities. This can lead to substantially different
time-of-travel capture zone lengths in the various layers. If it is considered important to include strati-
fication in the definition of a wellhead protection zone, a time-of-travel capture zone can be calculated
using a porosity value for the aquifer that incorporates both the effect of the highest conductivity
and the associated porosity. For instance, for an aquifer layer with five times the average hydraulic
conductivity the model should be given a porosity that is five times smaller than the average porosity
in order to get the correct time-of-travel capture zone length in that layer. This is assuming that the
porosity is the same in all layers.

5. Aquifer Heterogeneities. Small inclusions of clay or gravel with different hydraulic conductivities
than the average regional value are inconsequential, except if their size is on the order of the capture
zone width or larger.

6. Transient Flow. Highly permeable unconfined aquifers, and most confined aquifers, will exhibit
“summer conditions” and “winter conditions”. These can be approximated by steady-state solutions
using recharge rates and surface water levels observed in the summer and winter, respectively. When
delineating capture zones for wells based on residence times greater than say 5 years, the detailed
day-to-day pumping records can safely be averaged over the observation period, provided there is no
6.2. CONCLUSION 67

Figure 6.1: Decision tree for complexity of capture zone delineation method, after [Haitjema,2002].

long term trend. Accurate delineation for transient non-community wells may require special attention;
projections into the future should include additional safety factors.

6.1.2 Choosing the Right Tool

The technique used to delineate time-of-travel capture zones for wellhead protection areas is a balance
between geohydrological complexity, data availability, time constraints, budget, and other factors. A decision
tree can assist in the justification of approach based on geohydrological complexity, as shown in Figure
6.1. Back-of-the-Envelope calculations are a good starting point for most problems [Haitjema, 2006]. If
Back-of-the-Envelope capture zones are not justified, then Simple WHPAs and geohydrological modeling
should be considered. The analytic solutions in the Simple WHPAs tool include well(s)-in-uniform-flow
solution with simple, linear-shape boundary conditions. Geohydrological modeling using WhAEM involves
solving the regional ground water flow differential equations using the analytic element method. Here the
conceptual possibilities include more complex river boundary conditions and piece-wise constant area recharge
distributions and hydraulic conductivity zones. Where can one get other analytic element models? Check-
out the community page (www.analyticelements.org/software.html). If the conceptual model requires multiple
aquifers, then the USGS finite difference model MODFLOW [USGS, 2017], [Aquaveo, 2017], or the analytic
solvers TimM L [Bakker, 2007] or MLAEM [Strack, 2006] are available. If a complex hydraulic conductivity
distribution is known, MODFLOW or a finite element solver may be the better choice.

6.2 Conclusion
WhAEM is a computer tool to support step-wise, progressive modeling and delineation of source water
areas for pumping wells. Our case study moved from arbitrary fixed radius, to calculated fixed radius, to
well-in-uniform-flow, to modeling with hydrologic boundary conditions only, to modeling with hydrologic
and geologic boundary conditions. Each solution was conceptually more sophisticated, and we assume
68 CHAPTER 6. SUMMARY

that the corresponding calculated capture zones were progressively more realistic. The emphasis of the
WhAEM project on “ease-of-use” and computational efficiency does not release you, the modeler, from
responsibilities in justifying the conceptual models, and defending the reasonableness of the solutions. We
have emphasized the uncertainties in conceptualization of the boundary conditions in our case study, but
uncertainties in parameterization are also important. WhAEM is public domain and no-cost software, and
while it is intended to capture the basics of ground-water modeling in support of State efforts in source
water assessments and wellhead protection, it is also intended to stimulate innovation in software design and
modeling practice in the private sector.
One more time, “Ground-water models should be as simple as possible, but not simpler.” Good modeling!
Bibliography

[Aquaveo, 2017] Aquaveo (2017). GMS groundwater modeling system v 10.3.5.

https://www.aquaveo.com/software/gms-groundwater-modeling-system-introduction.

[Bakker, 2007] Bakker, M. (2007). TimM L : a multiaquifer analytic element model.

https://github.com/mbakker7/timml.

[Bakker and Strack, 1996] Bakker, M. and Strack, O. (1996). Capture zone delineation in two-dimensional
groundwater. Water Resources Research, 32(5):1309–1315.

[Blandford and Huyakorn, 1993] Blandford, T. and Huyakorn, P. (1993). WHPA: a semi-analytical model
for the delineation of wellhead protection areas, version 2.2.

[Ceric, 2000] Ceric, A. (2000). Assessment of the applicability of simplified capture zone delineation techni-
ques for groundwater public water supply systems. Master’s thesis, School of Public and Environmental
Affairs, Indiana University-Bloomington.

[Ceric and Haitjema, 2005] Ceric, A. and Haitjema, H. (2005). On using simple time-of-travel capture zone
delineation methods. Ground Water, 43(3):408–412.

[Cole and Silliman, 1996] Cole, B. E. and Silliman, S. E. (1996). Estimating the horizontal gradient in
heterogeneous, unconfined aquifers: comparison of three-point schemes. Ground Water Monitoring and
Remediation, Spring:85–90.

[Dupuit, 1863] Dupuit, J. (1863). Etudes Theoriques et Pratiques sur le Mouvement des Eaux dans les
Canaux Decouverts et a Travers les Terrains Permeables, 2nd ed. unknown, Dunod, Paris.

[Feyen et al., 2001] Feyen, L., Beven, K., Smedt, F. D., and Freer, J. (2001). Stochastic capture zone
delineation within the generalized likelihood uncertainty estimation methodology: conditioning on head
observations. Water Resources Research, 37(3):625–638.

[Gray et al., 1970] Gray, H., Wayne, W., and Wier, C. (1970). Geologic map of the 1◦ x 2◦ vincennes
quadrangle and parts of adjoining quadrangles, Indiana and Illinois, showing bedrock and unconsolidated
deposits. Map, Indiana Department of Natural Resources, Geological Survey Division. scale 1:250,000.

[Haitjema, 1995] Haitjema, H. (1995). Analytic Element Modeling of Groundwater Flow. Academic Press,
San Diego. http://www.haitjema.com/howtoorder.html.

[Haitjema, 2002] Haitjema, H. (2002). Time of travel capture zone delineation for wellhead protection.
Environmental Science Research Center, Indiana University, Bloomington. Prepared for Drinking Water
Branch, Indiana Department of Environmental Management, Indianapolis, Indiana.

69
70 BIBLIOGRAPHY

[Haitjema, 2006] Haitjema, H. (2006). The role of hand calculations in ground water flow modeling. Ground
Water, 44(6):786–791.

[Haitjema, 2011] Haitjema, H. (2011). Brief analysis for including recharge in calculated radius time-of-travel
(TOT) capture zones. personal communication.

[Haitjema et al., 1995] Haitjema, H., Strack, O., and Kraemer, S. (1995). Demonstration of the analytic
element method for wellhead protection. EPA Project Summary EPA/600/SR-94/210, U.S. Environmental
Protection Agency, Office of Research and Development.

[Haitjema et al., 1994] Haitjema, H., Wittman, J., Kelson, V., and Bauch, N. (1994). WhAEM: Program
documentation for the wellhead analytic element model. EPA Project Report EPA/600/R-94/210, U.S.
Environmental Protection Agency, Office of Research and Development.

[Haitjema, 2005] Haitjema, H. M. (2005). Dealing with resistance to flow into surface waters.
http://haitjema.com/surface.html.

[Heath, 1983] Heath, R. (1983). Basic ground water hydrology. Water-Supply Paper 2220, U.S. Geological
Survey. republished in a 1984 edition by the National Water Well Association, Dublin, OH.

[Hill, 1998] Hill, M. C. (1998). Methods and guidelines for effective model cali-
bration. Water Resources Investigation Report 98-4005, U.S. Geological Survey.
https://water.usgs.gov/nrp/gwsoftware/modflow2000/WRIR98-4005.pdf.

[Kelson, 1998] Kelson, V. (1998). Practical advances in ground-water modeling using analytic elements. PhD
thesis, Indiana University-Bloomington.

[Kelson et al., 1993] Kelson, V., Haitjema, H., and Kraemer, S. (1993). GAEP: a geographical preprocessor
for ground water modeling. Hydrological Science and Technology, 8(1-4):74–83.

[Konikow and Bredehoeft, 1992] Konikow, L. F. and Bredehoeft, J. D. (1992). Ground-water models cannot
be validated. Advances in Water Resources, 15:75–83.

[Kraemer, 1990] Kraemer, S. R. (1990). Modeling of Regional Groundwater Flow in Fractured Rock Aquifers.
PhD thesis, Indiana University-Bloomington.

[Kraemer, 2007] Kraemer, S. R. (2007). Analytic element ground water modeling as a research program
(1980 to 2006). Ground Water, 45(4):402–408.

[Maas, 1999] Maas, C. (1999). A hyperbolic dispersion equation to model the bounds of a contaminated
groundwater body. Journal of Hydrology, 226:234–241.

[McDonald and Harbaugh, 1988] McDonald, M. and Harbaugh, A. (1988). A modular three-dimensional
finite-difference ground-water model. Techniques of Water-Resources Investigations ,Book 6, Chapter A1,
U.S. Geological Survey. 586p.

[Merkle et al., 1996] Merkle, J., Macler, B., Kurth, L., Bekdash, F., Solomon, L., and Wooten, H. (1996).
Ground water disinfection and protective practices in the United States. U.S. Environmental Protection
Agency Office of Ground Water and Drinking Water and Science Applications International Corporation
Report.
BIBLIOGRAPHY 71

[Molz, 2015] Molz, F. (2015). Advection, dispersion, and confusion. Groundwater, 53(3):348–353.
doi:10.1111/gwat.12338.

[Podgorney and Ritzi Jr., 1997] Podgorney, R. K. and Ritzi Jr., R. W. (1997). Capture zone geometry in a
fractured carbonate aquifer. Ground Water, 35(6):1040–1049.

[Reilly, 2001] Reilly, T. E. (2001). System and boundary conceptualization in ground-water flow simula-
tion. Techniques of Water-Resources Investigations Book 3, Applications of Hydraulics, Chapter B8, U.S.
Geological Survey. https://pubs.usgs.gov/twri/twri-3 B8/.

[Shedlock, 1980] Shedlock, R. (1980). Saline water at the base of the glacial-outwash aquifer near Vincennes,
Knox County, indiana. Technical Report 80-65, U.S. Geological Survey.

[Strack, 1989] Strack, O. (1989). Groundwater Mechanics. Prentice Hall, Englewood Cliffs, NJ. out of press.

[Strack, 2017] Strack, O. (2017). Analytical Groundwater Mechanics. Cambridge University Press, Cam-
bridge, UK. ISBN:9781107148833.

[Strack et al., 1994] Strack, O., Anderson, E., Bakker, M., Panda, W., Pennings, R., and Steward, D. (1994).
CZAEM user’s guide: Modeling capture zones of ground-water wells using analytic elements. EPA Project
Report EPA/600/R-94/174, U.S. Environmental Protection Agency, Office of Research and Development.
https://www.epa.gov/water-research/capture-zone-analytic-element-model-czaem.

[Strack, 2006] Strack, O. D. (2006). MLAEM - multilayer analytic element model, version 5.02. Strack
Consulting, Inc. http://www.strackconsultingLLC.com.

[USEPA, 1993] USEPA (1993). Guidelines for delineation of wellhead protection areas. Technical Report
EPA/440/5-93-001, U.S. Environmental Protection Agency, Office of Water Office of Ground Water Pro-
tection, Washington, DC.

[USEPA, 1994] USEPA (1994). Handbook: ground water and wellhead protection. Technical Report
EPA/625/R-94-001, U.S. Environmental Protection Agency, Office of Research and Development, Cin-
cinnati, OH.

[USEPA, 2017] USEPA (2017). Protect sources of drinking water. Office of Ground Water and Drinking
Water, https://www.epa.gov/sourcewaterprotection. updated.

[USGS, 1997] USGS (1997). The universal transverse mercator (utm) grid. Fact Sheet 142-97, U.S. Geolo-
gical Survey.

[USGS, 1999] USGS (1999). Map accuracy standards. Fact Sheet 171-99, U.S. Geological Survey.

[USGS, 2017] USGS (2017). MODFLOW and Related Programs. https://water.usgs.gov/ogw/modflow/.

[Vassolo et al., 1998] Vassolo, S., Kinzelbach, W., and Schafer, W. (1998). Determination of a well head
protection zone by stochastic inverse modeling. Journal of Hydrology, 206:268–280.
72 BIBLIOGRAPHY

This page is blank.

Appendix A

Conversion Factors

Multiply By To obtain
Length f eet(f t) 0.3048 meter(m)
mile(mi) 5280 ft
mi 1609.344 m
yard(yd) 0.9144 m
Area f t2 0.092903 m2
acre 4046.8564 m2
mi2 2.589988 x 106 m2
hectare(ha) 10, 000 m2
Rate inches/year(in/yr) 2.2815 x 10−4 f t/day(d)
in/yr 6.9541 x 10−5 m/d
f t/d 0.3048 m/d
gallons per day(gpd)/f t2 0.0407 m/d
Volume rate f t3 /second(s) 2, 446.5754 m3 /d
f t3 /s 86400 f t3 /d
gal/min 5.4510 m3 /d
gal/min 192.5 f t3 /d
gal/d 0.0038 m3 /d
gal/d 0.13368 f t3 /d

73
74 APPENDIX A. CONVERSION FACTORS

This page is blank.

Appendix B

WhAEM Base Maps and the Internet

B.1 Vector Basemaps

Geographically registered electronic maps that include points, lines, and polygons are very useful basemaps
for representation of surface water (hydrography), transportation, topography contours (hypsography), soils
and geology. WhAEM can import and display vector maps created from binary base maps (*.bbm) or
shapefiles (*.shp).

B.1.1 Getting Binary Base Maps (BBM) from the EPA Web Server
A Binary Base Map server is supported by the EPA Center for Exposure Assessment Modeling (CEAM). To
access the server point your web browser to ??www.epa.gov/exposure-assessment-models/whaem2000-bbm-files-
us. The web page is shown in Figure B.1. Left click on the State, then the 30x60 minute map, then the 15
minute series. Save the self extracting file to your local drive and data folder, and uncompress (by right-click
and unzipping the file). The binary base maps (bbm) can be loaded into your project database by using the
Project and Project Settings, and Add Map commands.
The EPA-supplied BBMs are a binary form of the USGS Digital Line Graph (DLG) maps (optional
format, 1:100,000 scale, UTM projection, most often NAD27 but sometimes NAD83 datum, check with
USGS). The BBMs are saved in files associated with their location and category, and the name reflects
this information. The 1:100,000 scale map is divided into eight 15-minute quads (See Table B.1). Each of
these quads contains the DLG for the associated category: hypsography, hydrography, roads, railroads, and
miscellaneous(See Table B.2). The EPA server has only limited availability of hypsography coverages. As
an example of the file naming convention, the hydrography quad for Vincennes, IN, is VI1HYF06, which
refers to the sixth quadrant (F06) of the Vincennes 1:100,000 map (VI1), containing hydrography (HY).

Table B.1: BBM index by Quads

F01 F02 F03 F04
F05 F06 F07 F08

B.1.2 USGS Digitial Line Graphs (DLG) on EarthExplorer

The USGS DLGs are available for download at their EarthExplorer webpage. Point your web browser to
earthexplorer.usgs.gov(See Figure B.2).

75
76 APPENDIX B. WHAEM BASE MAPS AND THE INTERNET

Figure B.1: The EPA binary base map server.

Figure B.2: USGS EarthExplorer web page.

B.2. RASTER BASEMAPS 77

Table B.2: BBM by Type

Category Name Abbreviation
hypsography hp
hydrography hy
boundaries bd
roads and trails rd
railroads rr
miscellaneous transportation mt

The DLG maps are available at two scales: large scale 1:24K or intermediate scale 1:100K. The maps
are stored in the original optional format or the newer SDTS format. Complete coverage is available of the
USA at the 1:100K scale, while coverage at the 1:24K scale is not complete. You will need to use mapping
software to convert to a shapefile format for import into WhAEM . An example of a large scale 7.5 minute
quad map showing the elevation contours (hypsography) for Vincennes IN is shown in Figure B.3.

B.1.3 USEPA BASINS maps (Shapefile)

The USEPA watershed modeling system BASINS has data downloads organized by USGS hydrological unit
code or HUCs www.epa.gov/exposure-assessment-models/basins. The vector maps are available in shapefile
format. The BASINS user has the opportunity to customize the map projection (you can save as UTM
to be consistent with your WhAEM project). Maps relevant for source water protection may include land
use/land cover and soils maps.
Shapefiles can be added to the WhAEM basemap through the sequence Project > Project Settings
> Add Map.

B.2 Raster Basemaps

A Raster basemap is a georeferenced electronic map that is composed of a matrix of cells (or pixels). Common
examples are aerial photographs, imagery from satellites, and scanned maps. WhAEM is able to import and
view georeferenced raster basemaps in geotiff format (*.tif).

B.2.1 USGS digital raster graphic (DRG) maps on National Map Viewer
The USGS has scanned their historical paper topographic maps for use on a computer. WhAEM is capable
of using the USGS digital raster graphics (DRG) in geotif format (extension tif) as a base map.
The USGS digital raster graphics (DRG) can be selected and downloaded from viewer.nationalmap.gov.
See Figure B.4. Also note that NHD hydrography vector maps are also available.
When dealing with raster maps, it is useful to have access to a software package that allows clipping and
projecting, such as Global Mapper (www.globalmapper.com/).
Projected GeoTif files can be added to the WhAEM basemap through the sequence Project > Project
Settings > Add Map. Be sure to include the complete set of files in the folder storing the map files (e.g,
projection file *.tfw).
78 APPENDIX B. WHAEM BASE MAPS AND THE INTERNET

Figure B.3: USGS 7.5 minute quad map showing topographic elevation contours for Vincennes, IN.
B.3. USGS NATIONAL WATER INFORMATION SYSTEM 79

Figure B.4: USGS National Map Viewer, a source for DRG maps with historical topography.

B.2.2 GoogleEarth maps (tif )

It is possible to save a GoogleEarth map as a GeoTiff file and add to the WhAEM basemap. You will probably
need to register a GoogleEarth screen capture using tic marks of known location (e.g., latitude-longitude,
WGS84 datum), and process in a mapping package like Global Mapper, and export as a georeferenced tif
file in the projection of your basemap (e.g., UTM Zone 16, NAD27 datum, meters). See Figure B.5

B.3 USGS National Water Information System

The USGS online database has information of stage and discharge at thousands of stations on rivers throug-
hout the USA. The web link is waterdata.usgs.gov/nwis. The USGS stage information at stream gages can be
used as an alternative and more accurate method than using topo maps for creating a location of a known
water elevation (e.g., a hydro text label on the basemap).
80 APPENDIX B. WHAEM BASE MAPS AND THE INTERNET

Figure B.5: GoogleEarth map of Vincennes, IN, showing registration tic marks.