U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2007–5259

Prepared in cooperation with

Barnstable County and the Cape Cod Commission

Use of Numerical Models to Simulate Transport of 

Sewage-Derived Nitrate in a Coastal Aquifer, Central and 

Western Cape Cod, Massachusetts
Cover. View to the east from Salt Pond Outlet over Salt Bay Pond, Eastham, Massachusetts. Photo by Andrew 

Massey, U.S. Geological Survey
Use of Numerical Models to Simulate 

Transport of Sewage-Derived Nitrate in 

a Coastal Aquifer, Central and Western 

Cape Cod, Massachusetts

By Donald A. Walter

Prepared in cooperation with Barnstable County and the Cape Cod Commission

Scientific Investigations Report 2007–5259

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

U.S. Geological Survey, Reston, Virginia: 2008

For product and ordering information:

World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

For more information on the USGS--the Federal source for science about the Earth, its natural and living resources, 

natural hazards, and the environment:

World Wide Web: http://www.usgs.gov

Telephone: 1-888-ASK-USGS

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Walter, D.A., 2008, Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer, 

Central and Western Cape Cod, Massachusetts: U.S. Geological Survey Scientific Investigations Report 2007–5259,

41 p.
iii

Contents

Abstract ...........................................................................................................................................................1

Introduction.....................................................................................................................................................2

Purpose and Scope ..............................................................................................................................2

Methods of Analysis .............................................................................................................................4

Hydrogeology.........................................................................................................................................4

Geologic Setting ...........................................................................................................................4

Hydrologic Setting .......................................................................................................................6

Nitrogen Transport in Ground Water ..........................................................................................................8

Numerical Models .........................................................................................................................................8

Regional Models .................................................................................................................................10

Regional Model Design .............................................................................................................10

Modifications to the Regional Models ...................................................................................10

Subregional Model .............................................................................................................................12

Model Design..............................................................................................................................12

Model Stresses ..........................................................................................................................12

Solute Transport .........................................................................................................................14

Simulation of Nitrate Transport .................................................................................................................14

Advective Transport and Discharge to Wells and Ecological Receptors .................................15

Discharge to Wells and Ecological Receptors .....................................................................15

Advective Transport in Different Hydrologic Settings .........................................................20

Solute-Transport Simulations ...........................................................................................................22

Comparison to Advective-Transport Simulations .................................................................31

Effects of Dispersion .................................................................................................................35

Limitations of Analysis ................................................................................................................................37

Summary........................................................................................................................................................37

References Cited..........................................................................................................................................40
iv

Figures

1. Map showing location of the Sagamore and Monomoy flow lenses and wastewaterdisposal facilities (WDF), Cape Cod, Massachusetts.............................................................3

2–4. Map and schematic section showing—

2. (A) The surficial geology of western and central Cape Cod, and (B) geologic 

section C-C’ through western Cape Cod, Massachusetts............................................5

3. (A) The regional water table and general flow paths for central and western 

Cape Cod and (B) generalized vertical ground-water flow, western Cape Cod, 

Massachusetts.....................................................................................................................7

4. (A) Finite-difference grids and simulated hydrologic boundaries, and (B) vertical 

discretization for regional models of Sagamore and Monomoy Flow lenses, 

Cape Cod, Massachusetts .................................................................................................9

5. Map showing current (2006) and potential sewered areas and areas of centralized 

wastewater disposal, central and western Cape Cod, Massachusetts............................11

6. Grid showing (A) hydrologic boundaries in the subregional model and simulated 

water-table contours from the regional and subregional models and (B) recharge 

zones in the subregional model, Barnstable and western Yarmouth, Cape Cod, 

Massachusetts............................................................................................................................13

7. Graph showing simulated hydrologic-budget components in relation to flow from 

the subregional model and the coincident parts of the regional model............................14

8. Grid and cross section showing (A) simulated water-table elevations, traveltimes, 

and general advective-flow paths, and (B) generalized vertical ground-water flow 

for unstressed conditions for Barnstable and western Yarmouth, Cape Cod, 

Massachusetts............................................................................................................................16

9. Grid showing simulated drawdowns and simulated particle tracks representing 

advective-flow paths downgradient from wastewater-disposal sites for a 

hypothetical scenario, Barnstable, Cape Cod, Massachusetts..........................................17

10. Graphs showing model-estimated wastewater discharge rates into (A) coastal water 

bodies, (B) ponds, streams, and (C) production wells for a hypothetical wastewaterdisposal scenario, Barnstable and Yarmouth, Cape Cod, Massachusetts .......................19

11. Boxplots showing model-estimated traveltimes to (A) coastal water bodies and 

(B) production wells, Barnstable and Yarmouth, Cape Cod, Massachusetts...................21

12. Boxplots showing (A) ground-water velocities and (B) traveltimes for advectivetransport paths from four wastewater-disposal sites for a hypothetical 

wastewater-disposal scenario, eastern Barnstable, Cape Cod, Massachusetts ...........23

13. Cross sections showing vertical particle tracks representing advective-transport 

paths under the condition of no pumping downgradient from (A) site C and (B) site B 

near the Barnstable Wastewater Disposal Facility, Cape Cod, Massachusetts .............24

14. Grids and cross sections showing (A-F) simulated nitrate concentrations 

downgradient from the current wastewater-disposal facility for six horizontal 

sections and (G-H) simulated nitrate concentrations along vertical column 147 

and row 304..................................................................................................................................26

15. Graph showing simulated time-varying nitrate concentrations at well 89, Fawcetts 

Pond, Stewarts Creek, and Nantucket Sound, Cape Cod, Massachusetts ......................28

16. Grid showing simulated nitrate concentrations downgradient from five wastewaterdisposal sites and simulated advective-flow paths in model layer 13 (50 to 60 ft 

below NGVD 29) near the Barnstable Wastewater-Disposal Facility, Cape Cod, 

Massachusetts............................................................................................................................29
v

Tables

 1. Wells, ecological receptors, and corresponding codes in modeled area of western 

and central Cape Cod, Massachusetts. ..................................................................................18

 17. Cross sections showing nitrate concentrations and particle tracks representing 
vi

Conversion Factors, Abbreviations, and Datum

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

square mile (mi2) 2.590 square kilometer (km2) 

Flow rate

foot per day (ft/d) 0.3048 meter per day (m/d)

million gallons per day (Mgal/d) 0.04381 cubic meter per second (m3/s)

inch per year (in/yr) 25.4 millimeter per year (mm/yr)

Hydraulic conductivity

foot per day (ft/d) 0.3048 meter per day (m/d)

Transmissivity*

foot squared per day (ft2/d) 0.09290 meter squared per day (m2/d) 

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as follows:

°C=(°F-32)/1.8

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29).

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83).

Elevation, as used in this report, refers to distance above the vertical datum.

*Transmissivity: The standard unit for transmissivity is cubic foot per day per square foot times 

foot of aquifer thickness [(ft3/d)/ft2]ft. In this report, the mathematically reduced form, foot 

squared per day (ft2/d), is used for convenience.

Concentrations of chemical constituents in water are given either in milligrams per liter (mg/L) 

or micrograms per liter (μg/L).
vii

Acronyms and Abbreviations

BWDF Barnstable Wastewater-Disposal Facility

FDM finite-difference method

MassDEP Massachusetts Department of Environmental Protection

MEP Massachusetts Estuaries Project

MMR Massachusetts Military Reservation

MOC method of characteristics

SMAST School of Marine Science and Technology

WDF wastewater-disposal facility

TMDL total maximum daily load

USEPA U.S. Environmental Protection Agency

USGS U.S. Geological Survey

mg/L milligrams per liter

kg/L kilograms per liter

kg/d kilograms per day
viii
Use of Numerical Models to Simulate Transport of 

Sewage-Derived Nitrate in a Coastal Aquifer, Central and 

Western Cape Cod, Massachusetts

By Donald A. Walter

Abstract

The unconsolidated glacial sediments underlying Cape 

Cod, Massachusetts compose a regional aquifer system that 

is used both as a source of drinking water and as a disposal 

site for wastewater; in addition, the discharge of clean ground 

water from the aquifer system is needed for the maintenance 

of freshwater and marine ecosystems throughout the region. 

Because these uses of the aquifer conflict with one another 

in many areas of the Cape, local and regional planners have 

begun to develop sustainable wastewater plans that will facilitate the disposal of wastewater while protecting water supplies 

and improving the health of aquatic ecosystems. To assist 

local and regional planners in these efforts, the U.S. Geological Survey conducted a 2-year investigation to (1) assist local 

and regional planners in the evaluation of potential wastewater 

scenarios, (2) use results and interpretation from these analyses to develop hydrologic concepts transferable throughout the 

region, and (3) establish and test methods that would be of use 

in future evaluations.

Wastewater-disposal scenarios need to be evaluated in 

the context of the regional ground-water-flow system. For 

a given rate of disposal, wastewater from sites at or near a 

regional ground-water divide is transported in a wider arc of 

flow directions, flows deeper in the system, and contaminates 

a larger part of the aquifer than does wastewater discharged 

from sites farther from the divide. Also, traveltimes of wastewater from sites near a ground-water divide to receptors are 

longer (as much as several hundred years) than traveltimes 

from sites farther from the divide. Thus, wastewater disposal 

at or near a divide will affect a larger part of the aquifer and 

likely contribute wastewater to more receptors than wastewater disposal farther from a divide; however, longer traveltimes 

could allow for more attenuation of wastewater-derived nitrate 

from those sites.

Ground-water-flow models and particle tracking can 

be used to identify advective-transport patterns downgradient from wastewater-disposal sites and estimate traveltimes; 

however, these tools cannot predict the distribution of mass or 

concentrations of wastewater constituents, such as nitrate, in 

the aquifer. Flow-based particle-tracking analyses can be used 

to estimate mass-loading rates and time-varying concentrations at wells and ecological receptors by the accounting of 

mass-weighted particles discharging into the receptor of interest. This method requires no additional development beyond 

the flow model; however, post-modeling analyses are required. 

In addition, the method is based on the assumption that no 

mass is lost during transport, an assumption that likely is not 

valid in many systems. Solute-transport models simulate the 

subsurface transport of nitrate through the aquifer and predict 

the distribution of the mass of a solute in the aquifer at different transport times. This method does require additional model 

development beyond the flow model, but can predict timevarying concentrations at receptors. Estimates of mass-loading 

rates require minimal post-modeling analyses.

Time-varying concentrations and mass-loading rates 

calculated for wells in eastern Barnstable by the two methods generally were in reasonable agreement. Inherent in the 

flow-based particle-tracking method is the assumption that 

mass is conserved along a given flow line and that there is no 

spreading of mass in the aquifer. Although the solute-transport 

models also incorporate a system-wide conservation of mass, 

these models allow for a spreading of mass in the aquifer, and 

mass is not conserved along a given flow line. As a result, 

estimates of concentrations and mass loading rates generally 

were higher in particle-tracking analyses than in solute-transport simulations. Results from the two types of simulations 

agreed best for wells that receive large amounts of wastewater 

with short traveltimes (less than 10 years) because insufficient 

transport time likely had elapsed to make the assumption of 

conservation of mass along flow lines problematic. The results 

agreed less closely for wells that receive wastewater with long 

traveltimes because dispersion may have attenuated simulated 

concentrations. In these cases, the conservation of mass was 

a more problematic assumption. Given particle densities that 

are practicable, particle-tracking simulations may not represent advective transport in divergent flow fields with enough 

resolution to adequately estimate time-varying concentrations 

and mass-loading rates in those areas. Also, particle tracking 

is more likely to overpredict mass-loading rates for wells in 

areas where dispersion may be important, such as near the 

edge of a wastewater plume, at large transport distances from 

a source, or downgradient from a pond.
2 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

Introduction

Cape Cod, Massachusetts, (fig. 1) is undergoing rapid 

urbanization and development as the region transitions from 

a semirural region strongly dependent on agriculture and 

tourism to a suburban region with a large year-round population and substantial industry. Communities on Cape Cod are 

among the fastest growing in Massachusetts. The year-round 

population of the Cape has more than doubled since 1970, 

and some parts of the region have had population increases of 

more than 70 percent since 1990 (Walter and Whealan, 2004). 

Although population has increased substantially in most areas, 

large additional increases are still possible because many 

towns are well below projected buildout development.

The unconsolidated glacial sediments underlying Cape 

Cod constitute an aquifer system that is the sole source of 

potable water to these communities. Most of the population 

of Cape Cod resides over the Sagamore and Monomoy flow 

lenses, the largest of the six lenses that compose the regional 

aquifer system (fig. 1) (LeBlanc and others, 1986). Pumping 

from these lenses totaled about 25 Mgal/d in 2004 and is 

projected to increase to more than 39 Mgal/d by the year 2020. 

On average, about 7 percent of total recharge to the aquifer is 

withdrawn for water supply (Walter and Whealan, 2004); most 

of this pumped water (about 85 percent) is returned to the 

aquifer as wastewater, resulting in a total consumptive loss of 

about 1 percent. Most naturally recharged ground water (about 

93 percent) discharges into ecological receptors, including 

ponds, streams, and coastal water bodies; ground-water discharge from the aquifer system supports important freshwater 

and marine ecosystems.

Most wastewater enters the aquifer as untreated septicsystem effluent from residential areas or as treated wastewater 

from five wastewater-disposal facilities (WDFs) (fig. 1). The 

discharge of wastewater into the aquifer can adversely affect 

water supplies. Nitrogen, in the form of nitrate, is a common 

constituent in wastewater and can adversely affect human 

health. The discharge of sewage-derived nitrate into production wells is of concern on Cape Cod, and water suppliers 

monitor nitrate concentrations in production wells. The U.S. 

Environmental Protection Agency (USEPA) has established a 

drinking-water standard of 14 mg/L.

In addition to adverse effects on water supplies, wastewater can adversely affect critical habitats and ecosystems in 

both freshwater and saltwater. Wastewater contains nutrients 

such as phosphorus and nitrogen that can cause eutrophication in some waters. Several coastal water bodies on Cape 

Cod are considered eutrophic or mesotrophic with the primary 

cause attributed to the subsurface discharge of sewage-derived 

nitrogen into the coastal waters. Eutrophication in most coastal 

waters is nitrogen-limited, and the discharge of excess nitrogen into these waters can lead to the growth of harmful algae, 

decreases in water clarity, depletion of oxygen, and the loss 

of indigenous flora and fauna (School of Marine Science and 

Technology, 2005).

The Massachusetts Department of Environmental Protection (MassDEP) has initiated the Massachusetts Estuaries Project (MEP) to estimate the total maximum daily load 

(TMDL) of nitrogen for estuaries on Cape Cod. (A TMDL 

represents the maximum load that an individual estuary can 

receive without adverse effects on water quality or estuarine 

ecosystems.) At present (2005), several communities on Cape 

Cod are developing wastewater-management plans to (1) limit 

the discharge of nitrogen into the watersheds of sensitive 

coastal waters in accordance with the TMDLs established by 

the MEP and (2) protect current and potential water supplies 

from the adverse effects of wastewater disposal. Wastewatermanagement plans include the sewering of residential areas 

and the disposal of treated sewage effluent at centralized 

disposal locations. The redistribution of wastewater disposal 

can change the amount and distribution of nitrogen discharge 

into wells and ecological receptors and affect the hydrologic 

system by changing water levels, hydraulic gradients, and 

streamflows. Owing to the fact that communities share a 

regional aquifer system, the hydrologic effects of wastewatermanagement scenarios on individual communities are better 

understood in the context of the regional ground-water flow 

system. A regional analysis, as opposed to a town-by-town 

analysis, allows hydrologic analyses of water supply and 

wastewater-management strategies for an individual town to 

be done in a context that takes into account similar actions by 

neighboring towns. The approach also allows a sustainable 

balance between competing uses of the aquifer—water supply 

and wastewater disposal—to be achieved on a regional basis.

In 2003, the U.S. Geological Survey (USGS), in cooperation with Barnstable County and the Cape Cod Commission, 

began an investigation into the potential hydraulic effects of 

wastewater disposal on the regional aquifer system (water 

levels, hydraulic gradients, and streamflows) as well as on the 

wells and ecological receptors within the Sagamore and Monomoy flow lenses. Specifically, the objectives of the investigation were to use numerical models to (1) evaluate regional 

environmental effects of potential wastewater-disposal scenarios, (2) assess advective transport of sewage-derived contaminants in different hydrologic settings within the regional flow 

system, and (3) evaluate the use of solute-transport modeling 

to simulate the transport of nitrate, the wastewater constituent 

of most concern, in the Cape Cod aquifer system.

Purpose and Scope

This report documents hydrologic analyses of wastewater 

management on central and western Cape Cod. Examples are 

presented that illustrate the use of ground-water-flow models 

in evaluating the effects of wastewater disposal on the hydrologic system and the advective transport of sewage-derived 

contaminants from disposal areas to wells and ecological 

receptors. Case studies that compare and contrast advective transport of sewage-derived contaminants in different 

hydrologic settings are presented. The report also documents 
Introduction 

3

Figure 1. Location of the Sagamore and Monomoy flow lenses and wastewater-disposal facilities (WDF), Cape Cod, Massachusetts.

Massachusetts

Military

Reservation

EXPLANATION

Wastewater-disposal facility WDF

Massachusetts

Military

Reservation

Phlashes Br.

Hawes Run

Mill Cr.

Santuit R.

Long Pond

Bourne R.

Herring R. Pocasset Brook

Marston Mills R.

Lower  Mill Pond

Little R. Mashpee R.

Backus R.

Squateague Brook

Upper  Mill Pond

WhitePond

Hinckley Pond

Seymour Pond

Hamblin Pond Wequaquet 

 Lake

Johns Pond

