United States	Office of Research and	EPA/600/R-00/022
Environmental Protection	Development	April 2000
Agency	Washington DC 20460
v»EPA Working with WhAEM
Source Water
Assessment for a
Glacial Outwash Wellfield,
Vincennes, Indiana

revised March 2003
Working with W/iAEM2000
Source Water Assessment for a Glacial Outwash
Wellfield, Vincennes, Indiana
Stephen R. Kraemer
U.S. Environmental Protection Agency
Ecosystems Research Division
Athens, GA 30605
Hcnk M. Haitjcma
Haitjema Consulting. Inc.
Bloomington, IN 47401
Victor A. Kelson
Wittman Hydro Planning Associates
Bloomington, IN 47404
Office of Research and Development
U.S. Environmental Protection Agency
Washington, DC 20460

The U.S. Environment,al Protection Agency (EPA) through its Office of Research and Development (ORD)
partially funded and managed the research described here under contracts (P.O. 8L-0425-NTSA, P.O. 0L-
2294-NASA) to Haitjema (Consulting, Inc. It. has been subjected to the Agency's peer and administrative
review and has been approved for publication as an EPA document.
Mention of trade names or commercial products does not constitute endorsement or recommendation for
This document was typeset using WinEdtS of Aleksander Sirnonic (www.winedt.org) and MikTgX(including
Larlj?JX2e) software of Christian Shenk (www.rniktex.org). Tp]X. is a trademark of the American Mathematical

The purpose of this document is to introduce through a case study the use of the ground water geohy-
drology computer program WAAEM for Microsoft Windows (32-bit), or WMEM2000. W/iAEM2000 is a
public domain, ground-water flow model designed to facilitate capture zone delineation and protection area
mapping in support of the state's and tribe's 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 scries of progressively more complex representations of the well field tapping
a glacial outwash aquifer for the city of Vincennes, Indiana. W/i,AEM2000 provides an interactive computer
environment for design of protection areas based on radius methods, well in uniform flow solutions, and
geohydrologic modeling methods. Protection areas are designed and overlaid upon US Geological Survey
Digital Line Graph (DLG) or other electronic base maps. Base maps are available for download from the
EPA Center for Exposure Assessment Modeling web site. Geohydrologic modeling for steady pumping wells,
including the influence of hydrological boundaries, such as rivers, recharge, no-flow boundaries, and inho-
mogeneity 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/iAEM2000
has on-line help and tutorials.

The National Exposure Research Laboratory's Ecosystems Research Division (ERD) in Athens, Georgia,
conducts process, modeling, and field research to assess the exposure risks of humans and ecosystems to
chemical and non-chemical stressors. This research provides data, modeling tools, and technical support to
EPA Program and Regional Offices, state and local governments, and other customers, enabling achievement
of Agency and ORE) strategic goals in protection of human health and the environment.
ERD research includes studies of the behavior of contaminants, nutrients, and biota in environmental
systems, and the development of mathematical models to assess the response of aquatic systems, watersheds,
and landscapes to stresses from natural and anthropogenic sources. ERD field and laboratory studies support
process research, model development, testing and validation, and characterize variability and prediction
Leading-edge computational technologies are developed to integrate core sciences into multi-media (air,
surface water, ground water, soil, sediment, biota), multi-stressor, and multi-scale (organism, population,
community, ecosystem; field site, watershed, regional, national, global) modeling systems that provide pre-
dictive capabilities for complex environmental exposures.
Exposure models are distributed and supported via the EPA Center for Exposure Assessment Modeling
The development and release of the Wellhead Analytic Element Model (WAAEM2000) 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. WAAEM2000 is a general purpose ground water hydrology tool based on the innovative
modeling technique known as the analytic element method and a standard Windows graphical user interface.
WMEM2000 assists in the delineation of the area contributing source water to public water supply wells.
The development of W7iAEM2000 has been supported by the EPA Office of Research and Development and
the EPA Office of Ground Water and Drinking Water.
Rosemaric C. Russo, Ph.D.
Ecosystems Research Division
Athens, Georgia

1.1	Source Water Protection Programs				 						1
1.2	Modeling Philosophy		2
1.3	The Emergence of the WHPA Model . 					 						4
1.4	The Evolution of WftARM		5
1.5	The Base Map 			 										6
1.6	Exercise: Exploring the USGS quad map Vincennes, Indiana 		7
2.1 Setbacks		11
2.1.1 Exercise: Vincennes, Indiana . 					 						11
3.1	Calculated Fixed Radius 		19
3.1.1 Exercise: Vincennes, Indiana		20
3.2	Well in Uniform Flow 		21
3.2.1 Exercise: Well in Uniform Flow, Vincennes, Indiana 					23
4.1	Introduction to the Analytic Element Method		31
4.2	Rivers and Line-sinks		32
4.2.1 Exercise: Vincennes and Line-sinks . 								32
4.3	Geologic Contacts and No-Flow Elements		40
4.3.1 Exercise: Vincennes and no flow boundaries) 					42
4.4	Aquifer Inhornogeneities		44
4.4.1 Exercise: Vincennes and Inhornogeneities 			 						44
4.5	Riverbed Resistance		47
4.5.1 Exercise: Vincennes and Resistance Linesinks		48
4.6	Discussion		49
5.1	Good Modeling Practice		53
5.1.1	Rules of Thumb 		53
5.1.2	Simplified methods 		55
5.2	Conclusion 		56

A Conversion Factors	59
B W/tAEM2000 Binary Base Maps and the Internet 61
B.l Getting BBMs from the EPA Web Server	 61
B.2 Creating BBMs from USGS DLGs maps from the USGS Web Server	 61

Thanks to all participants in the testing of the code, including the participants in the EPA Regional Training
Courses (R5 Chicago, May, 1999),(R3 Philadelphia, June, 1999, August 2001), (RIO Seattle, January, 2000),
(HQ, March, 2000), (R8 Denver, February 2001), (R9 San Francisco, April 2002). Thanks to our 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 ERD-Athens for his review. Thanks for support from the EPA
Office of Ground Water and Drinking Water and Mr. Jim Hamilton, Ms. Denise Coutlakis, and Dr. Marilyn
Ginsberg. The development of W/iAEM benefited from comments from the former EPA Ground Water
Protection Technical Forum.


Chapter 1
The objective of this document is to introduce the geohydrology computer model W/iAEM2000. The pre-
sentation will cover both basic model operations and modeling practice. We will use a source water area
assessment associated with the wellfiekl of the city of Vincennes, Indiana, to emphasize a step-wise and pro-
gressive modeling strategy. The case study is not an approved component of the Indiana wellhead protection
1.1 Source Water Protection Programs
Ground wafer is a valued source of public drinking water in many parts of the country. Not. only Is the ground
water resource renewable if properly managed, its quality can be excellent, due to the natural cleansing
capabilities of biologically active soil cover and aquifer media. It has come to the public attention in the
past few decades that ground water wells can be threatened by over-drafting and pollution. Source water
protection is based on the idea of protecting public water wells from contamination up front, and working
to maintain the quality of the resource. Source water protection means taking positive steps to manage
potential sources of contaminants and doing contingency planning for the future by determining alternate
sources of drinking water. The Safe Drinking Water Act has required the USEPA to develop a number of
programs that involve the source water protection concept for implementation at the State level, such as the
Wellhead Protection Program and the Source Water Assessment Program [USEPA,2003].
The main steps of the source water protection process involves the assessment of the area contributing
water to the well or wellfield, a survey of potential contaminant sources within this area, and an evaluation
of the susceptibility of the well to these contaminants. This includes the possibility of contaminant release
and the likelihood of transport through the soil and aquifer to the well screen.
The designation of the wellhead protection area is then a commitment by the community to source
area management. Delineation of the wellhead protection area is often a compromise between scientific and
technical understanding of geohydrology and contaminant transport, and practical implementation for public
The EPA Office of Ground Water and Drinking Water 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. The methods applying the criteria
range from the simple to the complex, including setback mapping, calculated radius, and geohydrologic
modeling (see Table Li). This document will review selected criteria and their methods using WL4EM2000

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









residence time




assimilative capability

