Reactive-Transport Simulation of 

Phosphorus in the Sewage Plume at the 

Massachusetts Military Reservation,

Cape Cod, Massachusetts

By DAVID L. PARKHURST, KENNETH G. STOLLENWERK, and JOHN A. COLMAN

Northborough, Massachusetts

2003

TOXIC SUBSTANCES HYDROLOGY PROGRAM

In cooperation with the

AIR FORCE CENTER FOR ENVIRONMENTAL 

EXCELLENCE

Water-Resources Investigations Report 03-4017
For additional information write to:

Chief, Massachusetts-Rhode Island District

U.S. Geological Survey

10 Bearfoot Road

Northborough, Massachusetts 01532

or visit our Web site at

http://ma.water.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.

Copies of this report can be purchased 

from:

U.S. Geological Survey

Branch of Information Services

Box 25286, Federal Center

Denver, CO 80225-0286

U.S. DEPARTMENT OF THE INTERIOR

GALE A. NORTON, Secretary

U.S. GEOLOGICAL SURVEY

Charles G. Groat, Director
Contents III

CONTENTS

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

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

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

Methods................................................................................................................................................................ 6

Reactive-Transport Model .......................................................................................................................... 6

Computer Resources ................................................................................................................................... 7

Column Experiments................................................................................................................................... 7

Sorption/Desorption Column Experiments........................................................................................ 7

Fitting Surface-Sorption Parameters from Sorption/Desorption Experimental Data ........................ 8

Phosphorus Desorption Experiments................................................................................................. 9

Analytical Methods............................................................................................................................ 10

Reactive-Transport Simulations........................................................................................................................... 10

Base-Case Model ........................................................................................................................................ 10

Geometry of the Model Grid.............................................................................................................. 10

Hydraulic Properties .......................................................................................................................... 12

Hydrologic Boundary Conditions ...................................................................................................... 12

Chemical Initial and Boundary Conditions........................................................................................ 13

Chemical Reactions in the Sewage Plume......................................................................................... 13

Sorption Reactions.................................................................................................................... 14

Iron and Manganese Minerals................................................................................................... 14

Organic Decomposition Reactions ........................................................................................... 14

Iron Phosphate Minerals ........................................................................................................... 16

Irreversible Denitrification........................................................................................................ 16

Fitting Surface-Sorption Parameters to Sorption/Desorption Experimental Data............................. 16

Sensitivity Simulations ...................................................................................................................... 16

Revised Model ............................................................................................................................................ 21

Modifications to the Base-Case Model.............................................................................................. 22

Sorption and Subsequent Reaction of Organic Carbon ............................................................ 22

Irreversible Denitrification........................................................................................................ 23

Non-Desorbable Phosphorus .................................................................................................... 23

Infiltration at the Surface of the Aquifer................................................................................... 23

Leaky Boundary Condition for the Southern Boundary........................................................... 24

Number of Sorption Sites.......................................................................................................... 24

Oxygen Concentration in Sewage Effluent............................................................................... 26

Comparison of Measured and Simulated Phosphorus Concentrations .............................................. 26

Phosphorus Distribution in 1993 .............................................................................................. 26

Phosphorus at Site S469............................................................................................................ 26

Transport of Phosphorus to Ashumet Pond ....................................................................................... 28

Limitations of the Revised Model...................................................................................................... 30

Summary and Conclusions................................................................................................................................... 31

References Cited .................................................................................................................................................. 32
IV Contents

FIGURES

1–4. Maps showing: 

1. Location of disposal beds in the study area at the Massachusetts Military Reservation 

sewage-treatment plant near Ashumet Pond, Massachusetts, and extent of the sewage 

plume in Ashumet Valley as of 1993–94................................................................................... 3

2. Areal distribution of maximum dissolved phosphorus concentrations in ground water 

near Ashumet Pond, August to November 1993........................................................................ 4

3. Altitude of the water table and approximate direction of ground-water flow in January 

1994 near Ashumet Pond ........................................................................................................... 5

4. Model grid with 50-meter spacing for reactive-transport model of the area near 

Ashumet Pond ............................................................................................................................ 11

5–7. Graphs showing: 

5. Results of column experiments of phosphorus sorption and desorption on Cape Cod 

sediments and results of reactive-transport simulations based on fitted parameters.................. 17

6. Sensitivity of phosphorus load transported to Ashumet Pond to variation in selected 

reactive-transport model parameters .......................................................................................... 20

7. Measured and simulated phosphorus concentrations for a column desorption 

experiment on a sediment core from the sewage-contaminated zone of the aquifer ................. 25

8. Maps showing (A) measured and (B) simulated phosphorus concentrations in ground 

water from August to November 1993 near Ashumet Pond.............................................................. 27

9, 10. Graphs showing:

9. Measured phosphorus concentrations in July 1996 and July 1999 in multilevel sampler 

S469M1 and simulated phosphorus concentrations for the vertical set of nodes 

nearest to the location of disposal beds 5–8............................................................................... 28

10. Simulated load of phosphorus transported to Ashumet Pond and estimated load of phosphorus to 

disposal beds during the period 1936 to 2055, calculated with revised model parameters and 

revised model parameters excluding vivianite formation ..........................................................29

TABLES

1. Properties used in reactive-transport simulation of the phosphorus sorption/desorption 

column experiments........................................................................................................................... 9

2. Composition of synthetic ground water used in phosphorus desorption column experiment ........... 9

3. Base-case parameters for three-dimensional reactive-transport modeling ........................................ 12

4. Estimated loading of disposal beds during the period of sewage disposal........................................ 13

5. Solution compositions for uncontaminated ground water, treated sewage effluent, and 

rainwater ............................................................................................................................................ 14

6. Major chemical reactions used in reactive-transport simulations...................................................... 15

7. Parameters derived from fitting sorption/desorption column experiments........................................ 18

8. Base-case and perturbed values for sensitivity simulations with the three-dimensional 

reactive-transport model .................................................................................................................... 19

9. Comparison of parameters and processes for base-case and revised models .................................... 22
Conversion Factors, Vertical Datum, Abbreviations, and Acronyms V

CONVERSION FACTORS, VERTICAL DATUM, ABBREVIATIONS,

AND ACRONYMS

CONVERSION FACTORS

VERTICAL DATUM

Sea level: Vertical coordinate information is referenced to the North American Vertical Datum of 

1929 (NGVD29).

ABBREVIATIONS AND ACRONYMS

g/L grams per liter

kg/yr kilograms per year

MMR Massachusetts Military Reservation

μm micrometer

μmol/L micromoles per liter

m2/g square meter per gram

m2/L square meter per liter

USGS U.S. Geological Survey

Multiply By To obtain

centimeter (cm) 0.3937 inch

gram (g) 0.03527 ounce, avoirdupois

kilometer (km) 0.6214 mile

meter (m) 3.281 foot

meters per day (m/d) 3.281 foot per day

meters per second (m/s) 3.281 foot per second

milliliter (mL) 0.0002642 gallon

millimeter (mm) 0.03937 inch

million gallons per day (Mgal/d) 0.003785 cubic meters per day

Abstract 1

Reactive-Transport Simulation of 

Phosphorus in the Sewage Plume at the 

Massachusetts Military Reservation, 

Cape Cod, Massachusetts

By David L. Parkhurst, Kenneth G. Stollenwerk, and John A. Colman

ABSTRACT

The subsurface transport of phosphorus 

introduced by the disposal of treated sewage 

effluent to ground-infiltration disposal beds at 

the Massachusetts Military Reservation on 

western Cape Cod was simulated with a threedimensional reactive-transport model. The 

simulations were used to estimate the load of 

phosphorus transported to Ashumet Pond during 

operation of the sewage-treatment plant—from 

1936 to 1995—and for 60 years following 

cessation of sewage disposal. The model 

accounted for spatial and temporal changes in 

water discharge from the sewage-treatment plant, 

ground-water flow, transport of associated 

chemical constituents, and a set of chemical 

reactions, including phosphorus sorption on 

aquifer materials, dissolution and precipitation of 

iron- and manganese-oxyhydroxide and iron 

phosphate minerals, organic carbon sorption and 

decomposition, cation sorption, and irreversible 

denitrification.

The flow and transport in the aquifer were 

simulated by using parameters consistent with 

those used in previous flow models of this area of 

Cape Cod, except that numerical dispersion was 

much larger than the physical dispersion estimated 

in previous studies. Sorption parameters were fit to 

data derived from phosphorus sorption and 

desorption laboratory column experiments. Rates 

of organic carbon decomposition were adjusted to 

match the location of iron concentrations in an 

anoxic iron zone within the sewage plume. The 

sensitivity of the simulated load of phosphorus 

transported to Ashumet Pond was calculated for a 

variety of processes and input parameters. 

Model limitations included large 

uncertainties associated with the loading of the 

sewage beds, the flow system, and the chemistry 

and sorption characteristics in the aquifer. The 

results of current model simulations indicate a 

small load of phosphorus transported to Ashumet 

Pond during 1965–85, but this small load was 

particularly sensitive to model parameters 

that specify flow conditions and the chemical 

process by which non-desorbable phosphorus is 

incorporated in the sediments. The uncertainties 

were large enough to make it difficult to determine 

whether loads of phosphorus transported to 

Ashumet Pond in the 1990s were greater or less 

than loads during the previous two decades.

The model simulations indicate substantial 

discharge of phosphorus to Ashumet Pond after 

about 1965. After the period 2000–10 the 

simulations indicate that the load of phosphorus 

transported to Ashumet Pond decreases 

continuously, but the load of phosphorus remains 

substantial for many decades. The current 

simulations indicate a peak in phosphorus 
2 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

discharge to Ashumet Pond of about 1,000 

kilograms per year during the 1990s; however, 

comparisons of simulated phosphorus 

concentrations with measured concentrations in 

1993 indicate that the peak in phosphorus load 

transported to Ashumet Pond may be larger and 

moving more quickly in the model simulations 

than in the aquifer.

The results of the three-dimensional 

reactive-transport simulations are consistent 

with the loading history, experimental laboratory 

data, and field measurements. The results of the 

simulations adequately reproduce the spatial 

distribution of phosphorus concentrations 

measured in 1993, the magnitude of changes in 

phosphorus concentration with time in a profile 

near the disposal beds following cessation of 

sewage disposal, the observed iron zone in the 

sewage plume, the approximate flow of treated 

sewage effluent into Ashumet Valley, and 

laboratory-column data for phosphorus sorption 

and desorption. 

INTRODUCTION

Secondarily treated sewage effluent was 

discharged from 1936 to 1995 to sand beds at the 

Massachusetts Military Reservation (MMR), Cape 

Cod, Massachusetts, which is on a broad, glacial 

outwash plain on western Cape Cod (LeBlanc, 1984; 

Walter and others, 1996). The disposal of treated 

sewage effluent has created a plume of sewagecontaminated ground water in the local aquifer that has 

moved southward to Ashumet Pond and into Ashumet 

Valley (fig. 1). Disposal of treated sewage effluent 

ended in 1995, but effects of the long-term disposal 

operations remain (LeBlanc and others, 1999).

At the time of cessation of sewage disposal, 

phosphorus concentrations were greatest in that part of 

the plume just downgradient from the southeastern 

disposal beds—between the beds and Ashumet Pond 

(fig. 2, Walter and others, 1996). The elongation of the 

phosphorus-concentration contours toward Ashumet 

Pond indicates flow from the disposal beds toward 

Ashumet Pond. Ground-water flow directions in the 

aquifer (Walter and others, 1996) also show that ground 

water flows from the southeastern disposal beds to 

Ashumet Pond. However, the flow directions indicate 

that ground-water flow from the western disposal beds 

skirts the Pond and moves down Ashumet Valley 

(fig. 3). Flow directions in the vicinity of Ashumet 

Pond may vary by as much as 15o (Walter and others, 

1996). 

The magnitude of the phosphorus load 

transported to Ashumet Pond directly affects the 

health of the pond ecosystem. Estimates of the past 

and future loads of phosphorus are thus important 

information for planning a management strategy 

for Ashumet Pond. To address these needs and 

other environmental issues, the U.S. Geological 

Survey, in cooperation with the Air Force Center for 

Environmental Excellence, has studied the distribution 

and transport of phosphorus and other contaminants in 

the sewage plume at the MMR.

PURPOSE AND SCOPE

The purpose of this report is to describe a 

reactive-transport model that was developed to 

simulate phosphorus in the sewage plume at the MMR. 

Loads of phosphorus transported to Ashumet Pond 

were simulated over the course of the sewagetreatment-plant operation at MMR and over a 60-year 

time period after cessation of sewage disposal. The 

simulated phosphorus loads transported to Ashumet 

Pond are needed to formulate effective management 

strategies for Ashumet Pond. The report describes the 

selection of a set of chemical reactions for phosphorus 

and cation sorption, sorption and degradation of 

organic carbon, sequestration of non-desorbable 

phosphorus, and mineral equilibria; the evaluation of 

parameters for phosphorus and cation sorption to 

aquifer sediments by using data from laboratorycolumn experiments; the estimation of the rate of 

organic degradation reactions in the sewage plume; the 

testing of the sensitivity of the model to model 

parameters; the comparison of measured phosphorus 
Purpose and Scope 3

Ashumet

Pond

EXPLANATION

SEWAGE PLUME

SEWAGE DISPOSAL BEDS AND IDENTIFIERS

