EPA/600/11-05/151
                                                         November 2005
Capture Zone Delineation for a City Wellfield in a
  Valley Fill Glacial  Outwash Aquifer Supporting
                    Wellhead Protection
                                By
                        Stephen R. Kraemer. Ph.D.
                    U.S. Environmental Protection Agency
                       Ecosystems Research Division
                           Athens, GA 30605
                        Henk M. Haitjema, Ph.D.
                        Haitjerna Consulting, Inc.
                         Bloomington, IN 47401
                         Victor A. Kelson. Ph.D.
                     Wittrnan Hydro Planning Associates
                         Bloomington, IN 47404
                     Office of Research and Development
                    U.S. Environmental Protection Agency
                         Washington, DC 20460

-------
Notice

The information in this document has been funded (wholly) or (in part) by the United States Environmental
Protection Agency under contracts (P.O. 8L-0425-NTSA, P.O. OL-2294-NASA) to Haitjcma Consulting. Inc.
and under contract (68W01032) to Computer Sciences Corporation. It has been subject to the  Agency's
peer and administrative review, and it has been approved for publication as an EPA document.  Mention of
trade names of commercial products does not constitute endorsement or recommendation for use.
   This document was typeset using WiriEdtb of Aleksandcr Simonic (www.wincdt.org) and MikTpX'.(including
LaTfjX2e) software of Christian Shenk (www.niiktex.org). Tf]X is a trademark of the American Mathematical
Society.

-------
                                                                                                Ill

Abstract

The purpose of this document is to introduce through a case study the use of the ground water geohydrology
computer program W/iAEM for Microsoft Windows (32-bit), or W/iAEM2000. W/iAEM2000 is a public do-
main, ground-water flow model designed to facilitate capture zone delineation and protection area mapping
in support of the state's arid tribe's Wellhead Protection Programs (WIIPP) 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 a well field tapping a glacial
outwash aquifer. WMEM2000 provides an interactive computer environment for design of protection areas
based on simple WHPAs (e.g., 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),
Digital Raster Graphic (DRG) or other electronic base maps. Geohydrologic modeling for steady pumping
wells, including the influence of hydrological boundaries, such as rivers, recharge, no-flow boundaries, and in-
homogeneity 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 import and export utilities for DXF files and Shapefiles. W/&AEM2000 has on-line help and tutorials.
Install scripts  and base maps are available  for download  from the EPA Center for Exposure Assessment
Modeling web  site (www.epa.gov/athens/software/whaem/index.html).

-------
IV

Foreword

The National Exposure Research Laboratory's Ecosystems Research Division (ERD) conducts process, mod-
eling,  arid 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 Office of
Research and Development (ORD)  strategic goals in the 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
uncertainty.
   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
predictive capabilities for complex environmental exposures.
   To support this mission, ERD is comprised of dedicated scientific and administrative professionals who
deliver high-quality, timely, effective and understandable products and applications that meet our customers
needs, with a commitment to exceed their expectations.
   Exposure models are distributed and supported via the  EPA  Center for Exposure Assessment Modeling
(CEAM).
   The  development and release  of the Wellhead Analytic  Element Model (W/iAKM2000)  is in support of
the State's  and Tribe's Source Water Assessment and Wellhead Protection Programs as required by the Safe
Drinking Water Act.  W/iAEM2000  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.
W/iAKM2000 assists in the delineation of the area contributing source  water to public water supply wells.
The development  of W/iAEM2000 has been supported by the EPA Office of Research and Development and
the EPA Office of Ground Water  and Drinking Water.
                                                  Eric J. Weber, Ph.D.
                                                  Acting Director
                                                  Ecosystems Research Division
                                                  Athens. Georgia

-------

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

2  PROTECTION ZONES I - DISTANCE CRITERIA                                     15
   2.1   Setbacks	L5
        2.1.1  Exercise: Glacial outwash wellfield	   15

3  PROTECTION ZONES II — RESIDENCE TIME CRITERIA                          19
   3.1   Calculated Fixed Radius  	   19
        3.1.1  Exercise: Vincennes,  Indiana	   20
   3.2   Well In Uniform Flow	   22
        3.2.1  Exercise: Well in Uniform Flow, Vincennes, Indiana	   23

4  PROTECTION ZONES III — FLOW BOUNDARIES CRITERIA                      31
   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: No flow boundaries	   40
   4.4   Aquifer  Tnhomogeneities	   44
        4.4.1  Exercise: Inhomogeneities	   44
   4.5   Riverbed Resistance	   45
        4.5.1  Exercise: Resistance Linesinks  	   48
   4.6   Discussion	   50

5  SIMPLE WHPAs                                                                         53
   5.1   Introduction	   53
   5.2   Back-of-the-Envelope  WHPAs	   54
   5.3   Exercise: Simple WHPAs arid no ambient flow   	   55
        5.3.1  Exercise: Simple WHPAs with ambient flow	   55

-------
vi                                                                                CONTENTS

   5.4  Exercise: Simple WHPAs with stream boundary	  58
   5.5  Exercise: Simple WHPAs with multiple wells	  58
   5.6  Choosing the Right Tool	  63

6  SUMMARY                                                                           65
   6.1  Good Modeling Practice	  65
       6.1.1   Rules of Thumb   	  65
   6.2  Conclusion	  67

A  Conversion Factors                                                                    69

B  W/iAEM2000  Base Maps and the Internet                                              71
   B.I  Getting Binary Base Maps from the EPA Web Server  	  71
   B.2  Creating BBMs from USGS DLGs maps from the USGS Web Server	  71
   B.3  USGS raster maps (DRG)	  74
   B.4  USEPA BASINS maps (Shapcfilc)  	  74

-------
CONTENTS                                                                                 vii

Acknowledgements

Thanks to all participants in the testing of the code, including the participants in the EPA Regional Train-
ing Courses: (R3 Philadelphia, June, 1999, August 2001, December 2003), (R4 Atlanta, March 2004),(R5
Chicago, May,  1999),  (R8 Denver, February 2001),  (R9 San Francisco, April 2002,  September 2005),  (RIO
Seattle, January, 2000), (HQ, March, 2000). Thanks to our initial external peer reviewers, Dr. Steve Silli-
rnan 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 Ms. Joan TIarrigan-Farrelly and Mr. Roy Simon. The development of WMEM benefited
from comments from the former EPA Ground Water Protection Technical Forum.

-------
Vlll
                                                                                    CONTENTS

-------
Chapter  1
BACKGROUND
The objective of this document is to Introduce and demonstrate the the use of the geohydrology computer
model W7tAEM2000 for delineating capture zones for pumping wells. The presentation will cover both basic
model operations and modeling practice.  We will use a capture zone delineation case study to illustrate a
step-wise and progressive modeling strategy.


1.1    Source Water Protection Programs

Ground water 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 arid 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 requires the United States Environmental
Protection  Agency (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,2005].
   The main steps of the source  water protection process involves the assessment of the area  contributing
water to the well or wcllficld, 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
safety.
   The USEPA Office of  Ground Water and Drinking Water established guidance on the criteria and meth-
ods for delineating protection areas [USEPA,1993], [USEPA,1994].  The criteria  include: (1) distance; (2)
drawdown; (3) residence time; (-1) flow boundaries; and (5) assimilative capacity. The methods  applying the
criteria range from the simple to the complex,  including setback mapping, calculated radius, and geohydro-
logic modeling (see Table 1.1). This case study document will review selected criteria  and their methods

                                                 1

-------
                                                                    CHAPTER 1.   BACKGROUND
             Table 1.1: Wellhead protection area delineation criteria and methods.
METHODS


CRITERIA

distance
drawdown
residence time
geohydrologic
boundaries
assimilative capability
arbitrary
fixed
radius


X





calc.
fixed
radius


X
X
X



well in
uni-
form
flow
field

X
X



geo-
hydro-
logic
mapping




X


geo-
hydrologic
modeling



X
X
X


transport/
transfor-
mation
modeling


X
X
X

X
using W/iAEM2000. 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  alternative modeling
approaches based on stochastic theory [Vassolo et al., 1998] ,  [Feyen et al., 2001].
    Modeling of flow and transport through geohydrologic systems is a complex undertaking.  Our presenta-
tion 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 abstrac-
tion 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 water
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 (Thiem's solution for a well with an image recharge well across
from the equipotential line).  Of course, the price paid is the lack of realism of this conceptual model when
compared to the complex reality of a partially penetrating well near a meandering stream  with a silty stream
bottom in a heterogeneous aquifer of varying thickness.
    A realistic conceptual model is often thought of as a representation of reality that includes most if not
all of the geologic and hydrologic complexities of an aquifer or aquifer system. The more  realistic (read
complex) a conceptual model is the more involved the modeling becomes, mostly for the  following reasons:

-------
1.2.  MODELING- PHILOSOPHY                                                                 3

   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 believed to have offered the following advice when  talking
about models of the universe:  "'Things should be  made as simple as possible ... but not simpler". In other
words, good modeling practice requires a balance between the realism of the conceptual  model and the
practical constraint of parameterizing the resulting ground-water flow model.
   When deciding on an appropriate conceptual model it is important to 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)chcmical 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.
   A goal of the ground-water flow modeling is to find a proper balance between a simple conceptual model
and a realistic one. The most economical way to judge whether or not a concept.ua! 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 the cndpoint of the
modeling. 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

-------
4                                                                   CHAPTER 1.   BACKGROUND

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-thumb'' 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 in the 1990s

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 were 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) was 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
programs:

   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 ground water 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  takes ambient ground-water flow into account (approximately) and
offers a significantly better approximation of the actual timc-of-travel capture zone.  That approach was
implemented in a computer code WHPA, made public  by the USEPA  [Bl an ford  and IIuyakorn.1993].  In
addition, WHPA allowed for wells in  flow fields  generated  by enclosed (rectangular) boundaries, such as
provided  by combinations  of canals (constant head) or geologic contacts (no  flow).  WHPA used particle
tracking to draw bundles of streamlines toward the well, including the boundary streamlines that define the
capture zone for the well. 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 was menu driven,  had a windows-like
graphical user  interface under DOS, and required minimal computer resources.

-------
1.4.  THE EVOLUTION OF          TO THE PRESENT                                        5

1.4   The  Evolution  of W/iAEM to the Present

W/iAEM (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 cal-
culates 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. W/iAEM currently lacks continuously varying aquifer bottom eleva-
tions, three-dimensional  elements, multi-aquifer flow, transient flow, etc.. which are supported in full-feature
ground water models. While this limits the applicability of VV/iAEM. 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, timc-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 accu-
rate) definition of the tirne-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 model user 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 ground water flow solver 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.1989], Haitjema,1995|.
The two-dimensional analytic element models invoke the Dupuit  assumption  jDupuit, 1863 , meaning that
resistance to vertical flow is negligible and that the head is constant with depth (zero vertical head gradients).
The Dupuit assumption has been shown to be very powerful for modeling flow in shallow aquifers.  The Dupuit
assumption is reasonable for capture zone delineation when the capture zone is an order of magnitude wider
than the saturated thickness of the aquifer. The significance of the Dupuit assumption is 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]. CZAEM is a DOS program with a command
line interface.  CZAEM  is available for download from the  EPA  Center for Subsurface Modeling Support
(http://www.cpa.gov/ada/csmos/models/czaem.html)
   A graphical user interface (GUI) and geographic analytic element  preprocessor, GAEP [Kelson et al.,
1993],  was added to solution engine CZAEM" in the first release of the modeling system W/iAEM  [TIaitjema
et.  al, 1995],  [Haitjema et.  al, 1994].  This version ran under Windows  3.1.  The development  of GAEP
opened the door for on-screen conceptual model building, where the model elements (e.g. wells,  line-sinks)

-------
6                                                                   CHAPTER  1.   BACKGROUND

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 Ground water Modeling System  USACOE,2005 .
    WMEM for Windows 95 was released in 1997.
    The release of W/&AEM2000 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,1998] that included arbitrarily shaped no-flow  boundaries. The release of WMEM2000
version 1 supported delineation by fixed radius,  calculated  radius, well  in uniform flow, and steady flow in
confmed/unconfmed aquifers with constant aquifer properties and recharge.
    W/&AEM2000  version  2 incorporated the GFLOW1 solution engine  [Haitjcma, 1995], and included new
elements: inhomogeneities and resistance linesinks.
    The current release. W/iAEM2000 version 3.  has additional support  for raster base maps, such as USGS
digital raster graphics (DRG) and shapefiles (SHP). Also, WMEM2000 has simple WHPA tools to include
the functionality of the EPA WHPA code. WM.EM2000 runs under Windows 98/NT/2000/XP. WMEM2000
is available for download from the EPA Center for Exposure Assessment
(http://www.epa.gov/athens/software/whaem/index.html).

1.5   The         Map

Ground-water models are applied to regular Gartesian 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
or computer screen. 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/iAEM2000can use any Cartesian coordinate system, including universal transverse mercator (UTM),
State Plane, local site coordinates, etc. The graphical user interface (GUI) for W7tAEM2000 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 electronic vector and raster maps.  The  Digital  Line
Graph series (DLG) 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 arc available for download
on the Web (See Appendix T3.2).  The Digital  Raster Graphics (DRG) are scanned bitmaps of the paper
topographic maps, and are also available at intermediate and large scales.  DRG maps are available for order
(See Appendix B.3).
    The USGS maps use the UTM projection [USGS,1997]. In this grid  system, the globe is divided into 60
north-south zones, each covering a strip 6 degrees wide in longitude. The  zones are numbered from 1 to 60,
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/&AEM2000  reads BBM files from the disk each time it redraws a map on the screen.

-------
1.6.  EXERCISE: EXPLORING  A  USGS  QUAD MAP
                  126°   120°   114°   108°   102°   96°   90°   84°   78°   72°   66°
                           Figure 1.1: UTM zones of the conterminous U.S.

W/iAEM2000 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.
   W/iAEM2000 has a Reproject Maps under the Tools command  that allows the user to reproject maps
from one UTM zone to another, or from UTM to State Plane.
   Maps in any projection stored in digital exchange format (DXF) may be translated to BBM format using
the Tools .
1.6   Exercise: Exploring  a USGS quad map

We will use a source water area delineation for a city well field in a glacial outwash valley of a large river.
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.
   The floodplain of the river is characterized by valley fill from glacial outwash, and consists of coarse sand
and fine gravel on top of underlying bedrock. The area 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.8 MGD in 2005.  You will learn more details
about the wellfield as you progress through the case study.
   In this exercise, we will place the wellfield in its  geographic context. You will locate the wellfield on the

-------
                                                                       CHAPTER  1.   BACKGROUND
                                  Office of Research and Development
                                  Office of Ground Water and Drinking Water
                                                   Version 3 2 0 - 29 April gOOB
                     .  Btnk Haitjema  <

MBisGahfttt aephwAfe^iJ*, Sins'. .,,;  -,, -v -- -
                                                             , W»rtH*«l*m« •;*'"/!  •
                                   fmtftti rt Republic domain, wth somai tftStt^J^hs, "I He
                             software nnd documcnlalion may be freely copied and dishitjuted
                             The graphical aMf-lMJlMb'C* «hd Sn^reMOFLOWt') «« Copyi njhted
                            by Haitjema Soft^j£»r^^sten'WHirfB<*^ii^^ttng« Int , and
                              are open source tpitiff* JSWfkjfr.few^JniWlJprAt- ^   rh "
                             utilities "gap/gun-zip" and "tar" are free software and available from
                              '' wwwgnu oig/dn^ctoiyfGNtf/' lies tihi'tt*COPYINGtxt>
                                Figure 1.2:  W/iAEM2000 splash screen.
digital U.S. Geological Survey 7.5 minute DRG topographic map (vincennes_24k.tif). The DRG 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-
tion.
  Step 1   Start W/iAEM2000  and familiarize yourself with the graphical user interface (GUI).
   From the Start button, go to Programs, start W/iAEM2000 for Windows. The "About W/iAEM2000
for Windows" splash screen will appear briefly, and credits are given. Figure 1.2. Select the "Start WhAEM"
option on the welcome screen. Figure 1.3.
   The workspace for W/iAEM2000  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 1.4.
  Step 2   Create a new project database.
   From 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 1.5.  Click Open.
   The "New Database Wizard" will prompt you for a project description and the units associated with
the input data.  Figure 1.6.  Enter "Glacial   outwash  wellhead   case  study a project description.
W/iAEM2000 requires consistent units, either meters and days or feet and days.  The "Distance Units in the
Basemap Files" will be dictated by your source files. Select meters for the "Distance Units in Basemaps
Files" by clicking on the 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.

-------
1.6.  EXERCISE:  EXPLORING  A  USGS  QUAD  MAP
                           Welcome  to  WhAEM
                          pVhAEMZQOO is a public domain ground water flow modeling system     I
                            iased on the Analytic Element Method. The program solves ste
                           Now in a single aguifer applying the Dupuit-Forchheirner assumption (2D
                           eqtns). WhAEM is specifically designed for the delineation of time-of-travel
                           capture zones for wellhead and source water assessment programs. The
                           software is provided "as is", without warranty of any kind, express or
                           implied, including but not limited to the warranties of merchantibility, fitness
                           for a particular purpose and noninfringement. In no event shall the authors
                           or copyright holders be liable for any claim, damages or other liability,
                           whether in an action of contract tort, or otherwise, arising from, out of or in
                           connection with the software or the use or other dealings in the software.
                                Tutorial
                               IJJIIIIJII

                             Start WhAEM
                                   Don't show this dialog box in the future
                                        Figure 1.3: Start W/iAEM2000.
                                   WhAEM far Windows
                                 Project  Edit  View  Ejement Model Joels  Window  Help
                                 D  e
                                     Figure 1.4:  W/iAEM2000  menu bar.