LongPond

AshumetPond

Mashpee Pond Mary Dunn Pond

Shawme Pond

Spectacle Pond Lawrence Pond

Triangle Pond

Coonamesset Pond

ElbowPond

Town Cove

River

Bass

Quash

n

e

t

R.Ch

i

lds

R.

Co

o

n

a

m

e

s

s

e

tt

R.

unn

ma ed br.

Monomoy Flow Lens

WDF

WDF

WDF

WDF

WDF Baker

Pond

SnakePond

Monomoy Flow Lens Sagamore 

Flow Lens

MASSACHUSETTS

CAPE COD

0

0 10 KILOMETERS

10 MILES 5

5

Orleans

Chatham

Harwich

Brewster

Dennis

Yarmouth

Mashpee

Barnstable

Sandwich

Bourne

Falmouth

 Hatchville precipitation  gage

 USGS  ToxicsResearch  site Nantucket Sound

Cape Cod Bay

Buzzards Bay

Vineyard Sound

Atlantic Ocean

Pleasant Bay

Cape

Canal Cod

Base from U.S. Geological Survey digital data,Universal Transverse Mercator projection,Zone 19, 1:24,000, 1991

41o48'

41o33'

70o30' 70o00'
4 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

the development of a subregional-flow and solute-transport 

model of an area of central Cape Cod and presents examples 

of the use of this tool in simulating the conservative transport 

of nitrate to wells and ecological receptors. Finally, the report 

also compares and contrasts the use of different modeling tools 

(flow-based particle tracking and solute-transport models) in 

simulating the conservative transport of sewage-derived nitrate 

in the Cape Cod aquifer system.

Methods of Analysis

 Numerical ground-water-flow models were used to 

simulate the hydrologic systems of the Sagamore and Monomoy flow lenses on central and western Cape Cod. These 

models were developed in 2004 by the USGS, in cooperation with MassDEP. The models use the USGS software 

program MODFLOW-2000 (Harbaugh and others, 2000) to 

simulate three-dimensional ground-water flow in the region. 

Detailed documentation of these models, including model 

design, hydrologic boundaries, hydrologic stresses, aquifer 

characteristics, and model calibration is presented in Walter 

and Whealan (2004). The advective transport of contaminants 

in the aquifer was simulated by the USGS particle-tracking 

software program MODPATH (Pollock, 1994). Graphic 

display of simulation results was done by using a version of 

the GIS software application MODTOOLS (Orzol, 1997) 

modified to work with MODFLOW-2000. A subregional 

model was developed to simulate the transport of nitrate in the 

aquifer by use of the USGS software program MODFLOW-96 

(Harbaugh, 1996) and the linked solute-transport code MT3D 

(Zheng, 1990). Lateral boundaries for the subregional model 

were derived from the regional model of the Sagamore flow 

lens (Walter and Whealan, 2004).

Hydrogeology

The glacial sediments beneath Cape Cod generally are 

sandy, highly permeable, and can exceed 500 ft in thickness in 

the central and western parts of the Cape. The region receives 

substantial rainfall, and the unconsolidated sediments are 

an important source of potable water for local communities. 

The geologic history and hydrology of Cape Cod has been 

documented in numerous publications, including LeBlanc and 

others (1986), Masterson and others (1997a), Oldale (1992), 

and Uchupi and others (1996).

Geologic Setting

The unconsolidated sediments underlying Cape Cod 

consist of gravel, sand, silt, and clay that was deposited during 

the Pleistocene epoch between 15,000 and 16,000 years ago 

(Oldale and Barlow, 1986). Sediments were deposited at or 

near the edges of a retreating continental ice sheet through 

direct deposition from the ice and proglacial deposition from 

meltwater. The surficial geology of western and central Cape 

Cod is characterized by broad, gently sloping glacial-outwash 

plains and hummocky terrain associated with glacial moraines 

(fig. 2A).

Three general types of glacial deposits occur on Cape 

Cod—moraine, ice-contact, and outwash. Moraines were 

deposited in low-energy environments at or near the edge 

of the ice sheets and are either ablation moraines that were 

deposited in place by melting ice, such as the Buzzards Bay 

moraine, or tectonic moraines, such as the Sandwich moraine, 

that consist of reworked outwash sediments pushed into place 

by local readvances of the ice sheets (Uchupi and others 1996; 

Oldale, 1992; Oldale and O’Hara, 1984). Kames are icecontact deposits that were deposited in high-energy meltwater 

environments within holes in the ice. Ice-contact deposits 

also include sediments that were deposited by meltwater in 

high-energy fluvial environments near the ice margin. Outwash sediments, which compose most of the unconsolidated 

sediments underlying Cape Cod, were deposited by meltwater 

streams in depositional environments associated with proglacial lake deltas; these depositional environments are analogous 

to those found in present-day fluvial deltas (Oldale, 1992; 

Uchupi and others, 1996). Deltaic outwash sediments can be 

divided into three general facies: topset, foreset, and bottomset 

deposits (fig. 2B) (Masterson and others, 1997a). These sediments generally become finer with depth; in general, coarsegrained fluvial sand and gravel (topset sediments) overlie 

fine to medium sand deposited in nearshore lacustrine environments (foreset beds), which in turn are underlain by silty 

sand and clay deposited in offshore lacustrine environments 

(bottomset beds). The glacial sediments are underlain by basal 

till in most places (fig. 2B); basal till consists of fine-grained 

material that was produced by mechanical erosion of bedrock 

by movement of the overlying ice sheet. The unconsolidated 

glacial sediments are underlain by relatively impermeable 

crystalline bedrock.

The elevation of the bedrock surface on Cape Cod ranges 

from about 50 ft below NGVD 29 near the Cape Cod Canal 

to more than 900 ft below NGVD 29 beneath the outer part of 

Cape Cod. The thickness of the glacial deposits on central and 

western Cape Cod ranges from 70 ft near the Cape Cod Canal 

to more than 500 ft along Nantucket Sound (fig. 2B).

The lithology of the glacial sediments, which include 

moraines, kames and other ice-contact deposits, and stratified 

outwash, differs according to the environment in which 

the sediments were deposited. Deltaic outwash sediments 

generally are characterized by two grain-size trends: (1) a 

fining-down trend in which grain sizes decrease, and silt and 

clay content increases, with depth; and (2) a trend in which 

sediments become finer-grained, and coarse-grained sand 

and gravel deposits become thinner, with increasing distance 

from sediment sources (fig. 2B). Sediment sources generally 

were close to the present-day shore of Cape Cod Bay near the 

apex of the Sandwich and Buzzards Bay moraines and near 

the Dennis and Nauset Heights ice-contact deposits (fig. 2A), 
Introduction 5

Figure 2. (A) The surficial geology of western and central Cape Cod, and (B) geologic section C-C’ through western Cape Cod, 

Massachusetts.

Nantucket Sound

0

0 10 KILOMETERS

5 10 MILES

5

A

Falmouth 

 deposits

 (Qfd)

Buzzards Bay moraine (Qbm)

Sandwich moraine (Qsm)

C'

C

 Mashpee

plain (Qmo)

 Harwich

 plain (Qhw)

Barnstable

 plain 

 (Qbn)

Buzzards Bay

plain (Qbo) Cape Cod Bay 

Lake deposits (Qld) 

Nantucket ice-contact

 deposits (Qnd)

 Dennis ice-contact

 deposits (Qdi)

Nauset Heights 

 ice-contact

deposits (Qnh) Qdi

EXPLANATION

Buzzards Bay lobe deposits Cape Cod Bay lobe deposits

Cape Cod Bay

Qbm

Qfd

Qbm Buzzards Bay Moraine

Qfd Falmouth ice-contact deposits

Qsm

Qnd

Qmo

Qbn

Qbo

Qdi

Qhw

Qnh

Qld Qld Cape Cod Bay Lake deposits

Qnh Nauset Heights ice-contact deposits

Qhw Harwich Plain deposits

Qdi Dennis ice-contact deposits

Qbo Buzzards Bay Plain deposits

 Qbn Barnstable Plain deposits

Qmo Mashpee Outwash Plain deposits

Qnd Nantucket Sound ice-contact deposits 

Qsm Sandwich Moraine

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator projection,

Zone 19, 1:24,000, 1991

Qdi

41o48'

41o

33'

70o30' 70o

00'

Wakeby Pond

Mashpee Pond

Peters Pond

CAPE COD BAY

NANTUCKET SOUND

400

0

ALTITUDE, IN FEET RELATIVE

-500 TO NGVD 29

C C'

B

Bend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in sectionBend in section

Topset deposit Moraine

Bottomset deposit

Foreset deposit

Basal till

EXPLANATION

Depositional units

Kame

Lacustrine deposit
6 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

and sediments generally become finer to the south. Moraine 

sediments, which consist of gravel, sand, silt, and clay, are 

generally poorly sorted with a highly diverse lithology. Icecontact deposits and kame deposits are well-sorted coarsegrained sediments, and basal till is fine-grained, generally 

homogenous, and consists primarily of clay.

Grain-size and degree of sorting determine the watertransmitting properties of aquifer sediments. The trends in 

hydraulic conductivity of outwash sediments are parallel to the 

trends in grain size; the hydraulic conductivity of sediments 

generally decreases with increasing depth and with increasing distance from sediment sources, or generally southward 

(Masterson and others, 1997a). Previous investigations have 

identified general relations between sediment grain size and 

hydraulic conductivity, as determined from aquifer tests 

(Masterson and others, 1997a; Masterson and Barlow, 1994). 

Medium to coarse sand and gravel deposits have hydraulic 

conductivities that typically range from 200 to 350 ft/d. Fine 

to medium sands have hydraulic conductivities typically ranging from 70 to 200 ft/d. The hydraulic conductivities of very 

fine sand and silt typically range from 30 to 70 ft/d, and silt 

and clay deposits have hydraulic conductivities of between 

10 and 30 ft/d. Ice-contact and kame deposits consist generally of medium to coarse sand and gravel and have hydraulic 

conductivities similar to those of coarse-grained outwash 

deposits. The lithology of moraine deposits ranges from gravel 

and sand to silt and clay, and moraines generally have lower 

average hydraulic conductivities than outwash deposits. Most 

areas, including moraines, have trends of decreasing hydraulic 

conductivity with increasing depth.

Hydrologic Setting

The unconsolidated glacial sediments underlying Cape 

Cod compose an unconfined aquifer system that is surrounded 

by saltwater: Cape Cod Bay to the northeast, Cape Cod Canal 

to the northwest, Buzzards Bay to the west, and Vineyard and 

Nantucket Sounds to the south. The Sagamore and Monomoy 

flow lenses on central and western Cape Cod are the largest 

and southernmost of six separate ground-water-flow lenses 

that underlie Cape Cod (LeBlanc and others, 1986); the two 

flow lenses are hydraulically separated by the Bass River 

(fig. 3A). The Sagamore flow lens is hydraulically separated 

at its northwestern extent from mainland Massachusetts by 

the Cape Cod Canal. The Monomoy flow lens is hydraulically 

separated from an adjacent flow lens by Town Cove at its 

northeastern extent. Each flow lens represents a distinct 

aquifer system that is hydraulically separate from adjacent 

flow lenses. The aquifer systems are bounded below by 

impermeable bedrock and at the top by the water table across 

which recharge enters (fig. 3B). Recharge from precipitation 

is the sole source of water to the aquifer system. About 45 

in. of precipitation falls annually on Cape Cod; slightly more 

than half of the precipitation recharges the aquifer (LeBlanc 

and others, 1986). The remainder is lost to evapotranspiration; 

surface runoff generally is negligible owing to the sandy soils 

and low topographic relief of the area.

Ground water flows outward from regional ground-water 

divides towards natural discharge locations at streams, estuaries, and the ocean (fig. 3A). Most recharge flows through 

shallow sediments and discharges to streams and estuaries; 

ground-water recharging the aquifer near the ground-water 

divides flows deeper in the aquifer and discharges to the ocean 

(fig. 3B). Most ground water (about two-thirds) discharges 

into saltwater bodies. About 25 percent of ground water is 

discharged into freshwater streams and wetlands, and a small 

amount (less than 10 percent) is removed from the system for 

water supply (Walter and Whealan, 2004). Maximum watertable elevations of the Sagamore and Monomoy flow lenses 

are more than 65 and 30 ft above NGVD 29, respectively 

(fig. 3A). Water-table contours and ground-water-flow patterns 

are strongly affected locally by kettle-hole ponds because the 

ponds are areas in the aquifer with no effective resistance to 

flow. Kettle-hole ponds are flow-through ponds characterized 

by ground-water-flow paths converging in areas upgradient from the ponds where ground water discharges into the 

ponds, and diverging in downgradient areas where pond water 

recharges the aquifer. Some ponds have surface-water outlets 

where ponds drain into freshwater streams. Streams generally 

are areas of ground-water discharge (gaining streams) and 

receive water from the aquifer over most of their length. Some 

stream reaches may lose water to the aquifer (losing streams), 

particularly in areas downgradient from pond outflows.

Owing to high recharge rates (over 25 in/yr) and 

the generally high permeability of the aquifer sediments, 

advective transport likely is the dominant component of 

contaminant transport in the aquifer (LeBlanc, 1984). Some 

contaminant plumes around the Massachusetts Military 

Reservation (MMR) on western Cape Cod (fig. 1) have 

migrated more than 4 mi downgradient since the mid-1950s. 

Ground-water velocities of more than 1.5 ft/d have been 

observed at the USGS Toxics Substances Research Site at 

the former sewage-treatment facility at the MMR (fig. 1) 

(LeBlanc and others, 1991); ground-water velocity is a function of aquifer porosity, which is about 0.35 near the MMR 

(Garabedian and others, 1985).

The rate of advective transport in the aquifer is a function 

of the locations of source areas relative to regional groundwater divides; water recharged near divides, where groundwater flow is more vertical, moves more slowly than water 

recharged farther from divides, where horizontal flow dominates (fig. 3B) (Walter and Masterson, 2003; Walter and others, 2004). Traveltime, defined as the total time it takes water 

to move from recharge locations at the water table to natural 

discharge locations, is longest near regional ground-water 

divides and ranges from essentially zero adjacent to discharge 

boundaries to hundreds of years near ground-water divides 

(Walter and others, 2004).

About 7 percent of the water recharging the Sagamore 

and Monomoy aquifer system is removed for water supply 

(Walter and Whealan, 2004). Most of this water is returned to 
Introduction 7

Figure 3. (A) The regional water table and general flow paths for central and western Cape Cod and (B) generalized vertical groundwater flow, western Cape Cod, Massachusetts.

Town Cove

Massachusetts

Military

Reservation 20

30

10

30

10

20

10

20

30

40

50

60

20

30

10

30

40

50

60

10

20

Nantucket Sound

Cape Cod Bay

Buzzards Bay

Vineyard Sound

Pleasant

Atlantic Ocean

 Bay

Town Cove

Cape

Canal

Cod

o

A

General direction of ground-water flow

Line of equal water-table elevation—In feet

above NGVD 29. Contour interval is 5 feet.

EXPLANATION

Saline 

ground 

water

Bedrock

Freshwater-saltwater 

interface

0

-100

-200

-300

 200

 100

North South

Vineyard

 Sound

-400

-500

-600

ELEVATION, IN FEET RELATIVE TO

NGVD 29

Fresh ground water

Cape Cod

 Bay

Areal Recharge

B

C C'

C ́

C

0

0 10 KILOMETERS

5 10 MILES

5

0

0 10 KILOMETERS

5 10 MILES

5

Bass River

Sagamore Flow Cell

Monomoy Flow Cell

Massachusetts

Military

Reservation

30

30

Base from U.S. Geological Survey

digital data,Universal Transverse

Mercator projection,

Zone 19, 1:24,000, 1991

41o48'

41o

33'

70o30' 70o00'
8 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

the system as disposed waste-water, either as dispersed septicsystem discharge or as point discharges to the aquifer at sewage-treatment plants. Although most of the pumped water is 

returned to the aquifer—a process that creates a regional mass 

balance—the water is usually disposed of at some distance 

from where it was withdrawn, particularly in areas served by 

public water supply. Large-capacity production wells decrease 

ground-water levels and can affect natural resources by drying 

vernal pools, drawing down ponds, and decreasing streamflows by changing hydraulic gradients and either intercepting 

water that would have discharged to a surface-water body 

or by inducing infiltration from the stream (Masterson and 

others, 2002). In the vicinity of large wastewater-disposal 

facilities, ground water can form a mound, and wastewater 

disposal can increase streamflows.

Nitrogen Transport in Ground Water

Nitrogen occurs in natural waters primarily in the form 

of nitrate (NO3

-

), ammonia (NH4

+

), dissolved nitrogen gas 

(N2(g)), or organic nitrogen; the predominant species of 

inorganic nitrogen in ground water is determined by the pH 

and redox conditions of the water (Hem, 1981). Ammonia 

is a reduced form of nitrogen that strongly sorbs to aquifer 

sediments, whereas nitrate, an oxidized form of nitrogen, is an 

anion that does not sorb to aquifer sediments and is transported conservatively in ground water (Hem, 1981). LeBlanc 

(1984) reported that nitrate was the predominant form of 

nitrogen in treated sewage effluent at the MMR wastewatertreatment facility; Weiskel and Howes (1992) found that most 

nitrogen emanating from septic systems in a sandy coastal 

aquifer in mainland southeastern Massachusetts also was in 

the form of nitrate.

Nitrate can be lost or attenuated in an aquifer through 

the process of denitrification in which nitrate is reduced to 

nitrogen gas (Hem, 1981). In most aquifers the denitrification 

process is coupled with the oxidation of organic carbon present in aquifer sediments. The sandy aquifer underlying Cape 