WELL-CLUSTER SITES AND LOCAL IDENTIFIERS

1

1

2

3 4

5

6

7

8

9

10

11

12

13

14

15

16

17 18 19 21

22

24 23

25

F230

F231 F232

F236

F239

F300

F347 F383

F388

F422

F432

F564

F567

S316

S344

S434

F587

S593

S590

F589 F588 F424

S313

S317

S318 S314

S467

S423

S395

S436

S437

F343 F378F379

F570

S438; S440

F235

F238

F254

F254

F267 F348

S315

F582

F583

F584

F592 F237 F244

S469

Ashumet Valley

70o33'00

41o38'30"

70o32'00"

41o37'30"

Modified from Savoie and LeBlanc (1998)

41 ̊30'

42 ̊30'

73 ̊

71 ̊

0

0

50

50

MILES

KILOMETERS

MASSACHUSETTS STUDY

AREA

0

0

500

1,500

1,000 METERS

3,000 FEET

Figure 1. Location of disposal beds in the study area at the Massachusetts Military Reservation sewage-treatment plant near Ashumet 

Pond, Massachusetts, and extent of the sewage plume in Ashumet Valley as of 1993–94. 
4 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

Ashumet Valley

Ashumet

Pond

3

30

120

150

90

120

150 180

30

3

60

90

3

EXPLANATION

0

0

500

1,500

1,000 METERS

3,000 FEET

LINES OF EQUAL MAXIMUM DISSOLVED 

 PHOSPHORUS CONCENTRATION—Dashed 

 where approximately located. Interval, in 

 micromoles per liter, is variable.

WELL-CLUSTER SITES

MODELED DISPOSAL BEDS

AND IDENTIFIERS

SEWAGE DISPOSAL BEDS 

Modified from Walter and others (1996)

70o33'00" 

41o38'30"

70o32'00" 

41o37'30"

Figure 2. Areal distribution of maximum dissolved phosphorus concentrations in ground water near Ashumet Pond, Massachusetts, 

August to November 1993.
Purpose and Scope 5

Ashumet Valley

Ashumet

Pond

70o33'00"

41o38'30"

70o32'00"

41o37'30"

SEWAGE DISPOSAL BEDS

WATER-TABLE CONTOURS—Show altitude of water table. Dashed where 

 approximately located. Contour interval is 0.15 meter. Datum is NGVD29.

APPROXIMATE DIRECTION OF GROUND-WATER FLOW 

Modified from Walter and others (1996)

EXPLANATION

14.65

14.50

14.35

14.20

14.05

13.90

13.75

13.60

13.45

13.30

13.15

13.45

0

0

500

1,500

1,000 METERS

3,000 FEET

Figure 3. Altitude of the water table and approximate direction of ground-water flow in January 1994 near Ashumet Pond, 

Massachusetts.
6 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

concentrations to simulated phosphorus 

concentrations; and the estimation of the loads of 

phosphorus transported to Ashumet Pond.

A base-case model was developed on the basis of 

chemical analyses of samples taken from the aquifer 

in 1993, ground-water data from 1996, and results 

from laboratory-column experiments. A revised model 

was then developed by the inclusion of additional 

chemical processes and by the adjustment of the basecase model parameters on the basis of sensitivity 

simulations and comparisons with additional 

laboratory-column experiments. The reactivetransport simulations produced results for dissolved 

concentrations of each major element plus nitrogen, 

phosphorus, dissolved organic carbon, pH, and redox 

conditions; sorbed concentrations of phosphorus 

and major cations; and concentrations of solid 

constituents—iron and manganese oxyhydroxides, 

iron-phosphate minerals, and sorbed organic carbon. 

For the purposes of this report, only the simulation 

results for phosphorus are discussed in detail. 

METHODS

This section describes a three-dimensional 

reactive-transport model used to simulate transport of 

phosphorus downgradient from the disposal beds at 

the MMR. Model definitions include geometry and 

boundary conditions, initial conditions, and selection 

of chemical reactions. A description of computer 

resources used for the simulations is provided. Finally, 

this section describes experimental methods for 

laboratory experiments to determine phosphorus 

sorption and desorption from sediment columns and 

numerical methods for the determination of 

equilibrium constants from the column-experiment 

results. 

Reactive-Transport Model

The modeling described in this report expands 

on previous work by Stollenwerk and Parkhurst (1999) 

in which PHREEQC, a one-dimensional model, was 

used to simulate reactive transport in the sewage 

plume. Because of the simple mineralogy and the 

rapidity of the flow system in the aquifer in the vicinity 

of the MMR, few reactions need to be considered in a 

reactive-transport model of phosphorus. Carbonate 

minerals are absent from the aquifer sediment, and the 

primary minerals that compose the sediments—quartz 

and feldspars (Barber, 1990)—are relatively inert on 

the time scale of the flow system. Minor amounts of 

glauconite and clay minerals are present (Barber, 

1990), which allow for a small amount of cation 

exchange or sorption. Iron and manganese 

oxyhydroxide minerals and coatings of the aquifer 

sediments can be assumed to react rapidly under 

anoxic conditions. Iron and aluminum oxyhydroxide 

minerals provide large surface areas that can reversibly 

sorb phosphate and other ions, and can buffer the pH 

of the ground water to a significant extent. Phosphate 

can be sequestered in iron and aluminum oxyhydroxide 

minerals as solid solutions or in authigenic iron and 

aluminum phosphate minerals. 

The chemical reactions in the contaminated part 

of the aquifer are driven primarily by the loading of 

treated sewage effluent. The phosphorus in the effluent 

is the source of the phosphorus plume in the aquifer. 

Organic carbon introduced with the effluent enters the 

ground-water system and is readily reactive, but this 

organic carbon also may sorb on aquifer sediments 

and continue to react over an extended period of time. 

All of the reactions that occur rapidly can affect the 

transport of phosphate, in large part because they 

affect the pH of the ground water. The buffering of 

pH by these reactions directly affects the sorption of 

phosphorus and the formation of phosphate minerals.

Reactive-transport models are a recent 

development and modelers do not have a large body of 

work from which to draw. The combined capability to 

model flow, transport, and chemical reactions provides 

a systematic approach for studying ground-water 

processes. The model developed in this report attempts 

to include all of the significant processes affecting the 

transport of phosphorus to Ashumet Pond and to 

provide a comprehensive description of the evolution 

of ground-water and sediment chemistry. Even though 

large uncertainties are associated with the modeling 

results, a reactive-transport model is the only 

systematic method available to estimate the time 

dependency of the loads of phosphorus transported to 

Ashumet Pond and to assess the sensitivity of the load 

estimate to various chemical and physical processes.

The three-dimensional reactive-transport model 

used to develop the base-case and revised models was 

PHAST (Parkhurst and others, 1995; Parkhurst and 

Kipp, 2002). PHAST couples capabilities to simulate 

ground-water flow and associated solute transport with 
Methods 7

the capability to simulate chemical reactions. The 

ground-water flow and transport parts of PHAST are 

derived from HST3D (Kipp, 1997) and the chemical 

reaction part of the model is derived from PHREEQC 

(Parkhurst and Appelo, 1999). In this study, PHAST 

was used (1) to estimate sorption properties of the 

aquifer sediments, (2) to estimate the rate constant for 

organic decomposition reactions in the aquifer, and 

(3) to simulate the evolution of the sewage plume over 

time, specifically, to calculate loads of phosphorus 

transported to Ashumet Pond.

The hydraulic parameters and boundary 

conditions for the three-dimensional reactive-transport 

model were taken from models developed for previous 

studies (LeBlanc, 1984; LeBlanc and others, 1991; 

Garabedian and others, 1991; Walter and others, 

1996; Masterson and Walter, 2000). Although the 

USGS computer code MODFLOW (McDonald and 

Harbaugh, 1988; Harbaugh and McDonald, 1996) was 

used to simulate ground-water flow in these studies, 

PHAST solves the same partial differential flow 

equations as MODFLOW and uses the same hydraulic 

parameters. Simulated ground-water heads for the 

model domain used in this study were not adjusted 

by modifying parameters in a formal calibration 

procedure, but the simulated heads used in this study 

approximate heads simulated for previous studies. The 

flow system simulated by PHAST approximates the 

flow system in the aquifer in the study area and is 

appropriate for the reactive-transport simulations 

described in this report.

Computer Resources

Calculations for the reactive-transport 

simulations described in this report were run with a 

parallel-computing version of PHAST on a Beowulf 

cluster with up to 14 dual-processor personal 

computers. Individual processor speeds ranged from 

0.5 to 1.8 gigahertz and each computer had 0.75 

or 1.0 gigabytes of memory. Local-area-network 

connection was by dual 100-megabit Ethernet with an 

effective data-transfer rate of about 150 megabits per 

second. The parallel-computing version of PHAST 

(Parkhurst and Kipp, 2002) used a single processor to 

run the flow and transport calculations, whereas the 

chemical calculations for the nodes were distributed 

among the available processors with a load-balancing 

strategy. Calculations for the base-case model required 

about two hours of elapsed time on a cluster of 20 

processors.

Column Experiments

Sorption/desorption column experiments 

(Stollenwerk, 1996) were performed with pristine 

sediments unaffected by sewage effluent to evaluate 

phosphorus and cation-sorption parameters, which 

included equilibrium constants for sorption of 

phosphorus, cations, and hydrogen ion, and the 

concentration of sorption sites for phosphorus and 

cations. Desorption column experiments in which 

phosphorus was leached from contaminated sediment 

were simulated with PHAST to test the validity of the 

sorption parameters for simulating desorption of 

phosphorus from sewage-contaminated sediments in 

contrast to sorption and desorption from pristine 

sediments.

Sorption/Desorption Column 

Experiments

The sorption/desorption column experiments 

were designed to evaluate sorption and desorption of 

phosphorus from uncontaminated aquifer sediment. 

Sediment was collected 1.3–1.6 m below land surface 

from a trench excavated in the unsaturated zone at a 

gravel pit about 200 m south of the sewage disposal 

beds. Previous studies have shown that sorption 

properties of this sediment are similar to sediment 

from the saturated zone (Stollenwerk, 1995).

Uncontaminated ground water was collected 

from a well that was screened above the sewage plume. 

Sewage-contaminated ground water was collected from 

a well that was screened in the suboxic zone of the 

sewage plume. Unfiltered ground water was collected 

in polyethylene bottles that were stored in ice during 

shipment. All water was filtered in the laboratory 

through a 0.4-μm pore-size, polycarbonate-membrane 

filter and refrigerated until needed. The suboxic ground 

water became saturated with oxygen during the 

experiments; however, because some oxygen was 

already present in the ground water and the 

experiments remained oxic throughout, the additional 

dissolved oxygen introduced by air equilibration is 

expected to have had little effect on the experiment.
8 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

Replicate Plexiglas columns, approximately 

30 cm in length with 2.5-cm inside diameter, were 

filled with the homogenized fraction of the sediment 

sample less than 2 mm in size. This size fraction 

accounted for 90 percent by weight of the sediment. 

A peristaltic pump was used to drive solutions through 

the columns at a velocity of about 0.4 m/d, which is 

approximately the velocity of ground water in this part 

of the aquifer. Leachate was collected in a fraction 

collector that was sealed and humidified to minimize 

evaporation. Pre-sewage-plume conditions were 

established by rinsing columns with uncontaminated 

ground water until the chemistry of the column 

leachate stabilized at that of the eluent. A bromide 

tracer was then added to estimate longitudinal 

dispersivities of the columns for transport simulations. 

After the bromide was eluted from the 

columns, the eluent solution was changed to sewagecontaminated ground water. One set of experiments 

used the phosphorus concentration of 16 μmol/L that 

was present in the sewage-contaminated ground water 

collected for the experiment. In a second set of 

experiments, this ground water was spiked to 196 

μmol/L phosphorus. The two concentrations represent 

the range of phosphorus concentrations measured in 

contaminated ground water in the plume. Sewagecontaminated ground water was eluted through the 

columns until the leachate concentration of phosphorus 

approached the eluent concentration. The eluent 

solution was changed to uncontaminated ground water 

after 115 pore volumes for the 16-μmol/L experiment 

and after 48 pore volumes for the 196-μmol/L 

experiment to evaluate desorption. 

Fitting Surface-Sorption Parameters from 

Sorption/Desorption Experimental Data

PHAST was used to simulate flow and 

chemical reactions for the sorption/desorption column 

experiments. A version of the parameter-estimation 

code, UCODE (Poeter and Hill, 1999), modified 

to use multiple computer processors, was used to fit 

chemical parameters for phosphorus and cation 

sorption on aquifer sediments.The pH and phosphorus 

concentration in the column leachate were fit by 

adjusting surface-complexation parameters. The 

data included 432 phosphorus and 180 pH analyses 

from column 1 (16-μmol/L phosphorus eluent) and 

135 phosphorus and 31 pH analyses from column 3 

(196-μmol/L phosphorus eluent). The phosphorus data 

were weighted by using a standard deviation equal to 

the square root of the concentration (μmol/L) for each 

point, or 0.1 if the analysis was below the detection 

limit of the chemical analytical method; the pH data 

were weighted by using a standard deviation of 0.2 for 

each point.

Table 1 lists properties of the column used in 

phosphorus sorption/desorption experiments. Length 

and porosity were measured. Mass of solid per liter of 

