Response Curves for Phosphorus Plume 

Lengths from Reactive-Solute-Transport 

Simulations of Onland Disposal of 

Wastewater in Noncarbonate Sand and 

Gravel Aquifers

By John A. Colman

In cooperation with the 

Massachusetts Department of Environmental Protection

Scientific Investigations Report 2004-5299

U.S. Department of the Interior

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

Gale A. Norton, Secretary

U.S. Geological Survey

Charles G. Groat, Director

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

For sale by U.S. Geological Survey, Information Services

Box 25286, Denver Federal Center

Denver, CO 80225

For more information about the USGS and its products:

Telephone: 1-888-ASK-USGS

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

Any use of trade, product, or firm names in this publication 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:

Colman, J.A., 2005, Response curves for phosphorus plume lengths from reactive-solute-transport simulations 

of onland disposal of wastewater in noncarbonate sand and gravel aquifers: U.S. Geological Survey Scientific 

Investigations Report 2004-5299, 28 p.
iii

Contents

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

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

Approach to Computer Simulations of Phosphorus-Plume Length ................................. 7

Conceptual Model ................................................................................. 7

Processes Affecting Phosphorus Transport ................................................... 8

Physical Properties ........................................................................... 8

Biogeochemical Processes .................................................................. 9

Three-Dimensional Reactive-Solute-Transport Model for the Saturated Zone ............ 11

One-Dimensional Reactive-Solute-Transport Model for the Unsaturated Zone ............ 14

Combined Unsaturated- and Saturated-Zone Models .......................................... 16

Simulation Results for the Saturated- and Unsaturated-Zone Models ........................... 17

Base-Case Saturated-Zone Model and Stages of Plume Development ..................... 17

Sensitivity Analysis ................................................................................. 20

Unsaturated-Zone Model ........................................................................ 21

Response Curves of Phosphorus Plume Length as a Function of Unsaturated-Zone 

Thickness and Effluent Phosphorus Concentration .................................... 22

Model Verification ................................................................................. 22

Evaluation of Simulation Results ........................................................................ 24

Application of Simulation Results ............................................................... 25

Model Strengths, Weaknesses, and Reliability ................................................ 25

Efficient Use of Phosphorus Removal in the Subsurface ....................................... 26

Summary and Conclusions .............................................................................. 26

References ................................................................................................ 27

Figures

 1–3. Maps showing:

 1. Measured dissolved phosphorus concentrations in the aquifer 

downgradient from sewage infiltration beds at the Massachusetts 

Military Reservation, Cape Cod, Massachusetts, 1993 .............................. 4

2. Distribution of measured concentrations of dissolved phosphorus in ground 

water approximately 2 feet below the bottom of Ashumet Pond, 1999 ............ 5

3. Extent of nitrate plume emanating from the septic system of the Walden 

Pond State Reservation ............................................................... 6

4. Photograph showing sewage-infiltration beds at a new (2002) treatment plant in 

Acton, permitted for a disposal rate of 0.25 million gallons per day ..................... 7

5. Schematic showing conceptual model of wastewater infiltration through the 

unsaturated zone from an infiltration bed and transport downgradient through 

the aquifer .................................................................................... 8

6. Map showing comparative transport of phosphorus and other sewage-derived 

constituents in the ground-water contaminant plume from wastewater disposal 

at the Massachusetts Military Reservation, Cape Cod ................................. 10
iv

7. Model grid showing A, Saturated-zone PHAST model grid used in this study; 

and B, The modeled extent of the dissolved phosphorus plume at year 70 for 

base case—50 years of continuous discharge—with dissolved phosphorus 

concentrations shown on a log scale ...................................................... 12

 8, 9. Graphs showing:

8. Base-case unsaturated-zone PHAST model boundaries and grid ............... 16

9. Solid-phase phosphorus concentrations beneath the infiltration beds at the 

Massachusetts Military Reservation, Cape Cod .................................... 16

 10–16. Graphs showing:

10. Base-case simulations of the stages of phosphorus plume development in 

the saturated zone. ..................................................................... 18

11. Simulated stages of development of the base-case phosphorus plume 

length (defined as the extent of dissolved phosphorus equal or greater than 

0.015 milligram per liter) in the saturated zone. .................................... 19

12. Comparison of flux of effluent phosphorus into the infiltration bed and 

output through the model boundaries from 1,700 to 2,500 feet 

downgradient ........................................................................... 20

13. Model sensitivity analysis showing base-case phosphorus plume length 

and plume lengths altered in response to perturbation of the input 

parameters .............................................................................. 21

14. Phosphorus concentration of effluent entering and exiting the unsaturated 

zone for simulations of the Walden Pond leachfield and for the base-case 

conditions ................................................................................. 22

15. Simulated unsaturated-zone concentrations of phosphorus minerals in the 

solid phase at 10-year intervals during active discharge of wastewater, and 

of sorbed phosphorus after 10 years ................................................ 23

16. Simulated dissolved phosphorus-plume length, as a function of unsaturatedzone thickness and phosphorus concentration in the discharge .................. 24

Tables

1. Solution compositions for uncontaminated ground water and treated sewage 

effluent. ....................................................................................... 13

 2. Major chemical reactions used in reactive-transport simulations ..................... 14

3. Base-case values for the three-dimensional reactive-transport model and the 

unsatured zone model. ........................................................................ 15
v

Multiply By To obtain

foot (ft) 0.3048 meter (m)

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

foot per year (ft/yr) 0.3048 meter per year (m/yr)

gallon per day (gal/d) 3.785 liter per day (L/d)

gallon per square foot per day 

(gal/ft2/d)

40.7 liter per square meter per day 

(L/m2/d)

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

mile (mi) 1.609 kilometer (km)

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

square foot (ft2) 0.09290 square meter (m2)

square mile (mi2) 2.590 square kilometer (km2) 

ton per year (ton/yr) 0.9072 metric ton per year (ton/yr)

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

°F=(1.8×°C)+32

Vertical coordinate information is referenced to the North American Vertical Datum of 1929 

(NGVD 29).

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

kg/yr kilogram per year

mg/g milligram per gram

mg/L milligram per liter

mL milliliter

μmol/L micromoles per liter

Conversion Factors, Datums, and Abbreviations

Abstract

Surface-water resources in Massachusetts often are 

affected by eutrophication, excessive plant growth, which has 

resulted in impaired use for a majority of the freshwater ponds 

and lakes and a substantial number of river-miles in the State. 

Because supply of phosphorus usually is limiting to plant 

growth in freshwater systems, control of phosphorus input to 

surface waters is critical to solving the impairment problem. 

Wastewater is a substantial source of phosphorus for surface 

water, and removal of phosphorus before disposal may be 

necessary. Wastewater disposed onland by infiltration loses 

phosphorus from the dissolved phase during transport through 

the subsurface and may be an effective disposal method; 

quantification of the phosphorus loss can be simulated to 

determine disposal feasibility. In 2003, the U.S. Geological 

Survey, in cooperation with the Massachusetts Department 

of Environmental Protection, initiated a project to simulate 

distance of phosphorus transport in the subsurface for plausible conditions of onland wastewater disposal and subsurface 

properties. A coupled one-dimensional unsaturated-zone and 

three-dimensional saturated-zone reactive-solute-transport 

model (PHAST) was used to simulate lengths of phosphorus 

plumes. Knowledge of phosphorus plume length could facilitate estimates of setback distances for wastewater-infiltration 

sites from surface water that would be sufficient to protect 

the surface water from eutrophication caused by phosphorus 

transport through the subsurface and ultimate discharge to 

surface water. 

The reactive-solute-transport model PHAST was used 

to simulate ground-water flow, solute transport, equilibrium 

chemistry for dissolved and sorbed species, and kinetic regulation of organic carbon decomposition and phosphate mineral 

formation. The phosphorus plume length was defined for the 

simulations as the maximum extent of the contour for the 

0.015 milligram-per-liter concentration of dissolved phosphorus downgradient from the infiltration bed after disposal 

cessation. Duration of disposal before cessation was assumed 

to be 50 years into an infiltration bed of 20,000 square feet 

at the rate of 3 gallons per square foot per day. Time for the 

maximum extent of the phosphorus plume to develop is on 

the order of 100 years after disposal cessation. Simulations 

indicated that phosphorus transport beyond the extent of the 

0.015 milligram-per-liter concentration contour was never 

more than 0.18 kilogram per year, an amount that would likely 

not alter the ecology of most surface water. 

Simulations of phosphorus plume lengths were summarized in a series of response curves. Simulated plume lengths 

ranged from 200 feet for low phosphorus-concentration 

effluents (0.25 milligram per liter) and thick (50 feet) unsaturated zones to 3,400 feet for high phosphorus-concentration 

effluents (14 milligrams per liter) discharged directly into the 

aquifer (unsaturated-zone thickness of 0 feet). Plume length 

was nearly independent of unsaturated-zone thickness at phosphorus concentrations in the wastewater that were less than 

2 milligrams per liter because little or no phosphorus mineral 

formed at low phosphorus concentrations. For effluents of 

high phosphorus concentration, plume length varied from 

3,400 feet for unsaturated-zone thickness of 0 to 2,550 feet for 

unsaturated-zone thickness of 50 feet.

Response Curves for Phosphorus Plume Lengths from 

Reactive-Solute-Transport Simulations of Onland Disposal 

of Wastewater in Noncarbonate Sand and Gravel Aquifers

By John A. Colman
Model treatments of flow and equilibrium-controlled 

chemistry likely were more accurate than rates of kinetically 

controlled reactions, notably precipitation of iron-phosphate 

minerals; the kinetics of such reactions are less well known 

and thus less well defined in the model. Sensitivity analysis 

indicated that many chemical and physical aquifer properties, 

such as hydraulic gradient and model width, did not affect the 

simulated plume length appreciably, but duration of discharge, 

size of infiltration bed, amount of dispersion, and number of 

sorption sites on the aquifer sediments did affect plume length 

appreciably. 

Because simulation of plume length in carbonate-mineral 

sediments indicated that the plume would be substantially 

longer than in noncarbonate-mineral sediments, the application of the response curves in locations with carbonate-mineral 

sediments would be inappropriate. The effect of carbonate minerals in sediments is to increase pH, which causes 

decreased sorption of phosphorus on aquifer sediments.

Phosphorus removal from solution by precipitation onto 

aquifer sediments is more efficient at high concentrations of 

disposed phosphorus than at low concentrations. At very low 

phosphorus concentrations, the solubility product of phosphorus minerals is not exceeded and no phosphorus mineral 

forms. An important consequence is that removal of dissolved 

phosphorus from the plume by processes in the subsurface is 

decreased the more that removal efforts are applied in treatment before wastewater is disposed.

Model simulations indicate that removal of phosphorus 

from wastewater disposed through septic systems would have 

the advantage of efficient phosphorus removal in the subsurface because phosphorus concentrations are high in septicsystem effluent. Short plume lengths result from wastewater 

disposal through septic systems because of the efficient phosphorus removal and because of the low volume of wastewater 

involved. The simulation results for small-volume systems are 

not quantitative, however, because wastewater-infiltration rates 

are much lower than those of the higher-volume system that 

was used to calibrate the model and to create the plume-length 

response curves.

The response curves for phosphorus plume lengths, as 

defined by the maximum extent of the 0.015 milligram-perliter concentration contour, is clearly defined in the model 

simulations, although the relation between simulated plume 

length and protective setback distance is subject to interpretation. Phosphorus does move beyond the point at which the 

simulated 0.015 milligram-per-liter concentration contour has 

stopped, so that a determination of protective setback distance 

must include a consideration of whether that continued flux, 

or some other flux amount, is appropriate. Also, simulations 

indicate that phosphorus plumes do not reach their full extent 

until 50 to 200 years after disposal cessation, depending on 

concentration of phosphorus disposed. No phosphorus plume 

has been monitored for that long after cessation, so there is no 

way to verify the long-term simulation results.

Introduction

Onland disposal of wastewater has been practiced for 

many years in Massachusetts at large and small wastewatertreatment facilities underlain by porous sediment deposits. 

Publicly owned facilities treating millions of gallons per day 

use onland disposal to stabilize waste further and lessen the 

effect of discharge on surface waters as compared to disposal 

through pipe outfalls. Smaller facilities, disposing on-site 

through septic-system leachfields, rely on waste stabilization 