Cod generally is oxic and poor in organic carbon (LeBlanc, 

1984), and therefore, nitrate attenuation through denitrification may be limited by available carbon. Weiskel and Howes 

(1992) found that a total of about 26 ±10 percent of nitrogen 

(primarily as nitrate) in septic-system effluent was lost during 

transport over small transport distances (less than about 600 

m). The mechanism for the loss, however, is unknown; much 

of the nitrate may have been lost within the septic systems or 

in the underlying unsaturated zone and not in the aquifer. Most 

data suggest that the oxic, carbon-poor conditions in the Cape 

Cod aquifer generally would be favorable for the persistence 

of nitrate in the subsurface.

The one-dimensional subsurface transport of a solute can 

be described as (modified from Fetter, 1988):

∂

∂ = ∂∂ + ∂∂ − ∂∂

C

t

D C

x D Cx v Cx L d x

2

2

2

2 − ∂

∂ + ∂∂

B

n

C

t

C

t

d s

rxn

⎧

⎪

⎩

⎫

⎭

⎪ ⎧⎪

⎩

⎫

⎭

⎪ ⎧⎪

⎩

⎫

⎭

⎪ ⎧⎪

⎩

⎫

⎭

⎪ ⎧⎪

⎩

⎫

⎭

⎪ ⎧⎪

⎩

⎫

⎭

⎪

(Dispersion) (Diffusion) (Advection) (Sorption) (Reactions)

(1)

where

C = concentration of dissolved solute,

t = time,

DL = longitudinal dispersivity,

Dd = diffusivity,

vx = linear ground-water velocity,

Bd = bulk density of aquifer sediments,

n = porosity,

Cs = concentration of sorbed solute, and

rxn = subscript indicating chemical or biological 

reactions.

Advection refers to transport of a solute with average groundwater flow; dispersion, a process that occurs in all natural 

aquifers, refers to a spreading of solute mass with respect 

to the path of average ground-water flow owing to aquifer 

heterogeneity and preferential flow. Ground-water-flow 

models are mathematical representations of the aquifer that 

account for advective transport only. In highly permeable, 

well sorted aquifers, such as that which underlies Cape Cod, 

advection generally is the dominant transport mechanism, 

and ground-water models can be used to simulate this part of 

the transport process. Solute-transport models, such as those 

applied in this investigation, supply the additional terms for 

sorption, diffusion, and simple one-component reactions 

like decay. In the case of nitrate, which is conservative and 

nonreactive, the last two terms of the equation (sorption and 

reactions) likely are negligible, and diffusion generally is 

negligible in advection-dominated systems. The result is a 

transport equation in which the change in concentration with 

time is a function of advection and dispersion. It should be 

noted that the assumption of conservative transport not only is 

generally supported by field data, but also results in worst-case 

estimates of mass-loading rates; therefore, this assumption 

is favorable in formulating strategies for adequate resource 

protection.

Numerical Models

The set of modeling tools used in the analyses documented in this report include recently (2004) developed 

regional models of the Sagamore and Monomoy flow lenses 

as well as a subregional model of parts of Barnstable and 

Yarmouth on central Cape Cod (fig. 4A). The numerical 

models are discretized mathematical representations of the 

aquifer system that can simulate hydrologic conditions (heads 

and flows) in the aquifer and the movement of water through 

the aquifer.
Numerical Models 9

 Sagamore model row 119

60

0

-100 Sagamore active model area

Inactive model area

Sagamore model Monomoy model

Monomoy active model area

-240

-140

C C'

Additional layers

-503

ALTITUDE, IN FEET

RELATIVE TO NGVD 29

 Monomoy model row 79

Orleans

Chatham

Brewster

Harwich

Dennis

Yarmouth

Mashpee

Barnstable

Sandwich

Bourne

Falmouth

0

0 10 KILOMETERS

10 MILES 5

5

AB

C C'

Buzzards Bay

Cape Cod Bay

Nantucket Sound

Vineyard Sound

Finite-difference grids

Active model boundary

Active model boundary

Coastal boundary (GHB)

Estuary boundary (DRN)

Saltwater boundary

Freshwater boundary

PondStream (STR)

EXPLANATION

Well

Subregional model

grid

Subregional model grid

Bass River

Town 

Cove

Cape Cod 

 Canal

Base from U.S. Geological Survey digital data,Universal Transverse Mercator projection, Zone 19, 1:24,000, 1991

41o48'

41o33'

70o30' 70o00'

Figure 4. (A) Finite-difference grids and simulated hydrologic boundaries, and (

B) vertical discretization for regional models of Sagamore and Monomoy 

Flow lenses, Cape Cod, Massachusetts.
10 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

Regional Models

The regional models of the Sagamore and Monomoy 

lenses are documented in detail in Walter and Whealan (2004), 

Appendix 1. These models were previously used to delineate 

recharge areas to wells and ecological receptors on Cape Cod 

(Walter and others, 2004) as well as to evaluate the effects of 

time-varying pumping and recharge on the regions’ hydrologic 

system.

Regional Model Design

The models extend from the Cape Cod Canal in Bourne 

to Town Cove in Orleans; the boundary between the active 

parts of the two model domains is the Bass River in Yarmouth 

(fig. 4A). The finite-difference ground-water modeling 

program MODFLOW-2000 (Harbaugh and others, 2000; 

McDonald and Harbaugh, 1988) was used to simulate the 

ground-water systems of the two flow lenses. The Sagamore 

model has 246 rows and 365 columns with a total active 

modeled area of 246 mi2, and the Monomoy flow model has 

164 rows and 220 columns with a total active modeled area of 

106 mi2. Both models have a uniform horizontal discretization 

of 400 by 400 ft; the grids are coincident where the models 

overlap (fig. 4A). Both models have 20 layers with thicknesses 

ranging from a maximum of 10 ft in the top 17 model layers to 

more than 200 ft in the deepest layer (fig. 4B).

Estuaries, open coastal waters, and streams are 

represented as head-dependent flux boundaries in the models 

(fig. 4A). Open coastal waters, saltwater estuaries, and some 

freshwater wetlands were simulated by using the General 

Head Boundary (GHB) and Drain (DRN) packages in 

MODFLOW-2000 (Harbaugh and others, 2000; McDonald 

and Harbaugh, 1988). Freshwater streams and some freshwater wetlands were simulated by using the Streamflow Routing 

Package (STR) (Prudic, 1989) which allows for an accounting 

of streamflow and pumping-induced streamflow depletion. 

Freshwater ponds, which were represented as regions of high 

hydraulic conductivity (100,000 ft/d), were active parts of the 

aquifer, and, therefore, pond levels fluctuate in response to 

changing hydraulic stresses.

The only input of water into the model is recharge from 

precipitation. A natural recharge rate to aquifer sediments 

of 27.1 in/yr was specified; this value was obtained from 

long-term records from the precipitation gage at Hatchville, 

Mass. (fig. 1), was adjusted during model calibration, and is 

consistent with recharge values used in previous investigations on Cape Cod and southeastern New England (Barlow 

and Dickerman, 2000; Desimone and others, 2001; Masterson and others, 1997b; Walter and Masterson, 2003). Natural 

recharge to surface-water bodies was adjusted according 

to the receptor: 16 in/yr to ponds representing net recharge 

after pan evaporation (Farnsworth and others, 1982) and no 

recharge to wetlands, which likely are areas of net groundwater discharge. Recharge was further adjusted to account 

for septic-system return flow in residential areas. Because 

parcel-scale water-use data were not available, the volume of 

generated wastewater was determined from the total volume 

of ground water withdrawn by each town; a consumptive loss 

of 15 percent was assumed, and the remaining volume was 

evenly distributed across model cells representing nonsewered 

residential areas. Recharge also was adjusted in areas of 

current wastewater-disposal facilities (WDFs) in Barnstable, 

Chatham, Falmouth, Orleans, and at the MMR (fig. 1). In 

2004, about 156 production wells operated on central and 

western Cape Cod (Walter and Whealan, 2004). Withdrawal of 

water at production wells was simulated using the Well (WEL) 

Package (McDonald and Harbaugh, 1988). The regional models of the Sagamore and Monomoy flow lenses simulate current ground-water withdrawals of about 17.3 Mgal/d and 7.6 

Mgal/d, respectively. The models also simulate future (2020) 

ground-water withdrawals of about 27.4 and 11.6 Mgal/d, 

respectively, for the two flow lenses.

Aquifer properties were estimated from lithologic logs 

and from previously developed depositional models of central 

and western Cape Cod (Masterson and others, 1997b; Byron 

Stone, U.S. Geological Survey, written communication, 2002). 

Hydraulic conductivities varied spatially with depth and 

ranged from 350 ft/d for coarse sand and gravel to 10 ft/d for 

silt and clay (Walter and Whealan, 2004). The steady-state 

regional models were calibrated on the basis of measured 

long-term water levels and streamflows as well as delineated 

contaminant plumes as indicators of advective flow paths. 

Initial estimates of recharge and intrinsic aquifer properties 

were adjusted during the calibration process to achieve an 

acceptable match between observed and simulated hydrologic 

conditions at the calibration points (Walter and Whealan, 

2004). Steady-state models were used in the analysis of advective transport. Although recharge and pumping vary through 

time and versions of the regional models that incorporate 

these time-varying stresses have been developed, steady-state 

models are adequate for the simulation of advective flow 

because the time scale of the variability of these stresses is 

small compared to the time scale of transport in the aquifer 

(Walter and Masterson, 2003).

Modifications to the Regional Models

The hydrologic boundaries, intrinsic aquifer properties, 

and natural recharge rates described in Walter and Whealan 

(2004) also were used in this investigation. Modifications 

were made to human-influenced recharge stresses simulated 

in the regional models. Changes in pumping rates were minor 

and were based on additional information made available from 

local water suppliers during the investigation.

The amount and spatial distribution of recharge representing septic-system return flow also was changed. The additional 

recharge representing septic-system return flow was removed 

from some areas to simulate the diversion of water from these 

areas to areas of centralized sewage disposal (fig. 5). Recharge 

was enhanced in areas of the model grid that represent possible locations of centralized wastewater disposal, simulat-
Numerical Models 11 Figure 5. Current and potential sewered areas and areas of centralized wastewater disposal, central and western Cape Cod, Massachusetts.

Buzzards Bay Nantucket Sound

Vineyard Sound

Site F Site C Site D

Site B

Barnstable

wastewater-disposal

Cape Cod Bay

facility

Cape

Canal Cod

0

0 10 KILOMETERS

10 MILES 5

5

Orleans

Chatham

Harwich

Brewster

Dennis

Yarmouth

Mashpee

Barnstable

Sandwich Bourne

Falmouth

Current (2006)Potential

Wastewater-disposal facility

Potential

EXPLANATION

Sewered areas

Towns not participating in study

Current (2006)

Base from U.S. Geological Survey digital data,Universal Transverse Mercator projection,Zone 19, 1:24,000, 1991

41o48'

41o33'

70o30' 70o00'
12 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

ing the disposal of wastewater onto sand-infiltration beds at 

these locations (fig. 5). Each modeled scenario represented 

an adjustment of model recharge according to this method. 

Owing to the fact that septic-system return flow was estimated 

from general residential land use as opposed to actual parcelscale water use, the volume of water made available through 

the conversion of on-site septic disposal to sewering often did 

not equal the volume of water specified for disposal at centralized locations. In these cases, return flow in nonsewered areas 

was adjusted to maintain a mass balance of water (adjusted for 

consumptive loss) within each town. In all cases, the discrepancy was less than 1 percent of total recharge in the model, 

and this simplification did not affect simulated advective flow 

paths in the aquifer.

Subregional Model

A subregional model was developed as part of this 

investigation to simulate ground-water flow and transport in 

the eastern part of Barnstable (fig. 4A). This area was chosen 

because the area is the most densely populated part of Cape 

Cod, and wastewater disposal is an issue of great concern 

in Barnstable. The town is in the process of developing a 

wastewater-management plan; at the time of this investigation 

(2006), the community had already identified four potential 

wastewater-disposal sites as wells as additional areas to be 

sewered.

Model Design

The model grid, which includes the eastern part of 

Barnstable and the western part of Yarmouth, encompasses 

an area of about 43.2 mi2 on central Cape Cod (fig. 4A). The 

model has 376 rows, 320 columns, and 26 layers; horizontal 

discretization is a uniform 100 ft on a side. Layer thicknesses 

from the water table to 100 ft below NGVD 29 do not exceed 

10 ft; this area of the model is of most interest because wells 

and hydrologic boundaries as well as most ground-water 

flow are within the upper part of the aquifer. Six layers were 

added to the subregional model below the elevation of 100 ft 

below NGVD 29: layer 18 in the regional model (100 to 140 

below NGVD 29) was equally subdivided into two layers, 

layer 19 (140 to 240 ft below NGVD 29) was subdivided into 

three layers, and layer 20 (240 to 516 ft below NGVD 29) in 

the regional model was subdivided into two layers. The finer 

vertical discretization was designed to limit vertical numerical dispersion at depth in the aquifer to the extent practicable. 

Most wells and all hydrologic boundaries are within the top 17 

layers where vertical discretization is 10 ft or less.

The modeled area includes outwash and ice-contact 

deposits to the south and moraine and lacustrine deposits in 

the north, along Cape Cod Bay (fig. 2). Hydraulic conductivities of the outwash sediments range from 220 ft/d in the 

shallow part of the aquifer to 10 ft/d in the basal sediments, 

and hydraulic conductivities of ice-contact deposits decrease 

from 260 to 10 ft/d with depth; the finer-grained moraine and 

lacustrine deposits have hydraulic conductivities ranging from 

100 to 10 ft/d at depth. The hydrologic-boundary geometries 

used in the regional model were incorporated into the subregional model, so that the two models were consistent (fig. 6A). 

Boundary leakances also were the same: 0.2 ft/d for open 

coastal waters (GHB boundaries), 0.1 ft/d for estuaries and 

salt marshes (DRN boundaries), and 1.0 ft/d for streams (STR 

boundaries) (Walter and Whealan, 2004).

The subregional model was linked to the regional model 

of the Sagamore flow lens by constant-head boundaries along 

the eastern, western, and northern edges of the subregional 

model; the version of the regional model used to generate 

boundary conditions was modified to have the same 26 layers 

as in the subregional model. For each subregional model run, 

hydraulic stresses were modified in the regional model to 

match those to be simulated in the subregional model. The 

resulting regional head distribution was used as a specified 

head boundary in the subregional model. On the basis of 2004 

stresses, simulated heads in the subregional model closely 

agreed with heads simulated by using the regional model 

(fig. 6A); the two sets of heads were within 0.3 ft of one 

another (or less than 1 percent of the total hydraulic gradient) 

within the subregional-model domain. The hydrologic budget 

of the subregional model closely matched the hydrologic budget of the coincident part of the regional model (fig. 7). Most 

budget components—recharge to the aquifer and discharge 

into streams, estuaries, and open coastal waters—were within 

1 percent of their regional-model equivalents; flows across 

the constant-head boundaries differed by about 5 percent from 

flows across the coincident boundaries in the regional model.

Model Stresses

The distribution of recharge was the same in both the 

regional and subregional models; recharge zones included 

aquifer sediments, ponds, and wetlands (fig. 6B). Recharge 

values were set to match the values used in the regional 

model: 27.1 in/yr to aquifer sediments, 16 in/yr to surfacewater bodies, no recharge to wetlands, and additional recharge 

to aquifer sediments in residential areas representing septicsystem return-flow. Twenty-five production wells are within 

the subregional-model domain, 17 in Barnstable and 8 in 

Yarmouth (fig. 6A). In 2006, the wells withdrew a total of 

4.81 Mgal/d and are projected to withdraw about 6.59 Mgal/d 

in the future (2020). Simulated pumping rates in 2006 were 

5 percent higher than the 2004 pumping rates reported in 

Walter and Whealan (2004); the projected 2020 pumping rate 

is about 5 percent lower than the 2020 rates reported in Walter 

and Whealan (2004).

The recharge distribution includes areas of current and 

potential sewering and five locations of possible centralized 

wastewater disposal, including the Barnstable wastewater-disposal facility (BWDF) currently in operation and four candidate sites (sites B, C, D and F) for disposal as identified by 
Numerical Models 13

Figure 6. (A) Hydrologic boundaries in the subregional model and simulated water-table contours from the regional and subregional 

models and (B) recharge zones in the subregional model, Barnstable and western Yarmouth, Cape Cod, Massachusetts.

Yarmouth

Barnstable

Yarmouth

Barnstable

Site F

Site C

Site B

Site D

BWDF

Barnstable (2020 sewered)

Yarmouth (2006 sewered)

Pond

Freshwater wetland

Water body

EXPLANATION

Natural aquifer only

Septic-system return flow

Saltwater body

Barnstable

Barnstable (2006 sewered)

Potential

Current (2006)

Wastewater disposal

Coastal boundary (GHB)

Estuary boundary (DRN)

Saltwater boundary

Pond

Stream (STR)

Freshwater boundary

EXPLANATION

Well

Constant-head boundary

Simulated water-table contour—

 Number is feet above NGVD 29

Regional model

Subregional model

Potential wastewater-disposal site

Wastewater-disposal facility (BWDF)

70°20' 70°16'

41°

42'

41°

38'

BWDF

34

14

22

26

30

18

10

6

36

0

0 1

1/2 2 MILES

1/2 2 KILOMETERS

1

4

4

A B 70°20' 70°16'
14 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

Recharge

Estuaries

Coast

Streams

Constant head

Wells

0.0

1.0

2.0

3.0

4.0

5.0

6.0

FLOW, IN CUBIC FEET PER DAY (X 106

)

HYDROLOGIC-BUDGET COMPONENT

