June 2018
Working with W/iAEM
Demonstration of Capture Zone Delineation for a
City Wellfield in a Valley Fill Glacial Outwash
Aquifer for Wellhead Protection
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

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 W/iAEM 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
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 (68W0f032) to Computer Sciences Corporation.
This document was typeset using WinEdt of Aleksander Simonic (www.winedt.com) and MikT^(including
LaT|5X2e) software of Christian Shenk (www.miktex.org). is a trademark of the American Mathematical

The purpose of this document is to demonstrate through an example study the use of the ground water
geohydrology computer program W/iAEM™ for Windows, which was developed by the U.S. Environmental
Protection Agency (EPA). The Wellhead Analytic Element Model (W/iAEM) 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. W/iAEM 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. W/iAEM has import and export utilities for DXF files and Shapefiles. W/iAEM 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).

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 (W/iAEM) 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. W/iAEM is a general purpose ground water hydrology desktop modeling tool based
on the innovative computational solution technique known as the analytic element method. W/iAEM assists
in the delineation and mapping of the area contributing source water to pumping wells, such as public
water supply wells. The development of W/iAEM 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.
Systems Exposure Division
Cincinnati, Ohio

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 W/iAEM to the Present	 5
1.5	The Base Map 	 6
1.6	Exercise: Exploring a USGS quad map	 7
2.1 Setbacks	 15
2.1.1 Exercise: Glacial outwash wellfield	 15
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.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.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

5.5.1 Exercise: Resistance Linesinks 		61
5.6 Discussion		62
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 W/iAEM Base Maps and the Internet	75
B.l Vector Basemaps		75
B.l.l Getting Binary Base Maps (BBM) from the EPA Web Server		75
B.l.2 USGS Digitial Line Graphs (DLG) on EarthExplorer		75
B.l.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

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.

This page is blank

Chapter 1
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 (W/iAEM) 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,

Table 1.1: Wellhead protection area delineation criteria and met
well in







draw down


residence time




assimilative capability

and geohydrologic modeling (see Table 1.1). This case study document reviews all of the criteria and their
methods using W/iAEM 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,

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

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

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
1.4 The Evolution of W/iAEM to the Present
The W/iAEM (pronounced warn) 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, W/iAEM 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
W/iAEM they also make the code easier to operate and thus easier to learn. Furthermore, W/iAEM requires
less input data. Time-of-travel capture zone delineation with W/iAEM 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 W/iAEM 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 W/iAEM must be employed. It is the responsibility of the
model user to determine whether W/iAEM 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 W/iAEM , a graphical user interface (GUI) and geographic ana-

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 W/iAEM 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]. W/iAEM for Windows 95
was released in 1997. The release of W/iAEM 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 W/iAEM 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. W/iAEM version
2 incorporated the GFLOW1 solution engine [Haitjema, 1995] and included new elements: inhomogeneities
and resistance line-sinks.
The current release, W/iAEM version 3, has additional support for raster base maps, such as USGS
digital raster graphics (DRG) and shapefiles (SHP). Also, W/iAEM has simple WHPA tools to include the
functionality of the EPA WHPA code. W/iAEM runs under Windows 7 and 10. W/iAEM is available for
download from the EPA Center for Exposure Assessment web site (https://www.epa.gov/exposure-assessment-
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.
W/iAEM can use any Cartesian coordinate system, such as Universal Transverse Mercator (UTM), State
Plane, or local site coordinates. The W/iAEM 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
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

Degrees of west longitude
126 120 114 108 102 96 90 B4 76 72 66
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. W/iAEM 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 W/iAEM 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 W/iAEM ; 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

EPA	wellhead analytic element model
United States
Environmental Protection
" WMmM
Principal Investigator/Solver: Dr. Henk Haitjema' Versjon 3 3 0 . August 30, 2015
Graphical User Interface and base map support: Dr. Vic Kelson2™"3, Mr. Mike Galvin4"™8,
Mr. Stephen Alberty4, Mr. Jim Evans4"01"5, Mr. Henry Helgen6, Dr. Mart Haitjema'rK""7
Project Officer: Dr. Stephen Kraemer*
'Haitjema.,Consulting, Bloomington. IN, VMtman Hydro, Bloomington, IN;
3Layne Christen sen, Bloomington, IN; 'Computer Sciences Corp., Athens, GA;
sSRA Athens, GA; ^Computer Sciences Corp., Duluth, MN; 'Skyport Systems, Mountain View, CA
sU.S. Environmental Protection Agency, Office of Research and Development
W/iAEM"™ is placed in the public domain, with some restrictions. The
software and documentation may be freely copied and distributed.
The graphical user interface and solver (GFLOW1) are copyrighted
by Haitjema Software, a subdivision of Haitjema Consulting, Inc., and
are open source (see the Artistic License in the file ALtxt). The
utilities "gzip/gunzip" and "tar" are free software and available from
WiA(w.gnu.org/directory/GNU/ (see thefileCOPYING.txt).
Figure 1.2: W7/AKM 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 (T.AM web data download area (wwwiepa,gov/exposure-assessment-models/wliaem2000)).
We also assume user familiarity with the Microsoft Windows operating system, and basic mouse operation.
Step 1 Start W/iAEM and familiarize yourself with the graphical user interface (GUI):.
From the Start button, go to All Programs, start WMEM , The "About W/iAEM" splash Screen
(See Figure 1.2) will appear briefly, and credits are given.
The workspace for WMEM 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.1). 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).
W/iAEM requires consistent units,, either meters and days or feet and days- The "Distance I'nits 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