processes that occur in the subsurface. Although generally 

environmentally beneficial to surface water, as compared to 

outfall discharge (Jantrania, 2000), onland disposal results in 

ground-water transport of wastewater constituents, including 

phosphorus, that may affect surface freshwater upon discharge 

from the aquifer.

Transport of phosphorus in the subsurface from onland 

disposal likely depends on physical and chemical properties of the aquifer that could remove phosphorus from the 

dissolved (mobile) fraction, and on volume and phosphorus concentration of the disposed wastewater. Quantitative 

knowledge of how these factors affect phosphorus would 

facilitate licensing of new onland disposal facilities. In 2003, 

the U.S. Geological Survey (USGS), in cooperation with 

the Massachusetts Department of Environmental Protection, 

initiated a project to simulate distance of phosphorus transport 

for plausible conditions of onland disposal and subsurface 

properties. The simulation model was based on data and a 

model developed for the sewage plume at the Massachusetts 

Military Reservation (MMR) (Parkhurst and others, 2003) and 

data from the septic-system leachfield at Walden Pond State 

Reservation, Massachusetts (Colman and Friesz, 2001). The 

simulation results are likely applicable to other noncarbonate sand and gravel aquifers in addition to the glacial outwash 

deposits of Cape Cod and eastern Massachusetts on which the 

model was based.

2 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
The multiple concerns regarding wastewater disposal 

include removal of solids, biological oxygen demand (BOD), 

bacteria, nitrogen, phosphorus, trace metals, and emerging 

contaminants (for example, pharmaceuticals and endocrine 

disrupters). Depending on degree of treatment, many of these 

materials are decreased before disposal at treatment plants. 

Onland disposal of wastewater through correctly functioning septic systems can also be effective in eliminating the 

polluting constituents of wastewater, because these materials degrade in the unsaturated zone or are immobilized by 

adsorption and precipitation in the unsaturated zone (the 

zone of the subsurface that is above the water table) and the 

saturated zone (the aquifer) (Robertson, 2003; Weiskel and 

others, 1996; Kent and others, 2000, Kolpin and others, 2003). 

Although less investigated than other constituents, the emerging contaminants may also be degraded effectively or removed 

from ground water onto the solid phase with onland disposal, 

because many of these compounds are biodegradable organic 

compounds and may also be surface reactive. Phosphorus, a 

cause of eutrophication in freshwater, and nitrogen, a potential 

contaminant of drinking-water supplies and cause of eutrophication of coastal marine ecosystems, are only partially 

removed in treatment plants and from the dissolved phase 

during transport in the subsurface. Subsurface transport of 

phosphorus, in addition to nitrogen, where drinking-water 

supplies or coastal marine resources are involved, is a concern 

for locating onland disposal systems.

Phosphorus is usually the limiting nutrient for plant 

growth in freshwater; therefore, supply to surface water can 

control eutrophication (Scheffer, 1998). Of the more than 

three-quarters of the lake and pond acreage surveyed in 

Massachusetts for an annual report to the U.S. Environmental 

Protection Agency and Congress described as “not supporting 

all beneficial uses,” the most common reason for the unsupported use was excessive plant growth (U.S. Environmental 

Protection Agency, 1998). Introduction of nonnative species 

may have an effect; but overall, plant growth can continue 

only to the extent that plant nutrients, such as phosphorus, 

are available. Disposal of wastewater can result in substantial 

eutrophication of receiving waters because of the phosphorus 

contained in the effluents (Moore and others, 2003). Onland 

disposal of wastewater may partially or wholly mitigate the 

effects of phosphorus on surface water by phosphorus removal 

processes during transport through the unsaturated zone and in 

the saturated zone. 

Quantifying the absorptive capacity of the unsaturated 

and saturated zones for phosphorus would allow this capacity 

to be used most effectively. The capacity of the subsurface to 

absorb phosphorus clearly depends on the phosphorus load 

applied, and there are examples of onland disposal systems 

that were overloaded resulting in substantial amounts of 

phosphorus transported to surface water (figs. 1 and 2), as well 

as examples where no appreciable transport of phosphorus 

has occurred (fig. 3). In addition to load, physical, chemical, 

and biological processes affect phosphorus retention in the 

subsurface.

Time-dependent effects of these processes affect absorptive capacity and must be considered because plumes of 

dissolved phosphorus may persist in the subsurface long after 

active discharge has ceased. Because persistence of phosphorus plumes in the aquifer is long in comparison to the lengths 

of investigations that have measured their effects, long-term 

effects must be determined by computer model simulations 

based on the known geochemistry and hydrology of the subsurface. Until recently, models have been inadequate to simulate the interacting three-dimensional flow (3-D) and chemical 

and biological reactions that determine phosphorus transport. 

However, PHAST, a USGS-developed 3-D reactive-solutetransport model, has become available (Parkhurst and others, 

2004). By using PHAST simulations, engineering-response 

curves can be prepared that show how the length of the phosphorus plume in the aquifer relates to chemical, biological, 

and aquifer properties including phosphorus concentration in 

the effluent and thickness of the unsaturated zone.

This report describes the hydrology and geochemistry that affects phosphorus retention in the unsaturated and 

saturated zones of the subsurface. Computer models are 

described that simulate transport of phosphorus for onland 

disposal of wastewater primarily for large installations, with 

20,000 ft2 infiltration beds, infiltration rates of 3 gal/ft2/d, 

and disposal duration of 50 yrs. Simulation results are summarized in graphical engineering response curves that show 

the relations among effluent phosphorus concentrations, 

thickness of the unsaturated zone, and length of the dissolved 

phosphorus plume in the aquifer. The dissolved phosphorus 

plume (hereafter called the phosphorus plume) is the zone 

in the aquifer characterized by above-background dissolved 

phosphorus concentrations. Use of the response curves can 

facilitate determination of the disposal location and the degree 

of phosphorus-removal treatment needed to meet low-impact 

disposal objectives. Sensitivity analysis is included to quantify 

the effects on the response curves of the number of phosphorus binding sites on the aquifer solids, duration of disposal, 

and other aquifer and disposal conditions.

Introduction 3
Ashumet

Pond

Base from U.S. Geological Survey digital data

1:24,000, 1991

Universal Transverse Mercator projection

Zone 19

Modified from Walter and others (1996)

0

0

500

250

1,000 2,000 FEET

500 METERS

70 ̊32'50"

41 ̊38'30"

41 ̊38'

SANDWICH ROAD

CURRIER ROAD

TOWN OF FALMOUTH

TOWN OF MASHPEE

TOWN OF FALMOUTH

TOWN OF SANDWICH

0.1

1.0

2.0

3.0

4.0

5.0

10.0

3.0

4.0

5.0

6.0

1.0

0.1

0.1 LINE OF EQUAL MAXIMUM DISSOLVEDPHOSPHORUS CONCENTRATION—Dashed 

where approximately located. Interval, in 

milligrams per liter, is variable.

SAMPLING SITE—Sites include observation 

 wells, multilevel samplers, and screened 

 auger borings

EXPLANATION

ACTIVE SEWAGE-DISPOSAL BED

ABANDONED SEWAGE-DISPOSAL BED

70 ̊32'20"

CAPE COD 

BAY

NANTUCKET

SOUND

MMR

Ashumet Pond

Study Site

70o 70 00' o 71 30' o

00'

42o

00'

41o

30' 0 10

0 10 20 KILOMETERS

20 MILES

Figure 1. Measured dissolved phosphorus concentrations in the aquifer downgradient from sewage infiltration beds at the 

Massachusetts Military Reservation (MMR), Cape Cod, Massachusetts, 1993. Disposal ranged from 0.1 to 1.5 million gallons per 

day during 60 years of use.

4 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
Figure 2. Distribution of measured concentrations of dissolved phosphorus in ground water approximately 2 feet below 

the bottom of Ashumet Pond, 1999. The phosphorus derives from the sewage plume emanating from sewage infiltration 

beds at the Massachusetts Military Reservation, Cape Cod, Massachusetts.

EXPLANATION

U.S. Geological Survey digital data, Universal Transverse Mercator projection, 

Zone 19, 1:24,000, 1991; State plane coordinate system datum is NAD 83 in meters.

0.1

0.1

0.1

1.0

0.5

2.0

0.1 0.1

0.1

 3.0

0.5

Sandwich Road

Currier Road

0

0

250 500 FEET

50 100 METERS

Ashumet

Pond

Boat Ramp

TOWN OF MASHPEE

TOWN 

OF FALMOUTH

F624

1.9

F388

5.5

F565

1.6

F619

4.6

F622

2.4

F621

3.1

F300

1.6

F564

0.83

F424

1.2

F239

0.15

F239

0.15

Modified from McCobb and others (2003)

AREAL EXTENT OF DISSOLVED PHOSPHORUS 

 IN POND-BOTTOM GROUND WATER (1999)—In 

 milligrams per liter as P. Contour is dashed where inferred. 

 Interval is varied

MULTILEVEL-SAMPLER SITE, IDENTIFIER, AND 

 MAXIMUM PHOSPHORUS CONCENTRATION IN 

 GROUND WATER—In milligrams per liter as P

TEMPORARY POND-BOTTOM 

 DRIVE-POINT SAMPLER Greater than 3.0

2.0 to less than 3.0

1.0 to less than 2.0

0.5 to less than 1.0

0.1 to less than 0.5 

Less than 0.1

41o38'05"

70o 70 32'20" o32'25"

41o38'00"

280,000 280,100 280,200

821,000

820,900

820,800

Introduction 5
160

170

180

190

200

210

220

210

190

170

190

170

200

ESTIMATED EXTENT OF PLUME

BUILDINGS AND ROADS

WATER-TABLE CONTOUR—Contour interval 

is half-foot, dashed where inferred

TOPOGRAPHIC CONTOUR—Contours

are in 10-foot increments, NGVD 1929

MONITORING WELL—Alphanumeric code

is local identifier

BEACH PERIMETER TEMPORARY

DRIVE-POINT SAMPLER—Number is specific 

conductance, in microsiemens per centimeter at 25oCelsius

NITRATE CONCENTRATION—Number indicates

nitrate concentration, in milligrams nitrogen per liter;

bar height is proportional to nitrate concentration; <, less than

From Colman and Friesz, 2001

42o30'

73o00' 72o

00' 71o

00'

41o30'

0 50 KILOMETERS

0 50 MILES

300,000

200,000 100,000 950,000

850,000

MASSACHUSETTS

Walden Pond

213,400 213,600

910,000

909,800

Figure 3. Extent of nitrate plume emanating from the septic system of the Walden Pond State Reservation, Massachusetts. 

Dissolved phosphorus concentrations were at background values and are not shown. Essentially no phosphorus from the 

septic system disposal site has been transported to the pond during the 25 years of system operation. Disposal rate averages 0.001 million gallons per day.

6 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
Approach to Computer Simulations of 

Phosphorus-Plume Length

The approach to simulations uses a conceptual model to 

define conditions of wastewater disposal, phosphorus plume 

formation, and safe setback distances for locations of disposal 

from surface water. The conceptual model is quantified by use 

of computer models based on chemical and physical principles 

of subsurface flow and transport.

Conceptual Model

Typically, onland disposal of wastewater at treatment 

plants involves application of treated sewage effluent to 

infiltration beds excavated into the ground. A recent installation of infiltration beds for onland disposal of treated sewage 

has been completed in Acton, MA. The beds were set next to 

a 0.25 Mgal/d treatment facility, initially running 0.1 Mgal/d; 

effluent will be treated to decrease phosphorus concentrations to 0.15–0.2 mg/L; monitoring wells were installed 50 

to 100 ft downgradient from the beds (Acton Department of 

Health, oral commun., 2003; fig. 4). Alternatively, onland 

disposal of sewage can be through septic-system leachfields 

without pretreatment other than settling in septic tanks with 

effluent phosphorus concentrations of 10–15 mg/L. Common 

for the small scale of individual homes with effluent volume 

of 300–400 gal/d, septic systems can also be used to dispose 

of sewage at larger facilities such as the Walden Pond State 

Reservation, Concord, MA (up to 5,500 gal/d; Colman and 

Friesz, 2001; fig. 3).

Wastewater disposed to the ground infiltrates vertically 

through the unsaturated zone to the water table where it flows 

as ground water through the aquifer and eventually discharges 