water was calculated by using an estimated density 

of 2,650 g/L for solids and porosity of 0.33 to give 

5,380 g solid/L water. Surface area per liter of water 

was calculated by using an average measured specific 

surface area of 0.3 m2/g solid (Stollenwerk, 1995).

Iron- and manganese-oxyhydroxide minerals 

were specified to be present in the uncontaminated 

aquifer in concentrations representative of those 

measured by chemical extraction of the solids, 

80,000 μmol/L of iron oxyhydroxide and 700 μmol/L 

of manganese oxide (Coston and others, 1995). Ferric 

hydroxide was chosen as the representative iron phase, 

which was assumed to react to equilibrium with the 

aqueous phase throughout the simulations. The 

equilibrium constant for amorphous ferric hydroxide 

was taken from the phreeqc.dat database of PHREEQC 

[Fe(OH)3(a)] (Parkhurst and Appelo, 1999). Pyrolusite 

was chosen as the representative manganese oxide and 

was assumed to react to equilibrium with the aqueous 

phase, unless the mineral was completely removed by 

dissolution. The equilibrium constant for pyrolusite in 

the database phreeqc.dat was used for all simulations 

(Parkhurst and Appelo, 1999). 

Phosphorus sorption was modeled as a surfacecomplexation process (Dzombak and Morel, 1990; 

Parkhurst and Appelo, 1999, p. 28–29) that proceeds 

by the reaction:

, (1)

where represents a phosphorus-sorption site on 

the aquifer solids. The surface complexation massaction equation includes an electrostatic term to 

account for the work done in moving a charged species 

to the charged surface (Parkhurst and Appelo, 1999, 

p. 14–15). The equilibrium constant for this reaction 

was fit to the data from the sorption/desorption column 

experiments. Although many different combinations of 

chemical sorption reactions and physical parameters 

were tried, the final fits were made by using equation 1. 

No competing anion sorption reactions were 

considered for phosphorus sorption, except for 

protonation and deprotonation of the sorption site.

SiteOH PO3-

4 2H+ SiteHPO4

- + + = + H2O

Site
Methods 9

Phosphorus Desorption Experiments

In the desorption column experiment, the 

desorbable fraction of phosphorus associated with 

sewage-contaminated sediment was measured. The 

purpose was two-fold: (1) to evaluate desorption of 

phosphorus from sediment that had been exposed to 

sewage effluent for many years, in contrast to 

desorption measured during the sorption/desorption 

experiments, which involved pristine sediments 

exposed to sewage-contaminated water for less than 90 

days, and (2) to test the accuracy of the model 

parameters for simulating phosphorus desorption from 

sewage-contaminated aquifer sediments. Cores were 

collected in April 2000 from contaminated parts of the 

aquifer downgradient from the disposal beds. A total of 

17 cores was collected in plastic liners from locations 

between the disposal beds and Ashumet Pond. The 

cores were stored on ice during shipment and 

refrigerated until needed. The only desorption column 

experiment described in this report is for the core 

collected from site S469 at the southern edge of the 

disposal beds (site S469 in fig. 1). This core was 

collected from a depth of 15.5–15.9 m above 

NGVD29. The 33-cm-long core was fitted with 

endcaps containing inlet and outlet ports. Synthetic 

ground water that was representative of the major-ion 

composition of uncontaminated ground water in the 

aquifer (table 2) was eluted through the core. Measured 

porosity was 0.3 and an average flow velocity of 0.4 

m/d was used (range = 0.3–0.5 m/d). Initially the eluent 

pH was set at 5.78; after 50 pore volumes, the eluent 

pH was increased to 7.42; and after 150 pore volumes, 

the eluent pH was decreased back to 5.78.

The reason for the pH changes was to evaluate the 

changes in phosphorus concentrations in response to 

changes in eluent pH. Phosphorus concentration and 

pH of the leachate were monitored throughout the 

experiments.

Desorption of phosphorus in the experiment with 

sewage-contaminated sediment is more representative 

of processes occurring in the aquifer after cessation of 

sewage disposal than the sorption/desorption column 

experiment with pristine sediment. However, the 

desorption experiments require several months to run. 

To obtain data from a larger number of core samples in 

a reasonable time frame, batch desorption experiments 

were conducted at pH of 10. A pH of 10 is high enough 

to maximize phosphorus desorption, yet low enough to 

minimize dissolution of iron and aluminum oxides, the 

principal sorbents for phosphorus in this aquifer.

Total desorbable and acid-extractable 

phosphorus were measured in samples from all 17 

cores. For total desorbable phosphorus at pH of 10, 

15 g of wet core material was equilibrated with 

25 mL of uncontaminated ground water in 40-mL 

polyethylene centrifuge tubes, and pH was adjusted to 

10 with NaOH. Experiments were run in triplicate. The 

suspensions were equilibrated for 2 days on an endover-end rotating mixer and then centrifuged for 20 

minutes at 15,000 revolutions per minute. The 

centrifugate was analyzed for phosphorus and pH. This 

process was repeated until phosphorus desorption 

became negligible, usually after eight extractions. The 

total concentration of phosphorus desorbed from each 

Table 1. Properties used in reactive-transport simulation of the 

phosphorus sorption/desorption column experiments

[g/L, grams per liter; m, meter; μmol/L, micromoles per liter; m2/g, square 

meter per gram]

Column properties Value

Length, m 0.3

Porosity, unitless .33

Density of solid, g/L solid 2,650

Mass of solid per liter of water, g/L water 5,380

Surface area, m2/g solid .3

Surface area, m2/L water 1,614

Fe(OH)3(a), μmol/L water 80,000

MnO2, μmol/L water 700

Table 2. Composition of synthetic ground water used in phosphorus 

desorption column experiment

[Concentrations are in micromoles per liter, except for pH, in standard units]

Constituent Concentration

pH 5.78

Calcium 72

Magnesium 33

Sodium 240

Potassium 62

Chloride 370

Sulfate 45

Bicarbonate 80

Phosphorus <.6

Iron <.05
10 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

core was then calculated. Total acid-extractable 

phosphorus was analyzed on separate splits of the core 

samples, and non-desorbable phosphorus was 

calculated by difference. For acid-extractable 

phosphorus, 15 g of wet core material was equilibrated 

with 25 mL of 1N HCl for 24 hours on an end-over-end 

rotating mixer and then centrifuged for 20 minutes at 

15,000 revolutions per minute. These experiments were 

run in triplicate and involved only one extraction.

Analytical Methods

Phosphorus was analyzed by the 

phosphomolybdate colorimetric procedure (Fishman 

and Friedman, 1985). pH was measured with an Orion 

pH meter and Ross pH electrode. Cations were 

analyzed by inductively coupled atomic-emission 

spectroscopy; anions were analyzed by ion 

chromatography.

REACTIVE-TRANSPORT 

SIMULATIONS

This section describes the base-case 

model, results of fitting model parameters to the 

sorption/desorption experimental data, and a series of 

sensitivity simulations. A revised model is developed 

by making adjustments to the base-case model. The 

revised model was used to make estimates of the load 

of phosphorus transported to Ashumet Pond.

Base-Case Model

The MMR site and the underlying aquifer have 

been extensively studied over the last 20 years, and 

many of the hydraulic parameters for the aquifer are 

well known. Most of the hydraulic properties for the 

aquifer, chemical compositions of ground water and 

sewage effluent, and characteristics of the flow system 

were selected from previous studies of the aquifer 

(LeBlanc, 1984; LeBlanc and others, 1991; Garabedian 

and others, 1991; Walter and others, 1996; Masterson 

and Walter, 2000). 

Chemical reactions and reaction parameters 

were estimated through interpretation of the column 

experimental data and concentration data from the 

aquifer.

Geometry of the Model Grid

In the vicinity of the MMR, water in the aquifer 

flows radially outward from a ground-water high about 

5 km north of the sewage disposal beds (Masterson and 

Walter, 2000). The model for PHAST was chosen to 

simulate flow and transport in a region of the aquifer 

that includes the area of the sewage disposal beds, the 

northeast corner of Ashumet Pond, and the northern 

end of Ashumet Valley. An area 1,600 m in length that 

widens from 800 m at the northern end to 1,600 m at 

the downgradient southern end was chosen for the 

model domain (fig. 4). The widening area allows for 

spreading of the flow system, which is apparent in the 

water-table contours of figure 3 and is consistent with 

regional flow (Masterson and Walter, 2000). The 

eastern and western edges of the domain were chosen 

to align with the predominant flow directions (fig. 3). 

The eastern boundary is a straight line, whereas the 

western edge is stair-stepped progressively farther 

westward from north to south. The model domain is 

45 m thick and extends from 30 m below NGVD29 to 

15 m above NGVD29. 

The model grid is a point-distributed mesh, such 

that the edges of the model are defined by nodes. Cell 

faces are placed halfway between nodes in each 

direction; for cells that border the exterior of the 

domain, cell faces coincide with the node. The model 

grid was discretized with 33 nodes at 50-m spacing in 

the X and Y directions, and 19 nodes at 2.5-m spacings 

in the Z direction (table 3). Each model simulation 

calculates a ground-water head and constituent 

concentrations at each of the nodes and a ground-water 

flux and constituent fluxes across each cell face at each 

time step.

The model grid is rotated counterclockwise 

relative to north and south by about 40 degrees. 

However, north, south, east, and west will be used to 

refer to the top, bottom, right, and left boundaries as 

shown on figure 4. The coordinate axes of the model 

grid are oriented so that positive X is south (the general 

direction of flow) and positive Y is east.
Reactive-Transport Simulations 11

Figure 4. Model grid with 50-meter spacing for reactive-transport model of the area near Ashumet Pond, Massachusetts. 
12 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

Hydraulic Properties

Hydraulic properties are listed in table 3.The 

hydraulic conductivity of the aquifer was specified to 

be 1.06 × 10-3 m/s in the horizontal direction and 0.3 ×

10-3 in the vertical direction; the porosity was specified 

to be 0.39 (Garabedian and others, 1991). It is 

estimated that numerical dispersion in the model is far 

greater than the 1-m dispersivity measured in the 

tracer-test experiments carried out at locations within 

the area of interest (Garabedian and others, 1991). In 

the reactive-transport model, dispersivity was set to a 

small value (0.0001 m), so that numerical dispersion 

from the relatively large grid cells and time steps was 

the major cause of dispersion in the results. The 

numerical dispersivity is approximated by the 

following equation:

, (2)

where is the numerical dispersivity (m), is 

the grid spacing in the direction (m), is the 

time step (day), and is velocity (m/day). The 

numerical dispersivity in the direction (the 

predominant direction of flow) is approximately

 m, 

for 50-m grid cells and half-year (180-day) time steps. 

The effects of numerical dispersion on phosphorus 

loads transported to Ashumet Pond are discussed in the 

“Sensitivity Simulations” section.

Hydrologic Boundary Conditions

Boundary conditions are listed in table 3. The 

northern boundary is assigned a leaky boundary 

condition and uses the approximate location of the 

60-ft water-table contour from Masterson and Walter 

(2000) for the fixed head. This boundary condition 

specifies a fixed head of 18.3 m at a distance of 

1,700 m from the boundary (X = -2,300 m). Heads on 

the leaky boundary will depend on flow conditions 

within the model domain, and the flux of water through 

the boundary will depend on (1) the head gradient 

between the boundary and the fixed head at the 

specified distance and (2) a specified hydraulic 

conductivity of 1.06 × 10-3 m/s.

The southern boundary also is assigned a leaky 

boundary condition and uses the approximate location 

of the 40-ft contour from Masterson and Walter (2000) 

for the fixed head. For the base-case model, the fixed 

head is 12.19 m at a distance of 800 m downgradient 

from the southwest corner of the grid (head = 12.19 m 

at X = 1,800 m, Y = -1,200 m), varying linearly to a 

distance of 1,300 m downgradient from the boundary at 

Y = 0 m (head = 12.19 m at X = 2,300 m, Y = 0 m). 

Thus, the fixed head is specified downgradient at an 

angle relative to the boundary of the model grid from 

Y = -1,200 m to Y = 0 m. The fixed head is parallel to 

the boundary of the model grid from Y = 0 m to Y = 

400 m, and the distance to the fixed head is defined to 

be 1,300 m from the boundary (X = 2,300 m). The 

hydraulic conductivity associated with this boundary is 

1.06 × 10-3 m/s.

Ashumet Pond is represented as a set of 

constant-head nodes in the southeast corner of the 

model grid. Constant-head nodes are defined at an 

elevation of Z = 7.5 m to Z = 15 m in a zone extending 

from X = 600 m to X = 1,000 m, and from Y = 0 m to 

Table 3. Base-case parameters for three-dimensional reactive-transport 

modeling

[m, meters; m/d, meters per day; m/s, meters per second; yr, year]

Parameter Value

Northern leaky boundary

Distance to fixed head, m 1,700

Fixed head, m 18.3

Hydraulic conductivity for boundary, m/s 1.06 × 10-3

Southern leaky boundary

Distance to fixed head, m 800–1,300

Fixed head, m 12.19

Hydraulic conductivity for boundary, m/s 1.06 × 10-3

Specified head for Ashumet Pond, m 13.42

Recharge, m/d 1.8 × 10-3

Porosity 0.39

Hydraulic conductivity, m/s

X - direction 1.06 × 10-3

Y - direction 1.06 × 10-3