WhAEM Version 3.3.1
_ |n|x||
~ | B
Edit Hew Element Model lools ffiindow Help
Figure 1.3: W/iAEM user interface.
| ^ Create New Project Database

Libraries ~ Documents ~
My Documents ~
WhAEM projects ~ Vincennes
if) | Search Vincennes
Organize ~ New folder

s -
H ©
1 WhAEM projects
) verifications
Documents library
Arrange by: Folder ~


I WinEdt9


; Public Documents
My Music


Public Music

.«£. Muse

^ Pictures

Select a file
H Videos

to preview.
& Kraemer, Stephen

, .orade_jre_usage

Jl .swt


fc. Contacts

III Desktop

it Downloads


New Database Wizard
- Project Description —
Glacial outwsh wellhead case study
-Distance units in Basemap Files —
(• Meters	C Feet
Units for Computations (Elevations, Properties, Pumping Rates) —
C Meters and Days	<• jFeet and Daysi
Create Database
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.

Project Settings
Project Description
Glacial outwsh wellhead case study
Currently Used Base Maps
c:\users\skraerner\documents\whaern projects\vincennes\vincennes_24k.tif
Remove Map
Add Map...
Intermediate Model Files
Base Filename (DOS)
| SI c: [OSDisk]
Distance units in Basemap Files Meters

Units for Computations:
a Documents
OWhAEM projects
! Vincennes
Change Units..
Figure 1.6: Project; settings.

C:\Users\SKRAEMER\Documents\WhAEM projects\Vlncennes\viribase.whm - WhAEM Version 3.3.1
Project Edit
~ & | i3
View Element Model Tools Window Help
x|^|t|ar| *|«|^±|h| %|1 ft|#| ¦ |n..\ \o.\m|E31 f 141

: : .V . ¦
i-i i	/ < -• , V.
a.. ..,./v:r n;'; -r-
¦¦'W -	¦'

^ . iy ¦
Click mouse button to select element
441,364.7 m 4,282,440.4 m
Figure 1.7:: The WMEM base map view of USGS DRG.
S3 |
t |
4% 1
Figure 1.8: The View smart icons.

J	I »|g|^|B '¦ | s?| £ |0 | B|Pt| ®|Q|B|B|
Driv6-in .Theater
Lincoln Kpr
C:\Users\SKRAEMER\Pocuments\WhAEM projects\Vincennes\vinbase.whm - WhAEM Version 33.1
| Eroject Edit View Element Model loots Window Help
I ~
Figure 1.9: Zoom-in view of city wellfield.

This page is blank

Chapter 2
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
Open (or restore to the view) the project database vinbase, and make a duplicate database "setbackl":
(Project , Make Duplicate Database, setbackl). W/iAEM 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' '.

Well Properties
(F1 for Help)
General Other I
C Head-Specified
(* Discharge-Specified
X coordinate
Y coordinate
452650 meters
4280665 meters
0 H"3/dav
1 feet
Figure 2.1: Well dialog box.
Well Properties
(F1 for Help)
Label	|Well 1
General Other |
Setback radius
I- Trace particles from well
Number of particles |
Starting elevation
r Choose Pathline Color
Pathline Color	Color...
200 feet
0 feet
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 ft setbacks.
Draw the wellhead protection area (Edit , Draw wellhead protection area). Name the protection
area "200 ft 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-t.he-ground realities. And as in this case, when multiple
capture zones are generated using various conceptual models, a "wellhead protection area" can be drawn

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.

This page is blank.

Chapter 3
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:
The first term represents the inflow due to aquifer recharge. Since the area ranges from 0 to ttR2 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:
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
^NnR2t + nnR2H & Qt

Isochrone t
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
t^RdivideN = Qt	(3-3)
so that
Rdivide = y	(3.4)
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:
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.

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 W/iAEM 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 setbackl.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 vincfrl 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 ft3/d , the areal
recharge to the water table N = 0.00276 ft/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 tt = 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 ft3/d, aquifer saturated thickness
H = 67.5 ft, 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.

• / OtiCAL DM 'vprivb^'n JmBaier } -
i^OirHl^|TR|C^FR (t=5 years)

Lincoln Mjrt
'¦fc. J - ^iewaee •
Figure 3.2: Calculated fixed radii for city wellfield based on recharge (R=6532 feet) and volumetric, methods
(R=3992 feet),

To determine this capture zone (e.g. by use of W/iAEM ), 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,f983](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.

[head, 2S-.ro ni
v* i n i
I fittid i 26,5 £ 
, I <# iri i
[In.J, !6.'J7wl
0	51)
	1	i	I	
[bj (£6.26- £S.?0) (Z6.S6-it5.01')
* ¦ 6H n — -
ffi ?Gffl
£ i i «« ( ?
* L ZD m
le] IL-JL Z6_S!
ht Q.I3 n»
i " !M «
,u> (t	>
\. ~\\J 90" f
_ j Direction oP
{i* ~C Jr rn b « # «n # n I
* Lyi
tB.07 n
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 vincfrl. From
the Menu, Project , Make Duplicate Database, and enter vinunil.
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 e veins.
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 ! foot
contour lines, and 5 feet for 10 foot contour lines [USGS, 19!)!) . The spatial accuracy of the topographic
lines are accurate approximately plus or minus 40 feet [US.GS, 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 ft); (2) Otter Pond (403 ft); and (3) the; 400 ft topographic line
crossing of the Wabash River.
Step 2
Use the uniform flow element in W/iAEM to fit the estimated uniform flow.


•^^HULvon S;i.


IflFM^iCnF.. E
.N.VI PO y.Ms-Jtiif
fLttAfr i
A .f!

Figure 3,4: W/iAEM view of test points used to calculate the uniform flow vector.

Model Settings
Aquifer | Contouring | Tracing | Solver]
-Aquifer Properties —
Base Elevation	|33Q
Hydraulic Conductivity
(F1 for Help)
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 ft. 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 ft, the thickness to 100 ft, the hydraulic conductivity to 350 ft/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 ft, the
maximum elevation of 430 ft, and a contour interval of 1 ft. 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 (403/t and utm 452325, 4277311) and the 400 ft 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.

Model Settings
Aquifer Contouring | Tracing] Solver]
Piezometric Contour Settings —
(F1 for Help)
Show Contours	<• Coarse C	Detailed C Custom
Minimum Contour	I390	feet
Maximum Contour	p30	feet
Contour Interval	Fj	feet
Figure: 3.6: Contouring settings dialog box.
Model Settings	(F1 for Help)
Aquifer | Contouring ] Tracing Solver |
-Iteration Settings —
Number of Iterations:	[2"
W View error log file after solve
F* View message log file after solve
Line-sinks with resistance and inhomogeneities with a change in base elevation
may result in non-linear equations that require several iterations for an accurate
solution. The solver employs" iterative refinement" to improve accuracy of the
solution. Three (3) iterations will suffice for most problems, but check the
"Runtime Message File" (Model>View Model Run Files) for errors at boundary
conditions. Errors should generally be less than 1 %. Iterations may be added
using the "Continue Solve" option on the Model menu.
Figure 3.7: Solver settings dialog box.



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 a.-terix 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 Messa.ge.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

Well Properties
(F1 for Help)
Label	|Well 1
General | Other |
<~ Head-Specified
(• Discharge-Specified
X coordinate
Y coordinate
4526E0 meters
4280665 meters
370000 ft~3/day
4 feet
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 ft3/d,. Enter a radius of
4 ft 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 W/iAEM 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 fe (that
is 5 ft 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(Jyyr). 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.

Figure 3.11: Capture zone for 5 years residence time based on t.racelines for well in uniform flow field. The
previous CFR solutions, volumetric and recharge, are also shown for scale.

A+4 46
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 ft). 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.

This page is blank.

Chapter 4
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 W/iAEM 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.

Plan view
q=Q/A =
Q0™ = Q
-k^ = -ki
Cross sectional area
Q„ = -
Figure 4.1: The magnitude of the discharge vector of the uniform flow field is calculated knowing the hydraulic
conductivity (k), the hydraulic gradient (?'), 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 (f).
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 Qa [L'2/T], and can be estimated from the hydraulic
gradient i [-] and the aquifer transmissivity kH (hydraulic conductivity k times saturated aquifer thickness
H) [L'2/T]. See Figure 4.1. The magnitude of the uniform flow rate is calculated as
Qo = kiH
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= —
where T is the time-of-travel [T] and Ta [T] is a reference time defined as:
T0 =

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
B=L1543/il	(4-4)
When 0.1 < T < 1, the CFR is given by
R = Ls[1.161 + ln(0.39 + f)]	(4.5)
where Ls is the distance from the well to the stagnation point down gradient from the well given by
Ls = Q/(2ttQ0)	(4.6)
and where the eccentricity 5 is given by
5 = L .,[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/2Q0  Project Settings > Add Map); to select the wellhead protection area
for display (View > Base Maps > Vector Graphics).

T - —
> i
> o.i & < 1
Ecentric circular capture
Centric circular capture
capture zone
Figure 4.2: Simplified delineation techniques for a well pumping at rate Q, in an ambient flow field Qa,
with aquifer saturated thickness H and porosity n, and time of travel capture zones of time T, after [Hai-
1 Calculated WHPA
F1 for Help

Selected Well
Time of Travel
1826.25 days
Aquifer Porosity

Aquifer Thickness
Ambient Flow

67.5 ft
0 ft"2/day


This can only be done for one well at atime. Based on the
data, a dimensionless travel time parameter T is calculated.
used to decide on a centric circle, excentric circle, or pathlines
for well in ambient flow.

Figure 4.3: Dialog box for calculated WHPA.

irrfplej^HPA (red)
\ * Kr
, ?V\ jr ' . \^r
¦ _ , tk'estpm^
CFR Volumetric 5 yr
-- sf

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.

1 Calculated WHPA
F1 for Help

Selected Well
Time of Travel


1826.25 days
Aquifer Porosity

Aquifer Thickness

67.5 ft
Ambient Flow

23.625 ft"2/day



This can only be done for one well at atime. Based on the
data a dimensionless travel time parameter T is calculated.
used to decide on a centric circle, excentric circle, or pathlines
for well in ambient flow.

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 ft, thickness at 1000 ft,
hydraulic conductivity at 350 ft/d, and porosity at 0.2), contouring (390 to 400 step 1 ft), and tracing
(1825 d).
Step 3 Add the wellfield as before (single pumping center well: utm — x = 452650, utm — y = 4280665,
discharge=370000 ft3/d, radius 4 ft). Add other properties, such as 20 tracking particles starting from
elevation 335 ft. 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 traeelines 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

Well in Uniform Flow
F1 for Help
T = 1.2822 (T > 1)
Direction of Flow
150 Degrees
The direction of flow is measured in counterclockwise
direction from the horizontal x-axis. For example, north is
90* and south is -90*
Make sure that the "trace particles" on the "other" tab of
the Well Properties dialog is set and a sufficient "number
of particles" are entered.
Figure 4.6: Dialog box simple well in uniform flow.
Stream Boundary
Do you want to use a Stream Boundary
<® No
Head at first point
Hydraulic Conductivity
0 ft
0 ft/day
Stream Boundary:
Explanation: You will define the (infinitely long) straight stream
boundary by entering 2 points on the map with the mouse cursor
(left click). The specified head will be applied to the first point you
enter. The head along the stream boundary may vary depending
on the orientation and strength of the uniform flow and the hydraulic
conductivity you entered. Both pathlines and head contours will be
drawn. Make sure that the "Show Contours" is selected on the
"Model>Settings>Contouring" tab and the proper contouring
parameters are selected.
No Boundary:
Only pathlines will be drawn to define the WHPA(s).
Figure 4.7: Dialog box simple stream boundary.

Figure 4.8: Simple WHPA including ambient flow traeelines compared to well in uniform flow envelope.

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

Multiple Wells WHPA
Ambient Flow
Aquifer Porosity
Aquifer Thickness
Direction of Flow
F1 for Help
Select Well(s)
23.625 ftA2/day
67 5 ft
150 Degrees
To ensure proper treatment of well interference, capture zones will
be represented by sets of particle traces from each well selected
Make sure that for each well selected the "Trace particles from
wells'1 on the well properties dialog is checked and a sufficient
"number of particles" are entered.
Figure 4.10: Dialog box Simple WHPAs with multiple wells.
Figure 4.11: Simple WHPAs including multiple wells in ambient flow.

Chapter 5
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 W/iAEM . 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
W/iAEM 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 W/iAEM 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 W/iAEM 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. W/iAEM 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 W/iAEM 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 W/iAEM — the head and average ground-water flow rate can be calculated at any

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 W/iAEM . 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,f989] and [Haitjema,f995].
Analytic elements are defined in W/iAEM 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 W/iAEM 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 W/iAEM 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
W/iAEM 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 W/iAEM 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 W/iAEM by selecting the Element menu, and selecting New and then Line-
5.2.1 Exercise: Study area and Line-sinks
Step 1 Locate the pumping center.
Create a new project database called hydro 1 with basemap units meters and computational units ft and d.
Add the regional series of USGS 7.5 minute raster basemaps
Vincennes_24K_utml6n27.tif ,Bicknell_24K_utml6n27.tif, Decker_24K_utml6n27.tif,
EMtCarmel_24K_utml6n27.tif, Frichton_24K_utml6n27.tif, Heathsville_24K_utml6n27.tif,
Iona_24K_utml6n27.tif, Lawrenceville_24K_utml6n27.tif, MonroeCity_24K_utml6n27.tif,

0akton_24K_utml6n27 .tif, Patoka_24K_utml6n27. tif, Plainville_24K_utml6n27 .tif,
Russellville_24K_utml6n27. tif, SandyHook_24K_utml6n27 .tif, St Franc isville_24K_utml 6n27 .tif,
Washington_24K_utml6n27.tif, Wheatland_24K_utml6n27.tif to the project, and view.
Position the pumping center as described earlier, and create the discharge specified well (Qw = 370,000 ft3/d,
radius = 4 ft, 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 ft, hydraulic conductivity 350 ft/d, porosity 0.2, thickness 500 ft).
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 ft 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 ft
contour and click the left mouse button. Type the contour elevation 400 ft 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 W/iAEM 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 ft
elevation and 400 ft 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 WabashROl.
Treat the linesinks as "far field" and the click the "along stream centerline" radio button. Next fill in the
Starting Head of 410 ft and Ending Head of 400 ft (See Figure 5.3). Click OK. Move the cursor to



Wabash River

£' < Jrj. .
White River (
¦ !	r/	: y ;;
Mh	^ ' "
/ o nVS i ' f "	. 380
/' r - " , • 380 . ,

Figure 5.1: Regional XJSGS DRG basemaps with hydro elevation labels.
* nl^
Figure 5.2: Elements smart icons: wells, linesinks, inhomogeneities.,

Linesink String Properties (F1 for Help)
Label	E
General |
(• Head-Specified
C Discharge-Specified
P" Treat as "far-field"
-Linesink Location
C Unknown
(• Along stream centerline
Along surface water boundary

Starting Head
41 ojlfeet
Ending Head f
Am ^et
Figure 5.3: Linesink properties dialog box.
center of the Wabash at the 410 ft 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 ft
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 W/iAEM based on a linear interpolation of the starting head and ending head along the line-sink
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 ft and 397.5 ft elevations
(estimated half way between the 400 ft and 395 ft 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 ft contour, and continue to the estimated location
of 397.5 ft.
Return to Wabash River centerline linesinks from 397.5 ft to 395 ft and between 395 ft and 390 ft, and
between 390 ft and 380 ft.
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

X I Sj V I
	 !	i 	j
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
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 ft/d) [Shedlock,
1980]. To do this, from the Menu, select Element > New > Inhomogeneity. Type in the recharge rate of
0.00276 ft/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 ft and maximum contour of 500 ft and interval 1 ft. 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 ft 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

iash River
Figure 5.5: Hydrography labels and line-sinks representing the Wabash River and White River.

Inhomogeneity Properties	(F1 for Help)
General |
I Change hydraulic conductivity from default
Conductivity	I	35q|| ft/day
p" Change recharge rate from default
Added Recharge Rate | 0.00276 ft/day
| Change porosity from default
Porosity	I	g~2
| Change base elevation from default
Base Elevation | 330
I Provi d e e sti m ate of ave rag e h e ad
Average Head | 830
Color... |	Inhomogeneity Color
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

>ash River
Figure 5.7: Hydrography labels, line-sinks, and inhomogeneity (recharge) element (blue rectangle).

Figure 5.8: 5 year capture zone including the influence of hydrologie boundaries.
~ +0.62
¦7 8S
~ +1.73
-0 95
~ -1 34
~ -0 12
~ +1.67
I5 58
Figure 5.9: Modeled head contours compared to field observations at monitoring wells.

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 W/iAEM 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 W/iAEM 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 W/iAEM 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.
W/iAEM 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 nof lowl. You may find it helpful to add the USGS DLG hypsography
(topography) maps (VIlhp06.bbm, VIlhp07.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.

Viefcc ditcl
Geology from H.H. Gray, i.i. *ayne and C.E. W[«r (1970)
Alluvial silt,
sand, and gravel
Wind-blown sand
and some silt
Wind-blo*n silt and
'fine sand and clay
Outwash gravel,
tend, and silt
Figure 5.10: Surtieial geology of the study region. jShedlock, 1980][Gray et al.,1970].

Alluvial silt, sand, and gravel
Geology nodified from H.H. Cray,
W.J. Wayne, and C.E. W|ar(1970)
Wind-blown silt, fine sand, and clay
0utwash gravel, sand, and silt
Glacial till, mostly clay
Bedrock, sandstone, and some shale
Figure 5.11: Generalized geologic cross section of the Wabash River Valley near the study region [Shed-
lock, 1980] [Gray et al.,1970].

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 W/iAEM in-
terface. These domains can have different hydraulic conductivities, aquifer bottom elevations, effective poro-
sities, and/or areal recharge due to precipitation. W/iAEM 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 Fl). See Figure 5.13. It may be useful to edit the vertices that make up the inhomogeneity

Drawing the inhomogeneity
Enter vertices.
Right-click on
last vertex closes
the domain.
Do not create sharp
Incorrect: line-sink
should not intersect
Correct: line-sink vertex is
on line-doublet.
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 W/iAEM to represent the outwash.
We will build an inhomogeneity domain to represent the outwash of the Wabash River with a hydraulic
conductivity of 350 ft/d, in contrast to the rock uplands area with a hydraulic conductivity of 10 ft/d.
The first thing to do is to make a duplicate database called inhomol. Now go to Model , then Settings,
and enter a new hydraulic conductivity as the background value of 10 ft/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
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 VIlhp06.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 ft/d. Add the recharge
N = 0.00276 ft/d. Click OK. Use the hypsography basemaps (vilhpf06, vilhpf07) 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.

Inhomogeneity Properties	(F1 for Help)
Label	| outwash
[7 Change hydraulic conductivity from default
Conductivity	|	350 ft/day
[~ Change recharge rate from default
Added Recharge Rate | 0 00276 ft/day
(~~ Change porosity from default
Porosity	|	Q2
Change base elevation from default
Base Elevation	I	33Q feet
Provide estimate of average head
Average Head	|	630 'eet
Cclci... :|	Inhomogeneity Color
OK |	Cancel
Figure 5.14: Dialog window for entering inhomogeneity properties.
The vertices of linesinks and inhomogeneit.ies should overlap where these elements cross (See online
help). Therefore, refine the linesink vertices such that the linesinks are completely inside and outside the
inhomogeneit.ies. Also note that inhomogeneit.ies 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).

Figure: 5.15: Layout of linesinks, wellfleld, 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.



Uf t
A 5
d 6 °~K
! V

o = wq_
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 (6) divided by the hydraulic conductivity (kc) of the layer.
c = f	(5-2)
The vertical specific discharge qz[L/T} is equal to the head in the aquifer (
Liresink String Properties (F1 for Help)
(* Head-Specified
C Discharge-Specified
(Treat as "far-field1 ]
Linesink Location
C Unknown
C Along stream centerline
* Along surface water boundary
Starting Head
Ending Head
400 feet
397 5 feet
14 3 days
SI 5 feet
1 feet
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
Calculate the resistance of the bottom materials in the Wabash River.
Step 1
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 ft/d, and of average thickness of 1 ft. Let's work
with the value of kc = 0.07 ft/d and 5 = 1 ft, 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 LTSGS
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 inhomo 1. whm. Create a duplicate
database resistl. 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.

Figure 5.19: 5-year capture zones including resistance of Wabash River (c = 14.3 d, <5 = 1 /t, width B
=615 ft, 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 ft/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 (JFLOW [Haitjema, 19951 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

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]. W/iAEM 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,
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).

rfek trailer
Era'faf ' \
'Pa'M r/
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,1991]).

Chapter 6
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
by surface water when the distance of the well to the surface water feature is greater than 4\JQT/ (nHri),
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:
(a)	The well or well field is within a distance of 2y/kfl/kvH 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.
(b)	The width of the capture zone in a two-dimensional model is less than 2v/kh/kvH1 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.

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 4y/Tc, where T (ft2/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 (6) divided by its hydraulic conductivity (k) or c = b/k.
(c)	A hydrological boundary (stream or lake) is present within a distance of 4\/Tc from the well,
where T (ft2/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-t ravel 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-t ravel 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

that the corresponding calculated capture zones were progressively more realistic. The emphasis of the
W/iAEM 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. W/iAEM 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!

[Aquaveo, 2017] Aquaveo (2017).	GMS groundwater modeling system v 10.3.5.
[Bakker, 2007] Bakker, M. (2007).	TimML: a multiaquifer analytic element model.
[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.

[Haitjema, 2006] Haitjema, H. (2006). The role of hand calculations in ground water flow modeling. Ground
Water, 44(6): 786-79 f.
[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.
[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(l-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 Al,
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

[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-3J38/.
[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 Cliff's, 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.
[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.

This page is blank.

Appendix A
Conversion Factors

To obtain
meter (m)





2.589988 x 106

inches/year {in/yr)
2.2815 x 10"4

6.9541 x 10"5


gallons per day(gpd) / ft2
Volume rate
ft3 / second(s)







Appendix B
W/iAEM Base Maps and the Internet
B.l 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. W/iAEM can import and display vector maps created from binary base maps (*.bbm) or
shapefiles (*.shp).
B.l.l 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.l. 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.l). 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 (VII), containing hydrography (HY).
Table B.l: BBV
index by Quads
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).

WhAEM2000 BBM Files: U.S.
Binary base map (BBM) files are vectorized graphics derived from U.S. Geological Survey (USGS)
digital line graphs. BBM files can be imported directly into WhAEM2000 software and serve as the
geographic basis for groundwater modeling. To find BBM data for a specific geographic area, click on
the appropriate state below. A more detailed graphic of the state you selected will display.
Alternatively, you can choose from an alphabetical text list of available states.

Alabama | Arizona | Arkansas | California (North) | California (South) | Colorado | Connecticut-
Massachusetts. Rhode Island | Florida | Georgia | Idaho | Illinois | Indiana | Iowa [ Kansas | Kentucky |
Louisiana | Maine | Maryland. Delaware. Distr. of Columbia | Michigan (Lower Peninsula) | Michigan
(Upper Peninsula) | Minnesota | Mississippi | Missouri | Montana | Nebraska | Nevada | New Hampshire-
Vermont | New Jersey | New Mexico | New York | North Carolina | North Dakota | Ohio | Oklahoma [
Oregon | Pennsylvania | South Carolina | South Dakota | Tennessee | Texas (West) | Texas (East) | Utah |
Virginia | Washington [ West Virginia | Wisconsin | Wyoming
Figure B.l: The EPA binary base map server.
C A Secure https://earthexplorer.usgs.gov
4. Search Results
If you selected more than one data set to search, use
the dropdown to see the search results for each specific
data set.
Note: You must be logged in to download and order
Show Result Controls "
Dot a Set
Click here to export your results » C?
DLG 1:100k

« First * Previous | 1 r i Next > Last »
Displaying 1 -10 of 66 O
Entity ID:1J072207WTR02931-462850
1 Show Metadata Overlay:Roads
Revision Date:01-JAN-86
Entity ID: 1J072207WHP02928-462880
2 Show Metadata Overlay:Hypsography
Revision Date:01-JAN-86
'* /lO
1 At®
Entity ID:1J072207ETR02926-462836
Entity ID:1J072207WHY02929-462828
3 Show Metadata Overlay.Hydrography
Revision Date:01-JAN-86
i' 39'W N, 087* 30' 33* W) Options Overlays

* University

t 1
„ \%.*A
& \
* ^ c
©* \


/ ©

Map data ®2017 Google 500 m
, $ w
Terms of Ure Report a map er
• Google map is not for purchase or for download; it is to be used as a guide for reference and search purposes only.
Figure B.2: USGS BarthExplorer web page.

Table B.2: BBM by Type
Category Name
roads and trails
miscellaneous transportation
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 W/iAEM . 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 W/iAEM project). Maps relevant for source water protection may include land
use/land cover and soils maps.
Shapefiles can be added to the W/iAEM 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. W/iAEM 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. W/iAEM 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 W/iAEM 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).

Figure B.3: USGS 7.5 minute quad map showing topographic elevation contours for Vincennes, IN.

The National Map I
G A Secure I https://viewer.na:ionalmao.gov/basic/
Figure 15.1: USGS Hational 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 W/iAEM 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, XAl )27 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 waterdatauusgs.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).

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