to surface water (fig. 5). The required safe setback of the infiltration beds from downgradient surface-water bodies depends 

on whether, or at what point downgradient, the continued 

movement of dissolved phosphorus is effectively limited by 

interaction with the solid phase so that transport of phosphorus 

to a surface-water body is not a problem. Phosphorus transport 

in the unsaturated and the saturated zones is affected directly 

by processes of phosphorus precipitation and adsorption to 

the solid phase and indirectly by geochemistry that alters the 

rates or amounts of precipitation and adsorption. Simulation 

of phosphorus plume length to facilitate determination of 

protective onland disposal setback distance as a function of 

disposal and aquifer conditions is the objective of the modeling effort and can be quantified once the computer model has 

been developed. 

Figure 4. Sewage-infiltration beds at a new (2002) treatment plant in Acton, Massachusetts, permitted for a disposal rate 

of 0.25 million gallons per day.

Approach to Computer Simulations of Phosphorus-Plume Length 7
Processes Affecting Phosphorus Transport

Phosphorus transport in the subsurface depends on physical and biogeochemical processes. The complex interactions 

of these processes that affect phosphorus transport are best 

accounted for by computer simulations.

Physical Properties

Ground-water velocity, which refers to speed and direction of ground-water flow, determines where a plume will 

migrate from a given wastewater-disposal location. The effect 

of physical properties on flow direction can be simulated with 

computer models (MODFLOW) in a useful way for watersupply development. Many MODFLOW models for aquifers in Massachusetts are available (for example, Masterson 

and Barlow, 1997). These models incorporate the effects of 

recharge, withdrawals, porosity, hydraulic conductivity and 

simulate hydraulic head, ground-water velocities, and water 

flux. By using available flow-model results, determination of 

flow direction may be made, for example, in the vicinity of a 

proposed infiltration bed and a nearby kettle-hole pond. In the 

absence of a ground-water-flow model, determinations of flow 

direction may be derived from ground-water-level information, 

such as available water-table maps. 

Physical properties of the aquifer also affect chemistry. 

Flow rate and dispersion affect the contact time of water with 

aquifer solids and the volume of aquifer solids encountered by 

a ground-water phosphorus plume, and thus, the progress of 

kinetically controlled reactions and the size of the reaction surface area for reactions between dissolved phosphorus and the 

solid phase. The reaction surface area also depends on grain 

size—for a given weight of sediment, the smaller the grain 

size, the greater will be the surface area. The effect of physical properties on both ground-water flow and chemistry can 

be evaluated by simulation with PHAST, the computer model 

used in this study.

The range of possible values for physical properties 

investigated in this study is limited by practical considerations for wastewater disposal: at installations treating tens of 

thousands of gallons per day, a conductive aquifer is required 

so that infiltration rates will be in the range of 2–3 gal/ft2/d, 

and the area required for infiltration beds can be kept reasonably small. Aquifer sediment, such as a sand or fine sand, will 

accept 2–3 gal/ft2/d; these materials are present throughout 

Cape Cod and southeastern Massachusetts as glacial outwash deposits. Valley-fill deposits in central and western 

Massachusetts are also often relatively coarse-grained and will 

accept infiltration at a high rate. Commonly, infiltration is sufficiently high in the same kinds of geological deposits as are 

suitable for public ground-water-supply development; maps 

are available for these deposits (for example, Frimpter, 1980; 

Friesz, 1996).

Figure 5. Conceptual model of wastewater infiltration through the unsaturated zone from an infiltration bed and transport downgradient through the aquifer.

8 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater

Schematic only, not to scale Modified from Walter and others (2002)

Plume

Aquifer

Direction of ground-water flow

Water

table

Pond surface

Pond-bottom sediments

Pond

Wastewater infiltration bed

Unsaturated zone

Unsaturated zone
Biogeochemical Processes

Knowledge of geochemical controls and field observations indicate that phosphorus infiltrating the ground as a 

consequence of sewage disposal does not biodegrade, as 

organic contaminant compounds might, or change into benign 

forms, such as nitrate conversion into unreactive nitrogen 

gas. Rather, inactivation of phosphorus in the subsurface with 

respect to transport and discharge to surface water can occur 

only through phosphorus removal from the dissolved phase 

by adsorption on the aquifer solids or precipitation in mineral 

phases.

Phosphorus adsorption and precipitation can account for 

substantial amounts of phosphorus removal from the dissolved 

phase of wastewater plumes; these processes greatly retard 

the transport of phosphorus with respect to ground-water 

velocity. Retardation of phosphorus transport in the subsurface is clearly demonstrated at the decommissioned sewagedisposal site at the MMR, Cape Cod, MA, by comparison of 

the sewage-plume extent, determined by elevated conductance, 

and the phosphorus-plume extent, determined by elevated 

phosphorus concentration, downgradient from the infiltration 

beds (fig. 6). During 60 years of sewage disposal at the site, 

from 1936 to 1995, the estimated 516 metric tons of phosphorus disposed (Parkhurst and others, 2003) has been transported 

about 2,500 ft downgradient in dissolved concentrations 

greater than 0.1 mg/L, whereas more mobile constituents of 

the contaminant plume have progressed 15,000 ft, possibly to 

discharge at the wetland head of a coastal embayment (fig. 6) 

(Walter and others, 1996). Although the bulk of the sewage 

plume moved toward the coast, a portion of the phosphorus 

plume discharged to the kettle hole Ashumet Pond, only 

1,600 ft downgradient from the beds (fig. 2), and contributed 

to the eutrophication of that water body (McCobb and others, 

2003).

In another location, Walden Pond, Concord, MA, phosphorus attenuation associated with onland disposal of sewage 

has been so great as to render phosphorus concentrations in a 

sewage plume indistinguishable from background concentrations. At Walden Pond, where the unsaturated zone is 40 ft 

thick, untreated sewage from a visitor center and bathhouse 

has been discharged to a septic-system leachfield since 1975 

at an average rate of 0.001 Mgal/d. Phosphorus concentrations 

above background levels have not been detected in wells 

screened at the water table as close as 10 ft downgradient 

from the leachfield (Colman and Friesz, 2001). Attenuation 

of dissolved phosphorus plumes also has been measured in 

noncarbonate-mineral aquifers studied in Canada (Robertson, 

2003).

Phosphorus is present as an amphoteric oxyanion in water 

(dihydrogen or monohydrogen phosphate ion) that, depending 

on pH, may bind strongly with metal hydroxides of aluminum 

and iron at aquifer-solid surface sites. Unconsolidated deposits 

and deposit surface coatings in much of Massachusetts are 

formed from crystalline bedrock of metasediments and granite. 

As the deposits weather, ferrous iron is released, oxidized, and 

precipitated as iron oxyhydroxide coating on aquifer solids. 

This surface, and similar aluminum oxyhydroxide surfaces, 

interact electrostatically, because of the charged solids and 

charged anions, and through covalent bonding to adsorb 

anions like phosphate.

The descriptive sorption chemistry of oxyanions on pure 

metal hydroxide surfaces is well known (Dzombak and Morel, 

1990). Because aquifer-solid surfaces contain hydroxides of 

various metals, their sorption characteristics may be similar, 

but not the same, as measured for the pure systems. Several 

adsorption-desorption experiments have been completed on 

overburden sediments from eastern Massachusetts: column 

experiments involving sorption and desorption with uncontaminated sediments, and column and batch experiments using 

wastewater-contaminated sediments (Walter and others, 1996; 

Parkhurst and others, 2003). Estimates of the logarithm of the 

equilibrium constant (log K) for the reaction

 SiteOH + PO4

3- + 2H+

 = SiteHPO4

-

 + H2O, (1)

where Site represents a phosphorus-sorption site on the aquifer 

solids, are 25.39 for purely iron hydroxide sites (Dzombak 

and Morel, 1990), and 26.7 for sites present on aquifer solids 

(Parkhurst and others, 2003). The number of sites can be 

computed by fitting experimental data from column experiments with the double-layer model of Dzombak and Morel 

(Parkhurst and others 2003). The estimated number of phosphorus sorption sites measured for the aquifer sediments present at the MMR is 3,450 μmol/L (Parkhurst and others, 2003). 

Because hydrogen ion binds at both cation and anion 

sorption sites, phosphorus sorption can be affected by changes 

in pH and in cation concentration. These effects were observed 

during column experiments when a spike in phosphorus 

desorption occurred coincident with a change in the eluent 

from an ionic composition similar to treated sewage to that 

of background aquifer water (Walter and others, 1996). The 

results indicate that phosphorus release may occur from sorption sites in the aquifer when uncontaminated ground water 

moves into the plume after cessation of wastewater discharge. 

This effect is important to include in models that simulate 

phosphorus transport in the aquifer.

Less is known about the reactions of chemical precipitation that include phosphorus—the rates of formation and the 

types of materials that form—than is the case for sorption 

of phosphorus. Mineral-forming constituents of the ironphosphate minerals vivianite and strengite (Parkhurst and 

others, 2003) and the aluminum-phosphate mineral variscite 

(Robertson, 2003) typically are supersaturated in wastewater 

plumes, so formation is under kinetic rather than solubility 

control.

Approach to Computer Simulations of Phosphorus-Plume Length 9
EXPLANATION

EXTENT OF PHOSPHORUS PLUME FROM SEWAGE- 

 TREATMENT SITE NEAR ASHUMET POND

 (from McCobb, 1999)

EXTENT OF SEWAGE PLUME

MODEL-SIMULATED WATER-TABLE 

 CONTOUR—Shows altitude of simulated water table. 

 Contour interval is 10 feet. Datum is NGVD 29

 (from Savoie, 1995)

BUZZARDS 

BAY

VINEYARD SOUND

NANTUCKET SOUND

5 MILES

5 KILOMETERS

0

0

1

2

Ashumet

Pond

Base from U.S. Geological Survey digital data

Universal Transverse Mercator projection

Zone 19, 1:24,000, 1991

70o37' 30" 70o

27' 30"

2 3 4

1 3 4

Johns

Pond

Snake

Pond

SANDWICH

MASHPEE

FALMOUTH

BOURNE Peters

Pond

Mashpee

River

Santuit

Pond

Santuit River

Coonamessett

Pond

Wakeby Pond

Mashpee

Pond

60

50

40

30

20

10

10

Round

Pond

Long

Pond

Crooked

Pond

Jenkins

Pond

ATLANTIC

OCEAN

Green

Pond

Flax

Pond

Mystic

Lake

Middle

Lake

Hamblin

Lake

Moody

Pond

Coonamessett R.

Backus R.

Bourne R.

Child R.

Quashnet R.

Waquoit

Bay

MASSACHUSETTS 

MILITARY 

RESERVATION

41 ̊42'30"

41 ̊32'30"

815,000

270,000 280,000

42o30'

73o00' 72o

00' 71o

00'

41o30'

0 10 20 30 40 50 KILOMETERS

0 10 20 30 40 50 MILES

300,000

200,000 100,000 950,000

850,000

MASSACHUSETTS

820,850

Figure 6. Comparative transport of phosphorus and other sewage-derived constituents in the ground-water contaminant plume from wastewater disposal at the Massachusetts Military Reservation, Cape Cod, Massachusetts.

10 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
Sediment investigation has demonstrated that substantial 

retention of phosphorus occurs by processes other than sorption. During experimental work at the MMR, sediments from 

the contaminated zone were washed repeatedly with pH 10 

water. This washing released the sorbed fraction of sediment 

phosphorus. Total nondesorbable phosphorus that remained 

after the washing ranged from about 0.015 to 0.25 mg/g of 

sediment (Parkhurst and others, 2003), whereas background 

phosphorus concentration is consistently about 0.02 mg/g of 

sediment (Walter and others, 1996). Apparently, the nondesorbable phosphorus was deposited in mineral form and had 

accumulated to concentrations as much as 10 times greater 

than background concentrations. Evidence from another site, 

Walden Pond, indicates that substantial amounts of phosphate 

minerals form in the unsaturated zone as well. At Walden 

Pond, essentially all of the phosphorus discharged through the 

septic-system leachfields during 25 years has been sequestered 

in the unsaturated zone (Colman and Friesz, 2001), whereas 

mass-balance calculations using the number of sorption sites 

available indicate that breakthrough of phosphorus to the saturated zone would be expected if sorption alone were involved. 