Regional model

Subregional model

Figure 7. Simulated hydrologic-budget components in relation 

to flow from the subregional model and the coincident parts of the 

regional model.

regional and local planners (Cape Cod Commission, written 

commun., 2004) (fig. 6B). Sewered areas and areas of centralized wastewater disposal are represented as different recharge 

zones that allow for the simulation of different combinations 

of sewering and disposal of different volumes of wastewater at 

the designated sites. The total volumes of wastewater disposal 

simulated in the different scenarios for this area ranged from 

1.5 to about 7.0 Mgal/d.

Solute Transport

The subsurface transport of nitrate was simulated by 

using the three-dimensional transport code MT3D (Zheng, 

1990). MT3D simulates several components of transport, 

including advection, mechanical dispersion, and diffusion, 

and can also simulate some single-species chemical reactions, 

including contaminant decay and various types of linear and 

nonlinear sorption. In the case of nitrate, which is conservative 

and nonreactive (Hem, 1981), the most important components 

of transport are advection and, to a lesser extent, mechanical 

dispersion. MT3D is linked to MODFLOW and uses steadystate cell-by-cell flows and velocities from the flow model 

to simulate the movement of nitrate through the aquifer. The 

simulation outputs are dissolved and sorbed concentrations of 

the solute of interest for cells within the active model domain 

at different transport time steps. This information can be used 

to develop maps of solute concentrations as well as concentration breakthrough curves.

In addition to the physical boundary conditions in the 

flow model, concentration boundaries were required for the 

solute-transport simulations. Areas of wastewater disposal 

(fig. 6B) were simulated as a source of nitrate that was continuous over the period of simulation. Treated wastewater was 

assumed to have a nitrate concentration of 10 mg/L (Thomas 

Cambareri, Cape Cod Commission, written commun., 2005). 

Areas of septic-system return flow also were represented as 

a continuous source of nitrate. Septic-system effluent was 

assumed to have a nitrate concentration of 35 mg/L (Weiskel 

and Howes, 1992) and was adjusted to account for dilution 

by natural recharge. Inflow concentrations in residential areas 

(fig. 6B) ranged from about 3.7 mg/L in Barnstable to about 

6.9 mg/L in Yarmouth.

The transport equation can be solved by using either the 

finite-difference method (FDM) or particle-based methods 

such as methods of characteristics (MOC). The FDM was used 

for these simulations because the method provides good conservation of mass and is computationally efficient. In advection-dominated systems, particle-based methods generally 

have less numerical dispersion; however, conservation of mass 

can be poor. The use of a mass-balance method, such as FDM, 

in advection-dominated systems can result in some numerical 

dispersion; however, the problem can be minimized by using 

numerical models with finer discretization. Numerical dispersion can be approximated as a function of model discretization, velocity, and time-step length with the equation

, 2 2

x v t

α

Δ Δ = + (2)

where

	 α = dispersivity,

x = cell length,

v = velocity, and

t = the time-step length (Fletcher, 1991).

With a cell length of 100 ft, a velocity of 1 ft/d, and a time 

step of 4 d, parameters that are consistent with the subregional 

model, the approximate numerical dispersion is 52 ft. The 

Peclet number, which is defined as the ratio of cell length to 

dispersivity, is about 2. For advection-dominated systems, a 

Peclet number of less than about 4 indicates that the FDM is 

appropriate for use in solving the solute-transport equation.

Simulation of Nitrate Transport

The regional models were used to address two sets of 

questions regarding the regional effects of wastewater disposal: (1) the hydrologic effects of sewering and centralized 

wastewater disposal, including changes in water levels in 

wells and ponds and in streamflows, and (2) the discharge of 

sewage-derived contaminants, particularly nitrate, into wells 

and ecological receptors. The regional effects of wastewaterdisposal scenarios in 8 of the 11 towns on the Sagamore and 
Simulation of Nitrate Transport 15

Monomoy flow lenses were analyzed as part of this investigation (fig. 5). Because each town had multiple disposal sites 

as well as a range of possible wastewater-disposal volumes, 

multiple scenarios were analyzed as part of the investigation. 

In this report, a modeling analysis of a single hypothetical 

scenario of nitrate transport is presented to feature hydrologic 

concepts transferable to other scenarios and to document 

methods used in the investigation.

Advective Transport and Discharge to Wells and 

Ecological Receptors

Particle tracking was used to simulate the transport of 

sewage-derived contaminants through the aquifer from source 

areas to discharge locations at wells and hydrologic boundaries. Particle-tracking algorithms use model-calculated cell-bycell flows to track the movement of particles of water through 

the subsurface. Advective transport, which is the predominant 

component of transport in the aquifer, is adequate for representing the movement of conservative, nonreactive solutes 

with ground-water flow if no loss of mass can be assumed. For 

this analysis, the contaminant of most concern in wastewater 

is nitrate, which generally is conservative and nonreactive 

in ground water. As a result, it can be assumed that nitrate 

is transported primarily by advection, and therefore, flowbased particle tracking can be used to simulate the transport 

of nitrate through the subsurface and to estimate the potential 

of effects of sewage-derived nitrate on wells and ecological 

receptors. Nitrate can be attenuated by denitrification in the 

presence of organic carbon (Colman and others, 2004), and 

some nitrate attenuation is assumed in estimating TMDLs for 

estuaries on Cape Cod (School of Marine Science and Technology, 2005). The aquifer sediments on Cape Cod, however, 

are generally carbon poor (LeBlanc, 1984), and it is assumed 

that nitrate is transported conservatively. If some nitrate is 

attenuated, the estimates of mass-loading rates to receptors in 

this analysis would represent worse-case conditions.

The movement of wastewater through the aquifer was 

quantified by assigning particles at wastewater-disposal sites 

and tracking the particles to their discharge locations. An 

accounting of particles was used to estimate the volumes of 

wastewater discharging into wells and receptors; the specified volume of wastewater and the distribution of assigned 

particles were uniform within disposal areas so that an equal 

volume of wastewater could be assigned to each particle and 

the volumes of wastewater discharging to receptors could be 

estimated from particle counts. These estimates were used, 

in turn, to estimate mass-loading rates of nitrate by assigning 

concentrations to the wastewater and assuming that no mass is 

lost during transport.

The hypothetical scenario features full sewering of 

residential areas within Barnstable and the disposal of 

5.36 Mgal/d of the generated wastewater at the BWDF and 

at sites B, C, and F (figs. 5 and 6B). Disposal volumes at 

the four sites are as follows: 1.70 Mgal/d at the BWDF, 

similar to current flows, and 1.22 Mgal/d each at potential 

sites B, C, and F. It should be noted that this scenario was 

and is not under consideration by local planners. However, 

this hypothetical scenario is similar in terms of locations 

and volumes of disposal to actual scenarios analyzed and 

illustrates the types of analyses done as part of this investigation as well as concepts generally pertinent to wastewater 

disposal in the region.

Discharge to Wells and Ecological Receptors

The water table in east Barnstable and west Yarmouth 

features an east-west trending ground-water divide; groundwater flow, and therefore advective transport, is away from 

the divide. Under 2004 stresses, which include pumping and 

a loading rate of 1.7 Mgal/d of wastewater at the BWDF, the 

divide trends east-west and is closer to coastal boundaries to 

the north than those to the south (fig. 8A). The divide deviates 

to the south near the BWDF owing to mounding associated 

with wastewater disposal at the site. The areas on the water 

table that contribute recharge to wells and natural receptors—

ponds, streams, and the coast—extend upgradient from the 

receptors in the direction of the regional ground-water divide 

(fig. 8A). These contributing areas can change in response to 

changing hydraulic stresses such as centralized wastewater 

disposal; any wastewater within the contributing area to a 

receptor will contribute water to that receptor.

For the hypothetical scenario, the simulated disposal 

of an additional 3.66 Mgal/d of wastewater (1.22 Mgal/d at 

sites B, C, and F) changed water levels and hydraulic gradient 

directions and magnitudes in the aquifer as well as the position 

of the regional ground-water divide (fig. 9). The divide shifted 

to the north in response to disposal at sites B and F (fig. 9) 

relative to its position (assuming 2006 hydraulic stresses) 

(fig 8A). Sites B and F were near but to the north of the divide 

when the sites were not in operation and were on top of the 

divide when the sites were in operation (fig. 9). The hydraulic 

effects of disposal generally are greater near divides because 

hydraulic gradients are less steep than in areas closer to discharge boundaries (figs. 6A and 8). As an example, mounding 

at the water table exceeded 6 ft at sites B and F—near the 

ground-water divide—whereas the same amount of wastewater disposal at site C, an area of steeper hydraulic gradients 

farther from the divide, resulted in mounding of about 2 ft 

(fig. 9).

Sewage-derived contaminants were transported away 

from the source areas and discharged into wells, ponds, 

streams, and coastal water bodies (fig. 9). The names of these 

receptors and their corresponding codes are presented in table 

1. Of the 5.36 Mgal/d of disposed wastewater, most (about 

70 percent) discharged directly into 10 different coastal-water 

bodies (fig. 10A). These water bodies (fig. 9) represent estuaries and open coastal-water bodies for which TMDLs will be 

evaluated as part of the Massachusetts Estuaries Project (Brian 

Howes, School of Marine Science and Technology (SMAST) 

and Cape Cod Commission, written commun., 2001). About 
16 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

Figure 8. (A) Simulated water-table elevations, traveltimes, and general advective-flow paths, and (B) generalized vertical 

ground-water flow for unstressed conditions for Barnstable and western Yarmouth, Cape Cod, Massachusetts.

10

34

18

41°42'

70°20' 70°16'

 

EXPLANATION

Modeled receptors

 Freshwater boundary

 Pond

 Stream (STR)

 Saltwater boundary

 Coastal boundary (GHB)

 Estuary boundary (DRN)

Contributing areas to wells

 and freshwater receptors

 Ponds

 Streams

 Wells

 Constant head

 Traveltime, in years

 0 to10

 greater than 10 to 20

 greater than 20 to 30

 greater than 30 to 50

 greater than 50 to 100

 greater than 100

Barnstable wastewater-

 disposal facility (BWDF)

General advective

 transport directions

 Water-table elevations—

 Number is feet above 

 NGVD 29

Well 

BWDF

34

14

26

30

18

10

6

6

36

22

B

A

C'

C

C C'

0

0 1

1/2 2 MILES

1/2 2 KILOMETERS

1

0

0

1/2 1 MILE

1/2 1 KILOMETER

Traveltime of vertical

ground-water flow, in years 

0 to 10

greater than 30 to 50

greater than 20 to 30

greater than 10 to 20

greater than 50 to 100

greater than 100

Ground-water

 divide

Ground-water divide

41°38'

18

18
Simulation of Nitrate Transport 17

Figure 9. Simulated drawdowns and simulated particle tracks representing advective-flow paths downgradient from wastewaterdisposal sites for a hypothetical wastewater-disposal scenario, Barnstable, Cape Cod, Massachusetts.

0

0 1 KILOMETER

1/2

1/2

1 MILE

Column 65 Column 233

 Current 

ground-water 

 divide

 Current 

ground-water 

 divide

 New 

ground-water 

 divide

EXPLANATION

Modeled receptors—Numbers are identifiers, 

summarized in table 1

 Constant-head boundary

 Saltwater boundary

 Coastal boundary (GHB)

 Estuary boundary (DRN)

 Freshwater boundary

 Pond

 Stream (STR)

Wastewater-disposal site

Simulated drawdown

Simulated path of advective transport and

 total traveltime, in years

 Less than 20 years

 21 to 50 years

 51 to 100 years

 More than 100 years

 Well

-4

 2

-2

41°42'

70°20' 70°16'

41°38'

-6 -4 -4

-6 -8

Site D

Site C

Site F

BWDF

Site B

Site B

105

150

106

38

120

151

150

30

15

104

74

73

555

666

151

191

73

72

74

89

71

47

43

46

96

41

94

28

11

26

30

52

50

58

15

32

32

92 91

89

72

70

-2

-2
18 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

12 percent of the wastewater discharged into a stream and 

entered the coastal waters as surface-water inflow. Generally, 

the direct discharge of sewage-derived nitrogen into coastal 

waters is more problematic than wastewater that reaches the 

coast as streamflow (Colman and others, 2004). The remaining 

18 percent of the wastewater was intercepted by production 

wells. About 13 percent of the wastewater flowed through 

freshwater ponds prior to later discharge into coastal waters 

or wells. Ponds are considered to be effective at removing 

nitrogen from the system. Also, pond ecosystems generally are 

phosphorus-limited and, therefore, are not as sensitive to the 

adverse effects of sewage-derived nitrogen.

Given the location and volume of wastewater disposal for 

this scenario, Barnstable Harbor (code 104, fig. 9) received the 

largest input of sewage-derived nitrogen expressed as volume 

of discharging wastewater, about 1.7 Mgal/d. Discharging 

wastewater exceeded 0.5 Mgal/d in Stewarts Creek (code 120) 

and Nantucket Sound (code 666) (fig. 10A). Generally, open 

coastal water bodies, such as Nantucket Sound, can more readily assimilate nitrogen owing to greater circulation of water 

and, as a result, are more favored as receptors of wastewater 

than are more constricted inland estuaries (Edward Eichner, 

Cape Cod Commission, oral commun., 2005). In addition to 

total volumes, the source of discharging wastewater is important to local and regional planners in evaluating the viability 

of wastewater-disposal sites. Barnstable Harbor received large 

volumes (more than 0.5 Mgal/d) of wastewater from both sites 

B and F, whereas most other estuaries received wastewater 

from a single disposal site (fig. 10A). Also, some sites, such 

as site B, contributed large amounts of wastewater to a single 

coastal receptor, whereas other sites, such as the BWDF and 

site F, contributed smaller volumes to a larger number of 

receptors (fig. 10A).

Simulated wastewater discharge from the BWDF as 

well as from proposed sites B and F flowed through freshwater ponds prior to discharging into downgradient wells and 

ecological receptors (fig. 10B). The largest input of wastewater, about 0.22 Mgal/d, was into Fawcetts Pond (code 52) 

downgradient from the BWDF. The wastewater-discharge rate 

exceeded 0.1 Mgal/d to three of the seven ponds that received 

wastewater in this scenario. No wastewater from these three 

sites discharged into a freshwater stream. About one-half of 

the wastewater from site C (about 0.66 Mgal/d) discharged 

into Phlashes Brook (code 15). No wastewater from site C 

passed through ponds or discharged directly at the coast.

The discharge of wastewater into a well can adversely 

affect water quality in the well and can be detrimental to 

human health. Eleven production wells received wastewater 

from either the BWDF or from a proposed site; some wastewater from all four sites discharged into a well. Wastewaterdischarge rates exceeded 0.1 Mgal/d to 4 of the 11 wells: 

codes 94, 71, 73, and 74 (fig 10C; table 1). Well 74 received 

the most wastewater, about 0.25 Mgal/d. Three of these four 

wells are downgradient from site C; the remaining well 94 is 

downgradient from the BWDF. Wastewater from sites B and 

Table 1. Wells, ecological receptors, and corresponding 

codes in modeled area of western and central Cape Cod, 

Massachusetts.

[Ecological receptors are ponds, streams, and coastal water bodies]

Name Code

Coastal water bodies

Hyannis Inner Harbor, Barnstable 73 

Mill Creek, Barnstable 74 

Great Marshes, Sandwich/Barnstable Harbor 104

Maraspin Creek, Barnstable 105

Mill/Short Wharf,Barnstable 106

Stewarts Creek, Barnstable 120

Upper Bumps River 150 

Lewis Bay Proper 151 

Coast Cape Cod Bay 555

Coast, Nantucket Sound 666

Ponds

Hinckley Pond 11

Israel Pond 26

Lamson Pond 28

Wequaquet Lake 30

Long Pond 50

Fawcetts Pond 52

Simmons Pond 58

Streams

Phlashes Brook 15

Wells

HIGGINS CROW #3a,Yarmouth 41

HIGGINS CROW #1,Yarmouth 43

RT 132 WELL #3,BFD 46

Rt. 132 Well #4,BFD 47

MAHER #3,BWC 71

MAHER #1a,BWC 72

MAHER #1,BWC 73

MAHER #2,BWC 74

Straightway 2,BWC 89

SIMMONS POND W,BWC 94

HYANNISPORT W,BWC 96
Simulation of Nitrate Transport 19

Figure 10. Model-estimated wastewater discharge rates into (A) coastal water bodies, (B) ponds, streams, and (C) production wells 

for a hypothetical wastewater-disposal scenario, Barnstable and Yarmouth, Cape Cod, Massachusetts.

FRESHWATER-RECEPTOR CODE

Ponds Streams

Site B Site F Site C

A

B

0

0

0

0.2

0.4

0.6

0.8

1.0

WASTEWATER DISCHARGE, IN MILLION GALLONS PER DAY

0.2

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

C

COASTAL-RECEPTOR CODE

PRODUCTION-WELL CODE

89 94 96 47 46 71 72 73 74 41 43

52 58 11 30 50 26 28 15

73 120 151 666 104 150 666 74 104 105 106 151 555

BWDF

BWDF

BWDF

Site B

Site B

Site F

Site C Site F
20 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

F discharged into two wells each; however, the volumes of 

discharge were low (0.02 Mgal/d or less).

The use of numerical models to estimate wastewaterloading rates to wells and natural receptors provides several 

benefits in evaluating wastewater-disposal scenarios. Modeling analyses can differentiate wastewater that discharges 

into receptors of greater concern, such as inland estuaries 

and wells, from wastewater that discharges into receptors of 

less concern, such as ponds, streams, or open coastal waters. 

Modeling analyses with a sufficient level of detail also can 