' ''3
;S
'/ A';,1
jajgjL 3rnpfe,whCT

Fife name: >/iKi-,H jj Open j
Fteoiiype: riuieii C'diabj-eFles I' whm) jj Cancel J
                                       Figure 1.5: Create new database.

-------
10
            CHAPTER 1.  BACKGROUND
                       New Database Wizard
                        -' Project Deicriptiprr
                          NewWhftEM Project
                             • Distance units in Basemap Files
                                    (• Metets
'.Feet
                              Units for Computations [Elevations, Properties, Pumping Rates)
                                     (*" Meters and Days       (• [Feet and Days!
                                  .Create, Database
  Abort
                                   Figure 1.6: New database wizard.
   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  1.7. 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 the file vincennes_24k.tif. Then click Open. The file will be displayed
in the box "Currently Used Basemaps". Enter in the Base Filename (DOS) box "vinbase".  Click OK and
wait a few seconds. The map will be loaded and plotted on the project workspace.  Figure 1.8.
   From the smart icons bar, test the zoom in, zoom out,  zoom to extents, zoom to window icons, scrolling,
and re-centering. Figure 1.9.

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

-------
1.6.  EXERCISE:  EXPLORING  A  USGS QUAD  MAP
                                                         11
                   Project Settings
                       Project Description
                       Glacial out wash wellhead case study
                       Currently Used Base Maps
                       eAwhaemyrAvineenne'wineennes 24k.tif
                                          Remove Map
Intermediate Model Files
Base Filename (DOS)
                                                    Drive
                        Path
                      J9t5e25k             pile:


                        Distance units in Basemap Files    Meters
                aWHAEM
                •JIN
                        Units for Computations:
Feet
                                                       OK
                                         Figure 1.7: Project settings.

-------
12
                  CHAPTER 1.   BACKGROUND
                                  •> ' .  0 »  '« 4 31  * * t
                                                                  .
                                                                 •-*•-
                                                         '       * ,'/-' *'• \f *

                                                                •'«$1&'J*
                                                           •«    4.^-^ $f <*• c,
                                                             ,-..^>-^Si\%
                                                          _^.B^,' isHtC....^ • ">. .»
                                                                                 Vinoennes_24k.cdr
                       Figure 1.8: The W/iAEM2000 base map view of USGS DRG.
                                             Zoom to
                                   Zoom out    extents
Scroll
right
Scroll
down
                               Zoom in     Zoom to      Scroll     Scroll    Re-center
                                         window      left       up
                                    Figure 1.9:  The View smart icons.

-------
1.6.  EXERCISE:  EXPLORING A  USGS QUAD  MAP
13
                  ""jumping center
              ^jr.   w wmmmi, •     ; ••
                            Figure 1.10: Zoom in view of city wellfleld.

-------
14                                                CHAPTER 1.  BACKGROUND

-------
Chapter  2



PROTECTION  ZONES   I  - DISTANCE

CRITERIA



2.1   Setbacks

One of the basic criteria for defining a protection zone surrounding a pumping well or wellfield is distance.
The method of designating a setback based on a fixed radius is commonly used by the States [Merkle et
al.,1996]. For example, at one time, Indiana had a setback distance of 200 feet between potential sources and
public drinking water wells.  The fixed radius is the first line of defense against surface contaminants that
could reach the wellhead and travel down an improperly sealed borehole into the ground water near the well
screen. A fixed radius wellhead protection zone is also used as a setback for potential sources of microbial
pathogens.  The assumption is  that pathogens  entering ground water outside the setback area will become
inactive before entering the well. The actual survival time of pathogens in ground water, however,  is still
the topic of ongoing research.  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:  Glacial outwash wellfield
 Step 1  In this step you will locate a well in the center of the city wellfield, and draw the setback protection
area.

   Open (or restore to the view) the project database vinbase, and make a duplicate database "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 recommended that you make a new project
database anytime you start a new exercise.
   Zoom to the city wellfield. Position the cursor near the pumping center and read off the UTM coordinates.
Note that 452650m, 4280665m is about the  pumping center.
   To create a well, from the Menu select Element ,  then select  New and click on Well. Move the cursor
to the approximate location of one of the five wells, and click the  left mouse button. A dialog box will open
as in Figure 2.1. Type in a name for a label, e.g.  ' 'well  1'.'
   Click on the  Other tab, and fill in a setback radius of 200 feet. See Figure 2.2.

                                              15

-------
16
CHAPTER  2.   PROTECTION ZONES  I  - DISTANCE  CRITERIA
                                 Well Properties
                                  Label
                                    General ] Other 1
                                    T Head-Specified
                                    (* Discharge-Specified

                                   X coordinate
                                   Y coordinate
                                   Discharge
                                   Radius
             (452646.1550  meters
             _______
                      0 rT3/daji
                      1  feet
                                               OK
                  Cancel
                                          Figure 2.1: Well dialog box.
                                  Well Properties
                                  Label       (welll
                                    General  Other |
                                    Setback radius        |
                                    f"" Trace particles from well
                    200) feet
                                    |~" Choose Pathline Color
                                               OK
                   Cancel
                                          Figure 2.2:  The Other tab.

-------
2.1.  SETBACKS
17
                                     Figure 2.3:  200 foot setback.
   Click OK to add the well to the database.  Zoom into the well to see the dotted setback circle.  Now
create the four other wells.
   Draw the wellhead protection area (Edit , Draw wellhead protection area). Name the protection
area "200 ft setback". With the left mouse button, click a series of short segments along the outside envelope
of the dashed setback circles. The closing of the wellhead protection area is accomplished with a right-mouse
button click, and the wellhead protection area will show as a blue line. Figure 2.3.
   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. And as in  this  case, when 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.
   Minimize the view to continue to the next chapter, or exit the program by clicking File and Exit.

-------
18                          CHAPTER 2,  PROTECTION ZONES  I - DISTANCE  CRITERIA

-------
Chapter 3
PROTECTION  ZONES  II  —
RESIDENCE TIME  CRITERIA
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 flow.

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 H(m), the cylindrical boundary of radius R(m) is delineated by an isochrone of residence
time t(days), which means that any water particle that  enters the cylinder or is present in the cylinder will
travel no longer than t days before being pumped up  by the well  (See Figure 3.1). The pumping rate of
the well is Q(m3/day), the areal recharge to the water table rate due to precipitation is  N(m/day), and the
aquifer porosity is n(-}. A water balance for the period t yields:

                                     N7tR2t + mrR2H = Qt                                 (3.1)
   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:
                                                                                        <3'2>
   When t becomes infinitely large, the radius R represents the complete capture zone:


                                                                                        (3.3)
                                                K 1 V

                                             19

-------
20
CHAPTER  3.   PROTECTION  ZONES  II — RESIDENCE TIME CRITERIA
                                 N
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 TVvrt becomes small, due to a small value of t or N or both, equation 3.2 reduces to:
                                                                                              (3.4)
   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 H, 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  H in Equation  3.4.
3.1.1   Exercise: Vincennes, Indiana
 Step 1  Use W/iAEM2000  to plot the Calculated  Fixed  Radius protection areas  based  on the  recharge
method and the volumetric method.

-------
3.1.  CALCULATED  FIXED  RADIUS
21
               Project Edit View Element Model lools Window Help
                                    Recharge CFR (t=infinite)
                                           ~
                                    Volumetric CFR (t=5
   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.
   We cannot ignore the  effect of interference of the multiple pumping wells.  So for practical purposes we
will define a single pumping center, with the total pumping rate. We will assume the pumping center of the
wellfield is located at x=452650 m, y=4280665 m.
   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.00276 ft/day(l2.lin/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 IT = 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.
   Save these protection  areas for future overlay.  This is a two  step process.  First, export to shapefile:
Tools , Export, Shapefiles, Wellhead  protection areas,  and name the files.  Then, convert to BBM:
Tools , Convert Shapefile to BBM, and name the file vin_cfr5.bbm.
   Compare your results to Figure 3.2.  (3.1.1).

-------
22                        CHAPTER, 3.  PROTECTION ZONES  II  — RESIDENCE TIME CRITERIA

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
flow.
   To determine this capture zone (e.g. by use of W/iAEM2000) it is necessary to  determine:
   • direction of the ambient flow.

   • the hydraulic gradient.

   • the aquifer transmissivity.

   • the pumping rate of the well.

   • the desired maximum residence time.
   Given three water levels (average, static), it is possible to estimate the direction of the ambient flow and
the hydraulic gradient [Heath, 1983](See Figure 3.3):

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

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

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

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

   e. Divide the difference between the head of the well 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
http://www.epa.gov/athens/learn2model/part-two/onsite/gradient3ns.htm.
   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 unconflned 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 transmissivity is measured by  a pumping test. This  is  an  expensive procedure, and trans-
missivity is known to be scale dependent. Perhaps transmissivity measurements arc available from locations
not too far from the  well.  Estimates of hydraulic conductivity based on geology (e.g., rock type, unconsoli-
dated media type, grain size)  are available  in text books.  Check the USGS  Ground Water Atlas for regional
information about your aquifer system http://capp.water.usgs.gov/gwa/gwa.html.

-------
3.2.  WELL  IN  UNIFORM FLOW
23
         Will 2
     Ibtsi, SS.?0 HI
                                    WiEl I
      0   2$  50     100 METERS
       i   i    i  	i
                                      £6.Q7n)
                                       Heath.cdr
                                                                       (ZG.26-2S.Or)
                                                                                     J6 JB 91
Figure 3.3: Procedure A for determination of an equipotential contour and direction of ground water flow in
homogeneous,isotropic aquifer (after [Heath, 1983]).
   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.

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 vincfrl.  From
the Menu, Project , Make Duplicate Database, and enter vinunil.
   As a first estimate, we will  attempt to characterize the uniform flow field  by using the information
contained in the USGS DRG topographic map. The "lakes and ponds" in this area are actually abandoned
quarries, and the removed overburden gives a window into the elevation of the water table. We assume these
pits act like large piezometers, and field observations confirm that their water levels do in fact respond to
recharge 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.

-------
24
CHAPTER 3.   PROTECTION ZONES  II  — RESIDENCE  TIME  CRITERIA
                Project Edit View E.temenl Model Joels Window Help
                O|B|H| a| x|.<.|  |g-| -H»N '.
                                             '•/ test point - 400 ft topo
                               test point - Otter pond (403 ft)
                                                     flow function - Mirror Lake (407 ft)
                Select    Click mouse button to seteet element
                                                                         459,430.7m  4.281 
-------
3.2.   WELL   IN  UNIFORM FLOW
                                                                          25
                                    Model Sellings
                                     ...Aquifer j| Contouring] Tracing! Solver]
                                      Aquifer Properties
                                        Base Elevation
                                        Thickness
                                        Hydraulic Conductivity
                                        Porosity
J330
|100
|350
[02
                                        feet
                                        feet
                                        feet/day
                                        [dimensionless]
                                                             OK
                                                                           Cancel
                                            Figure 3.5: Aquifer settings dialog  box.
                                    Mode! Settings
Aquifer  Contouring | Tracing ] Solver
 Piezometric Contour Settings
  p1 Show Contours  <*  Coarse  C Detailed
  Minimum Contour
  Maximum Contour
  Contour Interval
J39Q
1430
                                                                             feet
                                                                             feet
                                                                             feet
                                                             OK    |        Cancel   j
                                          Figure 3.6:  Contouring settings dialog box.

-------
26                         CHAPTER  3.   PROTECTION ZONES  II  — RESIDENCE TIME  CRITERIA

                                            - ,uo Model
                                            Solve  Settings
                                                    n-
                                       Figure 3.7:  Solve icons.

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.

         Solve the model and plot the contour map.

   From the menu select Model , then Solve (Alternatively, use the F9 mapped function key, or the smart
icon Figure 3.7).  A black  box appears with some progress indicators for the solution process. When the
solver is  done the base map reappears with an overlay of piezometric contours (blue lines) that are labeled
with their  piezometric head elevations.  On top of the base map the two log files, specified on the Solver
tab, are opened in Notepad.  The file Error.log should have only an asterix on the start of the first line and
be empty otherwise.  The file is used by the Solver to report problems  (errors or warnings) with the  input
data.  Close the file after reviewing its contents.  The file Message.log contains Solver activity messages and
a report  of maximum errors  in boundary conditions. In general, errors should be less than 1%.  Close the
file after  reviewing its contents.
   The hydraulic head contour map of this solution is shown in Figure  3.8.
   Note:  since our model heads will be close to the test point heads,  the triangles designating the head difference will
be too small to see. Change the direction or magnitude of the uniform flow field and you will see the triangles grow in
size. You  may modify the test point appearance by selecting Test Points on the View  menu.
 Step 5  Enter the pumping rate and radius of the pumping center.
   Locate the Vincennes pumping center either 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 4/t
to represent our pumping center.  The Vincennes wellfield 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 water 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.

-------
3.2.  WELL IN UNIFORM  FLOW
                            27
                   Project Edit View Element Model JTools Window Help
                        l  a|   |  |"|  I -:-l;;N -.|-
Vinuni1.cdr
                   Select      Click mouse butlon to sdect element
                   Figure 3.8: Hydraulic head  contours of uniform flow field with test points.
                                     Well Piopert.es       (F1 foi Help)
                                     Label        jwell field

                                       General ] Other ]
                                       C Head-Specified