Other evidence for mineral formation has been obtained by 

microprobe analysis of phosphorus-aluminum-iron compounds 

in authigenic coatings on sand grains collected by coring 

beneath septic-system disposal tiles (Robertson, 2003).

Although phosphorus in natural waters has only one 

important oxidation state, other constituents that precipitate with phosphorus, such as iron, change oxidation state 

according to redox conditions in the aquifer. Because other 

redox-sensitive species, such as nitrogen, manganese, carbon, 

and oxygen, will determine the redox environment of iron, 

the biogeochemistry of all of these constituents needs to be 

included in a comprehensive phosphorus-transport simulation. 

Accounting for all the factors that could affect phosphorus 

transport, therefore, requires a model that simulates many 

chemical constituents and reactions simultaneously. 

Carbonate aquifers, such as those in the Housatonic 

and Hoosic River Valleys, may retard phosphorus less than 

the noncarbonate aquifers, which are present in the rest of 

Massachusetts. The higher pH in ground water from aquifers 

that are present in carbonate regions would cause phosphorus to sorb less readily than in noncarbonate aquifers. Also, 

comparison of septic-system plumes in carbonate and noncarbonate aquifers indicates that phosphorus is transported in 

the former but not the latter (Robertson, 2003). Phosphorus 

retardation in noncarbonate aquifers is attributed by Robertson 

(2003) to aluminum dissolution from gibbsite, under acidic 

conditions that develop in noncarbonate-aquifer plumes, and 

reprecipitation with phosphate. In contrast, less aluminum 

would dissolve under the near-neutral conditions that prevail 

in carbonate aquifers. 

Three-Dimensional Reactive-Solute-TransportModel for the Saturated Zone

The computer model of the saturated zone was based 

on previous computer models of phosphorus transport at 

the MMR on Cape Cod, MA. These models include a onedimensional (1-D) model (Stollenwerk, 1996; Stollenwerk 

and Parkhurst, 1999) and a 3-D model (Parkhurst and others, 

2003). The previous models were created to simulate transport of phosphorus in the aquifer after 60 years of wastewater 

disposal in infiltration beds at the MMR. Investigation of 

the aquifer downgradient from the infiltration beds has been 

conducted for 25 years with a network of hundreds of well 

clusters and multilevel samplers (Walter and others, 1996; 

LeBlanc, 1984). Data from the network were used to verify 

the models. The more recently developed 3-D model used new 

estimates of organic-matter decay rates, phosphorus sorption, 

and phosphorus precipitation and included the effect of 3-D 

ground-water flow in shaping the plume.

Despite the variable wastewater loadings associated with 

large fluctuations in population as military efforts changed at 

the MMR, the 3-D model accurately simulated distribution 

of 1993 and 1999 dissolved-phosphorus concentrations in the 

aquifer downgradient from the infiltration beds (Parkhurst 

and others, 2003). Because of the mechanistic geochemical 

and flow features built into the model for simulation of future 

phosphorus transport resulting from past variable disposal 

rates of wastewater, the model can be adapted for predicting 

phosphorus plume length as a function of loading rate in other 

similar situations.

The computer model used in this investigation to simulate 3-D transport in the saturated zone is a finite-difference 

code developed by the USGS called PHAST (Parkhurst and 

others, 2004). PHAST is a combination of the geochemical 

code PHREEQC (Parkhurst and Appelo, 1999) and the finitedifference flow-and-transport code HST3D (Kipp, 1997). 

PHAST is a general-purpose code that can accurately simulate 

flow, transport, and chemical processes controlled by chemical 

equilibrium or by reaction kinetics. 

The principal modifications to the previous 3-D 

saturated-zone PHAST model of Parkhurst and others 

(2003) involved simplifying the ground-water-flow-system 

geometry and decreasing the number of infiltration beds to 

one. Geometry of the base-case model grid includes 100-ft 

node spacing in the horizontal directions, extending 5,300 ft 

in the direction of flow by 1,100 ft in the lateral direction, and 

10-ft node spacing in the vertical direction, extending 150 ft 

in depth (fig. 7). At each time step (0.5 yr), the model solves 

for hydraulic head and constituent concentrations at each point 

of intersection in the grid (the model nodes), and for groundwater and chemical-constituent flux across each cell face. Cell 

faces extend from halfway between adjacent nodes to halfway 

to the next node or to the node if the node is on the boundary.

Approach to Computer Simulations of Phosphorus-Plume Length 11
Figure 7. A, Saturated-zone PHAST model grid used in this study. Boundaries are constant-flux on top, 

leaky boundaries on ends as indicated by yellow shading, and no-flow on other sides; and B, The modeled 

extent of the dissolved phosphorus plume at year 70 for base case—50 years of continuous discharge—with 

dissolved phosphorus concentrations shown on a log scale.

LEAKY BOUNDARY

BOUNDING BOX

FOR MODEL GRID

BOUNDING BOX

FOR SATURATED

OR PARTLY SATURATED

MODEL CELLS IN

CUTAWAY VIEW

DISTANCE DOWNGRADIENT, IN FEET

0

2,500

5,000

DISTANCE DOWNGRADIENT, IN FEET

0

2,500

5,000

1,100 FEET

150 FEET

A.

B.

 Area of infiltration bed

 Area of infiltration bed

EXPLANATION

DISSOLVED PHOSPHORUS CONCENTRATION,

 IN MILLIGRAMS PER LITER

10

1

0.1

0.01

0.001

12 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
Hydraulic boundaries of the model are no-flow boundaries on the long sides and the bottom, leaky boundaries on 

the short ends, and flux boundary on the top (fig. 7). The 

flux boundary allows recharge from rainfall, 0.0059 ft/d (26 

in/yr), except in the area of the infiltration bed, where recharge 

from wastewater flux was 0.40 ft/d (3 gal/ft2/d). The sewageinfiltration bed, 200 by 100 ft, was set transverse to the direction of flow with its center 200 ft from the upgradient leaky 

boundary and midway between the lateral no-flow boundaries 

(fig. 7). The leaky boundaries, which allow flow in or out 

of the model domain, are constrained by remote fixed heads 

outside of the model grid—50 ft above an arbitrary datum 

upgradient and 30 ft above the datum downgradient, 7,000 ft 

distant from, and parallel to, the leaky boundaries. The average 

hydraulic gradient established by the fixed heads is 0.001 ft/ft. 

Water generally flows in through the upgradient leaky boundary and out through the downgradient leaky boundary as 

determined by the fixed heads of the leaky boundaries and 

recharge through the top of the model. Hydraulic properties of 

the saturated-zone model include the hydraulic conductivity of 

300 ft/d in the horizontal directions and 85 ft/d in the vertical direction, and a porosity of 0.39. These model parameters 

result in an average ground-water velocity of about 1.5 ft/d. 

The hydraulic properties were taken from the previous MMR 

model (Parkhurst and others, 2003) based on measurements 

made on Cape Cod.

Although the model can simulate dispersion, it is effective only when dispersivity, the dispersion coefficient, is 

greater than the inherent numerical dispersivity of the model. 

Numerical dispersivity is a function of model-cell size and 

time step approximately according to the relation

D = ∆X/2 + V∆T/2, (2) 

where D is the numerical dispersivity, ∆X is the cell dimension 

in the X direction, ∆T is the time step and V is the velocity 

(Fletcher, 1991). In the saturated-zone model, dispersivity in 

the principal direction of flow is approximately

100/2 + 1.5(180)/2 = 185 ft. (3)

This dispersivity is much greater than the longitudinal dispersivity of 3 ft estimated in tracer-test experiments in the aquifer 

at the MMR (Garabedian and others, 1991). In the simulations, 

dispersivity was set to the small value of 0.0003 ft, so that 

the effective longitudinal dispersivity would be the numerical 

dispersivity alone.

Initial solution conditions for the saturated-zone model 

are the solute concentrations measured in uncontaminated 

ground water at the MMR (table 1). Initial conditions of 

the reactive surfaces were set to be in equilibrium with the 

uncontaminated (background) ground water. The composition 

of water moving into the model through the upgradient leaky 

boundary was set to uncontaminated ground water (table 1) 

as was the composition of water recharged as rainwater. 

Rainwater was assumed to be altered to the composition of 

uncontaminated ground water during passage through the soil 

and unsaturated zones. Chemical composition of wastewater 

effluent recharged through the infiltration bed was the same as 

that measured in wastewater at the MMR (table 1).

The chemical reactions simulated in the saturated-zone 

model were the same as those in the MMR sewage-plume 

investigation (Parkhurst and others, 2003). These included 

reactions relating directly to dissolved phosphorus concentration, such as sorption on metal oxyhydroxide surfaces and 

precipitation, and those relating indirectly through redox, pH 

or ionic-strength effects. Ionic-strength effects were included 

in all equilibrium reactions. 

Equilibrium was assumed for homogeneous aqueous 

reactions including protonation, ion-pair formation, and redox 

reactions as found in the database to PHREEQC (Parkhurst 

and Appelo, 1999). In addition, equilibrium and kinetic reactions involving sorption, mineral equilibria, organic carbon, 

and dissolved nitrogen gas were included (table 2). 

Table 1. Solution compositions for uncontaminated ground 

water and treated sewage effluent.

[Modified from Parkhurst and others, 2003. Concentrations are milligrams 

per liter, unless otherwise noted. TDIC, total dissolved inorganic carbon; 

°C, degrees Celsius]

Constituent or physical 

property

Uncontaminated 

ground water

Secondary 

treated 

sewage 

effluent

pH, standard units 5.6 6.0

Temperature, °C 14 14

Oxygen 4.0 4.0

Calcium 1.2 13

Magnesium .75 4.1

Sodium 4.6 48

Potassium .39 9.3

Sulfate, as SO4 8.3 27.9

Chloride 15.6 148

Nitrate, as N 0 15

Ammonium, as N 0 2.5

Phosphate, as P 0 5.9

Iron 0 .18

Manganese .035 .022

Total organic carbon, as C 0 19

Total dissolved inorganic 

carbon, as C

.34 14

1

Chloride was adjusted to achieve charge balance in the solution.

Approach to Computer Simulations of Phosphorus-Plume Length 13
Sites for the sorption reactions were of two types—

for cations (Cation_site, table 2) or for phosphate (Site, 

table 2)—except for hydrogen ion, which was allowed to react 

at both types of sites. Cations allowed to sorb were ammonium, calcium, magnesium, sodium, iron, and manganese. A 

single equilibrium constant was applied to all cation-sorption 

reactions because insufficient data were available to distinguish individual equilibrium constants for each cation (tables 2 

and 3). Phosphate and hydrogen ion were the only chemical 

species allowed to sorb at the phosphate sites. The surface 

complexation mass-action equation used to simulate sorption 

includes an electrostatic term to account for the work done in 

moving a charged species to a charged surface (Parkhurst and 

Appelo, 1999, p. 14–15).

Kinetics were assumed to control formation and dissolution of strengite and vivianite, decay of organic carbon, 

and removal of nitrogen gas. The rate or equilibrium-constant 

values used in the reactions, method of determination, and 

reference for calibration are summarized in table 3. The 

kinetic formulation for precipitation of minerals involved the 

saturation index so that the reaction rate would approach zero 

as equilibrium was established. Denitrification of nitrate was 

under thermodynamic control, and kinetic removal of nitrogen 

gas (reaction 13, table 2) was included to prevent reoxidation to nitrate. Sorption of organic carbon was under kinetic 

control during the active discharge period. After cessation of 

wastewater discharge, decay of the sorbed carbon was also 

expressed by kinetic reaction (reaction 10, table 2). 

One-Dimensional Reactive-Solute-Transport 

Model for the Unsaturated Zone

Flow in the unsaturated-zone model was simulated as 

1-D transport. The 1-D model is a simplification of the complex processes of unsaturated flow. Nevertheless, the time of 

transport, ratio of solution to reactive surfaces, and exposure to 

oxygen of actual unsaturated flow can be approximated with 

a 1-D model so that kinetic, equilibrium, and redox reactions 

are simulated as accurately as possible. For the base case, 15 

column cells were used to cover a vertical distance of 50 ft 

(fig. 8), so that each cell was 3.33 ft thick. Simulations of the 

unsaturated zone with different thicknesses were accomplished 

by adding or subtracting column cells to produce the required 