help determine the viability of proposed wastewater-disposal 

sites on the basis of possible ecological effects on sensitive 

receptors.

Advective Transport in Different Hydrologic 

Settings

Patterns of advective transport would be expected to 

differ in different hydrologic settings, particularly on the 

basis of the location of a potential disposal site relative to 

ground-water divides and hydrologic boundaries. In areas near 

divides, hydraulic gradients are small, and flow can diverge 

in different directions. Also, vertical flow is enhanced and 

velocities are smaller than in downgradient areas closer to 

hydrologic boundaries, where hydraulic gradients are steeper 

and flow directions more constrained.

Flow is outward from the regional divide towards 

hydrologic boundaries. Water enters coastal estuaries either as 

direct ground-water discharge or as streamflow. Some water 

is extracted for water supply, and recharge within contributing 

areas to ponds flows through the ponds prior to discharging 

into estuaries, streams, or wells (fig. 8A). Traveltimes, which 

are defined as the total transport time from recharge at the 

water table to discharge at a hydrologic boundary, are important considerations in the evaluation of wastewater disposal 

because nitrogen attenuation is assumed to increase with 

longer traveltimes (Colman and others, 2004; Valiela and 

others, 2000). Wastewater disposal that contributes water to 

an estuary with traveltimes longer than 10 years is assumed 

to pose less ecological risk than disposal of wastewater with 

traveltimes shorter than 10 years (School of Marine Science 

and Technology, 2005). Traveltimes to the coast are largest 

(longer than 100 years) for recharge that originates near the 

divide because velocities are smaller there, and these areas are 

farther from the discharge locations (fig. 8A); traveltimes for 

recharge that originates close to discharge locations are small 

(less than 1 year). Traveltimes to coastal discharge boundaries 

are shorter than 10 years in about 58 percent of the total watertable area that contributes water to the discharge boundaries. 

For ponds and streams, the percentages of contributing area 

for which recharge traveltimes are shorter than 10 years are 

82 and 72 percent, respectively. About 57 percent of water 

captured by production wells traveled less than 10 years from 

the recharge area at the water table.

Vertical flow paths and subsurface traveltimes under 

unstressed conditions (no pumping or return flow) are shown 

in figure 8B. In areas close to the divide, vertical gradients 

are enhanced relative to horizontal gradients, and advective 

flow paths are more downward than areas nearer the coast 

where vertical gradients are lower and flow paths are more 

horizontal. Flow is locally upward at coastal-discharge 

locations. This general pattern differs in the vicinity of local 

hydrologic features, such as ponds, streams, and wells; some 

perturbations are related to changes in lithology within the 

simulated aquifer. The particle paths shown in figure 8B 

were weighted with equal volumes of recharge water so that 

the volumes around the particles are analogous to stream 

tubes, and therefore, qualitatively represent ground-water 

fluxes with depth. Fluxes and velocities are higher nearer the 

coast as illustrated by the larger density of particles in those 

areas; conversely, fluxes and velocities near the ground-water 

divide are lower. In some areas, such as immediately below 

the divide, fluxes are small, and solutes move more slowly 

through that part of the aquifer; this pattern is reflected in the 

regional distribution of traveltimes (fig. 8A); long traveltimes 

(longer than 100 years) for recharge that originates near the 

regional ground-water divide are a result of transport through 

deep areas of the aquifer with small ground-water fluxes.

In the hypothetical scenario, wastewater disposal at three 

sites was on or close to the regional ground-water divide 

(BWDF and sites B and F) and at one site farther from the 

divide (site C) (fig. 9A). Hydraulic conditions at the three 

sites were reflected in the advective flow paths emanating 

from the source areas. Along the ground-water divide at sites 

B and F, simulated particle tracks fanned out in a wide arc of 

flow directions (fig. 9A). At site B, particle paths were both 

to the north into the Barnstable Harbor (code 104) estuarine 

system as well as to the south into Wequaquet Lake (code 30) 

and the coastal boundaries near Nantucket Sound (code 666). 

Similarly, advective flow from site F discharged to coastal 

boundaries to the north and south and to production wells 

to the east. In contrast, flow paths from site C had a smaller 

range of flow directions; wastewater discharged directly into 

Phlashes Brook (code 15) and into production wells directly 

downgradient from the site. Flow paths from the BWDF were 

to the south, but particle paths fanned out in a wider range of 

flow directions than particles from site C owing to the proximity of the regional ground-water divide to the BWDF source 

area (fig. 9A).

Traveltimes to coastal receptors differed by source area 

and receptor (fig. 11A). Median traveltimes from the BWDF 

source area to coastal receptors ranged from 22.7 years for 

Stewarts Creek (code 120) to 133 years for Nantucket Sound 

(code 666) (fig. 11A). The difference in traveltimes is a 

function of the setting of these hydrologic boundaries in the 

regional flow system. Stewarts Creek is an inland estuary that 

intercepts shallow ground water from a part of the aquifer with 

higher fluxes and velocities. Conversely, Nantucket Sound is 

the outer boundary of the flow system and captures ground 
Simulation of Nitrate Transport 21

Figure 11. Model-estimated traveltimes to (A) coastal water bodies and (B) production wells, Barnstable and Yarmouth, Cape Cod, 

Massachusetts.

N = 30

N = 1 N = 6

N = 2

N = 2

N = 30 N = 3

N = 40

N = 13 N = 15 N = 33

N = 4

Site B Site F

Median

90th percentile

10th percentile

25th percentile

75th percentile

Outlier

N = 13 Number of estimates

0

100

200

300

400

500

600

BWDF

BWDF

PRODUCTION-WELL CODE

TRAVELTIME, IN YEARS

0

400

800

1,200

1,600

2,000

2,400

Site C

COASTAL-RECEPTOR CODE

N = 259

N = 467

N = 5

N = 33

N = 443

N = 11

N = 446

N = 7

N = 255

N = 14

N = 9

N = 133

N = 30

A

B

Site B Site F

89 94 96 47 46 71 72 73 74 41 43

73 120 151 666 104 150 666 74 104 105 106 151 555

EXPLANATION
22 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

water from deeper in the flow system; this ground water 

recharged the aquifer through a part of the source area closer 

to the divide and traveled deeper in the aquifer at a slower 

velocity than water recharged near the coast (fig. 3B). Stewarts 

Creek and Nantucket Sound received similar volumes of 

wastewater from the BWDF; about 0.52 and 0.58 Mgal/d, 

respectively (fig. 10A). The traveltimes, however, indicated 

that Stewarts Creek may be more adversely affected by 

sewage-derived nitrate than Nantucket Sound owing in part 

to shorter traveltimes and, presumably, less subsurface nitrate 

attenuation.

Most wastewater from sites B and F was transported to 

the north and discharged into Barnstable Harbor (code 104) 

(fig. 10A). Median traveltimes for wastewater discharged 

from the two source areas to Barnstable Harbor were 26.8 and 

47.4 years, respectively. Both sites were at the ground-water 

divide and flow fanned out in a number of different directions 

(fig. 9A). As a result, traveltimes of wastewater from sites B 

and F to individual receptors also differed widely (fig. 11A). 

Wastewater from site B discharging into Nantucket Sound 

(code 666) had a median traveltime of more than 600 years; 

wastewater from site F discharging into Lewis Bay (code 151) 

had a median traveltime of about 900 years. In both cases, the 

wastewater flowed very deep into the system and discharged 

along the south shore. Most receptors had outliers with very 

large traveltimes (longer than 1000 years); these outliers represent particles that moved deeply through near-stagnant parts of 

the aquifer with little ground-water flux and slow velocities.

Traveltimes to production wells also differed by source 

area and individual well. Source areas at the BWDF and site 

C discharged large volumes of wastewater to wells (fig. 10C). 

For wastewater originating at the BWDF, median traveltimes 

ranged from 30.3 years (well 89) to 88 years (well 96) 

(fig. 11B). At site C, which discharged nearly half of its generated wastewater into four production wells, median traveltimes 

were all shorter than 5 years, owing to the shallower flow 

paths and faster velocities in that part of the aquifer as well as 

to the proximity of the wells to the source area (fig. 11B). As 

an example of the importance of traveltimes, well 94 received 

nearly 0.2 Mgal/d of wastewater from the BWDF, similar to 

the volumes received by wells 71 and 74 downgradient from 

site C; however, the median traveltime of wastewater received 

in well 94 is about 40 years compared to fewer than 5 years 

for wells 71 and 74 (fig. 11B). Based on the assumption that 

nitrate is attentuated over long transport times, the results of 

this simulation suggest that more nitrate would be attenuated 

prior to wastewater discharge into wells downgradient from 

the BWDF source area, and that wells downgradient from site 

C would be more adversely affected by wastewater disposal. 

Only small volumes of wastewater originating from sites B 

and F were captured by production wells. Traveltimes ranged 

from 86 to more than 400 years, indicating that the wells capture wastewater that had traveled along deep flow paths from 

the source areas.

Ground-water fluxes and velocities depend on location 

within the aquifer, and wastewater discharged in different 

hydrologic settings would be transported at different velocities. To generate a three-dimensional image of the velocity 

field within the aquifer, velocity vectors were computed for 

each segment of each particle track from the four simulated 

source areas to the discharge areas (fig. 12A). At the BWDF 

and at sites B and F near the ground-water divide, median 

velocities were less than 0.5 ft/d. At site C farther from the 

divide and in a shallower, more horizontal flow regime, the 

median velocity was 1.7 ft/d. Median total traveltimes of 

wastewater originating from the three source areas near the 

divide (the BWDF and sites B and F) were 66, 31, and 48 

years, respectively, whereas the median traveltime at site C 

farther from the divide was less than 10 years (fig. 12B).

Vertical flow paths along two columns in the model 

from sites in two different hydrologic settings—sites B and 

C—are shown in figure 13; the sites are at similar distances 

from hydrologic boundaries and receive the same volume of 

wastewater. For these simulations, no pumping was simulated 

at wells 71, 72, 73, 74, 46, and 47 so that advective transport 

would proceed to natural hydrologic boundaries from both 

sites; it should be noted that these stresses were different than 

in the scenario shown in figure 9, and therefore, flow lines in 

the two scenarios would not coincide. Flow paths from site B 

at the top of the ground-water divide were very deep—to near 

the bottom of the aquifer—and large volumes of the aquifer 

received wastewater. Flow was predominantly to the north; 

however, some flow was also to the south into Wequaquet 

Lake. Wastewater entering the lake also discharged from the 

lake and proceeded to the coast, but it was assumed that any 

associated nitrogen was attenuated within the pond. Some 

particles underflowed the lake and discharged into coastal 

boundaries along the south shore, but only after traveltimes 

longer than 100 years (fig. 13A).

From site C, flow paths were shallow and a smaller 

volume of aquifer received wastewater. Traveltimes, however, 

were significantly shorter (fig. 13B). Traveltimes of wastewater discharged from site B were longer and more variable 

than traveltimes of wastewater from site C; the minimum, 

maximum, and median traveltimes from site B were about 8, 

2095, and 32 years, respectively. The median traveltime of 

wastewater discharged from site C was about 10 years, and the 

range was about 3 to 72 years. This illustrates the importance 

of hydrologic setting in considering the downgradient proximity of receptors of concern and evaluating total transport time 

to the receptors from proposed wastewater-disposal sites.

Solute-Transport Simulations

The use of a solute-transport code has some advantages 

for evaluating the transport of wastewater-derived nitrate from 

source areas to receptors. Simulations with such a code can (1) 

estimate the distribution of nitrate in the aquifer, expressed as 

a concentration, (2) evaluate changes in nitrate concentrations 

over time at receptors, and (3) simulate the effects of 

mechanical dispersion on nitrate transport in the aquifer. 
Simulation of Nitrate Transport 23

Figure 12. (A) Ground-water velocities and (B) traveltimes for advective-transport paths from four wastewater-disposal sites for a 

hypothetical wastewater-disposal scenario, eastern Barnstable, Cape Cod, Massachusetts.

N = 320

N = 1408

N = 704

N = 10731

N = 55562

N = 640

DISPOSAL SITE

TRAVELTIME, IN YEARS

BWDF

Site B

Site C

Site F

0

100

200

300

400

500

EXPLANATION

Median

90th percentile

10th percentile

25th percentile

75th percentile

Outlier

N = 13 Number of estimates

0

1

2

3

4

5

A

N = 8554 

N = 10731

N = 55562

N = 21909

VELOCITY, IN FEET PER DAY

B
24 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

Figure 13. Vertical particle tracks representing advective-transport paths under the condition of no pumping downgradient from (

A) site C and (

B) site B near the Barnstable 

Wastewater Disposal Facility, Cape Cod, Massachusetts.

Barnstable Harbor

Wequaquet Lake

Barnstable Harbor

Upper Bumps River

Lewis Bay

Centerville  Harbor Long Pond

AB

Pond shadow

Site C

Site B

EXPLANATION

Simulated path of advective transport and total traveltime, in years 

 Less than 20 years 21 to 50 years51 to 100 yearsMore than 100 years

Particles flowing through ponds (Site B)

N

S

N

S

0

0 1 KILOMETER

1/2

1/2

1 MILE

0

-140-170-200-240-307

-40

-100-120

-40-70

0

-140-170-200-240-287

-40

-100-120

-40-70

VERTICAL EXAGGERATION X 5,000 

ALTITUDE, IN FEET RELATIVE TO NGVD 29
Simulation of Nitrate Transport 25

Particle tracking also can yield insight into the distribution of 

nitrate in the aquifer; however, the information from particletracking simulations is qualitative. As an example, the particle 

paths shown in figures 8 and 13 illustrate that a larger portion 

of the aquifer received wastewater from site B than from 

site C; wastewater from site B flowed across a wider arc of 

flow directions horizontally and much deeper in the aquifer 

vertically. The particle tracks, however, do not indicate the 

distributions of solute mass or of concentrations, which are of 

importance in evaluating the potential effects of wastewater 

disposal on ground-water quality. In addition, solute-transport 

models can more readily be used to evaluate time-varying 

concentrations and mass-loading rates at downgradient 

receptors. Although particle tracking can be used to quantify 

mass-loading rates at receptors by using mass-weighted 

particles, these methods often involve additional postmodeling 

analyses. The use of volume- or mass-weighted particles in 

modeling analyses for ground water on Cape Cod has been 

reported previously in Barlow (1997) and Walter and others 

(2002).

Solute-transport models can be used to represent the 

subsurface distribution of nitrate in the aquifer. Simulated 

nitrate concentrations downgradient from the BWDF at 

different horizontal sections through the nitrate plume are 

shown in figs. 14A-F. The concentration of nitrate in the 

discharged wastewater was simulated as 10 mg/L; after 

100 years of disposal, and assuming conservative transport, 

the maximum extent of nitrate concentration exceeding 1 mg/L 

was in model layer 18 (at 110 ft below NGVD 29) (fig. 14D). 

Maximum nitrate concentrations in this layer exceeded 8 mg/L 

at the coast. The distribution of nitrate is controlled in part 

by hydrologic boundaries. As an example, Fawcetts Pond 

(code 52) is a flow-through pond that focuses ground-water 

flow, and as a result, the patterns of nitrate concentrations are 

affected by the pond (figs. 14B–F).

Nitrate concentrations along vertical sections (column 

147 and row 304, fig. 14A) are shown in figures 14G and 

H. Concentrations exceeded 1 mg/L to an elevation of more 

than 300 ft below NGVD 29 (layer 24) beneath the disposal 

beds, near the bedrock surface; this pattern was consistent 

with the location of the BWDF near the regional ground-water 

divide where vertical flow paths were enhanced (fig. 14B). 

Nitrate concentrations along the coast exceeded 1 mg/L to a 

depth of about 250 ft below NGVD 29. Nitrate concentrations 

exceeding 1 mg/L were simulated at shallower depths in the 

aquifer near hydrologic boundaries than in the aquifer farther 

from the boundaries; this result indicates upward flow toward 

those boundaries. Examples are Fawcetts Pond (fig. 14G), 

Stewarts Creek, and Nantucket Sound (fig. 14H).

The solute-transport model also can be used to evaluate 

nitrate concentrations in receptors that receive wastewater 

discharge. Nitrate concentrations vary with time and with the 

location of a receptor relative to a source area (fig. 15). At 

Fawcetts Pond, concentrations increased to 8.1 mg/L within 

about 45 years after discharge from the BWDF; the pond 

is close to the BWDF source area (fig. 14A). At well 89, 

which is near the western edge of the nitrate plume (fig. 9), 

concentrations increased to about 1.3 mg/L within about 

65 years. Concentrations were constant after 100 years of 

transport; this result indicated that a sufficient transport time 

had elapsed to establish a steady-state condition with respect 

to nitrate concentration. At Stewarts Creek, concentrations 

increased to about 7.9 mg/L after 100 years of transport 

(fig. 15); although concentrations were still increasing after 

that time, the reduced rate of increase indicates a near-steadystate condition. At Nantucket Sound, concentrations increased 

to about 3.2 mg/L after 100 years of transport. The continuing 

increase in concentration indicates that a steady-state condition 

had not been reached and that concentrations would continue 

to increase for transport times beyond 100 years.

After 100 years of transport, the distribution of mass 

downgradient from the four wastewater-disposal sites is 

a result of advective transport patterns, and therefore the 

hydrologic setting of the disposal site. As an example, nitrate 

concentrations exceeded 0.1 mg/L over a large area downgradient from site F (fig. 16). Because site F is near the 

regional ground-water divide, discharge from this site generated a wide arc of advective flow paths from the site (figs. 9 

and 16). Beneath site B, also near the divide, concentrations 

exceeding 9 mg/L extended to near bedrock (the full saturated 

thickness of the aquifer) owing to enhanced vertical flow 

