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