unsaturated-zone thickness. Vertical flow was simulated with 

the TRANSPORT option of PHREEQC, included in PHAST, 

14 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater

Table 2. Major chemical reactions used in reactive-transport simulations.

[Modified from Parkhurst and others (2003); C, concentration; k, rate constant; SI, saturation index—log ion activity product divided by the solubility 

product; DOC, dissolved organic carbon; indicates absolute value]

Reaction

number

Reaction Reaction

type Equation

1 Phosphorus sorption Equilibrium SiteOH PO4

3- 2H+ SiteHPO4

- + + = + H2O

2 Phosphorus-sorption site protonation Equilibrium SiteOH H+ SiteOH2

+ + =

3 Phosphorus-sorption site deprotonation Equilibrium SiteOH SiteO- H+ = +

4 Cation sorption Equilibrium Cation_siteOH Catx+ Cation_siteOCatx – 1 H+ + = +

5 Cation-sorption site protonation Equilibrium Cation_siteOH H+ Cation_siteOH2

+ + =

6 Cation-sorption site deprotonation Equilibrium Cation_siteOH Cation_siteO- H+ = +

7 Fe(OH)3(amorphous) Equilibrium Fe OH ( )3 3H+ Fe+3+ = + 3H2O

8 Pyrolusite Equilibrium MnO2 4H+ 2e

- Mn+2

+ + = + 2H2O

9 Decomposition of dissolved organic carbon 

(DOC), reactant is CH2O Kinetic Rate kCDOC = –

10

Decomposition of sorbed organic carbon, 

reactant is CH2O (applies after cessation of 

wastewater disposal)

Kinetic Rate kO2CO2, aq – kNO3CNO3, aq = –

11 Vivianite precipitation and

dissolution, reactant is Fe3(PO4)2

Kinetic Rate kCPO4

SIvivianite

1 SI + vivianite

----------------------------------     = –

12 Strengite precipitation and

dissolution, reactant is FePO4

Kinetic Rate kCPO4

SIstrengite

1 SI + strengite

----------------------------------     = –

13 Removal of dissolved N2 (prevents N2

reoxidation to nitrate), reactant is N2

Kinetic Rate kCN2, aq = –
wherein water is shifted from one cell to the next down the 

column of cells. The shifts occur according to a specified timestep, and concentrations of all constituents are computed after 

each shift according to the chemical reactions specified in the 

model. At the base-case infiltration rate, 3 gal/ft2/d, downward 

velocity in the infiltration pond would be 0.401 ft/d, equivalent 

to 1.02 ft/d vertical velocity in the unsaturated zone, based 

on an aquifer porosity of 0.39. With 3.33-ft column cells, the 

number of cell volumes, or shifts, required per year to achieve 

the given infiltration rate is 114. 

The unsaturated-zone model generally is based on the 

same chemistry as that developed for the saturated-zone 

model. In the unsaturated-zone model, as in the saturatedzone model, phosphorus is retained by the solid phase through 

sorption and mineral precipitation. One difference between 

models, however, is that oxygen is introduced in the unsaturated-zone model to account for oxygen in soil gas to which 

the water would be exposed during infiltration. Thirty-percent 

saturation with oxygen was specified for a zone beginning 

6.7 ft below the surface and extending down to the water table. 

Thus, oxygen was available for reactions in the unsaturated 

zone and in the saturated zone because effluent entering the 

aquifer at the water table also contained dissolved oxygen.

Whereas the saturated-zone calibration can be based 

on that used for the MMR, unsaturated-zone transport needs 

to be calibrated in this study because no investigation of the 

unsaturated zone was included in the previous investigation. 

Calibration data are available from total amounts of solidphase phosphate measured in the sediment cores from the 

unsaturated zone under the sewage beds at the MMR (Walter 

and others, 1996). 

Background concentration of phosphorus in the solid 

phase was 0.02 mg/g (Walter and others, 1996), whereas solidphase phosphorus concentrations below an infiltration bed at 

the MMR ranged from 0.1 to 0.4 mg/g (fig. 9), with the highest concentrations immediately below the bed. Concentrations 

in the rest of the unsaturated zone and the top 20 ft of the 

saturated zone were about 0.1 mg/g. The maximum amount 

of sorbed phosphorus can be calculated from the number of 

anion sorption sites (Parkhurst and others 2003). Given anion 

sorption sites of 3,450 μmol/L of pore water, occupied entirely 

with phosphorus, the maximum concentration of sorbed 

phosphorus is 0.025 mg/g of aquifer solid, assuming pore 

volume of 0.39 and rock density of 2.7 (Parkhurst and others, 

2003). Thus, the sum of background and sorbed phosphorus 

concentrations is about 0.05 mg/g, which leaves 0.05 to 0.35 

mg/g that is attributable to mineral forms built up during 60 

yr of discharge. Such mineral formations are described by 

Robertson (2003) in his investigation of septic plumes with 

a conceptual model that iron is reduced near the beds and 

infiltrates to lower levels in the unsaturated zone, where it is 

oxidized and redeposited as Fe(OH)3 or FePO4. For the unsaturated-zone model to produce the appropriate concentrations 

of solid-phase phosphorus, the rate constant for formation and 

dissolution of strengite was increased from the saturated-zone 

value of 2.3 x 10-8 to 1.6 x 10-7 s-1 (table 3).

Table 3. Base-case values for the three-dimensional reactive-transport model and the unsatured zone model.

[Modified from Parkhurst and others (2003); P, phosphorus; μmol/L, micromole/liter; s-1, per second]

Property or process Saturatedzone model

Unsaturatedzone model Source Reference

P sorption constant (log K, eq. 1) 26.7 26.7 Column sorption/desorption 

experiments 

Walter and others, 1996; Parkhurst 

and others, 2003

P sorption sites, μmol/L 3,450 3,450 Column sorption/desorption 

experiments

Walter and others, 1996; Parkhurst 

and others, 2003

Cation-sorption constant (log K, 

reaction 4, table 2)

-1.8 -1.8 Column sorption/desorption 

experiments

Parkhurst and others, 2003

Cation-sorption sites, μmol/L 18,000 18,000 Column sorption/desorption 

experiments

Parkhurst and others, 2003

First-order rate constant for 

dissolved organic carbon, s-1

1.0 x 10-7 1.0 x 10-7 Fit to field data for O2 and NO3 Parkhurst and others, 2003

Rate constant for strengite, s-1 2.3 x 10-8 1.6 x 10-7 Fit to data for accumulation of 

nondesorbable phosphate

Parkhurst and others, 2003

Rate constant for vivianite, s-1 2.3 x 10-8 1.6 x 10-7 Fit to data for accumulation of 

nondesorbable phosphate

Parkhurst and others, 2003

y

Approach to Computer Simulations of Phosphorus-Plume Length 15
Combined Unsaturated- and Saturated-Zone 

Models

Because the reactive constituents—phosphorus, organic 

carbon, and nitrogen species—are altered as a function of time 

in the unsaturated zone as well as of the duration of disposal, 

the inputs to the saturated-zone model require a direct coupling of the two models. The output of the unsaturated-zone 

model becomes the time varying input to the saturated-zone 

model. The last concentration emerging from the unsaturated 

zone after each year was used as input concentration to the 

saturated-zone model for that year. The combined model was 

used to produce the plume-length response curves that are the 

product of this project. 

Combining the models could potentially alter the calibration relation between source loads and aquifer transport of 

phosphorus of the model constructed for this project, relative to the MMR model, because the MMR model did not 

include the unsaturated zone. At the phosphorus calibration 

concentration of the MMR model, however, the effect on phosphorus concentrations of the unsaturated-zone model through 

sorption is small because the sorption sites rapidly reach a 

steady state relative to the effluent phosphorus concentration, 

after which the input phosphorus concentration is not affected 

INFLOW FROM

INFILTRATION BED15 CELLS, 50 FEET

OUTFLOW TO

SATURATED-ZONE MODEL

Figure 8. Base-case unsaturatedzone PHAST model boundaries and 

grid.

0

70

10

20

30

40

50

60

0 0.2 0.4 0.5

Water table

(Unsaturated zone)

(Saturated zone)

EXTRACTABLE PHOSPHORUS 

(IN HYDROCHLORIC ACID), IN 

MILLIGRAMS PER GRAM

ALTITUDE, IN FEET ABOVE NGVD 29

Figure 9. Solid-phase phosphorus concentrations beneath the 

infiltration beds at the Massachusetts Military Reservation, Cape 

Cod, Massachusetts (Walter and others, 1996). Vertical bars 

represent cored interval. Horizontal bars represent variability 

among replicate extractions. Dashed line represents background 

extractable phosphorus concentrations collected from extractions at an upgradient sampling site.

16 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
by sorption. With respect to formation of phosphorus minerals, 

the unsaturated-zone chemistry results in continuous removal 

of about 20 percent of the phosphorus from the input under 

base-case conditions, such that the phosphorus concentration 

entering the saturated zone is constant at about 80 percent of 

the effluent phosphorus concentration.

Simulation Results for the Saturatedand Unsaturated-Zone Models

The saturated-zone model can be used to define the 

approach to estimating phosphorus plume length by indicating stages of plume development and phosphorus transport 

during and after disposal. The unsaturated-zone model output 

modifies the input to the saturated zone, but not the qualitative 

aspects of plume transport stages.

Base-Case Saturated-Zone Model and 

Stages of Plume Development 

The base-case saturated-zone model has one infiltration 

bed (100 by 200 ft) with a constant loading rate of 3 gal/ft2/d 

for 50 yr of wastewater with a composition of that of the 

secondary treated sewage at the MMR (table 1). This model 

simulates how phosphorus plumes might be distributed during 

discharge, immediately after discharge cessation, and long 

after discharge cessation. During a 200-yr simulation, progression of the leading edge of the phosphorus plume during 

active discharge (0–50 yr) is slow, about 12 ft/yr (fig. 10). 

After discharge cessation (at year 50), when the aquifer is 

beginning to receive clean water, plume progression almost 

triples to 34 ft/yr, and the highest concentration zone leaves 

the area under the infiltration beds and moves toward the leading edge of the plume (fig. 10). The counter-intuitive increase 

in plume speed following cessation of discharge results 

from desorbed phosphorus responding to changes in pH and 

ionic strength; the same effects were apparent in the column 

experiments that were used to calibrate the model (Parkhurst 

and others, 2003). The desorbed phosphorus increases the 

dissolved phosphorus concentration, which moves downgradient with ground-water flow to the front of the plumes where 

unoccupied sites are present. After the fast-progression period 

(year 50 to year 70), the plume slows again, and by year 100, 

forward movement has stopped, as best seen by plume crosssection plots of dissolved phosphorus concentrations in the 0 

to 0.015 mg/L range (fig. 11). Change after year 100 is a slow 

decrease in dissolved phosphorus concentration throughout the 

model cross section (fig. 11).

Even after downgradient progression of dissolved 

phosphorus in the 0 to 0.015 mg/L cross-section plot appears 

to stop, dissolved phosphorus transport continues, but at low 

rates. The continued phosphorus transport associated with low 

concentrations can be determined at an arbitrarily chosen contour of low concentration, for example, 0.015 mg/L, approximately the dissolved concentration of phosphorus measured at 

background locations in the aquifer of Walden Pond (Colman 

and Friesz, 2001). For the base-case simulation, the longest 

extent of the 0.015 mg/L phosphorus contour is 2,500 ft 

(indicated by the leading edge of the red zone, fig. 11), and 

transport of phosphorus through a plane at that location varies 

between 0 and 0.17 kg/yr during a period of 0 to 200 years 

after initiation of effluent discharge (fig. 12). For reference, 

the annual phosphorus input load from all sources for an oligotrophic lake, Walden Pond, was estimated as 32 kg/yr (Colman 

and Friesz, 2001). The phosphorus load to Walden Pond is 

greater by a factor of 190 than the highest rate of transport 

through the downgradient edge of the base-case plume, and 

greater by a factor of 670 than the average downgradient 

transport rate in the first 200 years. The simulations indicate 

that the association of phosphorus with the aquifer solids is 

substantial and that phosphorus is released to the dissolved 

phase slowly. If ground-water movement through an aquifer is 

considered as a nutrient-removal process for wastewater, the 

calculated phosphorus-removal efficiency for the base-case 

model over the first 200 years is greater than 99.9 percent.