(•" Pischafge-Specified
X coordinate 452650 meters
Y coordinate | 4280665 meters
Discharge
Radius

| 370000 r3/day
| 4 feet

OK | Cancel


                                             Figure  3.9:  Well dialog box.

-------
28
CHAPTER  3.   PROTECTION  ZONES II —  RESIDENCE TIME CRITERIA
               Project Edit View Element Model Joels Window Help
               DlBl'al'al  I M  I-:-lsUI • l-?kk-|
                                                              Uniform flow 5 yr
               s*d
                      Click rrouse button lo setecl etemeni
             Figure 3.10: Capture zone based on 5 years residence time in uniform flow field.
   From 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 335 feet
(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.
   Let's save the protection area for future overlay.  First, export to shapefile: Tools , Export, Shapefiles,
Wellhead  protection areas, and name the files. Then, convert to BBM: Tools , Convert Shapefile to
BBM, and name the file vin_ uniS.bbm
   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.

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

-------
3.2.   WELL  IN  UNIFORM  FLOW
29
                 Project £dit View Element Model look Window Help
                 g|B|a|a|  |  |"|  I •• l»l-l  - |-?KU I
                                                                               452.342.0m  4.282.191.2m
Figure 3.11: Capture zone based on 5 years residence time based on CFR volumetric and tracelines for well
in uniform flow field (zero gradient).

-------
30                    CHAPTER, 3.  PROTECTION ZONES II — RESIDENCE  TIME  CRITERIA

-------
                 4

PROTECTION  ZONES  III  —  FLOW
BOUNDARIES  CRITERIA
In the previous analysis we 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 W/&AEM2000. We will progressively refine the complexity of our  hydrogeological
conceptual model using the glacial  outwash example.  But before doing so we will  explain  briefly how the
analytic element methods work. The reader is encouraged to explore in more depth  the topic of system and
boundary conceptualization outside of this  tutorial [Rcilly, 2001].


4.1   Introduction to the Analytic Element Method

W7iAEM2000 uses the analytical element method to solve ground-water flow subject to hydrological boundary
conditions. Analytic elements arc mathematical functions that represent hydrogcologic 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 ground-water 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  W/iAEM2000 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/&AEM2000 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.f989]  and
[Haitjema,1995].
   Other examples of analytic  elements in W/iAEM2000 are  wells  (point sinks in two-dimensions),  line
doublets or "double layers" to represent aquifer irihornogeneities, 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 WMEM2000 	the head and average ground-water flow
rate can be calculated at any point by adding the influence of all the line-sinks, point-sinks, line doublets.
and areal elements. This is the principle of superposition.

                                               31

-------
32                     CHAPTER  4.   PROTECTION  ZONES III —  FLOW BOUNDARIES CRITERIA

   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 [Haitjema,1995].
   Analytic elements are defined in W/iAEM2000 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/iAEM2000  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 W/iAEM2000  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
water along its boundary only.  This appears to be a good  approximation  of reality when the size of the
lake  is large compared to the saturated thickness of the aquifer. For 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
Line-sinks.

4.2.1   Exercise:  Vincennes and  Line-sinks
         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.  Make  sure the raster  base
map vincennes_24k.tif and vincennes_100k.tif  are active  in the window. Find the location where the
400/t topographic contour crosses the Wabash River north of the city. Note below the UTM location of the
400/t contour:
                                utm — x(easting} =	                            (4.1)

-------
4.2.  RIVERS AND  LINE-SINKS                                                                33

                                                    Add
                                            Add      no-flow
                                            well      string
                                               Add
                                               line-sink
                                               string

                                  Figure 4.1: Elements smart icons.

                               •utm — y (northing) =	                           (4-2)

Zoom In to the location of the 400/t 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.
   It is often sufficient to add hydrography labels near the map boundaries and near confluences of streams.
When creating line-sink strings you will only be prompted  for a head at the start of the string and at  the end
of the string.  It is good practice to start a linesink string  at the higher elevation and work down. This will
ensure more meaningful data reporting in the graphical user interface (GUI). The W/iAEM2000 program
linearly interpolates  between the  starting head and ending head for the heads at  vertices of the linesink
stream.
   We have created a text file 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 fa = 1 for hydrography
labels (water levels in surface waters) and fa = 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
Hydro.  Example.txIn order to view the hydrography labels on the basemap, you may need to toggle them
on (View > show hydrography labels).

         Represent the River with a string of Line-sinks.

   At this point, it is recommended that you set your window to frame the region  including the 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.
   Note: It is recommended to enter line-sink strings along streams going down-stream (into  the direction of the

-------
34
CHAPTER 4.  PROTECTION  ZONES III — FLOW BOUNDARIES CRITERIA
                              tLine3ink String Pmperliet   (F1
                               Label      |Wabash River
                                 General |
                                £•  Head-Specified
                                f  Discharge-Specified
                                P7  iTreai"as"ivfar-fiidvi
                                Starting Head
                                Ending Head
                                  400 feet
                                  395 feet
                                            OK
                                 Cancel
                               Figure 4.2: Linesink properties dialog box.
stream flow). This is particularly important when querying line-sinks that are part of stream networks, as it will assure
that the stream flow reported actually occurs at the vertex that is highlighted.
    That completes your first line-sink string.  At the center of each line-sink in the string a head will be
specified by W/iAEM2000  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 hydrologic (line-sink) boundaries and solve.
   Position the pumping center as described in the previous Chapter, and create the discharge specified well
(Qw  = 370,000 ft3/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 wellfield area, those include City Ditch, Mantle Ditch, Kelso Creek, Ditch 2, Ditch 3, and the Wabash
River, of course [Shedlock, 1980].  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 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).

-------
4.2.  RIVERS AND LINE-SINKS
35
                Figure 4.3: Perennial streams in the study area, after [Shedlock, 1980].

-------
36
CHAPTER 4.   PROTECTION ZONES  III  — FLOW BOUNDARIES CRITERIA
                              i*h«»*9*ntiw PMpertte*      ffl tat;
                              Label      Jareal recharge 12.1 in per jir
                                General |
                               Added Recharge Rate
                            0.00276  ft/day
                                                 Color...| |
                               |~ Change base elevation from default
                               |~ Provide estimate of average head
                                        OK
                             Cancel
                          Figure 4.4:  Inhomogeneity element properties box.
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 wellfleld. 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.  We 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.
 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!

-------
4.2.  RIVERS AND LINE-SINKS
37

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

-------
                       CHAPTER 4.   PROTECTION ZONES  III  — FLOW BOUNDARIES  CRITERIA
               Project Edit View Element Model Jools Window Help
                      l  | |"|  I -:-|SUM-?kkl DlB^lQlBlDl
                                                                         455,188.0m  4,281,166.0m
            Figure 4.6: 5 year capture zones including the influence of hydrologic boundaries.
 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 study area for the period January
23-25, 1978  [Shedlock,  1980]. These water levels have been stored in a Test Point file. To import this file,
select Tools and then Import the file vin-mw.tp. 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.
   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

-------
4.2.  RIVERS  AND LINE-SINKS
                      39
                Project Edit View Element Model Joels Window Help
                DlBl'al'al  I  M I-:-lsUI • l-?kk-|
Twtpolnts1.cdr
                                                                 A/ / /  ,/ /  /  ,
               Setecl     Click mouse button to setecl etemeni
                                                                          454,144.4m  4,280 .^J ' rn
           Figure 4.7: Model head contours compared to field observations at monitoring wells.
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 W/iAEM2000  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 W/iAEM2000 model particularly if they show up as ephemeral
or intermittent steams on USGS topo maps.

-------
40                    CHAPTER 4.  PROTECTION ZONES III — FLOW BOUNDARIES  CRITERIA

4.3   Geologic  Contacts  and No-Flow Elements

The geology of the 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 and 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 wellfleld.
   Ideally,  we should represent  the high permeable rock and till  explicitly in our model.  As a first step,
we will represent in W/iAEM2000  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.
   W/iAEM2000 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.


4.3.1   Exercise: No flow boundaries


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

   First, create a duplicate database as nof lowl. 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 wellfleld.  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 wellfield and are irrelevant to the local capture zone solution. You should avoid crossing line-sink
elements and no-flow elements to avoid numerical difficulties.

         Solve and check the calibration.


         Plot the 5 year capture zone.

   Double click on the wellfield, 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 wellfield (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,
as we  shall see in the next  section.

-------
4.3.   GEOLOGIC  CONTACTS  AND  NO-FLOW ELEMENTS
                                                       41
                       Surf-geo.cdr

                           B
                                ft;
                                                   Geology (torn H.H. Gray,  1.1.  Wayne and C E. Kier (1970)


                                               EXPUN:AT ION
                              Alluvial silt.
                              sand,  antl fjra ve I
                              Wi lul-t] I own sand
                              anil so~e s-ll
Kind-liloKn sill and
< i tie sand and  clay
Ontwash gravel,
sand,  and silt

                      Figure 4.8:  Surficial geology of the study region.  [Shedlock,1980]

-------
42
CHAPTER 4.  PROTECTION ZONES  III  —  FLOW BOUNDARIES  CRITERIA
                           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

                   IERTICH SCM.E BSEtTLT
                                          S«olO(l cod: Hill 1ro» N.H. tci,,
                                          (.1. »irm, •«!« C.E. I In (H7C)
Figure  4.9:  Generalized geologic cross section of the Wabash River Valley near the study region [Shed-
lock, 1980]
                 Eroject £dit View EJement Model Idols Window
                      l a|
                                                                                   Llne-dlpole8.cdr
                 Select     Text: 460.0
                                                                               453,635.5m  4.283,334.3m
Figure  4.10:  Superposition of no-flow elements along  the geologic contact  inferred from the topographic

contours.

-------
4.3.   GEOLOGIC  CONTACTS  AND  NO-FLOW  ELEMENTS
                         43
                  Project Edit View Element  Model Jools Window Help
                       l a]  '|             '
                                                                                   448,694.8m  4,282,648.6m
            Figure 4.11: Hydraulic head contours and test points, including the no-flow elements.
                  Project Idil View Element Model Jools Window He^)
                   BlHlal  I  I  I  I -HHl^l • |t?k-UI nlfcl ltl*|-
Noflow5yr.cdr
                  Setect      Clicl mouse bi;i;i-:n o seled elerr.erpi
                                                                                   454,025.4m  4,273,536.0 rn
  Figure 4.12: 5 year capture zone solution including the influence of hydrologic  and no-flow boundaries.

-------
44
CHAPTER 4.  PROTECTION ZONES III
FLOW BOUNDARIES  CRITERIA
                      Drawing the inhomogeneity
                                  Enter vertices.
                   Right-click on ^^
                   last vertex closes
                   the domain.
               Do not create sharp
               corners!         \
                                        Incorrect: line-sink
                                        should not intersect
                                        line-doublet.
                                                Correct: line-sink vertex is
                                                on line-doublet.
                       Figure 4.13: Drawing inhomogeneities using the mouse.

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

         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
3.9).
   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
k = 350 ft/d. Label the inhomogeneity.  Check that the recharge rate is  correct. Click OK. Use the hypsog-
raphy basemaps (vilhpfOG, vilhpfOT) to guide your drawing the polygon, following the sharp transition from
outwash to  rock. Figure 4.15.

-------
4.5.  RIVERBED RESISTANCE
                                                45
                               Inhomogeneily Piopertie::
                               Label      joutwash
                                 General
                                Conductivity
                                Added Recharge Rate
                                Porosity
   350 ft/day
000276 ft/day
   0.2
                                                  Color..
                                f" Change base elevation from default
                                         OK
                                                     Cancel
                                                           lnhomo_properties.eps
                   Figure 4.14:  Dialog window for entering inhomogeneity properties.

          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.txThe 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
Arc View.  The near field linesinks representing the  center channel of  the Wabash River between elevations
400 ft and 395 ft should be entered as in a previous exercise.  See Figure 4.15.

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

-------
46
CHAPTER  4.   PROTECTION ZONES  III —  FLOW BOUNDARIES  CRITERIA
                                                                I              /
                                                                    Layout inhomo.eps
                                                      linesink string (ferfield)^
                                                 IrTK   J  ,,.-**«^i,  —eft  \
             Figure 4.15: Layout of linesinks, wellfleld, and aquifer inhomogeneity domain.
       Figure 4.16:  5 year capture zone solution including the influence of aquifer inhomogeneities.

-------
4.5.  RIVERBED  RESISTANCE
                                                    47
                                                  w
                                                        ?z=-
                                                               c
                                                           = wqz
                       Figure 4.17:  Conceptual model of stream/lake resistance.


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.
                                                   kc
                                                                                                (4.4)
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.

                                            qz = ^                                             (4-5)
The discharge of the stream segment a[L3/T] is equal to the specific discharge times the width.
a = wqz
                                                                                                (4.6)
   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.
   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 actual 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:

                                                                                                (4.7)
   There are some rules of thumb on setting the effective line-sink width parameter, given the stream width
and the leakage length.

-------
48                     CHAPTER  4.   PROTECTION  ZONES III — FLOW BOUNDARIES  CRITERIA
    Put line-sinks along both sides
    of the stream of width w:
    w = X                              A < O.IB
    w = Atanh(JB/2A)                O.IB < A < IB
    w = B/2                             \>2B

    Put line-sinks along the center
    of the stream of width w.
    w = B                               \>2B
4.5.1   Exercise: Resistance Linesinks
 Step 1  Calculate the resistance of the bottom materials in the Wabash River.
   The USGS report [Shedlock, 1980]  includes a discussion of the leakage from the Wabash River to the
outwash aquifer in the vicinity of the city pumping center.  They report a leakage layer of average hydraulic
conductivity 0.03 ft/d, and of local values of 0.07 and 0.7 ft/d, and of average thickness of 1 foot.  Let's work
with the value of kc = 0.07ft/d and 5 = Ift, and calculate the resistance using equation 4.4.
 Step 2  Calculate the characteristic leakage length for the Wabash River.
   Assume an aquifer thickness of H = 70ft and aquifer hydraulic conductivity of 350/t/d. Use equation
4.7 to calculate the characteristic length.
         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 raster base map
vincennes_24k.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 resistl. 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.

-------
4.5.  RIVERBED  RESISTANCE
49
Lrnesink String Properties (F1 for Help)

Label IKBSEffl
General |
f* Head-Specified
f" Discharge-Specified
P* Treat as "far-field"
Starting Head
Ending Head
Resistance
Width
Depth

OK

E


400 feet
337.5 feet
14.3 days
282 feet
1 feet
Cancel


                        Figure 4.18:  Dialog box for entry of resistance linesinks.

-------
50
CHAPTER  4.  PROTECTION ZONES  III
FLOW  BOUNDARIES CRITERIA
Figure 4.19: 5-year capture zones including resistance of Wabash River (c = 14.3 days, 5 = I ft, A = 592 ft,
B = 615 ft, and w = 282 ft
         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
sequence.
   Haitjema and  Kelson investigated the influence of complexity of the conceptual model on  the capture
zones for our study area [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, anisotropy,
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

-------
4.6.  DISCUSSION                                                                                51

the many capture zones generated by Haitjema 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 city wellfield. 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 city wellfield was for illustration purposes, and is not to be misunder-
stood as the approved State source water assessment.
   Now some other complexities.  We have not  introduced the concept of horizontal anisotropy in the
hydraulic conductivity distribution.  Horizontal anisotropy 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, [Kraerner, 1990],[Podgorney, Ritzi. 1997].  W/iAEM2000 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 concentration 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/athens/learn2model/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
(http://www.grac.org/pest.html).
   There are  other professional analytic element modeling packages,  including TWODAN,  GFLOW, and
MLAEM. The reader is encouraged to follow' developments of ArcAEM which is based on the solution engine
SPLIT by Igor Jankovic, State University of New York-Buffalo, and an ArcGTS interface. Three-dimensional
flow  analytic element 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/whaern/press.htiTil).

-------
52
CHAPTER 4.   PROTECTION ZONES  III —  FLOW BOUNDARIES  CRITERIA

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

-------
Chapter  5
SIMPLE WHPAs
5.1   Introduction

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

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

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

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

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

                                               53

-------
54
       CHAPTER  5.   SIMPLE  WHPAS
                                                Q
                dx
                                                        Q
                                                       wH
                   Plan view
                    w
                              tH
            Cross sectional area
= -kiH
                         Figure 5.1:  The magnitude of the discharge vector.

5.2   Back-of-the-Envelope  WHPAs

How does the Simple WHPAs tool work?  The calculations depend on several parameters, including the
magnitude and direction of the ambient flow near the well or well field, which is challenging to characterize.
The magnitude of the uniform flow is  denoted by Q0 [L2/T] and may be estimated  from the hydraulic
gradient i [-} and the aquifer transmissivity kH (hydraulic conductivity k times saturated aquifer thickness
H)  [L2/T]. See Figure 5.1. The magnitude of the uniform flow rate is calculated from
                                           Q0 = kHi
                                 (5.1)
The flow Qo is the total amount of water in the aquifer integrated  over the saturated thickness, per unit
width of the aquifer.
   The ambient  flow is challenging to estimate, and the reader is  referred  to the previous discussion of
Section 3.2.
   The shape and size of a simplified  time of travel capture zone may be related to a dimensionless travel
time parameter T. The term T is a dimensionless time-of-travel parameter defined as
                                                 T
                                                 -LO

where T is the time-of-travel [T] and T0 [T] is a reference time defined as:

                                                nHQ
                                          T =
                                            °
                                               2vrQ02
                                                                                            (5.2)
                                 (5.3)
where n is the aquifer porosity [-], and Q [L3/T] is the pumping rate of the well.

-------
5.3.  EXERCISE: SIMPLE WHPAS AND NO  AMBIENT FLOW                               55

   When T < 0.1, the calculated fixed radius  (CFR) centered on the well, including a safety factor for a
non-zero ambient flow  field, is
                                                                                            (5-4)
   When 0.1 < f < 1, the CFR is
                                   E = Ls[1.161-Mn(0.39 + T)]                                (5.5)
where the eccentricity is
                                     5  = Ls [0.00278 + 0.652T ]                                  (5.6)
given the distance from the well to the stagnation point down gradient from the well is

                                         Ls = Q/(2vrQ0)                                      (5.7)

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

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

given
                                     -QI1Q0 
-------
56
                                                                          CHAPTER  5.   SIMPLE WHPAS
                                 < 0.1
             Centric circular capture
                      zone
                                           Back-of-the-Envelope
                                                Technique
                                                     'i>O.J &< l
                                         Ecentric circular capture
                                                   zone
Boat-shaped
capture zone
Figure  5.2:  Simplified delineation techniques  for a well pumping at rate Q, in an ambient  flow field Q0,
with aquifer saturated thickness H and porosity n, and time of travel capture zones of time T, after [Hait-
iema.20021.
jema,2002].
                                Calculated WHMf
                                 Selected Well


                                 Time of Travel

                                 Aquifer Porosity

                                 Aquifer Thickness

                                 Ambient Flow
                                          Continue
                                 Explanation
                                                        well field
                                                              1826.25  days
                                                                 0.2
                                                                67.E) ft
                                                                  0 ff 2/day
                                                          Cancel
                                 This can only be done for one well at a time. Based on the
                                 data, a dimensionless travel time parameter T is calculated,
                                 used to decide on a centric circle, excentric circle, or pathlines
                                 for well in ambient flow.
                               Figure 5.3: Dialog box for calculated WHPA.

-------
5.3.  EXERCISE: SIMPLE WHPAS AND NO AMBIENT FLOW
                                     57
                                     recharge
                                            =  \.
                           ••'••' O'Nt^i. .'•  f-.M''.^ tlrf-ir. Ihtftar

                          •'-' .SiRi-THf ' '••  il'is"' "-:.?'-;.
                                          ?,,
                                CFR volumetric 5 yr "='%>-
                                                *.!»
   *   ».,i^-v   "• •' '' S r'\ :V*Vj, Jif \"';!H..
       •Ji-!;*x '» ,-^'' x ^.V"<--~---V,,i« S *: BT'^Vji ,.->«
        V.'N  ^^o^^H i*V^

   ,«= K,,^^^^Hf§|^₯i^S
     ..**"V <^V:S^^V^«^&i,
A-«-l«|F.'sW^:l'
f   ''i**jwi«Vv^.«\*-V •••*.  ^^-^"-^iw-S^?'1 -,
   .*-?"*  ""' '•.-.--••' *-;'""* ^- * '^- --=• ' :"J';- ?.?' v- •'^X.1 ". ",'i    /
                    x&,  E
                                                                      JL
Figure 5.4: 5-yr Simple WHPA with no ambient flow compared to calculated fixed radius solutions using

volumetric equation and recharge equation.

-------
58
           CHAPTER 5.  SIMPLE  WHPAS
                              Selected Well
                             1 Time of Travel
                             I
                             I Aquifer Porosity
                             i
                             j Aquifer Thickness
                             1
                             i Ambient Flow
                            ; Explanation
                                                   well field
 1825  days

  0.2

 67.5  ft

23.625  ff2/day
| Continue
.....
                                                      Cancel
                            | This can only be done for one well at a time. Based on the
                            j data, a dimensionless travel time parameter T is calculated,
                            I used to decide on a centric circle, excenttic circle, or pathlines
                            1 for well in ambient flow.
                   Figure 5.5: Dialog box calculated WHPA including ambient flow.
 Step 2  Add the wellfield.  Zoom in and add the single pumping center well as before (utm-x=452650,
utm-y=4280665, discharge=370000 ft3/d,  radius 4 ft).  Add other properties, such as 20 particles starting
 Step 3  Calculate the Simple WHPA with ambient flow. Fill in the dialog box shown in Figure 5.6.  Fill
in the next dialog box with a direction of flow of 150 degrees, as shown in Figure 5.7 Note the T value is
greater than 1.
   Click the no response to the issue of a stream boundary, see Figure 5.8  and Draw.
   The slight difference in the solutions for the Simple WHPA and the well in uniform flow function is that
the former is solved under the assumption of constant transmissivity, while the latter allows for unconfined
flow conditions. See Figure 5.9.
5.4   Exercise:  Simple WHPAs with stream boundary

The Simple WHPA tool can be used to create capture zones for the well field including the influence of the
stream boundary.  The reader should follow the instructions  in the dialog  boxes to create something like
Figure 5.10. Notice that  the 5-yr capture zone has shrunk in size since the well is now getting water from
the river.
5.5   Exercise:  Simple WHPAs with multiple wells

The Simple WHPA tool can be used to create capture zones for multiple wells within an ambient flow field.
The reader should follow the instructions in the dialog boxes  to create something like Figure 5.11 Notice
that the trace lines from each well can be assigned a different color.

-------
5.5.  EXERCISE:  SIMPLE   WHPAS  WITH  MULTIPLE  WELLS
                                                                                59

                                                   7 = 1.2813001557899 (TM)

                                       Direct!on of Flow               150  Degrees


                                                  Continue   '


                                       Explanation
                                       The direction of flow is measured in countei clockwise
                                       direction from the horizontal x-axis. For example, north is
                                       90' and south is -90*.
                                       Make sure that the "trace particles" on the "other" tab of
                                       the Well Properties dialog is set and a sufficient "number
                                       of particles" are entered.
                                  Figure 5.6:  Dialog  box simple well in uniform flow.

                                        Do you want to use a Stream Boundary
                                                  Yes
                                                                •* No)
                                                    Draw
                        Back
                                       Explanation
                                       Sir
Boundary:
                                       Explanation: You will define the (infinitely long) straight stream
                                       boundary by entering 2 points on the map with the mouse cursor
                                       (left click). The specified head will be applied to the first point you
                                       enter. The head along the stream boundary may vary depending
                                       on the orientation and strength of the uniform flow and the
                                       hydraulic conductivity you entered. Both pathlines and head
                                       contours will be drawn. Make sure that the "Show Contours" is
                                       selected on the "Model>Settings> Contouring" tab and the proper
                                       contouring parameters are selected.
                                       No Boundary:
                                       Only pathlines will be drawn to define the WHPA(s).
                                    Figure 5.7:  Dialog box simple stream boundary.

-------
60
CHAPTER  5.  SIMPLE WHPAS
          Figure 5.8: Simple WHPA including ambient flow compared to well in uniform flow.

-------
5.5.  EXERCISE:  SIMPLE  WHPAS  WITH MULTIPLE WELLS
61
                      f   <      /'
                           . ,?,/••
                            m
                          /

Figure 5.9: Simple WHPA including the influence of a stream boundary in comparison to capture zone for

a well in uniform flow field.

-------
62
CHAPTER  5.  SIMPLE  WHPAS
                                                           •

       .

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

-------
5.6.  CHOOSING  THE  RIGHT TOOL
                                                     63
                 SimpleWHPA
                    Analytic
  WhAEM
Analytic Element
                                             Back-of-the-
                                           Envelope based
                                                on  f
MODFLOW
Finite Difference
                                                                                    No
   Figure 5.11: Decision tree for complexity of capture zone delineation method, after [Haitjema,2002].

5.6   Choosing the Right  Tool

The technique  used  to delineate time-of-travel capture zones for wellhead protection  areas  is a balance
between geohydrological complexity, data availability, time constraints, budget, and other factors.  A decision
tree can  assist  in the justification of approach based on geohydrological complexity, as  shown in Figure
5.12.  If  Back-of-the-Envelope capture zones are not  justified,  then Simple WHPAs and geohydrological
modeling should be considered.  The analytic solutions  in the Simple WHPAs tool include well(s) in uniform
flow solution and simple  linear-shape boundary conditions. Geohydrological modeling using W/iAEM2000
involves solving the regional ground water flow differential equations using the analytic element  method,  and
the conceptual  possibilities include more complex river boundary conditions and piece-wise constant area
recharge distributions and hydraulic conductivity zones. Where can one get other analytic element models?
Check out the EPA frequently asked questions (#2) at http://www.epa.gov/athens/software/whaem/press.html.
If the conceptual model requires multiple aquifers, then the USGS finite difference model MODFLOW [USGS,
2005] or the analytic solvers TimML  [Bakker, 2004] or MLAEM [Strack, 2005] are  available.  If a complex
hydraulic conductivity distribution is  known, MODFLOW or a finite element solver might be the best choice.

-------
64                                                   CHAPTER, 5.  SIMPLE  WHPAS

-------
Chapter 6
SUMMARY
6.1   Good  Modeling Practice

Fortunately, we do not always have to repeat the procedure of hypothesis testing with increasingly complex
conceptual models for each individual modeling project. After a while some general conclusions 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.

6.1.1   Rules of Thumb
  1. River Proximity.  The residence time based protection zone, or capture zone, is not  likely influenced
     by surface water when the distance of the well to the feature is greater than 4v/QT/(irHn), 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 1^/kh/kvH from  a hydrological boundary, where k^
          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^kh/kvH, where kh and
          kv are the horizontal and vertical hydraulic conductivity, respectively, and where H is the average
          saturated aquifer thickness.

     Note that in the event that the (2D) capture zone width is 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.

                                                65

-------
66                                                                         CHAPTER 6.  SUMMARY

   3. Multi-Aquifer Flow.  The subsurface is usually stratified (layered) with each layer having its own
     hydraulic conductivity.  If the hydraulic conductivities of the  various formations are all within the
     same order of magnitude, the layers form one stratified aquifer, see item 4 in this list. If on the other
     hand one or more layers have a hydraulic conductivity that is more than an order of magnitude lower
     than the more permeable layers 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  hydrological  boundary (stream or lake) is present within a distance of  4\/Tc  from the well,
          where T (ft2/day} is the transmissivity of the aquifer and c (days) is the largest resistance of one
          of the adjacent  aquitards.

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

   4. Aquifer  Stratification.  An aquifer is considered  stratified  if the hydraulic  conductivities of the
     various geological formations (layers)  are different but within the same order of magnitude. 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 time-
     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 time-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
     layers.

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

   6. Transient  Flow.  Highly permeable  unconfined  aquifers,  and most confined  aquifers, will  exhibit
     "summer conditions" and "winter conditions" 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

-------
6.2.  CONCLUSION                                                                             67

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

6.2   Conclusion

W/&AEM2000 is a computer tool to support step-wise and progressive modeling and delineation of source
water areas for pumping wells. Our case study moved from arbitrary fixed radius, to calculated fixed radius.
to well In uniform flow, to modeling with hydrologic boundary conditions only, to modeling with hydrologic
and geologic boundary conditions. Each solution was conceptually more sophisticated, and we assume that
the corresponding capture zones were progressively more realistic. The emphasis of the W/iAEM2000 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 parameterizations are
also important. W/&AEM2000 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 arid wellhead
protection, it is also intended to stimulate innovation in software design  and modeling practice in the private
sector.
    One more time, "Ground water models should be as simple as possible, but not simpler." Good modeling!

-------
68                                                                CHAPTER 6.  SUMMARY

-------
Appendix A
Conversion Factors

Length



Area



Rate



Volume rate



Multiply
ft
miles
miles
yd
ft2
acre
mi2
hectare
inches/yr
inches/yr
ft/day
gpd/ft2
ft3/s
ft3/s
gal/min
gal/day
By
0.3048
5280
1609.344
0.9144
0.092903
4046.8564
2.589988 x 106
10,000
2.2815 x 10~4
6.9541 x 10~5
0.3048
0.0407
2,446.5754
86400
5.4510
0.0038
To obtain
m
ft
m
m
m2
m2
m2
m2
ft/day
m/day
m/day
m/day
m3 /day
ft3/day
m3/day
m3/day
                    69

-------
70                                           APPENDIX A.  CONVERSION FACTORS

-------
Appendix  B

W/iAEM2ooo   Base  Maps  and  the  Internet
B.I   Getting Binary Base Maps 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/iAEM2000, 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). The web page is shown in Figure B.I. 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). Click the unzip button in the Winzip program, and the bbms will
be saved to your hard drive (e.g., folder WHAEM, folder IN, folder Vincenne).  Next close the Winzip dialog
box. The binary base maps (bbm) can be 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 (optional format, 1:100,000 scale,
either NAD27  or NAD83 datum, check with USGS). The BBMs are saved in files associated with their
location and category, and the name reflects this information.  The 1:100,000  scale map is divided in eight
15-minute quads.  Table B.I.  Each of these quads contains the DLG  for the associated category: hyp-
sography, hydrography, roads, railroads,  and miscellaneous.  Table B.2.  (The EPA server has only lim-
ited availability of hypsography coverages ...  for more electronic maps check with USGS Earth Explorer
http://earthexplorer.usgs.gov). As an example, the hydrography quad for Vincennes, IN, is VI1HYF06,
which refers to the sixth quadrant (FOB)  of the  Vincennes 1:100,000 map (VII), containing hydrography
(HY).

                                        Table B.I:
F01
F05
F02
FOG
F03
F07
F04
F08
B.2   Creating BBMs from USGS DLGs maps from the USGS Web Server

W/iAEM2000 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
                                            71

-------
72
APPENDIX B.   WAEM2000 BASE  MAPS AND  THE  INTERNET
                                                                                 U.S. Environms
                     Exposure Assessment Models
         United States Binary Base Maps
          Binary base map (8BM) files ate vet-tensed graphics derrved from U S Geological Survey (USGS) digital
          hue qiaphs BBfvl files can be imported directly mto WNAEM2000 software and serve as the geographic
          basis foi gioundwatei modeling To find BBM data fos a specific geographic area, clack on the appropriate
          state below A more detailed graphic of the state you selected will display Alternatively, you can choose
          from an alphabetical %• • ; ; of available states
                            Figure B.I: The CEAM binary base map server.
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).
The  30 x 60 minute maps are stored in optional format  and SDTS format, while the 7.5 minute maps are
stored in SDTS format. 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
statgraph.cr.usgs.gov/.
   We will demonstrate the procedures for downloading  the 7.5 minute (1:24K scale)  DLGs for Vincennes,
Indiana from the USGS web site, and  convert  to BBMs. The 30x60 minute BBMs are available from the
EPA web server as described above.
   From the EROS Data center web page edc.usgs.gov/geodata/, select 1:24K DLG. You can select the DLG
maps by alphabetical List, by State, or by Graphics. Select "by State", then select Indiana.
   As an example, select the Vincennes quad, and download the hypsography DLG (1669766.HP.sdts.tar.gz)
to your project directory.
   From the W/iAEM2000 menu, select Tools  , then Convert.  For the Vincennes hypsography DLG, which
is  1:24K SDTS format, select SDTS to Basemap.  (If you were working with the USGS 1:100K optional
format DLG, you would select DLG to Basemap).  Select the compressed file (*.gz) you downloaded from
the Internet. Provide a name for the BBM. Be sure to limit  yourself to 8 characters.  W/iAEM2000 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. See Figure B.3.

-------
B.2.   CREATING BBMS  FROM  USGS  DLGS  MAPS  FROM  THE  USGS  WEB SERVER
                                                              73
                  I Products
                                 Science
                                                NASA IP (MAC
                                                             I Satellite Missions
                                                                            INSIRSDA
                                                                                          I AboilUiROS
                  USGS Geographic Data Download
For the text only page click here.
                            1:2SOK OEM 1:24K DEM 30 m NED 1:311 DLG 1:100K OLG  1:24K DLG  LULC   NtCD  NHD
                     1:24,000 Scale Digital Line Graphs (DLG) SDTS Format Only
                        •  Example
                        •  StatusJMajjs
                        •  Condensed User.Guide (Global Land Information System HTML)
                        •  Review the_ OOREADME before downloading!
                        •  SDTS format only.
                        •  Alpnabetical List
                        •  State
                        •  Graphics
                        •  SDTS DLGs require Master Data Dictionaiy
                        •  Back
                   U.S Department of the Interior | U.S Geological Survey (USGS) \ ERQS (An official website of the U.S.
                  government.)
                  URL t)tip'//e(ic.usg$.gov//geoclata/>nd8x.htmt
                  Maintainer: Contact Us via email at CustamB.LS£W£.BJsu.stsea^Msgs,ys.v.l / Wgtisils..Comraents
                  {emsweb&jsgs, gov)
                  Last Update: Fet> 16, 2005
                  Privacy Slstement   Disclaimer  Accessibility   FOIA
                   'TlRSTGOvl
                              Figure B.2:  USGS EROS Data Center Geodata portal.

-------
74
APPENDIX B.  WAEM2000 BASE MAPS AND  THE INTERNET
                                          Table B.2:
Category Name
hypsography
hydrography
boundaries
roads and trails
railroads
miscellaneous transportation
Abbreviation
hp
hy
bd
rd
rr
mt
              a . ~
                      3 ±l»l=s
        Figure B.3: Binary base map showing USGS 7-5 minute DLG for Vincennes hypsography.

B.3    USGS raster  maps  (DRG)

The USGS digital raster graphics  (DRG) and digital orthophoto quad (DOQ) maps can be selected and
ordered from earthexplorer.usgs.gov. W/iAEM2000 is now capable of using these as base maps.
B.4   USEPA BASINS maps  (Shapefile)

The USEPA watershed modeling system BASINS has data downloads organized by USGS hydrological unit
code http://www.epa.gov/OST/BASINS/.  The vector maps are available as shapefiles.  The user has the
opportunity to customize the map projection (you can save as UTM to be consistent with your W/iAEM2000
project). Interesting maps for source water protection include land use/land cover and soils maps.

-------

[Bakker, 2004] Bakker, M. (2004).  TimML: a multiaquifer analytic element model version  2.1.  University
  of Georgia, http://www.engr.uga.edu/people/mbakker/timnil.html.

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

[Blandford and Huyakorn, 1993] Blandford, T. and Huyakorn, P. (1993). WHPA: a semi-analytical model
  for the delineation of wellhead protection areas, version 2.2. http://www.epa.gov/ada/models.html.

[Ceric,  2000]  Geric, A. (2000). Assessment of the applicability of simplified capture zone delineation tech-
  niques for ground water public water supply systems. Master's thesis, School of Public and Environmental
  Affairs, Indiana University-Bloomington.

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

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

[Dupuit, 1963] Dupuit, J.  (1963).  Etudes Theoriques et Pratiques sur le Mouvement des Eaux dans les
  Canaux Decouverts et a  Trailers les Terrains Perm-cables, 2nd  ed. unknown,  Duriod, Paris.

[Feyen  et al., 2001] Feyen,  L., Beven, K.,  Smedt. F.  1).. 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 ct al., 1995] Haitjema,  H.,  Strack,   O.,   and  Kracmcr,   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.
  http://www.epa.gov/ada/download/project/wellhd.pdf.

                                                75

-------
76                                                                               BIBLIOGRAPHY

[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/whaem.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.
  http://wwwbrr.cr.usgs.gov/projects/GW.Subsurf/mchill/pubs/method/index.shtml.

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

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

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

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

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

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

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

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

-------
BIBLIOGRAPHY                                                                               77

[Straek et al., 1994]  Strack, ()., 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.
  http://www.epa.gov/ada/csmos/models/czaem.html.

[Strack. 2005] Strack. O. 1). (2005).  MLAEM - multilayer analytic element model  version 5.02.  Strack
  Consulting, Inc. http://strackconsulting.com/home/.

[USACOE, 2005]  USACOE (2005).  CMS ground water modeling system.  US Army Corps of Engineers
  http://chl.erdc.usace.army.mil/CirL.aspx7p  s&a  Software;!.

[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, 2005  USEPA  (2005).   Source water protection.  Office of  Ground Water  and Drinking Water
  http://www.epa.gov/safewater/protect.htnil. updated.

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

[USGS.  2005] USGS (2005). MODFLOW2000. http://water.usgs.gov/nrp/gwsoftware/modnow.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.

-------