United States Office of Ground-Water March 1991
Environmental Protection Protection
Agency
= Office of Water
&EPA WHPA
A Modular Semi-Analytical
Model for the Delineation of
Wellhead Protection Areas
Version 2.0
-------
WHP A
A Modular Semi-Analytical
Model for the Delineation of
Wellhead Protection Areas
Prepared by:
T. Neil Blandford
Peter S. Huyakorn
of:
HydroGeoLogic, Inc.
Herndon, VA 22070
For:
U.S. Environmental Protection Agency
Office of Ground-Water Protection
Washington, DC 20460
Under subcontract to
ICF. Incorporated
Contract No. 68-08-0003
March 1991
-------
Acknowledgements
>:-w»K'«»w»<*-fr«/^>w , ~-f " r f s
^^S^CvwrfvAft X* w< ^ ^ ^X-sV. ^ /X-ft-v -.S-SfjjtvM, > v/
-------
Disclaimer
The work presented in this document has been funded by the United States
Environmental Protection Agency. It has not been subject to the Agency's peer
and administrative review, and has as yet not been approved as an EPA
document. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use by the U.S. EPA. While EPA believes
that the WHPA model may be very useful for WHP program implementation,
EPA does not select, endorse or approve its use over any other approach.
In order to demonstrate use of the WHPA model, and to provide example
problems so that users may verify that the code is working correctly, several
WHPA delineation examples based upon actual hydrogeological studies in
Coming, NY, Albuquerque, NM and Seattle, WA are included in this document.
These examples are presented for instructional purposes only; due to a variety of
simplifying assumptions and inferred data, these case study examples may not be
indicative of actual site conditions.
-------
Abstract
WHPA is a modular, semi-analytical ground-water flow model designed to assist
State and local technical staff with the task of Wellhead Protection Area
(WHPA) delineation. The model is PC-based and very user-friendly.
The WHPA model consists of four independent computational modules that may
be used to delineate capture zones. Three of the modules contain semi-analytical
capture zone solutions; they are applicable to homogeneous aquifers that exhibit
two-dimensional, steady ground-water flow in an areal plane. Multiple pumping
and injection wells may be present. Barrier or stream boundary conditions which
exist over the entire aquifer depth may be simulated. Available aquifer types
include confined, leaky-confined and unconfined with area] recharge. One of the
three modules is a Monte Carlo module that provides an uncertainty analysis
capability. This module allows the effects of uncertain input parameters on the
delineated capture zone(s) to be assessed quantitatively.
The fourth module is a general particle tracking module that may be used as a
postprocessor for two-dimensional numerical models of ground-water flow.
Because this module uses the hydraulic heads output from a numerical model to
delineate capture zones, the hydrogeological scenarios that may be investigated
are only limited by the capabilities of the numerical code.
-------
Table of Contents
1.0 INTRODUCTION !-!
2.0 GETTING STARTED 2'1
2.1 WHPA Model Installation 2-1
2.2 A Brief Overview of the WHPA Model 2'3
2.3 A Brief Example 2'8
3.0 PROBLEM DESCRIPTION 3'1
3.1 Introduction 3-1
3.2 Capture Zone Types 3-1
3.2.1 Steady-State Capture Zones 3-1
3.2.2 Time-Related Capture Zones 3-4
3.23 Hybrid Capture Zones 3-4
4.0 SOLUTION TECHNIQUES 4'!
4.1 Particle Tracking 4'1
4.2 Pathline and Capture Zone Delineation 4-4
5.0 OVERVIEW OF THE WHPA MODEL 5'!
5.1 Objectives 5-1
5.2 Structure and Organization of the WHPA Model 5-1
5.2.1 Computational Modules 5-1
5.2.2 Preprocessor 5-4
5.2.2.1 Computational Modules Main Menu 5-4
5.2.2.2 Input Screens 5-6
5.2.2.3 Options Menu 5-8
5.2.2.4 Input Files 5-9
5.23 Postprocessor 5-9
5.2.3.1 GRAF Module Options 5-10
5.2.3.2 HPGL File Option 5-13
5.2.3.3 ARC/INFO File Option 5-14
5.2.3.4 General Output Files 5-14
5.3 Selection of Appropriate Modeling Options 5-15
iv
-------
6.0 RESSQC MODULE 6-1
6.1 Introduction 6-1
6.2 Capabilities 6-1
6.3 Assumptions and Limitations 6-1
6,4 Input Requirements 6-2
6.5 Example Applications 6-2
6.5.1 Hypothetical Example 6-2
6.5.2 Coming Example 6-6
7.0 MULTIPLE WELL CAPTURE ZONE MODULE (MWCAP) 7-1
7.1 Capabilities 7-1
7.2 Assumptions and Limitations 7-1
7.3 Input Requirements 7-2
7.4 Example Applications 7-2
7.4.1 Time-Related and Steady-State Capture Zone Comparison 7-2
7.4.2 Boundary Effects Example 7-5
7.4.3 Albuquerque Example 7-8
7.4.4 Seattle Example 7-11
8.0 GENERAL PARTICLE TRACKING MODULE (GPTRAC) 8-1
8.1 Introduction 8-1
8.2 Semi-analytical Option 8-1
8.2.1 Capabilities 8-1
8.2.2 Assumptions and Limitations 8-2
8.2.3 Input Requirements 8-2
8.2.4 Example Applications 8-5
8.2.4.1 Albuquerque Example 8-5
8.24.2 Seattle Example 8-7
8.3 Numerical Option 8-9
8.3.1 Capabilities 8-10
8.3.2 Assumptions and Limitations 8-10
8.3.3 Input Requirements 8-10
8.3.4 Example Applications 8-14
8.3.4.1 Hypothetical Examples 8-14
8.3.4.2 Coming Example 8-21
-------
8.4 Comparison of Semi-Analytical and Numerical Options 3-29
8.4.1 Strip Confined Aquifer Example 8-29
8.4.2 Unconfined Aquifer Examples 8-29
8.4.3 Leaky Aquifer Example 8-36
9.0 UNCERTAINTY ANALYSIS MODULE (MONTEC) 9-1
9.1 Introduction 9-1
9.2 Uncertainty Analysis Objectives 9-1
9.3 The Monte Carlo Technique 9-2
9.4 Application of Monte Carlo Technique to Capture Zone Analysis 9-3
9.4.1 Methodology 9-3
9.4.2 Input Requirements 9-5
9.4.2.1 Distribution Types 9-6
9.4.2.2 Distribution Bounds 9-8
9.4.2.3 Statistical Correlations Between Input Parameters 9-9
9.4.2.4 Data Analysis 9-12
9.4.3 output 9-12
9.5 Example Applications 9-13
9.5.1 Hypothetical Example 1 9-13
9.5.2 Hypothetical Example 2 9-16
10.0 PROBLEM CONCEPTUALIZATION vs. REALITY
POTENTIAL FOR ABUSE 10-1
10.1 Introduction 10-1
10.2 Effects of Input Parameters on Capture Zone Morphology 10-1
10.3 Effects of Boundary Conditions on Capture Zones 10-5
10.4 Analysis of Simulation Results 10-8
REFERENCES
VI
-------
Appendices
A THEORETICAL DEVELOPMENT AND IMPLEMENTATION
OF RESSQC MODULE A-l
A.I Introduction A-l
A.2 Analytical Solutions Used in RESSQC A-3
A. 2.1 Uniform Flow A-3
A.2.2 Point Sources and Sinks A-3
A.2.3 Finite Radius Source in Uniform Flow A-6
A2.4 Combination of Uniform Flow with Point Sources and Sinks A-7
A.3 Pathline and Capture Zone Delineation Algorithms A-9
A.3.1 Pathline Computation A-9
A.3.2 Front Tracking and Capture Zone Delineation A-12
A.4 Code Structure and Dimension Limits of RESSQC Module A-12
A.5 Differences Between RESSQC and the Original RESSQ Code A-12
B THEORETICAL DEVELOPMENT AND IMPLEMENTATION
OF MWCAP MODULE B-l
B. 1 Introduction B-l
B.2 Analytical Solutions Used in MWCAP B-2
B.2.1 Aquifer Without Lateral Boundary B-2
B.2.2 Aquifer With Stream Boundary B-4
B.2.3 Aquifer With Barrier Boundary B-9
B.3 Capture Zone Delineation Algorithms B-12
B.3. 1 Steady-State Capture Zones B-12
B.3.2 Time-Related Capture Zones B-16
B.3.3 Hybrid Capture Zone B-17
B.4 Code Structure and Dimension Limits of MWCAP Module B-20
vn
-------
C. GPTRAC FORMULATION C-l
C.I Introduction C-l
C.2 Theoretical Background C-l
C.2.1 Analytical Solutions used in GPTRAC C-2
C .2.1.1 Unconfined Aquifers with Recharge C-2
C.2.1.2 Leaky Semi-Confined Aquifers C-5
C.2.1.3 Strip and Three-sided Aquifers C-6
C.2.2 Darcy's Pathline Equations and Time Integration C-9
C.2.2.1 Darcy's Equations C-10
C.2.2.2 Numerical Integration C-10
C.2.2.3 Analytical Integration C-12
C.3 Numerical Implementation C-21
C.3.1 Grid Generation C-21
C.3.2 Velocity Computation C-22
C.3.3 Pathline and Capture Zone Delineation Procedure C-26
C.3.4 Finite Element Interpolation of Velocities C-30
C.4 Code Structure and Dimension Limits of GPTRAC Module C-30
D. SAVING TO WORDPERFECT D-l
D. 1 Retrieving HPGL Plotter Files Into WordPerfect 5.0 D-l
D.2 Retrieving HPGL Plotter Files Into WordPerfect 5.1 D-3
E. HEDCON UTILITY USER'S GUIDE E-l
E.I Introduction and Features E-l
E.2 Grid File Option E-3
F. DEVELOPMENT HISTORY OF WHPA CODE F-l
REFERENCES
vm
-------
List of Figures
Figure Page
3.1 Terminology for Basic Capture Zone Analysis 3-2
3.2 Generic Steady-State Time-Related, and Hybrid Capture Zone Shapes 3-3
4.1 Pathline Tracking to Delineate Time-Related (a) and
Steady-State (b) Capture Zones 4-5
5.1 Typical WHPA Model Input Screen (a) and Typical Input
Screen With The Options Menu Invoked (b) 5-7
5.2 Typical WHPA Model Plot Displayed on CRT Monitor Using
the GRAF Module 5-11
6.1 Contaminant Fronts Delineated Using RESSQC for Hypothetical Example 6-5
6.2 Capture Zones Delineated Using RESSQC for Hypothetical Example 6-7
6.3 General Site Map of Chemung River Valley in the Vicinity of Coming, N.Y 6-8
6.4 Cross Sections F-F and G-G' In the Vicinity of Corning
Reproduced from Ballaron (1988) 6-9
6.5 Predevelopment Water Table Map for Coming Region (Obtained From Numerical
Simulation) and Locations of Three Surficial Aquifer Pumping Wells 6-11
6.6 Five-Year Capture Zones Delineated Using RESSQC and Corresponding
Input Parameters for Coming Example 6-12
6.7 Five-Year Capture Zones for Corning Example Delineate Using RESSQC
With One Capture Zone Time Specified 6-13
7.1 Steady-State, and 10 and 25-Year Time-Related Capture Zones
Delineated Using MWCAP For The Time-Related Capture Zones the
Number of Pathline s is 50 7-4
7.2 Ten-Year Capture Zone For the Cases of No Boundary (a)
A 'Stream Boundary (b) and A Barrier Boundary (c) 7-7
7.3 General Ground-Water Flow System for Albuquerque MWCAP Example 7-9
7.4 Steady-State Capture Zone for Albuquerque Municipal Well Field
Delineated Using MWCAP 7-10
IX
-------
7.5 General Location Map for Seattle's Highline Well Field Example 7-12
7.6 Local Site Map With Intermediate Aquifer Well Locations for Highline
Well Field Example 7-13
7.7 Hydrostratigraphic Cross Section B-BThrough Highline Well Field Area 7-14
7.8 Contour Map of Piezometric Head for Highline Well Field Intermediate
Aquifer 7-16
7.9 Five-Year Hybrid Capture Zones Delineated Using MWCAP For
Riverton Heights and Boulevard Park Wells 7-18
7.10 Evaluation of Well Interference. Effects and Drawdown At Ground-Water
Divide Due to Each Well Using Theis Solution 7-20
8.1 Twenty-Five-Year Capture Zones for Three Albuquerque Municipal
Wells Located Near the Rio Grande 8-6
8.2 Five-Year Hybrid Capture Zones (MWCAP) and Five-Year Time-Related
Capture Zones (GPTRAC) Computed for Highline Well Field 8-8
8.3 Schematic Representation of Head Data File Format Required by GPTRAC
For Finite Element or Mesh-Centered Finite Difference Model Output
With Nodes Numbered in the y-Direction (a), or for Block-Centered
Finite Difference Model Output With Nodes Numbered in the x-direction (b) 8-13
8.4 Hypothetical Aquifer and Hydraulic Properties Used for First Two
GPTRAC Numerical Examples 8-16
8.5 One-Hundred-Day Capture Zones for Two Pumping Wells in Hypothetical Aquifer 8-19
8.6 One-hundred-Day Capture Zones for Two Wells and Three-Hundred-Day
Reverse- and Forward-Tracked Pathlines for Hypothetical Aquifer 8-20
8.7 General Site Map and Modeling Area Boundary for Surficial Aquifer in
the Vicinity of Corning, NY 8-22
8.8 Finite Difference Model Boundary Conditions and Distribution of
Hydraulic Conductivity for Surficial Corning Aquifer 8-23
8.9 Steady-State Head Field for Surficial Coming Aquifer Using MODFLOW
For (a) No Pumping Wells and (b) Three Pumping Wells 8-24
8.10 Five-Year Capture Zones of Three Pumping Wells in the Vicinity of Coming, NY 8-28
-------
8.11 Strip Aquifer Simulations Using the Semi-Analytical GPTRAC Module
for Two Barriers (a) Two Streams (b), and a Barrier (c) 8-30
8.12 Strip Aquifer Simulations Using the Numerical GPTRAC Module for Two
Barriers (a) Two Streams (b), and a Stream and a Earner (c) 8-31
8.13 Capture Zones Computed Using the Semi-Analytical (a) and Numerical (b)
GPTRAC Options for Unconfmed Aquifer with Zero Recharge
8.14 Capture Zones Computed Using the Semi-Analytical (a) and Numerical (b)
GPTRAC Options for an Unconfmed Aquifer with Areal Recharge
Capture Zones Computed Using Semi-Analytical GPTRAC Module for an
Unconfmed Aquifer with Zero Recharge and Ambient Flow and a Hydraulic
Conductivity of (a) 20 ft/d, and (b) 1 ft/d
Capture Zones for Semi-Confined (Leaky) Aquifer Using the Semi-Analytical
(a) and Numerical (b) GPTRAC Modules 8-37
9.1 Monte Carlo Method Applied to Capture Zone Analysis Series of Radial
Lines Used to Save Capture Zone Boundary (a) and Cumulative Distribution
Function For One Radial Line (b) 9-4
9.2 Distorted Shape of MONTEC Capture Zone 9-14
9.3 5th - 90th, and 99th Percentile Capture Zones For The First
Hypothetical Example Computed Using MONTEC 9-17
9.4 The 50th, 80th, and 95th Capture Zone Percentiles Calculated Using
MONTEC For The Second Hypothetical Example 9-19
9.5 Capture Zones Delineated Using MWCAP for (a) Lowest Randomly
Generated Transmissivity Value and (b) Highest Randomly Generated
Transmissivity Value 9-22
9.6 The 50th, 80th and 95th Capture Zone Percentiles Computed Using MONTEC
for the Second Hypothetical Example with Vertical Leakage 9-23
10.1 Mass Balance Analysis For Simple Hytrogeological System 10-3
10.2 Miscellaneous Capture Zones 10-6
10.3 Stream Boundary As Represented by the WHPA Model (a), and Stream
Boundary As Often Encountered in Practice (b) 10-7
XI
-------
10.4 Barrier Boundary As Represented by the WHPA Model (a), and Barrier
Boundary That May Be Encountered In Practice (b) . 10-9
A.I Fundamental Flow Fields Uniform Flow Field (a), and Radial Flow Field (b) A-4
B.I Schematic Flow System For MWCAP Stream Boundary Formulation B-5
B.2 Typical Pumping Well Capture Zones for a Well Near A Stream in an
Aquifer With Ambient Flow a) Pumping Rate Below Critical Value
b) Pumping Rate at Critical Value; c) Pumping Rate Above Critical
Value (from Lee 1986) B-7
B.3 Schematic Flow System For MWCAP Barrier Boundary Formulation B-10
B.4 Typical Steady-state Capture Zones That Can Be Expected For The
Situation Of An Aquifer With A Stream Boundary And Various Angles
Of Ambient Flow (Newsom and Wilson 1988)' B-13
B.5 Iterative Predictor-Corrector Particle Tracking Procedure B-15
B.6 Capture Zone Boundary Truncation Error That May Occur Using
RESSQC (a) and Improved Algorithm Used by MWCAP (b) B-l8
B.7 Schematic Description of the Procedure Used to Locate the Capping
Boundary of a Hybrid Capture Zone B-19
B.8 Simplified Flow Chart for MWCAP Module B-22
C.I Rectangular Grids Used by GPTRAC; Mesh-Centered Grid (a)
and Block-Centered Grid (b) C-l 1
C.2 Principal Velocity Components of a Rectangular Grid Block or Element C-14
C.3 Schematic Representation of a Particle Path Through A Grid Block C-16
C.4 Possible Patterns of the Principal Velocity Components Along the x-Direction C-20
C.5 Zonal Representation of an Inhomogeneous Aquifer C-23
C.6 Notation Used In Velocity Computation For (a) Grid Block i, and (b) Element e C-25
C.I Algorithm To Determine If A Particle Has Been Intercepted By A Nearby Well
(a) Raise Flag and (b) Continue C-27
C.8 Capture Zone Delineation By Reverse Pathline Tracking C-29
C.9 Simplified Flow Chart For GPTRAC Module C-33
xii
-------
List of Tables
Table
2.1
2.2
2.3
5.1
5.2
5.3
6.1
7.1
7.2
7.3
8.1
8.2
8.3
8.4
9.1
9.2
9.3
A.I
A.2
B.I
B.2
C.I
C.2
E-l
Graphics Devices supported by WHPA
Description of WHPA Model Computational Modules
Required Input for WHPA Model Computational Modules
Description of WHPA Model Computational Modules
Required Input for WI-WA Model Computational Modules
GRAF Module Options
RESSQC Input Requirements For Delineation of Capture Zones (Input
Differs Slightly For Delineation of Contaminant Fronts)
Input Requirements for MWCAP Module
MWCAP Input Parameters for Second Hypothetical Example
MWCAP Input Parameters used for Highline Well Field Example
Input Requirements for GPTRAC Semi-Analytical Option
Input Requirements for GPTRAC Numerical Option
Input Parameters for GPTRAC Hypothetical Aquifer (Examples One and Two)
Input Parameters for GPTRAC Numerical Option for Coming Example Problem
Input Parameters for First MONTEC Example
Input Parameters For Second MONTEC Example
MWCAP Input Parameters For Two Monte Carlo End-Member Capture Zones
Description of RESSQC subroutines Subroutines Listed In Alphabetical Order
Dimension Limits For RESSQC Arrays
Description of Subroutines That Compose The MWCAP Module
(Subroutines Listed In Alphabetical Order)
Dimension Limits For MWCAP Arrays
Description Of Key Subroutines In GPTRAC Module (Subroutines
Listed In Alphabetical Order)
Dimension Limits Of Key Array Variables Used In GPTRAC
Input Requirements for HEDCON Program
Page
2-2
2-5
2-6
5-2
5-3
5-12
6-3
7-3
7-6
7-17
8-3
8-11
8-17
8-26
9-15
9-18
9-21
A-13
A- 14
B-21
B-23
C-31
C-34
E-2
xiii
-------
1.0 Introduction
The Amendments to the Safe Drinking Water Act (SDWA), which were passed in June
1986, established the first nationwide program to protect ground-water resources used for
public water supplies from all (anthropogenic) potential threats. Unlike previous Federal
programs, that have tended to focus on individual contaminant sources, this new effort
approaches the assessment and management of ground-water quality from a more compre-
hensive perspective. The SDWA seeks to accomplish this goal by the establishment of State
Wellhead Protection (WFP) Programs that "protect wellhead areas within their jurisdiction
from contaminants which may have any adverse effect on the health of persons." A WHP
Program is part of a State's Ground Water Protection Strategy.
One of the major elements of WHP is the determination of zones within which
contaminant source assessment and management will be addressed. These zones, called
Wellhead Protection Areas (WHPAs), are defined in the SDWA as "the surface and
subsurface area surrounding a water well or wellfield, supplying a public water system,
through which contaminants are reasonably likely to move toward and reach such water well
or wellfield."
The States are given flexibility in determining appropriate operational approaches to
WHPA delineation. The U.S. Environmental Protection Agency (EPA) is required by the
SDWA to provide technical guidance on the hydrogeologic aspects of this task. In June
1987, EPA published the technical background document "Guidelines for Delineation of
Wellhead Protection Areas." In this report, EPA outlined five criteria that can be used as
a technical basis for WHPA delineation:
Distance
Drawdown
Time of travel
Flow boundaries
Assimilative capacity
-------
1-2 Introduction
States may select one, or some combination, of the above criteria to form a technical basis
for their WHP program. A State's choice of a criteria will likely be based on a combination
of technical and nontechnical (e.g., administrative) considerations.
Delineation methods are used to translate the criteria selected by the States to actual,
mappable delineation boundaries. EPA based upon existing ground-water protection
programs in the United States and Western Europe, has identified six primary methods for
WHPA delineation. The methods are listed below in order of increasing technical
sophistication:
Arbitrary fixed radii
Calculated fixed radii
Simplified variable shapes
Analytical methods
Hydrogeologic mapping
Numerical flow/transport models
A detailed explanation of each of these methods may be found in Chapter 4 of the EPA
Guidelines document (U.S. EPA 1987).
The development of a user-friendly computer model to assist State and local technical
staff with the delineation of WHPAs is an outgrowth of EPA's Guidelines document. This
report documents the capabilities and proper use of the WHPA model, which is a modular
semi-analytical and numerical code for the delineation of WellHead Protection Areas. Two
of the five delineation criteria, time of travel and flow boundaries, may be addressed using
the model. Although one model option is based on numerical model results, the delineation
methods used by the model are primarily semi-analytical.
All users are strongly encouraged to read this manual thoroughly prior to applying the
WHPA model. Chapter 2 is designed for those of you who refuse to do so ~ if you are an
experienced ground-water flow modeler it should provide enough information to get you up
and running. Chapter 3 is a general overview of the capture zone delineation problem,
while chapter 4 is an introduction to one of the most commonly used techniques to delineate
capture zones (particle tracking). Chapter 5 is an in-depth overview of the WHPA model.
Chapters 6-9 describe in detail the capabilities and limitations of, and provide examples for,
the individual capture zone delineation options contained in the WHPA model. Chapter 10
is the last and perhaps most important chapter; it describes some of the errors that may
occur if the WHPA model is applied improperly.
-------
2.0 Getting Started
2.1 WHPA Model Installation
The WHPA model is designed to run on any standard, IBM or compatible XT, AT,
or 386 microcomputer (PC) operating under MS DOS version 2.1 or later. The machine
must have a minimum of 640K RAM (Random Access Memory), a hard disk, and CGA
EGA VGA or Hercules graphics capability (EGA or VGA is highly recommended). The
model will automatically detect and use the math coprocessor if one is present but one is
not required.
A hard copy of capture zone plots observed on the monitor may be obtained using
EPSON, OKIDATA NEC PinWriter and IBM Proprinter dot matrix printers and
compatibles, Hewlett-Packard (HP) LaserJet laser printers and compatibles, and Hewlett-
Packard and Houston Instruments pen plotters and compatibles. A summary of the
available graphics devices supported by WHPA is provided in Table 2.1. If you do not have
access to any of the supported hard copy output devices, you may save plot files in standard
HPGL format (see section 5.2.3.2) and subsequently use any number of software packages,
such as WordPerfect (see Appendix D), PrintAPlot, or just about any CAD package, to
output the plot on other output devices.
To install the WHPA model on your computer, perform the following steps:
1. Establish a sub-directory on the hard disk
MKDHt ffHPA
2. Change directory
CD\HHPA
3. Place diskette in drive A
-------
2-2 Getting Started
Table 2.1
Graphics Devices Supported by WHPA
Device Type
Video
Dot Matrix Printer
Laser Printer
Pen Plotter
Standard Models Supported
CGA
EGA
VGA
Hercules Graphics Card (HGC)
EPSON FX, MX, LQ800, LQ1000, LQ1500,
LQ2500 and compatibles
OKIDATA and compatibles
NEC PinWriter and compatibles
IBM Proprinter and compatibles
HP LaserJet and compatibles (such as HP
DeskJet)
Hewlett-Packard and other HPGL-compatibles
Houston Instruments and other DM/PL-
compatibles
4. Copy files on WHPA diskette to hard disk. Note that steps 3 and 4 must be
performed two times.
COPY A:*.*
Once step 4 is completed the following files should exist in your WHPA subdirectory.
WHPA£XE
GRAFJEXE
PRERQC^XE
RESSQC.EXE
PREMW.EXE
MWCAP3XE
PKEGP.EXE
GPTRAC.EXE"
PREMC.EXE
MONTECJEXE
SETUP.EXE
HEDCON.EXE
-------
Getting Started 2-3
The next step is to execute the setup program by typing
SETUP
The setup program will display a series of menus that prompt for information
concerning graphical output. The program first attempts to determine the type of graphics
card installed in your computer; if the program determines the correct graphics card type,
type "V to select the default and continue on to select a hard copy output device. If you
type "N", a menu will appear that allows selection of the appropriate graphics mode (CGA
EGA, VGA or Hercules). Once the correct graphics card type is determined, a menu will
appear from which any of the hard copy output devices listed in Table 2.1 may be selected.
You must also select the port that the device is connected to (LPT1, LPT2, PRN, COM1
or COM2) and possibly some additional parameters that are device specific (e.g., the
number of pens available for a given pen plotter). The communication parameters specific
to your output device, such as baud rate, should be set using the MS DOS "MODE
command. The MODE command will generally be located in your AUTOEXEC.BAT file.
Refer to your MS DOS Users Manual for more details.
To execute the WHPA mode, simply type "WHPA".
WHPA
2.2 A Brief Overview of the WHPA Model
The primary objective of the WHPA model is to assist State and local technical staff
with the task of WHPA delineation. The WHPA model is an easy-to-use, widely applicable
tool for WHPA delineation based on state-of-the-art technology of the ground-water
industry. The WHPA model can be divided conceptually into two major sections. The
computational modules section contains the Fortran programs that compute the capture
zone(s) for a given physical scenario. All of the "number crunching" is performed by these
computational modules. The remaining portion of the WHPA model is the user-interface.
The interface provides an efficient mechanism for data entry as well as the viewing of model
results.
-------
2-4 Getting Started
The WHPA model contains four major computational modules: RESSQC, MWCAP,
GPTRAC, and MONTEC. The capabilities of each of these modules is summarized in
Table 2.2. Some users may recognize RESSQC as a modified version of the computer code
RESSQ presented by Javandel et al. (1984). The remaining three modules were developed
specifically for the EPA Office of Ground-Water Protection (OGWP) using state-of-the-art
technology and some recently published studies available in the literature (e.g., Newsom and
Wilson, 1988 and Pollock, 1988). The capabilities, assumptions, limitations, and input
requirements for each of the computational modules are discussed in detail in Chapters 6-9.
A matrix of the input requirements for each of the computational modules is provided in
Table 2.3. The MONTEC module is not listed in Table 2.3; it has the same input
requirements as MWCAP and semi-analytical GPTRAC with the addition that the uncertain
aquifer parameters and their associated probability distributions must be specified. Note
that each of the computational modules operate entirely independent of one another.
There' are two major assumptions common to all of the computational modules; 1)
flow in the aquifer is at steady state, and 2) flow in the aquifer is horizontal (two-
dimensional in areal view). The first assumption implies that the aquifer is under
equilibrium conditions, and therefore temporal variations in sources and sinks (including
pumping) are not considered. The WHPA model is therefore most applicable to
continuously used water-supply wells. For RESSQC, MWCAP, and MONTEC, the second
assumption implies that the aquifer is confined, or unconfined if the drawdown-to-initial
saturated thickness ratio is small (approximately less than 0.1). This assumption is also
applicable for the confined aquifer option in GPTRAC, but this module also has a special
unconfined aquifer option that may handle drawdown ratios much larger than 0.1 with
minimal error. None of the modules simulate vertical flow of water within the aquifer
explicitly.
The preprocessor is very straightforward to use. The user is prompted, through a
series of pop-up windows, for input required by the selected computational module. The
following key sequences will be useful when using the preprocessor:
Menu Commands:
H - Invoke a series of 1-4 pop-up help screens that define model input
parameters and provide guidance on proper model option selections.
M - Return to the main model menu. Any model input parameters that
were entered, or any changes made to an existing data set will be
saved. May also use to Home or End keys.
-------
Getting Started
2-5
Table 2.2
Description of WHPA Model Computational Modules
Module Name
Description
RESSQC
MWCAP
GPTRAC
MONTEC
Delineates time-related capture zones around pumping wells, or
contaminant fronts around injection wells, for multiple pumping
and injection wells in homogeneous aquifers of infinite areal extent
with steady and uniform ambient ground-water flow. Well interfer-
ence effects are accounted for.
Delineates steady-state, time-related or hybrid capture zones for
pumping wells in homogeneous aquifers with steady and uniform
ambient ground-water flow. The aquifer may be infinite in areal
extent or the effects of nearby stream or barrier boundaries can be
assessed. If multiple wells are examined, the effects of well inter-
ference are ignored.
Semi-analytical Option: Delineates time-related capture zones for
pumping wells in homogeneous aquifers with steady and uniform
ambient ground-water flow. The aquifer may be of infinite areal
extent, or it may be bounded by one or two (parallel) stream
and/or barrier boundaries. The aquifer may be confined, leaky
confined or unconfined with areal recharge. Effects of well inter-
ference are accounted for.
Numerical Option: Delineates time-related capture zones about
pumping wells for steady ground-water flow fields. Since this op-
tion performs particle tracking using a head field obtained from a
numerical (finite difference or finite element) ground-water flow
code, many types of boundary conditions as well as aquifer hetero-
geneities and anisotropies may be accounted for.
Performs uncertainty analysis for time-related capture zones for a
single pumping well in homogeneous aquifers of infinite areal
extent. The aquifer may be confined or leaky confined.
-------
2-6
Getting Started
Table 23
Required Input for WHPA Model Computational Modules
Required Input
Units used
Aquifer type*
Study area limits
Maximum step length
No. of pumping wells.
No. of recharge wells
Well locations
Pumping/injection rates
Aquifer transmissivity
Aquifer porosity
Aquifer thickness
Angle of ambient flow
Ambient hydraulic gradient
Areal recharge rate
Confining layer hydraulic
conductivity
Confining layer thickness
Boundary condition type
Perpendicular distance from
well to boundary
Orientation of boundary
Capture zone type
No. of pathlines used to
delineate capture zones
Simulation time
Capture zone time
Rectangular grid parameters
No. of forward/reverse pathlines
Starting coordinates for
forward/reverse pathlines
Nodal head values
No. of heterogeneous
aquifer zones
Heterogeneous aquifer properties
RESSQC
MWCAP
GPTRAC
Semi-
analytical
Numerical
Confined, unconfmed or leaky-confined.
-------
Getting Started
2-7
Menu Commands (cont'd):
N
P
S \
Q I
Fl:
Input Field Commands:
=
;
Move to the next input screen. May also use the Page
Up key.
Move to the previous input screen. May also use the
Page Down key.
Obtain a summary of the input parameters.
Quit editing/input session and return to main model
menu. Entered model parameters or changes made to
an existing file will not be saved.
Exit to-DOS. To return to the WHPA model, type
EXIT at the DOS prompt.
Clear the entire input field and reposition the cursor to
the beginning of the field. .
< ........ .....
Standard destructive backspace. \
Accept the current value; in an input field, and position
the cursor at the start of the next input field. If an
input field is left blank, a value of zero is assigned to
the corresponding modej parameter. One beep will
sound if inappropriate input is entered.
2.3 A Brief Example
The best way to introduce the WHPA model is to go through a step-by-step example.
This section is designed to do just that, so sit down at your computer and type "WHPA'.
The first screen that will appear is the disclaimer:
« DISCLAIMER **
Neither the U.S. EPAs Office of Ground-Water Protection, HydroGeoLogic,
Inc. , nor any person acting on behalf of either of these entities:
a) makes any warranty, express or implied, with respect to this
software; or
b) assumes any liabilities with respect to the use or misuse of
this software, or the interpretation or misinterpretation of
any results obtained from this software, or for damages
resulting from the use of this software.
-------
2-8
Getting Started
The disclaimer screen may be erased by pressing any key on the keyboard. The next screen
is the main menu for the WHPA code.
MAIN MENU FOR THE
WHPA CODE
RESSQC Modified version of the RESSQ code
MWCAP Capture zones for single or multiple,
(non-interfering) production wells
GPTRAC General particle tracking
MONTEC "Monte Carlo analysis
QUIT Halt execution of WHPA program
[t4.]move cursor, to select option
This screen allows the user to select one of five available options. Select the MWCAP
option by moving the cursor down one line using the down arrow key and depressing
"Enter." The next screen is the main menu for the MWCAP option.
-- MWCAP OPTION --=
(N)ew Problem
(C)ontinue Current Problem
(S)ave Current Problem
(L)oad Previous Problem
(R)un Model
(P)lot Results
(Directory
(H)elp
(E)xit
-------
Getting Started
2-9
Press "N" to input the parameters for a new problem. Next a series of input screens will
appear that prompt you for required model input. Type in the numeric values shown on
each screen.
MWCAP
Run Title: BRIEF WHPA MODEL EXAMPLE
Units to use for Current Problem: O
0 - meters and days
1 - feet and days
Number of Wells for which
Capture-Zones are desired: 1
Minimum X-Coordinate:
Maximum X-Coordinate: 3000
Minimum Y-Coordinate:
Maximum Y-Coordinate: 3000
Maximum Spatial Step Length: 50
Change Any Values On This Screen (Y/N)? i.
- select value - options menu = DOS shell
At the bottom of each input screen, the user may change values entered on the screen by
typing "Y". If the input values are all correct, type "N" to continue to the next screen.
AQUIFER PROPERTIES AND LOCATION FOR WELL # 1
X Coordinate (L): 500
Y Coordinate (L): 1500
Well Discharge Rate (L**3/T): 4000
Transmissivity (L**2/T): 1000
Hydraulic Gradient (dimensionless): 0.00150
Angle of Ambient Flow (degrees): 180
Aquifer Porosity (dimensionless): 0.25
Aquifer Thickness (L): 50
(H)ELP (M)AIN (N)EXT (P)REVIOUS (S)UMMARY (Q)UIT
- select value - options menu - DOS shell
-------
2-10 Getting Started
The options menu that is displayed at the bottom of the above input screen is obtained by
pressing the Escape key at any time. This menu allows you to obtain HELP information,
save the values that have been entered and return to the MAIN MWCAP menu, go to the
NEXT input screen, go to the PREVIOUS input screen, obtain a SUMMARY of
information already entered, or QUIT the current application (entered input data not
saved). Press "H" for help, and the following menus will appear:
MWCAP HELP
Well Location:
The location of each well is specified by a Cartesian coord-
inate pair (x,y). Well locations must reside within the study
area defined on the previous MWCAP input screen, or an error
will occur.
Well Discharge:
The pumping value for each well must be in ft**3/day or
m**3/day. Because the flow field is assumed to be at steady
state, only one discharge value per well is required. MWCAP
assumes that all wells are pumping wells. The sign of the
pumping value does not matter, it is set in the code.
Transmissivity:
The transmissivity (T) of an aquifer is a measure of the ease
with which water can travel through the porous media. T is
often computed from the equation T=Kb, where K-hydraulic
conductivity and b=aquifer thickness. The units of T are
ft**2/day or m**2/day.
MWCAP HELP
Gradient:
The hydraulic gradient (ft/ft or m/m = dimensionless) is most
commonly measured from a map of piezometric surface or water
table elevations. The average ambient gradient should be in-
put to the model, and therefore gradients prior to pumping, or
gradients not affected by the cone of depression should be used.
Direction of Ground-Water Flow:
Ground water flows from areas of high hydraulic head towards
areas of low hydraulic head; for homogeneous, isotropic
aquifers the direction of ground-water flow is perpendicular to
the hydraulic head contours. At a given site, the direction of
ground-water flow may be variable; in this case the average,
most dominant direction should be used. The direction of flow
may be 0-360 degrees, with 0=due east, 90=due north, etc.
Press any key to continue
-------
2-12 Getting started
CAPTURE-ZONE TYPE OPTION FOR WELL # 1
Capture-Zone Type Option: 2
0 - steady-state
1 - hybrid
2 - time-related
Travel Time: 3650
Number of Pathlines Desired: 20
(default- 20)
Plot Capture Zone Boundary ? 1
(O-No, 1-Yes)
Change Any Values On This Screen (Y/N)?
- select value - options menu - DOS shell
The above input screen is the last one for MWCAP. After entering the appropriate values
and typing "N", the main MWCAP menu will reappear.
MWCAP OPTION
(N)ew Problem
(C)ontinue Current Problem
(S)ave Current Problem
(L)oad Previous Problem
(R)un Model
(P)lot Results
(Directory
(H)elp
(E)xit
-------
Getting Started 2-11
MWCAP HELP
Porosity:
Porosity (dimensionless) is defined as the volume of the
voids within the aquifer divided by the total volume of the
aquifer, it must always be less than one by definition, and
values of 0.15-0.30 are characteristic of most aquifers.
Thickness:
The aquifer thickness has units of ft or m. If the aquifer
has a variable thickness, an average value for the aquifer
(generally in the vicinity of the pumping well) should be
used.
Press any key to continue
At any point while these help screens are displayed, pressing the Escape key will return you
to the appropriate model input screen, while pressing any other key will advance you to the
next help screen. Once you have entered the required values for the present input screen,
you will move through two more input screens in sequence
BOUNDARY CONDITION INPUT FOR WELL # 1
Boundary Type: 0
0 -
1 -
2 -
no boundary
stream boundary
barrier boundary
Change Any Values On This Screen (Y/N)?
= select value = options menu = DOS shell
-------
Getting Started 2-13
At this point type 'R" in order to execute the MWCAP module. The message will appear
** MBCAP
**
When the MWCAP run is completed, the main MWCAP menu will appear again. At this
point, type "P" to plot the results of the MWCAP run. Hit "Enter" to plot the current graph,
and the following figure will appear
(H)
3000
2400
1800
1200
BOO -
(H)
0 600 1200 1800 2400 3000
(S)ave Plot (H)ard Copy (O)verlay (R)etrieve Plot (M)ap Scale (E)xit
At this point the plotting module (GRAF) options may be used to save the plot file in
ASCII, HPGL or ARC/INFO format, to obtain a hard copy if you have access to one of the
output devices listed in Table 2.1, to overlay the results of up to fifteen different WHPA
model runs (this can not be done at this point because only one plot file has been created
thus far), to retrieve a previously created plot file, to scale the plot arbitrarily or to
correspond to some map scale, or you may exit the plotting module.
To exit the brief WHPA model example, type "E' to exit from the plotting module,
type another "E' to exit from the main MWCAP menu, and finally move the cursor to the
"Exit" option on the main WHPA model menu. Users are strongly encouraged to read the
remainder of the user's guide prior to applying the WHPA model for capture zone
delineation purposes.
-------
3.0 Problem Description
3.1 Introduction
A capture zone is defined as the zone surrounding a pumping well that will supply
ground-water recharge to the well. The zone of contribution (ZOC) of a well is identical
to the capture zone of a well as defined in OGWP's "Guidelines for Delineation of Wellhead
Protection Areas" (U.S. EPA 1987). For two-dimenensional areal ground-water flow
problems, the capture zone corresponds to the area of contribution surrounding the well.
A capture zone for a single well in a homogeneous, isotropic aquifer with an ambient
uniform flow is shown in Figure 3.1. Note that the extent of the capture zone in the
downgradient direction is defined by the location of a stagnation point. Water particles
between the well and the stagnation point travel towards the well in a direction opposite to
the regional hydraulic gradient. Water particles downgradient of the stagnation point travel
in the direction of regional flow, even though they may be located within the cone of
depression of the pumping well. The stagnation point itself is defined mathematically as a
point of zero ground-water flow velocity.
3.2 Capture Zone Types
The WHPA model may be used to delineate three types of capture zones; steady-state,
time-related, and hybrid. Steady-state and hybrid capture zones can only be obtained using
the MWCAP option, while the remaining options as well as MWCAP may be used to obtain
time-related capture zones. This section describes the three types of capture zones.
3.2.1 Steady-State Capture Zones
A steady-state capture zone is the surface or subsurface area surrounding a pumping
well that will supply ground-water recharge to the well over an infinite period of time. The
typical outline of a steady-state capture zone for a single pumping well is depicted in Figure
3.2. The open-ended shape of the capture zone is due to the fact that given enough time,
any particle of water upstream of the well, within the capture zone
-------
3-2 Problem Description
Figure 3.1
Terminology for Basic Capture Zone Analysis
Ground-water
divide
Cross section
Plan View,
Capture Zone
-------
Problem Description 3-3
Figure 3.2
Generic Steady-State, Time-Related, and Hybrid Capture Zone Shapes
Capture Zone Types
Steady-State
Time-Related
AMBENT
FLOW
Hybrid
-------
3-4 Problem Description
boundaries, will eventually travel to the well. In practice, the upstream end of the capture
zone would be "capped" in some manner due to physical and/or managerial restrictions. For
example, a steady-state capture zone may terminate at a ground-water flow divide.
3.2.2 Time-Related Capture Zones
A time-related capture zone is the surface or subsurface area surrounding a pumping
well that will supply ground-water recharge to the well within some specified period of time.
When calculating time-related capture zones, OGWP generally recommends that time
periods of 10-25 years be considered.
A typical outline of a time-related capture zone for a single pumping well is depicted
in Figure 3.2. A time-related capture zone is always represented by some closed shape. In
general, time-related capture zones are less conservative (enclose smaller areas) than steady-
state or hybrid capture zones. As the specified time increases, however, differences between
the three capture zone types in the proximity of the pumping well become negligible.
Note that time-related capture zones may be calculated when the ground-water flow
field is at steady-state. A steady flow field implies that the direction and magnitude of the
ground-water flow velocity at any point within the aquifer is constant for all time; this
concept should not be confused with the fact that it takes some finite period of time for a
water particle within a capture zone to travel to a pumping well located within a steady-state
flow field.
3.2.3 Hybrid Capture Zones
As the name implies, a hybrid capture zone is a combination between a steady-state
and a time-related capture zone. A typical outline of a hybrid capture zone is shown in
Figure 3.2. The hybrid capture zone is identical to the steady-state capture zone in all
respects except that it is "capped' on the upstream end (Figure 3.2). A particle of water
released from the mid-point of the capping segment will reach the pumping well within some
specified time. Therefore, the cap on the hybrid capture zone approximates a segment of
a time-related capture zone. The hybrid capture zone can be viewed as an implementable
alternative to the steady-state capture zone.
-------
4.0 Solution Techniques
The WHPA model delineates capture zones about pumping wells using the particle
tracking technique. The term "particle" is used only for conceptual purposes. One may view
a particle as an individual water molecule or an individual molecule of a conservative tracer
that moves through the aquifer coincident with the bulk movement of ground-water flow
dispersion and diffusion do not affect the particle location.
To obtain steady-state or hybrid capture zones, particles are released from the
stagnation point(s) of the system. Time-related capture zones are obtained by tracing the
pathlines formed by a series of particles placed around the well bore of the pumping well.
The code uses both forward and reverse particle tracking depending upon the option(s)
selected. Forward tracking involves the tracking of particles in the direction of ground-water
flow, while reverse tracking involves the tracking of particles in the direction opposite to
ground-water flow. The following two sections provide an introduction to the delineation
of capture zones using the particle tracking technique. For a more detailed explanation,
refer to Appendix A.3.
Although the term pathline is used throughout this document, it should be noted for
completeness sake that pathlines correspond to streamlines for the case of steady ground-
water flow. The term streamline, rather than pathline, is used often in the relevant
literature.
4.1 Particle Tracking
The particle tracking method requires knowledge of the ground-water flow velocity at
any point within the aquifer. The flow velocities are expressed in terms of Darcy's Law,
which may be written as:
Q = KiA (4-1)
where Q is the volumetric flow rate, K is the hydraulic conductivity of the porous medium,
i is the hydraulic gradient (change in hydraulic head over some specified horizontal distance)
and A is the cross-sectional area of flow. The specific discharge (or Darcy velocity) is
defined as:
-------
4-2 Solution Techniques
q = Q/A = Ki (4-2)
The average pore-water velocity for an individual fluid particle moving through the porous
medium may be written as
v = q/fl (4-3)
where v is the seepage velocity and 9 is the effective porosity of the medium. Equation (4-
3) may be generalized to describe the x and y components of seepage velocity for two-
dimensional, horizontal flow,
vx = qJ0 v = q/0 (4-4)
x y
Several methods are available to obtain the seepage velocity components, vx and v«
for a given flow field. For the RESSQC, MWCAP, MONTEC and semi-analytical GPTRAC
options of the WHPA model, the required velocities are obtained analytically. That is, there
are exact mathematical solutions for see page velocity programmed into the Fortran code.
Using these solutions, the code solves for vx and v at any specified location C^i) within the
aquifer.
The numerical option of the GPTRAC module uses a slightly different method to
calculate velocities. This option requires that hydraulic head values at the nodes of a
rectangular grid be supplied to the code. The grid may be a finite element or a finite
difference grid. In the latter case, the nodes may be located at the grid-block centers (block-
centered grid) or at the intersection of the grid lines (mesh-centered grid). Velocities at any
location are then calculated using a simple analytical solution within each grid block (see
Appendix C.3). The computed velocities are dependent upon the nodal hydraulic head
values, and the grid block hydraulic conductivity and effective porosity values.
Once velocities can be determined, pathlines (the route that an individual particle of
water follows through an aquifer) may be delineated using particle tracking. Particle
tracking delineates pathlines by calculating the distance d£, that is traversed in a given time
dt. The distance a particle travels during a given period of time is defined by
d£ = (dx2 + dy2)'72 (4-5)
-------
Solution Techniques
where
dx = vxdt = qxdt/0 (4-6a)
dy = vydt = qydt/0 (4-6b)
where dx and dy are the projections of d£ on the x and y axis, respectively. In practice,
equations 4-5 and 4-6 are solved using some form of numerical integration; i.e. the
differential dt is approximated by a finite time step At, and the differential d£ is approximat-
ed by a finite spatial increment A£. Therefore, the path of a water particle may be traced
using
= *i + AX = *i + vxAt (4-7a)
= y* + Ay = yi+ vyAt (4-7b)
where (x^) is the location of the water particle at time t, and (Xj+ 1}yi+ j) is the position of
the particle at time t + At.
The terms forward and reverse particle tracking are used frequently throughout this
document. Forward tracking refers to the procedure of tracking water particles in the
direction of ground-water flow, while reverse tracking refers to the procedure of tracking
water particles in the direction opposite to that of ground-water flow. Since ground water
flows towards a pumping well, reverse tracking is used when particles are released about the
circumference of a well bore. If particles are released upgradient of a pumping well,
forward tracking is used to determine whether or not the particle will be captured by the
well. Forward tracking is based upon equations (4-7a) and (4-7b), and reverse tracking is
based upon these equations with the sign of the velocity terms vx and v reversed (multiplied
by negative one).
Forward tracking can be used to determine whether or not a pumping well will be
contaminated by a particular contaminant source. For example, particles released at the
edge of a waste disposal facility may be forward-tracked for a specified time to determine
if they will enter a pumping well. If a production well has been contaminated, reverse
tracking may be used to determine potential sources of contamination. To do this, particles
would be released at the contaminated well and reverse-tracked through time to identify
potential sources of the contaminated water.
-------
4-4 Solution Techniques
4.2 Pathline and Capture Zone Delineation
Using the particle tracking method, pathlines within a region can be delineated for any
specified travel distance or time. The general procedure is illustrated in Appendix B
(Figure B.5).
Time-related capture zones are delineated by placing a series of water particles
(generally about 20-50) at sequential locations along the perimeter of a small circle
representing the well boundary. Individual pathlines for each of these particles are then
traced using reverse tracking. Pathlines are terminated when the assigned travel time value
is. reached or when a plotting boundary is encountered (Figure 4.1a). The capture zone
consists of the entire region enclosed by the delineated pathlines.
To delineate the boundaries of steady-state and hybrid capture zones, it is convenient
to release particles from the stagnation point(s). Pathlines that are forward tracked will
terminate at the pumping well, while pathlines that are reverse tracked will end at a study
area boundary (Figure 4.1b). The pathlines that emanate from the stagnation point(s) form
the steady-state capture zone boundary.
Note that in the case of Figure 4.1b where the capture zone intercepts the stream
boundary, the stream itself forms a segment of the capture zone boundary, and the two
pathlines extending from the stagnation points to the well partition the well recharge
attributable to the stream and the ambient aquifer flow respectively. The portion of the
capture zone attributable to the stream is referred to as the "stream capture zone".
The perimeter of a well bore and stagnation points are convenient locations to release
particles for capture zone delineation purposes; however, it is often desirable to release
particles from other locations within the study area (see previous section on reverse and
forward particle tracking). The RESSQC and GPTRAC options of the WHPA model allow
the user to specify arbitrary starting particle locations. The particles may be either forward
or reverse tracked.
-------
Solution Techniques 4-5
Figure 4.1
Pathline Tracking to Delineate Time-Related (a)
and Steady-State (b) Capture Zones
Y
X
Stream
Capture Zone
Boundary
\
Reverse Tracked
Pathline
Stagnation
Point
Ambient Row
(a)
Reverse Tracking
Stagnation
Point
Stream Capture
Zone
Forward Tracking
Aquifer Capture
Zone
Ambient Flow
(b)
-------
5.0 Overview of the WHPA Model
5.1 Objectives
The primary objective of the WHPA model is to assist State and local technical staff
with the task of WHPA delineation. The WHPA model is an easy-to-use, widely applicable
tool for WHPA delineation based on state-of-the-art technology of the ground-water indus-
try. The model is designed for widespread use among technical staff who may have only a
basic understanding of ground-water flow processes and computer model applications on the
IBM PC.
5.2 Structure and Organization of the WHPA Model
The WHPA model can be divided conceptually into two major sections. The computa-
tional modules section contains the Fortran programs that compute the capture zone(s) for
a given physical scenario. All of the "number crunching" is performed by these com-
putational modules. The remaining portion of the WHPA model is the user-interface. The
interface provides an efficient mechanism for data entry as well as the viewing of model
results. The composition of, and the relationship between, these two portions of the model
is discussed in the next two sections.
5.2.1 Computational Modules
The WHPA model contains four major computational modules: RESSQC, MWCAP,
GPTRAC, and MONTEC. The capabilities of each of these modules is summarized in
Table 5.1. Some users may recognize RESSQC as a modified version of the computer code
RESSQ presented by Javandel et al. (1984). The remaining three modules were developed
specifically for OGWP using state-of-the-art technology and some recently published studies
available in the literature (e.g., Newsom and Wilson, 1988). The capabilities, assumptions,
limitations, and input requirements for each of the computational modules are discussed in
detail in Chapters 6-9. A matrix chart of the input requirements for each of the computa-
tional modules is provided in Table 5.2. The MONTEC module is not listed in Table 5.2;
it has similar input requirements as MWCAP and semi-analytical GPTRAC with the addition
that the uncertain aquifer parameters and their associated probability distributions must be
specified. Note that each of the computational modules operate entirely independent of one
another.
-------
5-2
Overview of the WHPA Model
Table 5.1
Description of WHPA Model Computational Modules
Module Name
Description
RESSQC
MWCAP
GPTRAC
MONTEC
Delineates time-related capture zones around pumping wells, or
contaminant fronts around injection wells, for multiple pumping
and injection wells in homogeneous aquifers of infinite areal extent
with steady and uniform ambient ground-water flow. Well interfer-
ence effects are accounted for.
Delineates steady-state, time-related or hybrid capture zones for
pumping wells in homogeneous aquifers with steady and uniform
ambient ground-water flow. The aquifer may be infinite in areal
extent or the effects of nearby stream or barrier boundaries can be
assessed. If multiple wells are examined, the effects of well inter-
ference are ignored.
Semi-analytical Option: Delineates time-related capture zones for
pumping wells in homogeneous aquifers with steady and uniform
ambient ground-water flow. The aquifer may be of infinite areal
extent, or it may be bounded by one or two (parallel) stream
and/or barrier boundaries. The aquifer may be confined, leaky
confined or unconfined with areal recharge. Effects of well inter-
ference are accounted for.
Numerical Option: Delineates time-related capture zones about
pumping wells for steady ground-water flow fields. Since this op-
tion performs particle tracking using a head field obtained from a
numerical (finite difference or finite element) ground-water flow
code, many types of boundary conditions as well as aquifer hetero-
geneities and anisotropies may be accounted for.
Performs uncertainty analysis for time-related capture zones for a
single pumping well in homogeneous aquifers of infinite areal
extent. The aquifer may be confined or leaky confined.
-------
Overview of the WHPA Model
5-3
Table 5.2
Required Input for WHPA Model Computational Modules
Required Input
Units used
Aquifer type*
Study area limits
Maximum step length
No. of pumping wells
No. of recharge wells
Well locations
Pumping/injection rates
Aquifer transmissivity
Aquifer porosity
Aquifer thickness
Angle of ambient flow
Ambient hydraulic gradient
Areal recharge rate
Confining layer hydraulic
conductivity
Confining layer thickness
Boundary condition type
Perpendicular distance from
well to boundary
Orientation of boundary
Capture zone type
No. of pathlines used to
delineate capture zones
Simulation time
Capture zone time
Rectangular grid parameters
No. of forward/reverse pathlines
Starting coordinates for
forward/reverse pathlines
Nodal head values
No. of heterogeneous
aquifer zones
Heterogeneous aquifer properties
RESSQC
MWCAP
GPTRAC
Semi-
analytical
Numerical
* Confined, unconfined or leaky-confined.
-------
5-4
Overview of the WHPA Model
There are two major assumptions common to all of the computational modules; 1) flow
in the aquifer is at steady-state, and 2) flow in the aquifer is horizontal (two-dimensional in
areal view). The first assumption implies that the aquifer is under equilibrium conditions,
and therefore temporal variations in sources and sinks (including pumping) are not
considered. The WHPA model is therefore most applicable to continuously used water-
supply wells. For RESSQC and MWCAP, the second assumption implies that the aquifer
is confined, or unconfmed if the drawdown-to-initial saturated thickness ratio is small
(approximately less than 0.1). This assumption is also applicable for the confined aquifer
option in GPTRAC, but this module also has a special unconfined aquifer option that may
handle drawdown ratios much larger than 0.1 with minimal error. None of the modules
simulate vertical flow of water within the aquifer explicitly.
5.2.2 Preprocessor
The desired computational module (RESSQC, MWCAP, GPTRAC or MONTEC) is
selected from the main WHPA model menu that appears immediately after the disclaimer
screen when the code is run. Each of the modules have preprocessors that provide for the
convenient, interactive input of parameters. Parameter values are entered through a series
of input screens that appear on the monitor. Each input screen has one or more associated
help screens that define the input parameters and guide the user in the selection of
appropriate modeling options. The user is prompted only for those parameter values that
are required for a given application.
5.2.2.1 Computational Modules Main Menu
The "main menu" for each of the computational modules will appear after the desired
module is selected from the main WHPA model menu. The available options are as follows:
(N)ew Problem
(C)ontinue Current Problem
(S)ave Current Problem
(L)oad Previous Problem
(R)un Model
(P)lot Results
(D)irectory
(H)elp
(E)xit
The desired option is invoked by depressing the letter key marked by parenthesis.
-------
Overview of the WHPA Model 5.5
The "New Problem" option initiates a new run of the computational module. In this
case, all of the problem input parameters will need to be entered at the appropriate input
screen.
The "Continue Current Problem" option may be used when a problem has been run,
and the user wishes to rerun the code using the same, or similar, simulation options and/or
input parameters. When this option is chosen, the input parameters (e.g., well location) and
simulation options (e.g., capture zone type) will be echoed to the screen as each input
window appears. The user may then selectively change any input value, and input
parameters for the entire problem need not be reentered.
The "Save Current Problem" option is useful for saving the input associated with a
particular capture zone delineation scenario. The input is saved in a user-specified file.
Users should not confuse this option with the "Save Plot" option available in the GRAF
module. The "Save Plot" option saves the actual capture zone plot file output by a
computational module; the "Save Current Problem" option saves only the input parameters
required to calculate the capture zone(s).
The "Load Previous Problem" option is used to retrieve applications that were saved
previously using the "Save Current Problem" option.
Note that the "Save Current Problem" and "Load Previous Problem" options support
the use of full DOS pathnames. Therefore, files may be saved to, or retrieved from, drives
and directories other than the current drive and the current directory.
The "Run Model" option initiates execution of the selected computational module (e.g.
GPTRAC). The computed pathlines and capture zones are automatically written to the
WHPA.PLT file for subsequent plotting using the "Plot Results" option. Another, more
extensive output file that contains an echo of the input parameters as well as the model
output is written for each run (section 5.2.3.4).
Once an application has been run, the "Plot Results" option will produce a plot of the
calculated capture zone(s) on the screen. When this option is selected, the file WHPA.PLT
is plotted by default, or other plot files may be retrieved or overlayed.
-------
5-6 Overview of the WHPA Model
The "Directory" option may be used to display a listing of the files that reside in any
directory. When this option is selected the user is prompted for a pathname for which a
directory is to be displayed. The default (Enter key) is to display the contents of the current
directory. The directory option also supports use of the DOS wildcard character (*) in file
names.
The "Help" option will display a brief description of the other options, and the "Exit"
option returns the user from the selected computational module to the main menu of the
WHPA code.
5.2.2.2 Input Screens
If the "New Problem", "Continue Current Problem" or "Load Previous Problem" are
selected from the main menu of the executable module, the user is prompted for the
required model parameters through a series of input screens. A typical input screen is
shown in Figure 5.la; this particular screen was selected from the MWCAP module
preprocessor. For the application shown, the user has already input the location of pumping
well number one, and the next value that must be entered is the well discharge rate. At this
point the user simply types in the discharge value. Only numerals and decimal points will
be accepted - letters or other symbols will not appear on the screen. The arrow keys may
be used to move to the left or right of the current cursor position; they may not be used to
move the cursor up or down. The input field width is 10 spaces for real values and 1-5
spaces for integer values. If the cursor is moved past the end of this invisible "box", it will
appear at the first space again. Real and integer values are distinguished automatically by
the program, and therefore a decimal point does not need to be typed if a value is real. If
a mistake is made typing in the value, the Backspace or the Delete key may be used. The
Delete key clears the entire input field and positions the cursor at the beginning of the field.
Finally, the Enter key will enter the typed value into the program and position the cursor
at the beginning of the next input field. If an input field is left blank, a zero value will be
assigned to the corresponding parameter if a non-zero default is not specified.
Certain input values are not admitted if they are inappropriate for the current
parameter. For example, WHPA will not accept negative values for transmissivity or aquifer
thickness, and porosity values must be between zero and one.
-------
Overview of the WHPA Model
5-7
Figure 5.1
'Typical WHPA Model Input Screen (a) and
Typical Input Screen With The Options Menu Invoked (b)
AQUIFER PROPERTIES AND LOCATION FOR WELL # 1
X Coordinate (m) : 500
Y Coordinate (m): 1500
Well Discharge Rate (m**3/d): _
Transmissivity (m**2/d):
Hydraulic Gradient (dimensionless):
Angle of Ambient Flow (degrees):
Aquifer Porosity (dimensionless):
Aquifer Thickness (m):
- select value - options menu = DOS shell
(a)
AQUIFER PROPERTIES AND LOCATION FOR WELL # 1
X Coordinate (m): 500
Y Coordinate (m): 1500
Well Discharge Rate (m**3/d):
Transmissivity (m**2/d):
Hydraulic Gradient (dimensionless) :
Angle of Ambient Flow (degrees) :
Aquifer Porosity (dimensionless) :
Aquifer Thickness (m) :
(H)ELP (M)AIN (N)EXT (P)REVIOUS (S)UMMARY
(Q)UIT
= select value = options menu = DOS shell
(b)
-------
5-8 Overview of the WHPA Model
If the "Load Previous Problem" or the "Continue Current Problem" is selected, the
values of each input parameter and model option that exist in the current data set will
appear automatically on the screen. These values may be changed either by typing a new
value, in which case the previous value is erased, or by using the arrow keys to selectively
change one or more digits in an input field. In any case, the value that exists on the screen
when the Enter key is depressed will be the value that is read by the code.
Another useful option available by typing the Fl key is the ability to exit to DOS.
When this option is invoked, the WHPA model resides in memory, and you may return to
the model at any time by typing "EXIT" (this procedure is the same as that used by
Lotus 123 and other commercial software packages). Exiting to DOS may be useful if you
want to view certain files or directories at the DOS level, but you do not wish to exit the
WHPA model to do so. If this option is invoked, remember to return to the WHPA model
at some point and exit the code in a normal fashion.
Once all of the required parameters on a given screen have been entered, the prompt
"Change Any Values on This Screen (Y/N)?" will appear. If a "Y' is typed, the cursor will
be repositioned at the beginning of the first input field so that one or more parameter values
can be changed. The Enter key can be used to move down the screen to each parameter
input field. If a "N" is typed, the next input screen will appear.
5.2.23 Options Menu
The options menu can be accessed from any input screen by depressing the Escape
key. Figure 5.1b shows the same input screen as in Figure 5.la after the Escape key has
been depressed. The options available are as follows:
Help - Invoke a series of i-4 pop-up help screens to assist the user in
selecting appropriate parameter values and modeling options.
Main - Exit the current application and return.to the main menu of the
current computational module. Changes (if any) to the current data
set will be saved. May also use the Home or End keys.
Next - Move to the next input screen. May also use the Page Down key.
Previous - Move to the previous input screen. May also use the Page Up key.
Summary - Obtain a summary of the input values associated with the current
application and each well.
Quit - Stop execution of the current computational module. Changes (if
: any) to the current data set will not be saved.
-------
Overview of the WHPA Model
5-9
The desired function is selected by typing the first letter of the corresponding word (e.g.,
type 'H' to access the help screens). Depressing the Escape key again will clear the options
menu. Note that you can perform some of the functions on the options menu by using the
Home, End, Page Up or Page Down keys directly.
5.2.2.4 Input Files
Each time one of the WHPA code modules is used to run a problem, the parameter
values input via the preprocessor are saved in a file in the current directory. The default
file names are as follows: ,
Executable Module
Input File Name
RESSQC
MWCAP
GFTRAC
MONTSC
RESSQCSAV
MWCAESAV
GPTRACSAV
MONTECSAV
These files are used by some of the main menu options. For example, if the user has run
an application using MWCAP, and then depresses "S" for "Save Current Problem", the file
MWCAP.SAV would be copied to the user-specified file name. If the user wishes at a later
time to rerun the problem, the file could be retrieved using the "Load Previous Problem"
option, and the file that the input parameters were stored in would be copied to
MWCAP.SAV.
5.2.3 Postprocessor
The postprocessor of the WHPA code is a plotting routine named GRAF. When any
of the executable modules has been successfully run, the GRAF module may be invoked to
plot the capture zone results on the monitor. The capture zones for multiple wells will be
plotted in different colors if a color monitor is used, or in different shades of amber (or the
default display color) if a monochrome monitor is used. If, in addition to capture zone
analysis, forward- and/or reverse-tracked pathlines are specified, all of the forward-tracked
pathlines and/or all of the reverse- tracked pathlines will be plotted in the same color. A
hard copy of the plot can be obtained from any of the printer or plotter output devices listed
in Table 2.1.
-------
5-10 Overview of the WHPA Model
5.23.1 GRAF Module Options
The plot of an example application of MWCAP is shown in Figure 5.2 as it would
appear on a monitor with EGA resolution (colors not reproduced). Once the capture zones
are plotted, the user may select one of the six listed options to continue execution of the
code. The function of each option is documented in Table 5.3. The options are invoked by
typing the first letter of the description (e.g., type "S' to save the current plot file).
The default plot file generated by WHPA is WHPA.PLT. This file is copied to a user-
specified file name if the "Save Plot" option is invoked. The file can be replotted using the
new file name and the "Retrieve Plot" option. The "Overlay option retrieves up to 15 plot
files consecutively, and plots them one on top of the other. When the overlay option is used,
the plotting information for the overlay plot is saved in the file WHPABLD. Therefore, the
"Save Plot" option may be used subsequent to the "Overlay" option to save a series of related
plots that have been overlayed. The overlay feature is particularly useful for sensitivity
analysis, where the effect of varying input parameter values (e.g., transmissivity) on the size
and shape of a capture zone is investigated.
The output from all four computational modules is stored in a compatible format, and
therefore plots produced using different modules may be overlayed. However, prior to
overlaying two plots, GRAF checks if the length units (ft or m) for each plot are compatible.
If they are not, an error message is displayed and the user must choose another GRAF
option. If two plot files contain different values for the minimum and maximum x and y
coordinates, the smallest values will be used for the minimum coordinates and the largest
values will be used for the maximum coordinates when the plots are overlayed. This
procedure ensures that a portion of a plot is not erased inadvertently.
The "Map Scale" option allows arbitrary scaling of the plot or automatic scaling of the
plot to one of five standard USGS topographic map scales. The map scales available are
listed in Table 5.3. This option enables the user to transfer capture zone results directly
from WHPA model hard copy output to project base maps. Note that the physical area
available for plotting varies with the selected output device, and therefore capture zones that
are large relative to the selected map scale may be truncated in the x and/or y plotting
directions. Most of the available output devices (standard EPSON, OKIDATA NEC and
IBM dot matrix printers and HP LaserJet laser printers) are limited by a maximum plot size
of 10 x 8 inches. The wide carriage version of the dot matrix printers have a maximum plot
size of 13.6 x 10 inches. The available plot size for the Hewlett-Packard and Houston
Instruments pen plotters varies widely with the model type.
-------
Overview of the WHPA Model 5_\ \
Figure 5.2
Typical WHPA Model Plot Displayed on CRT Monitor Using the GRAF Module
(M)
3000
2400
1800
1200 -
600 -
600
1200
1800
2400
(M)
3000
(S)ave Plot (H)ard Copy (O)verlay (R)etrieve Plot (M)ap Scale (E)xit
-------
5-12 Overview of the WHPA Model
Table 53
GRAF Module Options
Option
Function
Save Plot
Hard Copy
Overlay
Map Scale
Retrieve Plot
Exit
Saves the current plot file (WHPA.PLT) to a user-specified
file name. File may be saved in ASCII, HPGL, or
ARC/INFO format.
Generates a hard copy of the capture zone plot on any of the
output devices listed in Table 2.1.
Plots up to 15 plot files, one on top of the other. The user is
prompted for the names of the files to plot. Plot will not be
generated if the units specified for each plot file are not
compatible.
Allows the user to scale the plot (up or down) to any arbitrary
ratio or to one of 5 standard USGS topographic map scales
(7.5 minutes, 7.5 x 15 minutes, 15 minutes, 30 x 60 minutes, 1
x 2 degree).
Plots a single, user-specified plot file.
Exit the GRAF module and return to the main menu of the
current computational module.
-------
Overview of the WHPA Model 5.13
The "Hard Copy" option allows the user to obtain a capture zone plot on EPSON,
OKIDATA NEC Pinwriter or IBM Proprinter (or compatible) dot matrix printers; HP
LaserJet or compatible laser printers; and Hewlett-Packard and Houston Instruments pen
plotters (or compatibles). The output device and the appropriate port specification (LPT1,
LPT2, COM1 or COM2) are selected using the SETUP program (see section 2.1). The dot
matrix and laser printers are monochrome output devices and only one color is used for
plotting. Most pen plotters have six or more pens available. If a monochrome plotting
mode was selected in the setup routine, only pen 1 of the plotter will be used to plot the
WHPA model results. If the color plot option was selected, the axis and borders of the plot
will be plotted using pen 1, but the capture zones and pathlines will be plotted using the
remaining pens in sequential order (i.e., pen 2 for the first capture zone, pen 3 for the next
one, etc.).
It should be noted that when the GRAF module is first invoked the six options
described above appear in menu form prior to the actual plot of the capture zone(s). To
proceed with the plot at this point, simply depress the Enter key. However, you may also
invoke any of the options listed prior to plotting the graph. For example, you may use the
"overlay" option to enter two or more file names, and then the overlayed plot will appear
when the capture zone(s) are plotted. This procedure is useful because the GRAF module
options may be invoked prior to viewing the plot.
The GRAF module is designed in such a way that scale distortions will not occur when
capture zones are plotted: 1 ft or m along the x-axis will always equal 1 ft or m along the
y-axis. For example, if the x-axis represents a distance of 6,000 ft when plotted, and the y-
axis represents a distance of 4,000 ft, the length of the plotted y-axis will be equal to two-
thirds the length of the plotted x-axis. The scaling is precise when model results are plotted
on an output device, but only approximate when results are plotted on the monitor.
5.2.3.2 HPGL File Option
The "Save Plot" option allows the user to save plot files in Hewlett-Packard Graphics
Language (HPGL) format. HPGL is a commonly used standard in the computer graphics
industry, and consequently many word processing, desktop publishing, computer aided design
(CAD) and specialized computer graphics packages have the capability to import and
process HPGL files. Most of these software packages also have the capability to output
graphics images (capture zone plots) on numerous types of laser and dot-matrix printers.
Step-by-step instructions for incorporating a HPGL file into WordPerfect are presented in
Appendix D. This option was included to aid users that do not have access to any of the
-------
5-14 Overview of the WHPA Model
output devices listed in Table 2.1, or who may wish to incorporate capture zone plots directly
into reports or documents.
Users are cautioned that the original plot scale may not be maintained by all software
packages. If secondary software is used to import and subsequently output WHPA model
HPGL files, the size of the resulting plot should be carefully checked.
5.233 ARC/INFO File Option
A geographic information system (GIS) can be a valuable tool in a WHP program.
GIS technology is commonly used to integrate and analyze numerous types of spatial data.
A particularly common use of a GIS is to develop maps or "coverage". A coverage may
consist of the spatial representation of a single attribute, such as land use, or of multiple
attributes, such as land use, aquifer type, and municipal well locations. Schooley (1989) used
a GIS in a pilot project to eliminate wells from the WHP program that were not
immediately susceptible to contamination to prioritize the remaining wells for WHPA
delineation; and to target problem sites for more intensive monitoring, enforcement, and
remediation efforts.
The "Save Plot" option provides a facility for saving WHPA model results in an ASCII
file that may be used as input to the ARC/INFO proprietary GIS developed by the
Environmental Systems Research Institute (ESRI). The ARC/INFO file option is designed
to construct files compatible with the ARC GENERATE LINES command. Note that once
an ARC/INFO input file is obtained using WHPA, additional processing may be required
to convert the WHPA model Cartesian coordinates to the appropriate coordinate system
(perhaps latitude and longitude or state plane) used for the GIS application. Such
processing should be performed by a GIS specialist.
Although the WHPA model supports output specifically for the ARC/INFO software,
it is a relatively simple matter to process the standard WHPA model plot files for input into
other GIS packages.
-------
Overview of the WHPA Model 5-15
5.2.3.4 General Output Files
Each time an executable module of the WHPA code is run, a default ASCII output
file is generated. The output file names are as follows:
Executable Module
RESSQC
MWCAF
GPTRAC
MONTEC
Output File Name
RESSQCOUT
MWCAKOUT
GFTRACOUT
MOIfTECOUT
Each output file contains an echo of the problem-specific input parameters, and the output
(capture zone and pathline coordinates) calculated by the code. Refer to Chapter 9 for a
description of additional files created when MONTEC is used.
53 Selection of Appropriate Modeling Options
Selection of the appropriate modeling options depends on the problem at hand.
Tables 5.1 and 5.2 can be used as guides in determiningg the appropriate computational
module to use. For example, if a user wishes to delineate capture zones for multiple, closely
spaced pumping wells in a well field, the RESSQC or GPTRAC semi-analytical modules
should be used because the effects of well interference are accounted for. If a user wishes
to delineate steady-state or hybrid capture zones, the MWCAP module must be used. If
recharge wells exist and their effect on the flow field is deemed important, the RESSQC or
GPTRAC semi-analytical modules should be used. If observed hydraulic head values in the
field or output heads from a two-dimensional numerical model are available, the GPTRAC
numerical option may be used.
To better understand the capabilities and limitations of each computational module,
as well as simplifying assumptions that may or may not be appropriate, users are strongly
encouraged to work through the practical examples presented for each module at the end
of Chapters 6-9.
-------
6.0 RESSOC Module
6.1 Introduction
The RESSQC module is a slightly modified version of the RESSQ code presented by
Javandel et al., 1984. The "C" was added to the program name because the code was
specifically modified for "C"apture zone delineation purposes. RESSQ was originally
designed to delineate contaminant fronts about injection wells, and this capability is retained
in RESSQC.
6.2 Capabilities
RESSQC can be used to delineate time-related capture zones for a system of one or
more pumping/injection wells that fully penetrate a homogeneous aquifer. The presence of
a stream or barrier boundary can be simulated using image well theory, but the image wells
need to be supplied explicitly to the code. If a stream or barrier boundary exists, the user
may save time and avoid mistakes by using the MWCAP or semi-analytical GPTRAC option.
The effects of well interference are accounted for through superposition of solutions. A
maximum of 50 pumping wells and 20 injection wells may be used in RESSQC.
The number of pathlines reverse tracked from each pumping well may be defined
interactively by the user. In addition, particles can be released at any point within the
system to be subsequently forward or reverse tracked (depending upon whether contaminant
fronts or capture zones are delineated).
6.3 Assumptions and Limitations
Capture zones and contaminant fronts delineated using RESSQC are valid for fully
penetrating pumping and injection wells screened in aquifers that are essentially homogeneo-
us. Ground-water flow must be two-dimensional in the x-y plane, and therefore the aquifer
may be confined or unconfined if the drawdown-to-initial saturated thickness ratio is small
(less than approximately 0.1). Steady ground-water flow is assumed.
-------
6-2 RESSQC Module
If a stream or a barrier boundary is implemented using image well theory, the
boundary is assumed to be a straight line and fully penetrate the aquifer. The latter
assumption is often violated in cases where stream boundaries exist. The effect of a partially
penetrating stream may be an important one and each application should be examined on
a site-by-site basis. In general, the greater the depth and breadth of the stream in relation
to the aquifer thickness, the more valid the fully penetrating stream assumption. Also,
stream boundary partial penetration effects decrease as the distance from the stream to the
well increases. The stream and the aquifer are assumed to be in perfect hydraulic
connection, the effects of a "clogging layer" between the streambed and the aquifer are not
considered.
If, in actuality, the stream is partially penetrating and/or there is a clogging layer of fine
grained material that lines the streambed, the capture zone obtained using RESSQC will be
smaller than the "true" capture zone. The amount of error incurred will be dependent upon
the degree to which the above assumptions are violated.
6.4 Input Requirements
The input requirements for the RESSQC option are outlined in Table 6.1. The well-
specific parameters must be input for each well specified in the study area.
6.5 Example Applications
Two examples are presented that illustrate use of the RESSQC module. The first
example is a hypothetical one taken from Javandel et al. (1984). The second example is a
site application in the vicinity of Corning, New York.
6.5.1 Hypothetical Example
Because the RESSQC module is an amended version of the RESSQ code (Javandel
et al. 1984), the example problems presented in Javandel et al. were run using RESSQC to
ensure that the changes made to the code did not alter its original capability. The original
RESSQ code has been verified, field validated and used extensively by numerous users (van
der Heijde and Beljin, 1988).
The second example problem presented by Javandel et al. is depicted in Figure 6.1.
In this example there is one injection well and one pumping well with equal rates of injection
and pumping. Uniform ambient flow within the aquifer occurs at an angle of 45°
-------
RESSQC Module
6-3
Table 6.1
RESSQC Input Requirements For Delineation of Capture Zones
(Input Differs Slightly For Delineation of Contaminant Fronts)
Program Variable
Description
For each problem
IUNIT:
NWP:
NWI:
XMIN:
XMAX:
YMIN:
YMAX:
TRANSM:
GRADNT:
ALPHA:
FOR:
HEIGHT:
DL
TMAX:
NCAPZ:
NRPATH:
Default units of input parameters (feet and days or meters
and days)
Number of pumping wells within the study area
Number of recharge (injection) wells within the study area
Minimum x-coordinate of study area (ft or m) '
Maximum x-coordinate of study area (ft or m)
Minimum y-coordinate of study area (ft or m)
Maximum y-coordinate of study area (ft or m)
Transmissivity of the aquifer (ftVd or mVd)
Regional hydraulic gradient (ft/ft or m/m)
Angle of ambient ground-water flow (0-360°)
Aquifer porosity (dimensionless)
Aquifer saturated thickness (ft or m)
Largest allowable step length, dt (ft or m) see section 4.1
Maximum amount of time for calculating the trace of a pathline
(days)
Number of time-related capture zones to be calculated for each
pumping well (maximum = 7)
Number of reverse-tracked pathlines started at arbitrary locations
within the study area
For each capture zone (1=1, NCAPZ) ~
DATE(I): Time value for capture zone (days)
-------
6-4
RESSQC Module
Table 6.1 (cont'd)
RESSQC Input Requirements For Delineation of Capture Zones
(Input Differs Slightly For Delineation of Contaminant Fronts)
Program Variable
Description
For each pumping well (1=1, NWP)
XW(1,I): x-coordinate of well (ft or m)
XW(2,I): y-coordinate of well (ft or m)
QPWELL(I): Well discharge rate2/ (ftfyd or m3/d)
RADW(I): Well radius (ft or m) - default = 0.3 ft or 0.1 m
ITRW(I): Ratio of number of pathlines to the number plotted (=1 to plot
all pathlines, =2 to plot every other pathline, etc.)
NPATH(I): Number of pathlines to be computed to delineate time-related
capture zone (default = 20)
For each injection well (I=NWP + 1, NWP + NWI) -
XW(1,I): x-coordinate of well (ft or m)
XW(2,1): y-coordinate of well (ft or m)
QPWELL(I): Well recharge rate^ (ft3/d or m3/d)
RADW(I): Well radius (ft or m) - default = 0.3 ft or 0.1 m
For each forward tracked pathline (1=1, NRPATH) ~
RSTART(I,1): x starting coordinate (ft or m)
RSTART(I,2): y starting coordinate (ft or m)
The sign (+,-) of the discharge or recharge rate need not be specified.
-------
RESSQC Module
6-5
Figure 6.1
Contaminant Fronts Delineated Using RESSQC for Hypothetical Example
(M)
1000
600 -
200 -
-200 -
-600 -
-1000
P = PUMRNG WELL
I = INJECTION WELL
-1000
-600
(M)
-200
200
600
1000
Simulation
Options
Units - m and days
Step Length = 10
Simulation time = 73,000
No, Fronts = 3
@ 183, 730 and 1,460 days
Aquifer
Properties
T
b
:B
i
a
- 100
= 10
= 0.25
= 0.00343
= 45
Injection
Well
X = -300 .
Y= 300
Q = 1,200
r = 1.0
No. Pathlines = 40
Plotting Interval « 4
Pumping
Well
X=- 300
y = -300
Q = 1,200
r = 1,0
-------
6-6 RESSQC Module
with respect to the x-axis. Note that the values of transmissivity (T) and hydraulic gradient
(i) agree with the 0.14 m/d pore water velocity specified by Javandel et al.
Because the original RESSQ code was designed to track the location of contaminated
water emanating from injection wells, option 2 of RESSQC (delineate contaminant fronts
about injection wells) was used to obtain Figure 6.1. In addition, 40 pathlines (plotted at
intervals of four) were specified. A comparison of Figure 6.1 with Figure 19 in Javandel et
al. (1984) shows a perfect correspondence. The three concentric rings about the injection
well in Figure 6.1 represent the areal extent of the injected water front for the time periods
of 0.5, 2 and 4 years.
To test the capture zone delineation capability of the RESSQC code, the problem
specified above was rerun with the only change being that simulation option 1 was used
(delineate capture zones about pumping wells). The results of this run are shown in Figure
6.2. As expected, Figure 6.2 is a rotated image of Figure 6.1. In Figure 6.2, the concentric
rings about the pumping well represent the 0.5-, 2- and 4-year capture zones for the
pumping well. Note that when delineating contaminant fronts the forward tracking method
is used (particles are tracked in the direction of ground-water flow), while when delineating
capture zones the reverse tracking method is used (particles are tracked in the direction
opposite to that of ground-water flow).
6.5.2 Coming Example
For the next example, RESSQC was used to delineate the capture zones for three wells
withdrawing water from the surficial aquifer in the vicinity of Corning, New York. The data
for this example was taken from Ballaron (1988).
Figure 6.3 is a general site map of the study region. The valley sediments in the
vicinity of Coming consist of stratified glacial drift deposits that are primarily interbedded
silty to clean sands and gravels. Relatively thin deposits of lacustrine clay, silt and fine sand
exist over much of the valley and separate a surficial, unconfmed aquifer from a confined
to semi-confined aquifer at depth. Two cross sections through the study area are presented
in Figure 6.4.
For the current example three wells screened in the surficial aquifer are considered.
Recharge from the Chemung River, local recharge from precipitation, and leakage between
the surficial and lower aquifer units were neglected. Hence, the capture zones computed
-------
RESSQC Module
6-7
Figure 6.2
Capture Zones Delineated Using RESSQC for Hypothetical Example
(M)
1000
600 -
200 -
-200 -
-600 -
-1000
P » PUMPING WELL
I = INJECTION WELL
-1000
-600
-200
200
600
(M)
1000
Simulation
Options
Units - m and days
Step Length = 10
Simulation time = 73,000
No. Fronts = 3
@ 183, 730 and 1,460 days
Aquifer
Properties
T = 100
b = 10
0 = 0.25
i = 0.00343
a = 45
Injection
Well
X=-30Q .
Y= 300 '
Q - 1,200 :
T = 10
No. Pathlines - 40
Plotting Interval - 4
Pumping
Well
X- 300
Y = -300
Q = 1,200
r:= 1.0
-------
6-8 RESSQC Module
Figure 6.3
General Site Map of Chemung River Valley in the Vicinity of Corning, New York
Corning, New York
1 875 1,750 FEET
-------
RESSQC Module
6-9
Figure 6.4
Cross Sections F-F' and G-G' in the Vicinity of Corning
Reproduced from Ballaron (1988)
1050 -
a
UJ
§
o
UJ
UL
ul
o
t
850 -
750
VERTICAL EXAGGERATION «10X
1000ft
Ul
I
ui
u.
Ul
o
F
G
1150 -
1050 -
950 -
850 -
G'
750
VERTICAL EXAGGERATION * 10X
10001
EXPLANATION
H FLOOD - PLAIN SILT AND SAND
. OUTWASH AND ICE - CONTACT
^ SAND AND GRAVEL
[el LACUSTRINE SILT AND SAND
g WELL OR TEST BORING AND
IDENTIFICATION NUMBER
|dj TILL
(] BEDROCK
[D FILL
-------
6-10 RESSOC Module
using RESSQC are conservative estimates (they are larger in areal extent than they would
be if all sources of recharge to the well were incorporated into the analysis).
The locations of the three pumping wells are superimposed on a pre-purnping steady-
state head map (refer to Section 8.3.4.2 for details on how the head map was obtained) in
Figure 6.5. From this map the hydraulic gradient («0.0019) and the direction of ambient
ground-water flow (-45°) was obtained. Additional required information such as reasonable
pumping rates for the wells, an average hydraulic conductivity of 150 ft/d, and an average
saturated thickness of 25 ft were taken from Ballaron (1988).
The five-year capture zones for the three pumping wells delineated using RESSQC,
as well as the input data used, are illustrated in Figure 6.6. Note that in this example
interference effects between the three wells are significant; the two wells farthest upgradient
(about in the middle of the study area) intercept a major portion of the ambient aquifer flow
that would otherwise supply recharge to the third well. The capture
zone of the third well, therefore, assumes a complex contorted shape because it must draw
water from regions of the aquifer that are not enclosed by the capture zones of the first two
wells.
The three pumping wells examined in this exercise were assumed to be screened over
the entire saturated thickness of the surficial aquifer. Because the aquifer is unconfined, its
saturated thickness decreases as water is withdrawn by pumping wells due to dewatering of
the effective pore space. The assumptions inherent in the RESSQC, MWCAP, and
MONTEC modules are only valid for this type of system if the ratio of the drawdown at (or
near) the pumping well(s) is small compared to the initial saturated thickness of the aquifer.
In general this ratio should not exceed 20%, or 10% if one desires to be conservative. As
a GPTRAC example in Chapter 8, the Corning case study was modeled using a numerical
ground-water flow code. The results of this simulation indicate a maximum drawdown to
initial saturated thickness ratio of approximately 15-20%. This ratio is small enough that
errors incurred due to using the confined aquifer ground-water flow equations are acceptable
for demonstration purposes.
Finally, the Coming example was rerun using the input parameters from Figure 6.6,
only for this run one capture zone for 1,825 days (5 years) was specified. The results are
shown in Figure 6.7. The crescent shapes are the 5-year capture zones that RESSQC
delineates for each well.
-------
RESSQC Module
6-11
Figure 6.5
Predevelopment water Table Map for Corning Region (Obtained From Numerical
Simulation) and Locations of Three Suficial Aquifer Pumping Wells
-------
6-12
Figure 6.6
Five-Year Capture Zones Delineated Using RESSQC and
Corresponding Input Parameters for Corning Example
(FT)
9000
7200
5400
3600
1800 -
(FT)
2100
4200
6300
8400
10500
Simulation
Options
Units = ft and days
Step Length = 25
Simulation time = 1,825
No. Capture Zone times 0
Aquifer.
Properties
T 3y75Q
U = 25:
6 -0.22,
i =0.0019
a =-45*
WeU#l
X = 8,000
Y = 2,500
Q =30,000
ffo, Patalines = 10 .
Plotting Interval = 1
Well#2
6,500
4,500
30,000
10
1
Well #3
4,500
5,000
25,000
10
1
-------
RESSQC Module
(FT)
Figure 6.7
Five-Year Capture Zones for Corning Example Delineated Using
RESSQC With One Capture Zone Time Specified
9000
7200
5400
3600
1800
(FT)
2100
4200
6300
8400
10500
-------
6-14 KKSSOr, Module
These odd shapes often occur when capture zones are delineated using RESSQC for
long periods of time - they obviously do not enclose the entire capture zone of the well.
They occur because, in order to delineate time-related capture zones, the code traces a
series of reverse pathlines from each well for a specified period of time. When a capture
zone time is reached, the corresponding point on each pathline is saved. These points are
then connected in sequential order (i.e. the point on pathline 1 is connected to the
corresponding point on pathline 2, then the point on pathline 2 is connected to pathline 3,
etc.) When the last pathline is reached, its capture zone point is connected to the capture
zone point on pathline 1; this final connection causes the line to be drawn that cuts across
the pathlines and truncates the capture zone. The simplest way to avoid this shortcoming
is to set the simulation time equal to the capture zone time desired, and then specify zero
capture zone time values as was done for the first Corning run. This method produces clear
plots and there is no ambiguity as to what area is enclosed by the capture zone(s).
-------
7.0 Multiple Well Capture
Zone Module (MWCAP
7.1 Capabilities
MWCAP is designed to provide efficient delineation of steady-state, time-related and
hybrid capture zones for one or more pumping wells in homogeneous aquifers. Each well
specified may be operating in an aquifer without a lateral boundary (an areally infinite
aquifer), or in an aquifer with a stream or a barrier boundary (semi-infinite aquifer). If a
stream or barrier boundary is present, the angle of ambient flow in relation to the boundary,
as well as the orientation of the boundary itself, may be completely arbitrary. MWCAP
requires that stream or barrier boundaries be represented by straight lines in plan view.
Although multiple wells within a study area may be specified, MWCAP assumes that
the wells operate independently of one another. Therefore, physical processes such as
increased drawdown due to well interference effects are ignored.
MWCAP is very efficient due to the small number of pathlines required to delineate
steady-state or hybrid capture zones. If a stream boundary is present and the capture zone
intersects the stream, the zone of induced recharge from the stream to the well will be
delineated automatically. MWCAP can also be used to delineate time-related capture zones.
7.2 Assumptions and Limitations
Capture zones delineated using MWCAP are valid for fully penetrating pumping wells
screened in aquifers that are essentially homogeneous. Ground-water flow must be two-
dimensional in the areal x-y plane, and therefore the aquifer may be confined or unconfined
if the drawdown-to-initial saturated thickness ratio is small (less than about 0.1). A steady-
state ground-water flow field is assumed.
If a stream or a barrier boundary is present, the boundary is assumed to be linear and
fully penetrating. The latter assumption is often violated in cases where stream boundaries
exist. The effect of a partially penetrating stream may be an important one and each
application should be examined on a site-by-site basis. In general, the greater the depth and
-------
7-2 Multiple Well Capture Zone Module (MWCAP)
breadth of the stream in relation to the aquifer thickness, the more valid the fully
penetrating stream assumption. Also, stream boundary partial penetration effects decrease
as the distance from the stream to the well increases. The stream and the aquifer are
assumed to be in perfect hydraulic connection: the effects of a "clogging layer" between the
streambed and the aquifer are not considered.
If, in actuality, the stream is partially penetrating and/or there is a clogging layer of fine
grained material that lines the streambed, the capture zones obtained using MWCAP will
be smaller than the "true" capture zones. The amount of error incurred will be dependent
upon the degree to which the above assumptions are violated.
Capture zones for multiple pumping wells within a study area maybe delineated with
one run of MWCAP, but each well is assumed to operate independently of every other well.
Therefore, each well may have a potentially unique set of input parameters. The effects of
well interference (increased drawdown due to overlapping cones of depression) are
neglected.
7.3 Input Requirements
The input requirements for MWCAP are outlined in Table 7.1. Note that the well-
specific parameters must be input for each well specified in the study area.
7.4 Example Applications
In this section, four examples of capture zones delineated using MWCAP are
presented. The first example compares steady-state and time-related capture zones; the
second example illustrates the effects of boundary conditions (stream or barrier) on capture
zone morphology and the third and fourth examples demonstrate the application of
MWCAP to actual field sites near Albuquerque, New Mexico and Seattle, Washington,
respectively.
7.4.1 Time-Related and Steady-State Capture Zone Comparison
Figure 7.1 shows the steady-state, 10-year and 25-year capture zones delineated using
MWCAP for a hypothetical set of input parameters. The time-related capture zones are
enclosed entirely by the steady-state capture zone. However, as the specified time
-------
Multiple Well Capture Zone Module (MWCAP)
7-3
Table 7.1
Input Requirements for MWCAP Module
Program Variable
Description
For each problem ~
IUNIT:
NWELL:
XMEV:
XMAX
YMIN
YMAX
DLMAX
For each well (1=
XWELL(I):
YWELL(I):
QWELL(I):
TRAN(I):
GRAD(I):
ANGLE(I):
POR(I):
THICK(I):
IBOUND(I):
DSW(I):
THETA(I):
ICZTYP(I):
TMCZ(I):
NSTLIN(I):
ICZPLT(I):
Default units of input parameters (feet and days or meters and days)
Number of pumping wells for which capture zones are to be delineat-
ed
Minimum x-coordinate of study area (ft or m)
Maximum x-coordinate of study area (ft or m)
Minimum y-coordinate of study area (ft or m)
Maximum y-coordinate of study area (ft or m) .-
Largest allowable step length, dt (see section 4.1)
=1, NWELL) -
x-coordinate of well (ft or m)
y-coordinate of well (ft or m)
Well discharge rate2/ (ft3/day or m3/d)
f\ A
Transmissivity of the aquifer (ft /d or m /d)
Regional hydraulic gradient ft/ft or m/m)
Angle of ambient ground-water flow (0-360°
Aquifer porosity (dimensionless)
Aquifer saturated thickness (ft or m)
Associated boundary type (no boundary, stream boundary, or barrier
boundary)
Perpendicular distance from stream or barrier boundary to the well
(ft or m)
Orientation of stream or barrier boundary (0-360°)
Capture zone type option (steady-state, time-related, or hybrid) .
Time value associated with capture zone (days); time-related and
hybrid capture zones only
Number of pathlines to be computed for the well in addition to
pathlines delineated automatically by the code
Flag indicating if capture zone boundary is to be plotted
The sign (+,-) of the discharge rate does not need to be specified.
-------
7-4
Multiple Well Capture Zone Module (MWCAP)
Figure 7.1
Steady-State, and 10- and 25-Year Time-Related Capture Zones Delineated
Using MWCAP. For The Time-Related Capture Zones the Number of Pathlines is 50.
(M)
3000
2400
1800
1200
600
(M)
600
1200
1800
2400
3000
Simulation Options
Well No. 1
Units - m and days : X - 500
Step Length = 50 - Y = 1,500
XMIN « 0.0 XMAX * 3000 Q ~ 4rOQO
YMIN = 0.0 7MAX=3QOO T= 1,000
No. Famping Wells - 1 = i - 0,0015
; . . ..! .a = 180°
No boundary
Steady-State Capture Zone
No. of Pathlines = 0
-------
Multiple Well Capture Zone Module (MWCAP)
increases, the time-related capture zones "expand" toward the steady state solution. This
makes intuitive sense, because the steady state condition can be viewed as the case where
time is infinitely large Figure 7. 1 also indicates that the time-related capture zones most
closely resemble the steady-state solution in the region close to the well. The user may find
it beneficial to check the steady- state capture zone depicted in Figure 7.1 using the analytical
solution presented in Todd (1980) and Javandel and Tsang (1986).
Figure 7.1 was constructed by making three separate runs of MWCAP and using the
"Save Plot" and "overlay" options of the postprocessor to save the plotting information for
each run and then to subsequently overlay the three plots. Once the input for the first run
(the steady-state example) was typed in, the subsequent two scenarios were run very quickly
using the "Continue Current Problem" option in the main MWCAP menu.
7.4.2 Boundary Effects Example
The second hypothetical example is designed to illustrate the effects that aquifer
boundary conditions may have on capture zone morphology. The input parameters for this
example are summarized in Table 7.2.
Figure 7.2 illustrates the MWCAP output for this example. The first capture zone
(Figure 7.2a) was obtained for a pumping well that was not influenced by an aquifer
boundary condition. Figures 7.2b and 7.2c portray the effects of a stream and a barrier
boundary respectively on the shape and areal extent of the same capture zone. Each
respective boundary condition was oriented north to south (parallel to the y-axis) and was
located 100 m due west of the pumping well (Table 7.2).
Because the stream boundary acts as a source of water to the well, the capture zone
in Figure 7.2b is smaller than that in either Figure 7.2a or c. Since the stream boundary is
considered to be fully penetrating, the capture zone expands in width along the stream, but
it does not extend beyond the stream. When a capture zone intersects a stream boundary,
the pumping well induces flow from the stream to the well.
Figure 7.2c shows the effect of a barrier boundary near the well. Because the barrier
boundary limits the portion of the aquifer from which the well can withdraw water, the
capture zone "fans out" to capture a greater proportion of ambient flow upstream of the
well. The areas enclosed by the capture zones in Figures 7.2a and 7.2c are equal because
the only source of water to the well in each of these cases is the ambient ground-water
-------
7-6
Multiple Well Capture Zone Module (MWCAP)
Table 7.2
MWCAP Input Parameters for Second Hypothetical Example
Simulation Options
units m and days
XMIN = 0 YMIN = 0
XMAX = 4,500 YMAX = 4,500
Step Length = 50
No. of Pumping Wells = 3
Input Parameter
X Coordinate
Y Coordinate
Discharge (Q)
Transmissivity (T)
Hydraulic Gradient (i)
Angle of Ambient Flow (a)
Porosity (0)
Aquifer Thickness (b)
Boundary Type
Distance from Well to Boundary
Boundary Angle
Capture Zone Type
Travel Time Value
No. of Pathlines
Capture Zone Plotting Option
Well No. 1
1,000
3,500
4,000
1,000
0.0015
180"
0.25
50
None
N/A
N/A
Time-Related
3,650
20
Yes
Well No. 2
1,000
2,300
4,000
1,000
0.0015
180"
0.25
50
Stream
100
0.0
Time-Related
3,650
20
Yes
Well No. 3
1,000
1,000
4,000
1,000
0.0015
180°
0.25
50
Barrier*
100
0.0
Time-Related
3,650
20
Yes
In practice, the direction of ambient ground-water flow will generally not be directly
toward a barrier boundary.
-------
Multiple Well Capture Zone Module (MWCAP)
7-7
Figure 7.2
Ten-Year Capture Zone For the Cases of No Boundary (a),
A Stream Boundary (b) and A Barrier Boundary (c)
(M)
4500
3600
(a)
2700
(b)
1800
900
(c)
(M)
900 1BOO 2700 3600 4500
-------
7-8 Multiple Well Capture Zone Module (MWCAP)
flow within the aquifer. The enclosed by the capture zone that intercepts the stream
is smaller than the other two because the stream boundary, in addition to the ambient
ground-water flow, acts as a source of water to the well.
7.4.3 Albuquerque Example
The third MWCAP example involves application of the code to an actual municipal
well field in New Mexico. The background information for this example was taken
primarily fromNewsom and Wilson (1988). Figure 7.3 presents the general hydrogeological
scenario; there are three closely spaced municipal supply wells located approximately half
a mile due east of the Rio Grande in the South Valley of Albuquerque, New Mexico.
The Albuquerque Basin lies within the Basin and Range physiographic province of the
United States. The basin fill (often referred to as "bolson fill") is primarily alluvial in nature
and consists of layered gravels, sands, silts and clays; it forms a very prolific aquifer. In the
vicinity of Albuquerque water-bearing sediments may have thicknesses on the order of
several thousands of feet. Leakage from the Rio Grande and mountain front recharge from
the Sandia Mountains to the east are the primary sources of recharge to the aquifer. The
climate is semi-and and therefore the assumption of negligible areal recharge due to rainfall
is a good one.
The direction of ground-water flow is from north to south (a = 270° or a = -90°) in
the ambient flow field west of the river which is unaffected by pumping. The sinuosity of
the Rio Grande was neglected and the river was approximated as a straight line boundary.
Because the three pumping wells are in close proximity to one another, and since MWCAP
does not account for well interference effects, the pumping from the three wells was
represented by an imaginary model well (marked M on Figure 7.3). The model well has a
discharge equal to the average sum of the discharges of the three municipal wells.
The steady-state capture zone delineated using MWCAP is presented in Figure 7.4.
Note that the y-axis in the figure corresponds to the stream (Rio Grande) boundary, arid
therefore a portion of the y-axis forms part of the capture zone boundary. The pathlines
that extend from the well to the stagnation point and from the well to the northern boundary
of the study area partition the recharge to the well; ground water to the left (west) of these
pathlines originates at the Rio Grande, while ground water to the right (east) of these
pathlines originates from ambient aquifer flow. All of the pathlines in Figure 7.4 were
generated automatically by MWCAP (i.e. the number of pathlines desired was set to zero).
If a non-zero number of pathlines had been specified, the additional pathlines would have
been started at the well bore and reverse tracked until they reached a study area boundary.
-------
Multiple Well Capture Zone Module (MWCAP)
7-9
Figure 7.3
General Ground-Water Flow System for Albuquerque MWCAP Example"
2 Miles
SCALE
Row Line
Water Level Contour,
Static Water Levels
Measured July and
August, 1983
Well Location, Wells
Completed Below
4600 Feet Above
Sea Level (After Kues, 1986)
P Real Pumping Well
O M Well Used for Model
From Newsom and Wilson, 1988. P's represent actual municipal wells, and M
represents the imaginary model well. The x and y coordinate axis used by MWCAP
are superimposed upon the general flow field.
-------
7-10 Multiple Well Capture Zone Module (MWCAP)
Figure 7.4
Steady-State Capture Zone for Albuquerque Municipal
Well Field Delineated Using MWCAP
(FT)
25000
20000
15000
10000
5000
STREAM CAPTURE
ZONE
WELL
N
STAGNATION POINT
_L
(FT)
5000
10000
15000
20000
25000
Simulation Options
Well No. 1
Units = ft and days
XMIN ~ 0.0; XMAX - 25000
YMIN = 0.0: YMAX = 25000
Step Length - 50
No. Pumping Wells = 1
X =
2,600
10,000
712,247
6,690
0.0013
T -
4-
I
a
S «= 0.25
b = 200
Stream Boundary
Distance to Boundary = 2,600
Boundary Orientation 0.0
Steady-State Capture Zone
No. of Patblines » 0
-------
Multiple Well Capture Zone Module (MWCAP) 711
The actual capture zone for the three Albuquerque wells probably differs somewhat
from that computed by MWCAP because the Rio Grande only partially penetrates the
aquifer (see Section 10.3 and Newsom and Wilson 1988). The depth of water in the Rio
Grande is probably no more than 5 to 10 ft at any given time, while the wells are pumping
from depths of up to 200 ft below land surface. Therefore, the capture zone in Figure 7.4
is not conservative because the recharge to the wells from the river is overestimated. A
conservative approach to this problem would be to ignore the presence of the stream
boundary.
Because the three wells are closely spaced relative to the size (especially the width) of
the capture zone, the simulation approach of lumping the discharge of the three wells to one
equivalent model well causes no significant errors. This point is illustrated in Section 8.2.4
where the Albuquerque example was run again using the semi-analytical GPTRAC module.
7.4.4 Seattle Example
The Final MWCAP example concerns the Highline well field in Seattle, Washington.
The background information for this example was taken from Seattle Water Department
(1986). A general location map is presented in Figure 7.5, and a local site map, along with
the locations of two pumping wells (Riverton Heights and Boulevard Park) is presented in
Figure 7.6.
The Highline area is underlain by three hydraulically interconnected aquifers (shallow,
intermediate, and deep) that occur within a thick sequence of unconsolidated glacial and
interglacial deposits. The aquifers are composed of cobbly sand and gravel deposits, while
the interbedded low-permeability aquitards that separate the major aquifer units are
composed of silt and clay. A hydrostratigraphic cross section through the Highline area is
presented in Figure 7.7.
The shallow aquifer is unconfined and ranges between 50 and 100 ft in saturated
thickness. This aquifer is not used as a major source of ground-water supply, but it does
supply significant recharge to the aquifers at depth.
-------
7-12
Multiple Well Capture Zone Module (MWCAP)
Figure 7.5
General Location Map for Seattle's Highline Well Field Example
SEATTLE
HIGHLINE
WELL FIELD
AREA
PIERCE
COUNTY
KING COUNTY
\
SNOHOMISH COUNTY
Supply Lines
N
-------
Multiple Well Capture Zone Module (MWCAP)
Figure 7.6
Local Site Map With Intermediate Aquifer Well Locations
for Highline Well Field Example
Legend:
Exploration Location and Number
W-1 ^ Test/Production Wei
B-10 O Boring/Piezometer
«*0 Spring*
Shallow Aquifer Monitoring
Well Location and Water
Quality Sampfing Point
iA A'i Hydrastraflgraphic Cross Section
I I Locaten and Designation
s* Boundary of Intermediate
_^> Aquiler (Approximaia)
8000
-------
7.14 Multiple Well Capture Zone Module (MWCAP)
Figure 7.7
Hydrostratigraphic Cross Section B-B' Through Highline Well Field Area
B North
W-2
Boulevard Pa*
400
OW-3
Glacier
South B*
W-1
Crestview
-400
Legend:
W-3 Wen Number
Well Location
Screened Section
Piezometric Surface
^_ Shallow Aquifer
I -2- Intermediate Aquifer
n JZ_ Deep Aquifer
Vertical Exaggeration x 10
Horizontal Scale in Feet
0 2000
C
0 200
Vertical Scale in Feet
4000
400
March 1986
HART-CROWSWER & associates inc.
-------
Multiple Well Capture Zone Module (MWCAP) 7
The intermediate aquifer is the target of most municipal ground-water development
and is the aquifer of interest for this example. It is confined, with static water levels typically
reaching 100 ft above the top of the aquifer. Most of the confinement is provided by the
upper and lower clay units (see cross section) that bound the aquifer above, below, and
along the east and west margins of the uplands. This aquifer is generally about 100 ft thick.
The deep aquifer is also confined and averages about 100 ft in thickness. It consists
of coarse cobbly sands and gravel deposits with localized silt and clay interbeds. This aquifer
is targeted for more limited ground-water development.
Figure 7.8 is a contour map of piezometric head for the intermediate aquifer. This
map indicates that the aquifer is recharged within the central upland area by the overlying
shallow aquifer, and that the general direction of ground-water flow is from the recharge
area toward regional discharge boundaries such as the Green River and Puget Sound.
Because the two production wells are relatively far apart (over 1 mile), an initial
assumption of negligible well interference effects was made and the five-year hybrid capture
zone for each well was delineated using MWCAP. The input for this example is provided
in Table 7.3. Proposed pumping rates and average transmissivity values in the vicinity of
each well were taken from Seattle Water Department (1986). The direction and magnitude
of the ambient hydraulic gradient was taken from Figure 7.8. Each well fully penetrates the
intermediate aquifer.
The MWCAP results were scaled using the GRAF "Map Scale" option so they could
be directly overlaid onto the base map to form Figure 7.9. Note in Figure 7.9 that the
capture zones extend across the regional ground-water flow divide in the intermediate
aquifer. This was permitted because drawdown due to the production wells will most likely
reach the divide, and alter its configuration. Ground-water flow divides may sometimes be
treated as barrier boundaries if pumping wells are sufficiently distant so as not to affect the
ground-water flow field in the vicinity of the divide.
The capture zones computed by MWCAP are conservative estimates since two major
sources of recharge to the well, leakage through the overlying aquitard due to pumping and
decreased leakage through the lower aquitard due to pumping, were neglected. Also, a
uniform ambient ground-water flow in the aquifer from southwest to northeast was assumed
to exist over the entire study area. Obviously, once the capture zone boundary
-------
7-16
Multiple Well Capture Zone Module (MWCAP)
Figure 7.8
Contour Map of Piezometric Head for Highline Well Field Intermediate Aquifer
Ground-Water w_2 A
Divid. >o \
\
W-1
Legend:
Exploration Location and Number
TesVProdocflonWei
Springs
Boundary o( IntermeoTate
Aquifar (ApprorimalB)
8000
-------
Multiple Well Capture Zone Module (MWCAP) n_\ n
Table 7.3
MWCAP Input Parameters used for Highline Well Field Example
Simulation Options
Units = feet and days
XMIN=0 YMIN =-15,000
XMAX = 15,000 YMAX = 10,000
Step Length = 100
No. of Pumping Wells = 2
Input Parameter
X Coordinate
Y Coordinate
Discharge
Transmissiviry
Hydraulic Gradient
Angle of Ambient Flow
Porosity
Aquifer Thickness
Boundary Type
Capture Zone Type
Travel Time Value
No. of Pathlines
Well No. &
12,428
-222
596,748
44,573
0.00385
45°
0.25
100
None
Hybrid
1,825
0
Well No. 2^
11,556
5,170
346,499
20,888
0.00385
45"
0.25
100
None
Hybrid
1,825
0
Riverton Heights
Boulevard park
-------
7-18
Multiple Well Capture Zone Module fMWCAP)
Figure 7.9
Five-Year Hybrid Capture Zones Delineated Using MWCAP
For Riverton Heights and Boulevard Park Wells
Legend:
Exploration Location and Number
W-1 ^ Test/Production We«
B-10 O Boring/Piezometer
Springs
A A
A*
» Shallow Aqurter Monitoring
We* Location and Water
Quality Sampling Point
i Hydrostratigraphic Cross Section
f Locaton and Designation
Boundary of Intermediate
Aquifer (Approximate)
8000
-------
Multiple Well Capture Zone Module (MWCAP) 7_1 Q
was extended across the ground-water flow divide this assumption became invalid. A
numerical ground-water flow model would have to be used to quantitatively assess the effects
of the reversal in ambient ground-water flow direction on the capture zones.
The Theis equation may be used to quickly assess some of the simplifying assumptions
made for this example. Jacob's "large time" approximation to the Theis equation may be
written:
2.3Q 2.25Tt
s = log (7.1)
4ni r S
where s is drawdown from an initial static water level at radial distance r from a well
pumping at rate Q, t is time since the start of pumping, and T and S are the aquifer
transmissivity and storativity respectively. The large time approximation is applicable
because only flow conditions that are approaching steady-state are of concern.
The radial distance from each well to the ground-water divide, and the radial distance
from each well to the mid-point between the wells were computed from Figure 7.8. Using
equation (7.1) the drawdown at each of these locations at various times was calculated; the
results of this analysis are presented in Figure 7.10.
Figure 7.10 indicates that significant drawdown may be observed at each of the selected
radial distance values. This analysis would suggest that the assumption of negligible well
interference was violated. It would also suggest that drawdown due to pumping will affect
the ground-water flow divide and therefore the divide should not be treated as a barrier
boundary. However, the simple Theis analysis presented above neglects two important
sources of recharge to the wells: ambient aquifer flow, and leakage from the confining beds.
To address these physical processes in an appropriate manner, a more sophisticated solution
for drawdown (analytical or numerical) would have to be used and perhaps some field
observations would be collected.
-------
7-20
Multiple Well Capture Zone Module (MWCAP)
Figure 7.10
Evaluation of Well Interference Effects and Drawdown At
Ground-Water Divide Due to Each Well Using Theis Solution
Jacob's Approximation to Theis Solution:
2.3Q 2.25Tt
s = log -=
r2S
(7.1)
Well No. 1-- Riverton Heights
Well No. 2~ Boulevard Park
Q = 596,748 f^/day
T = 44,573 ft2/day
S = 4.92 x ID"4
Q = 346,499 f^/day
T = 20,888
S = 1.04 x 10'3
Drawdown due to each well at mid-point of line segment connecting the wells
(r= 2,916 ft):
r(ft)
2,916
2,916
2,916
Well No. 1
t(days)
180
365
1,825
s(ft)
8.91
9.66
11,37
r(ft)
2,916
2,916
2,916
Well No. 2
t(days)
180
265
1,825
s(ft)
9.05
9.98
12.1
Drawdown due to each well at ground-water divide directly upgradient from well
location.
r(ft)
2,600
2,600
2,600
Well No. 1
t(days)
180
365
1,825
s(ft)
9.15
9.90
11.6
r(ft)
4,600
4,600
4,600
Well No. 2
t(days)
180
265
1,825
s(ft)
7.85
8.78
10.9
-------
8.0 General Particle Tracking
Module (GPTRAC
8.1 Introduction
The GPTRAC module consists of two major components. The first component,
referred to as the "semi-analytical option", provides the option for pathline and time-related
capture zone delineation for simple cases using analytical velocity computation techniques
similar to those of RESSQC and MWCAP. This option assumes a homogeneous aquifer but
is able to accommodate a wider range of aquifer and boundary conditions than the RESSQC
and MWCAP modeling options.
The second component, referred to as the "numerical option", provides pathline and
time-related capture zone delineation for general eases using numerical velocity computation
techniques that require input of hydraulic head values at nodal points of a rectangular grid.
This option is designed to be used as a postprocessor for numerical ground-water flow
models. Additionally, it can be used as a processor for interpolated head values derived
from a piezometric map (field data). Once the hydraulic head values corresponding to the
nodal points of a rectangular grid are read into the GPTRAC code, time-related capture
zone analysis and general particle tracking can be performed.
8.2 Semi-Analytical Option
8.2.1 Capabilities
The semi-analytical option of GPTRAC is capable of delineating time-related capture
zones for a system of pumping and injection wells that fully penetrate a homogeneous
aquifer. A steady flow field is assumed, and the aquifer may be confined, unconfmed or
semi-confined (leaky). For each aquifer type, a stream or barrier boundary can be specified
along any edge of the study area. The effects of well interference are accounted for through
superposition of solutions resulting from individual wells. If the aquifer is confined, "strip"
aquifers (aquifers with two parallel boundaries) may be simulated. A maximum of 50
pumping wells may be used in GPTRAC. For cases involving confined and semi-confined
aquifers, injection wells may also be used. The maximum number of injection wells allowed
by the code is 20.
-------
8-2 General Particle Tracking Module (GPTRAC)
The number of patties reverse tracked from each pumping well may be defined
interactively by the user. In addition, particles can be released at any point within the
system to be subsequently forward or reverse tracked.
8.2.2 Assumptions and Limitations
Capture zones delineated using the semi-analytical option of GPTRAC are valid for
fully penetrating wells screened in aquifers that are essentially homogeneous. Ground-water
flow in the aquifers must be two-dimensional in an areal x-y plane (the Dupuit assumptions
are used for unconfined flow cases). A steady-state ground-water flow field is assumed. In
the case of a leaky aquifer system, aquitard leakage is assumed to be vertical.
If a stream or a barrier boundary is present, the boundary is assumed to be linear and
fully penetrating. The latter assumption is often violated in cases where stream boundaries
exist. The effect of a partially penetrating stream may be an important one and each
application should be examined on a site-by-site basis. In general, the greater the depth and
breadth of the stream in relation to the aquifer thickness, the more valid the fully
penetrating stream assumption. Also, stream boundary partial penetration effects decrease
as the distance from the stream to the well increases. The stream and the aquifer are
assumed to be in perfect hydraulic connection; the effects of a "clogging layer" between the
streambed and the aquifer are not considered.
If, in actuality, the stream is partially penetrating and/or there is a clogging layer of fine
grained material that lines the streambed, the capture zone obtained using GPTRAC will
be smaller than the "true" capture zone. The amount of error incurred will be dependent
upon the degree to which the above assumptions are violated.
8.2.3 Input Requirements
The input requirements for the GPTRAC semi-analytical option are outlined in
Table 8.1. The well-specific parameters must be entered for each well in the study, area.
If a stream or barrier boundary is present, it is assumed to correspond to one edge (top,
bottom or either side) of the study area. If the strip aquifer option is selected, the two
boundary conditions must correspond to two opposite, parallel sides of the study area (top
and bottom, or left and right).
-------
General Particle Tracking Module (GPTRAC)
8-3
Table 8.1
Input Requirements for GPTRAC Semi-Analytical Option
Program Variable
Description
For each problem
IAQFR:
IUNIT:
NPWELL:
NRWELL:
XMIN:
XMAX:
YMIN:
YMAX:
DLMAX:
TRANSM:
GRADNT:
ALPHA:
POROS:
B:
KPRIM:
BPRIM:
CAPN:
CAPH:
RMAX:
IBOUND:
NFPATH:
NRPATH:
TMSIM:
TMCAPZ:
Aquifer type (confined, semi-confined, unconfined)
Default units of input parameters (feet and days or meters and
days)
Number of pumping wells within the study area
Number of recharge (injection) wells within the study area
Minimum x-coordinate of study area (ft or m) '
Maximum x-coordinate of study area (ft or m)
Minimum y-coordinate of study area (ft or m)
Maximum y-coordinate of study area (ft or m)
Largest allowable step length, d£ (see section 4.1)
Transmissivity of the aquifer (ft2/d or m2/d)
Regional hydraulic gradient (ft/ft or m/m)
Angle of ambient ground-water flow (0-360') "
Aquifer porosity (dimensionless)
Aquifer saturated thickness (ft or m)
Confining bed hydraulic conductivity^ (ft/d or m/d)
Confining bed thickness^ (ft or m)
Areal recharge rate^ (ft/d or m/d)
Aquifer saturated thickness prior to pumping ^ (ft or m)
Maximum radius of influence of pumping well^ (ft or m)
Boundary condition type (no boundary, stream or barrier
boundary on one side of study area, or strip aquifer)
Number of forward-tracked pathlines
Number of reverse-tracked pathlines
Time period for which GPTRAC will be executed^ (days)
Time value for time-related capture zones** (days)
-------
8-4
General Particle Tracking Module (GPTRAC)
Table 8.1 (continued)
Input Requirements for GPTRAC Semi-Analytical Option
Program Variable
Description
For each pumping well (1=1, NPWELL) ~
XPWELL(I):
YPWELL(I):
QPWELL(I):
NSTLIN(I):
x-coordinate of well (ft or m)
y-coordinate of well (ft or m)
Well discharge rate^ (f^/d or m3/d)
Number of pathlines to be computed to delineate time
related capture zone (default = 20)
For each injection well (1=1, NRWELL)
XRWELL(I): x-coordinate of well (ft or m)
YRWELL(I): y-coordinate of well (ft or m)
QRWELL(I): Well recharge rate^ (f^/d or m3/d)
For each forward-tracked pathline (1=1, NFPATH) ~
FSTART(I,1): x starting coordinate (ft or m)
FSTART(I,2): y starting coordinate (ft or m)
For each reverse-tracked pathline (1=1, NRPATH) ~
RSTART(I,1): x starting coordinate (ft or m)
RSTART(I,2): y starting coordinate (ft or m)
a/
b/
Sf
Only required for semi-confined (leaky) aquifers.
Only required for unconfmed aquifers.
A simulation time (TMSIM) different from the capture zone time (TMCAPZ) may be
useful when arbitrary forward- or reverse-tracked pathlines are desired. These
pathlines will represent the distance that a particle traveled for a time period equal to
TMSIM.
The sign (+,-) of the well discharge or recharge rate does not need to be specified.
-------
General Particle Tracking Module (GPTRAC) £.5
8.2.4 Example Applications
To demonstrate some of the capabilities of the GPTRAC semi-analytical module, the
New Mexico and Seattle case studies solved previously using MWCAP were re-examined
using GPTRAC. The confined aquifer type was used for each of these examples.
Additional examples are presented in Section 8.4.
8.2.4.1 Albuquerque Example
The ground-water system at the Albuquerque site was described in Chapter 7. Figure
7.3 portrays the hydrogeological setting of the Albuquerque municipal wells, and the aquifer
input parameters were discussed in section 7.4.3. However, because GPTRAC accounts for
the effects of well interference, the three pumping wells were not lumped to one equivalent
(imaginary) model well, but rather were treated individually. The total discharge of 712,247
ftVd was divided equally between the three wells so that each well had a pumping rate of
237,416 fr/d. The same coordinate system depicted in Figure 7.3 was used so that the
results of both modules could be overlayed. A stream boundary was specified along the left
side of the study area.
The 25-year capture zones for the three wells computed using GPTRAC are shown in
Figure 8.1. Each well draws water from the Rio Grande, but the capture zones of the top
and bottom wells force the middle well to obtain much of its discharge from captured
ambient flow. If these results are indicative of actual site conditions, one would expect water
quality analysis from the northernmost and southernmost wells to reflect river water
characteristics, while water quality for the middle well should reflect a combination of river
water and the ambient aquifer flow.
As expected, one can see by overlaying Figures 8.1 and 7.4 that the 25-year capture
zones lie within the steady-state capture zone delineated using MWCAP. Therefore, if the
objective of the capture zone delineation exercise was to designate a WHPA for the entire
well field, and the source of water for each individual well was not of concern, the approach
used in section 7.4.3 of analyzing one "equivalent" model well would be valid.
The same limitations addressed in section 7.4.3 apply to this analysis. Although the
Rio Grande probably contributes a significant portion of recharge to the pumping wells, the
amount is certainly overestimated by GPTRAC due to the fully penetrating stream
assumption. It is highly likely that the actual capture zones of the three wells extend to the
west underneath the river.
-------
8-6
General Particle Tracking Module (GPTRAC)
Figure 8.1
Twenty-Five-Year Capture Zones for Three Albuquerque
Municipal Wells Located Near the Rio Grande
(FT)
25000
20000
-< Rio Grande
15000
10000
5000
_L
J_
(FT)
5000 10000
15000 20000 25000
Simulation
Options
Units = ft and days
Aquifer Type = confined
Step Length - 200 ;
No. Pumping Wells = 3
Simulation Time - 9,125
Capture Zone Time = 9,125
Stream Boundary on Left
Side of Domain
Aquifer
Properties
T
fc
6
i
a
= 6,690
= 200
= 0,25
= Q.0013
= -90*
Well #1
X = 2,000
Y = 8,900
Q * 237,416
No. Pathlines = 10
Well #2
3,500
10,000
237,416
10
Well #3
2,100
11,000
237,416
10;
-------
General Particle Tracking Module (GPTRAC)
8.2.4.2 Seattle Example
The Seattle Highline Well Field example was also re-examined using GPTRAC. When
this example was presented in section 7.4 as an MWCAP problem, it was noted that a simple
Theis analysis indicated that well interference effects between the two pumping wells were
probably not negligible (i.e. the cones of depression for each well overlapped one another).
Because the semi-analytical GPTRAC module (unlike MWCAP) accounts for well
interference effects, the Seattle example was reanalyzed and differences in the capture zones
were compared.
Background information for the Highline well field is presented in section 7.4.4. Recall
that although the aquifer thickness remains fairly constant, transmissivity values reported in
the vicinity of the two wells are markedly different. Because the semi-analytical GPTRAC
module assumes a homogeneous aquifer, the average transmissivities for the Riverton
Heights area (T = 44,573 ftVd) and the Boulevard Park area (T = 20,888 ftVd) were
averaged to obtain an overall transmissivity of 32,730 ftVd. Although this average value of
transmissivity is sufficient for this example, several additional approaches would be worthy
of consideration if the actual WHPAs were to be designated. For example, a sensitivity
analysis would be in order where, in two separate model runs, the overall aquifer
transmissivity is set equal to the transmissivity values measured in the vicinity of each well.
The capture zone results obtained using MWCAP and GPTRAC are superimposed in
Figure 8.2. A transmissivity of 32,730 ft Yd, 10 pathlines for each well, the confined aquifer
type, and the additional parameters presented in Table 7.3 were used to obtain the
GPTRAC capture zones. The plot was not scaled to overlay on the base map for this
example so that the results of each module could be more easily compared.
The effects of well interference are evident in Figure 8.2 in that the GPTRAC capture
zones are not perfectly aligned with the angle of ambient ground-water flow (450), and there
also exists a slight curvature in the capture zone outlines (particularly the top well).
The time-related capture zones and the hybrid capture zones do not have compatible
lengths despite the fact that they both represent time periods of five years. This is because
different transmissivity values were used for each well in MWCAP, but a single
-------
8-8
Lrenerai Farticle i racking ivioauie (LTFIKAL)
Figure 8.2
Five-Year Hybrid Capture Zones (MWCAP) and Five-Year Time-Related
Capture Zones (GPTRAC) Computed for Highline Well Field
(FT)
10000
5000
-5000
-10000
-15000
(FT)
0 3000 6000 9000 12000 15000
-------
General Particle Tracking Module (GPTRAC)
"averaged" value was required by GPTRAC. This approach leads to an interesting
comparison of methods.
For the top well, the GPTRAC "average" transmissivity is higher than the MWCAP
value. Therefore the 5-year GPTRAC capture zone is narrower and longer than the
MWCAP 5-year hybrid capture zone. For this well, the MWCAP capture zone is more
consemative widthwise, but less conservative than the GPTRAC capture zone lengthwise.
For the bottom well the opposite is true. The "average" GPTRAC transmissivity is
lower for this well than was the MWCAP value. Therefore, the GPTRAC capture zone is
more conservative than the MWCAP capture zone widthwise, but less conservative
lengthwise.
83 Numerical Option
The numerical option of GPTRAC requires hydraulic heads at the nodes of a
rectangular mesh as input. If head values are observed in the field or read from a map, they
may be interpolated onto the nodal points of a grid using methods such as linear
interpolation or kriging. More commonly, however, the head values supplied to GPTRAC
will be the output of a finite difference or finite element ground-water flow model. In
addition to nodal heads, GPTRAC requires the aquifer transmissivity (T), porosity (o), and
thickness (b). If the aquifer is heterogeneous, it may be divided into multiple zones of
varying T, 6 and/or b, and if it is anisotropic, each zone may have directional transmissivities
TX and Ty.
GPTRAC uses the above information, in conjunction with linear finite element or finite
difference approximations, to determine the x and y velocity components of ground-water
flow at the edges of each rectangular element or grid block. Interlock or interelement
hydraulic conductivities are computed by taking harmonic averages of the block or element
hydraulic conductivities. If a finite element grid is used, the velocity components are also
computed at the element centroids for the elements that share the well nodes. Pathlines of
individual particles are delineated using particle tracking techniques based on semi-analytical
and numerical integration. The computational procedure is described in detail in
Appendix C.
-------
8-10 General Particle Tracking Module (GPTRAC)
8.3.2 Capabilities
The numerical option of GPTRAC can be used to delineate time-related capture
zones. If required, multiple pathlines starting at prescribed locations within the system can
also be delineated. The number of pathlines used to delineate capture zones may be
specified interactively by the user. Pathlines may be delineated using either forward or
reverse particle tracking.
8.3.2 Assumptions and Limitations
Because the numerical GPTRAC option is based upon the availability of an observed
or model calculated head field, the assumptions and limitations associated with it are
substantially less restrictive than those associated with the other WHPA model options. The
first assumption is that the ground-water flow field is at equilibrium (steady state). The
second limitation is that flow in the aquifer must be two-dimensional in the horizontal plane;
vertical flow components are neglected. Since GPTRAC does not obtain flow velocities
from an analytical solutiom the aquifer need not be homogeneous.
83.3 Input Requirements
The input requirements for the GPTRAC numerical option are outlined in Table 8.2.
Note that many of the variables may require iterative input depending upon the number of
pumping wells, recharge wells, aquifer material zones, forward-tracked pathlines, and
reverse-tracked pathlines.
The most substantial input requirement for the GPTRAC numerical option that is
different from the other WHPA model modules is the hydraulic head file. The structure of
this file is illustrated in Figure 8.3. The head file should be a standard ASCII file with the
format 5(I5,F10.3). The first five spaces of each input field contain a node number, and the
next ten spaces contain the head value associated with that node (head values may have a
precision up to three decimal points). This pattern is repeated five times across each line
of the input file. The utility program HEDCON, documented in Appendix E, is provided
to assist users with the construction of GPTRAC head files.
The node number associated with each hydraulic head value is not used explicitly by
the code, and therefore a dummy variable (e.g. set all node numbers to one) may be used
if desired. However, the scheme of associating a node number with each head value
provides a convenient method to assist users in visualizing the proper structure of the
hydraulic head input file (Figure 8.3).
-------
General Particle Tracking Module (GPTRAC)
Table 8.2
Input Requirements for GPTRAC Numerical Option
Program Variable
Description
For each problem--
IUNIT:
NPWELL:
NRWELL
XMIN:
XMAX:
YMIN:
YMAX
NZONES:
NROWS and NCOLS,
and XGRIDL(I),
YGRIDL(J) for
1=1, NCOLS and
J=l, NROWS
FILE1:
FILE2:
NFPATH:
NRPATH:
TMSIM:
TMCAPZ:
Default units of input parameters (feet and days or meters
and days)
Number of pumping wells within the study area
Number of recharge (injection) wells within the study area
Minimum x-coordinate of study area (ft or m)
Maximum x-coordinate of study area (ft or m)
Minimum y-coordinate of study area (ft or m)
Maximum y-coordinate of study area (ft or m)
Number of aquifer zones that have different" material
properties (if aquifer is nonuniform)
Number of grid-line rows and columns, and (if non-uniform
grid,) coordinates of each grid-line row and column
Input file name that contains head values for each node of
the finite element or finite difference grid.
Input file name that contains the number of grid-line
columns and rows and their coordinate values if file input
option is selected.
Number of forward-tracked pathlines
Number of reverse-tracked pathlines
Time period for which GPTRAC will be executed3' (days)
Time value assigned to time-related capture zones- (days)
For each pumping well (1=1, NPWELL) ~
XPWELL(I):
YPWELL(I):
QPWELL(I):
NSTLIN(I):
x-coordinate of well (ft or m)
y-coordinate of well (ft or m)
Well discharge rate^ (ft3/d or m3/d)
Number of pathlines to be computed to delineate time-
related capture zone (default = 20)
-------
8-12 General Particle Tracking Module (GPTRAC)
Table 8.2 (continued)
Input Requirements for GPTRAC Numerical Option
Program Variable
Description.
For each injection well (1=1, NRWELL) ~
XRWELL(I): x-coordinate of well (ft or m)
YRWELL(I): y-coordinate of well (ft or m)
QRWELL(I): Well recharge rate^ (f^/d or m3/d)
For each aquifer property zone (1=1, NZONES)^ ~
XMINZO(I):
XMAXZO(I):
YMINZO(I):
YMAXZO(I):
TRANSX:
TRANSY:
BZO(I):
PORZO(I):
Minimum x-coordinate of zone (ft or m)'
Maximum x-coordinate of zone (ft or m)
Minimum y-coordinate of zone (ft or m)
Maximum y-coordinate of zone (ft or m)
Transmissivity of zone - x direction^ (fr/d or m /d)
Transmissivity of zone - y direction (fr/d or m /d)
Saturated thickness of aquifer in zone (ft or m)
Porosity of aquifer in zone (dimensionless)
For each forward tracked pathline (1=1, NFPATH) ~
FSTART(I,1): x starting coordinate (ft or m)
FSTART(I,2): y starting coordinate (ft or m)
For each reverse tracked pathline (1=1, NRPATH) ~
RSTART(I,1): x starting coordinate (ft or m)
RSTART(I,2): "y starting coordinate (ft or m)
a/
£/
A simulation time (TMSIM) different from the capture zone time (TMCAPZ) may be useful
when arbitrary forward- or reverse-tracked pathlines are desired. These pathlines will represent
the distance that a particle traveled for a time period equal to TMSIM.
The sign (+,-) of the well discharge or recharge rate does not need to be specified.
The first zonal area specified should be the entire study area. Therefore, it is most convenient
to consider the largest zone to be zone 1, and then overlay the other zones on top of it. Zones
of inactive cells should be given zero values for transmissivity.
The calculations in GPTRAC use the hydraulic conductivity of each material zone. Therefore,
the transmissivity of each zone is divided by the thickness to obtain the conductivities,
CONXZO(I) and CONYZO(I) for each zone I.
-------
General Particle Tracking Module (GPTRAC)
Figure 8.3
Schematic Representation of Head Data File Format Required by GPTRAC
For Finite Element or Mesh-Centered Finite Difference Model Output With
Nodes Numbered in the y-Direction (a), or for Block-Centered Finite Difference Model
Output With Nodes Numbered in the x-Direction (b)
(4) 905.0
(3)
(2)
(1)
1
6
11
=16
882.0
880.0
876.0
877.0
(8) 895.0
(12) 885.0
(16) 877.0
894.0 (7)
887.0 (6)
882,0 (5)
888.0 (11)
880.0 (10)
876.0 (9)
876.0 (15)
872.0 (14)
868.0 (13)
871.0
865.0
862.0
(1) = node number
2 887.0 3 894.0 4 905.0 5 876.0
7 888.0 8 895.0 9 868.0 10 872.0
12 885.0 13 862.0 14 865.0 15 871.0
(a)
905.0
894.0
e
887.0
882.0
895.0
888.0
e
880. Q
876.0
885.0
876.0
s
872.0
868.0
877.0
871.0
«
865.0
862.0
864.0
861.0
e
860.0
858.0
1
2
1
6
11
16
882.0
887.0
894.0
905.0
2
7
12
17
876.0
880.0
888.0
895.0
3
8
13
19
868.0
872.0
876.0
885.0
4
9
14
19
862.0
865.0
871.0
877.0
5
10
15
20
858.0
860.0
861.0
864.0
-------
8-14 General Particle Tracking Module (GPTRAC)
Note that on the grid specification input screen, the user is prompted for the number
of grid-line columns and the number of grid-line rows. The number of grid-line columns is
equal to the number of columns in the grid plus one, and the number of grid-line rows is
equal to the number of rows in the grid plus one. The block-centered finite difference grid
in Figure 8.3b has 5 columns and 4 rows, but it has 6 grid-line columns and 5 grid-line rows.
Many finite difference codes (e.g., MODFLOW, McDonald and Harbaugh, 1988) use
formulations based upon a coordinate system where the y-axis is opposite that of a
conventional Cartesian system. These codes assume an origin in the upper left-hand comer
of the study area instead of the' conventional lower left-hand comer location. Some
postprocessing is required to rearrange the output from these codes to the format suitable
for WHPA.
If the MODFLOW code is used, the POSTMOD program supplied with the code (if
it was obtained from the International Ground Water Modeling Center) may be used to
process the binary hydraulic head output file. POSTMOD has the capability to rearrange
MODFLOW output to conform with a standard coordinate system where y-coordinates are
measured from the lower left-hand comer rather than the upper left-hand comer of the
modeled domain. The head file created using POSTMOD will have a x-coordinate, a y-
coordinate and a head value on each line of the file. The HEDCON program supplied with
the WHPA code may be used to convert the POSTMOD output file to a proper input
format for GPTRAC.
The HEDCON program may also be used to construct GPTRAC input head files from
data files that have nodal x- and y-coordinates and head values on each line, but which might
be organized in no consistent fashion. Refer to Appendix E for detailed instruction on the
capabilities and proper use of HEDCON.
83.4 Example Applications
Three examples that demonstrate the capabilities of the GPTRAC numerical option
are presented in this section. The first two examples are hypothetical; the third example re-
examines the Coming surficial aquifer previously analyzed using RESSQC.
-------
General Particle Tracking Module (GPTRAC)
8.3.4.1 Hypothetical Examples
The hypothetical aquifer used for the first two examples is shown in Figure 8.4. There
are three distinct material zones within the aquifer. Zone 1 is isotropic, and zones 2 and 3
are anisotropic (transmissivity is directionally dependent). Ambient ground-water flow is
from north to south. The left and right sides of the aquifer are no-flow (barrier) boundaries.
Two pumping wells discharging at 4,000 mVd and 3,000 mVd are screened over the entire
aquifer depth of 20 m.
A hydraulic head field for the described hydrogeological scenario was obtained using
a two-dimensional finite element ground-water flow and transport code called SAFTMOD
(Huyakom and Buckley, 1988). The steady-state head field computed by SAFTMOD, along
with the aquifer geometry and the zonal hydraulic properties, were input to GPTRAC and
the 100-day capture zone for each pumping well was delineated. The input parameters for
this example are presented in Table 8.3, and the capture zones are illustrated in Figure 8.5
The capture zones in Figure 8.5 have a complex shape and could not be delineated
using the semi-analytical module. The top well located in zone 1 preferentially draws water
from the region of high transmissivity, this is evidenced by the "clustering" of pathlines that
bend around the lower right-hand comer of zone 3. In general, pathlines that enter low per-
meability materials will exhibit increased spacing between them as opposed to pathlines
located in highly permeable materials. Note also that the pathlines refract at the material
property interfaces. The capture zone for the top well is slightly larger than the capture
zone of the bottom well primarily due to its larger pumping rate.
The results of the second example are presented in Figure 8.6. This example is
identical to the first example, except that the pumping rate of the second well has been
increased from 3,000 mVd to 5,000 mVd and its location was moved from zone 2 into zone
1 at x= 300 m andy = 333 m (Table 8.3). Again, it is readily apparent from Figure 8.6
that both wells tend to draw water from (and therefore have capture zones within) the zone
of highest transmissivity.
For the second example, two initial locations from which particles were forward- and
reverse-tracked for a period of 300 days were specified. The pathlines of these particles are
indicated in Figure 8.6. The forward-tracked particle entered the first pumping well at some
time less than 300 days. The reverse-tracked particle travelled quite slowly during the 300
day period (as indicated by its short pathline) because it was released and tracked in the
vicinity of the stagnation point formed by the bottom pumping well. If desired, additional
particles could have been released at arbitrary locations within the aquifer.
-------
8-16 General Particle Tracking Module (GPTRAC)
Figure 8.4
Hypothetical Aquifer and Hydraulic Properties Used
for First Two GPTRAC Numerical Examples
600 m
300 m
1000 m
h = 0
550 m
I
ZONE 3
ZONE 1
4000 m3/d
(Xw,Yw) = (233.3,500)
j^3000 m3/d
(Xw.Yw) = (500,200)
ZONE 2
300 m
h = -10 m
' 1000 m
'
Zone No. Tx '. Ty b
0
1 2,000 2,000 20 0.25
2 1,000 200 20 0.20
3 500 200 20 0.15
-------
General Particle Tracking Module (GPTRAC)
Table 8.3
Input Parameters for GPTRAC Hypothetical Aquifer
(Examples One and Two)
EXAMPLE ONE
GPTRAC Numerical Option Using Rectangular Finite Element Model
units = meters and days
Uniform Grid (Constant Spacings) - Yes
Uniform Aquifer Properties - No
Hydraulic Head Data File - TP8EX1
XMIN = 0.0 YMIN = 0.0
XMAX = 1000 YMAX = 1000
No. of Grid-line Rows =31
No. of Grid-line Columns = 31
Nodes are numbered along y - axis (option O)
No. of Pumping Wells = 2
No. of Injection Wells = O
No. of Material Zones = 3
Zone No.
1
2
0
XMIN
0
300
0
XMAX
1,000
1,000
550
YMIN
0
0
600
YMAX
1,000
300
1,000
Tx
2,000
1,000
500
Ty
2,000
200
200
b
20
20
20
8
0.25
0.20
0.15
Simulation Time =300 days
Capture Zone Time =100 days
Well No.
1
2
X
233
500
Y
500
200
Q
4,000
3,000
Capture Zone
Yes
Yes
No. of Pathlines
15
15
No. of Forward Pathlines = O
No. of Reverse Pathlines = O
-------
8-18 General Particle Tracking Module (GPTRAC)
Table 8.3 (continued)
Input Parameters for GPTRAC Hypothetical Aquifer
(Examples One and Two)
EXAMPLE TWO
Same as example 1 except for the following changes:
Hydraulic Head Data File - TP8EX2
Well No. X
1
2
233
300
Y
500
333
Q Capture Zone
4,000
5,000
Yes
Yes
No. of Pathlines
20
20
No. of Forward-Tracked Pathlines = 1
Starting Location: X = 450, Y = 800
No. of Reverse-Tracked Pathlines = 1
Starting Location: X = 175, Y = 175
-------
General Particle Tracking Module (GPTRAC) 8-19
Figure 8.5
One-Hundred-Day Capture Zones for Two Pumping Wells in Hypothetical Aquifer
(M)
1000
800 -
600
400 -
200 -
(M)
200 400 600 800 1000
-------
8-2O General Particle Tracking Module (GPTRAC)
Figure 8.6
One-hundred-Day Capture Zones for Two Wells and Three-Hundred-Day Reverse-
arid Forward-Tracked Pathlines for Hypothetical Aquifer
(M)
1000
800
600
400
200
2DNE3
FORWARD-TRACKED
FATHUNE
ZDNE1
REVERSE-TRACKED
PATHUNE
ZDNE2
I , I
(M)
200
400
600
800
1000
-------
General Particle Tracking Module (GPTRAC)
8.3.4.2 Corning Example
The final example is based on the results of a numerical ground-water flow simulation
for the surficial aquifer in the vicinity of Corning, New York. Some background information
for the Corning aquifer was provided in Section 6.5.2. The study area (Figure 8.7) is a
subdomain of the region studied by Ballaron (1988).
The valley sediments in the vicinity of Corning consist of stratified glacial drift deposits
that are primarily interbedded silty to clean sands and gravels. Relatively thin deposits of
lacustrine clay, silt and fine sand exist over much of the valley and separate a surficial,
unconfined aquifer from a confined to semi-cordined aquifer at depth. Two idealized cross
sections through the study area are presented in Figure 6.4.
The surficial aquifer has an average saturated thickness of about 25 ft. No-flow
boundaries were used where the surficial aquifer sediments pinch out or abut against the
older, low permeability valley walls. Values for the prescribed hydraulic head boundary
nodes (Figure 8.8) were interpolated from the steady-state head map for the surficial aquifer
presented in Ballaron (1988).
The approximate distribution of hydraulic conductivity for the surficial aquifer was also
taken from Ballaron and is presented in Figure 8.8. An areal recharge rate of 0.00297 ft/d
(13 in/yr) was used for the entire study area. Recharge from the Chemung River and
leakage between the surficial and lower aquifer units were neglected. Ground-water flow
in the aquifer was simulated using the USGS block-centered finite difference code
MODFLOW (McDonald and Harbaugh, 1988). The grid consisted of 36 rows and 42
columns (37 grid-line rows and 42 grid-line columns). Each cell was 250 ft on a side.
The results of the first simulation for steady-state conditions are illustrated in
Figure 8.9a. This is the steady-state head map that was used in the RESSQC examples
section to obtain the hydraulic gradient and direction of ambient flow information for the
Corning example.
For the second simulation, three pumping wells were assumed to fully penetrate the
surficial aquifer. The location and rates of pumping for these wells, along with the steady-
state head field predicted using MODFLOW are illustrated in Figure 8.9b.
-------
8-22 General Particle Tracking Module (GPTRAC)
Figure 8.7
General Site Map and Modeling Area Boundary for Surficial
Aquifer in the Vicinity of Corning, NY
70o Corning, New York
1 875 1,750 FEET
-------
General Particle Tracking Module (GPTRAC) 8-23
Figure 8.8
Finite Difference Model Boundary Conditions and Distribution of
Hydraulic Conductivity for Surficial Corning Aquifer
1 '
. ks1 FT/O I
ks30FT/D
1
k a 150 FT/D
k=400FT/D
Prescribed Head Boundary
No-Row Boundary
Hydraulic Conductivity interface
Pumping Well
1,750ft
Scale
-------
8-24 General Particle Tracking Module (GPTRAC)
Figure 8.9
Steady-State Head Field for Surficial Coming Aquifer Using MODFLOW
for-(a) No Pumping Wells and (b) Three Pumping Wells
-------
General Particle Tracking Module (GPTRAC) £-25
Once the head field for the second simulation was obtained using MODFLOW, the
WHPA model was applied to delineate the 5-year capture zones for the three pumping
wells. The required input for the GPTRAC numerical option is provided in Table 8.4, and
the delineated capture zones are presented in Figure 8.10. The capture zones for this
example are very complex due to the aquifer heterogeneities and well interference effects.
Note that since MODFLOW uses a coordinate system based upon an origin in the
upper left-hand comer of the study area, some manipulation of the output head file was
performed to obtain the head values ordered in a suitable fashion for input into WHPA.
The standard POSTMOD program supplied with the MODFLOW code (when it is obtained
from the International Ground Water Modeling Center) was used to process the binary head
file output by MODFLOW. The POSTMOD code has an option to produce an ASCII file
in which the head values are ordered according to a standard coordinate system (i.e. origin
in lower left-hand comer of study region). This file may then be used directly by numerous
software plotting packages such as SURFER (Figure 8.9 was constructed using the ASCII
file and SURFER). This ASCII file was further rearranged using the program HEDCON
supplied with the WHPA code. One HEDCON program option is to create an input head
file suitable for WHPA using an ASCII file output by POSTMOD.
By comparing the capture zones delineated using GPTRAC and RESSQC (Figures
8.10 and 6.6), the effects of areal recharge and aquifer heterogeneities on the Coming
capture zones may be qualitatively assessed. The capture zones for wells 1 and 2 in the
center of the study area are markedly similar. The RESSQC capture zones are larger in
areal extent because recharge from local precipitation was neglected.
The capture zones for the third well are quite different due to aquifer heterogeneities.
In Figure 8.10 this capture zone extends from the well almost due west; in Figure 6.6 this
capture zone is aligned more with the angle of ambient flow (northwest). This effect is do
primarily to the fact that well 3 is located within a zone of high aquifer permeability and its
effects on the capture zone could not be represented using RESSQC. For well 3, neglecting
aquifer heterogeneities has very significant consequences.
-------
8.26 General Particle Tracking Module (GPTRAC)
Table 8.4
Input Parameters for GPTRAC Numerical Option for Corning Example Problem
Block-Centered Finite Difference Option
units ' feet and days
Uniform Grid (Constant Spacings) - Yes
Uniform Aquifer Properties - No
Hydraulic Head Data File - CRNMOD1.HED
XMIN =0.0
XMAX = 10500
YMIN =0.0
YMAX = 9000
No. of Grid-line Rows = 37
No. of Grid-line Columns = 43
Nodes are numbered along.x - axis (option 1)
No. of Pumping Wells = 3
No. of Injection Wells = 0
No. of Material Zones = 20
Zone No.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
XMIN
0
0
0
1,750
1,750
3,250
4,250
4,750
4,250
1,750
3,750
5,250
6,750
8,250
8,250
8,750
9,750
9,750
1,750
2,750
XMAX
10,500'
1,750
1,750
3,250
3,250
4,250
4,750
5,250
7,250
3,750
5,250
6,750
8,250
10,500
8,750
9,750
10,500
10,200
8,250
3,250
YMIN
0
5,250
2,750
6,250
7,750
6,250
6,750
6,750
6,250
2,250
1,750
1,250
750
250
2,750
2,750
2,750
4,000
3,250
8,750
YMAX Tx
9,000 0
8,250 1,080
5,250 780
7,750 1,440
8,750 48
9,000 1,080
9,000 1,020
7,250 1,020
6,750 960
3,250 11,200
3,250 11,200
3,250 10,400
3,250 9,600
2,750 6,400
5,750 3,300
5,250 3,000
4,000 2,400
4,250 2,400
6,250 2,400
9,000 48
iy
0
1,080
780
1,440
48
1,080
1,020
1,020
960
11,200
11,200
10,400
9,600
6,400
3,300
3,000
2,400
2,400
2,400
48
b
1
36
26
48
48
36
34
34
32
28
28
26
24
16
22
20
16
16
30
48
6
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
0.22
, 0.22
Simulation Time = 1825 days
Capture Zone Time = 1825 days
-------
General Particle Tracking Module (GPTRAC) £-27
Table 8.4 (continued)
Input Parameters for GPTRAC Numerical Option for Corning Example Problem
Well No.
1
2
3
X
4,375
6,375
7,875
Y
4,875
4,375
2,375
Q Capture Zone
25,000
30,000
30,000
Yes
Yes
Yes
No. of Pathlines
10
10
10
No. of Forward Pathlines = 0
No. of Reverse Pathlines = 0
-------
8-28 General Particle Tracking Module (GTPRAQ
Figure 8.10
Five-Year Capture Zones of Three Pumping Wells in the Vicinity of Corning, NY
(FT)
9000
7200 -
5400 -
3600 -
1800 -
(FT)
0 2100 4200 6300 8400 10500
-------
General Particle Tracking Module (GPTRAC)
8.4 Comparison of Semi-Analytical and Numerical Options
In this section, - both the semi-analytical and numerical options of GPTRAC are
exercised with the intent of providing test examples that verify some additional analytical
solutions available in the code for various aquifer and boundary conditions. In each
example, capture zones obtained using the two modeling options are compared.
8.4.1 Strip Confined Aquifer Examples
This example involves strip confined aquifers bounded on two (parallel) sides. Three
cases with different combinations of barrier and stream boundaries are considered. For each
case, the aquifer transmissivity and thickness are 1,000 ft Yd and 50 ft, respectively, and the
effective porosity is 0.25. There is one well located at x=2,000 ft and y= 1,200 ft, and the
pumping rate of the well is 3,000 ft Yd. Ambient ground-water flow is assumed to be zero.
Figures S.lla through 8.lie show capture zones computed using the semi-analytical module
of GPTRAC for a very large time value, t= 106days (to ensure a steady-state condition).
Note that the three cases correspond to: (1) aquifer with two barrier boundaries, (2) aquifer
with two stream boundaries, and (3) aquifer with a stream and a barrier boundary. In each
case, the boundaries are parallel to the x-axis at y=0 and at y= 1,500 ft. The three cases
were also simulated using the numerical option of the code. For each case, a uniform finite
element grid with AX=AY= 100 ft was used, and zero drawdown conditions were imposed on
the left and right hand side boundaries of the flow domain. Figure 8.12 illustrates the
resulting capture zones. A comparison of the two sets of Figures (8.11 and 8.12) shows good
agreement between the semi-analytical and the numerical solutions for each case.
8.4.2 Unconfined Aquifer Examples
Four unconfmed flow examples with different values of pumping rate, recharge and
hydraulic conductivity are presented in this section. In each case, there is one pumping well
and the aquifer is assumed to be homogeneous and isotropic. The effective porosity of the
aquifer is 0.25.
For the first case, an unconfined aquifer with a well pumping near a stream is
considered. The well is located at x=350 ft, and y= 1,250 ft. The pumping rate is
5,000 ftYd. The hydraulic conductivity and initial saturated thickness of the aquifer are 20
ft/d and 48.12 ft, respectively. Ambient ground-water flow and area! recharge are assumed
to be zero. Figures 8.13a and 8.13b are the 20,000-day capture zones computed using the
-------
8-30 General Particle Tracking Module (GPTRAC)
Figure 8.11
Strip Aquifer Simulations Using the Semi-Analytical GPTRAC Module for Two
Barriers (a), Two Streams (b), and a Stream and a Barrier (c)
(FT)
1500
(FT)
1500
(FT)
1500
(FT)
BOO 1600 2400 3200 '4000
(a)
(FT)
800 1600 2400 3200 4000
(b)
(FT)
BOO 1600 2400 3200 4000
(c)
-------
General Particle Tracking Module (GPTRAC) g-31
Figure 8.12
Strip Aquifer Simulations Using the Numerical GPTRAC Module for Two
Barriers (a), Two Streams (b), and a Stream and a Barrier (c)
(FT)
1500
0
0
800
(FT)
1600 2400 3200 4000
(FT)
1500
0
(FT)
1500
_L
(a)
(FT)
800 1600 2400 3200 4000
(b)
(FT)
800 1600 2400 3200 4000
-------
8.32
General Particle Tracking Module (GPTRAC)
Figure 8.13
Capture Zones Computed Using the Semi-Analytical (a) and Numerical (b)
GPTRAC Options for Unconfined Aquifer with Zero Recharge
(FT)
2500
1666
833
(FT)
2500
1666 -
833 =.
(FT)
833 1666 2500
(a)
(FT)
2500
-------
General Particle Tracking Module (GPTRAC) £.33
semi-analytical and numerical options, respectively. For the semi-analytical option, a stream
boundary condition with a constant head value of 48.12 ft was imposed along the y-axis of
the plot. For the numerical option, a uniform rectangular grid with nodal spacings
AX=Ay=50 ft was used. Nodal head values were obtained using a finite element flow model
with a constant head boundary condition of h = 48.12 imposed on each side of the rectangular
domain. A comparison between Figures 8.13a and 8.13b shows an excellent agreement
between the capture zones delineated using the semi-analytical and numerical GPTRAC
options.
The second case concerns an unconfmed aquifer with areal recharge of 0.0023 ft/d
(10.0 in/yr), and a well pumping rate of 20,000 fr/d. The remaining data and boundary
conditions are the same as those for the previous case. Figures 8.14a and 8.14b illustrate
the 2,000-day capture zones obtained using the semi-analytical and numerical options
respectively. Again, the two solutions agree quite well.
For the third and fourth cases, the semi-analytical option of GPTRAC was used to
model a well pumping in an areally infinite unconfined aquifer with recharge equal to 0.0023
ft/d and no ambient flow. The two cases were designed to demonstrate the effect of the
reduction of hydraulic conductivity and saturated thickness of the aquifer on the extent of
the capture zone. Figures 8.15a and 8.15b show the 4,000-day capture zones predicted by
the GPTRAC semi-analytical model for hydraulic conductivity values of 20 and 2 ft/d,
respectively. Both capture zones correspond to the well pumping rate of 5,000 fr/d. Note
that the two capture zones are circular and the low-conductivity capture zone (K=2 ft/d) has
a slightly higher radius. The reason for this is that the decrease of hydraulic conductivity
from 20 ft/d to 2 ft/d led to a reduction in saturated thickness (or increase in drawdown) in
the inner portion of the flow region surrounding the well. This resulted in higher velocities
for the same pumping rate.
-------
8-34 General Particle Tracking Module (GPTRAC)
Figure 8.14
Capture Zones Computed Using the Semi-Analytical (a) and Numerical (b)
GPTRAC Options for an Unconfined Aquifer with Areal Recharge
(FT)
2500
1666
833
(FT)
2500
1666
833
(FT)
833 1666 2500
(a)
(FT)
833 1666 2500
(b)
-------
General Particle Tracking Module (GPTRAC)
Figure 8.15
Capture Zones Computed Using Semi-Analytical GPTRAC Module for an Unconfined
Aquifer with Zero Recharge and Ambient Flow and a Hydraulic Conductivity
of(a)20ft/d,and(b)lft/d
(FT)
2500
1666
833
(FT)
2500
1666
833
I
(FT)
833 1666 2500
(a)
(FT)
833 1666 2500
(b)
-------
8.36 General Particle Tracking Module (GPTRAC)
8.43 Leaky Aquifer Example
In this example, a semi-confined aquifer underlain by an impermeable (barrier
boundary) and overlain by a semi-permeable aquitard layer is considered. When such an
aquifer is pumped, vertical leakage occurs through the aquitard. The top of the aquitard
is assumed to be subject to a zero drawdown condition. For the case simulated, the
hydraulic conductivity and thickness of the aquifer are 20 ft/d and 50 ft, respectively, and the
hydraulic conductivity and the thickness of the confining bed are 0.2 ft/d and 40 ft,
respectively. The well is pumping near a stream at a rate of 20,000 ftVd. Depicted in
Figures 8.16a and 8.16b are 10,000-day capture zones obtained using the semi-analytical and
numerical GPTRAC options," respectively. Note that a zero drawdown condition was
assumed on the boundary of a uniform rectangular finite element grid (AX=Ay=50ft).~
comparison of Figures 8.16a and 8. 8.16b shows very good agreement between the semi-
analytical and numerical solutions to the problem.
-------
General Particle Tracking Module (GPTRAC) £.37
Figure 8.16
Capture Zones for Semi-Confined (Leaky) Aquifer Using the Semi-Analytical (a)
and Numerical (b) GPTRAC Modules
(FT)
2500
1666
833
(FT)
2500
1666
833
_L
(FT)
833 1666 2500
(a)
(FT)
833 1666 2500
(b)
-------
9.0 Uncertainty Analysis
ModulelMONTEC
' - . - -
9.1 Introduction
The computational modules described previously delineate various types of capture
zones for a variety of ground-water flow scenarios. Each module requires a number of user-.
specified input parameters such as transmissivity, porosity, etc. Typically, the values of these
parameters are not known exactly due to measurement errors, data limitations, and/or
inherent spatial and temporal variability. Often, no tests have been performed to determine
the aquifer hydraulic properties for a given site, and the values used will be an "educated
guess" on the part of the delineator.
If the uncertainties associated with one or more input variables are deemed significant,
it may be appropriate to express the variable(s) in a random context. Random variables
may be thought of as variables with more than one potential value, and therefore they are
described by probability distributions. The most likely value that a random variable will
assume is the mean of its associated probability distribution. The variance of a probability
distribution is a measure of the range of values that the corresponding random variable may
assume.
If one or more input variables of the capture zone model are considered random, then
Monte Carlo analysis may be used to assess the effect of the uncertain input parameters on
the model output (capture zones). For reasons discussed in Section 9.4.1, only time-related
capture zones for a single pumping well in a homogeneous, isotropic, confined or semi-
confined aquifer with steady and uniform ambient ground-water flow may be analyzed using
the Monte Carlo option.
9.2 Uncertainty Analysis Objectives
The objective of the Monte Carlo approach is to estimate the uncertainty in capture
zone size and shape given the uncertainty in the model input parameters. Alternatively, the
objective is to estimate the cumulative probability distribution of the capture zone boundary
-------
9-2 Uncertainty Analysis Module (MONTEC)
given the probability distribution of the input parameters. Thus, if R represents some
measure of the capture zone size and X represents the vector of all model inputs:
R = W(X) (9.1)
where W represents the WHPA model. Note that only some of the components of X may
vary in an uncertain way, i.e., they may be random variables defined by a probability
distribution. Thus, the goal is to calculate the cumulative distribution function FR(R') given
a probabilistic characterization of 2L FR(R') is defined as
FR(R') = Probability (R < R') (9.2)
where R is a given capture zone extent.
9.3 The Monte Carlo Technique
Given a set of values for each of the model input parameters, Xj, X2 . . . Xn, the
WHPA model computes the capture zone boundary R:
R = W(Xj, X2, X3 ... XJ (9.3)
Application of the Monte Carlo simulation procedure requires that at least one of the
input variables, Xj . . . X^, be uncertain and that the uncertainty be represented by some
probability distribution. The method involves the repeated generation of pseudo-random
values of the uncertain input variable(s) (drawn from the known distribution and within the
range of any imposed bounds) and the application of the model using these values to
generate a series of capture zones. The term pseudo-random (hereafter referred to as
random) number is used because truly random numbers can not be generated on the
computer; the computer is only capable of generating sets of numbers that exhibit the
characteristics of a given probability distribution. The capture zones are then statistically
analyzed to yield a cumulative probability distribution. Thus, the various steps involved in
the application of the Monte Carlo simulation technique are:
1. Selection of representative probability distributions for the input
variables that are to be considered uncertain.
2. Generation of a random number from each of the distributions
selected in (1). These values represent a possible set of values (a
realization) for the input variables.
-------
Uncertainty Analysis Module (MONTEC) n_i
3. Application of the WHPA model to compute a capture zone.
4. Repeated., application of steps (2) and (3) for a specified (large)
number of iterations.
5. Presentation of the series of output values generated in step (3) as a
cumulative probability distribution function (CDF).
6. Analysis and application of the cumulative probability distribution as
a tool for decision making.
The number of Monte Carlo iterations that must be conducted is problem specific, but in
general a minimum of at least several hundred iterations should be used (see section 9.4.2).
9.4 Application of Monte Carlo Technique to Capture Zone Analysis
9.4.1 Methodology
Figure 9.1 illustrates the methodology used to link the Monte Carlo method to the
WHPA model. Figure 9.1a represents one realization for a time-related capture zone about
a pumping well. The capture zone would be obtained by running MONTEC for a single
iteration subsequent to the generation of a random number for each of the uncertain input
variables.
Because MONTEC will calculate many capture zones during a Monte Carlo run, the
shape and areal extent of each capture zone needs to be saved so that a comparison of all
of the capture zones can be made. A series of radial lines, with the well as their common
origin, are used to accomplish this task. Each line has associated with it a unique angle, aj.
The length of each radial line segment and its corresponding angle compose a piecewise
constant representation of the capture zone. Therefore, for each capture zone generated,
the radial distance from the well to the capture zone boundary for n lines at each 04
(i = l,n) is saved by the code. The total number of values saved is N x n, where N is the
number of Monte Carlo runs.
-------
9-4
Uncertainty Analysis Module (MONTEC)
Figure 9.1
Monte Carlo Method Applied to Capture Zone Analysis
Series of Radial Lines Used to Save Capture Zone Boundary (a)
and Cumulative Distribution Function For One Radial Line (b)
LL
Q
O
(a)
= 20
Radial Distance
(b)
-------
Uncertainty Analysis Module (MONTEC) o
Once N capture zones have been generated by MONTEC, the N distances saved for
each radial line are sorted into ascending order. The resulting tabulation is the experimental
cumulative distribution function (CDF) for each line. A sample CDF, which corresponds
to the radial line defined by a^ in Figure 9.la, is shown in Figure 9.1b.
Once the CDF is obtained, any percentile of radial distance may be chosen. For
example, the 85th percentile radial distance implies that 85 percent of the radial distances
for the current line are less than or equal to that value. The 85th percentile capture zone
would be determined by finding the 85th percentile distance for each radial line, and then
connecting the endpoints of each line after plotting them at the appropriate angles.
One can see from Figure 9.la that the quality of the "preservation" of each capture
zone depends upon the number and spacing of the radial lines. In addition, the number of
lines required to preserve a capture zone increases with the complexity of the capture zone
shape. However, it is inconvenient to use large numbers of radial lines due to increased
storage and computational requirements. For these reasons, the Monte Carlo option of the
WHPA model may only be used for simple WHPA delineation scenarios. The following
restrictions apply to use of the MONTEC module:
Only the capture zone for a single pumping well in a confined or
semi-confined aquifer may be analyzed using the Monte Carlo
. technique. The capture zone type must be time-related.
The aquifer must be infinite in areal extent (effects of nearby stream
or barrier boundaries may not be assessed).
The portion of MONTEC that calculates capture zones is similar to
MWCAP, and therefore most of the restrictions applicable to
MWCAP are applicable to MONTEC as well.
9.4.2 Input Requirements
The input requirements for MONTEC are essentially the same as those of MWCAP
for the case of a single, time-related capture zone for a well pumping in an aquifer of
infinite areal extent (no boundaries). Additional input requirements are the number of
Monte Carlo runs (NRUN), the desired capture zone percentile(s) (PERCEN; up to five
may be specified), and the statistical attributes (mean and standard deviation) of the
probability density functions (PDF's) assigned to the uncertain input variables. Note that
MONTEC may also handle aquifers subject to leakage through a confining bed, whereas
MWCAP can not.
-------
9-6 Uncertainty Analysis Module (MONTEC)
A maximum of 1,000 Monte Carlo runs may be specified. In general, it is recommend-
ed that for all final analysis the maximum number of iterations be used. For screening
purposes, a smaller, number of Monte Carlo iterations may be sufficient. For several
example problems analyzed using MONTEC, it was observed that the Monte Carlo capture
zones obtained using 200-500 iterations were almost identical to those obtained using 1,000
iterations. Users that have access to a powerful PC such as the 80386 series machines will
be able to run a Monte Carlo problem with a large number of iterations within several
hours. Users with less computing power (e.g. an IBM XT class PC) may have to wait
overnight for results. Note that proper selection of the number of pathlines delineated and
the maximum spatial step length are very important when MONTEC is used.
The input variables that may be considered uncertain are
Well discharge rate
Hydraulic conductivity
Hydraulic gradient
Effective porosity
Aquifer thickness
Confining layer hydraulic conductivity
Confining layer thickness
Any one, some combination of, or all of the above parameters may be represented as a
probability distribution during a Monte Carlo run.
9.4.2.1 Distribution Types
If one or more input variables is considered uncertain, the probability distribution type
of the variable(s) must be specified. There are six distribution types to choose from in the
WHPA model:
Normal
Lognormal
Exponential
Uniform
Log 10 uniform
Empirical
-------
Uncertainty Analysis Module (MONTEC)
The normal, exponential and uniform distributions are commonly used in practice and a
description of each may be found in any introductory statistics text. The lognormal and
log 10 uniform distributions involve a mathematical transformation of the original data. The
result of the mathematical transformation, referred to as the transformed data, follows a
specified probability distribution that the untransformed (original) data does not. The two
transformations are defined as:
Lognormal: Y = in(X) (9-4)
LoglO uniform: Y = log10(X) (9-5)
where in denotes natural log, logjQ is logarithm to the base 10, X is an untransformed
random variable, and Y is a transformed random variable. The variable Y will have a
normal distribution when the lognormal transformation applies, and a uniform distribution
when Iog10 uniform transformation applies. The untransformed random variable X
represents the original data prior to any logarithm operations.
When a normal or uniform random number is generated from one of the logarithm
distributions in the Monte Carlo module, the generated value must be transformed from the
transformed space back to untransformed space. This is accomplished through the inverse
transforms defined by:
Lognormal: X = exp(Y) (9-6)
LoglO uniform: X = 10^ (9-7)
where exp denotes the exponential. Note that Y is the random value generated by
MONTEC, but X is the value that is used by the WHPA model.
The empirical distribution is a set of coordinates that define an observed cumulative
distribution function (CDF) for a random process. The y-axis of the CDF has a range of 0-1
and represents the probability that any value X will be less than or equal to y. The x-axis
of the CDF represents values of the random variable of interest.
The normal and lognormal distributions require the user to specify the mean and
standard deviation of the distribution. Note that for the lognormal distribution the mean
and standard deviation are entered in untransformed space, and the transformation of these
parameters to log space is made by the code using
-------
9-8 Uncertainty Analysis Module (MONTEC)
- 0.5 (oj) (9-8)
and
where jUy and ay are the mean and standard deviation in untransformed space, and MT and
0T are the mean and standard deviation in log space, respectively.
The uniform and log 10 uniform distributions require minimum and maximum values.
The log 10 uniform distribution bounds are input in untransformed space. The exponential
distribution requires only one parameter; the mean of the distribution. For the empirical
distribution, the user must input the coordinates of the cumulative distribution function
(minimum 2 pairs, maximum 50 pairs) which is subsequently treated as a piecewise linear
curve.
Once the distribution type for each uncertain input variable has been specified, the
MONTEC module generates a random number for each variable for each Monte Carlo run.
If the generated value lies outside any imposed bounds (see next section), the number is
regenerated. The random number generation routines used in MONTEC are documented
in McGrath and Irving (1973):
9.4.2.2 Distribution Bounds
The user may assign limiting upper and lower values that a random variable may
assume if the variable is represented by the normal, lognormal or exponential distribution
types. Upper and lower bounding values are inherent in the definitions of the other
distribution types. All upper and lower bounds are input in untransformed space.
Extreme caution should be used when specifying limiting values for distributions that
have none by definition. If the bounds lie too close to the mean of the distribution, the
available sampling space may be too limited and the nature of the statistical distribution will
not be preserved during the sampling (random number generation) process. Bounding
values should be used only to avoid physically unrealistic scenarios that could be caused by
an extreme value generated for a random variable. For example, if well discharge is
characterized using a normal distribution, an upper bound might be specified for Q such that
complete dewatering of the aquifer in the vicinity of the well would not occur.
-------
Uncertainty Analysis Module (MONTEC) no
The upper and lower bounding values imposed automatically by MONTEC are as
follows:
Well discharge must be greater than zero
Hydraulic conductivity must be greater than zero
Hydraulic gradient must be greater than or equal to zero
Porosity must be greater than zero and less than one
Aquifer thickness must be greater than zero
Aquitard hydraulic conductivity must be greater than zero
Aquitard thickness must be greater then zero
9.4.2.3 Statistical Correlations Between Input Parameters
If more than one input variable is considered uncertain, statistical correlations that may
exist between variables are neglected. Each random variable is assumed to be independent
of every other random variable during the random number generation process. This
assumption is based primarily upon practical rather than technical considerations. In
practice, it is generally extremely difficult due to data limitations to derive statistical
correlations (covariances) between input parameters.
However, statistical correlations between aquifer parameters certainly exist. For
example, it follows from Darcy's Law that hydraulic conductivity, aquifer thickness and
hydraulic gradient are correlated with one another since per unit width of aquifer
Q = Ti (9-10a)
where
T = Kb (9-10b)
Therefore high values of T are associated with small values of i, and high values of i (steep
gradients) are associated with small T values. Furthermore, the variable T is a direct
function of K and b.
These correlations may pose a problem if statistical independence between variables
is assumed. For a given Monte Carlo run, there may be a significant probability that a low
value of i will be generated in conjunction with a low value of the product Kb. Such a
situation might be unreasonable in light of the previous discussion.
-------
9-10 Uncertainty Analysis Module (MONTEC)
To avoid extremely unrealistic hydrogeological scenarios, an automatic check based
upon drawdown at the pumping well was implemented into MONTEC. If the maximum
allowable drawdown specified by the user is exceeded, the current set of random values
generated by the code is ignored and a new set is generated until the drawdown criterion
is met. Hydraulic parameter combinations that lead to excessive drawdowns at the well
during a Monte Carlo run are output to the file DD.LST (drawdown list).
Drawdown at the pumping well is calculated using the Theim equation:
5 2^5,"-(7^ <9-n>
where s is drawdown from the initial piezometric surface at radial distance r from the
pumping well, Reis the radial distance from the well at which there is no drawdown (the
radius of influence of the well), and Q, K and b are as defined previously. For each Monte
Carlo run the values of Q, K and b are known. Because drawdown at the well is being
calculated, r corresponds to the effective radius of the well bore. In order to calculates, the
only remaining unknown variable is Re.
Unfortunately, there is no sound quantitative method to estimate Re. This parameter
is often estimated from experience or actual field observation. Two approximate methods
for estimating Reare implemented in MONTEC. The method used by the code depends
upon the selected aquifer type (confined or leaky-confined).
Because Reis related to the recharge processes that supply ground water to the well,
the following method is used to estimate Rein MONTEC when the confined aquifer option
is selected. Using equation (9-10) and the relation
Q = qW (9-12)
where W is the width of aquifer required to supply the total recharge to a well pumping at
rate Q, Re is calculated from
Re = 0.5 W = 0.5 ^ (9-13)
Fortunately, Reand r appear in the form of a log function in equation (9-11), and therefore
even a large error in the estimation of one of these parameters does not appreciably affect
-------
Uncertainty Analysis Module (MONTEC) Q_J j
the drawdown computed using equation (9-11). The drawdown check is not made if there
is no ambient ground-water flow (i = O).
If the leaky-confined aquifer option is used in MONTEC Reis obtained from
Re = 1.23 A (9-14)
and
A* = (9.15)
K'
where A is a characteristic length of the leaky aquifer called the leakage factor (Bear, 1979),
b' is the thickness of the semi-permeable confining unit K' is the hydraulic conductivity of
the semi-permeable confining unit, and the other variables are as defined previously. In this
case, the pumping rate of the well, Q, will be equal to the volume of leakage per unit time
supplied to the aquifer within the circle defined by radius Reln this case the drawdown
calculation does not depend on the presence of ambient ground-water flow.
The maximum drawdown that may not be exceeded is input by the user; an
appropriate value for this parameter must be derived on a case-by-case basis. It might be
obtained from field observation (e.g. the maximum drawdown observed in a well in a
particular region), or from theoretical considerations (e.g. the piezometric surface can not
drop below the well bottom or some other datum such as the base of an overlying confining
unit). If the maximum allowable drawdown is set to zero, which is the default value, the
drawdown check is not conducted.
The procedure of eliminating certain variable sets due to excessive drawdown will
indirectly induce correlation between model parameters. For example, combinations of low
T and low i values create large drawdowns, and therefore the elimination of parameter sets
in which this condition occurs tends to allow only combinations of low T with high i and high
T with low or high i. Note that if the maximum allowable drawdown is set to a large value,
the condition of complete independence between parameters will again exist because the
drawdown criterion will never be exceeded.
-------
9.12 Uncertainty Analysis Module (MONTEC)
9.4.2.4 Data Analysis
The data required to derive the probability distributions for uncertain input variables
can be obtained from a variety of sources. The most obvious source of data is field tests
conducted in the aquifer under study. For example, if 50 aquifer tests have been conducted
for a region of interest, standard statistical techniques might be used to estimate a
probability distribution type for hydraulic conductivity. A lognormal distribution might be
found and this distribution type, along with its mean and standard deviation, could be input
into MONTEC (note that a lognormal distribution is commonly observed for the hydraulic
conductivity of an aquifer).
In general it is uncommon that enough field data exists to obtain estimates of PDF's
via traditional methods. Additional sources of data may become available if the WHPA
delineator is willing to make certain assumptions. For example, suppose that a delineator
would like to determine the lo-year capture zone for some municipal' water supply wells
located in a shallow, river-valley aquifer. The delineator would like to consider hydraulic
conductivity to be a random variable because no hydraulic conductivity values estimated
from field tests are available in the vicinity of the wells. However, there are other not too-
distant aquifers similar in morphology (stream valley settings, similar thicknesses, similar
aquifer materials, etc.) for which many hydraulic conductivity estimates are available. If the
assumption is made of some equivalence between aquifers due to similar physical
characteristics, the measurements made in the surrounding aquifers may be used to obtain
a PDF of hydraulic conductivity for the aquifer of interest.
Barring a reasonable quantitative means of obtaining sample PDF's, a distribution type
may be chosen using qualitative professional judgement. Often the uniform distribution is
used if a "guess" at the distribution type is made. A uniform random variable has an equal
likelihood of assuming any value between the specified upper and lower limits of the
distribution. The lower and upper distribution limits correspond to the minimum and
maximum parameter values thought to be reasonable for the given site. Other distributions
(e.g., normal) may be hypothesized as well.
9.4.3 Output
Output of the MONTEC module consists of a capture zone that corresponds to any
specified percentile value (85th, 90th, etc.). The 90th percentile capture zone implies that
there is a 90 percent chance that the "actual" capture zone boundary for the well will be less
than or equal to the areal extent of the delineated capture zone. Conversely, there is a 10
percent chance that the actual capture zone may exceed the bounds of the delineated
capture zone. Up to five percentiles may be specified in the MONTEC preprocessor, and
-------
Uncertainty Analysis Module (MONTEC) Q_"|3
the outline of each corresponding capture zone will be computed and plotted. The
percentile value used for WHPA delineation purposes is a policy decision; however,
generally values of 90 or 95 percent are used when uncertainty analysis is applied in a
regulatory framework.
The capture zones produced by MONTEC may have slightly distorted shapes due to
the piecewise linear representation used (Figure 9.2). The distortion will generally be minor,
and the actual shape of the capture zone is usually discernable.
In addition to the standard output file MONTEC.OUT, the MONTEC module
automatically creates three additional files that may be of interest to the user. These files
and their contents are listed below.
CAPTUR.VAR - Contains the MONTBC input parameter* Q, K, i, b; KV
''-. . b', and 0 generated for each Monte Carlo Iteration.
RANBND.LST . - Contains a listing of the generated random values that
. were rejected because they did not lie within any
i imposed bounds. Fife will be; empty if no lower and
; . upper bounds are set. . \
DDJLST . - Contains a listing of input parameter combinations that
' caused excessive drawdown at the pumping well (as
determined by the maximum permitted drawdown input
parameter).
9.5 Example Applications
Two hypothetical Monte Carlo problems were run to illustrate the capabilities of the
Monte Carlo option. Note that the MONTEC module is substantially more computationally
intensive than the other WHPA model modules. The first example problem ran in about
20 minutes on a Dell 80386-20 mgHz PC and in about 1.5 hours on an Everex 80286-12
mgHz PC. The second example ran in about 1 hour and 20 minutes on the Dell 80386
machine. Both machines were equipped with math coprocessor. Run times would be
substantially longer on an IBM XT model PC.
9.5.1 Hypothetical Example 1
For the first example problem, the input parameters well discharge, effective porosity
and aquifer thickness were considered uncertain. The input required by MONTEC is listed
in Table 9.1. Since this is only a hypothetical example, the maximum allowable drawdown
was set to a large value (40 m). The three percentiles specified were the 5th, 90th and 99th.
-------
9-14 Uncertainty Analysis Module (MONTEC)
Figure 9.2
Distorted Shape of MONTEC Capture Zone
Legend:
Actual Capture Zone
Capture Zone Outline Saved by MONTEC
-------
Uncertainty Analysis Module (MONTEC) 0-15
Table 9.1
Input Parameters for First MONTEC Example
Aquifer Type: Confined
Units: Meters and Days
Problem Boundaries: XMIN = YMIN = 0.0 m
XMAX=YMAX=5,OOOm
Maximum Step Length: 50 m
Well Location: x = 750 m; y = 2,500 m
Well Radius : 0.25m
Angle of Ambient Flow: 180°
Maximum Permissible Drawdown: 40m
No. of Monte Carlo Runs: 1,000
No. of capture zone percentiles to be solved for: 3
Percentile Values: 5, 90, 99
Monte Carlo Variables:
Distribution Standard Lower Upper
Parameter Type Mean Deviation Bound Bound
Well Discharge Normal 3,000 mVd 500 mVd 0 0
Hydraulic Conductivity Constant 15 m/d
Gradient Constant 0.0015
Porosity Uniform -- - 0.2 0.25
Thickness Empirical* ~
Empirical distribution for thickness
Thickness
40m
50m
62m
70m
0.0
0.5
0.8
1.0
Capture Zone Time = 3,650 days (10 years)
Number of Pathlines = 20
-------
9-16 Uncertainty Analysis Module (MONTEC)
Note in Table 9.1 that the well discharge has a normal distribution with the upper and
lower bounds on the distribution set equal to zero. When the lower and upper bounds of
a distribution are set-m zero, MONTEC does not check the random number generated prior
to each iteration to ensure that it lies within the specified bounds.
A uniform distribution with lower and upper bounds of 0.2 and 0.25 respectively was
specified for aquifer porosity. Therefore, for each Monte Carlo iteration, the random value
generated for porosity had an equal likelihood of lying anywhere within the range 0.2-0.25.
The empirical distribution used for aquifer thickness is listed at the bottom of
Table 9.1. Reading the table, one finds that the smallest possible value of aquifer thickness
is 40 m, there is a 50 percent chance that the aquifer thickness is less than or equal to 50
m, and so on. Such a distribution might be derived for a particular site using available well
logs for the region.
The 5th, 90th and 99th capture zone percentiles are plotted in Figure 9.3. The 5th and
the 99th percentile capture zones lie near the smallest and largest Monte Carlo capture
zones respectively. Obviously, the 90th percentile capture zone lies between these two
extremes.
9.5.2 Hypothetical Example 2
The second Monte Carlo example is loosely based on statistical distributions of flow
parameters derived for the Quadra Sand, British Columbia (Smith, 1981). The hypothetical
setting is a stratified sand aquifer with a single pumping well withdrawing water for public
supply at the rate of 1 mgd (1.34 x l(r ft /d). The aquifer ranges in thickness from 80-120
ft. The spatial variability of hydraulic properties is ignored and the statistical distributions
are viewed as representing the uncertainty associated with parameters describing a
homogeneous aquifer. The uncertain input parameters are hydraulic conductivity, effective
porosity and aquifer thickness. Table 9.2 is a summary of the MONTEC input for this
example.
The 50th, 80th and 95th capture zone percentiles for this example are presented in
Figure 9.4. The most striking feature of these capture zones is their atypical pear shape.
This shape results from the composite effects of averaging many different capture zones of
various shapes and sizes.
-------
Uncertainty Analysis Module (MONTEC) 0_1 7
Figure 9.3
5?, 90th, and 99th Percentile Capture Zones For The
First Hypothetical Example Computed Using MONTEC
(M)
5000
4000
.th
3000
2000
1000
(M)
1000
2000
3000
4000
5000
-------
9-18 Uncertainty Analysis Module (MONTEC)
Table 92
Input Parameters For Second MONTEC Example
Aquifer Type: Confined
Units: Feet and Days
Problem Boundaries: XMIN = YMIN = 0.0 ft
XMAX=YMAX=10,000ft
Maximum Step Length: 50 ft
Well Location: x =2000 ft; y = 5,000ft
Well Radius : 1 ft
Angle of Ambient Flow 180°
Maximum Permissible Drawdown: 100 ft
No. of Monte Carlo Runs: 1,000
No. of capture zone percentiles to be solved for: 3
Percentile Values: 0.5, 0.8,0.95
Monte Carlo Variables:
Distribution Standard Lower
Parameter Type Mean Deviation Bound
Well Discharge
Hydraulic Conductivity
Gradient
Porosity
Thickness
Constant 134,000^
Lognormal 142 ft/d 50 ft/d
Constant 0.002
Lognormal 0.43 0.03 0
Uniform - -- 80ft
Upper
Bound
0
120ft
Capture Zone Time = 3,650 days (10 years)
Number of Pathlines = 40
-------
Uncertainly Analysis Module (MONTEC) 9-19
Figure 9.4
The 50th, 80th, and 95th Capture Zone Percentiles Calculated Using
MONTEC For The Second Hypothetical Example
(FT)
10000
8000
6000
4000
2000
0
0
2000
(FT)
4000
6000
8000 10000
-------
9-20 Uncertainty Analysis Module (MONTEC)
To illustrate the cause of the pear shapes more clearly, two additional capture zones
were delineated using MWCAP (recall that much of the MONTEC module is based upon
the MWCAP module). The smallest and the largest transmissivity values generated during
the 1,000 MONTEC iterations were used as input for the first and second MWCAP runs
respectively. The values were extracted from the standard MONTEC output file
CAPTUR.VAR. All of the input parameters for each MWCAP run are presented in
Table 9.3.
Figure 9.5 portrays the MWCAP results. For the run that used the low value of
transmissivity as input, a nearly circular capture zone was delineated, while an oblong
capture zone resulted from using the high value of transmissivity as input. These capture
zones may be conceptualized as the two "extreme" or "end-member" capture zones
delineated during the MONTEC run. The origin of the pear shape becomes immediately
apparent when the circular and elongated capture zones are visually laid one atop the other.
Note that the pear shape will only occur when input parameters that effect the determina-
tion of ambient ground-water flow (i.e. hydraulic conductivity, aquifer thickness, and
hydraulic gradient) are considered uncertain.
To demonstrate the semi-confined aquifer option in MONTEC, the second example
was rerun using the leaky aquifer option. The input for this run was the same as that
presented in Table 9.2 with the exception that the hydraulic conductivity and thickness of
the confining bed were also specified as uniform random variables with bounds of 0.1-0.001
ft/d and 10-25 ft respectively. The results of this run are presented in Figure 9.6. A
comparison of Figures 9.4 and 9.6 shows the general effect of decreased capture zone size
due to leakage through the confining bed over the zone of influence of the well.
-------
Uncertainty Analysis Module (MONTEC) 9-21
Table 9.3
MWCAP Input Parameters For Two Monte Carlo End-Member Capture Zones
For Each Run:
Units = ft and days
XMIN = 0 YMIN = 0
XMAX = 10,000 YMAX = 10,000
Step Length =50
No. of Pumping Wells = 1
Input Parameter
X Coordinate
Y Coordinate
Discharge
Transmissivity
Hydraulic Gradient
Angle of Ambient Flow
Porosity
Aquifer Thickness
Boundary Type "
Distance from Well
to Boundary
Boundary Angle
Capture Zone Type
Travel Time Value
No. of Pathlines
Capture Zone Plotting
Option
Run No. 1
2,000
5,000
134,000
4,218
0.002
180"
0.43
87
None
N/A
N/A
Time-Related
3,650
20
Yes
Run No. 2
2,000
5,000
134,000
41,840
0.002
180°
0.43
109
None
N/A
N/A
Time-Related
3,650
12
Yes
-------
9-22 Uncertainty Analysis Module (MONTEC)
Figure 9.5
Capture Zones Delineated Using MWCAP for (a) Lowest Randomly Generated
Transmissivity Value and (b) Highest Randomly Generated Transmissivity Value
(FT)
10000
8CCG
6000
4000
2000
(FT)
0 2000 4000 6000 8000 10000 '
(a)
10000
8000
6000
4000
2000
(FT)
0 2000 4000 6000 BOOO 10000
(b)
-------
Uncertainty Analysis Module (MONTEC)
Figure 9.6
The 50th, 80* and 95th Capture Zone Percentiles Calculated Using MONTEC
for the Second Hypothetical Example with Vertical Leakage
(FT)
10000
8000
6000
4000
2000 -
0
I
(FT)
2000 4000 6000 8000 10000
-------
10.0 Problem Conceptualization vs.
Reality: Potential for Abuse
10.1 Introduction
In order to understand what the capture zones output by the WHPA model truly
represent, the effects that the input parameters and the underlying model assumptions have
on capture zone morphology must be understood. The effect of the model input parameters
and the underlying model assumptions on capture zone size and shape is the topic of this
chapter. The majority of the material presented herein is applicable to the RESSQC,
MWCAP, MONTEC and GPTRAC semi-analytical modules. The GPTRAC numerical
module is limited only by the capabilities of the numerical model used to obtain the input
hydraulic head field.
10.2 Effects of Input Parameters on Capture Zone Morphology
The effects of the input parameters transmissivity (T), aquifer thickness (b), hydraulic
gradient (i), aquifer porosity (0), well discharge (Qw), and the angle of ambient flow (a) on
capture zone morphology may be assessed using Darcy's Law and the mass balance
equation. Darcy's Law may be written
Q = KiA (10.1)
where Q is discharge, K is the hydraulic conductivity of the porous medium, and A is the
cross-sectional area of flow. If Q is assumed to represent the discharge per unit width of
aquifer, it holds that
(10.2)
or
Qa = Ti (10.3)
-------
10-2 Problem Conceptualization vs. Reality: Potential for Abuse
where Qa is the ambient ground-water flow per unit width of aquifer. The average pore
water velocity (seepage velocity) for a fluid particle moving through the porous medium may
be written as
v = q/0 (10.4)
where q = Q/A = Ki is called the specific discharge or Darcy velocity.
The equation of mass balance may be written
I-O = AS (10.5)
where
I = inflow per unit time
O = outflow per unit time
AS = change in storage per unit time
Equation 10.5 may be applied to any hydrologic system. For the hydrogeologic systems to
which the WHPA model is applicable, I represents ground-water recharge moving into the
aquifer (i.e., ambient ground-water flow, induced infiltration from streams, injection well
recharge, areal recharge, leakage from an overlying confining unit), and O represents
ground-water discharge from the aquifer (i.e., well discharge, ground-water discharge to
streams, ambient ground-water flow that is not captured by the pumping well(s)). Because
the flow fields considered are assumed to be at steady-state, the AS term is equal to zero.
Therefore, the mass balance equation for the WHPA model may be written
I = O (10.6)
Equation 10.6 implies that water discharged from the aquifer by the pumping well(s) must
be obtained from a combination of ambient ground-water flow within the aquifer, induced
infiltration from streams, and injection well recharge. Figure 10.1 is a schematic
representation of a mass balance analysis for a simple hydrogeologic system.
With equations 10.1, 10.3, 10.4 and 10.6 in hand, the effects of T, b, i, 0, Qw and a on
capture zone morphology cart be assessed. First of all, consider an aquifer of infinite areal
extent (no barrier or stream boundaries) with no recharge wells. If a pumping well exists
in the aquifer, its Qwmust be derived entirely from the ambient ground-water flow.
-------
Problem Conceptualization vs. Reality Potential for Abuse 10-3
Figure 10.1
Mass Balance Analysis For Simple Hydrogeological System
GENERAL SYSTEM
0
AS
= 0
I
I = 0
Aquifer with uniform ambient flow
Q,
Qi
-------
10-4 Problem Conceptualization vs. Reality Potential for Abuse
Therefore, as Qw increases the capture zone must expand in all directions to capture a
larger portion of the ambient aquifer flow, Qa. Conversely, as Qw decreases the areal extent
of the capture zone will decrease.
From equation 10.3 one can see that the ambient ground-water flow per unit width of
aquifer is the product of T and i. Therefore, for a constant Q^ increasing T and/or i will
decrease the capture zone size because the well can derive a larger portion of water from
each unit aquifer width contained within the capture zone. Decreasing T and/or i will
increase the areal extent of the capture zone. The capture zone for a given scenario will not
change if T and i are adjusted in opposite directions in such away that Qa remains constant.
The above generalizations hold if a stream or barrier boundary exists within the
aquifer. If a stream boundary is in close enough proximity to the well, it will act as a source
of water to the well and the areal extent of the capture zone will be smaller than it would
otherwise be. In addition, the stream may act as a limiting condition on the capture zone
extent. If a barrier boundary is in close enough proximity to the well, it will act as a barrier
to ground-water recharge to the well and the areal extent of the capture zone will be larger
than it would be otherwise. The orientation of the barrier boundary in conjunction with the
angle of ambient flow will dictate the direction in which the capture zone must expand.
Steady-state capture zones are not dependent upon porosity, 6, or thickness, b. For
the time-related capture zone solutions, 9 and b do not affect the position or orientation of
the pathlines; they only affect the velocity of the fluid flow along each pathline. Equation
10.4 indicates that fluid velocity is inversely proportional to porosity (i.e., decreasing porosity
increases fluid velocity). Therefore, when solving for time-related capture zones, decreasing
the porosity will have the same effect as increasing the time for which the pathlines are
delineated and vice versa; the same is also true for b. Note that the increased seepage
velocity due to decreasing 6 or b leads to more conservative capture zones because the
capture zone will become larger for a given time.
The aquifer parameters T and i will also effect the velocity of ground-water flow.
Recall that K = T/b. Equations 10.1 and 10.4 indicate that, if B remains constant, the
seepage velocity will increase or decrease in direct proportion with K and/or i. However,
unlike porosity, K and i will also effect the location and orientation of pathlines because they
determine the ambient flow rate, Qa.
-------
Problem Conceptualization vs. Reality Potential for Abuse 10-5
For steady-state capture zones, decreasing Tori will lead to more conservative capture
zones for a well pumping at a constant rate Q^ However, for time-related capture zones
simple generalizations should not be made regarding the effect of these parameters.
Although increasing Tori may lead to higher seepage velocities and therefore larger time-
related capture zones, the increased ambient flow rate will tend to limit the areal extent of
the capture zone.
The effect of one final input parameter, a, is readily apparent. A capture zone will
always orient itself so that it extends upgradient in the direction opposite to that of ambient
flow (i.e., if a =45° (NE), the bulk of the capture zone will extend in the direction of 225°
(SW)). Many interesting capture zone shapes may occur when a stream or barrier boundary
exists within the ambient flow field (Figure 10.2).
For the sake of simplicity, the effects of areal recharge and vertical leakage on capture
zones in unconfined and leaky-confined aquifers, respectively, were not addressed in the
preceding discussion. The general effect of these processes, obviously, is to decrease the
capture zone size. In the case of a leaky-confined aquifer, larger values of the confining bed
hydraulic conductivity (K'), and smaller values of the confining bed thickness (b'), lead to
greater quantities of leakage supplied to the aquifer through the confining bed and therefore
smaller capture zones.
103 Effects of Boundary Conditions on Capture Zones
The presence of a boundary condition in the proximity of a pumping well can
dramatically alter the shape and areal extent of the well's capture zone., The boundary
condition options incorporated in the WHPA model are designed to simulate those that
often occur in the field. However, in many cases encountered the boundary conditions as
formulated in the WHPA model may not adequately represent the natural system, and
therefore their use may lead to erroneous conclusions regarding capture zone size and
shape. Whenever the stream or barrier boundary conditions are applied, the user must be
aware of the underlying assumptions and how they affect capture zone morphology.
Of the two boundary types available, the stream boundary condition assumptions are
the most likely to be violated. Figure 10.3a is a schematic diagram of how the WHPA model
simulates a stream boundary the stream is assumed to fully penetrate the aquifer and
therefore can provide water to the aquifer over the entire saturated thickness, b. With this
conceptualization, the capture zone may expand along the stream boundary axis, but may
never extend across the stream.
-------
10-6 Problem Conceptualization vs. Reality Potential for Abuse
Figure 10.2
Miscellaneous Capture Zones
(M)
4000
3200
2400
1600
800 -
BARRER
BOUNDARY
STREAM
BOUNDARY
(M)
0
800
1600
2400
3200
4000
-------
Problem Conceptualization vs. Reality Potential for Abuse j 0
Figure 10.3
Stream Boundary As Represented by the WHPA Model (a), and
Stream Boundary As Often Encountered in Practice (b)
1 1
b
I
yy
^ J
1 *T
fc.
».
:--
: -«
: "*
Stream
(a)
-------
10-8 Problem Conceptualization vs. Reality Potential for Abuse
Figure 10.3b diagrams the stream boundary condition as it is often observed in
practice. Note that the stream is partially penetrating it does not exist over the entire depth
of the aquifer. In such a case the well may be pumping at a larger rate than that which can
be sustained by induced infiltration from the stream into the aquifer. When such a situation
occurs, the capture zone will extend across the stream boundary, and the pumping well will
draw a component of recharge from the opposite side of the stream. The occurrence, or
potential for occurrence, of this phenomena should be carefully examined whenever the
stream boundary condition is used. Morrissey (1987) has examined this problem in detail.
The assumption of a fully penetrating barrier boundary is less apt to be violated in
practice than is the fully penetrating stream assumption. Figure 10.4a is a schematic
diagram of how the WHPA model simulates a barrier boundary, the impermeable rock
formation is assumed to extend, and therefore act as a barrier to flow, over the entire
saturated thickness of the aquifer. With this conceptualization, the capture zone may never
extend across the boundary face.
In some instances, the barrier boundary may not fully penetrate the aquifer and the
capture zone may actually extend across the boundary face (Figure 10.4b). Detailed
geological and hydrogeological data would have to be available to identify this situation.
Whenever a stream or barrier boundary condition is partially penetrating, a potentially
significant vertical component of flow is introduced into the flow system. The WHPA model
is designed to simulate horizontal, two-dimensional flow and therefore the effects of partially
penetrating boundaries can not be assessed using the model. Even if the model is run
assuming that no boundary is present, the resulting capture zone may still not be
conservative. In order to rigorously assess the effects of partially penetrating boundaries,
a numerical model should be used.
10.4 Analysis of Simulation Results
All capture zones output from the WHPA model should be viewed in light of the
assumptions applied (e.g., steady flow, homogeneous aquifer, etc.). Although the model
assumptions may not allow an exact representation of many hydrogeological systems, the
assumptions will often be close enough to reality so that a reasonable capture zone solution
may be obtained. When one or more model input parameters are thought to be significantly
uncertain, the Monte Carlo module (MONTEC) can be used, or a sensitivity analysis can
be performed by altering the input parameter values in a series of WHPA model runs and
overlaying the resulting capture zones.
-------
Problem Conceptualization vs. Reality Potential for Abuse i n_g
Figure 10.4
Barrier-Boundary As Represented by the WHPA Model (a), and
Barrier Boundary That May Be Encountered In Practice (b)
fir-*
- i - i -
i-i-i-i-i-i-i-
i-i-i-i-i-i-i-i
-i-i-i-i-i-i-i-
i-i-i-i-i-i-i-i
-i-i-i-i-i-i-i-
xxxxxxxxxxxxxxxxxx x x x x x x
Barrier
(a)
i - i - r=
-i-i-i-i-i-i
i-i-i-i-i-i-i
-i-i-i-i-i-i-
i-i-i-i-i-i-
XXXXXXXX
-------
References
Ballaron, P.B., 1988. Ground-Water Flow Model of the Corning Area, New York.
Publication No. 116 of the Susquehanna River Basin Commission, Harrisburg, PA.
Huyakom, P.S. and J.E. Buckley, 1989. SAFTMOD - Saturated Zone Flow and Transport
Two-Dimensional Finite Element Model. HydroGeoLogic, Inc., Herndon, VA.
Javandel, I., C. Doughty, and C.F. Tsang, 1984. Groundwater Transport: Handbook of
Mathematical Models. Water Resources Monograph 10, American Geophysical
Union, Washington, DC.
Javandel, I., and C.F. Tsang, 1986. Capture-Zone Type Curves: A Tool for Aquifer Clean
Up. Ground Water, v. 24, no. 5, pp. 616-625.
Lee, K.H.L. and J.L. Wilson, 1986. Pollution Capture Zones for Pumping Wells in Aquifers
with Ambient Flow. EOS, Transactions of the American Geophysical Union, v. 67,
p. 966.
McDonald, M.G. and A.W. Harbaugh, 1988. A Modular Three-Dimensional Finite
Difference Ground-Water Flow Model. Techniques of Water-Resources Investigations
of the USGS, Book 6 Chapter Al.
McGrath, E.J., and D.C. Irving, 1973. Techniques for Efficient Monte Carlo Simulation.
Technical report prepared for Department of the Navy, Office of Naval Research,
Arlington, VA.
Morrissey, D.J., 1987. Estimation of the Recharge Area Contributing Water to a Pumped
Well in a Glacial-Drift, River-Valley Aquifer. U.S. Geological Survey Open-File
Report 86-543.
Newsom, J.M., and J.L. Wilson, 1988. Flow of Ground Water to a Well Near a Stream -
Effect of Ambient Ground-Water Flow Direction. Ground Water, v. 26, no. 6, pp.
703-711.
-------
References
Pollock, D.W., 1988. Semianalytical Computation of Path Lines for Finite-Difference
Models. Ground Water, v. 26, no.6, pp. 743-750.
Schooley, J., 1989. Application of a Geographic Information System in New Jersey's
Wellhead Protection Program. Proceedings of the 9th Annual ESRI User Conference,
Palm Springs, CA
Seattle Water Department, 1986. Proposal for Groundwater Recharge Demonstration
Program at Highline Well Field. Prepared by CH2M Hill, Hart-Crowser & Associates,
Economic and Engineering Services Inc. and Herrera Environmental Consultants.
Seattle, WA.
Smith, L. 1981. Spatial Variability of Flow Parameters in a Stratified Sand. Mathematical
Geology v. 13, no 1, pp. 1-21.
Todd, D.K 1980. Groundwater Hydrology. John Wiley & Sons, New York, NY.
U.S. EPA Office of Ground-Water Protection, 1987. Guidelines for Delineation of
Wellhead Protection Areas. Washington, DC. EPA 440/6-87-010.
van der Heijde, P.K.M., and M.S. Beljin, 1988. Model Assessment for Delineating Wellhead
Protection Areas. Final report for Office of Ground-Water Protection, U.S.
Environmental Protection Agency, Washington, DC. EPA 440/6-88-002.
-------
Appendix A
Theoretical Development and
Implementation of RESSQC Module
A.I Introduction
In this appendix, analytical solutions used in the RESSQC module are presented.
These solutions are derived from simple applications of the complex potential theory of
fluid dynamics. The analytical solutions implemented in RESSQC may be constructed
using three fundamental component solutions. These components are: (1) the uniform
flow solution (2) the point source/sink solution; and (3) the finite radius source/sink in
uniform flow solution. The various components can be combined or superposed to derive
a more complex overall solution of the ground-water flow system. The basic assumptions
used in developing the analytical solutions in RESSQC are as follows:
The aquifer is homogeneous, isotropic, and of constant saturated thickness.
The flow of ground water in the aquifer is two-dimensional in a horizontal plane
and is at a steady state.
The fundamentals of complex potential theory as applied to ground-water flow
problems are outlined below. The notation used here is consistent with that used by
Javandel et al. (1984) in their presentation of the original RESSQ code.
Let 0 and $ denote two functions known as the velocity potential function and stream
function, respectively. These functions satisfy the following pair of Laplace's equations:
32
+ = 0 (A.I)
-------
A-2 RESSQC Formulation
+ = 0 (A.2)
ax2 ay2
where x and y are Cartesian coordinates in a horizontal areal plane. Equations (Al) and
(A.2) can be used to describe steady-state two-dimensional ground-water flow if the
following definitions are adopted:
d + i\l/ (A. 5)
where i
With the above definition, the complex potential function, W, can be used to provide
simple representations of various two-dimensional flow patterns. The real part of W,
which is 0, can be used to produce equipotential lines (or lines of constant hydraulic head
values) in the flow field. Using Darcy' s law, the velocity potential, 0, and the hydraulic
head, h, are related by 0 = Kh. On the other hand, the imaginary component of W,
which is 0, can be used to produce streamlines (or flow lines). Each streamline
corresponds to a constant value of 0.
-------
RESSQC Formulation
A.2 Analytical Solutions Used in RESSQC
A.2.1 Uniform Flow
A uniform ambient ground-water flow in a direction making an angle a with the
positive x-axis (Figure Ala) may be represented by
w = -UZ(cosa-isina) (A.6)
where U is the magnitude of the Darcy velocity, and Z is a complex variable defined as
z - x + iy (A.7)
Combining Equations (A.6), (A.5), and (A.7) yields
+ i\l/ = -U (cosa-isina) (A. 8)
from which the following expressions for and $ are obtained:
0 - -U (A.9)
^ = U (A. 10)
Components of the Darcy velocity are obtained by differentiation of (A.9). Thus,
d
-------
A-4
RESSQC Formulation
t
Figure A.I
Fundamental Flow Fields.
Uniform Flow Field (a), and
Radial Flow Fields (b)
X
(a)
Sink
X
Source
-------
RESSQC Formulation A-5
W = - in (Z-Z0)
2ffb
(A. 13)
where
Z0 = x0 + iy0
(A. 14)
Substituting for complex numbers Z and Z0 in (A. 13), the velocity potential and stream
function for the specified source may be obtained as follows:
W = +
0 + ill/ =
2TTb
tn (Z-Z0)
,
in [(x-x0)2 + (y-y0)2]
(A. 15)
Thus
Q
-
4rrb
£n [(x-x0)
(y-y0)2]
(A. 17)
tan"1
The Darcy velocity components, q* and q_y, are given by
(x-x0)
(y-y0)
(A.19)
27rt>
+ (y-y0)
(A.20)
The above analytical solution describes the radial flow field emanating from the source
located at fa, y0). Replacement of Q by -Q has the effect of reversing the flow direction.
The resulting flow field then corresponds to that of a converging radial flow field toward
a point sink located at (x,,, y0).
-------
A-6 RESSQC Formulation
A.2.3 Finite Radius Source in Uniform Flow
In a situation where a source of contamination covers a large area, representing it
by a point in the two-dimensional flow field may be inaccurate. A more realistic
representation of this case is to treat it as a circular source whose radius is computed such
that the area of the circle equals the area of the actual contaminant source. Physically,
the idealized circular source may represent either a surface impoundment a waste pond
or a lagoon that discharges leachate into the ground-water system.
Consider a particular case where such a source operates at a constant head in an
aquifer with ambient flow in the positive x-direction. The center of the source is located
at the origin of the coordinate axes, and the source is assumed to be a fully penetrating
source of radius I0. The complex velocity potential describing the resulting flow field may
be expressed as
Ur? Z Qp
W = -UZ + --- £n Z + C (A. 21)
|z|2 2nrb
where Qp is the rate of outflow from the finite-radius source, z and z are the conjugate
and the modulus of Z, respectively, C is a constant, and the remaining symbols are as
defined previously.
The expressions for potential and stream functions, 0 and $, can be obtained by
separating the real and imaginary parts of W. The resulting expressions are given by
0 = -Ux + - - - £n (x2+y2) + Cl (A. 22)
x2 + y2 47rb
n
tan
where c2 and c2 are constants. C2 may be set equal to zero and Cj is determined by
imposing that # satisfies the boundary condition of $ = H^ at the boundary of the source
where r = r0. Thus, c7 is evaluated as being equal to HO + ^ to* . The expression for 0
then becomes
-------
RESSQC Formulation
A-7
Uric
x2 + y2 4Trb
x2+y2
(A.24)
The Darcy velocity components, q, and q^, are obtained by differentiating (A.24) and
reversing the signs:
U H
n*.a
ur.2
r 9 9 "i
x2 - y2
*\
2xy
L (x2+y2) ~ j
2wb x2 + y2
QP y
27rb x2 + y2
(A.25)
(A.26)
A.2.4 Combination of Uniform Flow with Point Sources and Sinks
In a general situation where there are groups of discharging and recharging wells
operating simultaneously in an aquifer with ambient ground-water flow, the resulting flow
field can be constructed by means of superposition. The overall complex velocity potential,
W, of such a system is given by adding the uniform flow component to the source and sink
components. The general expression of W may be written as
N
W = -UZ (cosa-isina) + Z
in (Z-Zj)
M Qk
Z
k=l 2rrb
£n (Z-Zk)
(A.27)
where Qy is the flow rate of discharging well j (located at Z=Z^), Qk is the flow rate of
recharging well k (located at Z=ZJfc), N and M are the numbers of discharging and
recharging wells, respectively, and the remaining symbols are as defied previously.
-------
A-8 RESSQC Formulation
A before, the expressions for potential and stream functions, 0 and ft can be
obtained by separating the real and imaginary parts of W. These expressions are given
by
N
0 = -U (xcoso+ysina) +
M Qk
- s
k=l
_ =1 4rrb
£n [(x-xk)2 + (y-yk)2]
£n [(x-Xj)2 + (y-yj)2]
(A.28)
N Q,
U (xcosa+ycosa) + S tan
-i
x-x
J J
M
S
k=l
tan
"1
x-xk
(A.29)
The Darcy velocity components, q* and q^ are obtained by differentiating (A.28):
d N
= =s Ucosot Z
(x-x.,)2 + (y-y.,)
M
2
k=l
(x-xk)
(x-xk)2 + (y-yk):
(A.30)
<90 N Q.,
- = Usina - S
dy j=l
(x-x.,)
M Qk (y-yk)
+ s
k=i 2?rb (x-xk)2 + (y-yk):
(A.31)
Additional composite potential fields can be generated through various combinations
of the fundamental component solutions given in the previous sections.
-------
RESSQC Formulation
A.3 Pathline and Capture Zone Delineation Algorithms
The RESSQC module perform; the computation of pathlines by tracing the
movement of individual fluid particles through the ground-water flow system. (Because
steady-state flow is assumed, the term pathline and streamline can be used interchangeably.
The term "streamline" is normally defined as a line that is tangent at a given instant to the
velocity vector at each point upon it.) By tracing multiple streamlines, contaminant
frontal positions and capture zones may be delineated. The computational procedures are
described in the following subsections.
A.3.1 Pathline Computation
For a given particle, the pathline computation begins with the evaluation of particle
velocity at the starting position. For water particles, the seepage (or average pore water)
velocity components in the x and y directions are evaluated using the following formulas:
vx = q,/* (A.32)
vy = qy/0 (A.33)
where Y,. and vy are the seepage velocity components,
-------
A-10 RESSQC Formulation
where pB is the bulk density of the porous medium (aquifer material), and kd is the
distribution coefficient.
Consider the situation involving a water particle moving through the flow system.
The pathline of the particle can be traced by solving the following pair of differential
equations:
dx
dt = v* = V* (A- 37a)
dy
= Vy = qy/0 (A.37b)
Subject to the initial conditions:
x(t=0) = xc (A. 38a)
y(t=0) = yc (A. 38b)
where (xoi y0) is the initial location (starting position) of the particle.
The solution of the above equations is obtained through numerical integration. In
RESSQC, the Euler integration scheme is used to determine a new position of the particle
for each specified time increment, At. According to this scheme, given that at time t the
particle is at the point (^, y;-), then its new position at time t+At is determined from
xj+1 = Xj + (q^/0) At (A.39a)
Yj+i = Yj + (qy/0) At (A.39b)
Thus, starting from the initial position of the particle at t=O, one can use equations
(A.39a) and (A39b) repetitively to determine successive locations of the particle at later
time values. The accuracy of such a forward integration procedure depends on the size
of the time step, At.
In RESSQC an automatic control of the time step value is made by limiting Ai, the
incremental displacement length, along the pathline. The value of A* is determined in the
calculation process such that the following criteria are satisfied:
-------
RESSQC Formulation
1. At must not exceed the maximum prescribed step length input by the user.
2. The velocity direction change over the displacement length must not be greater
than a prescribed tolerance. More precisely, the magnitude of angular
acceleration must not exceed 1. When this criterion is violated, A« is reduced
by 1/2. The angular acceleration, w, is computed from
+ [(vyN/KI) - (vyp/|vp|)32 (A. 40)
where y^ and v^ are the seepage velocity components at the new point, y^,
and vyP are the seepage velocity components at the previous point, and | VN |
and vp are the magnitude of seepage velocity at the new and previous points,
respectively.
3. The selected value of A 2 must satisfy a velocity convergence criterion within 5
iterations. The iteration is done by checking the average of the velocities at
the previous and tentative points on the pathline against the velocity used to
compute the tentative point. If this differs by more than 1 percent of the
magnitude of the averaged velocity, then that average is used to compute a new
tentative point and the iteration is repeated. If convergence is not achieved
within 5 iterations, the value of A* is reduced by 1/2.
Using the algorithm described above, the code draws the flow pattern in the aquifer
by tracing the pathlines that emanate from specified starting locations. Each starting
location may correspond to the location of an injection well, a production well, or simply
a point from which a single pathline (or streamline) is to be traced. In the case of an
injection well, the flow rate and radius of the well need to be specified in conjunction with
the number of pathlines to be calculated for the well. If the flow rate is not specified, the
code assumes a zero default value and only one pathline is traced from the specified
starting point. Production wells are distinguished from injection wells via different input
specifications. In the original RESSQ code presented by Javandel et al. (1984), neither
-------
A-12 RESSQC Formulation
pathline nor capture zone delineation is performed for production wells. This limitation
is removed in RESSQC, which is the modified code presented in this report. In RESSQC,
pathlines and capture zone boundaries can also be delineated for production wells by
reversing the velocity field before particle tracking is performed.
A.3.2 Front Tracking and Capture Zone Delineation
Tracking of contaminant fronts surrounding an injection well (or an areal contaminant
source modeled as a circle of finite radius) can be made by keeping record of the travel
times of contaminant particles released from the boundary of the well (or circular source).
RESSQC requires the user to specify as input the travel time values at which the fronts
are to be calculated. At each prescribed travel time value, the terminal points of all
pathlines that originate from the contaminant source are identified. The frontal position
can then be drawn by joining these terminal points in a sequential order. Time-related
capture zones for production wells are also delineated by RESSQC using a similar
approach. In this case, reverse pathlines of water particles are traced from the well bore
of each production well. The travel times taken by these particles to move along their
pathlines are recorded. At each prescribed capture zone time, the terminal points of all
pathlines are identified. The capture zone boundary is drawn by joining these terminal
points sequentially.
A.4 Code Structure and Dimension Limits of RESSQC Module
Table A.I is a list of program segments (subroutines) that compose the RESSQC
module. The dimension limits of RESSQC arrays are listed in Table A.2.
A.5 Differences Between RESSQC and the Original RESSQ Code
Three major changes were incorporated into the original RESSQ code of Javandel
et al. 1984:
1. The code was modified to allow the delineation of capture zones as well as
contaminant fronts. If reverse tracking is to be performed, the velocity
components v, and vy are multiplied by -1. Use of the adsorption coefficient
is not permitted with this new option.
-------
RESSQC Formulation A-13
Table A.1
Description of RESSQC Subroutines
(Subroutines Listed In Alphabetical Order)
Subroutine Name
Description
FLOW Traces the path of a streamline. Streamlines terminate at pumping wells
(contaminant fronts), at injection wells (capture zones), or at the model
boundary.
PRINTWL Prints the well information to output file.
READW Reads the input parameters associated with each recharge/ pumping well
from the input file.
RESSQ2 Continuation of main program. Reads some input parameters from file,
initializes variables, calls subroutines to calculate streamlines, and writes
results to output file.
SORT1 Sorts the values of the time-related capture zones or contaminant fronts
into ascending order.
SORT2 Sorts the array TARR and NWARR first by the number of the well of
arrival, then by the number of the well of departure, and finally by the
time of arrival.
TESTAR Checks if the current particle location is within distance DL (maximum
allowable step length) of a production well (contaminant fronts) or
injection well (capture zones).
VXVY Computes the ground-water flow velocity vector at any given point (x,y).
WEED Weeds out every other point in the trace of a streamline starting with
the fourth point in order to obtain enough room in the XL and YL
arrays to save the entire streamline.
-------
A-14 RESSQC Formulatiom
Table A.2
Dimension Limits For RESSQC Arrays
Array Name
BETA
DATE
INDW
ITRW
NAMEW
NPATH
NWARR
QW
RADW
TARR
XI?'
XL
XW
YL
Dimension
140
7
140
140
140
140
140
140
140
140
(2,140,140)
7,000
(2,140)
7,000
Description
Angle of first pathline to leave the well.
Contaminant front or capture zone time values.
Contains flag for each well to indicate whether
or not contaminant fronts or capture zones
should be plotted.
Ratio of number of pathlines to number plotted.
Character array containing the name of each
pumping and injection well.
Number of pathlines to be computed to
delineate contaminant fronts or capture zones.
Well of arrival of pathline.
Flow rate into/from the well(s).
Well radii.
Time of arrival of pathline in pumping well
(contaminant fronts) or injection well (capture
zones).
Coordinates of contaminant fronts or capture
zone boundaries.
X-coordinates of pathline.
Well coordinates.
Y-coordinates of pathline.
-------
RESSQC Formulation A-15
2 . The code was "undirnensionalized". The units for all model input parameters
used by RESSQC must be consistent (e.g. all length values might be in feet,
and all time values might be in days). RESSQ has internal conversion factors
so that parameters must have a "practical" or "COS" system of units.
3. The common block structure of RESSQ was changed in RESSQC so that the
storage space for CHARACTER INTEGER, and REAL data is separated and
Fortran EQUIVALENCE statements are no longer used.
In addition to these major changes, several smaller changes were incorporated so that
RESSQC would conform to the overall WHPA model structure. For example, the plotting
output for pathlines is written to the file WHPA.PLT instead of TAPE7, and several
format statements have been altered.
-------
Appendix B
Theoretical Development and
Implementation of MWCAP Module
B.I Introduction
The MWCAP module is designed especially for the efficient delineation of capture
zones of individual production wells. Three types of capture zones can be delineated.
These include (1) steady-state capture zones (open-ended capture zones that correspond
to infinite travel time), (2) time-related capture zones (closed capture zones that
correspond to specified values of travel time), and (3) hybrid capture zones (combined
capture zones formed by closing a steady-state capture zone with a time-related upstream
boundary). MWCAP uses potential function analytical solutions to determine the locations
of stagnation points. The stagnation points are used as control points in delineating
capture-zone boundaries. The computational procedure is based on the fundamental
analytical assumption of steady-state two-dimensional flow in a homogeneous aquifer of
constant saturated thickness. In addition, MWCAP assumes that well interferences due
to simultaneous pumping have negligible effects on capture zones of the individual wells.
Three major cases are considered by MWCAP. These cases concern (1) an aquifer
without a lateral boundary (infinite areal-extent aquifer), (2) an aquifer with a stream
boundary, and (3) an aquifer with a barrier boundary. Much of the material presented
here follows the work of Newsom and Wilson (1988) and Lee and Wilson (1986). It
should be noted, however, that the present notation differs somewhat from the previous
works in terms of the definitions of potential and stream functions and the convention for
identifying the ambient flow angle.
-------
B-2 MWCAP Formulation
B.2 Analytical Solutions Used in MWCAP
B.2.1 Aquifer Without Lateral Boundary
For this simple case, the expression of the potential function, 0, for a pumping well
in uniform ambient flow may be written as
Q
-U(xcosa + ysina) + - £n [(x-x0)2 + (y-y0)2] (B.I)
where U is the Darcy velocity of the ambient flow, a is the angle of ambient flow
measured counter-clockwise from the positive x-axis, Q is the well discharge rate, b is the
aquifer thickness and ^, and y0 are the x and y coordinates of the well respectively.
If the coordinate system is oriented such that the origin is at the well and the x-axis
is parallel to the flow direction, then (B.I) reduces to
Q
0 = -Ux + £n [x2+y2] (B.2)
Using equation (B.2), the following expressions for Darcy velocity components are
obtained:
Q
= u -
(B'4)
From equations (B.3) and (B.4), it follows that there is one stagnation point (where q* =
q^ = 0), and its coordinates (x^y,) are given by
xa = Q/(2wbU) (B.5a)
y. = o (B.5b)
-------
MWCAP Formulation
B-3
Two streamlines pass through this stagnation point. These are the dividing streamlines
which separate the capture zone of the well from the rest of the aquifer. Both dividing
streamlines can be traced by releasing two particles in the neighborhood of each stagnation
point. The coordinates of these neighborhood points, points A and B, are set such that
(B.6a)
and
(xB,yB) =
(B.6b)
where is a small positive number. This approach is used by MWCAP to locate the two
bounding streamlines representing the boundary of the steady-state capture zone. An
analytical expression for the two dividing streamlines can be obtained as follows.
First the expression for stream function for this case is written as
-Uy +
Q
tan
MS]
(B.7)
which can be rearranged in the form
2?rb
2?rbU
y + tan
-i
y
X
(B.8)
where riD is a dimensionless stream function. The values of r)D corresponding to the two
dividing streamlines are K and-ff, respectively. Substituting these values into (B.8) gives
the equations for these streamlines. Thus one obtains
y - ±
Q
tan
-i
2bU 2ffbU
y
x
(B.9)
'
x -> oo, --> 0 and tan'O)
-------
B-4 MWCAP Formulation
Equation (B.9) indicates that the water-dividing streamlines approach asymptotically (i.e.
as x -> ») the lines y = ± Q/2bU.
B.2.2 Aquifer With Stream Boundary
Although the previous analytical solution is fairly simple, it is limited to the situation
where the well is far enough removed from any physical boundaries so that the effects of
these boundaries on the capture zone are negligible. A more general solution that
accounts for the presence of a nearby stream boundary is also available and implemented
into MWCAP. For this case, the potential function, $, is obtained using image well theory.
Consider the flow system shown in Figure B.I. To satisfy the stream boundary
condition (zero drawdown along the stream), an image pumping well is introduced at the
location (x,y) = (-d,0). The resulting expression of potential function 0 is given by
-U(xcosa + ysina) +
47Tb
(x-d)2 + y2
(x+d)2 + y2
(B.10)
Differentiating (B.10) gives the following expressions for Darcy velocity components:
d Q dF
= Ucosoc +
dx 4ffb ax
(B.ll)
qy = -
ay
Q 3F
Usina +
4rrb dy
(B.12)
where
ax
2 (x-d)
2 (x+d)
(x-d)2 + y2 (x+d)2 + y2
(B.13)
ar
ay
2y
. (x-d)2 + y2
2y
(x+d)2 + y2
(B.14)
-------
MWCAP Formulation B-5
Figure B.I
Schematic- Flow System For MWCAP Stream Boundary Formulation
IMAGE
WELL
STREAM BOUNDARY
AMBIENT
FLOW
-------
B-6 MWCAP Formulation
A general expression for stagnation-point coordinates can be derived by differentiating the
complex velocity potential, W, (referred to in Appendix A) and setting it to zero. For the
particular case considered here, the expression of W is given by
Q I Z-d
W = -U(cosa - isina)Z + £n
Z+d
(B.15)
where Z = x + iy, and i = j-\. Since W is an analytic function its derivative is given by
dW Q Q
= -U(coso - isinot) + - (B.16)
dZ 27rb(Z-d) 2wb(Z+d)
Setting dW/dZ = 0 at Z = Z, and solving for Z, gives
Z8 = ± d [1 + /9(cosa + isina)]^ (B.I 7)
where ft is a dimensionless parameter defined as
Q
and
Z8 - xs + iys (B.19)
where (x^y,) are stagnation point coordinates.
For the special case where the ambient ground-water flow is parallel to the negative
x-axis (a = 180°), and equation (B.I7) reduces to
Zs = ± d [1 - £]* (B-2°)
which, upon combination with (B.19), leads to three following possibilities (see Figure B.2).
(1) When ft < 1, only one stagnation point exists on the x-axis and its x-coordinate
is given by
xs = d [1 - 0]* (B.21)
-------
Figure B.2
Typical Pumping Well Capture Zones for a Well Near A Stream in an Aquifer With Ambient Flow:
a) Pumping Rate Below Critical Value; b) Pumping Rate at Critical Value; c) Pumping Rate Above Critical Value
(from Lee, 1986)
a) Small QQ0
STREAM
U
STREAM
U
'*S*&v*&pv33^^
^mM^mrn
u
d -H
LEGEND
CAPTURE ZONE
SWEPT ZONE
-------
B-8 MWCAP Formulation
Note that the other (negative) solution for x, is unrealistic because it is outside the flow
domain. This case of ft < 1 corresponds to Figure B.2a. As can be seen, the capture
zone does not intercept the stream because the pumping rate Q is smaller than the critical
value Qc (Qc = ffdtJb). Thus the stream does not contribute recharge to the well.
(2) When 0 = 1, there is also one stagnation point but it is located at the origin
of the coordinate axis. This value of ft is the critical value above which the stream will
begin to contribute to the flow of the well. The case of /3 = 1 corresponds to Figure
B.2b. As can be seen, the stagnation point contacts the stream when Q = Qc.
(3) When (I > 1, two stagnation points exist. They are both located on the y axis.
Their y-coordinates are given by
ys = ± d [ft - 1]* (B.22)
This case corresponds to Figure B.2c. As depicted, when the pumping rate Q exceeds the
critical value, Q., the stream contributes to the well flow. The shaded part of the capture
zone is called the stream capture zone; within this region the flow lines emanate from the
stream.
In a general situation where the direction of ambient flow is such that both cosa and
sina are nonzero, the stagnation point coordinates (x,,y,) must be determined from
Equation (B.I 7) by taking the complex root. It can be shown that the general expressions
for Xj and y, are given by
xs = ()' [7 ± (1 + /Scosa)]"2 (B.23a)
2
S = ± <)M7 ± (flcosa + 1)] (B.23b)
2
where
7 = (l-20cosa + 02)* (B.24)
-------
MWCAP Formulation
B-9
Physical consideration leads to the conclusion that the number of stagnation points cannot
exceed 2. Furthermore, the values of x, and y, that are acceptable must satisfy the
following constraints: (1) both x, and y, must be real and x,> 0, and (2) at (Xj,yJ the
velocity must be zero. These constraints are used to determine the correct solution for
x, and yf
B.23 Aquifer With Barrier Boundary
Another scenario that may be addressed using MWCAP is the case of an aquifer
with a barrier boundary (Figure B.3). For this case, the expression of the potential
function is obtained once again using image well theory. An image discharging well is
introduced at (x,y) = (-d,0). The resulting expression of the potential function, 0, is given
by
(f> = -U (xcosa + ysina) + (£n[(x-d)2 + y2]
4TTb
+ tn [(x+d)2 + y2]] (B.25)
Differentiation of Equation (B.25) gives the following expressions for Darcy velocity
components:
d0 Q dF
qx = = Ucosa +
dx 4trb dx
(B.26)
= - = usina +
where
Q dF
4?rb dy
(B.27)
dF
dx
2(x-d) 2(x+d)
+
(x-d)2 + y2 (x+d)2 + y2
(B.28)
2y
2y
(x-d)2 + y2 (x+d)2 + y2
(B.29)
-------
B-10
MWCAP Formulation
Figure B.3
Schematic Flow System For MWCAP Barrier Boundary Formulation
ON
IMAGE
WELL
BARRIER BOUNDARY
Q
AMBIENT
FLOW
-------
MWCAP Formulation B-ll
Note that the given potential function, 0, yields a nonzero velocity at the barrier boundary
(x=0). Indeed, the values of q* and qy correspond to the x and y components of ambient
flow velocity, respectively. This indicates that the barrier boundary simulated is not an
impermeable boundary but rather a semi-permeable boundary across which there is
ambient flow.
As before, a general expression for the stagnation point coordinates can be derived
by differentiating the complex velocity potential, W, which is given by
Q
W = -U(cosa - isina)Z + in [Z2-d2] (B.30)
2?rb
Differentiation of equation (B.30) gives
dW Q
= -U(cosa - isina) +
dZ
2Z
(B.31)
Z2 - d2
Setting dW/dZ = 0 at Z = Z, and solving for Z,, one obtains
/3dexp(ia ± [/92d2 exp (2ia) + 4d2]^
Z8 = (B.32)
2
where ft is defined by Equation (B.I8), and Z, = x, + iy,. From Equation (B.32), it can
be shown that the stagnation point coordinates are given by
ca = - [jSdcosa ± (7/2) * (1 +
2
ys = - [jSdsina ± (7/2)? (1 - cosO^] (B.34)
2
where
- (/32d2cos2a + 4d2) /7 (B.35)
-------
B-12 MWCAP Formulation
For the special case where the ambient flow is parallel to the negative x-axis (a = 180°),
Equations (B.33) and (B.34) reduce to
1 ,
xs - - [-/3d ± (/32d2 + 4d2)*] (B.36)
2
Y. - 0 (B.37)
Since x, must not be negative, there is only one feasible solution corresponding the
stagnation point located on the positive x-axis at
1 ,
xa = - [-^d + (^2d2 + 4d2)^] (B.38)
2
In a general situation where the direction of ambient flow is not necessarily parallel
to the x-axis, the feasible solutions for x, and y, can be determined from Equations (B.33)
and (B.34) using the set of constraints described previously for the case of an aquifer with
a stream boundary.
B.3 Capture Zone Delineation Algorithms
General algorithms for delineating three types of capture zones have been
incorporated into the MWCAP module. In a general situation involving multiple wells or
a well field, these algorithms are formulated in terms of local coordinates and applied to
the individual wells. The local coordinate system is set up for each well and related to
the global coordinate system of the overall problem by means of a coordinate transforma-
tion matrix. The procedure for setting up a local coordinate system is simply to orient the
local x-axis in the direction opposite to that of ambient ground-water flow. The local y-
axis is oriented at 90° and counterclockwise to the local x-axis.
B.3.1 Steady-State Capture Zones
The delineation of steady-state capture zones is performed in a straightforward
manner by performing pathline tracking of particles released from the neighborhood of
each stagnation point. Typical patterns of capture zones that can be expected for the
situation of an aquifer with a stream boundary and various ambient flow angles are
depicted in Figure B.4. Note that in case b of the figure where a = 240° and the capture
-------
MWCAP Formulation
B-13
Figure B.4
Typical Steady-state Capture Zones That Can Be Expected For The
Situation Of An Aquifer With A Stream Boundary And Various Angles
Of Ambient Flow (Newsom and Wilson, 1988).
LEGEND
Pumping WrtI
*A
*STAO
|"|y~,| Aquifer Capture Zon*
fjtfi-iH Slr**m Caplur* Zone
V///A Zont of Throughfloo
-------
B-14 MWCAP Formulation
zone boundary becomes tangential to the stream boundary, it is necessary to identify
point A (the critical tangent point). Once this point is located two additional particles
can be released from its neighborhood to trace the boundary segments of the capture
zone. "The location of A can be determined using the following expression of stream
function derived from (B.I5):
= Uxsina -Uycosa +
2trb
tan
-i
y
x-d
- tan"1
y
x+d
(B.39)
At point A, XA = 0 and Equation (B. 39) reduces to
y
ss -Uycosa -
]
(B.40}
Differentiation of Equation (B.40) with respect to y and setting the 'derivative to zero
yields the following equation for yA:
Q d
- Ucosa ( )
Trb d2 + y2.
A
which can be rearranged in the form:
(B.41)
(B.42)
coscr
where ft = QArdUb. To obtain a feasible solution from (B.42), we must have
COSOC
> 1
(B.43)
The algorithm for tracing a pathline from a given starting position is now presented.
This algorithm, which is used in MWCAP, is different from that used in the RESSQC
module. Refer to Figure B.5 as a guide for the following discussion.
Suppose the particle is currently located at point i. Its tentative (future) position is
first calculated from
VlTAt =
(B.44a)
-------
MWCAP Formulation B-15
Figure B.5
Iterative Predictor-Corrector Particle Tracking Procedure
Streamline
-------
B-16 MWCAP Formulation
= Yi + viyAt = Yl + Ayt (B.44b)
where v^ and v^ are seepage velocity components evaluated at point i. The initial value
of At is calculated from At = Aij^n/jvj, where A*,^ is the maximum prescribed spatial
step length and v | is the magnitude of the velocity at point i.
Next, a modified velocity vector is computed. Its components are determined from
= vx (Xi + AXJ2, Yi + Ayt/2) (B.45a)
- vy (Xi + AXJ2, y± + Ayt/2) (B.45b)
The modified velocity is then used to compute a corrected tentative position of the
particle.
* *
xi-n = xi + vi+%,x At = Xi + AXi (B.46a)
Yi+i = Yi + vi+%,y At = Yi + Ayi (B.46b)
The next step is to check if the corrected particle position is sufficiently close to the
previous predicted position. The checking is performed by evaluating the magnitude of
incremental displacement vector 3s (Figure B.5).
ds » As* - ASi (B.47)
If 13s | /1 As | is greater than a prescribed tolerance, then the value of At is reduced
by " and the calculation is repeated with the reduced time increment. The iteration is
continued until the error is within the specified tolerance.
B.3.2 Time-Related Capture Zones
A time-related (transient) capture zone is constructed by MWCAP using a procedure
similar to that used by RESSQC. In essence, a series of water particles are placed at
sequential locations along the perimeter of a small circle representing the well boundary.
These particles are then allowed to move along their pathlines until the assigned travel
time value is reached or until a physical or plotting boundary is encountered. The
terminal point of each particle is saved in two storage vector arrays. Finally, the capture
zone boundary is drawn by joining the particle terminal points in sequential order.
-------
MWCAP Formulation B-17
The algorithm in MWCAP is designed to circumvent the potential problem of an
overly truncated capture zone that may occur using the simpler scheme of RESSQC
(Figure B.6a). Note that in the RESSQC scheme, the initial particle locations are
uniformly spaced along the well boundary. Consequently, one may need to use an
excessive number of pathlines (streamlines) to obtain an accurate delineation of the
capture zone for a situation where the travel time is large or the shape of the capture
zone is elongated.
The improved scheme used in MWCAP is shown schematically in Figure B.6b; it
utilizes the control feature (the stagnation point) of the capture zone boundary. First, the
location of each stagnation point is identified, and then two controlling pathlines that make
sharp turns near the stagnation points are traced (Figure B.6b). After their sharp turn,
the points on these controlling pathlines are recorded. Coordinates of the selected points
and the stagnation point are stored in the vector arrays used to store coordinates of the
terminal points of the standard regularly spaced pathlines. When all the points are joined
in their sequential order, an accurate delineation of the transient capture' zone is obtained.
B.3.3 Hybrid Capture Zone
A hybrid capture zone is constructed by closing the steady-state capture zone with
a time-related capping boundary (Figure B.7). For a specified value of travel time, the
capping boundary is drawn by locating its control points (A, B and C) and passing a
quadratic curve through these control points. Control point B is located by performing
reverse particle tracking along the pathline emanating from the well and proceeding in the
upstream direction of ambient flow. The radial distance from the center of the well to
this control point is denoted by r. This radius is then used to locate the remaining control
points (A and C). The algorithm for determiningg the coordinates of A and C is presented
in the following paragraphs.
Since both points A and C are tentative points on the dividing streamlines, their
coordinates can be readily computed. For point A, one obtains (Figure B.7):
As
vix At = xt + v^ (B.48)
-------
B-18 MWCAP Formulation
Figure B.6
Capture Zone Boundary Truncation Error That May Occur Using RESSQC (a)
and Improved Algorithm Used by MWCAP (b)
TRUNCATED PORTION OF CAPTURE ZONE
\
(a)
SELECTED POINTS ON
CONTROL PATHLINE
CONTROL PATHLINE
NCONTROL PATHLINE
STAGNATION POINT
STREAM BOUNDARY
AMBIENT FLOW
(b)
-------
MWCAP Formulation
B-19
Figure B.7
Schematic-Description of the Procedure Used to Locate the Capping
Boundary of a Hybrid Capture Zone
WATER-DIVIDING
CAPTURE-ZONE
BOUNDARY
AMBieRT
FLOW
TIME RELATED
CAPPING BOUNDARY
-------
B-20 MWCAP Formulation
As
Y = Yi + vly At = Yi + viy (B.49)
I Vl
where v^ and v^ are seepage velocity components at point i (the previous point), vi is the
magnitude of the seepage velocity at i, and As is the incremental length along the
streamline (distance between points A and i).
Let x,* be defined as x,* = ^ - d, where d is the distance from the stream to the
pumping well. The prescribed radial distance from the well to point A may be expressed
as
r2 = [Xi + ( )]2 + [yi + ( ) As]2 (B.50)
which may be rearranged in the form
(As)2 + 2Ax [ - + - 3 ~ (r2 - r42) = 0 (B.51)
Vi Vi
where
rt2 = (x*.)2 + y2 (B.52)
Equation (B.51) can readily be solved for As. Since it is a quadratic equation, there are
two roots. However, the feasible root must be a positive value. Once As is determined,
the coordinates of point A can be obtained from equations (B.48) and (B.49), respectively.
The coordinates of point C can be determined in a similar manner using the same
calculation procedure.
B.4 Code Structure and Dimension Limits of MWCAP Module
Table B.I is a list of program segments (subroutines) that compose the MWCAP
module. A simplified flow chart for MWCAP is presented in Figure B.8. The dimension
limits of MWCAP arrays are listed in Table B.2.
-------
MWCAP Formulation
B-21
Table B.I
Description of Subroutines That Compose The MWCAP Module
(Subroutines Listed In Alphabetical Order)
Subroutine Name
Description
AKTRAK Traces steady-state capture zone boundaries using particle tracking.
Particles are released from the stagnation point(s).
BKTRAK Delineates streamlines using reverse particle tracking. Particles are
released at the pumping well and are adsorbed at the stagnation point(s)
or the aquifer boundary.
CHECK Checks the accuracy of the projected particle location iteratively. If the
error is too great, the time step (DT) is halved until the error is less
than the tolerance (TOL) or the maximum number of iterations
(NITMAX) is reached. Currently, TOL=0.05 and NITMAX=20.
CINTER Performs quadratic finite element interpolation of coordinates of hybrid
capture zone capping segment.
CZINTP Traces bounding streamlines for time-related capture zone boundary.
MWCAP Executes the MWCAP computational modules iteratively for 1=1,
NWELL where NWELL=number of wells for which capture zones are
to be delineated. Initializes all problem-specific variables.
SIAQXA Computes pumping well capture zones and swept zones using particle
tracking for an areally infinite aquifer or a semi-infinite aquifer bounded
by a gaining stream or a barrier boundary. Any angle of ambient flow
toward the stream or barrier boundary may be specified.
STAGPT Determines the coordinates of the stagnation points (XS,YS) and the
bounding point (YA) of a stream capture zone. Applicable to aquifers
with no boundary or a stream boundary only.
VSIAQ Computes particle velocity components (VX,VY) as a function of
location (X,Y) for an aquifer with no boundary, a stream boundary, or
a barrier boundary. Ambient flow is at angle ALPHA measured
counter-clockwise with respect to the negative x-axis.
XYTRAN Converts any set of local coordinates to global coordinates based on the
angle of rotation, ALPHA.
-------
B-22 MWCAP Formulation
Figure B.8
Simplified Flow Chart for MWCAP Module
C Start )
Input problem variables
using pre-processor
I
MWCAP: Initialize problem-
specific variables for
current well
N
BKTRAK: Delineate NSTRL
streamlines that define
time-related capture zone
AKTRAK: compute
steady-state and
hybrid capture
zones
Output results
to plot file
WHPA.PLT
N
-------
MWCAP Formulation
B-23
Table B.2
Dimension Limits For MWCAP Arrays
Array Name
Dimension
Description
ANGLE
DSW
GRAB
IBOUND
ICZTYP
NPSTKP
NSTLIN
FOR
QAMB
QWELL
STRANG
THETA
THICK
50
50
50
50
50
50
50
50
50
500
50
50
Angle of ambient flow (0-360')
Perpendicular distance from well to associated
barrier or stream boundary (=0 if no boundary
exists)
Ambient hydraulic gradient
Array indicating the existence or non-existence
of a boundary condition for a given well
Array that specifies capture zone type (steady-
state, time-related, or hybrid)
Array to save the number of points retained
for each streamline stored in the XSTLKP and
YSTLKP arrays
Number of streamlines to be delineated using
reverse particle tracking
Aquifer porosity
Ambient rate of ground-water flow (equal to
transmissivity * gradient)
Well discharge
Array of angles at which streamlines emanate
from the pumping well
Orientation of boundary condition (0-3607
Aquifer thickness
-------
B-24 MWCAP Formulation
Table B.2 (continued)
Dimension Limits For MWCAP Arrays
Array Name
TMCZ
IRAN
XCZ
XG
XL
XSTAG
XSTL
XSTLKP
XTFRNT
XTMCZ
XWELL
YCZ
YG
YL
Dimension
50
50
(500,6)
500
500
2
500
(500,4)
100
40
50
(500,6)
500
500
Description
Time value associated with time-related capture
zones
Transmissivity
x-coordinates for 6 possible capture zone
boundary segments
Global x-coordinates of current pathline
Local x-coordinates of current pathline
x-coordinates of stagnation points
x-coordinates of current streamline or capture
zone boundary segment
x-coordinates of special streamlines tracked
close to stagnation points for time-related
capture zone (maximum of 4 streamlines with
500 points per line)
x-coordinates defining the boundary of a time-
related capture zone
x-coordinates of capping segment for hybrid
capture zone
x- coordinate of pumping well
y-coordinates for 6 possible capture zone
boundary segments
Global y-coordinates of current pathline
Local y-coordinates of current pathline
-------
MWCAP Formulation
B-25
Table B.2 (continued)
Dimension Limits For MWCAP Arrays
Array Name
Dimension
Description
YSTAG
YSTL
YSTLKP
YTFRNT
YTMCZ
2
500
(500,4)
100
40
y-coordinates of stagnation points
y-coordinates of current pathline or capture
zone boundary segment
y-coordinates of special pathlines tracked close
to stagnation points for time-related capture
zone (maximum of 4 streamlines with 500
points per line)
y-coordinates defining the boundary of a time-
related capture zone
y-coordinates of capping segment for hybrid
capture zone
YWELL
50
y-coordinate of pumping well
-------
Appendix C
Theoretical Development and
Implementation of GPTRAC Module
C.I Introduction
GPTRAC is a general-purpose particle tracking module designed to supplement the
RESSQC and MWCAP modules. GPTRAC can handle site-specific conditions that may
invalidate certain analytical assumptions mentioned in Appendices A and B. The GPTRAC
module offers both semi-analytical and numerical options to perform time-related pathline
and capture zone delineation. The semi-analytical option is resewed for cases involving
homogeneous aquifers with relatively simple boundary conditions. It is designed to
accommodate various aquifer types including confined and leaky semi-confined aquifers as
well as unconfined aquifers subject to recharge or infiltration. Additionally, aquifers
bounded on two or three sides by barrier and stream boundaries can be handled. The
GPTRAC semi analytical approach also accounts for well interferences. The numerical
option is intended for more complex situations involving heterogeneous and anisotropic
aquifers subject to various boundary conditions. It utilizes velocity computational schemes
that require knowledge of a steady-state areal distribution of hydraulic (piezometric) head
in the aquifer. The required hydraulic head data can be determined using a numerical
(finite element or finite difference) model to solve the governing partial differential equation
of ground-water flow. Alternatively, the head data may be derived from field mapping of
potentiometric or water level data.
C.2 Theoretical Background
In GPTRAC, computation of ground-water velocities is performed using two alternative
approaches. The first approach is via analytical formulas derived from the potential function
analytical solutions of well flow problems. This approach is used in the semi-analytical
option for pathline delination. The second approach is via application of Darcy's law to the
hydraulic head distribution. This approach is used in the numerical option.
-------
C-2 GPTRAC Formulation
C.2.1 Analytical Solutions Used in GPTRAC
The well flow analytical solutions used in GPTRAC include those involving confined
aquifers described previously in Appendices A and B and three other solutions described in
this section. The additional analytical velocity computation formulas correspond to the
following situations:
well flow in an unconfined aquifer with areal recharge
well flow in a leaky semi-confined aquifer
well flow in strip aquifers with two parallel boundaries
C.2.1.1 Unconfined Aquifer with Recharge
The analytical solutions presented in Appendices A and B are intended for nonleaky
confined aquifer conditions. They may be inadequate for use in some unconfined flow
situations. An improved solution that accounts for the water-table and recharge conditions
of the unconfined system is presented in this section.
Consider a case involving a fully penetrating well pumping at a rate Q from an
isotropic unconfined aquifer replenished by a constant rate of accretion N reaching the water
table. The flow problem is nonlinear due to the variation of the aquifer saturated thickness.
To linearize the problem, the following definition of the potential function,*, is introduced:
* = Kh2 (Cl)
2
Where K is the hydraulic conductivity, and h is the hydraulic head above the base of the
aquifer. By assuming that the flow in the aquifer is essentially horizontal (i.e., using the
Dupuit assumptions), the governing equation for a steady-state case may be expressed in the
form:
d2* + d2* + N = 0 (C.2)
dx2 dy2
-------
GPTRAC Formulation £.3
Equation (C.2) may be solved for a general situation involving an areally infinite aquifer with
a finite recharge applied over a circular area bounded by 0< rK (C.3b)
where
r2 = (x-xj1 H- (y-yj2 (C.4a)
To ensure continuity of hydraulic gradient (i.e., zero gradient) at r=R, the following
condition is imposed:
^ = ^recharge
thus
Q = 7rR2N
yielding
Jt - (-2-f (C.5)
M
-------
C-4 GPTRAC Formulation
The analytical solution given in (C.3a) and (C.3b) does not incorporate ambient flow.
The effect of ambient flow may be taken into account in an approximate manner by
superposition. This leads to the following modification of the solution:
$ = -Ub (XCOSCL + ysinct) +
-^ [1-ln (Rlr)] + - (R2-r*) for r z R (C.6a)
2it 4
O = -Ub (xcosa + ysintt) + 0 for r>R (C.6b)
Where b is the ambient saturated thickness, assumed to correspond to H.
Once the distribution of the potential function is known, the hydraulic head and Darcy
velocity components can be determined from:
h = 2*/X(1/2> (c 7)
q
,-! fj*]
h ( obc J
(C.8a)
q = -K *- -I fi*) (C.8b)
ay h(dy)
For the case involving a semi-infinite unconfmed aquifer with a pumpint well and a
stream boundary, the problem can be solved by means of an imaginary recharging well. At
some point (x,y) in the physical domain, the contribution of the imaginary source to the
potential field can be evaluated using the fundamental solution given in (C.3a) and (C.3b).
Similarly, the case involving a semi-infinite unconfmed aquifer with a pumping well and
a barrier boundary can also be handled using an imaginary pumping well.
-------
GPTRAC Formulation £.5
C.2.1.2 Leaky Semi-Confined Aquifer
A leaky semi-confined aquifer is an aquifer that can lose or gain water through an
underlying or overlying formation. When pumping takes place in such an aquifer, the
resulting drawdown produces leakage into or out of the aquifer through the confining
(semipervious) layers.
Consider a situation involving a uniform leaky confined aquifer with an overlying
aquitard that is subject to constant head at the top boundary (Bear, 1979, p. 313). Assuming
that flow in the aquifer is horizontal and leakage in the aquitard is vertical, a steady-state
well flow potential solution ,$, for an areally infinite aquifer may be expressed as
where the potential function * is defined as § = Kh
$0 is a reference potential at r=<», Q is the well pumping rate, b is the aquifer
thickness, K0 is the modified Bessel function of the second kind and zero order, r is given
by equation (C.4a) and A is a leakage factor defined as
k2 = Kbb'lK' (c-10)
where K is the hydraulic conductivity of the aquifer, and K' and b' are the hydraulic
conductivity and thickness of the aquitard, respectively.
Equation (C.9) is a fundamental analytical solution that can be extended in the same
manner as described previously (Appendices A and B) to account for ambient flow and the
presence of a barrier or a stream boundary. The Darcy velocity components in the x and
y directions can be obtained by differentiating equation (C.9) with respect to x and y,
respectively.
-------
C-6
GPTRAC Formulation
C.2.1.3 Strip and Three-Sided Aquifers
A strip aquifer is defined as an aquifer bounded by two parallel straight boundaries of
infinite length. Four possible cases are of practical interest. The first and second cases
correspond to an aquifer with two barrier boundaries and an aquifer with two stream
boundaries. The third and fourth cases correspond to strip aquifers with a stream boundary
on one side and a barrier boundary on the other side.
The analytical solution to steady-state well flow in a strip aquifer may be derived using
the conformal mapping theory. For the case of strip aquifer with two barrier boundaries
located at y = 0 and y = a, the complex potential resulting from a well at (x0, y0) is given
by
W =
(In
(C.ll)
a
where Q is the flow rate (> O for pumping and <0 for recharge), b is the aquifer thickness
and Z0 is the complex conjugate of Z0. The velocity field can be obtained by differen-
tiating the complex potential.
Re =
sinho
sno
dx) 4ab lcosha-cosp_ cosha-cosp+
(c
and
,
dx
4ab I cosh a -cos P _ cosh a -cos P
where
a =
; p_ =
; p+ =
-------
GPTRAC Formulation
C-7
Similarly, a strip aquifer that is bounded by two stream boundaries at y=0 and y=a can be
mapped into a half-plane with a stream boundary. The complex potential for this case is
W = -₯- In
2-Kb
e1tzla-evz'la
e e
(C.14)
and the gradient of the potential is
Q
4ab
tf-COSp_ a-COSp+
cosha-cosp_
cosha-cosp_
(C.15)
Im
4ab
sinp.
cosh a -cos P _ cosh a -cos p _
(C16)
whOere a, ^8. and /3+ are defined previously.
Next, we consider a strip aquifer bounded by a barrier boundary at y=0 and a stream
boundary at y=a. For this case the complex potential can be found using the theory of
image wells. The complex potential is given by
2Kb
(CM)
-------
C-8 GPTRAC Formulation
and the gradient of the potential is
where
dW _ Q
dx Sab
*- cosp_ e*+ cosp+
cosha- cosp_ cosha + cosp+
e"+ cos
8-
cosha + cosp_ cosha- cosp+
(C.I 8)
By Bab
sinp_
1
1
cosha - cosp_ cosha + cosp_'l
-sinpt
1
1
cosha
cosha -
(C.19)
« =
2a
; P_
2a
n
; P+ =
2a
The solution for a strip domain with a barrier at y=0 and a stream at y= a can be
found from the previous case using the following coordinate transformation:
x = x'
y = -y' + a
(C.20)
-------
GPTRAC Formulation £.9
Thus, a solution at (x', y') in the z' plane can be obtained by using
x = x' x. - x. (c 21)
y = -y + a y. = -y. + a
in the solution (C.17)-(C.19). However, the sign of the y-velocity so obtained must be
reversed to obtain a physically valid solution.
The analytical solutions of the four cases of strip aquifers can be readily extended to
corresponding cases involving aquifers bounded on three sides by barrier boundaries, stream
boundaries, or combinations of barrier and stream boundaries. The presence of a third
boundary is handled in the usual manner by the use of image well theory. 'Treating the third
boundary as a mirror, a discharging or recharging well may be introduced depending on
whether the boundary is a barrier or stream boundary.
C.2.2 Darcy and Pathline Equations and Time Integration
Given an initial location (x,0, y0) of the fluid particle, its pathline can be determined
by integrating the following pair of differential equations:
= vx (C.22)
dt
= vy (C.23)
dt
where (x(t),y(t)) are coordinates of the pathline at time t.
In GPTRAC, the integration is performed using either an analytical or a numerical
scheme. The appropriate scheme is selected automatically by the code. The analytical
scheme, which yields an exact integration, is used wherever possible (i.e., when its limiting
analytical assumptions are satisfied).
-------
C-10 GPTRAC Formulation
C.2.2.1 Darcy's Equations
Given a distribution of hydraulic (or piezometric) head in the aquifer, the components
of Darcy velocity (specific flux) in the x and y directions may be obtained from
(C-24)
f
where K^ and K^ are the principal components of hydraulic conductivity along the x and
y directions, respectively, and h is the hydraulic head. The components of ground-water
seepage velocity are given by
(C.26)
vy= qyIQ (C.27)
where 9 is effective porosity of the aquifer.
Equations (C.26) - (C.27) are used in GPTRAC for the numerical computation of
ground-water velocities. The distribution of hydraulic head is supplied to the code by a data
file containing steady-state values of h at the nodal points of a rectangular grid. GPTRAC
can accommodate both mesh-centered and block-centered grids. These two types of grids
are illustrated in Figure C. 1. The nodes are located at the intersections of grid lines for the
mesh-centered grid, and at the centers of grid blocks for the block-centered grid.
C.2.2.2 Numerical Integration
The numerical integration scheme implemented into the code is an Euler predictor-
corrector scheme with successive iterations to accommodate spatial variability of velocities
and curvature of the flow path. This scheme is described in Appendix B.
-------
GPTRAC Formulation
Figure C.I
Rectangular Grids Used by GPTRAC Mesh-Centered Grid (a)
and Block-Centered Grid (b)
(a)
©
<5) = Pumping Well
No. of nodes=90
No. of elements=72
0 =Pumping Well
No. of nodes=72
No. of grid blocks=72
-------
C-12 GPTRAC Formulation
An alternative integration scheme that is applicable, but which is not used in
GPTRAC, is the fourth-order Runge-Kutta scheme. This scheme has been used without
iterative corrections in a particle tracking code called GWPATH (Shafer, 1987). It is based
on the following RUNge-Kutta formulas:
(C 28)
^ J
where (x(n),y(n)) are coordinates at the previous location of the particle, (x(n+l),y(n+l))
are coordinates at the new location, and
ax = At * f[x(n),y(n)]
bx - At * g[x(n),y(n)]
a2 = At * f[x(n)+a1/2,y(n)+b1/2]
b2 - At * g[x(n)+a1/2,y(n))b1/2]
a3 = At * f[x(n)+a2/2,y(n)+b2/2]
b3 = At * g[x(n)+a2/2,y(n))b2/2]
a4 - At * f[x(n)+a3/y(n)+b3]
b4 = At * g[x(n)+a3/y(n)+b3]
Equations (C.30) and (C.31) yield a higher order numerical approximation than the
Euler integration formulas. However, it should be noted that the accuraq of the pathline
delineation also depends on the manner in which time stepping is performed. Thus, the
Euler scheme with iterative corrections and good control of time steps may be more
accurate than the fourth-order Runge-Kutta scheme used without iterative corrections to
control the time step size.
C.2.2.3 Analytical Integration
If a steady-state hydraulic head distribution is obtained using a numerical model, an
efficient analytical integration scheme can be derived for those grid blocks (or rectangular
elements) that do not contain the well nodes. Such a scheme has been successfully applied
by Pollock (1988) to velocity fields obtained using the block-centered finite difference
method to solve the ground-water flow equation. The scheme can also be adapted to the
rectangular finite element solution of the ground-water equation. It is based on two key
assumptions: (1) the x-velocity component, v^ within a grid block or element is a linear
-------
GPTRAC Formulation £-13
function of the x-coordinate and is independent of the y-coordinate, and (2) the y-velocity
component, v« within the grid block or element is a linear function of the y-coordinate and
is independent of x-coordinate. Via these assumptions, the principal velocity components
inside the grid block can be expressed in the form (Figure C.2):
vx = Ax(x-x1)+vxl (C.SOa)
vy = Ay(y-y1)+vyl (C.SOb)
where Ax and Ay are constants that correspond to the components of the velocity gradient
within the cell,
AX - (vX2-vxl)/Ax (C.Sla)
Ay = (vy2-vyl)/Ay (C.Slb)
vxl and Vtf, are values of the x-velocity component at the edges of the block or element sides
that are parallel to they-direction, vyl and v^ are values of the y-velocity components that
are parallel to the x-direction, and AX and Ay are the dimensions of the block or element in
the respective coordinate directions.
For each grid block or element the representative (principal) values of the velocity
components, vxl, v,^, vyl and v^, are computed using a simple linear finite difference
approximation for the head values at the block centers or element centroids. The centroidal
values for finite element grid blocks are obtained by arithmetic averaging of head values at
the corner nodes. Interlock or interelement conductivity values are obtained by harmonic
averaging of individual block or element conductivities.
Linear interpolation of velocity components using equations (C.30) and (C.31) has two
desirable properties. First, a continuous velocity vector field within each individual grid
block (or element) that satisfies the differential equation of mass balance everywhere within
the block is produced. Second, an analytical solution of the pathline equations, (C.22) and
(C.23), can be derived. The first property can be shown by substituting equations (C.SOa),
(C.SOb), (C.Sla) and (C.Slb) into the following differential mass balance equation:
fc
-------
C-14 GPTRAC Formulation
Figure C.2
Principal Velocity Components of a Rectangular Grid Block or Element
(x2,y2)
Side y2
Side xl
Side z2
Side yl
-------
GPTRAC Formulation C-15
where b is the saturated thickness of aquifer, and W is the volumetric rate of recharge or
discharge per unit surface area of aquifer. Assuming that b is constant over the grid block,
equations (C.30a)-(C.31b) may be substituted into (C.32) to obtain
(C.33)
Ax Ay bAxAy
which is indeed the incremental form of the mass balance equation for a grid block or
individual rectangular finite element.
Note that equation (C.33) is a good representation of (C.32) provided that any internal
sources or sinks may be regarded as being uniformly distributed within the grid block or
element. This is usually a valid assumption as long as the grid block does not contain a well
node. For such a grid block, the analytical equations describing the pathline of a given
particle can be derived as described below.
Consider the situation depicted in Figure C.3, where a particle enters the block at point
(Xp,yp) and at time t . At a later time, t, the particle has moved along the pathline to point
(x,y), which may be determined by integrating equations (C.22) and (C.23) using (C.30a) and
(C.3Ob) to substitute for vx and v as follows:
X t
[ ^ = fdt = t-t = At (C.34a)
* Ax* + B* {
and
y j_ '
' : = Af (C.34b)
/-
where Bv and EL are defined as
A y
Bx = VXI-A^! (C.3 5 a)
and
By = Vyi-AyY! (C.35b)
-------
C-16 GPTRAC Formulation
Figure C.3
Schematic Representation of a Particle Path Through A Grid Block
VX1
(Xp ,Yp
X-distance travelled during
time iiiterval
Y- distance travelled during
time interval
-------
GPTRAC Formulation C-17
The integrals in (C.34a) and (C.34b) can be evaluated to give
(Axx+Bx) / (AxXp+Bx) ] = AxAt (C.36a)
and
By)] = AyAt (C.36b)
which may be arranged in the form
x = xx + ^_ [v exp(AxAt)-vxl] (C.37a)
Ax
Y = Yi + 1_ [Vypexp(AyAt)-vyl] (C.37b)
Ay
As shown in Figure C.3, the particle leaves the grid block at an exit point (xe,ye). The
coordinates of the exit point can be determined from (C.37a) and (C.376b), provided that
the time increment, Ate, required by the particle to travel from point p to point e is known.
The value of Ate can be determined as follows. First, consider the time increment,
-------
C-18 GPTR AT Formulation
At^ that would be required by the particle to travel from point p to side x2 of the grid block.
This is obtained using equation (C.36a), which may be rewritten as:
Atx - - ln[Axx2+Bx)/(Axxp+Bx)]
or
Atx = J- In[vx2/vxp] (C.38a)
similarly, the time increment that would be required by the particle to travel from point p
to side y2 of the grid block is given by
Aty = J- In[vy2/Vyp] (C.38b)
If A^ is less than Aty, the particle will exit at side x2. Conversely, if Aty is less than At^ the
particle will exit at side y2- A third possibility is that At^ = Aty, in which case the particle
will exit the grid block at the corner where sides x2 and y2 intersect. In any case, the general
expression for Ate is given by,
Ate = min(Atx,Aty) (C.39)
where A^ and Aty are determined from (C.38a) and (C.38b) respectively.
Knowing Ate, the coordinates (xe,ye) of the exit point can be determined from
Xe = Xl
(C.40a)
Ye = yl + AT [Vypexp(AyAte)-vyl] (C.40b)
Once the exit point is known, intermediate points of the pathline within the grid block are
easily determined by substituting values of At that are within the range 0 < At < Ate.
-------
GPTRAC Formulation C-19
The analytical procedure presented is based on the assumption that neither Axnor Ay
is zero (i.e. the velocity gradient of the grid block is non zero). In the case where Ax = 0,
the x-velocity component is constant in the grid block and A^ and xe are obtained from the
simplified relations
Atx = (x2-Xp)/vxl (C.41a)
and
(C.41b)
Similarly, if A« = 0, then Atand ye are given by
Aty - (Y2-yP)/vyi (C.42a)
and
Ye = Yp+vy]Ate (C.42b)
The preceding example focused on a specific velocity pattern where the sense of
direction of all velocity components is positive. Other possible scenarios must be considered
as well. For each coordinate, there are 4 possible velocity patterns. Figure C.4 illustrates
the possible patterns for the x-velocity components. Pattern 2 (IP ATX = 2) has been dealt
with. Pattern 1 (IP ATX = 1) is simply the reverse of pattern 2. In this case the particle
enters the grid block at side x2 and exits at side Xj. Pattern 3 (IP ATX = 3) is an interesting
one. In this case, a local flow divide exists where vx = 0. If a particle resides in the grid
block and its v is greater than 0, then it will exit the grid block at side x2. On the other
hand, if v is less than zero the particle will exit at side Xj. The fourth pattern (IP ATX =
4) also requires special attention. In this case the particle can not leave the grid block in
the x-direction. When this condition exists, a flag is set to indicate the possbility of particle
entrapment in the grid block. Particle entrapment will occur if a second flag is raised for
the velocity component in the y-direction. This indicates that a strong sink (pumping well)
is present within the grid block and that the particle is captured by the sink. The equivalent
4 patterns of velocity components in the y-direction are treated in the same manner as the
corresponding patterns in the x-direction.
-------
C-20 GPTRAC Formulation.
Figure C.4
Possible Patterns of the Principal Velocity Components Along the x-Direction
IPATX=1
(a)
V
X1.
IPATX=2
(b)
IPATX=3
a
(c)
IPATX=4
(d)
-------
GPTRAC Formulation C-21
C.3 Numerical Implementation
The delineation of water-particle pathlines and well-capture zones is performed in
GPTRAC using two approaches, referred to as the semi-analytical and the numerical
options. The semi-analytical option uses analytical formulas to compute ground-water
velocities at a series of points along each pathline. Each pathline is traced using the
numerical integration procedure described in Appendix B.
In the numerical option, ground-water velocities are determined numerically from a
steady-state head distribution that is supplied to the code via an input data file. The
computational procedure is then applied in three major steps: (1) generation of a
rectangular grid, (2) computation of the principal velocity components, vxl, v^, v j and v^
at sides xl, x2, yl and y2 of each grid block, and (3) tracking of particle pathlines using
analytical integration or numerical integration when analytical integration is inappropriate.
Each of these steps is described in detail in the next section.
C.3.1 Grid Generation
The generation of the rectangular grid required in the numerical option is performed
automatically by GPTRAC using simple control parameters specified by the user. The two
key parameters are NROWS and NCOLS, which denote the number of row grid lines
parallel to the x-axis, and the number column grid lines parallel to the y-axis, respectively.
Figure C.I illustrates a sample rectangular finite element grid and a sample block-centered
finite difference grid each with NROWS=9 and NCOLS= 10. Locations of discharging and
recharging wells are assumed to correspond to nodal locations.
In the case of the finite element grid (Figure C.la), the nodes correspond to points of
intersection of row- and column-grid lines, and each rectangular element is surrounded by
4 nodes. Thus the grid consists of 90 nodes and 72 elements, and a pumping well is
represented by an internal node shared by 4 elements. In the case of the block-centered
finite difference grid (Figure c.lb), the nodes are located at the centers of rectangular grid
blocks. Thus the grid consists of 72 nodes and 72 grid blocks, and a pumping well is
represented by a node associated with one grid block.
In addition to the above two cases, GPTRAC also admits the use of a mesh-centered
finite difference grid in determining the hydraulic head distribution. This case is treated by
GPTRAC in the same manner as the rectangular finite element grid.
-------
C-22 GPTRAC Formulation
In any case, the user must specify what type of grid has been used to determine the
hydraulic head data required for velocity computation (This is done through an input control
parameter, IHDATA Head values need to be input sequentially starting from the first
through the last nodes. For the sake of convenience, the code admits either a column-wise
or row-wise node numbering scheme to be used for the input of head data. Grid spacings
in both directions may be uniform or nonuniform. The spacings of a uniform grid are
determined automatically by the code via a specification of the minimum and maximum
values of the x and y coordinates. The spacings of a nonuniform grid are computed by the
code once the user has input the x-coordinates of the column grid lines and the y-
coordinates of the row grid lines.
The numerical option of GPTRAC can delineate pathlines in homogeneous or
inhomogeneous aquifers. If the aquifer is inhomogeneous, a zonal representation of the
material properties is used (Figure C.5). Each zone is represented by a rectangular
subdomain and assigned one set of hydraulic properties. A subdomain 'overlay technique
(block overlay) is used to identify material zonal numbers for the individual blocks in the
grid. The largest zone should always be defined first (most commonly the first zone will be
the entire modeled domain), and then the smaller zones are subsequently "overlayed" upon
the first zone. Portions of the grid for which there is no zone specified will be assigned
hydraulic property values used for the first zone. Each zone is identified by the coordinates
of its bottom-left and top-right comers. Using this information, GPTRAC identifies the
material zonal numbers for each element or grid block by checking its centroidal location.
C.3.2 Velocity Computation
The semi-analytical option for pathline delineation in GPTRAC computes the ground-
water seepage velocity components using analytical formulas based on the potential function
analytical solution of the steady-state well flow problem. This option is limited to cases
involving a homogeneous aquifer. However, the aquifer may be bounded on one or two
sides by a stream and/or barrier boundary condition. The available cases include: (1) wells
in an areally infinite aquifer, (2) wells in a semi-infinite aquifer with a stream boundary, and
(3) wells in a semi-infinite aquifer with a barrier boundary. Each boundary is assumed to
fully penetrate the aquifer. In each case, GPTRAC computes values of velocity components
at a series of points along each pathline. Each pathline is traced using the numerical
integration procedure described in section C.2.2.2.
-------
GPTRAC Formulation C-23
Figure C.5
Zonal Representation of an Inhomogeneous Aquifer
A
D
Zone no.
1
2
3
4
B
E
F
Bottom Left Corner
A
B
D
E
Top Right Corner
C
C
F
F
-------
C-24 GPTRAC. Formulation
In the numerical option for pathline delineation, the velocity computation is performed
using values of hydraulic head at the centers of grid blocks or the centroids of rectangular
finite elements. In the latter case, the head value at the element centroid is determined by
averaging the head values at the four comer nodes.
For a typical grid block or element I, the principal velocity components vxl, v^, v j and
Vy2 are computed as follows (Figure C.6a):
vxi = -KXIW(hw-hI)/AxIH (C.43a)
Vx2 " -KxIE(hE-hl)/AxIE (C.43b)
vyl = -KyIS(hs-hI)/AyIS (C.43c)
(C.43d)
where Axiw = xw~xi * AXIE = XE~XI -'
AYis = Ys-Yi / AYiN = YN-YI /
and Kjjw, K^, K,IS and KLjjj are harmonically averaged effective seepage conductivities
between blocks in the west, east, south and north directions respectively. K^^ and K^yp are
given by
KYTKYW(AxT+Axw) ,^..-.
(C.44a)
KXIAxw+KxWAxI
KVIE
KXIAXE+KXEAXI
where K^, KxW and K^ are seepage conductive for grid blocks I, W and E, respectively.
and K are defined in a similar manner.
If the grid used corresponds to a finite element grid, velocity component values are
also determined at the element centroids. The centroidal values of the velocity components
are desirable for those elements that share well nodes for which pathline computation is
performed using the numerical integration approach. The formulas for computing the
centroidal velocity components of a typical element e (Figure C.6b) are:
-------
GPTRAC Formulation £-25
Figure C.6
Notation Used Tin Velocity Computation For (a) Grid Block i, and (b) Element e
N
T
L
4
1
S
(a)
X
-------
C-26 GPTRAC Formulation
vx = -Kxe(h2+h3-h1-h4)/2Axe (C.45a)
vy = -Kye(h3+h4-h1-h2)/2Aye (C.45b)
Where K^ and Kye are element seepage conductivities in the x and y direction, respectively,
and Axe and Aye are the x and y dimensions of the element.
C.3.3 Pathline and Capture Zone Delineation Procedure
GPTRAC delineates the pathline of a given particle using either the forward or reverse
tracking method. Starting from the initial particle location, the code uses analytical or
numerical integration to determine successive locations of the particle at increasing time
values. The pathline computation continues until one of several termination critera is met.
These criteria include: (1) checking to determine if the assigned value of travel time has
been exceeded, (2) checking for boundaries of the flow region or plotting boundary limits,
(3) checking if a stagnation point has been encountered, and (4) checking for particle
entrapment by a well.
When the particle enters a grid block or element having a well node, a series of checks
are made to determine if the particle will be intercepted by the well. Either pumping or
injection wells are considered in the termination criterion depending upon whether forward
or reverse tracking is being performed.
If a block-centered finite difference grid is used, the pattern of principal velocity
components of the well-node block is first identified by checking values of IP ATX and
IPATY, referred to in Figure C.4 and Section C.2.2.3. If both IP ATX and IPATY are equal
to 4, then the particle is entrapped in the grid block and its terminal point is at the well
node location.
If a rectangular finite element grid is used, the algorithm depicted in Figure C.7
provides an efficient way to determine well interception of a fluid particle. GPTRAC
employs this algorithm in conjunction with a numerical integration and finite element
interpolation scheme. During any particular time step, the distance from the previous
location, PJJ.J, to the current location, Pn is compared with the distance from the same
previous location to the next location Pn+ j (Figure C.7). If the distance from P^ to Pn+1
is less than the distance from Pn+ j to Pn, a flag is raised to indicate that a turning or
reversal in the direction of flow of 90 degrees or more has probably occurred. Next, a
-------
GPTRAC Formulation
C-27
Figure C.7
Algorithm To Determine If A Particle Has Been Intercepted
By A Nearby Well, (a) Raise Flag and (b) Continue
abs(p_1 p^ ) < abs(p^( pm) ; raise nag
(a)
-------
C-28 GPTRAC Formulation
second check is made to determine if the element where the particle currently resides shares
a well node with other neighboring elements. If so, identification of the well is made to
ensure that it is correct type (pumping well for forward tracking and injection well for
reverse tracking). If this check is positive, iterative projections of particle positions continue
with successive reductions of time step (1/2 reduction each time) until the tentative position
of the particle is sufficiently close to the well. At this point, the terminal position of the
particle pathline is taken to be the well location, and the pathline computation for that
particle is completed.
If pumping wells exist, GPTRAC allows the user to select an option to delineate the
capture zone of any individual well using the reverse pathline delineation technique. First,
an effective radius of the well node (Rw) is estimated. If the velocity field is determined
analytically, the following formula is used:
1^ = 0.0025*Min(XL,YL) (C.46)
where XL and YL are the x and y plotting dimensions of the flow domain.
If the velocity field is determined using gridded head data, the following formula is
used:
(C.47)
where Aj is the effective grid area covered by node I, assumed to be the well node, and f
is a factor set equal to 1 if a block-centered finite difference grid is used and equal to 0.5
if a rectangular finite element grid is used.
Having determined the effective radius, the next step is to place a series of particles
along the circular boundary of the pumping well node and perform reverse pathline tracking
for each particle (Figure C.8). The number of particles released may be specified by the
user - the default number is 20. The reverse pathline tracking starts at an initial time value,
to, and finishes at the specified value of travel time, T. The initial time value corresponds
to the time taken by the fluid particle to move from its location at the circular boundary to
the well node. Assuming radial flow within the circular region, the value of to is determined
from
t0 = Trtb/Q (C.48)
where Q is the pumping rate of the well, and b is the aquifer thickness.
-------
GPTRAC Formulation C-29
Figure C.8
Capture Zone Delineation By Reverse Pathline Tracking
LEGEND
WELL LOCATION
REVERSE PATH LINE
INITIAL PARTICLE LOCATION
> TERMINAL POINT
-------
C-30 GPTRAC Formulation
C.3.4 Finite Element Interpolation of Velocities
When a finite element grid is used, pathlines passing through elements that share a
well node are computed by numerical integration. The velocity components inside such
elements are determined using a simple one-dimensional, quadratic interpolation scheme.
The interpolation formulas for this scheme may be expressed as follows:
vx = vxlNi(^) + vx2N2(O + vxN3(O (C.49a)
y yi- + vy2N2(rj) + VyN3(r?) (C.49b)
v =
where vxl, v^, vyl and v^ are the principal, side velocity components of the element, and
vx and vy are the centroidal velocity components, Nl5 N2 and N3 are one-dimensional,
quadratic finite element basis functions, and f and 77 are dimensionless isoparametric
coordinates defined as
£ = 2(x-x1)/(x2-x1)-l (C.SOa)
77 = 2(y-y1)/(y2-y1)-l (C.SOb)
The general expressions for Nj, N2 and N3 are given by
= -^(l-tfO/2 (C.Sla)
= l-N!-N2 (C.Slb)
= l-^2 (C.Slc)
where ^ = | or 77.
C.4 Code Structure and Dimension Limits of GPTRAC Module
Table C. 1 is a list of program segments (subroutines) that compose the GPTRAC
module. A simplified flow chart for GPTRAC is presented in Figure C.9. The dimension
limits of the GPTRAC arrays are listed in Table C.2.
-------
GPTRAC Formulation
Table C.I
Description Of Key Subroutines In GPTRAC Module
(Subroutines Listed In Alphabetical Order)
Subroutine Name
Description
CPZONE Determines time-related capture zone boundary for pumping wells.
DTCOMP Computes the required time step for a particle to travel from its current
location to an exit point at an edge of the rectangular grid block.
ECHO1 Echoes input parameters to default output file GPTRAC.OUT.
ELPROP Identifies the material zone associated with each finite element or finite
difference grid block.
FRPATH Sets up the loop for forward or reverse pathline tracking.
GPTRAC Execution supervisor for the grid generation, capture zone delineation,
and forward and reverse particle tracking modules in GPTRAC.
GRIDGN Generates a rectangular grid for pathline computation.
HEDCTR Computes the nodal values of hydraulic head at the centroids of
rectangular elements.
HLTCHK Determines if pathline termination should occur due to presence of flow
domain boundary or exceedence of simulation time.
PCHECK Performs pathline computation using the Euler iterative predictor-
corrector algorithm for numerical integration.
PLIDEN Identifies the rectangular element in which the particle currently resides.
PTRACK Performs forward or reverse pathline tracking, using analytical or
numerical integration.
-------
C-32 flPTR AC Formulation
Table C.I (continued)
Description Of Key Subroutines In GPTRAC Module
(Subroutines Listed In Alphabetical Order)
Subroutine Name
Description
PTRAM Performs pathline computation using the analytical integration algorithm.
SWITCH Switches the node numbering of input head values from row-wise to
column-wise.
VCALBL Computes the principal velocity components of each rectangular grid
block or element, and identifies the velocity pattern.
VINTQ Computes values of velocity components at a point inside a rectangular
grid block or element using quadratic interpolation.
VREVRS Reverses the signs of velocity components for reverse pathline tracking.
WFANS Computes the analytical steady-state flow solution in a well field subject
to ambient ground-water flow.
WNIDEN Identifies the sequential node numbers corresponding to the specified
well locations.
-------
GPTRAC Formulation C-33
Figure C.9
Simplified Flow Chart For GPTRAC Module
START }
"^
Input problem variable*
using processor and
ASOI file
Read bead file.
generate grid and
compute velocities
from head data
I
1
dividu
Pathllne
Tracldn
Call fBPATH to preform
pathline computation
saisg nnmwioel or
integration
C«ll CPZOMX to perform
captor* lone delineation
mlng rerene traoktnf
Plot reenlU
STOP
-------
C-34 GPTRAC Formulation
Table C.2
Dimension Limits Of Key Array Variables Used In GPTRAC
Array Name
Dimension
Description
AXKEEP
AYKEEP
BZO
CONXZO
CONYZO
CORD
ETASTO
FSTART
HDCTR
HVAL
IAXCOD
IAYCOD
IPATX
IPATY
2916
2916
200
200
200
(3025,2)
4
(40,2)
2916
3025
2916
2916
2916
2916
Linear velocity interpolation coefficient for x-
direction for each grid block
Linear velocity interpolation coefficient for y-
direction for each grid block
Aquifer thickness for each zone
Hydraulic conductivity in the x-direction
for each zone
Hydraulic conductivity in the y-direction
for each zone
x and y nodal point coordinates
Gauss-Legendre quadrature points
x and y coordinates for starting location of each
forward tracked pathline (maximum = 40)
Hydraulic head values for the center of each
grid block or finite element
Nodal hydraulic head values
Flag for each grid block indicating if flow
across the block is uniform in the x-direction
Flag for each grid block indicating if flow
across the block is uniform in the y-direction
x-direction velocity pattern for each grid block
y-direction velocity pattern for each grid block
-------
GPTRAC Formulation C-35
Table C.2 (continued)
Dimension Limits Of Key Array Variables Used In GPTRAC
Array Name
IPWCZ
IRWCZ
IWCHK
IZOEL
NDWELL
NOP
NSTLIN
PLANG
PORZO
QPWELL
QRWELL
RSTART
TZISTO
VELX1
Dimension
50
20
2916
2916
50
(2916,4)
50
100
200
50
20
(40,2)
4
2916
Description
Indicates pumping wells that capture zones are
to be delineated for
Indicates recharge wells about which contami-
nant fronts are to be delineated (not used)
Flag array indicating if elements contain well
nodes
Material zone number associated with each
finite element or grid block
Contains node numbers that are pumping wells
Rectangular finite element connectivities
Number of pathlines to be traced for delinea-
tion of time-related capture zones
Angles at which pathlines emanate from the
well bore
Aquifer porosity for each heterogeneous zone
Pumping well discharge rate
Recharge well injection rate
x and y coordinates for starting location of each
reverse tracked pathline (maximum = 40)
Gauss-Legendre quadrature points
x- velocity at block edge Xj
-------
C-36 GPTRAf Formulation
Table (2.2 (continued)
Dimension Limits Of Key Array Variables Used In GPTRAC
Array Name
VELX2
VELY1
VELY2
VXEL
VYEL
XCTR
XGRIDL
XMAXZO
XMMZO
XPL
XPLFNT
XPLKP
XPWELL
XRWELL
YCTR
Dimension
2916
2916.
2916
2916
2916
2916
61
200
200
500
100
500
50
20
2916
Description
x-velocity at block edge x2
y-velocity at block edge yj
y-velocity at block edge y2
x-components of velocity at element or grid
block centers
y-components of velocity at element or grid
block centers
x-coordinates of finite element or grid
centroids
x-coordinates of gridline columns
block
Maximum x-coordinates of heterogeneous
aquifer zones
Minimum x-coordinates of heterogeneous
aquifer zones
x-coordinates of pathline
x-coordinates of pathline end points (not
x-coordinate of pumping wells
x-coordinates of recharge wells
y-coordinates of finite element or grid
centroids
used)
block
-------
GPTRAC Formulation C-37
Table C.2 (continued)
Dimension Limits Of Key Array Variables Used In GPTRAC
Array Name
YGRIDL
YMAXZO
YMMZO
YPL
YPLFNT
YPLKP
YPWELL
YRWELL
Dimension
61
200
200
500
100
500
50
20
Description
y-coordinates of gridline rows
Maximum y-coordinates for
aquifer zones
Minimum y-coordinates for
aquifer zones
y-coordinates of pathline '
heterogeneous
heterogeneous
y-coordinates of pathline end points (not used)
y-coordinate of pumping wells
y-coordinate of recharge wells
-------
Appendix D
Saving to WordPerfect
D.I Retrieving HPGL Plotter Files Into WordPerfect 5.0
Before trying to use this feature, make sure the version of WordPerfect 5.0 you are
running is dated no earlier than April 29, 1989 (your help key (F3) will 'indicate the version
you are running in the upper right hand comer of your "Help" screen). If you are running
an earlier version, call WordPerfect and ask them for the update.
With a blank screen in WordPerfect 5.0:
1. ALT F9 (Graphics Menu Screen)
1 Figure
2 Create
3 Filename: ENTER CORRECT PATH AND FILE NAME AND HIT the
ENTER KEY
WordPerfect 5.0 will prompt you in the lower portion of the screen ~ "Please Wait -
loading plotter file" message will appear. After it has loaded the plotter file, it will insert
your file name in No. 1 Filename with a "(graphic)". If you get a message "Incorrect File
Format", check your version of WordPerfect 5.0
-------
£)-2 Saving to WordPerfect
You can change the figure definition by: (Graphics Menu Definition Screen):
3 Type (Paragraph. Page, or Character)
4 Vertical Postilion - offset from top of paragraph, page or character (the same
selection as No. 3)
5 Horizontal Position - 1) Left, 2) Right 3) Center, 4) Both
6 Size
1 Wrap T&Around Box? Y or N
To view and/or edit your plotter file (Graphics Menu Definition Screen):
8 Edit
Your file is then displayed in a graphics box with a graphics edit menu at the bottom
of your screen: >
1 Move Use arrow keys to move
2 Scale Use PgUp/On keys to scale
3 Rotate
4 Invert On
5 Black & White
When you are finished editing and/or viewing your graphic, press F7 to get out of the
Graphics Menu. Hit the ENTER key several times, and you will see a "Fig 1" box appear.
Keep hitting the enter key until you get to the end of your graphics box to see the box
appear on your screen.
To save your document in WordPerfect 5.0
F7 Will prompt "Save Document? (Y/N)?" in lower left-hand corner Of screen
Y Will prompt "Document to be Saved:" Enter name of your document, ENTER
Will prompt "Exit WP: (Y/N)?" Say no.
Your document is now saved in WordPerfect 5.0.
-------
Saving to WordPerfect 1~)_3
To print your document in WordPerfect 5.0:
Shift F7 Will then give you the print menu screen
3 "Document on Disk" menu selection
Enter name of document to be printed>
You will then be prompted ("Page(s): (ALL)" hit the ENTER key.
Your document should start printing.
D.2 Retrieving HPGL Plotter Files Into WordPerfect 5.1
With a blank screen in WordPerfect 5.1:
1. ALT F9 (Graphics Menu Screen)
1 Figure
2 Create
3 Filename: ENTER CORRECT PATH AND FILE NAME AND HIT the
ENTER KEY
WordPerfect 5.1 will prompt you in the lower portion of the screen ~ "Please Wait -
loading plotter file" message will appear. After it has loaded the plotter file, it will insert
your file name in No. 1 Filename with a "(graphic)".
You can change the figure definition by: (Graphics Menu Definition Screen):
2 Contents (Displays "Empty" until you have loaded a graphic file - will then
say "Graphic")
3 Caption You can give your plotter file a caption name. Follow prompts on
screen to enable the caption mode.
4 Anchorr Type (Paragraph. Page, or Character)
5 Vertical Position - offset from top of paragraph, page or character (the same
selection as No. 4)
6 Horizontal Position - 1) Left, 2) Right. 3) Center, 4) Both
7 Size
-------
D-4 Saving to WordPerfect
8 Wrap Text Around Box? Y or N
To view and/or edit your plotter file (Graphics Menu Definition Screen):
9 Edit
Your file is then displayed in a graphics box with a graphics edit menu at the bottom
of your screen:
1 Move Use arrow keys to move
2 Scale Use PgUp/Dn keys to scale
3 Rotate
4 Invert On
5 Black & White
When you are finished editing and/or viewing your graphic, press F7 to get out of the
Graphics Menu. Hit the ENTER key several times, and you will see a "Fig 1" box appear.
Keep hitting the enter key until you get to the end of your graphics box to see the box
appear on your screen.
To save your document in WordPerfect 5.1:
F.7 Will prompt "Save Document? (Y/N)?" in lower left-hand corner of screen
Y Will prompt "Document to be Saved:" Enter name of your document, ENTER
Will prompt "Exit WP: (Y/N)?" Say no.
Your document is now saved in WordPerfect 5.1.
-------
Saving to WordPerfect £)-5
To print your document in WordPerfect 5.0
Shift F7 Will then give you the print menu screen
3 "Document on Disk" menu selection
Enter name of document to be printed>
You will then be prompted ("Page(s): (ALL)" hit the ENTER key.
Your document should start printing.
-------
Appendix E
HEDCON Utility
User's Guide
E.I Introduction and Features
The HEDCON program (CONvert HEaDs) supplied with the WHPA model is
designed to create hydraulic head input files in the format required by the GPTRAC
numerical option. The format required by GPTRAC is outlined in detail in Section 8.3.3.
The two main options available in HEDCON are:
1. Convert a POSTMOD output file in x,y,h format into a file suitable for input into
GPTRAC. This option should only be used if the USGS MODFLOW code
(McDonald and Harbaugh, 1988) was used to calculate a head field, and the
postprocessing code POSTMOD (available from the International Ground-Water
Modeling Center at Indianapolis, Indiana) was used to create a file in x,y,h
format with the origin of the coordinate system at the lower left-hand comer of
the grid.
2. Convert a randomly ordered file in x,y,h format into a file suitable for input into
GPTRAC.
The x,y,h format referred to above simply means that there is one x-coordinate, one
y-coordinate, and one hydraulic head value on each line of the file. This type of format is
commonly used in commercial software packages such as SURFER (Golden Software). The
x- and y-coordinates must correspond to the location of a grid node, and the head value (h)
must be the head value at that node. Either a mesh-centered or a node-centered grid may
be specified.
-------
E-2
HEDCON Utility User's Guide
The HEDCON program input requirements are outlined in Table E.I. The grid
parameters must be specified (either interactively or by reading a previously created file) so
that HEDCON can match the nodal coordinates of the grid with the coordinates read from
the input head file. HEDCON operates very similar to the WHPA model. The Escape key
may be used to display an options menu, and there is full Help support throughout the
program. To execute HEDCON, type HEDCON at the DOS prompt. HEDCON may not
be executed from within the WHPA code.
Table E.I
Input Requirements for HEDCON Program
Program Variable
Description
XMIN:
XMAX:
YMIrl:
YMAX:
Minimum x-coordinate of grid (ft or m)
Maximum x-coordinate of grid (ft or m)
Minimum y-eoordinate of grid (ft or m)
Maximum y-coordinate of grid (ft or m)
NROWS and NCOLS,
and XGRiDL(I),
YGRIDL(J) for
1=1, NCOLS and
J=l, NROWS
Number of grid-line rows and columns, and (if
non-uniform grid,) coordinates of each grid-
line row and column. Note that this data may
be read in using the grid file option (FILE2
below).
FILE1
FILE2
Input file name that contains head values for
each node of the finite element or finite differ-
ence grid. This file must be in x,y5h format.
Input file name that contains the number of
grid-line rows and columns and the x-coordi-
nate of each grid-line column and the y-coor-
dinate of each grid-line row (optional).
-------
HEDCON Utility User's Guide
E.2 Grid' File Option
HEDCON will automatically prompt the user for the number of grid-line columns and
rows, and the minimum and maximum x- and y-coordinates of the grid. If the grid is
uniform (the grid-line spacings in both the x and y directions is constant), no additional
information is required to generate the grid. If the grid is nonuniform, however, the x-
coordinates of the grid-line columns and the y-coordinates of the grid-line rows must be
supplied to HEDCON. These parameters may either be entered interactively, or they may
be read from a previously created grid file, the name of which may be entered into
HEDCON. If the grid file option is selected, the file must be constructed prior to executing
HEDCON, and it must conform to the following format:
Number of grid-line columns (NCOLS) on line 1
x-coordinate of grid-line column 1 on line 2
x-coordinate of grid-line column 2 on line 3
x-coordinate of grid-line column 3 on line 4
,,, and so on for each grid-line column
Number of grid-line rows (NROWS) on line NCOLS+ 1
y-coordinate of grid-line row 1 on line NCOLS+2
y-coordinate of grid-line row 2 on line NCOLS+3
y-coordinate of grid-line row 3 on line NCOLS+4
,,, and so on for each grid-line row
Refer to one of the ***.GRD files in the EXAMPL subdirectory on the WHPA model.
diskettes for an example of a grid imput file.
If you decide to enter the grid-line coordinates directly into HEDCON, the grid data
will automatically be written to the file GRID.SAV. This file may then be read by
GPTRAC, so the grid data does not have to be entered twice. Note that once a GPTRAC
input file has been saved, the necessary grid data is retained in the saved file and the
separate grid file is no longer required as input.
E-3
-------
Appendix F
Development History
ofWHPACode
WHPA Version 1.0 Released February, 1990
WHPA Version 2.0 Released March, 1991
Enhancements include:
New analytical solutions for unconfmed aquifers subject to area] recharge and
leaky-confined aquifers subject to vertical leakage were implemented into
GPTRAC.
New analytical solutions for strip aquifers with two parallel stream or barrier
boundaries or a stream and a barrier boundary were implemented into
GPTRAC.
The same leaky-confined aquifer analytical solution implemented into GPTRAC
was implemented into MONTEC.
The GRAF module (postprocessor) was substantially enhanced. The upgrades
included the ability to obtain hard copy output from a number of commonly used
devices other than the HP 7475A pen plotter; the ability to overlay up to fifteen
plot files at one time; and the ability to execute the GRAF module options (e.g.,
map scale and overlay) prior to viewing the current plot.
Numerous corrections, enhancements and upgrades were made to the prepro-
cessor, including the ability to exit to DOS without exiting the WHPA code.
-------
F-2 Development History of WHPA Code
The HEDCON utility program was revised to include menus, help screens and
error prompts. The ability to read in a file in x,y, head format in any order, and
to subsquently reorder the file in a format suitable for the GPTRAC numerical
option was also implemented.
The restriction of a math coprocessor was eliminated: WHPA 2.0 will run on
IBM compatible microcomputers with or without a math coprocessor.
-------
References
Golden Software, Inc. SURFER, Version 4 Reference Manual. Golden, Colorado 80402.
Javandel, I., C. Doughty, and C.F. Tsang, 1984. Groundwater Transport: Handbook of
Mathematical Models. Water Resources Monograph 10, American Geophysical
Union, Washington, DC.
Lee, Kenrick H.L., 1986. Pollution Capture Zones for Pumping Wells in Aquifers with
Ambient Flow. Masters Thesis, Hydrology Program, New Mexico Institute of Mining
and Technology, Socorro, New Mexico.
Lee, K.H.L. and J.L. Wilson, 1986. Pollution Capture Zones for Pumping Wells in Aquifers
with Ambient Flow. EOS, Transactions of the American Geophysical Union, v. 67,
p. 966.
McDonald, M.G. and A.W. Harbaugh, 1988. A Modular Three-Dimensional Finite
Difference Ground-Water Flow Model. Techniques of Water-Resources Investigations
of the USGS, Book 6 Chapter Al.
Newsom, J. M., and J.L. Wilson, 1988. Flow of Ground Water to a Well Near a
Stream - Effect of Ambient Ground-Water Flow Direction. Ground Water, v. 26, no.
6, pp. 703-711.
Pollock, D.W., 1988. Semianalytical Computation of Path Lines for Finite-Difference
Models. Ground Water, v. 26, no. 6, pp. 743-750.
Shafer, J.M. 1987. GWPATH: Interactive Ground-Water Flow Path Analysis. Bulletin 69,
State of Illinois Department of Energy and Natural Resources.
------- |