Phosphorus fluxes were also evaluated at distances 

closer than the 0.015 mg/L contour and increased sharply 

with decreasing distance to the infiltration beds (fig. 12). An 

8-percent decrease in distance (from 2,500 ft to 2,300 ft) corresponded to a 380-percent increase in the maximum phosphorus flux. Although a greater phosphorus flux might be deemed 

protective for surface water than that which continues past the 

0.015 mg/L contour, fluxes quickly become large (and likely 

not protective) as distances downgradient are decreased from 

that of the 0.015 mg/L contour.

In summary, the model simulates the stages of plume 

evolution during and after cessation of wastewater disposal. 

Although phosphorus migrates quickly immediately after 

cessation, the dissolved phosphorus plume eventually is characterized by low concentrations of dissolved phosphorus at the 

leading edge. The position of the contour for 0.015 mg/L stops 

advancing, and phosphorus flux at the leading edge is small by 

comparison with that entering oligotrophic lakes. Therefore, 

the distance from the beds to a selected low-concentration 

contour may help to indicate a protective setback distance 

from receiving surface water for an onland wastewaterdisposal location. The simulation effort of this project will 

focus on predicting the phosphorus plume length, defined as 

the distance from infiltration beds to the 0.015 mg/L phosphorus contour, for a range of wastewater loading and aquifer 

properties. 

Simulation Results for the Saturated- and Unsaturated-Zone Models 17
0

50

100

150

0

50

100

150

0

50

100

150

0

50

100

150

0

50

100

150

DEPTH BELOW WATER TABLE, IN FEET

0

50

100

150

0 1,000 2,000 3,000 4,000 5,000 0 1,000 2,000 3,000 4,000 5,000

DISTANCE DOWNGRADIENT FROM INFILTRATION BEDS, IN FEET

DISSOLVED PHOSPHORUS CONCENTRATION, 

IN MILLIGRAMS PER LITER—Log scale

3 2 1 0.1 0.03

10 YEARS

20 YEARS

30 YEARS

40 YEARS

50 YEARS

60 YEARS 120 YEARS

110 YEARS

100 YEARS

90 YEARS

80 YEARS

70 YEARS

Figure 10. Base-case simulations of the stages of phosphorus plume development in the saturated zone.

18 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
Figure 11. Simulated stages of development of the base-case phosphorus plume length (defined as the 

extent of dissolved phosphorus equal or greater than 0.015 milligram per liter) in the saturated zone.

Simulation Results for the Saturated- and Unsaturated-Zone Models 19

greater than or equal to 0.01500

0.01125

0.00750

0.00375

0.0

0

50

100

150

0

50

100

150

0

50

100

150

0

50

100

150

DEPTH BELOW WATER TABLE, IN FEET

DISTANCE DOWNGRADIENT FROM INFILTRATION BEDS, IN FEET

0 1,000 2,000 3,000 4,000 5,000 0 1,000 2,000 3,000 4,000 5,000

DISSOLVED PHOSPHORUS CONCENTRATION, 

IN MILLIGRAMS PER LITER—Log scale

50 YEARS

60 YEARS

70 YEARS

80 YEARS

90 YEARS

100 YEARS

120 YEARS

110 YEARS
Sensitivity Analysis

A sensitivity analysis was conducted to determine the 

relative dependence of phosphorus transport on aquifer properties and other model parameters. For sensitivity testing, only 

the base-case saturated-zone model was used. The unsaturated 

zone was omitted, in part, to save computer run time and 

also because the most sensitive variable in the unsaturatedzone model is unsaturated-zone thickness, which is one of 

the parameters varied in runs of the combined model. To test 

model sensitivities, parameters were usually increased or 

decreased by a factor of 2. Parameters increased by a factor 

of 2 were aquifer width (model width), duration of discharge, width of the infiltration bed, and hydraulic gradient. 

Parameters decreased by a factor of 2 were aquifer thickness, 

number of phosphorus sorption sites per liter of pore water, 

and effluent discharge rate. Other factors tested were elimination of iron-phosphate mineral precipitation, equilibration with 

dolomite (which buffers the pH at about 7.5), conversion of 

the model to two dimensions (2-D), and use of a finer mesh 

grid and shorter time step. The latter decreased dispersivity, in 

the longitudinal direction, for example, from 185 to 93 ft.

The sensitivity analysis of the base-case saturated-zone 

model indicated that the largest changes in plume length 

were in response to equilibration with dolomite, a carbonate 

rock (plume length reached the model length, 5,100 ft, at 260 

years and was still increasing). Equilibration with dolomite 

increased pH, which substantially decreased the sorption of 

phosphorus on the metal-hydroxide surfaces of the aquifer. 

This result was in agreement with observations in the field of 

extended transport of phosphorus plumes from small septic 

systems in carbonate as compared to noncarbonate aquifers 

(Robertson, 2003).

None of the factor-2 parameter changes caused a proportional change in plume length (fig. 13), although response 

to increased discharge duration was nearly so (plume length 

increased by 92 percent), and the response for decreased effluent discharge was also substantial (plume length decreased 

by 40 percent). The less-than-proportional response between 

parameter change and the size of the response was caused by 

chemical precipitation kinetics. When the dissolved phosphorus concentration is increased, such as by increased discharge 

or decreased sorption sites, the precipitation kinetics (reaction 

12, table 2) are faster. When the result of the parameter change 

is lower dissolved phosphorus concentrations, precipitation 

occurs at a slower rate. Thus, the effects on plume length of 

both increased and decreased dissolved phosphorus concentrations are compensated.

Increasing the width of the model decreased the plume 

length only 2 percent, indicating that the no-flow lateral 

boundaries of the model were not causing elongation of the 

plume. Changing the model to 2-D (in effect, to one cell 

width), however, caused the plume to lengthen by 48 percent. Increasing hydraulic gradient by a factor of 2 increased 

0

100

200

300

400

500

600

700

0 100 200 300 400

TIME, IN YEARS

PHOSPHORUS INPUT, IN KILOGRAMS PER YEAR 

0

2

4

6

8

10

12

14

PHOSPHORUS FLUX DOWN

GRADIENT,

IN KILOGRAMS PER YEA

R

Infiltration bed input

Flux 2,500 feet downgradient

Flux 2,300 feet downgradient

Flux 2,100 feet downgradient

Flux 1,900 feet downgradient

Flux 1,700 feet downgradient

Figure 12. Comparison of flux of effluent phosphorus into the infiltration bed and output through the model 

boundaries from 1,700 to 2,500 feet downgradient. Input and output are plotted at different scales.

20 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
plume length by only 13 percent. Decreasing the number 

of sorption sites by half caused a 33-percent increase in the 

plume length. When mineral precipitation of phosphorus was 

eliminated, plume length increased 33 percent; lack of greater 

response indicates that sorption, rather than chemical precipitation, is the dominant phosphorus-removal mechanism. 

Finally, decreasing the numerical dispersion by decreasing the 

grid mesh and time step caused the plume length to increase 

48 percent. Decreasing the dispersion caused the plume 

to become more narrow and to be exposed to less aquifer 

sediment. 

Unsaturated-Zone Model

Some of the constituents of wastewater change concentration during simulated percolation through the unsaturated 

zone, and the degree of change varies over the period of simulation. The concentration of dissolved phosphorus eluting 

from the unsaturated zone decreased from the base-case input 

concentration of 5.9 mg/L to nearly 0 during the first 1.5 

years. After 1.5 years, phosphorus breakthrough occurred. 

After breakthrough, concentrations eluting from the unsaturated zone increased to 4 mg/L but never reached the basecase input concentration (fig. 14). The difference between 

dissolved phosphorus concentrations after breakthrough and in 

the input is caused by precipitated phosphorus, which builds 

up continuously in the unsaturated zone (fig. 15). All sorption sites in the unsaturated zone are filled by the beginning 

of dissolved phosphorus breakthrough at 1.5 years (fig. 14). 

As the sorption sites fill, dissolved phosphorus concentrations increase, and precipitation of iron-phosphate minerals 

begins, continuing thereafter for the duration of the wastewater 

disposal. In the unsaturated-zone model, iron phosphate is 

formed by a kinetically controlled reaction (reaction number 

12, table 2) whose rate is dependent on the phosphorus concentration. Thus, the minerals form slowly, or not at all, while 

soluble phosphorus is filling sorption sites on the solid surfaces but at a greater rate after the sites are filled and dissolved 

phosphorus concentrations increase.

Generally, the unsaturated-zone model was run for just 

the 50 yr of active disposal. A model simulation continued 

past active disposal, however, indicated that, as in the saturated case, dissolved phosphorus concentrations are higher 

when disposal stops and background water (rainwater) 

moves through the unsaturated zone. In the unsaturated zone, 

however, the flux of water decreases after disposal cessation by a factor of 88, unlike in the aquifer, where water flux 

is controlled substantially by regional gradients in groundwater heads. In the unsaturated zone, water flux is limited to 

recharge only. The decreased flow of water in the unsaturated 

zone when disposal of wastewater ceases causes transport of 

all constituents, including phosphorus, to decrease greatly.

PARAMETER VARIED

Base case

2 x Base case

PLUME LENGTH, IN FEET

5,000

4,000

3,000

2,000

1,000

6,000

0

Base case

2 x width

2 x duration

2 x infiltration bed

2 x hydraulic gradient

1/2 aquifer thickness

1/2 sorption sites

1/2 discharge

No phosphate minerals

Dolomite equilibration

Two-dimensional model

1/2 discretization

Figure 13. Model sensitivity analysis showing base-case phosphorus plume length and plume lengths altered in 

response to perturbation of the input parameters. The arrow over the response for dolomite equilibration indicates that 

the plume was still increasing when the model boundary at 5,100 feet was reached.

Simulation Results for the Saturated- and Unsaturated-Zone Models 21
Response Curves of Phosphorus Plume 

Length as a Function of UnsaturatedZone Thickness and Effluent Phosphorus 

Concentration

Plume lengths, as represented by the 0.015 mg/L dissolved phosphorus contour, were simulated as a function of 

unsaturated-zone thickness and effluent phosphorus concentration using the combined model (fig. 16). Plume length 

ranged from 200 ft for the greatest unsaturated-zone thickness 

(50 ft) and lowest phosphorus concentration (0.25 mg/L), to 

3,400 ft at the lowest unsaturated-zone thickness (0 ft) and 

highest phosphorus concentration (14 mg/L). The relation was 

nonlinear with respect to both phosphorus concentration and 

unsaturated-zone thickness. At low phosphorus concentrations 

(less than 2 mg/L), plume length is essentially independent 

of unsaturated-zone thickness, whereas at higher concentrations, plume lengths are shorter for greater unsaturated-zone 

thicknesses. For a constant unsaturated-zone thickness, plume 

length increases less rapidly as phosphorus concentrations 

increase. The cause of the nonlinear response to unsaturatedzone thickness is the variable rate of precipitation of strengite 

(iron phosphate) in the unsaturated zone; the rates depend on 

the phosphorus concentration. Large amounts of strengite form 

in the unsaturated zone when the phosphorus concentration is 

high, but precipitation slows or even ceases when the phosphorus concentration is low and the solubility product for ironphosphate formation is not exceeded.

Model Verification

Verification of the combined model was difficult because 

no disposal system has been monitored long enough to provide 

accurate time histories of effluent phosphorus concentration 

and effluent flux rate, required inputs for the model. Likewise, 

phosphorus plumes from onland disposal have not been monitored for the length of time (greater than 50 years after cessation) necessary for maximum length of plumes to be developed. Comparison of simulation results with available data at 

the MMR and a small system at Walden Pond, however, is still 

of interest to determine if the model is generally in agreement 

with available data. 

At the MMR, a range of effluent disposal rates was 

used during the discharge period from 1936 to 1995, from 

0.32 to 3.1 gal/ft2/d. The number of infiltration beds and 

total discharge also fluctuated resulting in disposal areas that 

varied from 107,600 to 538,000 ft2 and total discharge from 

0.1 to 1.5 Mgal/d. The discharge in the base-case model, by 

comparison, is 3.0 gal/ft2/d at 0.06 Mgal/d, so that the basecase discharge was over a smaller infiltration area. Effluent 

0

2

4

6

8

10

12

14

0 10 20 30 40 50

TIME, IN YEARS

PH