Z - direction 0.3 × 10-3

Grid dimensions

X, m 50

Y, m 50

Z, m 2.5

Time step, yr 0.5

Estimated dispersivity, m 60

D ∆X

2

------- V∆T

2 = + ----------

D ∆X

X ∆T

V

X

50

2

----- 0.4 180 ( )

2

+ --------------------- ≅ 60
Reactive-Transport Simulations 13

Y = 400 m. All nodes with X > 600 m, Y > 0 m, and X 

> 7.5 m (the nodes exclusively within Ashumet Pond) 

were also defined to be inactive. All constant-head 

nodes are associated with Ashumet Pond and, 

therefore, flux to Ashumet Pond is equal to the flux 

through the active constant-head nodes of the model.

No-flow boundaries are applied to the eastern 

and the stair-stepped western edges of the model grid, 

including both X and Y faces of the stair steps. The 

bottom of the model grid is also a no-flow boundary. 

No water enters or leaves the system through the

no-flow boundaries.

The top of the model grid is a flux boundary, 

through which rainfall recharges the aquifer. The 

recharge rate from rainfall is specified to be 

1.8 × 10-3 m/d (0.66 m/yr).

A flux boundary condition accounts for recharge 

of sewage effluent through the disposal beds; during 

the period of use of a bed for sewage-effluent disposal, 

the flux of rainwater for the bed is ignored. The 

disposal beds, numbered 1 through 25 (figs. 1 and 4), 

have been loaded for different periods of time and at 

different rates over the course of operation of the 

sewage-treatment plant. The history of total plant 

discharge and usage of the individual disposal beds is 

uncertain. The total discharge from the plant is 

estimated to have varied from 400 to 6,000 m3/d (0.1 to 

1.5 Mgal/d) over the course of operation. For simplicity 

in the model, the disposal beds were defined as five 

100-m × 100-m areas. The current best estimate of the 

history of the plant discharge and loading for these five 

areas is listed in table 4 (D.R. LeBlanc, U.S. 

Geological Survey, written commun., 2001). 

Chemical Initial and Boundary 

Conditions

At the beginning of the simulation period (initial 

condition), the aquifer was assumed to be filled with 

uncontaminated ground water of the composition 

given in table 5. All water flowing into the model 

domain through the northern boundary also had the 

composition of uncontaminated ground water. For the 

base-case calculations, only dissolved oxygen and 

nitrate were considered in the recharge from rainfall to 

the aquifer within the model domain (table 5). Finally, 

the composition of treated sewage effluent (table 5) 

was assumed to be constant throughout the operation of 

the plant, except that dissolved oxygen was assumed to 

be absent prior to the sewage plant upgrade in 1984 

(base-case only) and present from 1984 through 1995.

Chemical Reactions in the 

Sewage Plume

Determination of loads of phosphorus 

transported to Ashumet Pond was the primary objective 

of this study. However, phosphorus sorbs strongly to 

aquifer sediments, and furthermore, phosphorus 

sorption is dependent on pH. Thus, when the fate and 

transport of phosphorus are simulated, it is important to 

consider not only phosphorus sorption, but any other 

reactions in the sewage plume that affect pH. The 

reactions for the base-case model of the sewage plume 

included phosphorus sorption, cation sorption, ironand manganese-oxyhydroxide mineral equilibria, and 

decomposition of dissolved organic carbon. Additional 

reactions were included in the revised model to account 

Table 4. Estimated loading of disposal beds during the period of sewage disposal

[Each set of beds was assumed to be 10,000 square meters in area for modeling. Beds 17 and 18 were assumed never to have been used. Mgal/d, Million 

gallons per day]

Time period

Total sewageeffluent load

(Mgal/d)

Flux (meter per day)

Beds 1–4 Beds 5–8 Beds 9–12 Beds 13–16 Beds 19–25

1936–40 0.1 0.037 0.0 0.0 0.0 0.0

1941–45 1.5 .114 .114 .114 .114 .114

1946–55 .3 .023 .023 .023 .023 .023

1956–70 1 .126 .126 .126 .0 .0

1971–77 .4 .05 .05 .05 .0 .0

1978–83 .3 .114 .0 .0 .0 .0

1984–95 .2 .0 .038 .038 .0 .0
14 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

for accumulation of non-desorbable phosphorus in 

the sediments, reaction of sorbed organic carbon, 

and irreversible denitrification. The major chemical 

reactions used in the base-case and revised models 

(excluding equilibrium among aqueous species) are 

listed in table 6.

Sorption Reactions

The equilibrium constants and sorption sites 

determined from results of the sorption/desorptioncolumn experiments were used to describe phosphorus 

and cation sorption for the three-dimensional reactivetransport model. The numbers of sorption sites for 

phosphorus and cations in the aquifer were adjusted 

according to the difference in porosity between the 

experimental columns and the aquifer as follows:

, (3)

where and represent the number of 

aquifer and column sorption sites, and and 

represent the porosities of the aquifer and the column.

Iron and Manganese Minerals

Amorphous ferric hydroxide and pyrolusite 

were assumed to be present in the same amounts as 

described for the column-experiment modeling above; 

the minerals react to equilibrium unless completely 

removed by dissolution.

Organic Decomposition Reactions

The decomposition of organic carbon derived 

from the sewage effluent is an important reaction in the 

aquifer. The decomposition reactions cause large 

suboxic and anoxic zones within the plume (Walter and 

others, 1996). The most reducing conditions in the 

aquifer are represented by an anoxic zone with large 

iron concentrations between beds 5-8 and Ashumet 

Pond. A constant concentration of dissolved organic 

carbon of 1,600 μmol/L was assumed to be present in 

the sewage effluent throughout the period of sewage 

disposal. This assumed concentration of dissolved 

organic carbon in the effluent is consistent with 

measured values from the 1980s (LeBlanc, 1984). In 

the simulations, the dissolved organic carbon was 

consumed by a first-order kinetic reaction with a rate 

constant of 1.0 × 10-7 s-1. The most thermodynamically 

favored electron acceptor available—dissolved oxygen, 

nitrate, pyrolusite, or iron oxyhydroxide—was reduced 

as the organic carbon was oxidized. Results of PHAST 

simulations with this rate constant and dissolved 

organic-carbon concentration show a zone of iron 

concentrations of approximately 300 μmol/L about 

Table 5. Solution compositions for uncontaminated ground water, treated 

sewage effluent, and rainwater

[Concentrations are in micromoles per liter, unless otherwise noted. TDIC, 

total dissolved inorganic carbon; TOC, total organic carbon; °, degree; --, 

constituent not included in rainwater]

Constituent or

physical property

Uncontaminated

ground water

Treated

sewage

effluent

Rainwater

pH, standard units 5.6 6.0 5.0

Temperature, °C 14 14 25

Dissolved oxygen 125 a125 220

Calcium 29 335 --

Magnesium 31 170 --

Sodium 200 2,100 --

Potassium 10 240 --

Sulfate 86 290 --

Chloride b158 b1,361 --

Nitrate 0 1,050 10

Ammonium 0 180 --

Phosphorus 0 190 --

Iron 0 4.5 --

Manganese .64 .4 --

TOC 0 1,600 --

TDIC 28 1,200 --

aTreated sewage effluent was assumed to have no dissolved oxygen 

prior to 1984. bChloride was adjusted to achieve charge balance in the solution.

sitesaq sitescol

Pcol

1 Pcol – ------------------ ⎝ ⎠ ⎛ ⎞ 1 Paq –

Paq

----------------- ⎝ ⎠ ⎛ ⎞ =

sitescol

0.33

0.67

---------- ⎝ ⎠ ⎛ ⎞ 0.61

0.39

---------- ⎝ ⎠ ⎛ ⎞ 0.77sites = = col

sitesaq sitescol

Paq Pcol
Reactive-Transport Simulations 15

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

[C, concentration; k, rate constant; SI, saturation index; DOC, dissolved organic carbon; indicates absolute value]

Reaction Reaction type Equation

Phosphorus sorption Equilibrium

Phosphorus-sorption site protonation Equilibrium

Phosphorus-sorption site deprotonation Equilibrium

Cation sorption Equilibrium

Cation-sorption site protonation Equilibrium

Cation-sorption site deprotonation Equilibrium

Fe(OH)3(amorphous) Equilibrium

Pyrolusite Equilibrium

Decomposition of dissolved organic carbon 

(DOC), reactant is CH2O

Kinetic

Decomposition of sorbed organic carbon, 

reactant is CH2O (post cessation)

Kinetic

Vivianite precipitation and

dissolution, reactant is Fe3(PO4)2

Kinetic

Strengite precipitation and

dissolution, reactant is FePO4

Kinetic

Removal of dissolved N2 (prevents N2

reoxidation to nitrate), reactant is N2

Kinetic

SiteOH PO4

3- 2H+ + + SiteHPO4- H2 = + O

SiteOH H+ + SiteOH2

+ =

SiteOH SiteO- H+ = +

Cation_siteOH Catx+ + Cation_siteOCatx – 1 H+ = +

Cation_siteOH H+ + Cation_siteOH2

+ =

Cation_siteOH Cation_siteO- H+ = +

Fe OH ( )3( ) a 3H+ + Fe+3 3H2 = + O

MnO2 4H+ 2e

- + + Mn+2 2H2 = + O

Rate kCDOC = –

Rate kO2CO2, aqk

NO3

C

NO3, aq

= – –

Rate kCPO4

SIvivianite

1 SIvivianite + ---------------------------------------

⎝ ⎠ ⎜ ⎟ ⎛ ⎞ = –

Rate kCPO4

SIstrengite

1 SIstrengite + ---------------------------------------

⎝ ⎠ ⎜ ⎟ ⎛ ⎞ = –

Rate kCN2, aq = –

150 m downgradient of disposal beds 5-8 in 1993, 

which is consistent with the field measurements shown 

in Walter and others (1996, p. 32). However, the 

simulated iron zone extends all the way to Ashumet 

Pond, which is contrary to field observations.

After cessation of sewage disposal, the 

introduction of dissolved organic carbon ceased, 

but sorbed organic carbon present in the aquifer 

continues to react. From measurements of oxygen 

demand of contaminated sediments, a concentration of 

41,500 μmol/L of sorbed organic carbon was estimated 

for the zone between the disposal beds and Ashumet 

Pond (Richard Smith, U.S. Geological Survey, written 

commun., 2001). After cessation of sewage disposal, 

sorbed organic carbon was allowed to react with 

dissolved oxygen and nitrate in the simulation until 

these electron acceptors were consumed. First-order 

kinetics were used for both electron acceptors with 

the rate constants 1.0 × 10-7 s-1 and 1.5 × 10-7 s-1, 

respectively. (The assignment of different values was 

done inadvertently but the difference is not significant.) 

This formulation of organic-carbon decomposition 

reactions for the period following cessation of sewage 

disposal consumes oxygen and nitrate, but does not 

allow for the consumption of iron oxyhydroxide and 

the maintenance of an iron zone in the aquifer; 
16 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

however, the effect of additional organic decomposition 

that maintains an iron zone was evaluated in the 

sensitivity simulations.

Iron Phosphate Minerals

Although not included in the base-case model, 

one of the sensitivity simulations of the base-case 

model and the revised model allowed for kinetically 

controlled precipitation and dissolution of ferric 

phosphate (strengite) and ferrous phosphate (vivianite) 

minerals. The rate equations for these two minerals are 

listed in table 6.

Irreversible Denitrification

Thermodynamic equilibrium between nitrate and 

dissolved nitrogen was specified for the base-case 

simulation and all sensitivity simulations. Dissolved 

nitrogen is relatively inert, however, so that 

denitrification is simulated more realistically as an 

irreversible process. Irreversible denitrification was 

implemented in the revised model and will be 

discussed in the “Revised Model” section. 

Fitting Surface-Sorption Parameters to 

Sorption/Desorption Experimental Data

Results of the sorption/desorption column 

experiments are plotted in figure 5 with PHAST 

simulation results. Parameter estimation was used to fit 

the number of phosphorus-sorption sites and the log K 

for the phosphorus-sorption reaction (eq. 1) that were 

used in this model. The log K for surface protonation 

also was fit using parameter estimation, but the model 

was insensitive to the log K for surface deprotonation, 

which was then arbitrarily assigned a value of -7.0. 

The initial plateau in the measured pH at small 

numbers of pore volumes (~pH 5, fig. 5) could not be 

fit with phosphorus sorption alone. It was assumed that 

a set of cation-sorption reactions must account for 

the shape of the pH curve at small numbers of pore 

volumes; however, the data were insufficient to fit the 

number of sorption sites and to discriminate among the 

major cations. Only two parameters were fit for cation 

sorption: (1) the number of cation-sorption sites and 

(2) a single equilibrium constant that applied to the 

sorption reactions for all cations—Ca, Mg, Na, K, 

NH4, Fe, and Mn—of the form:

, (4)

where represents a cation-sorption site, 

and is a cation with charge .

It was difficult to fit all parameters 

simultaneously; log K values were dependent on the 

number of sorption sites, such that increasing the 

number of sites and decreasing the log K tended to 

give a similar fit. The parameters listed in table 7 

were all fit simultaneously (except for the log K for 

deprotonation), and they produced the simulation 

results for both columns that are shown in figure 5. The 

parameter fit is not unique, and depends somewhat on 

starting values for the parameter estimation. Equally 