beneath the site (figs 13A and 17A). At site C, farther from the 

ground-water divide concentrations exceeded 0.1 mg/L over a 

smaller area than at sites B and F, reflecting more constrained 

advecive flow paths (figs. 9 and 16). Flow from this site was 

more horizontal and advective-transport paths were shallower; concentrations exceeding 9 mg/L extended to about 70 

ft below NGVD 29 (or the upper 25 percent of the saturated 

thickness) after 100 years of transport (figs. 13B and 17B). 

Although advective-transport paths did not extend deeper 

than about 140 ft below NGVD 29 (or the upper 28 percent of 

the saturated thickness) (fig. 13B), simulated concentrations 

beneath site C exceeded 0.1 mg/L to the depth of bedrock 

(fig 17B). It is likely, however, that there is no nitrate mass 

at depth in the aquifer and that the simulated concentrations 

were caused by numerical dispersion. Vertical discretization is 

larger at depth in the model, and therefore numerical dispersion, which is a function of model-cell size, would be expected 

to be enhanced there. Although mass fluxes at depth in the 

aquifer are small, this simulated scenario illustrates the need 

to consider numerical dispersion when evaluating simulated 

concentrations.
26 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

Figure 14. (A-F) Simulated nitrate concentrations downgradient from the current wastewater-disposal facility for six horizontal 

sections and (G-H) simulated nitrate concentrations along vertical column 147 and row 304.

Water table (layer 1) 0 feet relative to NGVD 29 (layer 8) 55 feet below NGVD 29 (layer 13)

185 feet below NGVD 29 (layer 21) 265 feet below NGVD 29 (layer 23)

Row 304

Column 147

A B C 

70°18'59" 70°16'21"

41°40'33"

41°37'34"

41°40'33"

41°37'34"

70°18'59" 70°16'21" 70°18'59" 70°16'21"

Site C

Site C

BWDF

Coastal boundary (GHB)

Estuary boundary (DRN)

Saltwater boundary

Pond

Freshwater boundary

EXPLANATION

Modeled receptors

Site C

Site C

Site C

Site C

Simulated nitrate concentrations—

 Number is milligrams per liter

1 to 2

Greater than 4 to 5

Greater than 3 to 4

Greater than 2 to 3

Greater than 5 to 6

Greater than 6 to 7

Greater than 10

Greater than 8 to 9

Greater than 7 to 8

Wastewater-

 disposal site

Fawcetts Pond

Nantucket

Sound

Stewarts

 Creek

89

Fawcetts Pond

0

0 1 KILOMETER

1/2 1 MILE

1/2

110 feet below NGVD 29 (layer 18)

D E F

Well
Simulation of Nitrate Transport 27

Figure 14. (A-F) Simulated nitrate concentrations downgradient from the current wastewater-disposal facility for six horizontal 

sections and (G-H) simulated nitrate concentrations along vertical column 147 and row 304.—Continued

Simulated nitrate concentrations—

Number is milligrams per liter

1 to 2

Greater than 4 to 5

Greater than 3 to 4

Greater than 2 to 3

Greater than 5 to 6

Greater than 6 to 7

Greater than 10

Greater than 8 to 9

Greater than 7 to 8

EXPLANATION

Layer 18

Layer 21

Layer 23

Layer 13

Layer 8

Layer 1

Layer 18

Layer 21

Layer 23

Layer 13

Layer 8

Layer 1

Vertical section along column 147

BWDF Fawcetts

 Pond

Nantucket

 Sound

Nantucket

 Sound

Stewarts

Well 89 Creek

Vertical section along row 304

0

0 1 KILOMETER

1/2 1 MILE

1/2

B
28 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

0

2

4

6

8

10

0 20 40 60 80 100

TRANSPORT TIME, IN YEARS

SIMULATED NITRATE CONCENTRATION,

IN MILLIGRAMS PER LITER

Wastewater concentration = 10 milligrams per liter

Well 89

Nitrate concentration

Fawcetts Pond

Nantucket Sound

Stewarts Creek

Figure 15. Simulated time-varying nitrate concentrations at well 89, Fawcetts Pond, Stewarts Creek, and Nantucket Sound, Cape Cod, 

Massachusetts.
Simulation of Nitrate Transport 29

Figure 16. Simulated nitrate concentrations downgradient from five wastewater-disposal sites and simulated advective-flow paths in 

model layer 13 (50 to 60 ft below NGVD 29) near the Barnstable Wastewater-Disposal Facility, Cape Cod, Massachusetts.

EXPLANATION

Modeled receptors—Numbers are identifiers

 Constant-head boundary

 Saltwater boundary

 Coastal boundary (GHB)

 Estuary boundary (DRN)

 Freshwater boundary

 Pond

 Stream (STR)

0

0 1 KILOMETER

1/2

1/2

1 MILE

Simulated nitrate concentrations—Number is 

 milligrams per liter

0.1 to 1

Greater than 1 to 2

Greater than 4 to 5

Greater than 3 to 4

Greater than 2 to 3

Greater than 5 to 6

Greater than 6 to 7

Greater than 10

Greater than 8 to 9

Greater than 7 to 8

Wastewater-disposal site

Simulated path of advective transport

Well

Site B

41°42'

70°20' 70°16'

41°38'

Column 65 Column 233

Site D

Site C

Site F

BWDF

Site B
30 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer VERTICAL EXAGGERATION X 5,000 

0

0 1 KILOMETER

1/2

1/2

1 MILE

Barnstable Harbor

Wequaquet Lake

Barnstable Harbor

Upper Bumps  River

Lewis Bay

Centerville

Harbor

Long Pond

AB

Pond Shadow

Site C

Site B

EXPLANATION

Simulated nitrate concentrations— Number is milligrams per liter

Simulated path of  advective transport

0.1 to 1Greater than 1 to 2Greater than 2 to 3Greater than 3 to 4Greater than 4 to 5Greater than 5 to 6Greater than 6 to 7Greater than 7 to 8Greater than 8 to 9Greater than 10

0

-140-170-200-240-307

-40

-100-120

-40-70

0

-140-170-200-240-287

-40

-100-120

-40-70

Figure 17. Nitrate concentrations and particle tracks representing advective-transport paths downgradient from (

A) site B and (

B) site C near the Barnstable WastewaterDisposal Facility, Cape Cod, Massachusetts.
Simulation of Nitrate Transport 31

Comparison to Advective-Transport Simulations

Both solute-transport and particle-tracking methods 

can simulate the movement of a conservative solute, such as 

nitrate, through the aquifer. Particle-tracking methods use 

mass-weighted particles to simulate the movement of a solute 

and can estimate mass-loading rates at receptors through an 

accounting of particles (fig. 18A). This method requires no 

additional model development beyond development of the 

flow model but does require postmodeling analyses to convert 

particle track information into useful estimates of mass and 

concentration. In addition, representing mass by discrete 

particles nay not yield sufficient resolution for some scenarios, 

and the assumption that there is no dispersion or loss of mass 

during transport is not valid in most systems. Conversely, 

solute-transport methods explicitly simulate concentration and 

mass at all points within the model domain at each transport 

time step and can simulate transport processes, such mechanical dispersion and sorption; however, this approach requires 

additional model development and the assembly of additional 

input-data sets. In addition, numerical dispersion needs to 

be considered when interpreting model results. Particletracking results in complete conservation of mass along flow 

lines (fig. 18A), whereas solute-transport methods allow for 

spreading of mass between flow lines (transverse and vertical 

dispersion) (fig. 18B). Solute-transport methods also allow 

for the spread of mass in the direction of flow (longitudinal 

dispersion). Although mass is conserved for the entire system 

in the case of nonreactive solutes such as nitrate, the inclusion 

of dispersion in the transport process can result in changes in 

mass along flow lines (fig. 18B). The choice of method should 

be based on the types of questions being addressed and the 

level of resolution and accuracy needed.

Particle tracks from the four wastewater-disposal sites 

and simulated nitrate concentrations in model layer 13 (50 to 

60 ft below NGVD 29) are shown in figure 16; it should be 

noted that the comparison is between mass in a single model 

layer and a projection of three-dimensional particle tracks 

onto a two-dimensional map surface. The spatial distribution 

of particle tracks generally corresponds to the distribution 

of mass simulated by the solute-transport model (fig. 16). 

Thus, the two methods delineate the same general areas of 

the aquifer that receive wastewater-derived nitrate; however, 

the particle tracks intrinsically do not include any quantitative 

information regarding the distribution of mass. The agreement 

between simulated vertical particle tracks and nitrate mass 

near sites B and C is also close (figs. 17A and B). Both also 

delineate the same general volumes of the aquifer that receive 

Figure 18. Mass-loading rates estimated from (A), flow-based particle tracking and (B) solute-transport modeling at two downgradient 

production wells, Cape Cod, Massachusetts.

6 6

8

3 3

Well 1

Well 2

Mass-weighting rate (kilograms per day) Mass-loading rate (kilograms per day)

Dispersion Dispersion

0 8 10 8 0 0 8 10 8 0

Well 1

Well 2

Well 1—26 kilograms per day

Well 2—0 kilograms per day

Well 1—20 kilograms per day

Well 2—3 kilograms per day

Loads Loads

8

10

0 0

8

No Dispersion

No Dispersion

A B

Nitrate Source Nitrate Source
32 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

wastewater-derived nitrate. The spreading of concentrations between 0.1 and 1 mg/L below the particle tracks to the 

bedrock surface beneath site C was likely caused by numerical 

dispersion (fig. 17B).

Both particle-tracking and solute-transport methods 

can provide estimates of mass-loading rates into wells and 

ecological receptors. Flow-based particle-tracking methods 

calculate mass-loading rates by using a set of discrete particles 

to represent the load of wastewater being discharged into 

the aquifer. From the known concentration of nitrate in the 

wastewater, a mass can be assigned to each particle, and the 

fate of the nitrate can be determined by an accounting of 

particles that discharge into each receptor. If the distribution 

of particles representing the loading at the site is uniform, 

the load associated with each particle is determined by 

multiplying the concentration by the wastewater-discharge rate 

and dividing that result by the number of particles. Important 

considerations to make when estimating mass-loading rates by 

particle tracking include (1) the need to use a sufficiently large 

number of particles to represent the advective-flow system 

adequately and minimize the amount of mass represented 

by each particle and (2) the possibly invalid assumption that 

mass does not change along flow lines because of dispersion. 

Solute-transport methods calculate mass-loading rates into a 

receptor by multiplying the simulated concentration by the 

ground-water flux into the receptor.

An important concern in Barnstable and western 

Yarmouth, the most densely populated area on Cape Cod, is 

the potential for the discharge of wastewater-derived nitrate 

into production wells from one current and several potential 

wastewater-disposal sites. For the hypothetical scenario 

described earlier, the particle-tracking method simulated 

a range of mass-loading rates to the wells of 10.1 kg/d 

(well 74) to 0.1 kg/d (well 46) (fig. 19). Well 74 is directly 

PRODUCTION WELL

MASS LOADING OF NITRATE, IN KILOGRAMS PER DAY

Site B BWDF Site F

0

2

4

6

8

10

12

EXPLANATION

Solute transport

Particle tracking

Solute transport

Particle tracking

Solute transport

Particle tracking

Solute transport

Particle tracking

BWDF Site B Site C Site F

Site C

89 94 96 47 46 71 72 73 74 41 43

Figure 19. Nitrate loads estimated by flow-based particle tracking and solute transport modeling in wastewater intercepted by 11 

production wells, Cape Cod, Massachusetts.
Simulation of Nitrate Transport 33

downgradient from proposed site C and would receive about 

half of its water from wastewater discharge from the site 

(fig. 9). Well 46, to the east of Site B (fig. 9), is not in the 

part of the aquifer that would receive significant amounts of 

wastewater, as indicated by advective flow paths; as a result, 

the mass-loading rate is small. Mass-loading rates estimated 

by the solute-transport model ranged from 8.8 kg/d at well 74 

to 0.1 kg/d at wells 41, 43, and 46 (fig. 19). Advective flow 

paths indicate that wells 41 and 43, to the east of site F (fig. 9), 

are in a part of the aquifer that would receive only small 

amounts of wastewater.

Mass-loading estimates calculated by the two methods 

generally were in agreement (fig. 19). Estimates by particletracking methods generally exceeded estimates by the solutetransport method; this was the case for 9 of the 11 wells. This 

likely was the result of the complete conservation of mass 

along flow lines inherent in the particle-tracking method. The 

numerical dispersion created by the solute-transport model 

would spread mass spatially and thus decrease concentrations 

in some areas; particle tracking does not account for 

dispersion and the resulting local change in mass. The results 

calculated by the two methods agreed best for the wells 

downgradient of site C. The median traveltimes for wastewater 

intercepted by these wells were small (about 5 years), and it 

is likely that, given the short transport time, little spreading 

of mass occurred prior to discharge into the wells. Thus, the 

assumption of complete conservation of mass along flow 

lines was not as important for the simulation of wastewater 

transport to these wells. The largest discrepancy was for 

well 94 where mass-loading estimates by particle-tracking and 

solute-transport methods were 7.3 and 4.8 kg/d, respectively. 

The well is directly downgradient from the BWDF (fig. 9), 

and the mass-loading rate into the well was represented by a 

larger count of particles. In this case, the discrepancy likely 

was a result of the presence of Simmons Pond (code 58) just 

upgradient from the well (fig. 9). In the particle-tracking 

simulation, particles flowed through the pond and into the 

well with no loss of mass along flow lines represented by the 

movement of particles. In the solute-transport simulation, 

mass decreased locally within the pond owing to enhanced 

dilution and dispersion. These processes resulted in lower 

concentrations and, therefore, lower mass-loading estimates 

downgradient at the well.

For wells 72 and 47, mass-loading estimates by the 

solute-transport model were higher. Well 72 is downgradient 

of site C in close proximity to three other wells. Because the 

well is screened 20 ft directly below well 73, the particletracking method may not represent advective flow paths with 

sufficient resolution to produce reliable estimates for each 

well. If the two wells are combined, mass-loading estimates by 

the particle-tracking and solute-transport methods are similar: 5.3 and 4.8 kg/d, respectively. Well 47 is east of site B in 

an area that would receive a small portion of the wastewater 

(fig. 9); the mass-loading rate to this well is represented by a 

small count of particles. In these cases, particle-tracking methods based on particle densities that are practicable may not 

represent advective flow paths with enough resolution, which 

would be particularly important in areas of strongly diverging flow paths, such as near sites B and F (fig. 9). As a result, 

particle-tracking methods may have limitations for wells with 

small mass-loading rates in strongly diverging flow fields 

given particle densities that are practicable.

Mass-loading rates estimated from mass-weighted 

particles represent a steady-state condition in which 

concentrations have reached a maximum and do not change 

further over time. In some cases, it is useful to determine the 

amount of time before wastewater disposal affects a well or 

ecological receptor. This determination is done by quantifying 

how concentrations change with time prior to reaching a 

steady-state condition. Solute-transport simulations yield 

nitrate concentrations for each time step directly from model 

output. At wells 71 and 74, concentrations began to increase 

within about 2 years of transport and reached a final steady 

state after about 13 and 24 years, respectively (figs. 20A and 

20B). Because wells 71 and 74 are about 2,100 ft directly 

downgradient from site C, they would receive a large amount 

of wastewater, and as a result, concentrations would reach a 

steady-state condition quickly. Well 94 is at a greater distance, 

about 7,800 ft, downgradient from the BWDF (fig. 9). Nitrate 

concentrations in wastewater intercepted by this well did not 

increase until after about 20 years of simulated transport from 

the BWDF and then increased gradually to a near-steady-state 

condition after about 100 years (fig. 20C). These results are 

consistent with the location of the well downgradient from 

the disposal site and near the edge of the simulated nitrate 

plume. At well 47 to the east of site B, concentrations began 

to increase after about 15 years and were still increasing 

steadily after 100 years of transport (fig. 20). These results 

also were consistent with traveltimes estimated by the flow 

model (fig. 11). The median traveltime for wastewater arriving 

at wells 71 and 74 is less than 10 years; thus, concentrations 

would increase to a steady-state condition more quickly at 

these wells than concentrations in wastewater intercepted by 

wells 94 and 47, which capture flow from deeper in the aquifer 

and have median traveltimes of 45 and 90 years, respectively 

(fig. 11).

Particle-tracking methods also can estimate timevarying concentrations by accounting for the traveltimes of 

individual mass-weighted particles discharging into each 

receptor; these traveltimes and time-specific particle counts 

can be converted to concentrations by dividing the calculated 

mass-loading rate by the water flux through the well. Because 

the process involves several postmodeling analyses, it is not 

as straightforward as the direct simulation of concentrations 

by the solute-transport method. Time-varying concentrations 

simulated by the particle-tracking method for wells 71 and 74 

were similar to the concentrations simulated by the solutetransport model (figs. 20A and 20B). For well 94, however, 

the steady-state concentration estimated by the two methods 

differed by about 2.5 mg/L owing to the assumption inherent 

in the particle-tracking method that no loss of mass occurs 

along flow lines. This assumption probably affected the results 
34 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

for well 94 because the well is far from the source area and 

downgradient of a flow-through pond; both of these factors 

likely caused a local decrease in mass due to enhanced dispersion. For well 47, the two methods gave similar results, but 

the solute-transport method predicted larger concentrations 

than the particle-tracking method. In this case, estimates of 

time-varying concentrations were based on a small number of 

particles; the results suggest that the particle-tracking method 

does not sufficiently define the advective-flow patterns to the 

well given particle densities that are practicable. Although the 

particle-tracking method does yield similar time-varying concentrations in some cases, the additional effort and possibly 

invalid assumptions involved in the approach make the choice 