OSPH

ORUS CONCENTRATION,

IN MILLIGRAMS PER LITE

R

Base-case breakthrough

Walden Pond leachfield input

Base-case input

Walden Pond breakthrough Walden 

Pond

leachfield

Base case

Figure 14. Phosphorus concentration of effluent entering and exiting the unsaturated zone for simulations of the Walden Pond, Massachusetts, leachfield and 

for the base-case conditions. Differences between input and exiting concentrations are caused by continuous phosphate mineral formation in the unsaturated 

zone. Delay of breakthrough is caused by adsorption. 

22 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
phosphorus concentration at the MMR was estimated to be 

11.8 mg/L before a treatment-plant upgrade in 1983 and 

5.9 mg/L after the upgrade until discharge cessation in 1996. 

Plume-length results from the combined model (fig. 16) for 

12 mg/L phosphorus, which applied for the first 50 years of 

discharge at the MMR, and 20 ft of unsaturated zone present 

at the MMR is 2,500 ft. The measured plume length at the 

MMR in 1993, plotted as the 0.1 mg/L contour, (fig. 1, Walter 

and others, 1996) also was 2,500 ft. Measured plume length as 

represented by the 0.015 mg/L contour would likely be longer 

than 2,500 ft and might have increased since 1993 because 

only 10 years had passed since cessation of wastewater discharge supplying the long axis of the plume. The longer MMR 

plume results are compatible with the simulation base-case 

results because wastewater was applied to a greater area in the 

measured case than in the base case. Simulations run for the 

sensitivity analysis indicated that the use of expanded infiltration bed size with the same flux resulted in longer plumes.

Data from the septic-system leachfield at Walden Pond 

provide information on transport through the unsaturated zone 

rather than information about plume length. After 24 years of 

operation of the leachfield, phosphorus had not arrived in the 

0

10

20

30

40

50

DEPTH OF UNSATURATED ZONE, IN FEET

0 0.1 0.2

SOLID PHASE PHOSPHORUS CONCENTRATION,

IN MILLIGRAMS PER GRAM 

FePO4 – 10 years

FePO4 – 20 years

FePO4 – 30 years

FePO4 – 40 years

Sorbed P – 10–50 years

FePO4 – 50 years

EXPLANATION 

Figure 15. Simulated unsaturated-zone concentrations of phosphorus minerals in the solid phase 

at 10-year intervals during active discharge of wastewater, and of sorbed phosphorus after 10 

years. The simulated sorbed concentration remains constant during the 10-to-50-year interval.

Simulation Results for the Saturated- and Unsaturated-Zone Models 23
saturated zone (Colman and Friesz, 2001). By comparison, the 

simulation for the Walden Pond leachfield with discharge of 

14 mg/L phosphorus concentration to the 40-ft thick unsaturated zone indicated that breakthrough would occur at 25 

years (fig. 14). Thus, the measured and simulated results are 

compatible, although the results imply that phosphorus breakthrough to the saturated zone might occur in the near future. 

Note that the plume lengths represented in figure 16 are not 

appropriate for the septic system at Walden Pond because the 

disposal rate at Walden Pond in terms of gallons per day is 60 

times less than that used in figure 16.

Evaluation of Simulation Results

The relation among unsaturated-zone thickness, effluent phosphorus concentration, and phosphorus plume length 

(fig. 16) can be used to help plan setbacks between onland 

disposal sites and surface-water bodies that are long enough 

to protect the surface water from phosphorus eutrophication. 

The simulated plume lengths can be used for planning and 

cost-benefit analysis of alternative disposal designs, such as 

treatment removal of phosphorus as compared to acquisition 

of land with characteristics favorable to phosphorus removal in 

the subsurface.

Application of Simulation Results

The model simulations can be used appropriately to 

simulate plume length where parameters other than phosphorus load and unsaturated-zone thickness (explicitly tested 

in this investigation) are similar to those used in the model, 

or differ only with respect to parameters determined in the 

sensitivity analysis not to affect plume length (fig. 13). The 

sensitivity analysis indicated that the tested parameters, other 

than those associated with load and with carbonate rocks 

(pH), changed the plume length less than 33 percent. Thus, 

although the calibration data used were for the MMR, the 

model results are expected to apply more generally to sand and 

gravel aquifers, with the exception of regions with carbonEFFLUENT PHOSPHORUS CONCENTRATION, IN MILLIGRAMS PER LITER

UNSATURATED-ZONE THICKNESS, IN FEET

500

1,000

1,500

2,000

2,500

3,000

Figure 16. Simulated dissolved phosphorus-plume length, as a function of unsaturated-zone thickness and phosphorus concentration in the discharge. Plume length (illustrated with the contours; contour interval is 250 feet) is defined 

by the extent of the 0.015 milligram per liter dissolved phosphorus contour after 50 years of discharge at the rate of 

3 gallons per square foot per year onto a 20,000 square foot infiltration bed. Dots indicate the values of unsaturated-zone 

thicknesses and phosphorus concentrations of effluent used in the model runs.

24 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
ate rocks. The simulations were made with a thick aquifer 

(150 ft), although plume length increased only 21 percent 

when the aquifer thickness was reduced by half (75 ft). The 

model was run for porous, hydraulically conductive material 

that is likely to allow infiltration rates of 3 gal/ft2/d. The model 

would not be applied for hydraulically unconductive material 

because onland disposal would not be selected at sites where 

infiltration would be much less than 3 gal/ft2/d.

Whether the application of simulation results (fig. 16) is 

for a specific location or for comparison of alternative sites, 

knowledge of the direction of flow, downgradient distance 

to surface water, and unsaturated-zone thickness is required. 

Ground-water-flow direction can be determined from results 

of previous regional-flow simulations or from hydraulic-head 

measurements for at least three wells screened at the water 

table (McCobb and others, 1999). Unsaturated-zone thickness 

can be determined as the difference between land-surface 

elevation and the water table by using methods established by 

Frimpter (1981). Once unsaturated-zone thickness and downgradient distance to surface water have been established, the 

effluent phosphorus concentration corresponding to a plume 

length equal to the downgradient distance at that site can be 

obtained from figure 16. The phosphorus concentrations and 

plume lengths shown in figure 16 represent model simulations; 

their use as protective combinations for onland disposal upgradient of surface water may require additional interpretation by 

the regulating agency.

Monitoring downgradient from onland disposal installations allows the predictions to be checked and knowledge 

of phosphorus mobility to be improved for further refinements of model-simulated estimates of phosphorus plume 

length. Monitoring wells installed downgradient of onland 

disposal installations that have intercepted plumes of effluentassociated constituents (conductance, nitrate) would be appropriate to monitor for eventual breakthrough of phosphorus. 

Effluent constituents in addition to phosphorus can be 

important at an onland disposal facility. Nitrogen transport in 

ground water is a concern, for marine areas, where nitrogen 

supply may control growth of aquatic plants, and for local 

drinking-water supply. Increasing the setback distance so that 

the aquifer is used efficiently in phosphorus removal may 

expose a larger zone of ground water to high nitrogen concentrations, rendering it unfit for water supply. In near-shore 

areas, a substantial fraction of nitrogen disposed onland may 

be transported to the coast. Some nitrogen removal in the aquifer or in the sediments at the point of discharge to the marine 

environment is likely, however. As a result, onland disposal of 

nitrogen-containing waste may be preferable to direct discharge of an equivalent waste to surface water. A satisfactory 

onland disposal balancing requirements of phosphorus and 

nitrogen pertains for the septic-system leachfield at Walden 

Pond. Setback is sufficient to eliminate phosphorus transport 

to the pond. Nitrogen is transported to the pond, but apparently 

without appreciable effect on water quality because phosphorus, not nitrogen, controls the growth of aquatic plants. No 

water-supply wells are present near in the plume between the 

leachfield and Walden Pond (Colman and Friesz, 2001).

Model Strengths, Weaknesses, and Reliability

The strengths (accuracy in simulating plume length) of 

the reactive solute-transport model derive from the coupling 

of full chemistry with flow in the 1-D unsaturated and in the 

3-D saturated zones—a model structure that is state-of-theart. The model accounts for the complex dilution that results 

as infiltrated water interacts with ground water in the aquifer, 

and for the deepening of the plume downgradient because 

of recharge (fig. 10). Chemical reactions are treated with the 

full set of equilibrium and kinetic reactions available from 

PHREEQC (Parkhurst and Appelo, 1999). The sequestration of phosphorus in both unsaturated and saturated zones, 

rather than just the latter can be important. Chemical effects 

relating to control of sorption by pH and control of pH by 

ion exchange and redox chemistry are included explicitly in 

the model simulations. These explicit reactions represent the 

important chemical effects that have been observed in column 

studies, where more phosphorus is released because of pH and 

ionic-strength effects when ground water enters the plume area 

after cessation of discharge.

Model weaknesses (uncertainties in simulating plume 

length) are related primarily to difficulty in simulating chemical reactions controlled by kinetics rather than by equilibrium 

conditions and in specification of sorption properties of the 

solids. Rates of reaction for organic-carbon degradation and 

phosphate-mineral formation are important, but only a few 

calibration points are available, from laboratory data from the 

MMR for carbon oxidation and from the observed buildup of 

phosphate minerals in the unsaturated and saturated zones at 

the MMR for phosphate minerals. The formation of ironphosphate minerals is not certain, but a similar result would 

pertain for other possible precipitates, such as aluminum phosphate. Dissolution rates for the phosphate minerals that likely 

occur after disposal cessation are also uncertain. Because rates 

in kinetically controlled processes may depend in a complex 

manner on pH, bacteria populations, temperature, and concentrations of reactive species, the limited calibration data 

available from the MMR studies are probably not sufficient to 

provide robust simulations. Nevertheless, the reactions have 

been confirmed to proceed and the relations used in the kinetics have been shown to be adequate in simulating phosphorus 

transport at the MMR.

The amount of iron-phosphate mineral precipitated is 

a larger percentage of the total phosphorus removed when 

onland wastewater discharge is low and phosphorus concentration in the effluent is high, such as is the case for septic 

systems (not covered in figure 16). Thus, the predictions for 

a septic system with a small plume length are more highly 

dependent on mineral formation than are predictions for the 

Evaluation of Simulation Results 25
base-case model, which apply to high disposal rates. Whereas 

phosphate-mineral formation has been observed in analysis of 

the unsaturated zone below septic systems (Robertson, 2003), 

the rate of formation has not been quantified.

Efficient Use of Phosphorus Removal in the 

Subsurface

The curvatures of the response curves in figure 16 

provide an important conclusion regarding the efficacy of 

treatment of phosphorus before disposal onto infiltration beds; 

phosphorus removal in the sewage-treatment plant before 

discharge onto the beds results in little or no removal by precipitation of mineral phases in the subsurface because of the 

low phosphorus concentration in the effluent. The reactions 

in the subsurface, precipitation of phosphorus with iron, are 

actually among those applied for phosphorus removal in treatment plants. If phosphorus removal is carried out in the plant, 

however, the resulting effluent is too dilute for precipitation 

reactions to occur in the subsurface.

Summary and Conclusions

Phosphorus loading to surface freshwater can be the 

principal environmental problem associated with disposal of 

domestic wastewater. Onland disposal of wastewater results 

in less phosphorus transport to surface water than does direct 

discharge through an outfall pipe because of geochemicalloss processes in the subsurface. Quantification of the loss 

processes, which may occur in the saturated and unsaturated 

zones, could be used to plan onland disposal systems. In 

2003, the U.S. Geological Survey, in cooperation with the 

Massachusetts Department of Environmental Protection, initiated a project to simulate distance of phosphorus transport 

for plausible conditions of onland disposal and subsurface 

properties.

The report summarizes the hydrology and geochemistry 

that controls phosphorus retention in the unsaturated and saturated zones of the subsurface. Computer models are described 

that simulate transport of phosphorus for onland disposal of 

wastewater primarily for large installations—3 gal/ft2/d to 

a 100 by 200 ft infiltration bed. Although disposal at lower 

rates, such as rates pertinent to septic systems, are discussed, 

plume lengths are not simulated because the septic-system 

disposal rates are too small by comparison with rates used to 

calibrate the combined model. Simulation results are summarized in graphical engineering response curves that show the 

relations among effluent phosphorus concentrations, thickness 

of the unsaturated zone, and length of the phosphorus plume 