good fits could be found with variations of the 

parameters. However, the number of cation sites from 

table 7 (23,000 μmol/L) corresponds to a cationexchange capacity of approximately 0.5 meq/100 g 

sediment, which is consistent with estimates of cationexchange capacity for these sediments (Ceazan and 

others, 1989). With the number of cation sites fixed at 

23,000 μmol/L, an optimum fit to the data could be 

found only if the number of phosphate sorption sites 

was within the range of 2,000 to 3,500 μmol/L. For this 

range in number of phosphate sorption sites, the 

corresponding variations in the fit of log K values for 

phosphate sorption and site protonation were both 

approximately 1.2 log units. The fit values for the 

phosphate parameters were not sensitive to the cation 

parameters. 

Sensitivity Simulations

The sensitivity of the three-dimensional reactivetransport model to selected parameters was investigated 

by individually perturbing model parameters or sets of 

model parameters that described a single process. The 

effects of 11 parameters and parameter sets were 

simulated, and for each perturbation, the load of 

phosphorus transported to Ashumet Pond was 

calculated for each year of the simulation. Only the 

relative effect of the parameter change on the load 

transported to Ashumet Pond should be considered 

because of the preliminary nature of the base-case set 

of model parameters.

Cation_siteOH Catx+ +

= Cation_siteOCatx – 1 H+ +

Cation_site

Catx+ x
17 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

0 50 100 150 200 250 300

NUMBER OF PORE VOLUMES

0

100

200

300

400

PHOSPHORUS, IN MICROMOLES

PER LITER

COLUMN 1, MEASURED

COLUMN 1, SIMULATED

COLUMN 3, MEASURED

COLUMN 3, SIMULATED

4

5

6

7

8

pH

ELUENT CHANGED TO UNCONTAMINATED GROUND WATER

ELUENT CHANGED TO UNCONTAMINATED GROUND WATER

ELUENT CHANGED TO UNCONTAMINATED GROUND WATER

ELUENT CHANGED TO UNCONTAMINATED GROUND WATER

Figure 5. Results of column experiments of phosphorus sorption and desorption on Cape Cod sediments and results of reactive-transport 

simulations based on fitted parameters. Influent phosphorus concentration in column 1 was 16 micromoles per liter and in column 3 was 196 

micromoles per liter.

Usually, a parameter value was increased and 

decreased to evaluate its effect; however, this pattern 

was not followed in the following four cases: (1) The 

sensitivity of the model to the maintenance of an iron 

zone following cessation of sewage disposal was 

investigated by adding a kinetic reaction that allowed 

sufficient sorbed organic carbon to react to generate 

approximately 150 μmol/L dissolved iron. This process 

was specified to occur for the 60 years following the 

cessation of sewage disposal in a 900-m × 300-m ×

25-m-thick zone that contained the disposal beds and 

extended to Ashumet Pond. (2) The sensitivity of the 

model to the fixed head of the southern leaky boundary 

was investigated by lowering the head, but not by 

raising the head. (3) The effect of numerical dispersion 

was investigated by halving the time step and X- and Ydirection grid size (Z direction unchanged), which 

should halve the horizontal numerical dispersivity to 

approximately 30 m. (4) The sensitivity of the model to 

sequestration of phosphorus in a non-desorbable form 

was investigated by including equilibrium and kinetic 

mineral reactions with iron-phosphate minerals. 

Iron-phosphate colloids have been identified 

near the disposal beds (Gschwend and Reynolds, 

1987), and phosphate sorption may ultimately lead 

to surface precipitation of solids or solid solutions 

(Dzombak and Morel, 1990). To examine the effects 

of incorporating phosphorus into the sediments in a 

non-sorbing form, strengite (ferric phosphate) and 

vivianite (ferrous phosphate) were included in the 
18 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

Table 7. Parameters derived from fitting sorption/desorption column experiments

[Site, phosphorus-sorption site. Cation-site, cation-sorption site. μmol/L, micromoles per liter]

Parameter Value

log K: (eq. 1) 26.7

log K: (eq. 4) -1.8

Sorption sites for PO4

3-, μmol/L 3,000

Sorption sites for cations, μmol/L 23,000

log K: and 4.1

log K: and -7.0

SiteOH PO4

3- 2H

+

+ + SiteHPO4

- H2 = + O

Cation_siteOH Catx++ Cation_siteOCatx – 1 H

+ = +

SiteOH H+ + SiteOH2

+ = Cation_siteOH H

+

+ Cation_siteOH2

+ =

SiteOH SiteO- H

+ = + Cation_siteOH Cation_siteO- H+ = +

model as equilibrium phases (log K values defined by 

wateq4f.dat, Parkhurst and Appelo, 1999). These 

phases were not present at the beginning of the 

simulation, but could form when the solution became 

saturated and subsequently continue to precipitate 

or redissolve to maintain equilibrium. A kinetic 

formation of these minerals was implemented with 

a rate expression that was formed as the product 

of the following: a rate constant (2.3 × 10-8 s-1), the 

concentration of dissolved phosphorus, and a Monod 

factor based on the saturation index of the mineral 

with half-saturation constant of 1.0 (unitless). This rate 

constant produced a distribution of strengite at the time 

of cessation that was greatest at the disposal beds and 

decreased downgradient to near zero in about 400 m. 

The same rate constant was used for vivianite. 

The base-case and perturbed values for 

parameters are listed in table 8 for 11 sensitivitysimulation cases. Figure 6 shows the phosphorus load 

transported to Ashumet Pond for simulations with 

perturbed parameters compared to the load calculated 

for the simulation with the base-case parameters. For 

most simulations the load transported to Ashumet Pond 

has a peak around 1970 and second peak around 2005.

The timing of the first peak is generally 

consistent among the calculations, but is delayed 

slightly by decreasing the cation-sorption log K 

(case 3) or increasing the number of cation-sorption 

sites (case 4). Decreasing the cation-sorption log K 

relative to the base case produces initial conditions 

with nearly twice as many uncharged OH sites that can 

participate in deprotonation reactions to maintain the 

pH at lower values for a longer period of time than the 

base case. The lower pH values cause stronger sorption 

and smaller dissolved phosphorus concentrations; 

similarly, increasing the number of cation-sorption 

sites allows a longer period of low pH initially, during 

which phosphorus is sorbed more strongly. 

The magnitude of the first peak decreases or 

increases with adjustment to most parameters. The 

following perturbations tend to decrease the magnitude 

of the first peak in phosphorus load transported to 

Ashumet Pond: increasing the phosphorus-sorption 

constant (case 1), increasing the number of 

phosphorus-sorption sites (case 2), decreasing the 

concentration of dissolved organic carbon in sewage 

effluent (case 5), decreasing the rate constant for 

decomposition of dissolved organic carbon (case 6), 

lowering the specified head for the southern leaky 

boundary (case 8), decreasing the phosphorus 

concentration in sewage effluent (case 9), and 

formation of iron-phosphate minerals (case 10). 

The magnitude of the first peak appears to be 

most sensitive to the concentration of dissolved organic 

carbon in the sewage effluent, the concentration of 

phosphorus in the sewage effluent, and the formation of 

iron-phosphate minerals. Doubling the concentration 

of dissolved organic carbon in the sewage effluent 

(case 5) causes intense iron dissolution and sulfate 

reduction in the simulation, which has not been 

observed in the plume. These redox reactions generate 

pH values as high as 9.5, which are also inconsistent 

with field measurements. The high pH values have the 

effect of diminishing the sorption of phosphorus and 

increasing the load of phosphorus transported to 

Ashumet Pond. Halving the concentration of dissolved 

organic carbon (case 5) causes the iron zone not to 

form in the simulation, and causes pH values of about 5 
Reactive-Transport Simulations 19

Table 8. Base-case and perturbed values for sensitivity simulations with the three-dimensional reactive-transport model

[X, X-direction (fig. 4); Y, Y-direction (fig. 4); Z, Z-direction (vertical); P, phosphorus; m, meter; μmol/L, micromoles per liter; s-1, per second; --, no second 

sensitivity simulation for this case]

Case Parameter or Process Base

case Value 1 Value 2

1 P sorption constant (log K, eq. 1) 26.7 27.2 26.2

2 P sorption sites, μmol/L a2,300 4,600 1,150

3 Cation-sorption constant (log K, eq. 4) -1.8 -0.8 -2.8

4 Cation-sorption sites, μmol/L a18,000 36,000 9,000

5 Dissolved organic carbon in sewage effluent, μmol/L 1,600 3,200 800

6 First-order rate constant for dissolved organic carbon, s-1 1.0 × 10-7 2.0 × 10-7 0.5 × 10-7

7 Generation of iron zone following cessation of sewage disposal no1 yes2 --

8 Specified head for southern leaky boundary, m 12.19 11.19 10.19

9 P concentration in sewage effluent, μmol/L 190 285 95

10 Non-desorbable phosphorus none Strengite and vivianite 

equilibrium

Kinetic formation of 

strengite and vivianite3

11 Numerical dispersivity, m ~60 ~30 --

Node spacing X, m 50 25 --

Node spacing Y, m 50 25 --

Node spacing Z, m 2.5 2.5 --

Time step, year 0.5 0.25 --

a Number of sites is adjusted from values in table 7 to account for the difference in porosity.

1 Sufficient organic carbon reacted to consume oxygen and nitrate in a 900 × 300 × 25-m zone containing the disposal beds and extending to Ashumet 

Pond.

2 200 micromoles of organic carbon per liter of water reacted kinetically to consume oxygen and nitrate, and generate approximately 150 μmol/L 

dissolved iron. This reaction was applied in a 900 × 300 × 25-m zone containing the disposal beds and extending to Ashumet Pond. 

3 Rate equations for strengite and vivianite were , where k = 2.3 × 10-8, is the total dissolved phosphorus 

concentration, and SI is the saturation index for the mineral.

Rate kCPO4

SI

1 + SI ----------------- ⎝ ⎠ ⎛ ⎞ = – CPO4

within the plume. These low pH values also are 

inconsistent with field measurements, but have the 

effect of increasing the sorption of phosphorus and 

decreasing the load of phosphorus transported to 

Ashumet Pond. The concentration of phosphorus in 

the sewage effluent has a direct effect on the loads of 

phosphorus transported to Ashumet Pond; increased 

phosphorus concentration leads directly to increased 

loads transported to Ashumet Pond and the reverse 

for decreased phosphorus concentration (case 9). 

Equilibrium with iron-phosphate minerals essentially 

sequesters all phosphorus in mineral phases within a 

short distance of the disposal beds (case 10); in this 

case, simulated small concentrations of phosphorus 

within the aquifer are inconsistent with large measured 

concentrations of phosphorus downgradient of the 

disposal beds. With the specified kinetic formation of 

iron-phosphate minerals, ferric phosphate accumulates 

near the disposal beds and ferrous phosphate 

accumulates in the large-iron-concentration zone, 

but phosphorus still discharges to Ashumet Pond. The 

load of phosphorus transported to Ashumet Pond is 

then sensitive to the kinetic-rate formulation—the 

simulated phosphorus load may vary from near zero 

for fast rates of iron-phosphate mineral formation 

(producing near-equilibrium conditions), up to the 

base-case results for slow rates (no mineral formation). 

Kinetic formation of iron-phosphate minerals has a 

greater effect on the first peak in phosphorus load than 

on the second peak.

The timing of the second peak in phosphorus 

load transported to Ashumet Pond does not appear 

to be affected by any of the parameter perturbations 

and occurs around 2005. The magnitude of the second 

peak appears to be most sensitive to four parameters 

or processes: (1) Increasing the number of phosphorussorption sites (case 2) tends to decrease the second 

peak in phosphorus load transported to Ashumet 

Pond. The number of sorption sites for phosphorus 

is uncertain by at least a factor of two, and additional 
20 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

1950 2000 2050

0

2,000

4,000

6,000

8,000

BASE CASE

SEWAGE DOC X 2

SEWAGE DOC X 0.5

5

0

2,000

4,000

6,000

8,000

PHOSPHORUS LOAD TO ASHUMET POND, IN KILOGRAMS PER YEAR

BASE CASE

CATION SORPTION LOG(K) + 1

CATION SORPTION LOG(K) - 1 

3

0

2,000

4,000

6,000

8,000

BASE CASE

P SORPTION LOG(K) + 0.5

P SORPTION LOG(K) - 0.5

1

1950 2000 2050

BASE CASE

DOC RATE CONSTANT X 2

DOC RATE CONSTANT X 0.5

6

BASE CASE

CATION SORPTION SITES X 2

CATION SORPTION SITES X 0.5

4

BASE CASE

P SORPTION SITES X 2

P SORPTION SITES X 0.5

2

YEAR

Figure 6. Sensitivity of phosphorus load transported to Ashumet Pond to variation in selected reactive-transport model parameters.

simulations in the “Revised Model” section 

indicate that the number of sites could be larger. 

(2) Maintenance of the iron zone following cessation 

of sewage disposal (case 7) tends to increase the second 

peak in phosphorus load transported to Ashumet 

Pond. Measurements of oxygen demand for the 

contaminated sediments indicate that hundreds of pore 

volumes of water are needed to consume the sediment 

oxygen demand. However, the extent to which iron 

oxyhydroxides will be reduced in the process of 

oxidizing the organic carbon (the cause of the oxygen 

demand) is unknown. (3) Allowing more water to flow 