through examination of the Vincennes, IN case study. The emphasis will be on the residence time criteria
and its methods. The evaluation of assimilative capability is an active area of research and will not be
discussed in detail in this document.
1.2 Modeling Philosophy
The philosophy presented below is for deterministic mathematical modeling, and a more complete discussion
can be found in [Konikow and Bredehoeft, 1992]. The reader is encouraged to explore an alternative modeling
approach based on stochastic theory [Vassolo et al., 1998] , [Feyen et ah. 2001].
Modeling of flow and transport through geohydrologic systems is a complex undertaking. Our definition of
good modeling practice revolves around step-wise 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 by 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 of constant head). Such a conceptual
model is a severe abstraction of a well near a stream or lake boundary. The actual meandering surface wafer
boundary is replaced by an infinitely long straight line, heads underneath it are assumed equal to the surface
water elevation, aquifer properties are assumed uniform and homogeneous throughout the flow domain, and
lastly the ground-water flow is assumed to be purely horizontal. The virtue of this conceptual model is the
ease with which a solution can be obtained (Thiern'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 when
compared to the complex reality of a partially penetrating well near a meandering stream with a siltv 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 may escalate when striving 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 may suggest an increasingly
accurate representation of the ground-water flow regime, the reality may well be different. Increased data
requirements and the associated data uncertainty, as well as the difficulty to interpret the results of complex
(multi-parameter models) can defeat the theoretical improvement in model realism. This is not only true
for ground-water flow modeling. Albert Einstein is known 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 not only take the complexity of
the geohydrology into account, but also the objective of the modeling exercise. In the context of wellhead
protection the purpose of the modeling is often the delineation of a "time-of-travel capture zone" for one or
more drinking water production wells. The delineation of the time-of-travel capture zone, for the purpose
of wellhead protection, 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 are staggering, while 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 to say the least. 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.
The essence of good 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 or not a conceptual model is not
too simple, is to gradually increase its complexity and investigate the effect on the purpose of the modeling;
the time of travel capture zone in this case. A very good example is the classical dilemma of modeling
two-dimensional flow or three-dimensional flow. Most modelers feel inadequate when using a horizontal flow
model (based on the Dupuit assumption), when realizing that the actual world is clearly three-dimensional.
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 aquifer heterogeneities, varying recharge rates, resistance to flow in or out of surface water
features, etc. can only 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 what we are
modeling for. Once progressively more complex conceptual models do not change the time-of-travel capture
zone in a meaningful way, the modeling may be considered complete. Sometimes, it may 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 may be obtained by making some assumptions that are known
to lead to a larger capture zone than the one that may result from a more detailed conceptual model
trading uncertainty into conservative/protective assumptions. Also, modeling experience will lead to general
"rules-of-thunib" in formulation of the conceptual model. These topics will be visited in the discussion in
the final chapter.
1.3 The Emergence of the WHPA Model
Nearly all ground water modeling projects are constrained by a finite budget and a finite amount of time.
This is particularly true for time-of-travel capture zone delineation in the context of wellhead protection.
Most States are confronted with the task of developing wellhead protection programs for thousands of public
drinking water supply systems in only a few years time. Conducting an extensive ground water modeling
campaign for each individual drinking water well (or well field) is out the question, both in view of 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
1.	arbitrary fixed radius method.
2.	calculated fixed radius method.
3.	simplified variable shapes.
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 latter (#7) all-out modeling exercise may completely satisfy the groundwater hydrologist,
the more approximate approaches are offered to get the job done. The rationale behind this thinking is
that it is better to define a less then 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, constraint of
the wellhead protection program.
The first two methods lead to circular wellhead protection zones and often over-protect some areas and
not protect other areas near the well.
The well in uniform flow method fakes ambient, ground-water How into account (approximately) and
offers a significantly better approximation of the actual time-of-travel capture zone. That approach has
been implemented in a computer code WUPA, made public by the USEPA [Blanford and ITuyakorn,1993].
In addition, WHPA allows 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 uses particle
tracking to draw bundles of streamlines toward the well, including the boundary streamlines that define the
capture zone for the well. By drawing streamlines with a constant ground water residence time, from starting
point to the well, a time-of-travel capture zone is delineated. WHPA is menu driven, has a windows-like
graphical user interface under DOS, and requires minimal computer resources.

1.4 The Evolution of W/iAEM
VV/tAEM (pronounced warn) offers an evolutionary but significant step. In addition to facilitating the sim-
plified techniques (e.g., fixed radius and uniform flow), the modeling environment allows the 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 (WHPA), but can be of arbitrary shape, further improving the realism of the conceptual
model. On the other hand, WMEM currently lacks varying aquifer bottom elevations, three-dimensional
elements, multi-aquifer flow, transient flow, etc., which are supported in full-feature ground water models.
While this limits the applicability of W/iAEM2000, it also makes the code easier to operate and thus easier to
learn. In addition, the relatively simple conceptual model of W/iAEM requires less input data, thus further
simplifying its use. Taken together, time-of-travel capture zone delineation with W/iAEM is considerably
cheaper than when using a full feature ground water flow model such as e.g. , MODFLOW [McDonald and
Harbaugh, 1988]. Obviously, the price for this is a less sophisticated (potentially less accurate) definition of
the time-of-travel capture zone. W/iAEM clearly offers more realism than the first three approaches in the
list, but also requires more work (more time and more money). In some cases this may be justified, in some
cases one of the simpler methods may be adequate, and in some other cases a more powerful ground water
model must be employed. It is the responsibility of the hvdrologist to determine whether or not 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 explained, should be based on technical merits, while taking into account the
practical limitations associated with the timely implementation of the entire wellhead protection program.
The release of the Capture Zone Analytic Element 'Model (CZAEM) [Strack et al.,1994] introduced a new
analytical solution technique to the wellhead 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, line
doublets (double layers) for geologic contacts | Strack, 19891, [Haitjema,1995|. 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 assmumption
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 to approximate a three-
dimensional flow field by a two-dimensional one. Implementation of the analytic element method in CZAEM
includes wells, line-sinks, uniform flow, and circular recharge elements, and is operated from a DOS command
line. CZAEM has sophisticated routines for drawing time-of-travel based capture zone envelopes and time
and source water sub-zones [Bakker and Strack,1996].
A graphical user interface (GUI) and geographic analytic element preprocessor, GAEP [Kelson et ah,
1993], was added to solution engine CZAEM. in the first release of W/iAEM [Haitjerna et. al, 1995], [Haitjema
et. al, 1994], This version ran under Windows 3.1. The development of GAEP opened the door for conceptual
modeling, 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 [USACOE,2002].

WMRM for Windows 95 was released in 1997.
The release of W/iAEM version 1 (April 2000) included an improved windows graphical user interface,
the use of U.S. Geological Survey (USGS) digital line graphs (DLGs) as the base map, a new solution
engine ModAEM [Kelson,f998] that included arbitrarily shaped no-flow boundaries. The release ofWMEM
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.
The most recent release, W/iAEM2000 version 2, uses the GPLOW1 solution engine [Haitjema, 1995],
and includes new elements: inhornogeneities and resistance linesinks. W&AEM2000 runs under Windows
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 the 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.
In addition, we need to define the origin of our ground-water model coordinate system in a standard way so
that we can overlap maps and spatial data from different sources.
W/iAEM.2000 can use any Cartesian coordinate system, including universal transverse mercator (UTM),
State Plane, local site coordinates, etc. The graphical user interface (GUT) for WftAKM'JOOO provides both
pre- and post- processing in the base map coordinates, while the computational engine operates in local
rectangular coordinates for reasons of numerical accuracy, but this is hidden from the user.
We will build our base maps using the USGS Digital Line Graph series. DLG maps include lines for
hydrography (rivers, lakes, etc), roads, hypsography (topography), 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). DLG maps are available for download on the Web (See Appendix 13.2).
The DLG 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 f to 60,
where by zones 10 to 19 cover the conterminous U.S. Figure 1.1. In each zone, coordinates are measured
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. Caution: Base maps do not smoothly join across UTM zones!
For performance reasons, we have developed a compressed binary form of the DLG called BUM (bi-
nary base map). W/tAEM2000 reads BBM. files from the disk each time it redraws a map on the screen.
WMEM2000 has tools to convert DLG to BBM. The EPA Center for Exposure Assessment (CEAM) oper-
ates a web server for download of BBMs. The CEAM web server has a graphical index, and when accessed
with a browser, you can select specific maps for your 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 also
variable; states that are elongate from east to west, such as New York, use a transverse Mercator projection,
while states that are elongate from north to south, such as California, use a Lambert conformal projection.
Maps in any projection stored in digital exchange format (DXF) may be translated to BBM format using
the Tools .

126° 120° 114° 108° 102° 96° 90° 84° 78° 72° 66°
Figure 1.1: UTM zones of the conterminous U.S.
1.6 Exercise: Exploring the USGS quad map Vincennes, Indiana
We will use a source water area delineation for the well field of Vincennes, Indiana, as a case study throughout
this document. The case study will be used to: (1) familiarize the reader with the basic operations of
W/iAEM2000; and (2) illustrate some basic ground-water modeling concepts.
Vincennes is the seat of Knox County and home to about 40,000 people. The city is located in the alluvial
flood plain of the Wabash River in the southwestern quadrant of Indiana. See Figure 1.2. The floodplain of
the Wabash River is characterized by valley fill from glacial outwash, and consists of coarse sand and fine
gravel on top of underlying bedrock. The area receives about 40 inches per year of precipitation. The areal
recharge to the aquifer is estimated to be 12.1 inches per year, and the average hydraulic conductivity of
the local aquifer is 350 feet per day, and an effective porosity is 0.2 [Shedlock, 1980]. The wellfield pumped
about 2.75 MGD (million gallons per day) in 1977, and about 3.5 MGD in 1999. You will learn more details
about the wellfield as you progress through the case study.
In this exercise, we will place the Vincennes wellfield in its geographic context. You will locate the well-
field on the digital U.S. Geological Survey 7.5 minute topographic map (vincenne.tif) using the DLG-Viewer
[USGS,2002]. The DLGV32 Pro program is available for download from the USGS webpage
(http://mcmcweb.er.usgs.gov/drc/dlgv32pro/). The 7.5 minute digital raster graphic (DRG) for the Vin-
cennes quad (vincenne.tif) is available at the CEAM web download area
(http: / / www.epa.gov/ceampubl/gwater/whaem/index.htm).
We will assume basic familiarity with the Microsoft Windows operating system, and basic mouse opera-
Step 1 Start the USGS DLG-View program and open the file vincenne.tif.
How is this accomplished? From Windows, use the mouse to Start, Programs, US Geological Sur-
vey, dlgv32. From the Menu bar, click on File, Open as New, look in the WhAEM folder for the file

Figure 1.2: Location of Vincennes, Indiana
vincenne.tif, click to select, and Open.
Size your view to fill the monitor screen. Figure 1.3.
Step 2 Test out the Smart Icons buttons: Zoom In, Zoom Out, Full View, Zoom Tool, Pan.
You will use the smart icons (see Figure 1.4) of the USGS DLG-View program to reposition the view.
Zoom In to the City of Vincennes, and find the city wellfield, Figure 1.5. Note the UTM coordinates of
the pumping center reported in the Status Message bar in the lower right corner of the Workspace (452650
m, 4280665 m). Pan up the Wabash River to the north and locate where the 400 foot topographic line
crosses the river.
Minimize the view to continue to the next chapter, or exit the program by clicking File and Exit.

Figure 1.3: USGS DLG-View of Vincennes, Indiana,
Out	Recenter
ig; Global Mapper - REGISTERED

File View Tools Searcl


^~|cnn]| [;^ |>-?| « |
] Atlas Shader _v|
Zoom Full Zoom
In	View to box
Figure 1.4: USGS DLG-View icons.

'avel Pit
'W»:UOD Hjfv.
' O'N'E^L '
pumping center
.	* V; \
rnpM IJTM («4U89! «.S'«S '51 X «'IS7r»'N BT Jl" 19IB"
Figure 1.5: Zoom in view of Vincennes wellfield.

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 by the States [Merkle et
al.,1996]. For example, at this time, Indiana has 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. And the setback is not likely protective from a gasoline spill containing the
additive MTBE (methy tertiary butyl ether). MTBE resists degradation in many subsurface environments.
2.1.1 Exercise: Vinceniies. Indiana
Step 1 Start W/iAEM2000 and familiarize yourself with the graphical user interface (GUI).
From the Start button, go to Programs, start WftAEM2000 for Windows. The "About. WMEM2000
for Windows" splash screen will appear briefly, and credits are given. Figure 2.1. Select the '"Start WhAEM"
option on the welcome screen. Figure 2.2.
The workspace for W/iAEM.2000 has standard windows features: a Menu bar on top: smart icons bar
underneath the menu bar; status bar on the bottom of the window. Figure 2.3.
Create a new project database.
the Menu Bar, drop) down the Project options. Select New Database, and select your project
folder (e.g., where you stored your bbms) on the hard drive, and type vinbase Into the dialog box for the
file. Figure 2.4. Click Open.
The "New Database Wizard" will prompt you for a project description and the units associated with
the input data. 2.5. WMEM2000 requires consistent units, either meters and days or feet and days. The
Step 2


(I	Office of Research and Development
Office of Ground Water and Drinking Water
Principal Investigator and Solver, Dr. Honk Haityama, Haitjema Consulting, www.haitjema.com
Graphical User Interface, VieKaSson, Wittman Hydro Planning Assoc., www.wittmanhydro.com
Mike Gaivin	Mart Haitjema
Project Officer, Steve Kraemer EPA Ecosystems Research Division, www.epaJgov/athens
This program uses the programs 'gzip' and 'tar', which were published by the
Free Software Foundation and are covered by the GNU Public License. See the
file COPYING for details or visit www.gnu.org.
Figure 2.1: WMFM2000 splash screen.
Welcome to WhAEM
WhAEM is a groundwater flow model specifically designed to assist in
we 11 h e ad p rote cti o ri are a d el ineati o n. If yo u are a fi rst-ti rn e WhAE M user,
we invite you to first enter the WhAEM tutorial by pressing the 'Tutorial'
button below.
Start the WhAEM Tutorial

Start WhAEM
Don't show this dialog box in the future
Figure 2.2: Start W/iAE'M2000.
WhAEM for Windows
Project Edit View EJement Model Tools Window Help
~ »l
Figure 2.3: WMEM2000 menu bar.

Create New Project Database
a c? m-
Figure 2.4: Create new database.
"Distance Units in the Basemap Files" will be dictated by your source files. Select meters for the "Distance
Units in Basemaps Files" by clicking on the radial 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. After you have made
these selections, click on Create Database.
Note: 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 a 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, which works in consistent units.
The next window facilitates the definition of base maps (in binary format). Figure 2.6. The window
"Currently Used Base maps" will be empty at this point. Click on Add BBM to display the available
*.bbm files. Do a shift click to select all of the bbm files. Then click Open. The list of bbm files will be
displayed in the box "Currently Used Basemaps". Click OK and wait a few seconds. The maps will be
loaded and plotted on the project workspace. Figure 2.7. See Appendix B.l for information about the USGS
naming convention for these files.
From the Menu Bar, drop down the View options. Go to Basemap to view the available layers. Play
with toggling these layers on and off and observe the change in the view.
From the smart icons bar, test the zoom in, zoom out, zoom to extents, zoom to window icons, scrolling,
and re-centering. Figure 2.8.
Step 3 In this step you will locate a well in the center of the Vincennes wellfield, and draw setback
protection area.
Make a duplicate database "setbackl". (Project , Make Duplicate Database, setbackl). W/iAEM2000
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 to the file. So it is useful
to make a new project database anytime you start a new exercise.
Minimize the W/iAEM2000 window. Restore (or open) the dlgv32 program and explore the Vincennes

New Database Wizard
Project Description
NewWhAEM Project
Distance units in Basemap Files
f* Meters	C Feet
Units for Computations (Elevations, Properties.. Pumping Rates)
C Meters and Days	(• [Feet and Days;
Create Database
Figure 2M New database wizard.
Project Settings
Project Description
IJ W 3 s i agrera
Currently Used Base Maps
Remove BBM
Intermediate Model Files
Base Filename (DOS)
~ c: [DRIVE-C]
ac a
\l Program Files
Figure 2.6: Project settings.

¦JM.I	u.iana
SSm* 1 T-i hwTTbT- l«'K'l 011.1 ^!S|G1 »H» ~I-I	whbasecdr
Figure 2.7: The base map view.
Zoom out
Zoom to
Zoom in
Zoom to

Figure 2.8: The View smart, icons.

Well Properties	(F1 foi Help)
jwell field
General Other |
Setback radius	|
P T race particles from well
Number of particles	|
Starting elevation	I
200 feet
0 feet
Figure 2.9: Well properties dialog box.
Well Properties	(F1 foi Help)
Label	| well field
General ) Other l
(• ID ischarge-S pecif ied
X coordinate
Y coordinate



452650	meters
4280665	meters
0	ft"3/day
1	feet
Figure 2.10: Well properties dialog box.
DRG map. Zoom to the Vincennes wellfield. Position the cursor near the pumping center and read off
the UTM coordinates. Note that 452650m, 4280665m is about, the pumping center. Minimize or exit the
program dlgv32.
Restore the W/iAEM2000 window. Zoom in to the city of Vincennes.
To create a well, from the Menu select Element , then select. New and click on Well. Move the cursor
to the approximate location of the well, based on your exploration of the Vincennes DRG map, and click the
left, mouse button. A marker for the location of the well is added to the base map and a. dialog box appears.
Type in a. name for a. label, e.g. * 'well field''. Under the Other folder enter the "Setback radius" of
200feet. Figure 2.9.
Under the General folder, re-adjust, the well location to be ut.m location x = 452650m., y = 4280665m by
typing into the X and Y boxes. Figure 2.10.
Click OK to add the well to the database.
Zoom in on the well and note the setback is plotted as a. dashed circle around the well. Figure 2.11.

3C:\Program Files\WHAEM\setback1 whm-WhAEM lor Windows Release Candidate 1.0.0
~ Ib|b| B| | | | | +|8|^|	D|S:U|Q|h|E3| «|->|t|4|»|
Figure 2.11: 200 foot setback.
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 dotted outline;
click with the right mouse to close the final segment.
The purpose of drawing this area by hand is separate the delineation of the capture zone (based on
hydraulics/hydrology) from the definition of the wellhead protection area (based on other considerations,
such as parcel or zoning boundaries). The draw utility allows the wellhead manager to fine tune the shape
of the wellhead protection area to correspond with on-the-ground realities. Also, if multiple capture zones
are generated using various conceptual models, a "wellhead protection area" may be drawn that envelops
the various capture zones for maximum protection. More on this later.


Chapter 3
The justification of the residence time criteria as protective of ground water source water protection 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
Vincennes case study: (1) calculated fixed radius; and (2) well in uniform flow7.
3.1 Calculated Fixed Radius
The fixed radius can be calculated based on a simple two dimensional static water balance analysis, assuming
negligible ambient flow in the aquifer. If we assume radial flow toward a well in an aquifer with a constant
saturated thickness ll(rri). the cylindrical boundary of radius R(rri) is delineated by an isochrone of residence
time I (days), which means that any water particle that, enters the cylinder or is present in the cylinder will
travel no longer than I days before being pumped up by the well (See Figure 3.1). 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 ??.(-). A water balance for the period I yields:
The first term represents the inflow due to aquifer recharge, 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 for the pumping period. The radius R can be expressed as:
NttRH + mrR2H = Qt
When t becomes infinitely large, the radius R represents the complete capture zone:

N ^
isochrone t
Figure 3.1: Water balance for radial flow to a well in a domain bounded by an isochrone of residence time t.
Use of this equation is called the "recharge method" [USEPA,1993].
If the term NttI becomes small, due to a small value of I or N or both, equation 3.2 reduces to:
Use of this equation is called the "volumetric method" [USEPA,1993].
The equations 3.2 through 3.4 only apply if the Dupuit assumption applies, that is 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 may become important, and the approximation
may lead to an underestimation of the capture zone isochrone radius R. Such a non-conservative result may
be avoided by replacing the saturated thickness H in Equation 3.4 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 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 which extends in the upgradient direction well beyond the circular
representation of that isochrone, as follows from 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 discussion so far assumed a constant saturated thickness //, as 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 II in Equation 3.4.
3.1.1 Exercise: Vincennes, Indiana
Step 1 Use WMEJV12000 to plot the Calculated Fixed Radius protection areas based on the recharge
method and the volumetric method.

Figure 3.2: Calculated fixed radii for Vincennes wellfield based on recharge and volumetric methods.
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.
For the recharge method, use equation 3.3, and assume the well field pumps Q = 370, 000ft3/day , and
the areal recharge to the water table is N = 0.00276ft/da,y(12.1in/yr), and calculate R. If you do not have
a portable calculator handy, we suggest you use the Microsoft Windows calculator (Start > Programs >
Accessories > 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. 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 previous chapter).
For the volumetric method, use equation 3.4, assume Q = 370,000ft3/day, aquifer saturated thickness
H = 67.5ft, the aquifer porosity is n = 0.2, and the pumping period is five years (1826.25 days). Draw the
wellhead protection area.
Compare your results to Figure 3.1.1.
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, may be 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 into the direction of uniform

To determine this capture zone (e.g. by use of WMEM20W) it is necessary to determine:
•	direction of the ambient flow.
•	the hydraulic gradient.
•	the aquifer transrnissivity.
•	the pumping rate of the well.
•	the desired maximum residence time.
Given three water levels (average, static), it is possible to estimate the direction of the ambient flow and
the hydraulic gradient [Heath, 1983]:
1.	Identify the well that has the intermediate wafer level.
2.	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.
3.	Draw a straight line between the intermediate well and the point identified in step 2. This line represents
a segment of the water level contour along which the total head is the same as that in the intermediate
4.	Draw a line perpendicular to the wrater level contour and through the well with the lowest (or highest)
head. This indicates the direction of ground water movement in an isotropic aquifer.
5.	Divide the difference between the head of the well and that of the contour by the distance between the
well and the contour. This gives the hydraulic gradient.
The reader is referred to EPA's on-line three-point gradient calculator created by Dr. Jim Weaver at
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 and sinks in between
[Cole, Silliman, 1996].
Be sure to do a search of the scientific literature (US Geological Survey, State Geological Survey) to find
any studies of the regional aquifer system and the regional hydraulic gradient.
The aquifer transrnissivity should be determined by a pumping test. This is an expensive procedure,
and transrnissivity is known to be scale dependent. Perhaps transrnissivity measurements are available from
locations not too far from the well.
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 phe-
nomenon. It is 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 to 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 protective area. In other words, use the "wellhead area" drawing tool to draw
an envelope around the set of capture zone realizations.

Will l
[hm i , Zfi.t 0 m )
Wtl[ I
( httid i 26.2 6 »li f
tfiri j
Ihtid, ?fi 07 m ]
0 35 50
L_	I	1	
(6 ) (26.26 - 26.20] ( 26 26 - it5 0? I
?fi ?B m
* L.=26 20 m "	T T<
te] 26. t -25.07
/ D irettion »f
{i* jf mowttatnl
EG. OT 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, Vincennes, Indiana
Step 1 Determine the magnitude and direction of a uniform flow field near the Vincennes pumping center.
First we will give our project a new name. Open the database from the previous exercise. 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 USGS topographic
map and the USGS DLGV32 software and the Vincennes DRG. The "lakes and ponds" in this area are
actually abandoned quarries, and the removed overburden gives a window into the elevation of the water
table. We assume these pits act like large piezometers, and field observations confirm that their water levels
do in fact respond to recharge events.
A note of caution is deserved 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 Vincennes maps was photo-revised in 1989. The elevation values are accurate
say plus or minus a foot. The spatial accuracy of the topographic lines are accurate; say plus or minus 5 feet.
These uncertainties are acceptable for this level of analysis.
Use the method of Heath or the OnSite Calculator discussed in the previous section,with the following
three points of known head: (1) Mirror Lake (407 feet); (2) Otter Pond (403 feet); and (3) the 400 foot
topographic line crossing of the Wabash River.
Step 2
Use the uniform flow element in WM.EM2000 to fit the estimated uniform flow.
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

Figure 3.4: W/iAEM2000 view of test points used to calculate uniform flow vector.
Lake. Enter in the elevation 407feet. Enter the gradient (for our case about 0.001). Enter the angle (for
our case about 150 degrees).
Step 3 Enter the aquifer and contouring properties and test points.
Erom the Menu, select Model, then Settings, and then select the Aquifer tab. Set the base elevation
to 330ft', the thickness to 100ft, the hydraulic conductivity to 350ft/day, and leave the porosity as
0.2. 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 390ft, the
maximum elevation of 430ft, and a contour interval of 1 ft.
Click on OK.
To test our conceptual solution, we may introduce so-called "test points" at locations of known head, e.g.
Otter Pond ( 103/7. and utm 452325, 4277311) and the 400ft topographic intersection futm 454315, 4283085)
with the Wabash River, by selecting Model , then add test points, locate, 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. 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 a test point).
Next, we will 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 may set. the number of iterations to 2
in order to achieve some "iterative refinement." of the solution.

Model Settings
Aquifer] j Contouring] Tracing ] Solver |
Aquifer Properties
Base Elevation	I330
Hydraulic Conductivity
(F1 for Help)
Figure 3.5: Aquifer settings dialog box.
Model Settings
Aquifer Contouring | Tracing | Solver |
Piezometric Contour Settings
W Show Contours (* Coarse C Detailed
Minimum Contour	I39Q
Maximum Contour	I43Q
Contour Interval	[-jj
(F1 for Help)
Figure 3.6: Contouring settings dialog box.


Figure 3.7: Solve icons.
1^ C:\Proqram Files\WHAEM\vinuni.whm - WhAEM for Windows Beta 3.1.1

Project Edit View EJement Model Tools Window Help
~ |b|b| a| < I sh* I I \|;?l
Well Pioperties	(F1 foi Help)
| well field
General | Other I
t* Discharge-Specified
X coordinate
Y coordinate
452650 meters
4280665 meters
370000 ft 3/day
4 feet
Figure 3.9: Well dialog box.
Locate the Vincennes pumping center either1 by double clicking on the existing well or creating a new one
(Element > New > Well). Enter in the well pumping rate of 370,000ft3/day. Enter a radius of 4ft
to represent our pumping center. The Vincennes wellftekl consists of a local cluster of at least 5 wells; we
will represent with a single pumping center. Too small an effective radius and the pumping center will draw
the wafer table below the aquifer base. The radius value does not effect, the flow solution, but will provide
the starting point for particle of particles during capture zone delineation. A radius of 4 feet will work fine.
Figure 3.9.
Enter OK.
Step 6 Plot a capture zone based on 5 years of residence time.
Erom the well properties box, click on the "Other" tab, on the well, and check the box Trace particles
from Well in the Other tab, and set the Number of pathlines to 20. The well is assumed fully penetrating
so we have an option to set the starting elevation of the trace particles. Set the Starting elevaton to 335feet
(that is 5feet 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 days (5 years). 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 , and 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.10. 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 parameterization.
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. Click on the well, enter the 5 year CFR volumetric previously calculated (R=3992 ft). Change

Files\WHAEM\vinuni1 .whm - WhAEM for Windows Release Candidate 1.0.0
Project Edit View EJement Model Tools Window Help
PlBrlal a		 :\x\e.\«\ a|n=l etNaM
Click mouse button to select element
flow 5 yr
453.G55.2m 4,280.289.9 m /y
Figure 3.10: Capture zone based on 5 years residence time in uniform flow field.
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. Figure 3.11.

Figure 3.11: Capture zone based on 5 years residence time based on CFR volumetric and t.racelines for well
in uniform flow field (zero gradient).


Chapter 4
In the previous analyses we have ignored the presence of hydrological boundaries. Instead, we lumped their
effects near the well into a uniform flow field. 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 power of WMEM2000. We will progressively refine the complexity of our hydrogeological
conceptual model using the Vincennes 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).
4.1 Introduction to the Analytic Element Method
WMEM2000 uses the analytical element method to solve ground-water flow subject to hydrological boundary
conditions. Analytic elements are mathematical functions that represent hydrogeologic features in a ground-
water flow model. For instance, a meandering creek or river may be represented by a string of straight line
elements, which approximates the geometry of the stream. In W/iAEM2000, each of these line elements is
associated with a "line-sink function", which represents mathematically the ground-water inflow, or outflow,
into the stream section. The sink density of the line-sink equals the groundwater inflow into the stream
section. Each stream section is likely to have a different ground-water inflow rate and thus a different sink
density. These sink densities are not known in advance. What is considered known are the heads underneath
the stream sections. For instance, in WMEM2000 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/iAEM2000 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
Other examples of analytic elements in W/iAEM2000 are wells (point sinks in two-dimensions), line
doublets or "double layers" to represent aquifer inhornogeneities, 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&AEM2000 	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 area! elements. This is the principle of superposition.

Note: the ground-water velocity is calculated by adding "derivative functions" for the analytic elements, which
have been obtained analytically and are programmed into W/iAEM2000. For details on the mathematical functions
and on the mathematical procedures to enable the superposition of solutions of both confined and unconfined flow, the
reader is referred to |Strack,1989] and |fiaitjema,1995].
Analytic elements are defined in W/&AEM2000 by clicking the Element menu, and selecting New.
4.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/iAEM2000 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/tAEM2000 this is handled by specifying a piezometric head at the
center of the line-sink, which 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 known 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 WMEM2000 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/iAEM2000 a lake receives or loses ground
wafer 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 streams that have a width that is large
compared to the aquifer thickness, the realism of the line-sink representation may be improved by positioning
line-sinks along both sides of the river.
The number of line-sinks required to represent a stream or lake depends on two factors: (1) the need to
follow 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 which may draw 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/iAEM2000 by selecting the Element menu, and selecting New and then
4.2.1 Exercise: Vincermes and Line-sinks
Step 1 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).
Open the database file vinbase.whm. Create a duplicate database hydrol. 'Minimize the W&AEM2000
window. Open the USGS 7.5 minute topographic map for Vincennes with the dlgv32 program. Find the
location where the 400ft topographic contour crosses the Wabash River north of the city. Note below the

UTM location of the 400/# contour:
utm — x(existing) =		(4.1)
utm — y(norihing) =		(4.2)
Now minimize the dlgv32 window, and restore the W/iAEM2000 window. Zoom In to the location of the
400ft contour. From the Menu, select Edit, and then Add text, move the text cursor to the proper location
(using feedback on the UTM location of the cursor from the bottom status line) on the stream and click the
left mouse button. Type the contour elevation 400 feet in the dialog box, check the Hydrography Label
option, and click OK. The label appears in the map and the text cursor reappears. You would need to repeat
this procedure to add all the rest of the hydrography elevations, i.e. elevations where topography contour
lines cross the surface water features.
Note: You could also eyeball the location of the hydro label on the WMEM2000 basemap by either having a paper
topo sheet in hand, or have the W/iAEM2000 and dlgv32 windows open at the same time on your computer screen.
Another option is to use a digitizing tablet and a paper topographic map, and saving the hydro labels to an ASCII
text file, as described below.
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 WMEM2000 program
linearly interpolates between the starting head and ending head for the heads at vertices of the linesink
We have created a text file for you with elevations to assist in the completion of your basemap. The label
file must be a comma delimited ASCII text file (filename.TX) in which each line has the following format:
utm — x, utm — y, h, label	(4.3)
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. Select the file
HydroJExample. tx. In order to view the hydrography labels on the basemap, you may need to toggle them
on (View > show hydrography labels).
Step 2 Represent the Wabash River with a string of Line-sinks.
At this point, It is recommended thai you set. your window to frame the region Including the 395 foot
elevation and 400 foot 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 Figure 4.1. You will be prompted
for the name of the line-sink string about to be entered. Let's start with Wabash River, hence type Wabash
River. Treat the linesinks as far field. Next fill in the Starting Head of 400 ft and ending head of 395 ft.
Figure 4.2. Click OK. Move the cursor to the east bank of the Wabash at the 400 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. You may elect to use shorter segments close to the wellfield. After clicking on the
395 ft location, click the right mouse button to end the linesink entry. Click OK.

+ 4-
¦t- t-
— I
Figure 4.1: Elements smart icons.
Linesink String Properties (F1 for Help]
Label	|Wabash River
General ]
(• Head-Specified
P Discharge-Specified
W iTreat as "far-field1:
Starting Head
Ending Head
| 335
1 o


Figure 4.2: Linesink properties dialog box.

Note: ft is recommended to enter line-sink strings along streams going down-stream (into 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 WAAEM2000 based on a linear interpolation of the starting head and ending head along the
line-sink string.
We recommend you represent the other side of the Wabash River with a string of line-sinks, particularly
near the wellfield. Start again at the of 400 foot contour, and continue to an estimated location of 397.5 feet
(halfway between 400 and 395 ft.
Step 3 Surround the wellfield with hvdrologic (line-sink) boundaries and solve.
Position the pumping center as described in the previous Chapter, and create the discharge specified well
(Qw = 370, 000ft3/d, radius = 4ft).
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 will include the perennial streams, or streams that are flowing all year round.
In the Vincennes wellfield area, those include City Ditch, Mantle Ditch, Kelso Creek, Ditch 2, Ditch 3, and
the Wabash River, of course [Shedlock, 1980J. See Figure 4.3. Perennial streams show up as solid blue lines
in USGS topographic maps, and dashed if ephemeral.
Create line-sinks representation of the perennial streams. Add the Vincennes well-field with the same
properties as before.
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.
You may change the location or properties of wells and line-sinks as follows. Click on the well or a
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 icon (the third
to the right of the printer icon). Acknowledge the warning. Notice the two added vertices (dots) on either
side of the vertex you clicked on. Repeat this for one or two other vertices of the Wabash River line-sink
string near the well.
Now let's add areal recharge based on the USGS suggested value of 12.1 inches per year (0.00276 ft/day)
[Shedlock, 1980]. To do this, from the Menu, select Element , New, Inhomogeneity. Type in the recharge
rate of 0.00276 ft/day. Click OK. Use the mouse and click to define the boundary vertices. Surround the
linesink elements with an inhomogeneity polygon. Wre will further the discussion of the inhomogeneity
elements later in this chapter.
Your model space should look something like Figure 4.5.
From the Menu, select Model, and reset your Settings, including Aquifer and Contouring parameters
as in previous exercises.

Figure 4.3: Perennial streams in the Vincennes area, after [Shedlock, 1980].

Inhomogeneity Properties
(F1 for Help)
Label	|areal recharge 12.1 in peryr
General |
|~~ [Change hydraulic conductivitiJ from defaults
Conductivity	|	~ ft/day
Added Recharge Rate | 0.0027E ft/day
Color... | |
I- Change base elevation from default
Base Elevation	I	q feet
I- Provide estimate of average head
Average Head	I	~ feet
Figure 4.4: Inhomogeneity element properties box.
Step 4 Solve and plot the 5 year capture zone for the well field.
In Model , Settings, select the Tracing tab and check the Compute Particle Paths. Double click
on the well and check the box Trace Pathlines Emanating From Well. Select Number of Pathlines
20. Click OK and then Solve. Draw the wellhead protection area. Figure 4.6.
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 5 Check the solution for integrity.
Select View Results Table from the Model menu, click on View and on Linesinks. You are looking
at the vertices of the line-sink strings. Look further to the right in the table to ensure that the "Specified
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 linesink string to the observed base
flow of a stream. The observation may 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 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.
It is useful to compare the model predicted heads to observations of heads at monitoring wells. Fortu-
nately, the USGS has reported a series of measured water levels in the Vincennes 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. The view will show the test point locations as triangles,
with the size and direction related to the difference between model head and observed head. Figure 4.7.
Click on the triangle to report the error on the bottom status bar. Notice we have some large errors in the
model to the south. The conceptual model will be refined in the next section.

Figure 4.5: Hydrography labels, recharge element, and line-sinks representing the Wabash River between
elevations 395ft and 400/t and perennial streams and ditches.

Figure 4.6: 5 year capture zones including the influence of hydrologic boundaries.
|^C:\Progtam Files\WHAEM\hydio0.whm - WhAEM for Windows Release Candidate 1.0,0
Project Edit View Element Model Iools Window Help
QlsLgj a] I I I I	vk?k:l^l Blw:| stlalislsl
Figure 4.7: Model head contours compared to field observations at monitoring wells.

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
triangles that point both up and down, whereby the differences between model head and observed head are
not too large. 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 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 may also indicate too low a model
transmissivity or too high an areal recharge rate (or both). When the actual aquifer is significantly stratified
or, in fact, consists of several vertically stacked interconnected aquifers, discrepancies between observed and
calculated heads may 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, may this trend be reversed.
Heads that are too high or too low may also be the result of improper representation of surface waters.
For instance, the headwaters of a creek may artificially raise the ground water table in WMEM2000 by
infiltrating large amounts of water (which these headwaters don't have to give). Such stream sections and
their line-sinks should be removed from the WMEM2000 model particularly if they show up as ephemeral
or intermittent steams on USGS topo maps.
4.3 Geologic Contacts and No-Flow Elements
The geology of the Vincennes study area deserves a closer look at this point as we refine the conceptual model.
The area is underlain by unconsolidated sediments shown in Figure 4.8. The sediments were deposited by
glacial processes in the Wabash River valley during the Pleistocene Epoch arid were reworked by wind
and water. The outwash pinches out a the margins or the river valley and generally rests directly on the
underlying bedrock, as shown in Figure 4.9. A couple of rock "islands" outcrop just south of the wellfield.
Ideally, we should represent the high permeable rock and till explicitly in our model. As a first step,
we will represent in W&AEM2000 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. This no-flow boundaries should be extended into the far field to ensure complete effect in the
near field.
WL4EM2000 conceptual models with and without these no-flow boundaries may be seen 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.

Geology from H.H. Gray, W.j. Wayne and C.E. Wier (1970)
AI luvial si 1t,
sand, and gravel
Wind-blown sand
and some silt
Wind-bI own silt and
fine sand and clay
Outwasti gravel ,
sand, and silt
Figure 4,8: Siirficial geology near Vincennes, Indiana, [Shedlock,1980]

alluvial silt, sand, and gravel
wind-blown silt, fine sand, and clay
outwash gravel, sand, and silt
glacial till, mostly clay
bedrock, sandstone,
and some shale
Figure 4.9: Generalized geologic cross section of the Wabash River Valley near Vincennes, Indiana [Shed-
lock, 1980]
4.3.1 Exercise: Vincennes and no flow boundaries)
Step 1 Create the no-flow boundary at the topographic contact between outwash and rock.
First, create a duplicate database as noflowl. From the menu, select Element , New, and Horizontal
Barrier. Name the feature "Outcrop" and click OK. Use the mouse create vertices along the boundary (left
click) and end (right click). Do the same for the island outcrops south of the wellfiekl. See Figure 4.10.
Remove the line-sinks that intersect or are to the east of the no-flow elements (click on string, then press
delete key). The link-sink elements to the east, of the no-flow elements are now hydrologically separated
from the wellfiekl and are irrelevant to the local capture zone solution. You should avoid crossing line-sink
elements and no-flow elements to avoid numerical difficulties.
Step 2 Solve and check the calibration.
Step 3 Plot the 5 year capture zone.
Double click on the wellfiekl, and release around 20 pathlines. Check the Model Settings, and Solve.
In View set the Pathlines to Show time tics to every 1 year. Draw the wellhead protection area.
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 wellfiekl (compare Figure 4.6 and Figure 4.12). If in reality some
water enters the outwash from the till and rock formations to the east, the capture zone will tend to extend
further to the east than for the case the rock and till are considered impermeable. This makes it reasonable to
explore solutions for the problem where the till and rock formations are given a finite hydraulic conductivity,

Figure 1.ID: Superposition of no-flow elements along the. geologic contact inferred from, the; topographic
Figure 4.11: Hydraulic head contours and test points, including the no-flow elements.

Select	Click mouse button to select element	454,025.4 m | 4,278,536.0 m >
Figure 4.12: 5 year capture zone solution including the influence of hydrologic and no-flow boundaries.
as we shall see in the next section.
4.4 Aquifer Inhomogeneities
Domains with differing properties are enclosed within polygons defining the inhomogeneity using the W/iAEM2000
interface. These domains may have a different, hydraulic conductivity, aquifer bottom elevation, effective
porosity, and/or areal recharge due to precipitation. W/iAEM2000 uses mathematical features called line-
doublets to simulate the polygon segments. Domains may be nested, but should not overlap. You should
avoid sharp corners when designing the domains. A line-sink should be either inside or outside the inhomo-
geneity domain. Extra vertices are needed near high capacity wells. See Figure 4.13.
4.4.1 Exercise: Vincennes and Inhomogeneities
Step 1 Create inhomogeneities using W/iAEM2000 to represent the outwash and rock outcrops.
We will build an inhomogeneity domain to represent the outwash of the Wabash River of hydraulic
conductivity 350 ft./d, in contrast, to the rock uplands area, of finite hydraulic conductivity, say 10 ft./d.
The first, thing to do is to open the vinbase project, database, and make duplicate database inhomol.
Now go to Model , then Settings, and re-enter the aquifer properties (Figure 3.5), but. enter a. new hydraulic
conductivity to a. background value of 10 ft./d. You should also re-enter the data, on the wellfield (Figure
To build the outwash inhomogeneity, select. Element , New, then Inhomogeneity. You will be pre-
sented with the opportunity to parameterize the inhomogeneity. Figure 4.14. Click the box to change
hydraulic condutivity from the default.. Enter the new hydraulic conductivity for the outwash domain

Drawing the inhomogeneity
Enter vertices.
Right-click on
last vertex closes
the domain.
Do not create sharp
corners!	x
Incorrect: line-sink
should not intersect
Correct: line-sink vertex is
on line-doublet.
Figure 4.13: Drawing inhomogeneities using the mouse.
Inhomogeneity Properties
(F1 for Help)
Label	|outwash
[~ [Change hydraulic conductivity from default!
Added Recharge Rate
350 ft/day
0.00276 ft/day
\~ Change base elevation from default
Base Elevation	|	33^ feet
I- Provide estimate of average head
Average Head	I	430 feet
I n horn o _properti es. eps
Figure 4.14: Dialog window for entering inhomogeneity properties.

Figure 4,15:: Layout of linesinks, wellfield, and aquifer Ijjhomogefteity domain.
k = 350 ft/d. Label the inhomogeneity. Check that the recharge rate Is correct- Click OK. Use the hypsog-
raphy basemaps (vilhpf06, vilhpf07) to guide your drawing the polygon, following the sharp transition from
outwash to rock. Figure 4.15.
Step 2 Enter the linesinks representing the near field and far field.
We will want to re-enter the linesinks to represent the Wabash River, and in addition, represent the White
River to the east and south. As before, bring in the: already prepared text label file by Tools > Import >
Text label file > Hydro_Example ,tx.: The far, field linesinks are entered in a very coarse representation of
the rivers. The elevations in the text label file were acquired using a USGS digital elevation model (DEM) and
ArcView. The near field linesinks representing the center channel of the Wabash River between elevations
400 ft and 39.5 ft should be entered as in a previous exercise. See Figure 4.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. A, capture zone solution is
shown in Figure 4.16. Notice that the tracelines now extend to the north under the Wabash River. The
interactions of the aquifer and the river deems further examination, as we shall do in the next section.
Advanced Use: When only the recharge differs, inhomogeneity domains may have arbitrarily large line-doublet sides,
the domain may overlap any other feature, including inhomogeneity domains and linesinks, and the domain may have
any shape, including sharp corners. When only porosity jumps, but not hydraulic conductivity or base elevation, the
inhomogeneity domains may not overlap other inhomogeneity domains, but the domains can take any shape, including
sharp corners and large sides, and the domains may intersect linesinks.

Figure 4.16: 5 year capture zone solution including the influence of aquifer inhomogeneit.ies.
4.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 may be 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 4.17. Resistance c[T] is equal to the
thickness of the layer divided by the hydraulic conductivity of the layer.
The vertical specific discharge qz[L/T] is equal to the head in the aquifer minus the head in the stream
divided by the resistance.
The discharge of the stream segment a[L3/T] is equal to the specific discharge times the width.
What values might one expect of the resistance parameter? Major rivers with sandy bottoms might have
resistance between 0 and 1 day. Small streams with silty bottoms might have resistance in the range between
1 and 10 days. Stream or lakes in till formations (overlying a deeper aquifer) might have resistance up to
100 days. Resistance greater than 100 days means almost no interaction with the aquifer. Resistance for a
specific site is dependent, on both the hydraulic conductivity of the riverbed and the thickness of the unit.
a = wqz

\ A
U* t d
\f \
_ 8
5 C~kc
a = wq:
Figure 4.17: Conceptual model of stream/lake resistance.
Discharge to resistance line-sinks in W/iAEM2000 is calculated given a layer of given resistance and
thickness, and effective width. Line-sink effective width w[L\ is calculated based on a characteristic leakage
length A[L] and the stream width B[L\. The characteristic leakage length is equal to the square root of the
product, of aquifer transmissivity (hydraulic conductivity times thickness) times resistance:
A = VkHc	(4.7)
There are some rules of thumb on setting the effective line-sink width parameter, given the stream width
and the leakage length.
Put line-sinks along both sides
of the stream of width w:
w = A	A < 0.1 B
w = Atanh(£>/2A)	0.1B < X < 2B
w = B/2	A > 2B
Put line-sinks along the center
of the stream of width w:
w = B	A > 2 B
4.5.1 Exercise: Vincennes and 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 Vincennes 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 foot..
Let's work with the value of kc = 0.07ft/d and 5 = 1 ft, and calculate the resistance using equation 4.4.
Step 2 Calculate the characteristic leakage length for the Wabash River.

Assume an aquifer thickness of II = 70ft and aquifer hydraulic conductivity of 350ft/d. Use equation
4.7 to calculate the characteristic length.
Step 3 Calculate the effective stream width for line-sink representation of the Wabash River.
Estimate the width B(ft) of the Wabash River near the pumping center using the USGS dlgv32 software
and the digital map vincennes.tif. Calculate an effective width w based on the table knowing B and
A. The Microsoft Windows calulator, under Accessories, can be set to scientific mode in the View options,
giving access to the hyperbolic tangent.
Step 4 Enter resistance linesinks to represent the Wabash River near field.
Let's build off of the previous database. Open the project database inhomol.whm. Create a duplicate
database resist 1. In the view, toggle on the hydrography labels. We are going to want to create resistance
linesinks on both sides of the Wabash river in the near field (the area of influence of the pumping well),
and we want to represent the far field Wabash on either side of the near field with single linesinks down the
center. In order to parameterize the heads in the near field, we suggest you query the vertices of the Wabash
river linesink for the elevations in the near field after a solve, and add those elevations to the basemap as
hydro text labels. You can then delete the Wabash river linesinks and build the new ones. We suggest you
create individual linesinks between each of the near field elevation labels, on both sides of the river channel,
starting from highest elevation to lowest. Enter the resistance, thickness, and width parameters as calculated
above. See Figure 4.18. Enter a string of far field (no resistance) linesinks down the center of the Wabash
River channel from elevation label 400 ft to the start of the near field linesinks. Do the same for the the
Wabash River to the west of the near field linesinks to the 395 ft contour.
Step 5 Solve the model and plot the 5 year capture zone.
4.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 fifty foot per day 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 may 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 may be
printed using the Printer setup and Print options on the Project menu. You may also export the graphics
in DXF format or export the piezometric contour and path line data in SURFER format for further post
processing. DXF maps can then be converted to binary base map format using the Tools and Convert
Haitjema, and Kelson investigated the influence of complexity of the conceptual model on the capture
zones for Vincennes [Haitjema, 1995]. They used the analytic element, code GFLOW [Haitjema,1995] to
include the influence of the low hydraulic conductivity of the rock outcrops. They used the finite difference
code MODFLOW [McDonald and Harbaugh, 1988] to represent localized three-dimensional flow, anisotropv,

398.25	^98
38	398.52
Wabash RjVer
F«u» 4.18: Dial
°S box for entiy of
distance iinestote
Wabash River
FL°w boundaries
1 on ZONES
Figure 4.19: 5-year capture zones
B - 615 ft, and w = 282 ft
deluding resistance of Wabash River (c
14-3 clays, S =

and resistance between the river and the aquifer. In carrying out this modeling, progressively more complex
conceptual models were tested. It appeared 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 4.20, "composite capture zones" are depicted, based on
the many capture zones generated by ITaitjema and Kelson. The capture zones were produced in GFLOW
using only 2D flow, but including resistance to flow from or into Wabash River.
Based on the exercises in this guide for the Vincennes wellfiekl, it appears that the 5 year capture zone
solution is insensitive to the exact value of the hydraulic conductivity of the rock outcrops — a no-flow
boundary (assuming the rock is impermeable) captures the essence. However, we did learn that the 5 year
capture zone is very sensitive to the resistance of the Wabash River channel — our capture zone extended
all the way to the north side of the river when resistance was given a conservative (high end) value. More
investment, in characterizing the riverbed resistance through field observations and modeling simulations
might be justified. Alternatively, one may assume the worst and accept the capture zone to the north side
of the river as realistic.
Remember, the case study for the Vincennes wellfield was for illustration purposes, and is not to be
misunderstood as the approved State source water assessment.
New some other complexities. Wc have not introduced the concept of horizontal an i sot ropy in the
hydraulic conductivity distribution. Horizontal anisotropv is known to influence the capture zone [Cole,
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],[Podgorncy, Ritzi, 1997]. WMEM2000 should
not be used in aquifers with a highly anisotropic hydraulic conductivity distribution.
We 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
represent the breakthrough time of the mean concent,ration of a conservative tracer at the well. If you wanted
to represent, the first arrival of the tracer in your delineations, and thus capture the initial breakthrough of
a potential contaminant, you would need to incorporate macroscale dispersion
(http://www.epa.gov/atliens/learii2model/part-two/onsite/conc.htm). Characterizing dispersion in capture zone
models is a topic of active research [Maas, 1999].
We have limited our guidance to manual model calibration. The reader is referred to [Hill, 1998] for
a discussion of automated methods of parameter estimation and model calibration. Automated non-linear
parameter estimation is implemented in UCODE
(http:/ /water.usgs.gov/software/ucode.html) and PEST
There are other professional analytic element modeling packages, including WTNFLOW, TWODAN,
GFLOW, and MLAEM. The reader is encouraged to follow developments of ArcElow which is based on the
solution engine SPLIT by Igor Jankovic, State University of New York-Buffalo, and the ARCVIEW interface
by the Minnesota, Department of Health [Blum, 2000]. Three-dimensional flow AE models include PhreFlow
by Randal Barnes of the University of Minnesota, and 3DFlow by David Steward at Kansas State University.
EPA maintains a web page with links to the above analytic element ground water models
(www.epa.gov/athens /'software/whaem/press.html).

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

Chapter 5
5.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 may be drawn
from comparing simple and complex conceptual models. In many cases the validity of certain assumptions
have already been investigated by others and are 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. We will, however, provide a few general rules of thumb that may be of help In making
some initial decisions regarding conceptual models, Figure 5.1 [ITaitjema,2002].
5.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 feature is greater than 4yQT / (tt tin), 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/kh/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 2*Jkh/kvH, where kh and
kv are the horizontal and vertical hydraulic conductivity, respectively, and where II is the average
saturated aquifer thickness.
Note that in the event that the (2D) capture zone width is in the same order than the aquifer thickness,
the three-dimensional capture zone is likely a slender three-dimensional stream tube that may occur
at some depth below the water table. A conservative (protective) time-of-travel capture zone may be
obtained by using a wellhead protection zone that is obtained from a horizontal flow model with a
confined aquifer of the same thickness as the length of the well screen.

Figure 5.1: Decision tree for complexity of capture zone delineation method, after [Uaitjerna.2002]. The
simplified techniques are discussed in section 5.1.2. Geohydro modeling involves solving the regional ground
water flow differential equations, including boundary conditions.
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 they are referred to as aquitards. Aquitards divide the subsurface into
more than one aquifer, forming multi-aquifer systems. In case of a multi-aquifer system a multi-aquifer
model (quasi three-dimensional model) may be necessary if all of the following circumstances apply:
(a)	The well occurs in a multi-aquifer zone and is not screened in all aquifers.
(b)	The length of the capture zone, based on flow in the aquifer(s) accessed by the well, is in the order
of (or smaller) than 4 \/Tc:, where T (ft2/day) is the transmissivity of the aquifer and c (days)
is the largest resistance of one of the adjacent aquitards. Resistance is equal to the aquitard
thickness divided by its hydraulic conductivity c = b/k.
(c)	A hvdrological boundary (stream or lake)is present within a distance of 4yTe from the well, where
T ( f t2/day) is the transmissivity of the aquifer and c (days) is the largest resistance of one of the
adj acent 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 tirne-of-travel capture zone will occur if interaction with
adjacent aquifers is ignored. If the direction of flow in the various aquifers are expected to be not too
different, the single aquifer capture zone may offer a conservative tirne-of-travel capture zone.
4.	Aquifer Stratification. An aquifer is considered stratified if the hydraulic conductivities of the
various geological formations (layers) are within the same order of magnitude. Aquifer 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 conductivity
and inversely proportional to their porosity. This can lead to substantially different tirne-of-travel
capture zone lengths in the various layers. If it is considered important to include this in the definition
of a wellhead protection zone, a tirne-of-travel capture zone may be produced 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
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" that may 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
long term trend. Accurate delineation for transient non-community wells may require special attention.
And projections into the future should include additional safety factors.
5.1.2 Simplified methods
The material presented in this section was sponsored by the Indiana Department of Environment,al Manage-
ment [Ceric, 2000J, [Haitjema,2002|.
How does one choose a simplified method? This decision depends on several parameters, including the
magnitude and direction of the ambient flow near the well or well field, which while important, is challenging
to characterize. The magnitude of the uniform flow is denoted by Qa and may be estimated from the
hydraulic gradient i and the aquifer transmissivity kH (hydraulic conductivity k times saturated aquifer
thickness H) Figure 5.2. The magnitude of the uniform flow rate is Qa — kHi. The flow Q0 is defined as
the total amount of water that flows in the aquifer per unit time per unit width. Techniques for estimating
the direction of the uniform flow field were discussed in Chapter 2.
The shape and size of a simplified time of travel capture zone may be related to a dimensionless travel
time parameter T/T, knowing the time of travel T for the capture zone, the uniform flow Qa, the saturated
aquifer thickness tl, the aquifer porosity n, and the pumping rate of the well Q, Figure 5.3.
T = ^	(5.1)
When T/T < 0.1, the calculated fixed radius (CFH) centered on the well, including a safety factor for a
non-zero ambient flow field, is
R = 1.154'Sy' (QT) / (nHn)	(5.2)
When 0.1 < T/T < 1, the CFR is
R - Ls[1.161 + /n(0.39 + j)\	(5.3)

Plan view
t H
Cross sectional area
q = Q/A
Qow = Q
Figure 5.2: The magnitude of the discharge vector.
where the eccentricity is
5 - Ls[0.00278 + 0.652-p]	(5.4)
given the distance from the well to the stagnation point downgradient from the well is
La = Q/(2ttQ0)	(5.5)
When T/T > 1, a uniform flow envelope can be defined as
x = y/tan(y/Ls)	(5.6)
and clipped at upgradient distance
Figure 5.3: Simplified delineation techniques for a well pumping at rate Q, in an ambient flow field Q0,
with aquifer saturated thickness II and porosity n, and time of travel capture zones of time T, after [Hait-

the corresponding capture zones were progressively more realistic. The emphasis of the W&AEM2000 project
on "ease-of-use" and computational efficiency does not release you, the modeler, from responsibilities in justi-
fying the conceptual models, and defending the reasonableness of the solutions. We emphasized uncertainties
in conceptualization of boundary conditions in our case study, but uncertainties in pararneterizations are
also important. WMEM2000 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
One more time, "'Ground water models should be as simple as possible, but not simpler." Good modeling!

Appendix A
Conversion Factors

To obtain





2.589988 x 106

2.2815 X nr4

6.9541 X nr5


Volume rate
2, 446.5754
rnA / day





Appendix B
W/1AEM2000 Binary Base Maps and the
B.l Getting BBMs 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 from within W/iAEM.2000, go to Tools select Base Map Browser, and select Use
Web. (Alternatively, from outside of the program you can point your web browser to
www.epa.gov/ceampubl/gwater/whaem/us.htm). 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 uncompress (by double
clicking the file) to your workspace. The binary basernaps (bbm) are now loaded into your project database
by using the Project and Project Settings, and Add BBM commands.
The EPA supplied BBMs are a binary form of the USGS DLG maps. The BBMs are saved in files
associated with their location and category, and the name reflects this information. The 1:100,000 scale
quad is divided in eight 15-minute quads. Table B.l. Each of these quads contains the DLG for the
associated category: hypsography, hydrography, roads, railroads, and miscellaneous. Table B.2. Therefore,
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:
B.2 Creating BBMs from USGS DLGs maps from the USGS Web Server
W/jAEM'2000 has utilities to convert the USGS digital line graph (DLG) maps into the binary base map
format (BBM").
The USGS DLGs are available for download on the Internet at their GeoData webpage. Point your web
browser to
edc.usgs.gov/geodata/. Figure B.2.
The DLG maps are available at two scales: 7.5 minute (or 1:24,000) and 30 x 60 minute (or 1:100,000).

¦¦it Ed* Mew Favorites Toob Met
>B«t • -J • J j] 3 ^Search _iJFavor*es QfMeda J -6* @ -

MAois |gj http: f/wwv. pea.  'i-wKfrrcter » feMi
United States Binary Base Maps
Binary base map (BBM) fifes are vectomed graphics derived from U.S. Geological Survey (USGS) diflal
line graphs BBM files can be imported directly into WhA£M2000 software and serve as Ifte geograph c
basis for giouicjrvater modeling To fid 88M data for a specific geographic area, click on the appropriate
stale below A mow dmailed graphic of th# siaie you solociod w>it display AMtmMivaly. you can choose
from an afflhabetical [e-i n of available states
ceam whaem.cdr
Figure B.l:
Table B.2:
Category Name
roads and trails
miscellaneous transportation
The 30 x 60 minute maps, are stored in optional format and SDTS format, while the 7.5 minute maps are
stored in SDTS formal. Complete coverage is available of the USA at the 30 x 60 minute scale, while coverage
at the 7.5 minute scale is not complete — you can check the status of maps for your state by going to the
web page
mcmc web .er.usgs. gov/status/dig _st at. html.
We will demonstrate the procedures for downloading the 7.5 minute (1:24K scale) DL.Gs for Vincennes,
Indiana from the USGS website, and convert to BBMs. The 30x60 minute BBMs are available from the EPA
web server as described above.
Erom the Eros Data center webpage, select 1:24K DLG. You can select the DLG maps by alphabetical
Listj by State, or by Graphics. Select. "FTP by Graphics", then select, conterminous 48 states, you will be.
presented with a. map of the USA. Locate your mouse, and click on the point on the map near west-central
Indiana that you want to zoom in to. Continue to point, and click until the quads are displayed as shown in
Figure B.3.
As an example, select the FrichtOn quad, and download the hypsography DLG (1669.820.HP.sdt.s.t.a.r.gz)
to your project directory.
From the W/tAEM2000 menu, select Tools , then Convert.. For the Erichton hypsography DLG, which
is l:24Iv SDTS format, select. SDTS to Basemap. (If you were working with the USGS l:100Iv optional

f$|ER05 Data Center, Sioux Falls, SD - Netscape
File Edit View Go Communicator Help
|| 4 3, $ st- [Si
Back Forward Reload Home Search Netscape
.J- >£
Print Security
m a

.Jf T B ookmarks Location: | http: //edc. usgs. gov/geodata/
Instant Message j§i Internet [_J Lookup L j New&Cool [{3
1 Support. Dell. Co

§41 ,k"i . •

,f; w/ui-:


| Products ^ j Science ^ J NASA DA AC
^ | Satellite Missions ^ | National Archive
1 About EROS
USGS Geographic Data Download
Data Sets: | 1:250KDEM 1.24KDEM 30 m NED 1:2M DLG 1:100KDLG 1:24KDLG LULC NLCD NHD
1:24,000 Scale Digital Line Graphs (DLG) SDTS Format Only

• Example

• Status Maps

• GLIS to check availability in Native or SDTS formats.

• Data User Guide (National Mapping Program FTP1

• Condensed User Guide (Global Land Information System HTML)

• Review the 00README before downloading!

• SDTS format only.

~ FTP via Alphabetical List

• FTP via State

• FTP via Graphics

• SDTS DLGs require Master Data Dictionarv

~ Back

( Water

U. S. Department of the Interior \ U. S. Geological Survey (USGS)
Please read this general Disclaimer
URL: http./fedc. usgs.gov//includes/footer_nomenu. html
Maintainer: EDC Web Master email at edcwebCcbusas. gov
Last Update:Monday, August 12, 2002
USGS Privacy Statement\Accessibility
Figure B.2:








(no data)
(no data)
(no data)
(no data)

(no data)
(no dat


(no data)
(no data)
(no data)
(no data)
(no data)
(no data)
(no dat

(no data)
(no data)
(no data)
(no data)
(no data)
(no data)
(no data)
(no dat

(no data)
(no data)
(no data)
(no data)
(no data)
(no data)
(no data)
(no dat

i 7

IL "
¦ IL


Figure B.3:

format DLG, you would select DLG to Basemap). Select the compressed file (*.gjs) you downloaded from
the Internet. Provide a name for the BBM. Be sure to limit yourself to 8 characters. W/iAEM.2000 will
automatically call a sequence of utilities (e.g., gzunzip, chop). The final utility will do the conversion from
DLG to BBM, and store the file in your workspace. You can now use this BBM in a new database, or add
it to an existing project with the Project , Base Map Settings, Add BBM sequence.

[Bakker and Strack, 1996] Bakker, M. and St rack, O. (1996). Capture zone delineation in two-dimensional
groundwater. Water Resources Research, 32(5):1309-1315.
[Blandford and Huyakorn, 1993] Blandford, T. and Huvakorn, P. (1993). Whpa: a semi-analytical model
for the delineation of wellhead protection areas, version 2.2. http://www.epa.gov/ada/models.html.
[Blum, 2000] Blum, J. L. (2000). ArcFlow User Manual (first draft). Minnesota Department of Health.
[Ceric, 2000] Ceric, A. (2000). Assessment of the applicability of simplified capture zone delineation tech-
niques for groundwater public water supply systems. Master's thesis, School of Public and Environmental
Affairs, Indiana Univcrsity-Bloomington.
[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.
[JDupuit, 1963] Dupuit, J. (1963). Etudes Theoriques et Pratiques sur le Mouvement des Manx dans les
Canaux Decouverts et a Travers les Terrains Permeables, 2nd ed. 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.
[Haitjema, 1995] Haitjema, H. (1995). Analytic Element Modeling of Groundwater Flow. Academic Press,
San Diego.
[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 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,
ht.tp:/ /www.epa.gov/ada/'download / project/wellhd.pdf.
[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,
http:/ / www.epa.gov/ada/download/models/whaern.pdf.

[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.
[Kelson, 1998] Kelson, V. (1998). Practical advances in ground-water modeling using analytic elements. PhD
thesis, Indiana University-Bloomington.
[Kelson ct ah, 1993] Kelson, V., Haitjerna, H., and Ivraemer, S. (1993). Gaep: a geographical preprocessor
for ground water modeling. Jlydrological 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] Kraerner, S. R. (1990). Modeling of Regional Groundwater Flow in Fractured Rock Aquifers.
PhD thesis, Indiana University-Bloomington.
[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, P., 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
[Podgorney and Robert W. Ritzi, 1997] Podgorney, R. K. and Robert W. Ritzi, J. (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. http://water.usgs.gov/pubs/twri/twri-3__B8/.
[Shedlock, 1980] Shedlock, R. (1980). Saline water at the base of the glacial-outwash aquifer near vincennes,
knox county, Indiana. Technical Report 80-65, U.S. Geological Survey.
[Strack, 1989] Strack, O. (1989). Groundwater Mechanics. Prentice Hall, Englewood Cliffs, NJ.
[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.
[USACOE, 2002] USACOE (2002). Groundwater modeling system. US Army Corps of Engineers
http:/ / chl.wes.arniy.mil/ software/gnis/.

[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, Cincin-
nati, OH.
[USEPA, 2003] USEPA (2003). Source water protection. Office of Ground Water and Drinking Water
ht.tp://www.epa.gov/safewater/protect.html. updated.
[USGS, 1997] USGS (1997). The universal transverse mercator (utm) grid. Fact Sheet 142-97, U.S. Geolog-
ical Survey.
[USGS, 2002] USGS (2002). Dlgv32 pro - windows display software for dig and drg data. US Geological
Survey http://memcweb.er.usgs.gov/drc/dlgv32pro/index.html.
[Vassolo et al., 1998] Vassolo, S., Kinzelbach, W., and Schafer, W. (1998). Determination of a well head
protection zone by stochatic inverse modeling. Journal of Hydrology, 206:268	280.