of solute-transport methods preferable.

The effect of dispersion on mass-loading rates and how 

those effects differ spatially are shown conceptually in figures 

18A and B, in which a flow field is represented as discrete 

flow lines. Figure 18A illustrates a hypothetical particletracking scenario in which mass-loading rates were estimated 

at two wells downgradient from a continuous source of nitrate. 

Particles within the source area are weighted with 8 to 10 

kg/d of nitrate, and particles from outside the source area are 

not weighted with any mass. An accounting of particles at the 

two wells shows that well 1 downgradient from the source 

receives a total of 26 kg/d of mass, whereas well 2 to the right 

of the source area receives no mass. If transverse dispersion 

allows for mass to spread between adjacent flow lines, mass 

decreases along the center flow line and increases along flow 

lines near the edges. The resulting mass-loading rates at wells 

1 and 2 are 20 and 3 kg/d, respectively (fig. 18B). Although 

mass is shifted within the system, the total mass is conserved 

(26 kg/d). This simplified example illustrates how dispersion 

can result in different distributions of estimated concentrations 

and mass within an aquifer.

0

1

2

3

4

5

6

TIME, IN YEARS

 NITRATE CONCENTRATION, IN MILLIGRAMS PER LITER

Solute transport

Particle tracking

B

D

A

Steady-state concentration 

Well-71 Well-74

Well-94 Well-47

C

0

1

2

3

4

5

6

0 20 40 60 80 100 0 20 40 60 80 100

Figure 20. Time-varying and steady-state nitrate concentrations in wastewater intercepted by four production wells, calculated by 

solute-transport and particle-tracking methods, Cape Cod, Massachusetts.
Simulation of Nitrate Transport 35

Effects of Dispersion

Advection, which refers to the movement of a solute with 

ground-water flow, is the dominant component of transport in 

most sandy coastal aquifers, including the Cape Cod aquifer 

system, where recharge rates and flow velocities are high; 

dispersion, however, is an additional component of transport 

that can affect the movement of nitrate through these aquifers. 

Dispersion refers to the differential movement of a solute 

through the aquifer as a result of heterogeneities within the 

aquifer matrix and, like heterogeneity, is scale dependent and 

is expressed in units of length. In a perfectly homogenous 

system, dispersivity is near zero; as the degree of heterogeneity increases, dispersivity increases. For one-dimensional 

transport through a homogenous sediment column, outflowing 

concentrations would increase to the initial inflowing concentration instantaneously at a time determined by the water 

velocity. If the sediments are heterogeneous, dispersion would 

cause outflowing concentrations to increase gradually to the 

inflowing concentration; the rate of the increase is a function 

of the dispersivity of the sediments. Longitudinal dispersion 

refers to a spreading of mass in the direction of ground-water 

flow; transverse and vertical dispersion refer to a spreading of 

mass orthogonal to the ground-water-flow direction. Longitudinal dispersion is typically much larger than the other two 

components.

Estimates of dispersivity increase with increasing 

transport distance. Horizontal transport distance, defined as 

the distance between the waste-disposal sites and discharge 

locations, ranged from about 1 to 3 mi for the nitrate plumes 

simulated in this analysis. This range of distances is referred 

to as plume-scale transport distances. On the basis of an 

empirical relationship derived from literature values (Luckner 

and Schestakow, 1991), longitudinal dispersivity ranges from 

42 ft to 55 ft for plume-scale transport distances of 7,500 

and 13,500 ft, respectively. Spitz and Moreno (1996) suggest 

a longitudinal dispersivity of about 90 ft for a plume-scale 

transport distance of about 13,000 ft. Longitudinal, transverse, 

and vertical dispersivities of 35, 3.5, and 0.35 ft best matched 

plume geometries at the MMR on western Cape Cod (U.S. 

Air Force Center for Environmental Excellence, 2000). The 

approximate numerical dispersivity generated by the solutetransport model—as estimated from equation (1)—is about 

55 ft. This similarity suggests that numerical dispersion alone, 

a modeling artifact, may reasonably represent the general scale 

of longitudinal dispersion over plume-scale transport distances 

in the aquifer.

Dispersion can affect simulated concentrations and 

estimated mass-loading rates of nitrate at wells and ecological 

receptors that receive wastewater. For a mass of solute 

released continuously from a source over time, the effect of 

dispersion on concentrations at a given time and receptor 

would depend on the location of the receptor relative to 

the source, including both the downgradient distance from 

the source and the location of the receptor relative to the 

center of advective transport from the source. The effects 

of dispersion also would differ if, at a given time, steadystate concentrations had been reached at a given location. 

For a continuous source of solute, the effect of longitudinal 

dispersion generally would be to decrease the time of transport 

to a receptor and thereby increase the rate at which concentrations increase prior to reaching steady state. The effects of 

transverse and vertical dispersion would be to spread mass 

over a larger volume of aquifer and thus to decrease steadystate concentrations at receptors directly downgradient from 

source areas (near the centers of nitrate plumes) and increase 

steady-state concentrations in areas near the edges of the 

plumes.

The effect of dispersion on time-varying concentrations 

in four wells is shown in figure 21. Longitudinal dispersivities 

of 20, 50, 100, 200, and 500 ft were simulated; these values 

represent mechanical dispersion in addition to numerical 

dispersion (about 50 ft). In all cases, transverse and vertical 

dispersivities were specified as being one and two orders of 

magnitude lower, respectively, than longitudinal dispersivity. 

In general, the results show that the simulated dispersivity 

does not appreciably affect the simulated concentration. 

It should be noted that a dispersivity of zero represents 

numerical dispersion only. At wells 71 and 74, concentrations 

during 100 years of transport did not change substantially with 

increasing dispersivity (figs 21A and B). For dispersivities 

between 0 and 200 ft, the steady-state concentrations in 

wells 71 and 74 differed only by about 0.04 and 0.08 mg/L, 

respectively. For a dispersivity of 500 ft, steady-state 

concentrations were about 0.2 and 0.3 mg/L lower than 

concentrations simulated for the other dispersivities in wells 

71 and 74, respectively. In these cases, sufficient time and 

transport distance may not have elapsed to allow dispersion 

to appreciably affect concentrations simulated at the wells. 

For well 94, nitrate concentrations for dispersivities between 

0 and 100 ft differed by about 0.46 mg/L after 100 years of 

transport, a difference of about 33 percent (fig. 21C). Timevarying concentrations in well 47 increased with increasing 

dispersivity for about the first 40 years of transport. For 

transport times of about 50 to 100 years, concentrations 

for dispersivities between 0 and 200 ft were similar. For 

dispersivities of less than 200 ft, concentrations in wells 94 

and 47 were still increasing after 100 years of transport.

For the largest dispersivity simulated (500 ft), steadystate concentrations in wastewater in wells 71 and 74 

were about 0.2 and 0.3 mg/L lower, respectively, than 

concentrations for dispersivities of less than 200 ft. The lower 

steady-state concentrations likely were a result of an enhanced 

spreading of mass at the larger transverse and vertical 

dispersivities (50 and 5 ft). For wells 94 and 47, a steady-state 

concentration was reached when a dispersivity of 500 ft was 

simulated (fig. 21C and D). For well 94, concentrations for 

longitudinal, transverse, and vertical dispersivities of 500, 50, 

and 5 ft, respectively, were similar to those for dispersivities 

of 200, 20, and 2 ft for about the first 50 years of transport. 

Between transport times of 50 and 100 years, concentrations 

for a longitudinal dispersivity of 500 ft neared a steady-state 
36 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

condition, whereas concentrations for the lower dispersivities 

continued to increase (fig. 21C). This effect is illustrated 

more clearly for well 47. Concentrations for a longitudinal 

dispersivity of 500 ft were larger than concentrations for 

lower dispersivities for the first 50 years of transport. After 

50 years of transport, concentrations assuming a dispersivity 

of 500 ft reached steady state, while concentrations for 

lower dispersivities were still increasing; the concentration 

assuming a dispersivity of 500 ft was about 0.4 mg/L lower 

than concentrations assuming lower dispersivities after 100 

years of transport (fig. 21D). Simulation results indicated 

that the effects of longitudinal dispersivity caused earlytime concentrations to increase with increasing dispersivity; 

however, at very high values of dispersivity, concentrations 

reached steady state. By that time in the transport process, 

the enhanced transverse and vertical spreading of mass may 

have depleted the mass transported along the central path 

directly downgradient from the source and thus lowered 

the concentration below the steady-state value reached at 

the lower dispersivities. The results indicate that the effects 

of dispersion on nitrate concentrations are enhanced in 

wells (1) where sufficient transport times have elapsed, 

(2) near the edges of a plume, or (3) downgradient from 

flow-through ponds.

TIME, IN YEARS

Well-71 Well-74

Well-94 Well-47

 NITRATE CONCENTRATION, IN MILLIGRAMS PER LITER 0

1

2

3

4

5

0

1

2

3

4

5

Numerical (0)

+ 20

Longitudinal dispersivity,

 in feet

+ 200

+ 100

+ 50

+ 500

0 20 40 60 80 100 0 20 40 60 80 100

B

D

A

C

Figure 21. Time-varying nitrate concentrations in wastewater intercepted by four production wells under simulated dispersivities of 0, 

20, 50, 100, 200, and 500 ft., Cape Cod, Massachusetts.
Summary 37

Limitations of Analysis

Several limitations apply to estimates of mass loading 

by particle-tracking (flow-based) and solute-transport 

methods. The analyses presented in this report were done 

with deterministic ground-water-flow models in which 

hydraulic conductivities, recharge, and boundary leakances 

were determined by trial-and-error calibration to observations 

of heads, streamflows, and advective flow, as indicated by 

contaminant plumes. The model is assumed to represent 

the flow system adequately because (1) the distribution of 

hydraulic conductivity is based on a model that incorporates 

an understanding of the depositional history of the glacial 

sediments in the region, (2) values of hydraulic conductivity, 

boundary leakance, and recharge are reasonable values 

consistent with previous investigations, and (3) the simulated 

results reasonably match observations of heads, streamflows, 

and advective transport (Walter and Whealan, 2004). The flow 

model produces a simulated flow field that is nonunique in 

the sense that different model configurations could produce 

similar head distributions. Particle-tracking results would be 

expected to differ as aquifer properties and recharge were 

changed; however, if the resulting head distribution was 

reasonably calibrated to observed conditions, particle tracks 

likely would be similar for different combinations of aquifer 

properties and recharge. In addition, flow lines in the aquifer 

are strongly affected by hydraulic boundaries, such as ponds, 

streams, and coastal waters; because boundary geometries 

are static, however, simulated particle tracks and estimates of 

advective flow likely would be similar for a range of model 

inputs.

Although advective-flow paths are likely to be reasonably 

robust for a reasonable range of model inputs, simulated 

traveltimes can be affected by several factors. If hydraulic 

conductivity and recharge rates were increased in a correlated 

manner relative to their values in a baseline simulation to 

produce a similar calibrated head distribution, ground-water 

fluxes would increase and traveltimes would decrease. 

Another factor affecting traveltimes is aquifer porosity. The 

uniform value of 0.35 used in these analyses was based on 

field-scale tracer tests on western Cape Cod (Garabedian and 

others, 1985); the porosity of glacial sand and gravel likely 

ranges between 0.25 to 0.4 (Freeze and Cherry, 1979). The 

linear velocity of ground water is inversely proportional to 

porosity; decreasing the porosity would increase velocity and 

decrease simulated traveltimes.

Additional considerations are necessary when particletracking methods are used to estimate concentrations and 

mass; uncertainties in simulated results can arise from the 

representation of volumes of water or mass by too small 

a number of discrete particles. If the number of particles 

is small, the particles may not represent the flow field at a 

resolution adequate to estimate traveltimes to a receptor. Even 

if adding particles does change the estimated mean traveltime 

to a receptor, the means should generally be similar because 

the position of the receptor in the flow system would have the 

most control on the age of water intercepted by the receptor. 

It is advisable to choose the largest number of particles that is 

practicable with the use of particle-tracking methods.

Another assumption in these analyses is that of steadystate flow based on constant recharge. Although actual 

recharge varies both monthly and from year to year, hydraulicgradient directions in the aquifer generally do not vary over 

time owing to the small time scale of recharge variability 

(months) compared to the much larger time scale of advective 

transport (multiple years). Walter and Masterson (2003) 

compared particle tracks on northwest Cape Cod through 

both steady-state and transient flow fields and found that 

the particle tracks were nearly identical. The results of this 

comparison indicated that it was reasonable to simulate 

steady-state advective transport even though recharge was 

known to vary through time. Ground-water withdrawals also 

change seasonally and over multiple years. These seasonal 

changes and longer-term trends also would be expected to 

alter hydraulic gradients; however, these changes likely would 

be limited to areas close to the wells.

The solute-transport methods used in these analyses also 

rely on intercell flows produced by the ground-water-flow 

model. As a result, estimated concentrations, mass-loading 

rates, and ground-water velocities are susceptible to the same 

uncertainties as the same quantities calculated from particle 

tracks produced by the flow model. Additional uncertainties 

include (1) the effects of dispersion on the transport of mass, 

(2) the possibility of processes affecting nitrate transport that 

are not accounted for in the solute-transport model, and (3) 

estimated nitrogen inputs, which will depend on nitrogentreatment levels at specific wastewater-treatment plants in the 

proposed source areas.

Summary

Cape Cod, a rapidly urbanizing region in southeastern Massachusetts, faces several water-resource issues. The 

unconsolidated sandy aquifer underlying the region is the 

sole source of potable water to local communities and is the 

primary receptor of the region’s wastewater. Ground-water 

discharge from the aquifer supports ecosystems in freshwater 

ponds and streams and in coastal waters.The disposal of 

municipal wastewater from the region into the aquifer has 

created the potential for adverse effects of wastewater disposal 

on production wells and ecological receptors. Nitrogen, in 

the form of nitrate, is a wastewater constituent of particular 

concern. Marine and estuarine ecosystems are sensitive to 

inputs of wastewater-derived nitrogen because excess nitrogen 

in marine waters can cause algal blooms, degradation of water 

clarity and, ultimately, loss of critical habitats, such as seagrass beds. In areas where wastewater-disposal sites are near 

production wells, the discharge of nitrate into the wells can be 

harmful to human health.

Many of the communities in the region are in the 

process of developing wastewater-management strategies 
38 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

that will allow for the sustainable use of the aquifer as both 

a source of potable water and a receptacle of wastewater and 

will maintain healthy aquatic ecosystems. To support these 

efforts, the Massachusetts Department of Environmental 

Protection is currently (2006) developing estimates of the total 

maximum daily loads of nitrogen into individual estuaries and 

coastal waters. In 2003, the U.S. Geological Survey began 

a cooperative investigation with Barnstable County and the 

Cape Cod Commission to describe the effects of wastewater 

disposal in the context of the regional-flow system. An 

important component of the investigation was the use of 

numerical models to simulate the transport of nitrate from 

areas of centralized wastewater disposal to wells and natural 

receptors, such as ponds, streams, and coastal waters. Possible 

wastewater-disposal scenarios, including sewering and 

centralized disposal, from individual towns were simulated 

using regional flow models of central and western Cape Cod. 

The dual purpose of this modeling support was to (1) assist 

communities in the initial evaluation of possible wastewaterdisposal scenarios and (2) develop a regionally applicable 

understanding of wastewater disposal in different hydrologic 

settings.

Numerical models of ground-water flow can be of use 

in the evaluation of wastewater-disposal scenarios by simulating both the hydraulic effects of wastewater disposal on 

water levels and streamflow and the advective transport of 

wastewater-derived nitrate through the aquifer. Of particular 

concern to local and regional wastewater planners are the 

potential mass-loading rates of nitrate into wells and ecological receptors—ponds, streams, and coastal waters. Particletracking methods use model-calculated cell-by-cell flows to 

track particles of water through the subsurface. These methods 

can quantify wastewater discharge to receptors by an accounting of volume-weighted particles discharging to the receptors. 

These estimated discharge volumes can be converted to massloading rates by multiplying these volumes by an assumed 

nitrate concentration in the treated wastewater.

The utility of flow models in evaluating wastewater-disposal scenarios is illustrated by a hypothetical scenario in eastern Barnstable. This area, which is the most urbanized part of 

Cape Cod, has an operating wastewater-disposal facility—the 

Barnstable Wastewater Disposal Facility (BWDF)—and four 

potential new wastewater-disposal sites (sites B, C, D, and F). 

The area also includes several production wells and sensitive 

ecological receptors. The hypothetical scenario includes the 

sewering of eastern Barnstable and the disposal of 1.7 Mgal/d 

of wastewater at the BWDF and 1.22 Mgal/d at sites B, C, and 

F. Although this hypothetical scenario was not analyzed as 

part of modeling support to the community, the volumes and 

locations of wastewater disposal are consistent with scenarios 

done as part of that modeling support, and the results illustrate 

the types of analyses that would be useful in evaluating actual 

scenarios.

The ground-water-flow system is characterized by an 

east-west trending regional ground-water divide; ground-water 

flow and advective transport is away from the divide toward 

downgradient discharge locations. Near the divide, vertical 

flow is enhanced and flow lines are deep in the system. Farther 

from the divide, ground-water flow is more horizontal and 

generally shallower. Traveltimes in the aquifer, referring to the 

elapsed time between recharge at the water table and discharge 

from the system, range from essentially instantaneous next 

to discharge locations to more than 100 years near regional 

ground-water divides.

According to the simulation results, wastewater would be 

transported by advection away from the current and proposed 

disposal sites toward discharge locations, and several sensitive 

receptors in this part of Cape Cod would receive some wastewater: 10 coastal water bodies, 1 stream, 7 ponds, and 11 production wells. The wastewater-disposal sites are in different 