out across the southern boundary (decreasing the 

specified head for the leaky boundary, case 8) causes

a decrease in the second peak in phosphorus load 

transported to Ashumet Pond. The relative flow of 

water out across the southern boundary compared to 

flow of water to Ashumet Pond is not well constrained, 

but additional simulations in the “Revised Model” 

section indicate that the plume in the Ashumet Valley 

may be better represented if more water is allowed 

to flow out across the southern leaky boundary. 

(4) Rapid formation of iron-phosphate minerals 

(equilibrium, case 10) nearly eliminates phosphorus 

loads transported to Ashumet Pond. 

In addition to the physical and chemical 

processes, the simulation results include the effects of 

numerical dispersion. Ideally, the grid and time step of 
Reactive-Transport Simulations 21

0

2,000

4,000

6,000

8,000

PHOSPHORUS LOAD TO ASHUMET POND, IN KILOGRAMS PER YEAR

BASE CASE

P CONCENTRATION x 1.5

P CONCENTRATION x 0.5

9

0

2,000

4,000

6,000

8,000

BASE CASE

IRON ZONE AFTER CESSATION

7

1950 2000 2050

YEAR

0

2,000

4,000

6,000

8,000

BASE CASE, DISPERSIVITY ~60 m 

REFINED GRID/TIME STEP,

DISPERSIVITY ~30 m 

11

1950 2000 2050

YEAR

BASE CASE

EQUILIBRIUM STRENGITE AND VIVIANITE

KINETIC STRENGITE AND VIVIANITE

10

BASE CASE

SOUTH LEAKY BOUNDARY -1 m

SOUTH LEAKY BOUNDARY -2 m

8

Figure 6. Sensitivity of phosphorus load transported to Ashumet Pond to variation in selected reactive-transport model parameters—Continued.

the model are refined until the effects of numerical 

dispersion are less than the physical dispersion that is 

simulated. Unfortunately, limitations in computer 

memory and computer processing time preclude 

refinement of the model grid and time step to the point 

of grid and time-step convergence, that is, to the point 

that no change in results occurs with further refinement 

of the grid and time step. The base-case model was run 

with a refined grid and decreased time step (table 8). 

The dispersivity for this refined grid and time step is 

estimated to be 30 m, compared to 60 m for the original 

mesh and time step. The results for the refined grid are 

similar to the results for the coarse grid (fig. 6, case 

11), but with somewhat increased loads for the last two 

peaks. Further refinement of the grid could change the 

peak loads substantially; however, the pattern of 

decreasing loads after about 2005 does not appear to be 

sensitive to numerical dispersion.

Revised Model

In this section, several of the base-case 

parameters are adjusted to develop a revised model that 

better simulates the available data. The results of the 

sensitivity simulations for the base-case model are used 

to provide qualitative interpretations for the changes in 

phosphorus load transported to Ashumet Pond that are 

simulated by the revised model.
22 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

Modifications to the Base-Case Model

By considering additional field and experimental 

data collected in 1996 and 1999, and reconsideration of 

chemical processes, the revised model was developed 

to account for additional processes that are evident in 

the aquifer, to refine estimates of model parameters, 

and to finalize estimates of the load of phosphorus 

transported to Ashumet Pond. Four processes that were 

not present in the base-case modeling are implemented 

in the revised model: (1) organic-carbon and oxygendemand accumulation on the aquifer sediments, 

(2) irreversible denitrification of nitrate, (3) nondesorbable phosphorus accumulation on the aquifer 

sediments, and (4) rapid reaction of aquifer sediments 

with rainwater that infiltrates the aquifer. This section 

also provides justification for adjusting some of the 

model parameters used in the base-case modeling, 

specifically, the head assigned to the southern model 

boundary, the number of phosphorus-sorption sites, and 

the oxygen concentrations in the sewage effluent 

during the first 50 years of disposal. Revisions to the 

base-case model are summarized in table 9.

Sorption and Subsequent Reaction of 

Organic Carbon

Over the period of sewage disposal, a large 

amount of organic carbon accumulated on the aquifer 

sediments. Measurements of total oxygen demand 

range from 2.1 to 12.1 μmol of O2 per g dry sediment 

in the area downgradient of the disposal beds (Richard 

Smith, U.S. Geological Survey, written commun., 

2001), with a general trend of decreasing oxygen 

demand with distance downgradient. Assuming 4,150 g 

of sediment per liter of water, each volume of aquifer 

containing one liter of water has from 9,000 to 

50,000 μmol of oxygen demand and requires 35 to 

200 pore volumes of air-saturated water 

(approximately 250 μmol/L of O2) to satisfy this 

oxygen demand. Thus, with flow time from the 

Table 9. Comparison of parameters and processes for base-case and revised models

[P, phosphorus; m, meter; μmol/L, micromoles per liter; s-1, per second]

Parameter or process Base-case

model

Revised

model

Sorb organic carbon during sewage disposal, react sorbed organic carbon following cessation of sewage 

disposal

no yes

Irreversible denitrification no yes

Kinetic formation of non-desorbable phosphorus (strengite and vivianite, table 6) no yes

Infiltration at top of aquifer rain background

water

Generation of iron zone following cessation of sewage disposal no no

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

P sorption sites, μmol/L 2,300 3,450

Cation-sorption constant (log K, eq. 4) -1.8 -1.8

Cation-sorption sites, μmol/L 23,000 23,000

Dissolved oxygen concentration in sewage effluent, 1936–83, μmol/L 0 250

First-order rate constant for conversion of dissolved nitrogen to non-reactive dissolved species, s-1 0 1.16 × 10-7

P concentration in sewage effluent, 1935–83, μmol/L 190 380

Rate constant for conversion of dissolved phosphorus to non-desorbable phosphorus (strengite and 

vivianite rate expressions, table 6), s-1

0 1.16 × 10-8

Dissolved organic carbon in sewage effluent, μmol/L 1,600 1,600

First-order rate constant for dissolved organic carbon, s-1 1.0 × 10-7 1.0 × 10-7

Sorbable organic carbon in sewage effluent, μmol/L 0 1000

First-order rate constant for sorbable organic carbon, 1936–95, s-1 0 1.16 × 10-8

First-order rate constant for consumption of electron acceptors following cessation of sewage disposal, 

1996–2055. Reaction consumes sorbed organic carbon, s-1

0 1.16 × 10-8

Specified head for southern leaky boundary, m 12.19 10.19
Reactive-Transport Simulations 23

disposal beds to Ashumet Pond on the order of 4 years, 

anoxic conditions are likely to persist downgradient 

from the disposal beds for many decades following 

cessation of sewage disposal.

To simulate the process of accumulation of 

oxygen demand within the aquifer, an additional pool 

of organic carbon was added to the sewage effluent, 

which represented that part of the organic carbon that 

sorbed on the sediment during disposal of sewage 

effluent from 1936 through 1995. The concentration 

of sorbable organic carbon was specified to be 

1,000 μmol/L in the sewage effluent, which sorbed by 

first-order reaction with rate constant 1.16 × 10-8 per 

second. This concentration and rate constant produced 

maximum sorbed organic carbon concentrations of 

16,000 μmol/L after 60 years of sewage-effluent 

disposal; the sediment oxygen demand corresponding 

to this accumulation of sorbed organic carbon is also 

16,000 μmol/L.

After cessation of sewage disposal, the model 

includes kinetic reactions that decrease the sorbed 

organic carbon concentration as the oxygen in the 

invading background water is consumed. For each 

mole of oxygen consumed, one mole of sorbed 

organic carbon is oxidized; the rate of consumption 

is first order with respect to the dissolved oxygen 

concentration. The presence of sorbed organic carbon 

is required for the process to continue.

Irreversible Denitrification

Thermodynamic equilibrium between nitrate 

and dissolved nitrogen was specified for the base-case 

modeling. Field data indicate that dissolved nitrogen 

is formed by denitrification within the plume as 

organic carbon is consumed; further reduction of 

nitrate to ammonium is not apparent. However, in the 

simulations, dissolved nitrogen within the plume mixed 

with dissolved oxygen outside the plume by the process 

of dispersion, and dissolved nitrogen reoxidized to 

nitrate, with a concomitant decrease in pH. Because 

dissolved nitrogen is relatively inert and unlikely to 

reoxidize, denitrification is simulated more realistically 

as an irreversible process. Irreversible denitrification 

was implemented in the revised model by removing 

dissolved nitrogen from the system as it formed from 

denitrification and adding an equal amount of a new 

dissolved nitrogen species that could not reoxidize to 

nitrate in the presence of oxygen. The removal of 

dissolved nitrogen was modeled as a first-order decay 

with rate constant 1.16 × 10-7 s-1; the new non-reactive 

dissolved nitrogen species was produced at an equal 

rate. The rate constant for these reactions was sufficient 

to remove all of the dissolved nitrogen essentially as 

soon as it was formed.

Non-Desorbable Phosphorus

Experiments with contaminated sediments from 

the aquifer show that not all of the phosphorus can be 

removed by desorption. Even after repeated rinsing 

with pH-10 water, most of the phosphorus is retained 

by the sediment. Total non-desorbable phosphorus 

ranges from about 0.5 to 8 μmol/g of sediment, or 

about 2,000 to 35,000 μmol/L of water in the aquifer. 

Formation of a large amount of non-desorbable 

phosphorus in the aquifer requires a greater input of 

phosphorus into the system, and accordingly, the 

concentration of phosphorus in the sewage effluent 

for the period 1936 through 1983 for the revised 

model was increased to 380 μmol/L from the basecase concentration of 190 μmol/L. This increase in 

concentration is justified, in part, by the large pool of 

sorbed phosphorus in the sediments, and, in part, by 

the use of high-phosphate detergents within this time 

period. 

The revised model implemented the kinetic 

reactions for ferric phosphate (strengite) and ferrous 

phosphate (vivianite) minerals used in the sensitivity 

simulations (case 10, fig. 6). The combination of the 

increased concentration of phosphorus in the 

sewage effluent through 1983 and the kinetic 

reactions produced concentrations of strengite up to 

7,000 μmol/L and vivianite up to 1,600 μmol/L. 

Simulation results indicate that most of the strengite 

forms immediately downgradient of the disposal beds, 

whereas vivianite tends to form in and downgradient 

from the iron zone.

Infiltration at the Surface of the Aquifer

The base-case model recharged dilute nitric acid 

(pH 5) into the aquifer to simulate the infiltration of 

rainwater (table 5). However, the reactive surfaces in 

the model were initialized to be in equilibrium with the 

uncontaminated ground-water composition (table 5). 

Reaction of the dilute infiltration water with the preequilibrated surfaces caused a pH higher than 6 in the 

shallow part of the aquifer throughout the simulation. 

The pH measured in the shallow aquifer is 
24 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

approximately 5.5, which is similar to the pH of the 

background water. Thus, for the revised model, the 

background-water composition was used for recharge 

in place of rainwater. It is assumed that rainwater reacts 

rapidly in the shallow part of the aquifer to achieve a 

composition similar to the uncontaminated background 

water (table 5).

Leaky Boundary Condition for the 

Southern Boundary

By using the specified concentration of chloride 

in the sewage effluent (1,361 μmol/L, table 5) and the 

prescribed loading history, the total mass of chloride 

applied through the disposal beds is 2,200,000 kg for 

all simulations. All of this chloride must eventually 

leave the model domain, either through the constanthead nodes that represent Ashumet Pond or through the 

leaky boundary at the southern edge of the domain. In 

the base-case model, 90 percent of the sewage-derived 

chloride left the domain through nodes that represent 

Ashumet Pond, and about 10 percent, or 200,000 kg, 

left the domain through the southern boundary. This 

small flux of chloride through the southern boundary 

appears to be inconsistent with the development of the 

large sewage plume in Ashumet Valley (fig. 1).

To estimate the mass of chloride contributed by 

sewage effluent to the Ashumet Valley plume, a mass 

accounting for chloride in Ashumet Valley was done. A 

total of 328 chloride-concentration data from 1994 

from sites throughout Ashumet Valley (Savoie and 

LeBlanc, 1998) were assembled. These data were 

interpolated to a regular, three-dimensional grid, from 

which the mass of chloride in the Valley was estimated 

by means of a zero moment routine (R.L. Naff, U.S. 

Geological Survey, written commun., 2002). Any 

excess chloride above the background concentration of 

10 mg/L was assumed to be contributed by the sewage 

effluent. The total mass of chloride in the valley 

contributed by the sewage effluent is estimated to be 

2,000,000 to 2,500,000 kg. This range of estimate still 

has a large uncertainty because of the sparsity of spatial 

coverage of the data, the selection of a constant 

background concentration, and because of the difficulty 

of accurately interpolating and integrating the data over 

a large region. However, the base-case calculation of 

only 200,000 kg of chloride entering Ashumet Valley 

appears to be an order of magnitude too small.

The parameter in the model that effectively 

distributes the flow between Ashumet Pond and the 

southern boundary is the specified head for the leaky 

southern boundary condition. By lowering the head 

for the leaky southern boundary condition, more 

water will flow through the southern boundary and less 

water will flow to Ashumet Pond. By lowering the 

head associated with the boundary by two meters (to 

10.19 m), the mass of chloride transported through the 

southern boundary is about 1,100,000 kg. Although 

still smaller than the estimate of sewage-derived 

chloride derived from the measured concentration data 