in the aquifer. The response curves can facilitate determination 

of the disposal location and the degree of phosphorus-removal 

treatment needed to meet low-impact disposal objectives. 

Sensitivity analysis is included to quantify the effects on the 

response curves of the number of phosphorus binding sites on 

the aquifer solids, duration of disposal, and other aquifer and 

disposal conditions.

The dependence of phosphorus plume length on unsaturated-zone thickness and wastewater phosphorus concentration 

was simulated by combining 1-D unsaturated zone and 3-D 

saturated-zone PHAST computer models. Grid-element sizes 

were 3.33 ft in the unsaturated zone and 100 by 100 ft by 10 ft 

in the saturated zone. The simulations gave estimates of phosphorus plume length for use in determining setback distances 

for onland disposal of wastewater required for protection of 

surface water from eutrophication by subsurface transport of 

phosphorus. Base-case disposal conditions were 3 gal/ft2/d 

of effluent disposed in a 100- by 200-ft infiltration bed for 

50 yrs. The length of the phosphorus plume was defined as the 

maximum extent of water with 0.015 mg/L or more dissolved 

phosphorus downgradient from the infiltration bed during 

discharge and after discharge cessation; the plume may take 

more than 100 years to reach its maximum length. For the base 

case, phosphorus transport beyond the extent of this contour 

did not exceed 0.18 kg/yr, an amount expected to have little or 

no effect on the aquatic plant growth of most surface water. 

Simulations of phosphorus plume lengths were summarized in a response-curves graph. Simulated plume lengths 

ranged from 200 ft for low-concentration (0.25 mg/L) effluent and thick (50 ft) unsaturated zones to 3,400 ft for highconcentration (14 mg/L) effluent discharging directly into 

the aquifer (unsaturated-zone thickness of 0 ft). Plume length 

was nearly independent of unsaturated-zone thickness at low 

values of phosphorus concentration in the wastewater because 

little or no phosphorus mineral formed at low dissolved 

phosphorus concentrations. At high concentrations of effluent, 

plume length varied from 3,400 ft for 0-ft unsaturated-zone 

thickness to 2,250 ft for 50-ft unsaturated-zone thickness. 

Plume-length dependence on the concentration of phosphorus in the wastewater was also variable, with diminished 

increase in plume length as wastewater phosphorus concentration increased. The diminished response was caused by the 

phosphate-mineral precipitation rate, which increases with 

increase in phosphorus concentration.

Model treatments of flow and equilibrium controlled 

chemistry were likely more accurate than were assumed 

properties of sorption surfaces and assumed rates of kinetically controlled reactions, notably precipitation of phosphate 

minerals. A sensitivity analysis indicated that many chemical 

and aquifer properties, such as hydraulic gradient and model 

width, did not affect the simulated plume length appreciably; 

but duration of discharge, size of infiltration bed, and number 

of sorption sites on aquifer sediments did appreciably affect 

plume length. 

An important result from the simulations is that phosphorus removal by aquifer sediment is more efficient at high 

effluent concentrations of phosphorus discharged than at low 

26 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater
concentrations. High-concentration phosphorus is removed 

more efficiently because the rate of phosphorus removal by 

precipitation of phosphorus minerals is a function of phosphorus concentration; at low phosphorus concentrations, the 

solubility product of minerals is not exceeded, so no mineral 

forms. A consequence of the dependence of mineral formation on phosphorus concentration is that extent of phosphorus 

removal processes in the subsurface is decreased as more 

complete removal efforts are applied in treatment before wastewater disposal. The relation between phosphorus concentration and phosphorus removal rate means that removal of phosphorus in the subsurface from small volumes of wastewater, 

such as are disposed through septic systems, may be highly 

efficient. Phosphorus concentrations and removal rates are 

high in septic-system effluent; volume of wastewater disposed 

is low; consequently plume lengths are short. Less efficient is 

phosphorus removal from effluents with phosphorus concentration that has been decreased by treatment before disposal. 

If phosphorus removal is carried out in the treatment plant, the 

resulting effluent may be too dilute for precipitation reactions 

to occur in the subsurface.

Whereas the results of this investigation are expected 

to have wide application to noncarbonate glacially deposited 

sediments of Massachusetts, simulation of plume length in 

carbonate-mineral sediments indicated that plumes would be 

much longer. Thus, the model results would be inappropriately 

applied to locations containing carbonate-mineral sediments. 

The simulations indicate that less phosphorus is removed in 

carbonate or other high-pH ground-water environments.

Even with the simulations of the response curves for 

phosphorus plume lengths, as defined by the 0.015 mg/L 

concentration contour completed, the relation between plume 

length and protective setback distance is subject to interpretation. Phosphorus does move beyond the point at which the 

0.015-mg/L-concentration contour has stopped, so that a 

determination of protective setback distance must include a 

consideration of whether that continued flux, or some other 

flux amount is protective. Also, model simulations indicate 

that phosphorus plumes do not reach their full extent until 50 

to 200 yrs after disposal cessation, depending on concentration of phosphorus disposed. No phosphorus plume has been 

monitored for that long a period after cessation so that there is 

no way to verify the long-term simulation results.

Acknowledgments

Guidance from David Parkhurst on computer modeling, 

and calibration results for phosphorus from Ken Stollenwerk 

and for organic carbon from Dick Smith were much used and 

appreciated.

References

Colman, J.A., and Friesz, P.J., 2001, Geohydrology and 

limnology of Walden Pond, Concord, Massachusetts: U.S. 

Geological Survey Water-Resources Investigations Report 

01-4137, 61 p.

Dzombak, D.A., and Morel, F.M.M, 1990, Surface complexation modeling—Hydrous ferric oxide: New York, Wiley 

Publishers, 393 p.

Fetter, C.W., 1993, Contaminant hydrology: New York, 

Macmillan, 458 p.

Fletcher, C.A.J., 1991, Computational techniques for fluid 

dynamics, v. 1 (2d ed.): Berlin, Springer-Verlag, 401 p.

Friesz, P.J., 1996, Geohydrology and recharge of stratified drift and streamflow in the Deerfield River Basin, 

Massachusetts: U.S. Geological Survey Water-Resources 

Investigations Report 96-4115, 49 p.

Frimpter, M.H., 1980, Map showing ground-water availability in the northern part of the Connecticut Valley 

urban area, central New England: Miscellaneous Geologic 

Investigations Map I-1074-I, scale 1:125,000.

Frimpter, M.H., 1981, Probable high ground-water levels in 

Massachusetts: U.S. Geological Survey Open-File Report 

80-1205, 2 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, p. 911–924.

Jantrania, A.R., 2000, Wastewater treatment in the 21st century—Technology, operation, management, and regulatory 

issues: Journal of Environmental Health, v. 63, p. 19–25.

Kent, D.B., LeBlanc, D.R., Abrams, R.H., and others, 2000, 

Modeling the influence of variable pH on the transport of 

zinc in a contaminated aquifer using semiempirical surface 

complexation models: Water Resources Research, v. 36, 

p. 3411–3425.

Kipp, K.L., Jr., 1997, Guide to the revised heat and solute 

transport simulator—HST3D: U.S. Geological Survey 

Water-Resources Investigations Report 97-4157, 149 p.

References 27
Kolpin, D.W., Zaugg, S.D., Barber, L.B., Buxton, H.T., 

Furlong, E.T., Meyer, M.T., and Thurman, E.M., 2003, 

Pharmaceuticals, hormones, and other organic wastewater 

contaminants in U.S. streams, 1999–2000—A national 

reconnaissance: Environmental Science and Technology, 

v. 36, p. 1202–1211.

LeBlanc, D.R., 1984, Sewage plume in a sand and gravel aquifer, Cape Cod, Massachusetts, in U.S. Geological Survey 

Water-Supply Paper 2218, p. 28. 

Massachusetts Department of Environmental Protection, 2003, 

Title 5 Septic Systems: accessed June 23, 2003, at http://

www.state.ma.us/dep/brp/wwm/t5about.htm

Masterson, J.P., and Barlow, P.M., 1997, 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 WaterSupply Paper 2447, 79 p.

McCobb, T.D., LeBlanc, D.R., and Hess, K.M., 1999, 

Determination of temporal and spatial variability of hydraulic gradients in an unconfined aquifer using three-point 

triangulation, Cape Cod, Massachusetts, in Morganwalp, 

D.W., and Buxton, H.T., eds., U.S. Geological Survey 

Toxic Substances Hydrology Program—Proceedings of the 

Technical Meeting, Charleston, South Carolina, March 8–

12, 1999—Volume 3 of 3—Subsurface Contamination from 

Point Sources: U.S. Geological Survey Water-Resources 

Investigations Report 99-4018C, p. 349–360. 

McCobb, T.D., LeBlanc, D.R., Walter, D.A., Hess, K.M., 

Kent, D.B., and Smith, R.L., 2003, Phosphorus in a groundwater contaminant plume discharging to Ashumet Pond, 

Cape Cod, Massachusetts, 1999: U.S. Geological Survey 

Water-Resources Investigations Report 02-4306, 69 p.

Moore, J.W., Schindler, D.E., and Scheuerell, M.D., 2003, 

Lake eutrophication at the urban fringe, Seattle region: 

USA, Ambio, v. 32, p. 13–18.

Parker, J.E., and Maberly, S.C, 2000, Biological response 

to lake remediation by phosphate stripping—Control of 

Cladophora: Freshwater Biology, v. 44, p. 303–309. 

Parkhurst, D.L., and Appelo, C.A.J., 1999, User’s guide to 

PHREEQC (Version 2)—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse 

geochemical modeling: U.S. Geological Survey WaterResources Investigations Report 99-4259, 312 p.

Parkhurst, D.L., Kipp, K.L., Engesgaard, Peter, and Charlton, 

S.R., 2004, PHAST—A program for simulating groundwater flow, solute transport, and multicomponent geochemical reactions: U.S. Geological Survey Techniques and 

Methods 6-A8, 154 p.

Parkhurst, D.L., Stollenwerk, K.G., and Colman, J.A., 2003, 

Reactive-transport simulation of phosphorus in the sewage 

plume at the Massachusetts Military Reservation, Cape Cod, 

Massachusetts: U.S. Geological Survey Water-Resources 

Investigations Report 03-4017, 33 p.

Robertson, W.D., 2003, Enhanced attenuation of septic system 

phosphate in noncalcareous sediments: Ground Water, v. 41, 

p. 48–56.

Scheffer, M., 1998, The ecology of shallow lakes: New York, 

Chapman and Hall, 357 p.

Stollenwerk, K.G., 1996, Simulation of phosphate transport in sewage-contaminated groundwater Cape Cod, 

Massachusetts: Applied Geochemistry, v. 11, p. 317–324.

Stollenwerk, K.G., and Parkhurst, D.L., 1999, Modeling the 

evolution and natural remediation of a ground-water sewage 

plume, in Morganwalp, D.W., and Buxton, H.T., eds., U.S. 

Geological Survey Toxic Substances Hydrology Program—

Proceedings of the Technical Meeting, Charleston, South 

Carolina, March 8–12, 1999, Volume 3 of 3, Subsurface 

Contamination from Point Sources: U.S. Geological Survey 

Water-Resources Investigations Report 99-4018C, 

p. 371–381.

U.S. Environmental Protection Agency, 1998, National Water 

Quality Inventory: 1998 Report to Congress, State Fact 

Sheets, Chapter 12 State and Territory Summaries, 

p. 318–319.

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., Rea, B.A., Stollenwerk, K.G., and Savoie, 

Jennifer, 1996, Geochemical and hydrologic controls 

on phosphorus transport in a sewage-contaminated 

sand and gravel aquifer near Ashumet Pond, Cape Cod, 

Massachusetts: U.S. Geological Survey Water-Supply Paper 

2463, 89 p.

Weiskel, P.K., Howes B.L., and Heufelder, G.R., 1996, 

Coliform contamination of a coastal embayment—Sources 

and transport pathways: Environmental Science & 

Technology v. 30, p. 1872–1881.

Zarriello, P.J., 2002, Effects of water-management alternatives 

on streamflow in the Ipswich River Basin, Massachusetts: 

U.S. Geological Survey Open-File Report 01-483, 30 p.

28 Response Curves for Phosphorus Plume Lengths from Onland Disposal of Wastewater