hydrologic settings within the ground-water-flow system. The 

BWDF and sites B and F are at or near the regional groundwater divide, and site C is farther from the divide. Wastewater 

originating from source areas near the ground-water divide 

would be transported over a wider arc of flow directions and 

deeper into the aquifer than wastewater emanating from the 

source area farther from the divide. These advective-transport 

patterns indicate that wastewater disposed near the divide 

would likely contaminate a larger volume of the aquifer. Of 

the 5.36 Mgal/d of discharged wastewater, about 70 percent 

discharged into coastal waters, 12 percent into streams, and 

18 percent into production wells; about 13 percent of the 

wastewater flows through freshwater ponds prior to discharge 

from the system. The volumes of wastewater that discharged 

into individual coastal water bodies ranged from about 0.01 

to nearly 1 Mgal/d. Discharge volumes into both ponds and 

wells ranged from about 0.01 to more than 0.2 Mgal/d. About 

0.65 Mgal/d of wastewater discharged into Phlashes Brook, 

the only stream that received wastewater. Wells, ponds, 

streams, and coastal waters differ in their capacity to assimilate wastewater-derived nitrate, and therefore calculating 

wastewater-discharge volumes to individual receptors from 

each individual source area is useful in evaluating the potential 

for adverse effects of wastewater disposal on the receptors. 

Because freshwater ecosystems are not as sensitive to nitrogen 

as are marine ecosystems, ponds and streams are receptors 

of less concern than are coastal receptors. Also, open coastal 

waters are of less concern than inland estuaries because open 

coastal waters have more tidal flushing and circulation. Wells 

are receptors of concern because of the potential for adverse 

effects on human health.

Nitrate may be attenuated over long periods of subsurface transport primarily by the process of denitrification, and 

traveltimes also are important in evaluating the potential for 

adverse effects on receptors that receive wastewater-derived 

nitrate. Traveltimes in the aquifer range from essentially 

instantaneous near discharge boundaries to hundreds of years 

near regional ground-water divides. Traveltimes of wastewater discharged from the source areas differ according to 

hydrologic setting. For the hypothetical scenario, wastewater discharged from source areas near ground-water divides 

were as large as several hundred years, whereas traveltimes 
Summary 39

for wastewater discharged from the source area farther from 

regional ground-water divides—where horizontal groundwater velocities are higher—were as small as 3 years. If 

pumping is removed from the scenario and all wastewater is 

allowed to discharge to natural receptors, median traveltimes 

of wastewater from the BWDF and sites B and F, near the 

regional divide, are 66, 31, and 48 years, respectively. The 

median traveltime of wastewater from site C farther, from the 

divide, to natural receptors is about 10 years. As an example 

of the importance of traveltimes in the analysis of regional 

wastewater disposal, both Stewarts Creek and Nantucket 

Sound receive a substantial amount of wastewater (more than 

0.5 Mgal/d) from the BWDF; however, the median traveltime 

of wastewater to Stewarts Creek is about 20 years, whereas the 

median traveltime to Nantucket Sound is about 130 years. As 

a result of the longer traveltime to the Sound, the potential for 

nitrate attenuation during transport to this receptor is greater.

A subregional solute-transport model of eastern 

Barnstable was developed and used to simulate nitrate 

transport resulting from the hypothetical wastewater-disposal 

scenario; results from the solute-transport simulations were 

compared with those from flow-based particle-tracking 

analyses to evaluate the utility of using particle-tracking 

and solute-transport methods for calculating mass-loading 

rates to wells and ecological receptors. Flow-based particle 

tracking can be used to simulate advective-flow directions 

and traveltimes, and the analyses require no further model 

development beyond the ground-water-flow model; however, 

the method does not yield substantial quantitative information 

regarding the distribution of mass in the aquifer. The approach 

represents movement of a mass of nitrate with discrete massweighted particles and, thus, is based on the assumption, likely 

to be invalid in most systems, that dispersion does not occur 

and mass is conserved along flow lines. Also, the use of flowbased particle tracking to estimate mass-loading rates and 

time-varying concentrations at receptors requires a number of 

additional postmodeling analyses. In contrast, solute-transport 

methods require additional model development but explicitly 

simulate the transport of nitrate through the subsurface and 

better estimate mass-loading rates and time-varying concentrations at receptors with no postmodeling analyses. Also, solutetransport models can simulate the effects of other components 

of transport, such as dispersion, on nitrate transport. Numerical dispersion, which is a function of model discretization, 

should be considered if solute-transport models are used to 

simulate the movement of nitrate. In advective-dominated 

aquifer systems such as that on Cape Cod, models with a sufficiently small discretization can be chosen to minimize the 

potential for numerical dispersion to be greater than actual 

dispersion in the aquifer.

Mass-loading rates to 11 production wells estimated from 

both particle-tracking and solute-transport methods generally 

were in good agreement. Mass-loading rates calculated by 

particle-tracking methods generally were higher than rates 

calculated by solute-transport methods owing to the assumption that no mass is lost during transport. Results for wells 

with high mass-loading rates and short traveltimes agreed best. 

More solute mass is conserved if nitrate concentrations reach 

steady state quickly with little time for mass to be lost during 

transport. In these cases, the invalid assumption that mass is 

conserved along flow lines is less important. In some cases, 

the two methods yielded substantially different results. Well 

94 is downgradient from a pond, and the enhanced dispersion 

of nitrate by the pond as represented by the solute-transport 

model causes concentrations and, therefore, mass-loading rates 

in wastewater intercepted by the well to be lower than those 

estimated by the particle-tracking method. The mass-loading 

rate calculated by the solute-transport model exceeded the 

particle-tracking estimate for wastewater received in well 47. 

The well is east of source area B in a diverging flow field near 

the ground-water divide. It is likely that the particle-tracking 

method does not represent advective flow, which was mostly 

to the north, with sufficient resolution to simulate nitrate 

transport in the diverging flow field and that the required 

particle density may be higher than that which is practicable. 

Flow-based particle tracking also can estimate time-varying 

concentrations by a cumulative accounting of particles arriving after different traveltimes at a well. For wells with large 

mass-loading rates and short traveltimes, the two methods 

yielded similar results. At wells 94 and 47, time-varying concentrations estimated by particle tracking were substantially 

different from those simulated by solute-transport modeling; 

however, the particle-tracking method did accurately represent the elapsed time before steady-state concentrations were 

reached.

Transport in the Cape Cod aquifer is dominated by advection, which refers to the movement of a conservative solute, 

such as nitrate, with ground-water flow. Another component of 

transport, however, is dispersion, which refers to the effect of 

aquifer heterogeneity on the differential transport of a solute 

through the aquifer and the resultant spread of mass in the 

aquifer; both aquifer heterogeneity and dispersion are scale 

dependent. Numerical dispersion is an artifact of the solutetransport model and is a function of ground-water velocities 

and the spatial and temporal discretization of the model. The 

subregional model of Barnstable has an approximate numerical dispersivity of about 50 ft; because this value is consistent 

with a range of reasonable plume-scale dispersivities for 

the Cape Cod aquifer, the model may reasonably represent 

dispersive transport with no additional simulated dispersivity. 

The sensitivity of simulated concentrations to dispersivity was 

evaluated for a range of values from 20 to 500 ft. In general, 

time-varying concentrations at wells were not greatly affected 

by dispersion. Concentrations in wastewater received in wells 

71 and 74 reached a similar steady-state value for dispersivities between 20 and 500 ft owing to the close proximity of 

the wells to the site C source area and the small traveltimes of 

nitrate from the source area to the wells. For wells 47 and 94 

farther and at greater transport distances from source areas, the 

effect of dispersion was more pronounced owing to the fact 

that the cumulative effect of dispersion was greater over the 

longer transport times and distances.
40 Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer

References Cited

Barlow, P.M., 1997, Particle-tracking analysis of contributing areas of production wells in simple and complex flow 

systems, Cape Cod, Massachusetts: U.S. Geological Survey 

Water-Supply Paper 2434, 66 p.

Barlow, P.M., and Dickerman, D.C., 2000, Numerical simulation and conjuctive-management models of the Hunt-Annaquatucket-Pettaquamscutt stream-aquifer system, Rhode 

Island: U.S. Geological Survey Professional Paper 1636, 

79 p.

Colman, J.A., Masterson, J.P., Pabich, W.J., and Walter, D.A., 

2004, Effects of aqufier travel time on nitrogen transport to 

a coastal embayment: Ground Water: Oceans Issue, v. 42, 

no. 7, p. 1069–1078.

DeSimone, L.A., Walter, D.A., Eggleston, J.R., and Nimiroski, M.T., 2002, Simulation of Ground-Water Flow and 

Evaluation of Water-Management Alternatives in the Upper 

Charles River Basin, Eastern Massachusetts: U.S. Geological Survey Water-Resources Investigations Report 02–4234, 

94 p.

Farnsworth, R.K., Thompson, E.S., and Peck, E.L., 1982, 

Evaporation atlas for the contiguous 48 states: National 

Oceanic and Atmospheric Administration Technical Report 

NWS 33, 26 p.

Fetter, C. W., 1988, Applied Hydrogeology (2d ed): Columbus, Ohio, Merrill, 592 p.

Fletcher, C. A. J., 1991, Computational techniques for fluid 

dynamics, v. 1 (2d ed): Berlin, Springer-Verlag, 401 p.

Freeze, R. A., and Cherry, J. A., 1979, Groundwater: Englewood Cliffs, N.J., Prentice-Hall, 604 p.

Garabedian, S.P., LeBlanc, D.R., Gelhar, L.W., and Celia, 

M.A., 1991, Large-scale natural-gradient tracer test in 

sand and gravel, Cape Cod, Massachusetts: 2. Analysis of 

spatial moments for a nonreactive tracer: Water Resources 

Research, v. 27, no. 5, p. 911–924.

Harbaugh, A.W., and McDonald, M.G., 1996 , User’s documentation for MODFLOW-96, an update to the U.S. Geological Survey Modular Finite-Difference Ground-Water 

Flow Model: U.S. Geological Survey Open-File Report 

94–4156, 56 p.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW-2000, The U.S Geological Survey 

modular ground-water model—User guide to modularization concepts and the ground-water flow process: U.S. 

Geological Survey Open-File Report 00–92, 121 p.

Hem, J.D., 1992, The study and interpretation of the chemical 

characteristics of natural water, 3rd ed.: U.S. Geological 

Survey Water-Supply Paper 2254, 225 p.

LeBlanc, D.R., 1984, Sewage plume in a sand and gravel 

aquifer, Cape Cod, Massachusetts: U.S. Geological Survey 

Water-Supply Paper 2218, 28 p.

LeBlanc D.R., Garabedian, S.P., Hess, K.M., Gellar, L.W., 

Quadri, R.D., Stollenwerk, K.G., and Wood, W.W., 1991, 

Large-scale natural-gradient tracer test in sand and gravel, 

Cape Cod Massachusetts—l. Experimental design and 

tracer movement: Water Resources Research, v. 27, no. 5, 

p. 895–910.

LeBlanc, D.R., Guswa, J.H., Frimpter, M.H., and Londguist, 

C.J., 1986, Ground-water resources of Cape Cod, Massachusetts: U.S. Geological Survey Hydrologic Atlas 692, 

4 pls.

Luckner , L., and Schestakow, W.M., 1991, Migration Processes in the Soil and Groundwater Zone: Chelsea, Mich., 

Lewis Publishers, 487 p.

Masterson, J.P., and Barlow, P. M., 1994, Effects of simulated 

ground-water pumping and recharge on ground-water flow 

in Cape Cod, Martha’s Vineyard, and Nantucket Island 

Basins, Massachusetts: U.S. Geological Survey Water-Supply Paper 2447, 79 p., 1 pl.

Masterson, J.P., Stone, B.D., Walter, D.A., and Savoie, Jennifer, 1997a, Hydrogeologic framework of western Cape 

Cod, Massachusetts: U.S. Geological Survey Hydrologic 

Atlas HA–741.

Masterson, J.P., Walter, D.A., and Savoie, Jennifer, 1997b, Use 

of particle tracking to improve numerical model calibration 

and to analyze ground-water flow and contaminant migration, Massachusetts Military Reservation, western Cape 

Cod, Massachusetts: U.S. Geological Survey Water-Supply 

Paper 2482, 50 p.

Masterson, J.P., and Walter, D.A., 2000, Delineation of 

ground-water recharge areas, western Cape Cod, Massachusetts: U.S. Geological Survey Water-Resources Investigations Report 00–4000, 1 pl.

Masterson, J.P, Hess, K.M., Walter, D.A., and LeBlanc, D.R., 

2002, Simulated changes in the source of ground water for 

public-supply wells, ponds, streams, and coastal areas on 

western Cape Cod, Massachusetts: U.S. Geological Survey 

Water-Resources Investigation Report 02–4143, 12 p.

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground- water-flow 

model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap. A1, 586 p.
References Cited 41

Oldale, R.N., and Barlow, R.A., 1986, Geologic map of 

Cape Cod and the Islands, Massachusetts: U.S. Geological 

Survey Miscellaneous Investigations Series Map I–1763, 

scale 1:100,000.

Oldale, R.N., 1992, Cape Cod and the islands—The geologic 

story: Parnassus Imprints, East Orleans, Mass., 205 p.

Oldale, R.N., and O’Hara, C.J., 1984, Glaciotectonic origin of 

the Massachusetts coastal end moraines and a fluctuating 

late Wisconsinan ice margin: Geological Society of America 

Bulletin, v. 95, p. 61–74.

Orzol, L.L., 1997, User’s guide for MODTOOLS—Computer 

programs for translating data of MODFLOW and MODPATH into Geographic Information Systems files: U.S. 

Geological Survey Open File Report 97–240, 86 p.

Pollock, D.W., 1994, User’s guide for MODPATH/ MODPATH-PLOT, version 3-a particle tracking post-processing 

package for MODFLOW, the U.S. Geological Survey 

finite-difference ground-water flow model: U.S. Geological 

Survey Open-File Report 94–464, 234 p.

Prudic, D.E., 1989, Documentation of a computer program to 

simulate stream-aquifer relations using a modular, finite-difference, ground-water flow model: U.S. Geological Survey 

Open File Report 88–729, 113 p.

School of Marine Science and Technology, University of 

Massachusetts at Dartmouth, 2005, Linked watershedembayment model to determine critical nitrogen loading 

thresholds for Popponesset Bay, Massachusetts: Massachusetts Estuaries Program Final Report.

Spitz, Karlheinz, and Moreno, Joanna, 1996, A practical guide 

to groundwater and solute-transport modeling: New York, 

John Wiley and Sons, 480 p.

Uchupi, Elazar, Giese, G.S., and Aubrey, D.G., 1996, The late 

Quaternary construction of Cape Cod, Massachusetts—A 

reconsideration of the W.M. Davis Model: Geological Society of America Special Paper 309, 69 p.

U.S. Air Force Center for Environmental Excellence, 2000, 

Final CS-10 in-plume remedial system design ground-water 

modeling report: AFC-J23-35S18405-M17-0005, prepared 

by Jacobs Engineering Group, Inc., for AFCEE.MMR, 

Installation Restoration Program, Otis Air National Guard 

Base, Mass., unpaginated.

Valiela, I., Geist, M., McClelland, J., and Tomasky, G., 2000, 

Nitrogen loading from watersheds to estuaries: Verification 

of the Waquoit May nitrogen loading model: Biogeochemistry, v. 49, p. 277–293.

Walter, D.A., Masterson, J.P., and LeBlanc, D.R., 2002, Simulated pond-aquifer interactions under natural and stressed 

conditions near Snake Pond, Cape Cod, Massachusetts: U.S. 

Geological Survey Water-Resources Investigations Report 

99–4174, 35 p.

Walter, D.A., and Masterson, J.P., 2003, Simulation of advective flow under steady-state and transient recharge conditions, Camp Edwards, Massachusetts Military Reservation, 

Cape Cod, Massachusetts: U.S. Geological Survey WaterResources Investigation Report 03–4053, 51 p.

Walter, D.A., Masterson, J.P., and Hess, K.M., 2004, Groundwater recharge areas and traveltimes to pumped wells, 

ponds, streams, and coastal water bodies, Cape Cod, Massachusetts: U.S. Geological Survey Scientific Investigations 

Map I–2857.

Walter, D.A., and Whealan, A.T., 2004, Simulated water 

sources and effects of pumping on surface and ground 

water, Sagamore and Monomoy flow lenses, Cape Cod, 

Massachusetts: U.S. Geological Survey Scientific Investigations Report 2004–5181, 85 p.; also available online at 

http://pubs.usgs.gov/sir/2004/5181

Weiskel, P.W., and Howes, B.L., 1992, Differential transport 

of sewage-derived nitrogen and phosphorus through a 

coastal watershed: Environmental Science and Technology, 

v. 26, no. 2, p. 352–360.

Zheng, C., 1990, MT3D: A modular three-dimensional 

transport model for simulation of advection, dispersion, 

and chemical reactions of contaminants in ground-water 

systems: U.S. Environmental Protection Agency, R.S. Kerr 

Environmental Research Laboratory, Ada, Okla, 219 p.
For additional information call or write to:

Director, U.S. Geological Survey

Massachusetts-Rhode Island Water Science Center

10 Bearfoot Road, Northborough, MA 01532

(508) 490-5000

http://ma.water.usgs.gov

Document prepared by the Pembroke and West Trenton Publishing Service Centers
Walter —Use of Numerical Models to Simulate Transport of Sewage-Derived Nitrate in a Coastal Aquifer, Central and Western Cape Cod, Massachusetts —Scientific Investigations Report 2007–5259

Printed on recycled paper Printed on recycled paper