in the Ashumet Valley, it is within a factor of two of the 

estimate, and it accounts for nearly half of the sewagederived load of chloride in the simulation. Further 

lowering the head of the boundary condition generates 

head gradients that are significantly greater than are 

currently measured in the flow system. The revised 

model uses a 2-m reduction in the head associated with 

the leaky southern boundary, which is a compromise 

between allowing sufficient sewage-derived chloride 

to flow into Ashumet Valley and maintaining head 

gradients that approximate water-level measurements.

Number of Sorption Sites

Model parameters for sorption processes were 

fit from sorption/desorption column experiments on 

uncontaminated sediments. A column experiment 

desorbing phosphorus from contaminated sediments 

also was run (see “Methods” section). This desorption 

column experiment was simulated by using the 

parameters derived from the sorption/desorption 

column experiments. 

The main difficulty in simulating the phosphorus 

desorption experiment was in the definition of initial 

conditions. The core material was obtained from site 

S469 at the southern edge of the disposal beds in April 

2000, approximately 4.5 years after cessation of 

sewage disposal. To account for the period of time 

between cessation of sewage disposal and sample 

collection, the initial composition of the reactive 

surfaces in the column were modified as follows: first, 

the surfaces were brought to equilibrium with treated 

sewage effluent (table 5), and second, the surfaces were 

flushed with 15 pore volumes of uncontaminated 

ground water (table 5). The simulated phosphorus 

concentrations in equilibrium with the surface 
Reactive-Transport Simulations 25

compositions following this procedure were similar to 

the initial phosphorus concentrations released during 

the phosphorus desorption experiment (fig. 7). 

The simulation results for the desorption column 

experiment generated with base-case parameters for 

sorption reactions are shown in figure 7. The general 

trend in phosphorus concentration in the leachate is 

captured by the simulation, although simulated 

phosphorus concentrations are consistently smaller 

than the experimental data. It is plausible that the 

number of sorption sites may vary within the aquifer 

and may vary as a result of reactions in the sewage 

plume, such as sorption of organic carbon or 

precipitation of additional iron oxyhydroxide surfaces. 

Surface areas measured on aquifer sediments varied by 

a factor of two within the same sediment core (Fuller 

and others, 1996), which implies that the number of 

sorption sites could also vary by at least this extent 

within the aquifer. Increasing the number of surface 

sites in the simulation of the desorption experiment 

tends to increase the phosphorus concentrations in the 

leachate because a larger pool of desorbable 

phosphorus is defined by the initial conditions. A 

factor-of-two increase in the number of sites, while 

maintaining the same log K for the phosphorussorption reaction, tends to produce simulated 

concentrations similar to measured concentrations for 

25 to 150 pore volumes, although simulated 

phosphorus concentrations are larger than measured 

concentrations for the first 25 pore volumes and 

smaller than measured concentrations for the last 50 

pore volumes (fig. 7). 

The number of sites for the revised model 

was increased by a factor of 1.5, splitting the 

difference between the number of sites fit from the 

sorption/desorption experiments and the apparent need 

for a greater number of sites in the desorption 

experiment. 

0 50 100 150 200

NUMBER OF PORE VOLUMES

0

10

20

30

40

PHOSPHORUS, IN MICROMOLES PER LITER

MEASURED

SIMULATED (BASE-CASE MODEL)

SIMULATED (PHOSPHORUS SITES X 2.0)

ELUENT pH CHANGED TO 7.42

ELUENT pH CHANGED TO 5.78

ELUENT pH INITIALLY 5.78

Figure 7. Measured and simulated phosphorus concentrations for a column desorption experiment on a 

sediment core from the sewage-contaminated zone of the aquifer.
26 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

Oxygen Concentration in 

Sewage Effluent

The oxygen concentration in the sewage 

effluent for the base-case model was arbitrarily set to 

zero for the period 1936–83. For the revised model, a 

constant oxygen concentration in the sewage effluent 

(250 μmol/L) was used for the entire period of sewage 

effluent disposal. The oxygen is consumed in the 

oxidation of organic carbon, which produces carbon 

dioxide. The additional input of carbon dioxide tends to 

buffer the pH at lower values for a longer period of 

time. The overall effect of the additional oxygen on the 

simulations is relatively small.

Comparison of Measured and 

Simulated Phosphorus Concentrations

Reactive-transport simulations provide 

concentrations of phosphorus at each active node 

in the model grid for each year of the period of 

simulation. Simulation results from the revised model 

are compared to measured phosphorus concentrations 

in space and time. 

Phosphorus Distribution in 1993

As a test of the revised model, simulated 

phosphorus concentrations were compared with a 

comprehensive set of measured phosphorus 

concentrations in the aquifer from August through 

November 1993 (Walter and others, 1996). Figure 8 

shows contours of the measured phosphorus 

concentration from 1993; maximum concentrations 

were used for sites that had samples from multiple 

depths. Figure 8 also shows contours of simulated 

phosphorus concentration for the end of 1993; for 

consistency, maximum concentration was used for 

each vertical set of nodes. The simulated phosphorus 

concentrations are similar to the measured phosphorus 

concentrations, both in location and in magnitude 

of the peak concentrations downgradient from the 

disposal beds. In the simulation, the easternmost large 

concentrations are derived from beds 1-4, which were 

the only beds used during the period 1978–83; those 

beds were not used thereafter. The westernmost large 

concentrations in the simulation are derived from beds 

5–12, which were the only beds used during 

the period 1984–95. The location of the observed 

phosphorus concentration peaks and the magnitude 

of the measured phosphorus concentrations are 

reproduced by the reactive-transport simulations, but 

the maximum simulated phosphorus concentrations 

southeast of beds 1–4 are larger (270 μmol/L) than 

maximum measured concentrations in this area 

(180 μmol/L). The simulated large phosphorus 

concentrations southeast of beds 1-4 extend all 

the way to Ashumet Pond, but the measured large 

concentrations decrease toward the Pond. The 

measured concentrations show smaller, narrower lobes 

of large phosphorus concentrations, which may be due 

to selective use of individual disposal beds (compared 

to the equal distribution of sewage effluent to groups of 

disposal beds simulated in the model), to the small 

dispersivity of the aquifer, or to inadequate spatial 

coverage of the ground-water sampling locations.

Phosphorus at Site S469

Decreases in phosphorus concentrations have 

been measured in wells near the disposal beds since 

the cessation of sewage disposal. Site S469 is 

at the edge of disposal beds 5-8, and a multilevel 

sampler at this location was sampled in July 1996 and 

July 1999. During that period, phosphorus 

concentrations decreased significantly 

at altitudes greater than -10 m. Measured phosphorus 

concentrations in multilevel sampler S469M1 and 

phosphorus concentrations simulated with the revised 

model for the vertical set of nodes nearest to the 

location of disposal beds 5-8 are shown in figure 9. The 

simulated concentrations are consistent with 

a marked decrease in measured phosphorus 

concentrations near the disposal beds in the time frame 

of 3 years. The magnitude of the simulated 

concentration change also is consistent with the 

magnitude of the measured concentration change.
Reactive-Transport Simulations 27 MODELED DISPOSAL BEDS AND IDENTIFIERS SEWAGE DISPOSAL BEDS MEASURED LINES OF EQUAL MAXIMUM DISSOLVED PHOSPHORUS CONCENTRATION—Dashed where approximately located. Interval, in micromoles per liter, is variable. MODELED DISPOSAL BEDS AND IDENTIFIERS SEWAGE DISPOSAL BEDS SIMULATED LINES OF EQUAL MAXIMUM DISSOLVED PHOSPHORUS CONCENTRATION— Interval, in micromoles per liter, is variable.

Ashumet

Pond

0.1

EXPLANATION

Ashumet

Pond

0.1

EXPLANATION

B. A.

60

30

150

120

90

180

240210 270

3

3

30

120

150

120

90 150180 30

3

60

90

70o33'00"

41o38'30"

70o32'00"

41o37'30"

70o33'00" 

41o38'30"

70o32'00" 

41o37'30"

00 500

1,500

1,000 METERS

3,000 FEET 0

0 500

1,500

1,000 METERS

3,000 FEET

Modified from Walter and others (1996)

Figure 8. (A) Measured and (B) simulated phosphorus concentrations in ground water from August to November 1993 near Ashumet Pond, Massachusetts.
28 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

0 50 100 150 200 250

PHOSPHORUS, IN MICROMOLES PER LITER

-30

-20

-10

0

10

20

ALTITUDE, IN METERS ABOVE NGVD29

1996 MEASURED

1996 SIMULATED

1999 MEASURED

1999 SIMULATED

MULTILEVEL SAMPLER S469M1

Figure 9. Measured phosphorus concentrations in July 1996 and July 1999 in multilevel sampler S469M1 and

simulated phosphorus concentrations for the vertical set of nodes nearest to the location of disposal beds 5–8.

Transport of Phosphorus to 

Ashumet Pond

The timing and sequence of the peaks in 

phosphorus load transported to Ashumet Pond are 

determined by a combination of effects related to 

phosphorus concentration in the sewage effluent, 

the sequence of loading the disposal beds, and the 

chemical response of the aquifer to the loading of 

sewage effluent and related chemical reactions. 

However, several results are consistent for all plausible 

simulations with parameters from the base-case, 

sensitivity, and revised-model simulations. First, the 

earliest substantial phosphorus loads transported to 

Ashumet Pond are simulated to have occurred in the 

1960s, nearly 30 years after the initial use of the 

disposal beds and 20 years after peak loading of the 

disposal beds. Second, the simulated phosphorus load 

transported to Ashumet Pond decreases with time after 

the 2000–10 time period; however, substantial 

phosphorus loads continue for decades. And third, the 

phosphorus load transported to Ashumet Pond is a 

small fraction of the total phosphorus loaded through 

the disposal beds. Most of the phosphorus is retained in 

the aquifer sediment in the sorbed and non-desorbable 

fractions. Consistent with the sorption/desorption and 

desorption column experiments, the simulations 

indicate that this large pool of sediment-bound 

phosphorus remains sequestered for a long period of 

time and ultimately maintains small concentrations of 

dissolved phosphorus in the ground water.

Loads of phosphorus at the disposal beds 

and the loads of phosphorus transported to Ashumet 

Pond simulated by using the revised model are 

plotted in figure 10. For the revised model, the 

simulated load of phosphorus transported to Ashumet 

Pond increased in the late 1960s and remained 

approximately constant until 1985. The constant load 

over this period of time is markedly smaller than the 

corresponding peak in load calculated for the base-case 
Reactive-Transport Simulations 29

1940 1960 1980 2000 2020 2040 2060

YEAR

0

10,000

20,000

30,000

40,000

50,000

60,000

PHOSPHORUS LOAD TO DISPOSAL BEDS,

IN KILOGRAMS PER YEAR

DISPOSAL-BED LOAD

0

500

1,000

1,500

2,000

2,500

PHOSPHORUS LOAD TO ASHUMET POND,

IN KILOGRAMS PER YEAR

ASHUMET POND LOAD, REVISED MODEL

ASHUMET POND LOAD, REVISED MODEL

WITHOUT VIVIANITE FORMATION

Figure 10. Simulated load of phosphorus transported to Ashumet Pond, Massachusetts, and estimated load of phosphorus 

to disposal beds during the period 1936 to 2055, calculated with revised model parameters and revised model parameters 

excluding vivianite formation.

model for this period of time (fig. 6), even though a 

larger concentration of phosphorus was introduced 

with the sewage effluent in the revised model over most 

of the period of sewage disposal. The larger input in 

phosphorus simulated in the revised model does not 

cause larger loads to be transported to Ashumet Pond 

compared to the base-case model because (1) the 

amount of phosphorus transported away from the 

Pond toward Ashumet Valley is greater because of the 

lower head for the southern leaky boundary, (2) the 

uptake of phosphorus by sorption is greater because 

of the increased number of phosphorus-sorption sites, 

and (3) the transfer of phosphorus to non-desorbable 

phosphorus sequesters phosphorus in the sediment. 

The load of phosphorus transported to Ashumet 

Pond is especially sensitive to the formation of 

vivianite. A simulation using the revised model, but 

without vivianite formation, shows larger loads of 

phosphorus transported to Ashumet Pond during the 

period 1965–85 (fig. 10) compared to smaller loads in 

the revised model with vivianite formation. In the 

absence of vivianite formation, however, the larger 

phosphorus loads transported to Ashumet Pond 

are sustained by larger phosphorus concentrations. 

Without vivianite formation, the western lobe of 

large phosphorus concentrations (>90 μmol/L near 

disposal beds 8–16) extends to Ashumet Pond and the 

northern part of Ashumet Valley (not shown), which is 

inconsistent with smaller measured phosphorus 

concentrations in this area in 1993 (see fig. 8 for 

location of western lobe of large measured phosphorus 

concentrations).

The revised model simulates a peak in 

phosphorus load transported to Ashumet Pond between 

1990 and 2010, and although the calculations are better 

suited to qualitative interpretations, the maximum load 

during this period is estimated to be about 1,000 kg/yr. 

The phosphorus in this peak is derived mostly from 
30 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

loading the eastern disposal beds (beds 1–4) prior to 

1984. The location of the maximum measured 

phosphorus concentrations was approximately 

350 m downgradient from beds 1–4 in 1993 (fig. 8). 

This downgradient location demonstrates the transport 

of phosphorus over this distance during a period of 

at least 10 years, the time since the beds were used. 

However, differences between the simulated and 

measured phosphorus concentrations downgradient of 

beds 1–4 in 1993 (fig. 8) indicate that the timing and 

magnitude of the peak in phosphorus load transported 

to Ashumet Pond (fig. 10) may not be accurately 

predicted by the revised model. The simulated 

phosphorus concentrations are larger than the measured 

concentrations, and the large simulated concentrations 

extend to a greater distance downgradient toward 

Ashumet Pond than the large measured concentrations. 

These differences between the simulated and measured 

phosphorus concentrations indicate that the simulation 

results describe a larger peak in phosphorus load 

transported to Ashumet Pond that occurs and dissipates 

more quickly than actually occurs in the aquifer.

Limitations of the Revised Model

Of the many uncertainties involved in developing 

the revised model, the following appear to have the 

greatest effects on the load of phosphorus transported 

to Ashumet Pond: 

• The loading history of the beds—The 

concentration of phosphorus in the sewage 

effluent directly affects the load of phosphorus 

transported to Ashumet Pond. In addition, the 

history of disposal-bed usage—the specific beds 

and rates of disposal, particularly for individual 

beds on the eastern side—is critical to the 

estimation of phosphorus loads transported to 

Ashumet Pond. 

• The mechanism by which phosphorus is 

incorporated into the non-desorbable sediment 

fraction—The chemical form of the nondesorbable phosphorus is largely unknown, but 

extremely important. The long-term stability of 

the non-desorbable phosphorus is uncertain, 

although column experiments indicate it is not 

readily released. The kinetic expressions 

describing the incorporation of the non-desorbable 

form into the sediment are plausible and result in a 

reasonable match between simulated and 

measured phosphorus concentrations, but are 

otherwise poorly constrained. In the revised 

model, the presence or absence of vivianite 

formation switches the period of peak phosphorus 

load transported to Ashumet Pond from the 1990s 

to the 1970s (fig. 10). Although the loads after 

1990 are not greatly affected, the uncertainty 

makes it difficult to assess whether current loads 

are greater or less than the loads transported in 

previous decades. 

• The description of the flow system—Small 

changes to the model flow system can distribute 

the phosphorus loaded on the western disposal 

beds (beds 5–25) either to Ashumet Pond or to 

Ashumet Valley. The flow system is much more 

dynamic than simulated in the revised model, with 

seasonal variations in water level in Ashumet 

Pond and ground-water flow direction in the 

aquifer; however, a model that accounts for the 

dynamics of the flow system is beyond the scope 

of this report.

• The extent of decomposition of organic carbon 

following cessation of sewage disposal—In the 

revised model, decomposition occurs only to the 

extent that dissolved oxygen and nitrate are 

consumed in the presence of sorbed organic 

carbon. Sulfate reduction and methanogenesis 

have not been observed during the period of 

sewage disposal and are unlikely to be initiated 

following cessation of sewage disposal, when 

the reactivity of sorbed organic carbon is expected 

to decline. However, it is possible that organic 

decomposition will proceed further than the 

revised model allows. Reductive dissolution of 

iron- (and manganese-) oxyhydroxide minerals 

has been observed in the core of the sewage 

plume, but no reaction is included in the revised 

model to maintain reductive dissolution of iron 

and manganese minerals following cessation 

of sewage disposal. Iron and manganese 

oxyhydroxide dissolution will tend to raise the 

pH and increase the mobility of phosphorus. 
Summary and Conclusions 31

Simultaneously, increased iron concentrations 

may tend to sequester phosphorus in ferrous 

phosphate mineral phases. Thus, the extent of 

reductive dissolution of iron-oxyhydroxide 

minerals following cessation of sewage disposal 

and related effects on phosphorus load transported 

to Ashumet Pond are difficult to assess. 

• Cation-sorption reactions—These reactions are 

poorly defined, but affect the buffering of pH in 

ground water, which affects the sorption and 

mobility of phosphorus in the aquifer.

SUMMARY AND 

CONCLUSIONS

The subsurface transport of phosphorus 

introduced by the disposal of treated sewage effluent at 

the Massachusetts Military Reservation on western 

Cape Cod was simulated with a three-dimensional 

reactive-transport model. The simulation was used to 

estimate the load of phosphorus transported to 

Ashumet Pond during sewage-treatment-plant 

operation—1936 to 1995—and for 60 years following 

cessation of sewage disposal. The model accounted for 

spatial and temporal changes in water discharge from 

the sewage-treatment plant, ground-water flow, 

transport of associated chemical constituents, and a set 

of chemical reactions, including phosphorus sorption 

on aquifer materials, dissolution and precipitation of 

iron- and manganese-oxyhydroxide and iron-phosphate 

minerals, organic carbon sorption and decomposition, 

and cation sorption.

The three-dimensional reactive-transport model 

combined flow and transport calculations with 

chemical reactions. The flow and transport in the 

aquifer was modeled by using parameters consistent 

with previous flow and transport models of this area of 

Cape Cod, except that numerical dispersion was much 

larger than the physical dispersion estimated in 

previous studies. 

A base-case model was developed by using 

measured phosphorus concentrations from the aquifer 

and results from laboratory column experiments. A 

series of sensitivity simulations was made by varying 

the parameters and processes that were used in the 

base-case model. A revised model was developed by 

consideration of sensitivity-simulation results, 

comparison with additional laboratory column 

experimental results, comparison with field data, and 

inclusion of additional chemical processes.

Chemical reactions were included in the revised 

model to simulate the major reactions in the aquifer. 

The surface-complexation equilibrium constants and 

number of sorption sites for phosphorus and cations 

were adjusted in a reactive-transport model to fit the 

pH and phosphorus concentrations measured in column 

experiments. Iron- and manganese-oxyhydroxide 

minerals were assumed to react to equilibrium. Rates of 

organic carbon decomposition were adjusted to match 

the location and iron concentrations in the anoxic iron 

zone within the sewage plume. Additional reactions 

were incorporated in the model to account for the 

accumulation of sorbed organic carbon and nondesorbable phosphorus in aquifer sediments and to 

model denitrification as an irreversible process. 

All plausible simulations with the base-case and 

revised models indicate that the phosphorus load 

transported to Ashumet Pond was minimal before 

1965; the load transported to Ashumet Pond decreases 

following the period 2000–10, but is substantial for 

many decades; and most of the phosphorus introduced 

through the disposal beds resides in the aquifer 

sediments. The current calculations indicate a peak in 

phosphorus discharge to Ashumet Pond of about 

1,000 kg/yr during the 1990s; however, comparisons 

of modeled phosphorus concentrations with measured 

concentrations from samples taken in 1993 indicate 

that the peak in phosphorus load transported to 

Ashumet Pond may be larger and moving more quickly 

in the model simulations than is evident in the aquifer.

Simulation results from the revised threedimensional reactive-transport model are consistent 

with the current understanding of the loading history, 

laboratory experimental data, and measured 

concentrations in water samples from the aquifer. 

However, large uncertainties are associated with the 

loading of the sewage beds, the flow system, and the 

chemistry and sorption characteristics in the aquifer. As 

an example of the large uncertainties, the revised model 

indicates a small load of phosphorus during 1965–85, 
32 Reactive-Transport Simulation of Phosphorus in the Sewage Plume at the Massachusetts Military Reservation, Cape Cod, Massachusetts

but this small load is particularly sensitive to model 

parameters that specify flow conditions and the 

chemical process by which non-desorbable phosphorus 

is incorporated in the sediments. The uncertainties are 

sufficient that it is difficult to determine whether loads 

of phosphorus transported to Ashumet Pond in the 

1990s were greater or less than loads during the 

previous two decades. The revised model adequately 

reproduces the spatial distribution of phosphorus 

concentrations measured in 1993, the magnitude of 

changes in phosphorus concentration in a vertical 

profile of the aquifer near the disposal beds after 

cessation of sewage disposal, the observed iron zone in 

the sewage plume, the approximate flow of treated 

sewage effluent into Ashumet Valley, and laboratory 

column data for phosphorus sorption and desorption. 

REFERENCES CITED

Barber, L.B., II, 1990, Geochemical heterogeneity in a 

glacial outwash aquifer—Effect of particle size and 

mineralogy on sorption of nonionic organic solutes: 

Boulder, University of Colorado, Department of 

Geological Science, Ph.D. dissertation, 238 p. 

Ceazan, M.L., Thurman, E.M., and Smith, R.L., 1989, 

Retardation of ammonium and potassium transport 

through a contaminated sand and gravel aquifer—The 

role of cation exchange: Environmental Science and 

Technology, v. 23, no. 11, p. 1402–1408.

Coston, J.A., Fuller, C.C., and Davis, J.A., 1995, Pb2+ and 

Zn2+ adsorption by a natural aluminum- and ironbearing surface coating on an aquifer sand: Geochimica 

et Cosmochimica Acta, v. 59, no. 17, p. 3535–3547.

Dzombak, D.A., and Morel, M.M., 1990, Surface 

complexation modeling: New York, John Wiley and 

Sons, 393 p.

Fishman, M.J., and Friedman, L.C., 1985, Methods for 

determination of inorganic substances in water and 

fluvial sediments: Techniques of Water-Resources 

Investigations of the U.S. Geological Survey, book 5, 

chap. 1, p. 477–478.

Fuller, C.C., Davis, J.A., Coston, J.A., and Dixon, Eleanor, 

1996, Characterization of metal adsorption variability 

in a sand and gravel aquifer, Cape Cod, Massachusetts, 

U.S.A: Journal of Contaminant Hydrology, v. 22, 

no. 3-4, p. 165–187.

Garabedian, S.P., LeBlanc, D.R., Gelhar, L.W., and Celia, 

M.A., 1991, Large-scale natural-gradient tracer test 

in sand and gravel, Cape Cod, Massachusetts—

2. Analysis of spatial moments for a nonreactive tracer: 

Water Resources Research, v. 27, no. 5, p. 911–924.

Gschwend, P.M., and Reynolds, M.D., 1987, Monodisperse 

ferrous phosphate colloids in an anoxic groundwater 

plume: Journal of Contaminant Hydrology, v. 1, no. 3, 

p. 309–327.

Harbaugh, A.W., and McDonald, M.G., 1996, User’s 

documentation for MODFLOW-96, an update to U.S. 

Geological Survey Modular Finite-Difference GroundWater Flow Model: U.S. Geological Survey Open-File 

Report 96-485, 56 p.

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.

LeBlanc, D.R., 1984, Sewage plume in a sand and gravel 

aquifer, Cape Cod, Massachusetts: U.S. Geological 

Survey Water-Supply Paper 2218, 28 p.

LeBlanc, D.R., Garabedian, S.P., Hess, K.M., Gelhar, L.W., 

Quadri, R.D., Stollenwerk, K.G., and Wood, W.W., 

1991, Large-scale natural-gradient tracer test in sand 

and gravel, Cape Cod, Massachusetts—Experimental 

design and observed tracer movement: Water Resources 

Research, v. 27, no. 5, p. 895–910.
References Cited 33

LeBlanc, D.R., Hess, K.M., Kent, D.B., Smith, R.L., Barber, 

L.B, Stollenwerk, K.G., and Campo, K.W., 1999, 

Natural restoration of a sewage plume in a sand and 

gravel aquifer, 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, Subsurface Contamination from Point 

Sources: U.S. Geological Survey Water-Resources 

Investigations Report 99-4018C, p. 245–259.

Masterson, J.P., and Walter, D.A., 2000, Delineation of 

ground-water recharge areas, western Cape Cod, 

Massachusetts: U.S. Geological Survey WaterResources Investigations Report 00-4000, 1 pl.

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground-water flow 

model: Techniques of Water-Resources Investigations 

of the U.S. Geological Survey, book 6, chap. A1, 

variously paged. 

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 Water-Resources Investigations Report 99-

4259, 312 p.

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

Coupling the geochemical model PHREEQC with a 3D 

multi-component solute transport model, in Program 

and Abstracts: State College, PA, V.M. Goldschmidt 

Conference, May 24–26, 1995, p. 77–78.

Parkhurst, D.L., and Kipp, K.L., 2002, Parallel processing 

for PHAST—A three-dimensional reactive-transport 

simulator, in Hassanizadeh, S.M., Schlotting, R.J., 

Gray, W.H., and Pinder, G.F., eds., Computational 

methods in water resources, v. 2, Developments in 

Water Science, n. 47: Amsterdam, The Netherlands, 

Elsevier, p. 711–718.

Poeter, E.P., and Hill, M.C., 1999, UCODE, A computer 

code for universal inverse modeling: Computer in 

Geosciences, v. 25, p. 457–462.

Savoie, Jennifer, and LeBlanc, D.R., 1998, Water-quality 

data and methods of analysis for samples collected near 

a plume of sewage-contaminated ground water, 

Ashumet Valley, Cape Cod, Massachusetts, 1993–94: 

U.S. Geological Survey Water-Resources Investigations 

Report 97-4269, 208 p.

Stollenwerk, K.G., 1995, Modeling the effects of variable 

ground-water chemistry on adsorption of molybdate: 

Water Resources Research, v. 31, no. 2, p. 347–357.

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, Subsurface Contamination from Point 

Sources: U.S. Geological Survey Water-Resources 

Investigations Report 99-4018C, p. 371–381.

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.