Groundwater and Surface-Water Interaction and Effects of 

Pumping in a Complex Glacial-Sediment Aquifer, Phase 2, 

East-Central Massachusetts

Prepared in cooperation with the Town of Framingham, Massachusetts

Scientific Investigations Report 2015–5174

U.S. Department of the Interior

U.S. Geological Survey
Cover. Pod Meadow Pond, Framingham, Massachusetts. Ice-free area of open water in the foreground is caused by warmer 

groundwater discharging to the pond. Photograph by Jack Eggleston, U.S. Geological Survey, February 2011.
Groundwater and Surface-Water 

Interaction and Effects of Pumping in 

a Complex Glacial-Sediment Aquifer, 

Phase 2, East-Central Massachusetts

By Jack R. Eggleston, Phillip J. Zarriello, and Carl S. Carlson

Prepared in cooperation with the Town of Framingham, Massachusetts

Scientific Investigations Report 2015–5174

U.S. Department of the Interior

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

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

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

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

resources, natural hazards, and the environment—visit http://www.usgs.gov or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod/.

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

U.S. Government.

Although this information product, for the most part, is in the public domain, it also may contain copyrighted materials 

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Eggleston, J.R., Zarriello, P.J., and Carlson, C.S., 2015, Groundwater and surface-water interaction and effects 

of pumping in a complex glacial-sediment aquifer, Phase 2, east-central Massachusetts: U.S. Geological Survey 

Scientific Investigations Report 2015–5174, 38 p., http://dx.doi.org/10.3133/sir20155174.

ISSN 2328-0328 (online)
iii

Acknowledgments

This study was supported by the Town of Framingham, Massachusetts. The authors extend 

thanks to Peter Newton (Bristol Engineering Advisors, Inc.), Eric Johnson (Town of Framingham, 

Mass.), and Diane Stokes (City of Cambridge, Mass.), for sharing their knowledge of the study 

area and collecting and providing field information used in the study.

v

.....

Contents

Acknowledgments........................................................................................................................................iii

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

Introduction.....................................................................................................................................................1

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

Study Area..............................................................................................................................................2

Previous Investigations........................................................................................................................2

Data .............................................................................................................................................................5

Seismic Survey......................................................................................................................................5

Specific Conductivity Measurements ...............................................................................................5

Streamflow Measurements ..............................................................................................................10

Observation Wells...............................................................................................................................11

Groundwater Elevations ....................................................................................................................11

Groundwater Model Modifications ..........................................................................................................13

Layering and Active Model Area .....................................................................................................14

Boundary Conditions and Surface-Water Body Representation................................................14

Hydraulic Parameter Assignment....................................................................................................17

Model Calibration................................................................................................................................20

Calibration Targets and Weights.............................................................................................20

Model Uncertainty and Sensitivity Analysis .........................................................................21

PEST Parameter Sensitivities .........................................................................................21

Numerical Instabilities.....................................................................................................21

Model Error........................................................................................................................24

Effect of Parameter Values on Streamflow Timing and on Induced Recharge 

from Lake Cochituate...................................................................................................25

Model Limitations.......................................................................................................................28

Effects of Pumping.......................................................................................................................................28

Steady State Simulations ..................................................................................................................30

Transient Simulations.........................................................................................................................32

Summary and Conclusions.........................................................................................................................37

References Cited..........................................................................................................................................37

Figures

1. Map showing glacial-sediment aquifer study area in east-central Massachusetts........3

2. Map showing groundwater levels in the lower aquifer for the study area in 

east-central Massachusetts.......................................................................................................4

3. Map showing seismic reflection transect lines made on May 30, 2013, 

Lake Cochituate, east-central Massachusetts........................................................................6

4. Example of marine transect profiles of Lake Cochituate, east-central 

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

5. Map showing locations of water specific conductance measurements made in

Lake Cochituate and near shore of Pod Meadow Pond in east-central 

Massachusetts..............................................................................................................................8
vi

6. Map showing locations of borehole drilling sites and wells installed in 2011 in the

Birch Road well area, east-central Massachusetts.............................................................12

7. Vertical cross sections of groundwater model layers A, north-south and 

B, east-west.................................................................................................................................15

8. Map showing boundary conditions in the Phase 2 groundwater model ..........................16

9. Maps showing horizontal hydraulic conductivity of groundwater model for 

A, Layer 1, B, Layer 2, and C, Layer 3.......................................................................................18

10. Graph showing observed water level compared to simulated steady-state

groundwater levels for the calibrated Phase 2 groundwater model .................................28

11. Map showing simulated groundwater levels under existing conditions in

model Layer 3 and error in simulated heads in model Layers 1 and 3 for the 

study area in east-central Massachusetts ............................................................................29

12. Graphs showing simulated monthly A, flow components under existing average

conditions and with additional pumping from Birch Road wells and B, changes in

flows..............................................................................................................................................33

13. Graph showing simulated streamflow response of the Sudbury River after 1 month

of pumping the Birch Road wells at 4.9 cubic feet per second ..........................................34

14. Graphs showing simulated A, monthly streamflows and streamflow depletions in

the Sudbury River at the model exit, with hypothetical groundwater withdrawals at

the Birch Road wells, under B, average, and C, dry conditions .........................................36

Tables

1. Specific conductivity profile measurements made on May 30, 2013, in Lake

Cochituate, east-central Massachusetts .................................................................................9

2. Specific conductivity measurements made on June 5, 2013, near the shoreline of

Pod Meadow Pond, east-central Massachusetts ................................................................10

3. Discharge measurements at the outflows from Pod Meadow Pond and Dudley

Pond, east-central Massachusetts .........................................................................................11

4. Observation wells installed in 2011 in the Birch Road well aquifer area, east-central

Massachusetts............................................................................................................................12

5. Assigned monthly stream inflow to the model, in cubic feet per second, under

average (25th percentile) and dry (10th percentile) climatic conditions...........................17

6. Calibrated hydraulic parameter values in the Phase 2 model, with comparison to

Phase 1 values ............................................................................................................................19

7. Calibrated streambed hydraulic conductivity values and comparison to Phase 1

values............................................................................................................................................19

8. Average observed and simulated groundwater levels and flow just downstream

from the Pod Meadow Pond outlet for the calibrated groundwater model of an

aquifer in east-central Massachusetts...................................................................................22

9. Observed and simulated maximum groundwater-level drawdown during 2006

aquifer pumping test, Phase 2 groundwater model of the aquifer in east-central

Massachusetts............................................................................................................................24

10. Calibration target groups used in automated parameter estimation, with

representative assigned weights by group............................................................................24

11. Parameter sensitivity values determined by PEST for the calibrated Phase 2

groundwater model of a shallow aquifer in east-central Massachusetts........................26

12. Summary of model errors..........................................................................................................30
vii

13. Sensitivity of pumping-induced recharge from Lake Cochituate to the aquifer 

and of streamflow response time to select steady-state and transient model 

parameters...................................................................................................................................31

14. Groundwater model scenarios used to evaluate the aquifer and streamflow

response to hypothetical withdrawals at the Birch Road wells, east-central

Massachusetts............................................................................................................................32

15. Steady-state simulated water budgets with and without hypothetical Birch Road

well withdrawals at a continuous rate of 4.9 ft3/s, east-central Massachusetts ............32

16. Simulated streamflow in the Sudbury River at the exit of the Phase 2 

groundwater model ....................................................................................................................35

Conversion Factors

[Inch/Pound to International System of Units]

Multiply By To obtain

Length

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

inch (in.) 2.54 centimeter (cm)

Area

square mile (mi2

) 2.590 square kilometer (km2) 

Flow rate

cubic foot per day (ft3

/d) 0.00001157 cubic foot per second (ft3/sec)

cubic foot per second (ft3

/s) 0.6463 million gallons per day (Mgal/d)

million gallons per day (Mgal/d) 1.547 cubic foot per second (ft3

/s)

foot per day (ft/d) 0.032854 inch per year (in/yr)

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

(NAVD 88).

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

‘Elevation’, as used in this report, refers to distance above the vertical datum (NAVD 88).

Specific conductance is given in microsiemens per centimeter at 25 degrees Celsius (μS/cm at 

25 °C).
viii

Abbreviations

EIR Environmental Impact Report 

EPA U.S. Environmental Protection Agency

GPR ground penetrating radar

GPS global positioning system

HSPF Hydrological Simulation Program—Fortran

Kh horizontal hydraulic conductivity (feet per day) 

Ks streambed hydraulic conductivity (feet per day)

LAK MODFLOW Lake package 

MassDEP Massachusetts Department of Environmental Protection

MassDCR Massachusetts Department of Conservation and Recreation

MEPA Massachusetts Environmental Policy Act

MODFLOW-NWT modular groundwater flow model with Newtonian solver

MWRA Massachusetts Water Resources Authority

NOAA National Oceanic and Atmospheric Administration

NPS National Park Service

NWIS National Water Information System 

PEST parameter estimation software

SFR MODFLOW Streamflow Routing package 

Ss specific storage (1/feet)

Sy specific yield (percent)

t

50 time for 50 percent of streamflow depletion to take place (day)

USFWS U.S. Fish and Wildlife Service

USGS U.S. Geological Survey
Groundwater and Surface-Water Interaction and Effects 

of Pumping in a Complex Glacial-Sediment Aquifer, 

Phase 2, East-Central Massachusetts

By Jack R. Eggleston, Phillip J. Zarriello, and Carl S. Carlson

Abstract

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

Town of Framingham, Massachusetts, has investigated the 

potential of proposed groundwater withdrawals at the Birch 

Road well site to affect nearby surface water bodies and 

wetlands, including Lake Cochituate, the Sudbury River, and 

the Great Meadows National Wildlife Refuge in east-central 

Massachusetts. In 2012, the U.S. Geological Survey developed 

a Phase 1 numerical groundwater model of a complex glacialsediment aquifer to synthesize hydrogeologic information and 

simulate potential future pumping scenarios. The model was 

developed with MODFLOW-NWT, an updated version of a 

standard USGS numerical groundwater flow modeling program that improves solution of unconfined groundwater flow 

problems. The groundwater model and investigations of the 

aquifer improved understanding of groundwater–surface-water 

interaction and the effects of groundwater withdrawals on 

surface-water bodies and wetlands in the study area. The 

initial work also revealed a need for additional information 

and model refinements to better understand this complex 

aquifer system.

In this second phase of the study, the original groundwater flow model was revised to improve representation of 

groundwater and surface-water hydrology, stabilize the model, 

and reduce model error. The model was simplified by reducing 

the number of layers from 5 to 3 and adding the MODFLOW 

lake package (LAK) to simulate Lake Cochituate and Pod 

Meadow Pond and better represent interaction between the 

lakes and the aquifer. Model revisions improved stability and 

shortened run times, allowing use of automated parameter estimation software (PEST) to further refine the model hydraulic 

parameters and reduce simulation errors. 

Model simulations indicate that under average baseflow conditions, the Birch Road wells have a small effect on 

flow in the Sudbury River during most months, even at the 

maximum pumping rate of 4.9 ft3/s (3.17 Mgal/d). Maximum 

percent streamflow depletion in the Sudbury River caused 

by simulated pumping takes place during simulated drought 

conditions, when streamflow decreased by as much as 

21 percent under maximum continuous pumping. Simulations 

also indicate that groundwater withdrawals at the Birch Road 

site could be managed so that adverse streamflow impacts 

are substantially ameliorated. Under the most ecologically 

conservative simulated drought conditions, simulated streamflow depletion was reduced from 21 percent to 3 percent by 

pumping at the maximum rate for 6 months rather than for 

12 months. Simulations that return 10 percent of the Birch 

Road well withdrawals to Pod Meadow Pond indicate a modest reduction in the Sudbury River streamflow depletion and 

provide a larger percentage increase to streamflow just downstream of the pond. The groundwater model also indicates that 

well locations can have a large effect on the sustainable pumping rate and so should be chosen carefully. The model provides 

a tool for evaluating alternative pumping rates and schedules 

not included in this analysis.

Introduction

The Town of Framingham in east-central Massachusetts 

operated several groundwater water-supply wells in the area 

known as the Birch Road site (fig. 1) from 1939 until about 

1979. In 2009, in accordance with the Massachusetts Environmental Policy Act (MEPA), the Town of Framingham filed an 

Environmental Impact Report (EIR) with the Massachusetts 

Department of Environmental Protection to reactivate these 

supply wells (SEA Consultants, Inc., 2009). The growing 

recognition of the interconnection between groundwater and 

surface-water resources and the potential effects of groundwater withdrawals on surface water led to concerns raised by 

the National Park Service (NPS), the U.S. Fish and Wildlife 

Service (USFWS), the U.S. Environmental Protection Agency 

(the EPA), state and local agencies, and environmental interest 

groups and citizens. The impacts from groundwater withdrawals on the nearby Sudbury River and the downstream Concord 

River are of particular concern to Federal interests because 

of the Great Meadows National Wildlife Refuge and the 

Minute Man National Historical Park, respectively. In addition, the Sudbury River is designated by Congress as “Wild 

and Scenic,” requiring special resource protection (U.S. Congress, 1999). In response to those concerns and notice by 
2 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

the U.S. Department of the Interior Office of the Solicitor to 

specifically address Federal interests (U.S. Department of Interior, 2009), the Town of Framingham and the U.S. Geological 

Survey (USGS) agreed to collaborate on a study of the aquifer 

system to better understand the potential effects of groundwater withdrawals on nearby surface-water resources.

In 2012, the USGS published a report (Eggleston and 

others, 2012) under the cooperative agreement with the 

Town of Framingham that characterized the complex glacialsediment aquifer system in the area of the proposed pumping 

wells. The report also described a groundwater flow model 

(MODFLOW–NWT) developed to better understand the connection between groundwater and surface water that was used 

to simulate the effects of pumping on the aquifer and nearby 

surface-water resources. Simulations indicated about one third 

of the proposed withdrawal is from induced infiltration of 

water from Lake Cochituate and the rate of streamflow depletion in the Sudbury River downstream from the oxbow (fig. 1) 

is about equal to the rate of pumping under steady state conditions, but streamflow depletion changed quickly in response 

to changes in pumping. The report identified the need for 

additional data collection and refinement of the model, which 

could affect the present understanding of groundwater/surfacewater interactions in the study area, and the report identified 

the need to reexamine withdrawal scenarios tested. 

Purpose and Scope

This report describes the second phase of the cooperative 

study to further improve the understanding of the complex glacial-aquifer system and the effects of groundwater withdrawals 

on surface-water resources in the vicinity of the Birch Road 

well site. This study was designed to improve understanding 

of the shallow groundwater system in the vicinity of the Birch 

Road well site and to build a groundwater model for assessing the effects of proposed pumping on nearby surface-water 

features including Lake Cochituate and the Sudbury River. 

This report describes the new information obtained to characterize the aquifer and its interaction with the surface-water 

system and the modifications made to the previously developed groundwater flow model (MODFLOW-NWT) of the 

glacial-sediment aquifer in the Birch Road area in east-central 

Massachusetts (Eggleston and others, 2012). This report also 

presents the revised model (herein referred to as Phase 2) and 

simulation results used to further quantify groundwater and 

surface-water interaction under present conditions and under 

various pumping scenarios at the Birch Road site to as much 

as a maximum of 3.17 million gallons per day (Mgal/d). The 

model scenarios include alternative pumping schedules and 

return of a portion of the withdrawal, which could reduce 

the effects of groundwater withdrawals during seasonal low 

streamflows when stream ecosystems are most sensitive. 

Study Area

The study area is 16 miles (mi) west of Boston in eastcentral Massachusetts mostly in the towns of Framingham 

and Wayland (fig. 1). The groundwater aquifer is a complex 

glacial-fill aquifer in a bedrock valley trending north to south 

with a large bedrock outcrop (bedrock island) near the center 

of the study area (fig. 1). Several large ponds, including Lake 

Cochituate, and numerous streams overlay the aquifer. The 

Sudbury River flows from the southwest toward the northeast 

through the study area. The active groundwater model area is 

about 5.5 square miles (mi2) and includes about 1.7, 3.4, 0.4, 

and 0.01 mi2 in the towns of Framingham, Wayland, Sudbury, 

and Natick, respectively. The active model area was reduced 

from about 6.1 mi2 in the Phase 1 study to 5.5 mi2 in the 

Phase 2 study to improve model stability as described in the 

Groundwater Model Modifications section.

The Sudbury River is the primary surface-water drainage feature and sets the natural base groundwater level in the 

study area (fig. 2). The Great Meadows National Wildlife 

Refuge includes lands adjacent to the Sudbury River north of 

the oxbow and wetland areas in the northern part of the study 

area, although most of the refuge is to the north of the study 

area (fig. 1). The north pond of Lake Cochituate covers about 

0.3 mi2 and is the largest surface-water feature in active model 

area. A mostly natural causeway forms a divide between the 

north pond and two ponds to the south. The northern pond of 

Lake Cochituate has a contributing drainage area of 17.5 mi2

and drains to the Sudbury River through the 1.4 mi long 

Cochituate Brook. Additional surface-water bodies include 

Dudley Pond, Pod Meadow Pond, and Heard Pond, which 

drain to the Sudbury River. 

The aquifer is a complex mix of stratified glacially 

deposited sediments, including melt water deltaic deposits 

and proglacial lake deposits that range in texture from clay to 

coarse gravel and boulders. Most of these deposits however 

are medium to fine sands, which form the primary aquifer. The 

glacial geomorphological sequences and depositional features 

are described in the Phase 1 study by Eggleston and others 

(2012).

Previous Investigations

Previous studies describing geology and hydrology of the 

area were summarized in the Phase 1 report by Eggleston and 

others (2012). The Phase 1 report provides additional information on the bedrock geology and depth to bedrock, interpretation of the stratigraphy of the surficial glacial and post-glacial 

deposits, surface-water features, and water use. The Phase 1 

report also documents the groundwater flow model development and parameterization, which is the basis for the model 

presented in this Phase 2 study, and the simulated effects of 

select groundwater withdrawals on select surface-water features, which were reexamined in the Phase 2 study.
Introduction 3

Figure 1. Glacial-sediment aquifer study area in east-central Massachusetts. Sudbury River 

 Cochituate Brook

Sudbury River 01098500 and

01098499

Dudley Pond

Lake Cochituate

(north pond)

Heard Pond

Pod Meadow 

Pond

GREAT 

MEADOWS 

NATIONAL

WILDLIFE 

REFUGE

Bedrock island Bedrock island

Wayland

Framingham

Natick

Sudbury

MA-WKW-2

71°24’ 71°22’

42°21’

42°19’

0 1,000 2,000 FEET

0 250 500 METERS

Base from U.S. Geological Survey and Massachusetts Geographic Information System data, variously dated, various scales

Massachusetts State Plane Coordinate System, Mainland zone

Oxbow

Pod Meadow

01098530

Birch Road wells

#

# #

#

Birch Road test wells #

##

MV–1

HH–2 HH–1

Study area

MASSACHUSETTS

Boston

Wetland

Inactive model area

Great Meadows National 

Wildlife Refuge

Town boundary

Birch Road test well

Wayland production well and number

U.S. Geological Survey (USGS) 

monitoring well and name

USGS streamgage and number

USGS partial record streamgage

EXPLANATION

#

#

01098500

MV–1

MA-WKW-2

Phase I model 

boundary
4 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Figure 2. Groundwater levels (simulated heads) in the lower aquifer (model Layer 3) for the study area in east-central Massachusetts.

90

Sudbury River 

 Cochituate Brook

Sudbury River

Lake Cochituate

(north pond)

Pod Meadow 

Pond

Wayland

Framingham

Natick

Sudbury

Oxbow

Pod Meadow

130

130

120

140

150

160

170

120

140

150

160

130

130

140

150

160

170

170

120

#

# #

#

#

##

Bedrock island

Birch Road

test wells

MV–1

HH–2

HH–1

Base from U.S. Geological Survey and Massachusetts Geographic Information System data, variously dated, various scales

Massachusetts State Plane Coordinate System, Mainland zone

71°24’ 71°22’

42°21’

42°19’

0 1,000 2,000 FEET

0 250 500 METERS

108

118

128

138

148

158

168

178

188

120

Inactive model area

Town boundary

Equal head contour 

and value—Contour 

interval 2 feet

General direction of 

groundwater flow

Birch Road test well

Wayland production 

well and number

EXPLANATION

#

#

MV–1

Simulated head, in feet
Data 5

Data

In response to data needs identified in Phase 1, data were 

collected to fill data gaps to the extent possible and to support 

work in Phase 2. Data collected in Phase 2 include a marine 

seismic survey, specific conductivity measurements of Lake 

Cochituate and Pod Meadow Pond, water-level measurements 

from newly installed wells, and flow from Pod Meadow Pond. 

In addition, groundwater levels used to calibrate the Phase 1 

model were adjusted to better reflect steady-state conditions.

Seismic Survey

Seismic reflection profiling was completed in Lake 

Cochituate on May 30, 2013, to better characterize lake bottom sediments and determine the extent of fine-grained lake 

bed sediment deposits (gyttja). The extent of gyttja deposits is 

important to know because gyttja has low hydraulic conductivity and can reduce the rate of flow between the lake and 

underlying aquifer. Previous attempts to characterize the lake 

sediments by use of ground penetrating radar (GPR), which is 

generally preferable for mapping coarse-grained lake bed sediments in fresh water systems, were unsuccessful because the 

high specific conductivity of the lake water quickly attenuated 

the radar signal. Marine (underwater) seismic waves scatter in 

response to gaseous fine-grained lake bed sediment deposits 

on the bottom of glacial ponds so they were used as a surrogate to map the extent of the gyttja deposits in the north pond 

of Lake Cochituate. 

Seismic response was measured with a chirp profiling 

system (EdgeTech XSTAR SB-216) that emits an acoustic 

signal in the water at regular intervals with a 4–20 kilohertz 

operating frequency. A towfish pulled behind a boat emits an 

acoustic energy, which is reflected from different boundaries 

and recorded by a hydrophone. The return signal is digitally 

recorded for later processing using DelphSeismic software. 

Nineteen transects made on May 30, 2013, generally traversed 

the lake from north to south and east to west (fig. 3). The 

profiles covered about 6 miles and a total of 18,732 seismic 

traces were collected. The position of the tow fish was tracked 

continuously during the profiling using a global positioning system (GPS). The profile transects were played back 

with Edge Tech Corp software (v. 3.52) from which areas of 

gaseous deposits were identified by the highly scattered seismic traces. Examples of the profile transects and interpreted 

gaseous deposits are shown for lines 11 (east-west line) and 

18 (north-south line) in figure 4. The extent of the gaseous 

deposits were mapped and used as an approximation of gyttja 

deposits extent.

The seismic profiles indicated an extensive area of 

gaseous gyttja deposits in the widest and deepest part of 

the north pond of Lake Cochituate (fig. 3). Smaller areas of 

gyttja deposits were interpolated in the northern part of the 

lake between underwater mounds or ridges. The mounds or 

ridges were most pronounced in an east–west direction where 

the lake shore narrows, however some mounds also appear 

to traverse in a north–south direction. One interpretation of 

how these underwater mounds formed can be made from the 

sequence of the glacial retreat (Clapp, 1904). Lake Cochituate 

formed as a kettle pond where a large block of ice persisted as 

the main glacier retreated and paused farther to the north. The 

ice block limited sediments from being deposited in proglacial 

Lake Charles and, when it eventually melted, it left a depression that formed Lake Cochituate (Gay, 1985). During this 

process the ice block that later created the lake was likely 

fractured where melt-water sediments collected. When the ice 

block finally melted these sediments formed the underwater 

mounds or ridges that are likely a mix of fine and coarse grain 

sediments, the locations of these deposits are inferred from 

sediment cores and geologic interpretation, but not known 

exactly. Hydraulically, the mounds or ridges may create preferential flow paths between the lake and underlying aquifer 

especially compared to the gyttja deposits that cover about 

30 percent of the lake bottom. Because gyttja deposits do not 

cover most of the lake bottom they are unlikely to substantially limit flow between the aquifer and the lake, so values of 

lakebed conductance were not constrained to low values in the 

model. 

Specific Conductivity Measurements

The Phase 1 study identified groundwater seepage into 

Pod Meadow Pond as a potentially important hydraulic 

exchange with Lake Cochituate. To further investigate possible groundwater flow from Lake Cochituate to Pod Meadow 

Pond, specific conductivity profile measurements were made 

on May 30, 2013, at four locations in Lake Cochituate at 5 to 

10 foot (ft) vertical intervals from near the water surface to 

near the lake bottom (fig. 5 and table 1). Values ranged from 

499 to 545 microsiemens per centimeter (μS/cm) at 25 degrees 

Celsius with a mean of 519 μS/cm that decreased slightly with 

depth, and were generally consistent at the four measurement 

locations.

Additional specific conductivity measurements were 

made along the edge of Pod Meadow Pond (table 2) on June 5, 

2013. The specific conductance values in Pod Meadow Pond 

were similar to those in Lake Cochituate although slightly 

lower at most locations (fig. 5 and tables 1–2). Anomalously 

lower specific conductance values were measured at sites 4 

and 5 (312 and 140 μS/cm, respectively) and were 30 and 

70 percent lower than the median specific conductance 

measurement of 435 μS/cm at Pod Meadow Pond. The lower 

specific conductance at sites 4 and 5 may provide another 

indication that groundwater discharge may be focused in the 

south-western part of the pond. Other indicators of focused 
6 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Figure 3. Seismic reflection transect lines made on May 30, 2013, Lake Cochituate (north pond), east-central 

Massachusetts.

90

Lake

Cochituate

0 1,000 2,000 FEET 

0 250 500 METERS 

71°22’

42°19’

42°20’

71°23’

Base from U.S. Geological Survey and Massachusetts Geographic 

Information System data, variously dated, various scales

Massachusetts State Plane Coordinate System, Mainland zone

L01

L02

L03 L04

L05

L06

L07 L08

L09

L10

L11

L12

L13

L14

L15

L16

L17

L18

L19

Interpreted gyttja deposits from marine seismic

Path of tow fish and line numbers for marine seismic

EXPLANATION

Pod Meadow

Pond
Data 7

Figure 4. Example of marine transect profiles (lines 11 and 18, figure 3) of Lake Cochituate (north pond), east-central 

Massachusetts.

0

10

20

30

Depth, in meters

0 100 200 300 400 500 600 700 800 900 1,000 1,100 1,200 1,300 1,400 1,500 1,600 1,700

Distance, in meters

NORTH SOUTH

Lake bottom

VERTICAL EXAGGERATION X 25.8

B. Line 18 (fig. 3)

Reflected waves are scatter indicative of gas 

in sediments characteristic of gyttja deposits 

0

10

20

30

WEST EAST

Depth, in meters

0 100 200 300 400 500 600

Distance, in meters

A. Line 11 (fig. 3)

Lake bottom

VERTICAL EXAGGERATION X 7.6
8 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

90

Lake

Cochituate

Lake elevation 138.5 feet

Pod Meadow

Pond

Pond elevation

125.8 feet

0 1,000 2,000 FEET 

0 250 500 METERS 

71°22’

42°19’

42°20’

71°23’

Base from U.S. Geological Survey and Massachusetts Geographic 

Information System data, variously dated, various scales

Massachusetts State Plane Coordinate System, Mainland zone

Vertical coordinate information is referenced to the North American 

Vertical Datum of 1988 (NAVD 88)

#

#

#

#

1

2

3

4

Bedrock

Island

#

3 Location and number of specific conductivity profile (tables 1 and 2)

EXPLANATION

1

2

3

4

13

10

12 11

9

8

7

6

5

EXPLANATION

Location and number of 

specific conductivity 

measurement

0 500 FEET

0 150 METERS

#

#

#

#

#

#

#

##

#

#

#

# #

Pod Meadow

Pond

13

#

Figure 5. Locations of water specific conductance measurements made in Lake Cochituate and near shore of Pod Meadow Pond in 

east-central Massachusetts.
Data 9

Table 1. Specific conductivity profile measurements made on May 30, 2013, in Lake Cochituate, east-central Massachusetts.

[ft, feet; WT, water temperature; °C, degrees Celsius, SC, specific conductivity; μS/cm, microsiemens per centimeter at 25 degrees Celsius]

Site 

(fig. 5) Latitude Longitude Water depth (ft)

WT 

(°C)

SC 

(μS/cm) Remarks

1 42°19ʹ33.7ʺ -71°22ʹ55.0ʺ 1 20.0 526 Northern end of lake.

7 17.8 535

17 13.0 512

27 8.9 499

2 42°19ʹ22.4ʺ -71°22ʹ52.8ʺ 1 22.0 528

7 18.0 534

13 16.6 531

19 11.8 511

25 8.9 505

31 8.2 506

3 42°18ʹ58.6ʺ -71°22ʹ48.9ʺ 1 21.0 532 Near lake outlet.

6 18.0 534

12 16.6 535

19 13.0 511

25 8.3 507

31 7.5 505

37 7.4 501

4 42°18ʹ56.2ʺ -71°22ʹ26.7ʺ 1 21.2 537 Not at lake bottom.

6 20.0 536

12 17.8 545

18 14.5 521

24 9.7 508

30 8.0 505

36 7.6 503

516 Median.
10 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

groundwater discharge are open water observed in the same 

area in February 2011, when other areas of the pond were 

frozen (Eggleston and others, 2012), and groundwater seepage noted along the shore of the pond during the June 2013 

measurements. If the south-western part of the pond is an area 

of focused groundwater discharge, there are several possible 

contributing causes. The horizontal distance between the north 

shore of the lake and the south shore of the pond is small, only 

about 700 ft, whereas the difference in lake and pond surface 

elevations, about 12 ft, providing a relatively steep hydraulic 

gradient to drive groundwater flow. The bedrock island just to 

the north of the pond (fig. 1) and the rise in the bedrock valley 

to the north may force groundwater flow upward in this area. 

The similarity between specific conductivity values 

measured at the bottom of Lake Cochituate near the north end 

(499) and at areas of observed groundwater discharge in Pod 

Meadow Pond (average of 466 μS/cm for sites 6–8, 10 in table 

2 further support the hypothesis of rapid flow of recharge from 

Lake Cochituate to the aquifer that subsequently discharges to 

Pod Meadow Pond.

Streamflow Measurements

To better understand rates of groundwater discharge 

to the Sudbury River, a series of streamflow measurements, 

referred to as seepage measurements, were made on July 

13, 2012, along a 1.3-mi reach of the Sudbury River. Flows 

in the Sudbury River at Saxonville (0109850) were then at 

about the 90-percent flow duration based on daily streamflow 

records, indicating base-flow conditions when groundwater 

is the primary source of water to the river. Four measurements were made in the Sudbury River (fig. 1) beginning at 

streamgage (01098530) at Saxonville (0 miles [mi]), about 

halfway between the streamgage and the oxbow (0.35 mi), at 

the oxbow (0.75 mi), and near where the river bends sharply 

to the east (1.3 mi). Streamflow at these locations were 16.1, 

15.5, 16.5, and 15.2 ft3

/s, respectively, and were rated good 

(plus or minus 5 percent) to fair-poor (plus or minus 8 percent) 

because of the very shallow depths. Seepage measurements 

indicate a loss between the most upstream and downstream 

sites ranging from -0.6 to 1.3 ft3/s but, given the potential 

measurement error, this may not be an accurate measure of 

groundwater discharge to the river.

Streamflow and groundwater-level measurements have 

continued to be collected since the Phase 1 study ended. Daily 

streamflow records have been collected for Cochituate Brook 

below Lake Cochituate at Framingham, Massachusetts (USGS 

station 01098500) and daily stage records have been collected 

for Lake Cochituate at Framingham, Mass., (USGS station 

01098499). Both stations are operated and maintained in a 

cooperative agreement with the Town of Framingham and the 

USGS. Streamflow exiting Pod Meadow Pond was measured 

on a monthly basis through December 2012 to extend the flow 

record and improve descriptive flow statistics (table 3). The 

additional Pod Meadow Pond outlet streamflow measurements 

were very important because average Pod Meadow Pond 

outlet flow value is the only flux observation used in model 

calibration, as described later (Model Calibration section).

Table 2. Specific conductivity measurements made on June 5, 2013, near the shoreline of Pod Meadow Pond, east-central 

Massachusetts.

[WT, water temperature; °C, degrees Celsius; SC, specific conductivity; μS/cm, microsiemens per centimeter at 25 degrees Celsius; ft, feet; -- not 

recorded]

Site 

(fig. 5) Time WT (°C)

SC 

(μS/cm) Remarks

1 11:22 22.5 424 shaded, 6 ft offshore.

2 11:26 22.6 428 shaded, 1 ft depth.

3 11:29 20.0 415 shaded.

4 11:38 13.2 312 shaded, very shallow, with algae.

5 11:40 14.2 140 in 1 ft depth hole dug at shore.

6 11:43 11.3 477 active groundwater discharge at pond edge.

7 11:50 11.4 501 active groundwater discharge, 6 ft offshore.

8 11:51 11.3 502 active groundwater discharge.

9 11:55 19.1 435 shaded, 8 ft offshore.

10 11:56 12.0 383 active groundwater discharge in hole dug 1 ft deep at shore.

11 11:58 15.5 450 6 ft offshore.

12 12:01 22.0 436 4 ft offshore.

13 -- 28.8 435 north shore of pond.

435 Median.
Data 11

Table 3. Discharge measurements at the outflows from Pod 

Meadow Pond and Dudley Pond, east-central Massachusetts.

[ft3/s, cubic feet per second]

Date Discharge (ft3

/s)

Pod Meadow Pond

03-23-2011 1.18

05-09-2011 0.79

06-14-2011 0.65

07-12-2011 0.43

08-01-2011 0.54

09-12-2011 0.74

10-07-2011 0.45

11-10-2011 1

0.14

11-29-2011 1

0.12

12-14-2011 0.58

01-09-2012 0.86

02-08-2012 0.68

03-08-2012 0.75

04-10-2012 0.63

05-07-2012 0.70

06-01-2012 0.65

07-05-2012 0.79

08-09-2012 0.60

09-11-2012 0.49

10-05-2012 0.44

11-09-2012 0.71

12-07-2012 0.54

Mean 2

0.66

Standard deviation 20.17

1

Discharge affected by upstream beaver activity.

2

November 2011 measurements excluded because of upstream beaver 

activity.

Observation Wells

Although many boreholes exist in the study area, more 

lithological information was needed in a few locations to 

better define the extent of a clay layer under Lake Cochituate 

and to determine depth to bedrock. Boreholes were drilled at 

five sites and seven observation wells installed at four of these 

sites by Bristol Engineering Advisors, Inc. (Bristol) on the 

behalf of the Town of Framingham for use in this study (fig. 6 

and table 4). The sites were strategically located in accessible 

areas around the Birch Road wells where additional lithology 

and water-level information were determined to be of greatest 

value. Lithologic results from the boreholes were incorporated 

into model layering and representation of low hydraulic conductivity clay in the model, as discussed later (Groundwater 

Model Modifications section). 

A series of seasonal water-level measurements at the 

newly installed observation wells were made by Bristol on 

October 1, 2011; February 15, 2012; June 19, 2012; and 

August 21, 2012. Observation well SB-3 was not measured 

because it was not screened and was used only for determining the depth to bedrock. The minimum water levels were 

measured near the end of August and ranged from 0.69 to 

1.69 ft lower than the maximum water levels measured in 

October or in February. Water-level measurements made 

at these wells had an average standard deviation of 0.44 ft 

reflecting relatively small seasonal variations. Differences in 

water levels at paired wells finished in the upper and lower 

aquifer (SB-1, SB-4, and SB-5) were small at SB-1 and SB-5 

(averaged 0.45 and 0.16 ft, respectively), but were relatively 

large at SB-4 (averaged 27.2 ft higher in the upper aquifer than 

the lower aquifer). The high groundwater level in the shallow 

well at SB-4 is similar to high levels in nearby shallow wells 

and likely indicates a perched upper aquifer that could be 

caused by the presence of the low hydraulic conductivity clay 

layer present in this area. Because a perched aquifer would 

be separated by an unsaturated zone from the deeper aquifer, 

high groundwater levels observed in shallow wells in this area 

were not used as calibration targets for the model. Water-level 

measurements from the other newly installed wells were used 

for model calibration.

Groundwater Elevations

Groundwater elevations from 63 observation wells were 

used to calibrate the Phase 2 groundwater model. The number 

of water-level measurements at these wells ranged from 1 to 

16 with a median of 3. The measurements spanned from 1946 

through 2012, but the time span varied among the observation 

wells and typically was grouped in clusters, such as around a 

2006 aquifer pumping test. For calibration of a steady-state 

model meant to represent average conditions, it is important 

that groundwater levels used as calibration targets not be 

strongly affected by seasonal or climatic variations lasting 

1–2 years or less. Typically, a mean groundwater level is compiled for each site having multiple measurements; however, 

when there are few measurements or when the measurements 

are temporally clustered around a period of relatively high or 

low water levels, a simple average may not be representative 

of average conditions over 10 years or more.

In the Phase 1 model calibration, a simple average of the 

water-level measurements at each observation well was used 

for calibration. To minimize bias in Phase 2 head calibration 

targets, the measured groundwater levels were adjusted by 

referencing them to water level measurements at USGS well 

MA-WKW-2 (421852071220501) (fig. 1), which has monthly 

data from 1965 to 2010 and continuous data thereafter. Deviations of the daily water level in MA-WKW-2 (interpolated 
12 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

SB-5

SB-1

SB-3

SB-2

SB-4

Pod Meadow

Pond

Lake

Cochituate

Sudbury

Riveer

0 0.1 0.2 0.3 0.4 KILOMETER

Base from Massachusetts Geographic Information System Ortho Imagery 15-centimeter resolution 0 0.1 0.2 0.3 0.4 MILE

Universal Transverse Mercator projection, zone 19N

71°24' 71°23'

42°20’

42°19’30”

EXPLANATION

Location and number of borehole 

drilling site (table 6)

SB-5

Figure 6. Locations of borehole drilling sites and wells installed in 2011 in the Birch Road well area, east-central Massachusetts.

Table 4. Observation wells installed in 2011 in the Birch Road well aquifer area, east-central Massachusetts.

[Lat., latitude, in decimal degrees; Lon., longitude, in decimal degrees; elevation (approximate) in feet NAVD 88; feet, feet below land surface; 

--, not applicable]

Site 

(fig. 6) Lat. Lon.

Land 

surface 

elevation 

(feet)

Well 

depth 

(feet)

Top of 

bedrock 

(feet)

Well screen elevation (feet) Aquitard elevation 

(feet) Clay 

present

Shallow well Deep well

Top Bottom Top Bottom Top Bottom

SB–1 42.3314 -71.3847 125 61.5 60 21.2 26.2 51.5 55.5 -- -- no

SB–2 42.3265 -71.3845 186 176 176 50 55 -- -- 55 80 yes

SB–3 42.3310 -71.3853 125 45 38 -- -- -- -- -- -- no

SB–4 42.3254 -71.3946 180 124 124 64 69 99 104 70 90 yes

SB–5 42.3351 -71.3827 125 170 170 43 48 90 95 48 58 no
Groundwater Model Modifications 13

when only monthly data were available) from mean water 

level for the full period of record were used to adjust waterlevel measurements at observation wells used in the Phase 2 

study to better account for seasonal water-level conditions. 

 MA-WKW-2 is located just outside the active model 

boundary to the south and about 1,060 ft east of Lake Cochituate (fig. 1) and is finished in coarse grain material at about 

32 ft depth below land surface. The groundwater-level 

adjustment at observation wells used to calibrate the Phase 2 

steady-state model assumes water-level fluctuations at MAWKW-2 are representative of fluctuations in the observation 

wells used in the steady-state calibration. The validity of this 

assumption was tested by comparing water-level fluctuations 

at MA-WKW-2 with the few observation wells in the study 

area with sufficient data to compare against. Observation wells 

WKW-119, WKW-117, and F1W-84, finished in the upper 

aquifer had about 1 year of monthly data and generally followed the same pattern of water-level fluctuations observed in 

MA-WKW-2. Observation wells MW-4 and MW-5, finished 

in the lower aquifer had less than 1 year of data and had more 

sporadic observations. Water-level fluctuations in these wells 

did not match the fluctuation in MA-WKW-2 as closely as 

wells finished in the upper aquifer. Although the adjustment 

method helps water level observations better reflect appropriate mean values, the limited groundwater level measurement 

data, particularly in the lower aquifer, are still a source of 

model calibration error and uncertainty.

Eight observation wells used in the Phase 1 model calibration were not included in the Phase 2 calibration because 

these wells each had only one measurement made before 1978 

when water levels may have been affected by the then active 

Birch Road supply wells. In addition, other nearby observation wells were available that had more measurements when 

the Birch Road wells were inactive. The total number of wells 

used in the Phase 2 calibration included 57 wells used in the 

Phase 1 calibration, plus 6 of the new wells installed in 2011 

by Bristol. The revised water-level calibration dataset has an 

average of 6.4 measurements per well and just 7 wells with 

only 1 measurement.

Observation wells in the USGS NWIS database (USGS, 

2014) (18 of the 63 wells) used in this study are in the 

National Geodetic Vertical Datum of 1929 (NGVD 29) datum 

and were corrected to the North American Vertical Datum of 

1988 (NAVD 88) (-0.774 ft correction). With the exception 

of the seven recently installed wells (2011) and the observation wells in NWIS, the datum of other observation wells 

is unknown. Well installation documentation refers to the 

well elevation in feet above mean sea level only, which was 

assumed to be the NGVD 29 datum. To bring all well coordinates into the NAVD 88 datum a -0.774 ft correction was 

applied to groundwater-level measurements suspected of being 

in NGVD 29 datum. 

Groundwater Model Modifications

In the Phase 1 study (Eggleston and others, 2012), a 

groundwater model was developed using MODFLOW-NWT 

(Niswonger and others, 2011) to simulate groundwater flow, 

stream base flow, and the effects of groundwater pumping 

on surface water. MODFLOW-NWT is a Newton-Raphson 

formulation for MODFLOW-2005 (Harbaugh, 2005) to solve 

problems involving drying and rewetting nonlinearities of the 

unconfined groundwater flow equation (http://water.usgs.gov/

ogw/modflow-nwt/). The original groundwater flow model 

was discretized into 360 rows and 280 columns of cells 50 ft 

on a side (Eggleston and others, 2012, fig. 9) and discretized 

into 5 layers of variable thickness including explicitly representing bedrock as a layer (Eggleston and others, 2012, 

fig. 11). The original model was surrounded by no-flow 

boundary conditions (Eggleston and others, 2012, fig. 9). 

Model cells representing Dudley and Pod Meadow Ponds, and 

Lake Cochituate were modeled with high hydraulic conductivity. The original model was run in both steady-state and 

transient modes; the transient model had 60 monthly stress 

periods representing 5 years of average monthly conditions. 

The original model was calibrated by trial-and-error, with 

the resulting hydraulic parameter distribution described in 

Eggleston and others (2012). The original model had some 

issues with numerical instability and was sensitive to hydraulic 

conductivity under Lake Cochituate and around the pumping wells (Eggleston and others, 2012). Several issues such 

as these were identified in the original report that, if resolved, 

could help refine simulation of the area. This section describes 

implementation of changes to the Phase 1 groundwater flow 

model. Changes made to the groundwater model in Phase 2 

include the following:

• Combining Phase 1 model Layers 1 and 2 (surface and 

near surface layers) and deleting Layer 5 (bedrock 

layer);

• Representing Lake Cochituate with the MODFLOW 

lake (LAK) package;

• Removing Dudley Pond from the model;

• Inactivating cells with shallow depths to bedrock; and

• Adjusting model parameterization using automated 

methods.

Changes to the model reduced run times (by about a factor of 30) and improved model stability, allowing automated 

parameter estimation to be applied.
14 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Layering and Active Model Area

Model layering was simplified from the 5-layer model 

developed in Phase 1 to a 3-layer model in Phase 2 (fig. 7). 

Layers 1 and 2 of the Phase 1 model were combined to form 

Layer 1 of the Phase 2 model. In Phase 1, Layer 1 represented 

a thin (10 ft thick or less) unconsolidated sediment layer that 

became mostly dry during steady-state simulation (fig. 12A

in Eggleston and others, 2012). Because this layer had little 

effect on the simulations it was combined with the underlying 

Layer 2 to form the uppermost layer (Layer 1) in the Phase 2 

model, representing surficial sediment and underlying sand. 

Layer 2 of the Phase 2 model represents silt and clay deposits 

present in the north and south of the study area. Where buried 

silt and clay deposits do not exist in the study area, model 

Layer 2 was assigned a thickness of 0.5 ft. Layer 3 of the 

Phase 2 model represents the deeper sand and gravel aquifer in 

which the simulated supply wells are screened. Layer 5 in the 

Phase 1 model (top 80 ft of bedrock) was removed as it was 

determined to be relatively unimportant in model simulations 

relative to the more permeable unconsolidated deposits and 

little is known about the hydraulic properties of bedrock in the 

study area.

The active model area was reduced from about 6.1 mi2

 in 

the Phase 1 model to 5.5 mi2 in the Phase 2 model by inactivating cells, mostly along the eastern model boundary and in 

other areas where bedrock was close to the surface (fig. 1). 

A notable change was the small bedrock island (fig. 1) just 

north of the Birch Road area, which was made inactive in the 

Phase 2 model. Other structural characteristics of the Phase 2 

model remained unchanged relative to the Phase 1 model. 

These structural characteristics include total depth of the 

aquifer above bedrock and a uniform spatial discretization of 

2,500 ft2 cells (50 by 50 feet) with 360 rows (north–south) and 

280 columns east-west. 

Boundary Conditions and Surface-Water Body 

Representation

High horizontal hydraulic conductivity (Kh) zones in 

the Phase 1 model aquifer that represent Pod Meadow Pond 

and Lake Cochituate were removed in the Phase 2 model and 

replaced with surface boundary lakes using the Lake Package 

(LAK) in MODFLOW (Merritt and Konikow, 2000). LAK 

boundary cells control the flux between the water body and 

the top layer of the aquifer, providing greater flexibility in 

assignment of lakebed conductivity between open water and 

the aquifer. The depth of Pod Meadow Pond was set to 9.5 ft, 

whereas the depth of Lake Cochituate varied spatially according to bathymetric data. Inflows to Lake Cochituate and Pod 

Meadow Pond were assigned with the stream package (SFR), 

but precipitation to and evapotranspiration from the two lake 

surfaces were not simulated. Dudley Pond was removed from 

the Phase 2 model because the extensive thick gyttja deposits 

on the pond bottom limit it as a source of water to simulated 

pumping of the Birch Road wells. The simulated stream cells 

through Dudley Pond remain unchanged. 

Specified flow boundaries representing pumping at Wayland water supply wells (HH-1, HH-2, and MV-1 on fig. 1) 

were unchanged from the Phase 1 model. Long-term average 

steady-state or monthly pumping rates were assigned at all 

3 wells for all simulations with the Phase 2 model. 

Streams represented by the streamflow-routing package 

(SFR2; Niswonger and Prudic, 2005) were assigned to appropriate cells in Layer 1 (fig. 8). SFR2 cells traversing Lake 

Cochituate or Pod Meadow Pond were removed to accommodate the LAK package in the Phase 2 model except for two 

SFR2 cells needed to interface with the LAK cells. A SFR2 

cell was specified at the southern end of Lake Cochituate to 

receive inflow from SFR2 cells to the south and another SFR2 

cell was specified at the outlet of Lake Cochituate to receive 

output from the lake at the head of Cochituate Brook. Similarly, SFR2 cells interface to LAK cells to account for indirect 

runoff to Pod Meadow Pond and output water from the LAK 

cells at the outlet. All other SFR2 cells in the Phase 2 model 

remain unchanged from the Phase 1 model. 

Inflows were specified at stream reaches at the model 

boundary to account for upstream drainage area contributions 

to streamflow. For the steady-state model, inflows assigned 

to the Sudbury River, Lake Cochituate, and the tributary to 

Cochituate Brook were the same as those specified in the 

Phase 1 model; 107.5, 11.6, and 1.4 ft3

/s, respectively (fig. 8). 

These values are based on the median daily flow (144 ft3/s) 

observed at the Sudbury River at Saxonville (01098530) 

streamgage from January 1980 through December 2010 

proportionally distributed based on the drainage area upstream 

from each boundary and the Saxonville streamgage, less a 

small amount (23.5 ft3/s) to account for the intervening area 

between the streamgage and the model boundaries.

Transient simulations specified monthly inflows at 

the stream boundaries. For transient simulations of average monthly conditions, the monthly inflows were held 

the same as those specified in the Phase 1 model based on 

the 25th-percentile daily flows during each month. The 

25th-percentile flows are equivalent to the monthly 75-percent 

flow duration, which is the flow value that is exceeded 

75 percent of the time. The monthly 25th-percentile flows 

were computed from observed flows at the Sudbury River 

at Saxonville (01098530), from November 1979 through 

November 2011, and were apportioned by drainage area at the 

Sudbury River boundary (80.2 percent) and the Lake Cochituate boundary (17.8 percent), less a small percent of the drainage area accounting for by the active model area between the 

streamgage and the stream cell boundary cells.

For transient simulations of dry conditions, monthly 

inflows at the stream boundaries were specified as 

10th-percentile daily flows for each month (table 5). The 

10th-percentile flows are equivalent to the monthly 90-percent 

flow duration, which is the daily flow value exceeded 90 percent of the time, and are representative of drought condition, 

especially when sustained for more than a 12-month period. 
Groundwater Model Modifications 15 Figure 7. Vertical cross sections of groundwater model layers 

A, north-south and 

B, east-west.

Layer 1

Layer 3

Sudbury River

Lake Cochituate

Fines

Inactive 

Inactive 

Till and bedrock

Layer 2

WEST EAST

VERTICAL EXAGGERATION X 10

0 2,000 4,000 6,000 8,000 10,000 12,000 14,000

300200100

-100-200

0

B B'

Distance along model grid, in feet

Elevation above North American

Vertical Datum of 1988, in feet

Elevation above North American

Vertical Datum of 1988, in feet

Inactive model areaLayer 1—sand and gravelLayer 2—silt and clayLayer 3—sand and gravel

EXPLANATION

A. Modeled layers along column 125

B. Modeled layers along row 250

VERTICAL EXAGGERATION X 10

Inactive Fines Layer 2

Layer 1

Layer 3

Till and bedrock

0 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000

0

200100

-100

A'

A

Layer 1

Layer 3

300

NORTH SOUTH

Layer 2

Birch Road well area Cochituate Brook

Sudbury River Pod Meadow

SECTION

B–B'

SECTION

A–A'
16 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Figure 8. Boundary conditions in the Phase 2 groundwater model.

Dudley Pond

Lake Cochituate

(north pond)

Pod Meadow 

 Pond

Oxbow

Pod Meadow 

wetland area

SECTION B–B' ( (fig. 7B)

SECTION A–A' (fig. 7A)

Birch Road wells

#

# #

#

#

#

#

MV-1

HH-2

HH-1

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources, variously dated, various scales

Massachusetts State Plane Coordinate System, Mainland zone

71°24’ 71°22’

42°21’

42°19’

0 1,000 2,000 FEET

0 250 500 METERS

#

#

#

#

#

1

2

3

4

6

7

8 10

11

12

13

14

15

16

18

21

20

19

17

1.4

0.2

107.5

11.6

5

9

0 to 0.5

Bedrock 

island

Existing Birch Road 

test wells

Model grid cells 50 feet by 50 feet

#2 #3 #4 Pod Meadow Pond

#1

Alternate #3

Model column number

Model row number

40

60

80

100

120

140

160

180

200

220

240

260

280

300

320

340

360

20

0

200 40 60 80 100 120 140 160 180 200 220 240 260 280

EXPLANATION

01098500

#

#

#

1

Inactive model area

Lake package cells

River cell and segment number 

(table 7)

U.S. Geological Survey (USGS) 

streamgage and number

USGS partial record site

Wayland production well

and number

 Birch Road test well

Specified stream inflow 

and value, in cubic feet 

per second

11.6

Bedrock island

Birch Road test wells

HH-2
Groundwater Model Modifications 17

The 10th-percentile daily flows for each month observed at 

the Saxonville streamgage (January 1980 through December 2014) were apportioned again by applying 80.2 percent 

of the flow to the Sudbury River boundary and 17.8 percent 

of the flow to the Lake Cochituate boundary. The monthly 

10th-percentile flows decreased by 34 to 68 percent compared 

to the 75th-percentile flows with an average monthly decrease 

of 56 percent. There were no inflows to the tributary to Cochituate Brook for the dry-condition simulations.

Recharge rates in the Phase 2 model were assigned as 

specified fluxes to the top of the model at rates equal to the 

Phase 1 model, except for formerly active model areas that 

were made inactive in the Phase 2 model had recharge reassigned to the nearest active model cells in proportion to the 

area made inactive. For steady-state simulations, recharge 

was assigned a spatially uniform value of 0.00502 feet per 

day (ft/d) (22 inches per year [in/yr]). For the transient model, 

recharge rates were varied monthly at the same rates as in the 

Phase 1 model (table 6, Eggleston and others, 2012), including 

for transient simulations of dry conditions in which recharge 

was turned off during July, August, and September. Recharge 

was not calibrated or treated as an adjustable parameter during 

the model calibration process because the rate is relatively 

well known and constrained by previous studies (Zarriello and 

others, 2010; DeSimone, 2004; DeSimone and others, 2002).

Hydraulic Parameter Assignment

A major objective of the Phase 2 model was to improve 

the representation of aquifer hydraulic parameters. To achieve 

this objective, the Phase II model was calibrated by automated parameter estimation, which provides an optimal set of 

parameter values for the model design and observation dataset, 

in contrast to the Phase I study in which the model was manually calibrated. Calibrated hydraulic parameter values are 

described in this section of the report and compared to Phase 

1 model values from Eggleston and others (2012). Automated 

parameter estimation techniques used to calibrate the Phase 2 

model are also described in this section of the report. Model 

error and sensitivity are discussed later in the report (Model 

Uncertainty and Sensitivity Analysis section).

Horizontal hydraulic conductivity (Kh) was adjusted 

using automated parameter estimation techniques only within 

2,500 ft of the Birch Road test wells (fig. 8) because the many 

boreholes in that area reveal a complicated alluvial sediment 

structure and provide numerous water-level observations for 

detailed model calibration. Beyond 2,500 ft from the Birch 

Road test wells, calibration data are sparse and automated 

calibration of Kh outside the 2,500-ft radius tended to produce 

unreasonably high and low values, based on previous studies 

of similar hydrogeologic settings in Massachusetts (DeSimone 

and others, 2002; Masterson and others, 2009; Eggleston and 

others, 2012). Hence, Phase 2 model Kh values outside of the 

2,500-ft radius were assigned according to the conceptualized 

hydrogeology described in Phase 1 of the study. One change to 

the conceptualized hydrogeology was made by adding a zone 

of sand and gravel (fig. 9) around the Happy Hollow wells 

(HH-1 and HH-2 in fig. 9). This zone encompassed 672 cells 

where Kh was increased from 10 or 20 ft/d to 200 ft/d in 

layers 2 and 3. The Kh increase prevented unrealistically 

high simulated groundwater-level drawdowns in response to 

pumping at those wells and is consistent with descriptions of 

highly conductive sand and gravel deposits in that area from 

Fortin (1981). 

Table 5. Assigned monthly stream inflow to the model, in cubic feet per second, under average (25th percentile) and dry 

(10th percentile) climatic conditions.

Month

Average conditions Dry conditions

Lake Cochituate Sudbury River TotalLake Cochituate Sudbury River Total

Stream segment 1 Stream segment 5 Stream segment 1 Stream segment 5

January 24.8 103.2 128.0 14.5 60.5 75.0

February 26.8 111.9 138.8 17.0 71.0 88.0

March 42.2 175.8 218.0 22.7 94.9 117.6

April 37.9 158.1 196.0 23.8 99.3 123.1

May 23.0 95.8 118.8 13.9 58.1 72.0

June 10.8 45.2 56.0 7.4 30.7 38.1

July 5.6 23.4 29.0 3.3 13.7 17.0

August 4.0 16.7 20.8 1.9 7.8 9.7

September 3.5 14.5 18.0 1.7 7.3 9.0

October 6.8 28.2 35.0 2.3 9.7 12.0

November 14.5 60.5 75.0 8.3 34.7 43.0

December 23.2 96.8 120.0 14.2 59.4 73.6
18 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Figure 9. Horizontal hydraulic conductivity (Kh) of groundwater model for A, Layer 1, B, Layer 2, and C, Layer 3.

D D

D

D

D D

D

D

D

D

D

D

D

D D

D

D

D

D D

D

D

D D

D

D

D

D

D

D

D

D D

D

D D D D

D

D

D D

D

D

D

D

D

D

D

D D

D

D

D

EXPLANATION

>0.01 to 30

>30 to 60

>60 to 90

>90 to 120

>120 to 150

>150 to 180

>180 to 210

>210 to 240

>240 to 270

>270 to 300

Hydraulic conductivity, in feet per day

Observation well—Upper aquifer

D

Observation well—Lower aquifer

Wayland production well and number

Pilot point

A. Layer 1 B. Layer 2 

C. Layer 3 

0 0.25 0.5 1 MILE

0 0.50.25 1 KILOMETER

Lake

Cochituate

Dudley Pond

Pod Meadow

Pond

Heard

Pond

Sudbury River

Lake

Cochituate

Dudley Pond

Pod Meadow

Pond

Heard

Pond

Sudbury RiverSudbury River

Lake

Cochituate

Dudley Pond

Pod Meadow

Pond

Heard

Pond

Sudbury RiverSudbury River

HH-1

MV-1 HH-2 HH-1

MV-1 HH-2 HH-1

MV-1 HH-2 HH-1
Groundwater Model Modifications 19

Parameterization was augmented by use of PEST software (Doherty, 2010) to minimize model error by automated 

adjustment of selected model parameters. Parameters calibrated using PEST include horizontal and vertical hydraulic 

conductivity, streambed hydraulic conductivity, lakebed 

conductance, specific storage, and specific yield (table 6). The 

initial value of each model parameter was assigned based on 

Phase 1 values and then adjusted by PEST. PEST repeatedly 

runs MODFLOW-NWT using iteratively adjusted parameter 

values until the objective function is satisfactorily minimized 

and further parameter value changes are small based on userspecified convergence criteria (Doherty and Hunt, 2010). The 

final optimized parameter values (tables 6 and 7) were then 

used to run transient model scenarios discussed later (Effects 

of Pumping section).

A pilot point process (Doherty and others, 2010) was 

used to assign Kh values within a 2,500-ft radius of the 

pumping wells. Instead of assigning Kh by zone or by layer, 

Kh was interpolated from Kh values assigned at 54 pilot 

points, 18 in each layer (fig. 9). Pilot point locations were 

selected subjectively with the goals of locating points in areas 

of observation wells, keeping pilot points separated from one 

another, and keeping the number of pilot points low enough 

to allow reasonable run times. Each pilot point Kh value is 

treated by PEST as a separate parameter to be fitted based 

on how changes to that pilot point Kh value affect simulated 

heads and streamflows. The calibrated pilot point Kh values 

are spatially interpolated with kriging to assign Kh values to 

model cells within a 2,500-ft radius of the Birch Road test 

wells. Model cells outside of the 2,500-ft radius were assigned 

Kh values equal to Phase 1 values (except in the zone around 

the Happy Hollow wells previously mentioned). Model cells 

from 2,000 to 2,500 ft from the pump test wells were assigned 

a distance-weighted average Kh value between the kriged 

Table 6. Calibrated hydraulic parameter values in the Phase 2 

model, with comparison to Phase 1 values.

[Ss, specific storage in 1/feet; E, exponent (for example, 1.5E3 means “1.5 X 

103”); Sy, specific yield in fractional percent; hydraulic conductivity in feet 

per day; --, not applicable]

Parameter

Phase Model layers

Upper Middle Lower

Phase 1 1 and 2 3 4

Phase 2 1 2 3

Hydraulic conductivity1

Horizontal Phase 1 53.5 51.6 139.2

Phase 2 57.5 57.1 82.6

Vertical 

(expressed as 

anisotropy2

)

Phase 1 10.0 20.0 10.0

Phase 2 52.2 1.0 19.2

Storage

Ss Phase 1 1.00E-05 1.00E-05 1.00E-05

Phase 2 4.95E-07 4.20E-07 2.84E-06

Sy Phase 1 0.15 0.15 0.15

Phase 2 0.14 0.10 0.10

Lakebed leakance3

 (1/feet)

Pod Meadow Pond 0.05 -- --

Lake Cochituate 0.20 -- --

1

Horizontal hydraulic conductivity values are average for model cells 

within 2,500 feet of the Birch Road test wells.

2

Anisotropy is the ratio of horizontal to vertical hydraulic conductivity.

3

Lake (LAK) package used only in Phase 2 model.

Table 7. Calibrated streambed hydraulic conductivity values and 

comparison to Phase 1 values.

[Ks, vertical streambed hydraulic conductivity in feet per day; Ratio, Phase 2 

Ks/Phase 1 Ks; Sensitivity, change in PEST objective function divided by 

change in model parameter value]

Segment 

(fig. 8)

Phase 1 Phase 2 Ratio 

(percent)

Sensitivity 

(10–5) Ks Ks

1 20 15.2 76 0.95

2 20 15.2 76 0.95

3 20 15.2 76 0.95

4 20 15.2 76 0.95

5 20 15.2 76 0.95

6 20 15.2 76 0.95

7 20 15.2 76 0.95

8 20 2.7 14 250

9 10 175.4 1,754 27

10 20 175.4 877 27

11 20 2.7 14 250

12 20 2.7 14 250

13 20 2.7 14 250

14 20 2.7 14 250

15 20 2.7 14 250

16 200 0.5 0. 25 0.84

17 200 0.5 0. 25 0.84

18 200 0.5 0. 25 0.84

19 2 3.5 177 0.84

20 2 0.5 25 0.84

21 200 0.5 0.25 0.84

Mean 52.1 22.8 164 8.2
20 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

values and Phase 1 Kh values. The pilot point process was not 

applied to parameters other than Kh. 

PEST calibrated Kh values within 2,500 ft of the test 

wells indicate a complex and spatially heterogeneous aquifer system. Spatial variation of Kh created by the automated 

calibration process is greater and more detailed than the spatial 

variation of Kh in the Phase 1 model resulting from assignment by hydrogeologic zone. Within the 2,500-ft radius, calibrated Kh values in Layer 1 (upper aquifer) ranged from 0.1 

to 156 ft/d with an area of high Kh near the simulated Birch 

Road wells and an area of very low Kh values just south of the 

high Kh values (fig. 9). This contrasts with the corresponding 

Phase 1 Kh values in the upper aquifer (Phase 1 Layer 2) that 

were 60 ft/d except in wetland areas where they were 10 ft/d. 

Layer 2 (semiconfining unit) Kh values within 2,500 ft of the 

test wells are 30–40 ft/d immediately around the Birch Road 

wells and increase to as much as 160 ft/d just to the south of 

the test wells. This contrasts with the corresponding Phase 1 

Kh values in the semiconfining unit (Phase 1 Layer 3), which 

were assigned lower Kh values of 10 ft/d in wetland areas to 

the north of the test wells and in the area south of the test wells 

that may contain less permeable silts and clays and a value of 

70 ft/d elsewhere. In Layer 3 (lower aquifer) of the Phase 2 

model, PEST calibrated Kh values are as high as 300 ft/d to 

the south of the test wells but as low as 2 ft/d elsewhere. In 

contrast, corresponding Phase 1 Kh values in Layer 4 (lower 

aquifer) were a uniform 140 ft/d. The spatial variability of Kh 

determined by PEST within 2,500 ft of the test wells is supported by variability in the lithology of many well logs in that 

area. The spatial variability suggests that simplified conceptual 

models of aquifer geology may not be sufficiently detailed to 

accurately assign aquifer properties to a groundwater model of 

a glacial sediment aquifer with the complexity and heterogeneity seen in this study. 

Hydraulic conductivity anisotropy (ratio of horizontal/

vertical hydraulic conductivity), specific storage, and specific yield were calibrated using PEST but were assigned by 

layer rather than by cell and were reassigned everywhere in 

the model rather than just within 2,500 ft of the test wells 

(table 6). Vertical hydraulic conductivity for each cell in 

the Phase 2 model was determined by dividing the cell’s 

Kh value by an anisotropy factor assigned by layer. Phase 

2 model anisotropy values indicate greater variability and 

are higher in the upper and lower aquifer and lower in the 

middle semiconfining unit than Phase 1 values. Specific storage and specific yield are lower in the Phase 2 model than 

in the Phase 1 model. Lake-bed conductivity, a parameter 

not used in the Phase 1 model, was uniform for each of the 

two lakes and these values, 0.20 (1/day) for Lake Cochituate and 0.05 (1/day) for Pod Meadow Pond, were calibrated 

with PEST. 

Streambed hydraulic conductivity (Ks) was recalibrated 

in the Phase 2 model and values varied considerably from 

Phase 1 values (table 7). Stream segments were grouped for 

calibration according to morphology, observed bed characteristics, and position relative to the Birch Road well site. For 

example, stream segments 1–7 (locations shown in fig. 8) are 

upstream from the Birch Road site, and generally have greater 

slopes and gravel beds, whereas stream segments 11–15 are 

in the Pod Meadow area downstream from the Birch Road 

well site and generally have lower slopes and mud beds. The 

greatest increases in streambed conductivity were for stream 

segments 9 and 10, representing the inlet to and immediate 

downstream drainage from Pod Meadow Pond. Model results 

were most sensitive to streambed conductivity of the reaches 

in the Pod Meadow wetlands area (segments 8 and 11 to 

15). The greatest decrease in streambed conductivity was for 

downstream parts of the Sudbury River; however, no measurements of streambed hydraulic conductivity were available to 

guide parameter estimation and streambed hydraulic conductivity had low sensitivity in the model and little effect over 

model error.

Model Calibration

The Phase 2 groundwater model was calibrated by 

adjusting model input parameter values, with the use of PEST, 

to minimize the difference between simulated and observed 

water levels and stream base flow (table 8). For automated 

parameter estimation using PEST, the model was run with 

60 stress periods, the first a steady-state stress period providing simulated values for comparison to steady-state calibration targets (water-levels and outflow from Pod Meadow 

Pond). After the initial steady-state stress period the model 

simulated 39 monthly stress periods to establish dynamic 

seasonal equilibrium (monthly heads and flows consistent 

from year to year), then 20 daily stress periods to simulate the 

April 26 through May 15, 2006, pumping test. All observation 

data including steady-state water levels, Pod Meadow Pond 

outflow, and transient drawdowns were used collectively in 

the objective function, rather than calibrating the steady-state 

and transient models separately. Parameters were calibrated 

sequentially with PEST to achieve reasonable calibration run 

times. The sequential calibration process was repeated until a 

stable suite of calibrated values was achieved that provided the 

desired match between simulated values and observations.

Calibration Targets and Weights

Targets for model calibration were observations seasonally adjusted average groundwater levels, outflow from Pod 

Meadow Pond, just downstream from the pond, and measured 

drawdowns during an aquifer pumping test (table 8). The 

average groundwater level values differed somewhat from 

the values used to calibrate the Phase 1 model (Eggleston and 

others, 2012, tables 2 and 8), because they included additional 

data and were adjusted to remove seasonal bias, as described 

earlier (Groundwater Elevations section). To calibrate aquifer 

storage properties, drawdown data from an aquifer pumping 

test (SEA Consultants, Inc., 2008) were included as calibration 

targets (table 9).
Groundwater Model Modifications 21

Each observation target is assigned a weight that is 

used by PEST in the parameter estimation process. Weights 

are useful for emphasizing some observations or groups of 

observations more than others and for normalizing observations of different units or magnitudes. Weights were adjusted 

for different PEST runs during the calibration process as 

different sets of parameters were calibrated, but each set of 

weights was assigned following similar guidelines. The first 

guideline was to emphasize Pod Meadow Pond outflow as an 

observation target relative to groundwater elevation observations, because it was the only flux observation. The second 

guideline was to assign roughly equal weights to three defined 

groups of groundwater elevation observations; groundwater 

level observations in Layer 1 (ssgwlvlslay1), groundwater 

level observations in Layer 3 (ssgwlvlslay3), and pumping 

test drawdown observations (pumpgwlevels). These groups 

of observations served different purposes in the calibration: 

ssgwlvlslay1 observations had the strongest control over 

hydraulic characteristics of the upper aquifer and interactions 

between groundwater and surface water, ssgwlvlslay3 observations had the strongest control over hydraulic characteristics in the lower aquifer and flow to pumping wells, whereas 

pumpgwlevels observations had the strongest control over 

aquifer storage coefficients and the timing of delays between 

groundwater pumping and streamflow depletion. A further rule 

was to adjust observation weights individually by well within 

observation groups ssgwlvlslay1 and ssgwlvlslay3. The individual weights were assigned to reflect confidence in the average observed groundwater elevation value and the relatively 

lower importance of spatially clustered observations. Relative 

weight of an average groundwater elevation was based on the 

number of observations for a given well (50 percent of variation), proximity to other observation wells (30 percent), and 

other sources of uncertainty such as accuracy of XYZ coordinate locations (20 percent). Representative weights used in 

later PEST runs are shown in table 10.

Model Uncertainty and Sensitivity Analysis

Groundwater models are inherently uncertain because 

they are simplifications of complex natural systems. To better understand the effects of uncertainty on the soundness of 

conclusions drawn from model results it helps to understand 

which model parameters have the largest effects on simulation 

outcomes. For this Phase 2 model uncertainty and sensitivity 

are analyzed in four separate modes, as follows: 

1. PEST determined parameter sensitivities, 

2. numerical instabilities, 

3. model error, and 

4. hydraulic response to changes in hydraulic parameter 

values. 

PEST Parameter Sensitivities

Sensitivity to model parameters is calculated by PEST as 

part of the optimization process. Any given parameter sensitivity is measured as the weighted change in the objective 

function divided by the change in that model parameter value. 

Doherty and Hunt (2010) provides detailed descriptions of the 

calculation and meaning of model parameter sensitivities. 

PEST run times were prohibitively long when all parameters in the Phase 2 model were calibrated simultaneously and 

in some cases when particular parameter value sets caused 

difficulty for MODFLOW-NWT to reach a solution. To avoid 

such problems, parameters were calibrated sequentially with 

PEST; sensitive parameters were calibrated first, whereas less 

sensitive parameters were fixed at likely values. Then the more 

sensitive parameters were fixed at calibrated values, whereas 

less sensitive parameters were calibrated. This process was 

iterated until parameter values indicated little change between 

iterations and simulation errors fell below target maximums. 

Some parameters were grouped for PEST calibration to reduce 

PEST run times. Storage parameters (Ss and Sy) and vertical hydraulic conductivity were assigned by layer and stream 

segments were lumped into five groups representing different 

reaches for determination of streambed conductivity.

Grouping like parameters together, the parameter groups 

in order of greatest sensitivity were specific storage, specific 

yield, streambed hydraulic conductivity, pilot point horizontal 

hydraulic conductivity, aquifer vertical hydraulic conductivity, 

and lakebed conductance. Sensitivity values shown (table 11) 

are for one of the later PEST runs, prun23, with sensitivities 

of streambed hydraulic conductivity included from prun26 

because streambed hydraulic conductivity was not an adjustable parameter in prun23. 

Numerical Instabilities

In some cases the calibrated Phase 2 model did not 

converge if groundwater withdrawals were set too high or if 

recharge or aquifer hydraulic conductivity were set too low. 

In these cases the model convergence problems are usually 

caused by model cells that contain pumping wells going 

drying and rewetting from one solver iteration to the next, 

creating a numerical instability. MODFLOW-NWT automatically reduces pumping rates to prevent cells from drying, 

which can cause simulated pumping rates to be lower than the 

target rates. 

Alternative pumping rates and locations were therefore used to help the model converge during assessment of 

model sensitivities. Hypothetical locations of future supply 

wells in the Birch Road area were varied from actual locations to improve model convergence. By trial and error it was 

determined that moving the model location of Test Well 3 

away from the other test wells into an area of greater aquifer 

transmissivity improved model stability. The best location 

is labeled “alternate #3” on the figure 8 inset map. Rates of 

pumping at four hypothetical supply wells also were varied 

to provide additional stability in the model under drought 
22 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Table 8. Average observed and simulated groundwater levels and flow just downstream 

from the Pod Meadow Pond outlet for the calibrated groundwater model of an aquifer in 

east-central Massachusetts.

[ft3/d, cubic feet per day; ft3/s, cubic feet per second; %, percent; Average groundwater level applies to 

observed only; feet MSL, feet below mean sea level]

Site Observed Simulated Error = simulated − observed

Stream base flow (ft3

/d)

Pod Meadow Pond outflow 57,024 

(0.660 ft3

/s) 

 57,395 

(0.6643 ft3/s)

371 (0.7%)

Average groundwater level in wells (feet MSL)

F1W–64 124.2 120.4 -3.8

F1W–84 138.7 134.9 -3.8

F1W–92 126.3 128.7 2.4

WKW–119 131.9 130.8 -1.1

MW–10 120.0 120.5 0.6

MW–12 117.0 118.4 1.4

MW–15 126.2 125.1 -1.0

MW–16 129.4 124.1 -5.3

MW–8 119.4 117.0 -2.5

MW–9 116.7 117.4 0.7

USGS–F1S 122.1 122.8 0.7

USGS–F2 140.4 135.9 -4.5

USGS–F5S 117.4 120.5 3.2

1–90 124.6 126.4 1.8

2–90 124.7 126.3 1.6

3–90 126.4 126.2 -0.2

5–90 126.8 126.0 -0.8

7–90 119.1 120.6 1.4

8–90 126.4 128.0 1.7

F1W–74 137.4 136.5 -0.9

F1W–88 126.5 130.7 4.2

F1W–89 126.6 130.7 4.1

F1W–90 132.1 133.4 1.3

F1W–91 133.0 133.4 0.5

F1W–93 126.2 128.9 2.7

F1W–94 126.2 128.9 2.7

WKW–123 138.0 135.9 -2.1

WKW–030 118.0 119.1 1.0

WKW–052 117.2 115.8 -1.5

MW–01 126.1 125.1 -1.0

MW–11 123.5 121.2 -2.3
Groundwater Model Modifications 23

Site Observed Simulated Error = simulated − observed

Average groundwater level in wells (feet MSL)—Continued

MW–13 126.3 124.3 -2.0

MW–14 125.8 124.2 -1.6

MW–14R 126.1 124.2 -1.9

MW–15D 126.5 124.6 -1.9

MW–02 125.5 122.7 -2.8

MW–2B 125.8 122.7 -3.2

MW–2D 125.3 122.7 -2.6

MW–03 124.0 121.9 -2.0

MW–04 123.1 122.5 -0.5

MW–05 124.7 122.5 -2.2

MW–5D 125.4 122.5 -2.9

MW–06 131.4 125.2 -6.2

MW–07 129.9 122.6 -7.3

MW–8D 120.6 119.0 -1.6

MW–9D 118.2 118.6 0.5

SEA–10 127.1 129.5 2.5

SEA–11 128.7 128.1 -0.6

SEA–12 121.0 123.6 2.7

SEA–13 127.0 124.1 -2.8

SEA–14 121.0 124.3 3.3

SEA–15 116.8 118.4 1.6

SEA–16 130.8 133.0 2.2

SEA–17 127.0 125.6 -1.4

SEA–18 126.8 127.7 0.9

SEA–02 124.7 125.9 1.2

SEA–03 121.5 122.8 1.4

SEA–04 126.7 126.6 -0.1

SEA–07 127.0 125.7 -1.3

SEA–08 126.1 126.2 0.2

USGS–F1RD 121.8 122.9 1.1

USGS–F4D 127.1 125.1 -2.0

USGS–F5D 117.2 121.1 3.9

Mean 125.5 125.1 -0.4

Table 8. Average observed and simulated groundwater levels and flow just downstream 

from the Pod Meadow Pond outlet for the calibrated groundwater model of an aquifer in 

east-central Massachusetts.—Continued

[ft3/yr, cubic feet per year; ft3/s, cubic feet per second; %, percent; Average groundwater level applies to 

observed only; MSL, mean sea level]
24 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Table 10. Calibration target groups used in automated parameter estimation, with representative assigned weights by group.

[--, not applicable]

Calibration target group Description Number of 

observations

Observation weights

Average Range

streamflow Pod Meadow Pond outflow, steady state 

(cubic feet per day)

1 0.00004 --

ssgwlvlslay1 Groundwater level observation, steady-state, 

Layer 1 (feet)

14 0.02786 0.0166–0.333

ssgwlvlslay3 Groundwater level observation, steady-state, 

Layer 3 (feet)

50 0.01894 0.0096–0.032

pumpgwlevels Observed groundwater level drawdown 2006 

pumping test (feet)

30 0.12100 0.121–0.121

Table 9. Observed and simulated maximum groundwater-level drawdown during 2006 aquifer pumping test, Phase 2 

groundwater model of the aquifer in east-central Massachusetts.

Observation well

Groundwater-level drawdown at end of test (feet)

Observed Simulated Error = simulated – observed

7–90 1.3 2.4 1.2

MW–1 5.9 4.3 -1.7

SEA–10 6.1 3.9 -2.2

SEA–11 2.9 3.1 0.1

SEA–15 0.4 1.0 0.6

SEA–3 0.9 1.7 0.8

Average 2.9 2.7 -0.2

conditions and higher pumping rates, unlike in the Phase 1 

model where all four wells were assigned equal pumping 

rates (Eggleston and others, 2012). Phase 2 simulations using 

a variety of hypothetical well locations and pumping rates 

indicated that 3.17 Mgal/d of total pumping is close to the 

maximum possible steady-state withdrawal from four wells 

in the Birch Road area. Model cells are less likely to go dry 

if (1) pumping wells are spaced farther apart, (2) pumping 

wells are placed away from the bedrock island (fig. 8), and 

(3) pumping rates are higher at well locations having relatively 

greater transmissivity. Any siting of future supply wells should 

take these findings into account so that greater flexibility in 

pumping rates and schedules is achieved. 

One set of well locations was determined (locations #1, 

#2, #4, and alternate #3 in fig. 8) that allowed a steady-state 

pumping rate of 3.17 Mgal/d. Note that the supply wells are 

unlikely to be operated at this full maximum pumping rate on 

a steady-state or constant basis, which could affect site selection for future groundwater supply wells. Alternative site analysis was not exhaustive and, with realistic pumping schedules, 

other well configurations and pumping rates would likely yield 

3.17 Mgal/d. In addition, the drilling of wells for the pumping 

tests in 2006 (SEA Consultants, Inc., 2008) indicated that the 

aquifer contains localized pockets of cobbles and boulders 

with exceptionally high hydraulic conductivity. Although 

there is insufficient field data to accurately support inclusion 

of these localized high-Kh deposits in the groundwater model, 

locating wells in such deposits would reduce well drawdown 

during pumping that could yield pumping rates greater than 

that simulated. 

Model Error

Simulated steady-state heads show good agreement 

with observed heads (figs. 10 and 11). Errors are lower for 

the Phase 2 model that for the Phase 1 model. Phase 2 model 

simulated steady-state water levels have an average error of 

-0.4 ft (table 12), average absolute error of 2.1 ft, and range 

of -7.3 to 4.2 ft. Simulated maximum drawdowns in six wells 

during the 2006 aquifer pumping test have a mean error of 

-0.3 ft (table 12). Simulated Pod Meadow Pond outflow in 

the Phase 2 model matched the observed outflow (table 12). 

All of these error measures improve on Phase 1 model 

calibration results. 
Groundwater Model Modifications 25

The Phase 2 model, like the Phase 1 model, provides 

an acceptable spatial representation of average groundwater 

levels, average Pod Meadow Pond outflow, and pumping test 

drawdowns. The Phase 2 model, however, provides better 

representation of the spatial heterogeneity in aquifer hydraulic 

conductivity and more realistic simulation of lake boundary 

effects, which improved the model calibration and stability 

compared to the Phase 1 model. Overall, the Phase 2 model 

better reflects the heterogeneity of the aquifer and provides an 

improved basis for comparing model results and uncertainties.

Effect of Parameter Values on Streamflow Timing 

and on Induced Recharge from Lake Cochituate

Hydraulic parameters in the model affect the timing of 

streamflow response to pumping and rates of recharge from 

Lake Cochituate to the aquifer induced by pumping. The sensitivity of these hydraulic responses was tested by individually 

varying a selection of model parameters and calculating the 

resulting change in the timing of streamflow response and in 

the rate of induced recharge from Lake Cochituate. Model 

hydraulic parameters were varied for steady-state simulations 

to determine sensitivity of induced recharge and for transient 

simulations to determine sensitivity of streamflow response 

times (table 13). Hydraulic parameters varied included Ks, Kh, 

hydraulic conductivity anisotropy, lakebed leakance, specific 

yield, and specific storage. For simplicity, the model parameter 

changes were made uniformly throughout the model by factors 

of 50 and 200 percent. 

Sensitivity of induced recharge to parameter values was 

determined for the parameters listed in table 13 using steadystate simulations. Induced recharge is the simulated difference 

in net recharge from Lake Cochituate to the aquifer, caused by 

pumping at the Birch Road wells, as was done for the Phase 1 

model sensitivity tests (Eggleston and others, 2012). A steadystate pumping rate of 423,459 ft3

/d (equal to 3.17 Mgal/d 

or 4.9 ft3/s) was simulated at the existing well locations (#1 

through #4 in the inset map on fig. 8), to evaluate induced 

recharge from Lake Cochituate. Variation of induced recharge 

caused by parameter value changes was calculated relative to 

the calibrated model induced recharge rate of 152,807 ft3/d 

(SsCalb102 in table 13). Induced recharge from Lake Cochituate was equal to 36 percent of the simulated pumping, 

as compared to 32 percent in the Phase 1 model. Results 

(table 13) show that induced recharge from Lake Cochituate 

is most sensitive to Kh, which decreased induced recharge 

by 19.8 percent when Kh was decreased by 50 percent and 

increased induced recharge by 9.1 percent when Kh was 

doubled (SsCalb105 and SsCalb106). Independently halving 

and doubling streambed hydraulic conductivity (SsCalb103 

and SsCalb104) and anisotropy (SsCalb107 and SsCalb108) 

had little effect on induced lake recharge, and lakebed leakance (SsCalb113 and SsCalb114) caused only minor changes 

in induced recharge from Lake Cochituate (changes ranged 

from -3.3 to +1.8 percent). Lakebed leakance has little effect 

on rates of flow between lake and aquifer induced by pumping, increased lakebed leakance slightly decreases the induced 

recharge rate. This somewhat unexpected result is likely due 

to lakebed conductance affecting not only rates of flow from 

the lake but also rates of flow into the lake. Under conditions 

of no groundwater pumping, Lake Cochituate has a simulated 

net groundwater inflow, while under conditions of 3.17 Mgal/d 

continuous pumping at the Birch Road site, Lake Cochituate 

has simulated net groundwater outflow. In the scenario with 

higher lakebed leakance (SsCalb114), discharge from the aquifer to the lake and recharge from the lake to the aquifer both 

increase, but the net groundwater inflow to the lake increased 

more. The sensitivity of lake recharge to Kh indicates that the 

complexity of aquifer sediment textures, reflected in the well 

log information, has large potential to affect model accuracy 

and uncertainty. In addition, the variable lake-bottom sediments indicated in the marine seismic survey also affects the 

lake-aquifer interaction. 

The sensitivity of surface-water response time was 

calculated using monthly average transient model simulations 

(TransCalb210–TransCalb222). Response time was defined 

as the number of days for 50 percent of streamflow depletion 

(t50) to take place following 1 month of pumping at 2.3 ft3

/s 

(1.5 Mgal/d), distributed equally among the existing Birch 

Road wells (#1 through #4; inset map on fig. 8). A lower rate 

of pumping, 2.1 ft3/s (1.3 Mgal/d), was imposed for scenario 

TransCalb213 to enable model convergence. The transient 

model was run for 5 years (60 months) with no pumping, 

except in January of year 4 (month 37), which simulated 

pumping of 2.3 ft3/s for 1 month and then shut off for the 

remainder of the 5-year simulation. The transient simulations 

allowed the models with different parameter values to stabilize 

during the first 4 years then continued through month 60 to 

ensure streamflow depletion, measured at the Sudbury River 

exit from the model, was fully accounted for. Depletion was 

calculated by running each scenario with and without pumping and subtracting the streamflow simulated by the two runs. 

Because the transient simulations had monthly stress periods, 

linear interpolation was used to interpolate between monthly 

responses and estimate the number of days to reach the t50

target. For example, if total streamflow depletion at the end of 

pumping in January was 40 percent, and at the end of February was 60 percent t50 was calculated to take place halfway 

through February, 45 days (31 +14 days) after pumping 

started. Some variation in t50 could be expected if pumping 

were simulated in a different month, other than January of 

year 4 in which groundwater recharge is relatively high. For 

example, if pumping were simulated in August when recharge 

is at an annual low, then t50 would be longer; however, all 

simulations are relative to the same pumping period (January 

of year 4) to assess model sensitivities. 

Timing of streamflow response was most sensitive to 

specific yield (Sy). Halving and doubling Sy decreased t50

by 32 percent and increased t50 by 74 percent, respectively. 

The parameters Kh and anisotropy also affected the timing of 

streamflow response. Halving and doubling Kh increased t50
26 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Table 11. Parameter sensitivity values determined by PEST for the calibrated Phase 2 groundwater model of a shallow aquifer in 

east-central Massachusetts.

[E, exponent (for example, 1.5E3 means “1.5 X 103”); +, plus]

Parameter group Parameter

Name Relative sensitivity Name Description Sensitivity

Specific storage 1.0E+00 sslay1 Layer 1 6.2E-01

sslay2 Layer2 4.2E-02

sslay3 Layer 3 9.2E-01

Specific yield 1.4E-02 sylay1 Layer 1 1.7E-02

sylay2 Layer2 2.6E-04

sylay3 Layer 3 5.0E-03

Streambed hydraulic conductivity 1.1E-02 strcon1 Segments 1–7 9.2E-05

strcon2 Segments 8, 11–15 1.4E-02

strcon3 Segments 9–10 7.4E-05

strcon4 Segments 16–18, 20–21 2.3E-03

strcon5 Segments 19 9.3E-04

Horizontal hydraulic conductivity 

pilot points

1.8E-04 ptl1kx1 Layer 1, pilot point 1 7.6E-06

ptl1kx2 Layer 1, pilot point 2 4.8E-06

ptl1kx4 Layer 1, pilot point 3 5.0E-06

ptl1kx5 Layer 1, pilot point 4 5.1E-06

ptl1kx6 Layer 1, pilot point 5 7.4E-06

ptl1kx7 Layer 1, pilot point 6 1.1E-05

ptl1kx8 Layer 1, pilot point 7 4.9E-06

ptl1kx9 Layer 1, pilot point 8 4.9E-06

ptl1kx10 Layer 1, pilot point 9 4.9E-06

ptl1kx11 Layer 1, pilot point 10 4.9E-06

ptl1kx13 Layer 1, pilot point 11 4.8E-06

ptl1kx14 Layer 1, pilot point 12 4.8E-06

ptl1kx15 Layer 1, pilot point 13 4.9E-06

ptl1kx21 Layer 1, pilot point 14 4.9E-06

ptl1kx22 Layer 1, pilot point 15 5.5E-06

ptl1kx24 Layer 1, pilot point 16 5.1E-06

ptl1kx30 Layer 1, pilot point 17 5.3E-06

ptl1kx31 Layer 1, pilot point 18 4.8E-06

ptl2kx1 Layer 2, pilot point 1 4.8E-06

ptl2kx2 Layer 2, pilot point 2 4.8E-06

ptl2kx4 Layer 2, pilot point 3 5.0E-06

ptl2kx5 Layer 2, pilot point 4 4.8E-06

ptl2kx6 Layer 2, pilot point 5 4.8E-06

ptl2kx7 Layer 2, pilot point 6 4.8E-06

ptl2kx8 Layer 2, pilot point 7 4.8E-06

ptl2kx9 Layer 2, pilot point 8 4.8E-06

ptl2kx10 Layer 2, pilot point 9 4.8E-06

ptl2kx11 Layer 2, pilot point 10 5.0E-06
Groundwater Model Modifications 27

Parameter group Parameter

Name Relative sensitivity Name Description Sensitivity

Horizontal hydraulic conductivity

pilot points—Continued

1.8E-04 ptl2kx13 Layer 2, pilot point 11 4.8E-06

ptl2kx14 Layer 2, pilot point 12 4.8E-06

ptl2kx15 Layer 2, pilot point 13 4.8E-06

ptl2kx21 Layer 2, pilot point 14 4.8E-06

ptl2kx22 Layer 2, pilot point 15 4.8E-06

ptl2kx24 Layer 2, pilot point 16 4.8E-06

ptl2kx30 Layer 2, pilot point 17 4.8E-06

ptl2kx31 Layer 2, pilot point 18 4.8E-06

ptl3kx1 Layer 3, pilot point 1 6.2E-06

ptl3kx2 Layer 3, pilot point 2 6.7E-06

ptl3kx4 Layer 3, pilot point 3 6.8E-06

ptl3kx5 Layer 3, pilot point 4 6.7E-06

ptl3kx6 Layer 3, pilot point 5 5.8E-06

ptl3kx7 Layer 3, pilot point 6 5.8E-06

ptl3kx8 Layer 3, pilot point 7 6.4E-06

ptl3kx9 Layer 3, pilot point 8 5.4E-06

ptl3kx10 Layer 3, pilot point 9 5.2E-06

ptl3kx11 Layer 3, pilot point 10 5.4E-06

ptl3kx13 Layer 3, pilot point 11 5.2E-06

ptl3kx14 Layer 3, pilot point 12 5.0E-06

ptl3kx15 Layer 3, pilot point 13 5.5E-06

ptl3kx21 Layer 3, pilot point 14 5.0E-06

ptl3kx22 Layer 3, pilot point 15 5.2E-06

ptl3kx24 Layer 3, pilot point 16 5.3E-06

ptl3kx30 Layer 3, pilot point 17 5.3E-06

ptl3kx31 Layer 3, pilot point 18 5.3E-06

Vertical hydraulic conductivity 3.2E-05 vklay1 Layer 1 2.6E-05

vklay2 Layer2 1.4E-05

vklay3 Layer 3 9.7E-06

Lakebed conductance 1.4E-16 lakek_coch Lake Cochituate 2.2E-16

lakek_pmp Pod Meadow Pond 1.1E-18

Table 11. Parameter sensitivity values determined by PEST for the calibrated Phase 2 groundwater model of a shallow aquifer in 

east-central Massachusetts.—Continued

[E, power of 10 exponent; +, plus]
28 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Observed water level, in feet

Simulated water level, in feet

115

120

125

130

135

140

115 120 125 130 135 140

1:1 line

Figure 10. Observed water level compared to simulated steadystate groundwater levels (heads) for the calibrated Phase 2 

groundwater model.

by 58 percent and decreased t50 by 29 percent, respectively. 

Halving and doubling anisotropy increased t50 by 28 percent 

and decreased t50 by 14 percent, respectively. Small changes 

in t50 take place when streambed hydraulic conductivity (1.3 

to 4.5 percent) or lakebed leakance (-4.4 to 7.7 percent) are 

varied. Varying specific storage has no effect on t50. 

Induced recharge from Lake Cochituate was most sensitive to Kh of the aquifer, reducing Kh by one-half caused a 

38-percent decrease in induced recharge. Streamflow response 

time was most sensitive to changes in Sy. Doubling Sy caused 

t

50 to increase by 23 days, whereas a 50-percent decrease in Sy 

caused t50 to decrease by 10 days. The t50 streamflow response 

time of the Phase 2 model (31 days) was nearly four times 

greater than that of the Phase 1 model (8.3 days). 

Model Limitations

The Phase 2 groundwater model provides some notable 

changes and improvements over the Phase 1 model but, like all 

models, representation of the aquifer system is still imperfect. 

As such, users should be aware of the model limitations. The 

model simulates only parts of Lake Cochituate’s hydrology, 

most importantly the exchange of lake water with the aquifer 

and the effects of groundwater pumping on this exchange. 

The groundwater model does not simulate the surface runoff 

that strongly affects Lake Cochituate’s water budget and so, 

by itself, cannot accurately simulate lake levels. Any changes 

simulated are therefore best viewed as a relative change compared to other conditions simulated by the groundwater model. 

Similarly, the groundwater model simulates stream base flow 

but does not simulate surface runoff and so cannot accurately 

simulate streamflow. Therefore simulated changes in streamflow also should be viewed as a relative change compared to 

other conditions simulated by the groundwater model.

The model has limitations in its ability to simulate the 

exchange of water between the lakes and the aquifer. Although 

the new marine seismic data improve understanding of where 

Lake Cochituate’s bed is covered with gyttja, the low hydraulic conductivity mud that likely reduces flow, the hydraulic 

conductivity of Lake Cochituate bed sediments is still not 

well known. Simulated net recharge of groundwater from 

Lake Cochituate to the aquifer is not sensitive to the lakebed 

leakance parameter because leakance simultaneously controls 

groundwater discharge to the lake and groundwater recharge 

from the lake. Limitations in the model’s ability to simulate aquifer-lake interaction also apply to Pod Meadow and 

Pod Meadow Pond. Although the Pod Meadow Pond stream 

outflow is accurately simulated under existing (no Birch Road 

well pumping) conditions, it is not known if simulations of the 

pond under pumping conditions will be accurate and whether 

the model can accurately answer questions such as what rates 

of pumping might cause Pod Meadow Pond to go dry.

Using PEST to automatically calibrate hydraulic conductivity and other key parameter values reduced model 

error; however, model error still exists because of the many 

unknown properties of the aquifer and the limitations of 

the head and flow data available for model calibration. Few 

groundwater elevation measurements are available outside the 

2,500-ft radius around the Birch Road well site, so hydraulic conductivity and storativity are likely not well calibrated 

there. The resulting uncertainties restrict the applicability of 

detailed head simulations. Outflow from Pod Meadow Pond 

and groundwater fluxes in the Pod Meadow wetlands also are 

affected by parameter uncertainty. 

Other important considerations are the starting conditions 

and the time-steps of the model. Stresses of greatest concern 

often take place at a shorter time scale than a month, such 

as the 7-day low flow. Likewise, pumping schedules in the 

model were simulated on a monthly basis. The Phase 2 model 

as described in this report is not designed to simulate finer 

temporal or spatial resolutions that may be important to future 

day-to-day operation and management of withdrawals. 

Effects of Pumping

After calibration, the Phase 2 groundwater flow model 

was used to simulate the effects of pumping the Birch Road 

wells on groundwater flow and streamflow. Pumping scenarios (table 14) are based on design criteria for the proposed 

Birch Road wells that varied pumping rates as much as 

4.9 ft3/s (3.17 Mgal/d) (Peter Newton, Bristol, written commun., December 2014). All scenarios simulated the Wayland 

wells pumping at their average reported annual or monthly 
Effects of Pumping 29

90

130

130

120

140

150

160

170

120

140

150

160

130

130

140

150

160

170

170

120

Bedrock island

Lake Cochituate

(north pond)

Pod Meadow 

Pond

Sudbury River 

 Cochituate Brook

Sudbury RiverSudbury River

Wayland

Framingham

Natick

Sudbury

#

# #

MV–1

HH–2 HH–1

71°24’ 71°22’

42°21’

42°19’

Base from U.S. Geological Survey and Massachusetts Geographic Information System data, variously dated, various scales

Massachusetts State Plane Coordinate System, Mainland zone

0 1,000 2,000 FEET

0 250 500 METERS

120

Inactive model area

Town boundary

Simulated head contour 

and value—Contour 

interval 2 feet

General direction of 

groundwater flow

Wayland production well

and number

EXPLANATION

# MV–1

Simulated head error, in feet

>2.0

 0.5 to 2.0

-0.5 to 0.5

-2.0 to -0.5

<-2.0 

Figure 11. Simulated groundwater levels (heads) under existing conditions in model Layer 3 and error in simulated heads (observed 

minus simulated) in model Layers 1 and 3 for the study area in east-central Massachusetts.
30 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

withdrawal rates. References to scenarios with or without 

pumping only refer to pumping at the Birch Road wells. 

Notable changes from the Phase 1 model scenarios include 

(1) constant pumping throughout the year was reduced slightly 

(from 4.9 to 3.9 ft3/s) during July through September in order 

for the model to converge to a solution, (2) more conservative 

drought conditions were imposed, and (3) an additional set of 

scenarios simulated 10 percent return flow from the pumped 

wells to Pod Meadow Pond. 

Steady State Simulations

Water balances for steady-state simulations of existing conditions (SsCalb101) indicate that without pumping at 

the Birch Road wells, 9.4 ft3/s of the total 10.4 ft3/s inflow to 

the aquifer is from areal recharge from surface precipitation, 

whereas the remaining 1.0 ft3/s is recharge from streams and 

lakes (table 15). With hypothetical pumping rates of 4.9 ft3/s 

(3.17 Mgal/d) from the Birch Road wells (SsCalb102), the 

same amount of inflow recharges from surface precipitation 

(9.43 ft3/s), but an additional 1.81 ft3/s of recharge is induced 

from streams and lakes and groundwater discharge to streams 

and lakes decreases by 3.11 ft3/s (table 15). Net induced 

recharge from Lake Cochituate caused by pumping, equal to 

inflow change minus outflow change, is 36.1 percent of the 

pumping (1.8 of 4.9 ft3/s). The Phase 2 model water balance is 

similar to the Phase 1 model water balance, but in the Phase 2 

model lake-aquifer exchange is explicitly simulated.

An estimate of how future pumping at the Birch Roads 

wells may affect Lake Cochituate water levels can be made by 

reconsidering surface-water modeling results from Zarriello 

and others (2010), in light of the groundwater modeling results 

presented here. The 2010 study used HSPF software (Bicknell 

and others, 2000) to simulate surface water in the SudburyAssabet basin and the reaction of Lake Cochituate to hypothetical pumping at a constant rate of 6.65 ft3

/s (4.3 Mgal/d) at the 

Birch Road well site. Four surface-water modeling scenarios 

were simulated assuming different proportions of water to the 

Birch Road wells came from induced infiltration from Lake 

Cochituate (table 14 and fig. 33 of Zarriello and others, 2010). 

The 3rd scenario (BIR-LC33), which assumes that one-third 

of the water (to hypothetical Birch Road supply wells) is from 

Lake Cochituate and two-thirds is from the Sudbury River, 

is the scenario that most closely matches the 36.1 percent 

induced recharge indicated by the steady-state groundwater 

modeling results. The HSPF scenarios simulated the 44-year 

period 1961–2004, using an hourly time step and average 

water use conditions for 1993–2003. The HSPF model, unlike 

the MODFLOW-NWT model presented in this study, performs 

detailed accounting of surface-water inflows and lake outflows 

(surface outflow, evaporation). The HSPF simulation results 

indicate that under these conditions (BIR-LC33 scenario) Lake 

Cochituate water elevations drop below the Cochituate Brook 

outfall elevation 32.3 days per year as compared to 17.0 days 

per year under existing conditions without Birch Road pumping (AVGWU scenario). 

The number of days that Lake Cochituate water elevations are below the Cochituate Brook outfall simulated with 

the HSPF model would be less if the pumping rate decreased 

(4.9 compared to 6.65 ft3

/s), or ceased, or some combination 

thereof during the summer months. The MODFLOW-NWT 

groundwater flow model developed in this study could be used 

to refine the lake-groundwater interaction for application in the 

HSPF surface-water model and simulate the effect of alternate 

pumping schedules that would likely be used to maintain 

desired streamflows during the normal low-flow period. 

Coupling of the MODFLOW-NWT groundwater model to the 

existing HSPF surface-water model was beyond the scope of 

this study but a future study could address these questions. 

Simulation scenarios SsCalb101 and SsCalb102 also indicate that hypothetical steady-state pumping at the Birch Road 

well site also decrease flows and water levels in Pod Meadow 

Pond, the downstream Pod Meadow wetland, and the Sudbury River. The question of how pumping schedules could be 

managed to reduce such surface-water effects during low-flow 

periods, is addressed using the transient flow model.

Table 12. Summary of model errors.

[ft3/s, cubic feet per second; ft, feet; Error = simulated − observed]

Calibration targets—average error Observed Simulated Error

Pod Meadow Pond outflow (ft3

/s) 0.66 0.66 0.00

Mean steady-state heads1

 (ft) 125.5 125.1 -0.4

Mean pumping test drawdown2

 (ft) 2.9 2.6 -0.3

1

Mean of all observed and simulated water levels used in Phase 1 and Phase 2 steady-state model calibrations. Note 

the number and the observed wells used differs slightly. 

2

 Mean maximum drawdown at six observations simulated in Phase 1 and Phase 2 models (table 9 of this report).
Effects of Pumping 31

Table 13. Sensitivity of pumping-induced recharge from Lake Cochituate to the aquifer and of streamflow response time to select steady-state and transient model parameters.[ft3/d, cubic feet per day; T50, time for 50 percent change in streamflow depletion; --, not computed]

Parameter 

Parameter 

change factor 

(percent)

Induced recharge rate (steady-state) Streamflow response time (transient)

Model run Volume (ft3/d)

Change 

relative to 

calibrated 

(ft3/d)

Percent 

change Model run T50 (days)

Change 

relative to 

calibrated 

(days)

Percent 

change

Calibrated model

-- SsCalb102 152,807 -- -- TransCalb210 31.5 -- --

Sensitivity runs

Streambed hydraulic conductivity (ft/d) 50 SsCalb103 153,238 431 0.3 TransCalb211 31.9 0.4 1.3

200 SsCalb104 152,320 -487 -0.3 TransCalb212 32.4 0.9 3.0

Horizontal hydraulic conductivity (ft/d) 50 1SsCalb105 122,534 -30,273 -19.8 2TransCalb213 49.8 18.4 58

200 SsCalb106 166,721 13,914 9.1 TransCalb214 22.5 -9.0 -29

Anisotropy—horizontal/vertical 

conductivity

50 SsCalb107 151,818 -989 -0.6 TransCalb215 40.2 8.7 28

200 SsCalb108 153,261 454 0.3 TransCalb216 27.2 -4.3 -14

Specific yield (unitless) 50 -- -- -- -- TransCalb217 21.5 -10.0 -32

200 -- -- -- -- TransCalb218 54.7 23.3 74

Specific storage (1/ft) 50 -- -- -- -- TransCalb219 31.5 0.0 0.0

200 -- -- -- -- TransCalb220 31.5 0.0 0.0

Lakebed leakance (1/d) 50 SsCalb113 153,900 1,093 0.7 TransCalb221 33.9 2.4 7.7

200 SsCalb114 150,944 -1,863 -1.2 TransCalb222 29.8 -1.7 -5.3

1Pumping reduced 50 percent to permit model convergence, rates then rescaled 2X for comparison.2One month of pumping at 90,000 ft3/d, as compared to 100,000 ft3/d for the other transient scenarios.
32 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Transient Simulations

The model was used to simulate monthly responses of 

groundwater and surface water to hypothetical Birch Road 

pumping under average and low-recharge conditions. In 

the first set of transient scenarios (TransCalb400 to TransCalb404, table 14), average recharge rates were assumed and, 

as described earlier (Boundary Conditions and Surface-Water 

Body Representation section), streamflow input to the model 

was assigned 25-percentile flow values to represent ecologically conservative base-flow conditions. 

Simulated monthly water balances show seasonal hydrologic changes (fig. 12A). Under existing conditions in the 

warmer months of June through September, inflow to the aquifer is mostly from storage, as the water-table drops, whereas 

outflows are mostly groundwater discharge to streams. The 

colder months of October through March indicate inflows 

dominated by recharge from precipitation, whereas outflows 

are again mostly to streams with some outflow going into 

aquifer storage as groundwater levels rise. The small amount 

of pumping outflow under existing conditions is from the 

Wayland well withdrawals. 

Table 14. Groundwater model scenarios used to evaluate the aquifer and streamflow response to hypothetical withdrawals at the 

Birch Road wells, east-central Massachusetts.

[Yellow shaded cells highlight months of maximum pumping. Highlighted cells indicate pumping at 4.9 cubic feet per second (ft3/s) (3.17 million gallons per 

day [Mgal/d]) or 3.9 ft3/s (2.50 Mgal/d)]

Scenario Hypothetical pumping of Birch Road wells (cubic feet per second)

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

Steady-state model, average recharge conditions

SsCalb101 0

SsCalb102 4.9

Transient model, average recharge conditions

TransCalb400 0 0 0 0 0 0 0 0 0 0 0 0

TransCalb401 4.9 0 0 0 0 0 0 0 0 0 0 0

TransCalb402 4.9 4.9 4.9 4.9 4.9 4.9 3.9 3.9 3.9 4.9 4.9 4.9

TransCalb403 4.9 4.9 4.9 4.9 4.9 4.9 0.1 0.1 0.1 4.9 4.9 4.9

TransCalb404 4.9 4.9 4.9 0.1 0.1 0.1 0.1 0.1 0.1 4.9 4.9 4.9

Transient model, dry (low recharge) conditions

TransCalb410a 0 0 0 0 0 0 0 0 0 0 0 0

TransCalb411a 4.9 4.9 4.9 4.9 4.9 4.9 3.9 3.9 3.9 4.9 4.9 4.9

TransCalb412a 4.9 4.9 4.9 4.9 4.9 4.9 0.1 0.1 0.1 4.9 4.9 4.9

TransCalb413a 4.9 4.9 4.9 0.1 0.1 0.1 0.1 0.1 0.1 4.9 4.9 4.9

Transient model, dry (low recharge) conditions, 10 percent return flow1

TransCalb411b 4.9 4.9 4.9 4.9 4.9 4.9 3.9 3.9 3.9 4.9 4.9 4.9

TransCalb412b 4.9 4.9 4.9 4.9 4.9 4.9 0.1 0.1 0.1 4.9 4.9 4.9

TransCalb413b 4.9 4.9 4.9 0.1 0.1 0.1 0.1 0.1 0.1 4.9 4.9 4.9

1

Scenarios return 10 percent of groundwater withdrawals at the Birch Road wells to Pod Meadow Pond.

Table 15. Steady-state simulated water budgets with and 

without hypothetical Birch Road well withdrawals at a continuous 

rate of 4.9 ft3

/s, east-central Massachusetts.

[ft3/s, cubic feet per second]

Flow component Average existing 

conditions

With continuous 

Birch Road pumping 

of 4.9 ft3

/s 

Inflows

Recharge 9.43 9.43

Recharge from streams 0.53 0.58

Recharge from lakes 0.48 2.24

Total 10.44 12.26

Outflow

Pumping 1.40 6.30

Discharge to streams 7.87 5.19

Discharge to lakes 1.21 0.78

Total 10.47 12.26

Numerical error -0.03 0.00
Effects of Pumping 33

Figure 12. Simulated monthly A, flow components under existing average conditions and with additional pumping from Birch Road 

wells and B, changes in flows.

A. Simulated monthly flow components to and from the aquifer

Inflows to aquifer

Outflows from aquifer

E

P

E

P E

E P

P

E

P

E

P

E

P

E

P

E

P

E

P

E

P

E

P

E Existing conditions (scenario TransCalb400)

P Maximum pumping Birch Road wells (scenario TransCalb402)

EXPLANATION

Additional inflows

to aquifer

Additional outflows

from aquifer 

B. Change in flow components in response to additional pumping

Flow, in cubic feet per second Flow, in cubic feet per second

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

Precipitation

Streams

Lakes

Storage

(water table falling)

0

5

0

5

0

-20

-15

-10

-5

0

5

10

15

20

Streams

Lakes

Storage

(water table rising) (

Pumping

Month

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

Month

Lakes

Storage

Streams

Pumping

Storage
34 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

The monthly aquifer water balance changes appreciably when the Birch road wells are simulated as pumping at 

maximum hypothetical rates (scenario TransCalb402) of 4.9 

and 3.9 ft3/s (fig. 12B). The additional outflow is balanced 

mostly by greater net inflows to the aquifer from lakes and less 

outflow to streams (fig. 12B). 

Simulated streamflow in the Sudbury River is assessed at 

the downstream exit of the active model area where all pumping effects can be seen (fig. 8). Streamflow depletion is calculated by subtracting the no-pumping scenario streamflow from 

the pumping scenario streamflow. For example, streamflow 

depletion in response to 1 month of pumping is calculated by 

subtracting simulated Sudbury River streamflow in scenario 

TransCalb400 from simulated streamflow in scenario TransCalb401. The timing of streamflow depletion in response to 

pumping in the Phase 2 model is similar to timing in the Phase 

1 model (fig. 13). Streamflow depletion as a percent of pumping for scenario TransCalb401, during the month of pumping, 

is slightly less in the Phase 2 model than in the Phase 1 model 

(37 and 43 percent, respectively), but for all succeeding 

months streamflow depletion is slightly greater in the Phase 2 

model than in the Phase 1 model. Within 1 month after pumping has stopped 66 percent of all resulting streamflow depletion has taken place, within 2 months 82 percent of streamflow 

depletion has taken place, and within 4 months 94 percent of 

depletion has taken place. The rapid response of surface water 

to pumping stresses seen in both models indicates good potential for pumping to be managed so that streamflow depletions 

are acceptable during periods of low flow.

Simulations of dry climatic conditions (scenarios 

TransCalb410a to TransCalb413b) were made by decreasing 

recharge and streamflow specified at the model boundaries. 

As in the Phase 1 study, recharge was decreased by setting 

recharge to zero for the months of July, August, and September. In Phase 2 simulations, additional dry conditions were 

imposed by reducing streamflow inflow at the model boundaries from the monthly 25-percentile to the 10-percentile daily 

flow values as described earlier (Boundary Conditions and 

Surface-Water Body Representation section). The monthly 

10-percentile flows that are low and are sustained more than 

12 months represent drought conditions. The reduced stream 

inflow to the model, combined with the assigned recharge 

reductions for scenarios TransCalb410a-413b, impose very dry 

conditions that are more ecologically conservative than those 

used in the Phase 1 study.

The very dry conditions of scenarios TransCalb410a413b made it difficult for the model to reach a solution 

because many cells dried. Scenarios TransCalb411a and TransCalb411b were not able to solve due to model cells going dry. 

To estimate streamflows for the conditions of these two scenarios (table 14), results from scenario TransCalb402, which 

applies the same pumping conditions, were instead used and 

reduced by the difference between the 25 percent and 10 percent duration inflows. This alternative estimation enforces the 

water balance and gives values that should be close to simulated values had the model been able to converge to a solution. 

All of the dry condition scenarios, TransCalb410a-413b, had 

some drier months with water mass balance errors of more 

than 1 percent. Various model adjustments, such as shortening time steps and adjusting convergence criteria, were made 

to reduce errors as much as possible. The resulting simulated 

values are generally representative of expected conditions but 

should be understood to have limitations in accuracy.

In scenarios TransCalb411a through TransCalb413b 

withdrawals from the Birch Road wells are assigned at various schedules and rates (table 14) to evaluate the effects of 

pumping on surface-water flows. Simulation results provide 

examples of managing and reducing streamflow depletion during average recharge conditions (TransCalb400 through TransCalb404) and during dry conditions (TransCalb410a through 

TransCalb413b). Scenarios TransCalb400 and TransCalb410 

set a baseline condition with no pumping at the Birch Road 

wells under average and dry climatic conditions, respectively. 

For scenarios with pumping, pumping rates were assigned 

by month at one of several fixed rates (table 14); minimal 

pumping (0.1 ft3

/s), reduced-maximum pumping (3.9 ft3/s), or 

maximum pumping (4.9 ft3/s). The 0.1 ft3/s pumping rate is the 

minimum required to keep a groundwater treatment facility 

operational (Peter Newton, Bristol Engineering Advisors, Inc., 

written commun., December 2014). The reduced maximum 

pumping rate of 3.9 ft3

/s was the maximum rate that the model 

could sustain without convergence failure. The maximum 

pumping rate of 4.9 ft3/s corresponds to a proposed pumping 

rate in the Framingham permit application (SEA Consultants 

Inc., 2009). During actual operation of a future well field, 

Figure 13. Simulated streamflow response of the Sudbury River 

after 1 month of pumping the Birch Road wells at 4.9 cubic feet 

per second (3.17 million gallons per day).

0

20

40

60

80

100

0

0.5

1.0

1.5

2.0

2.5

0 1 2 3 4 5 6 7 8 9 10 11 12

Cumulative streamflow depletion, as percent of pumping

Streamflow depletion, in cubic feet per second

Month

EXPLANATION

Phase 1 streamflow depletion

Phase 2 streamflow depletion

Phase 1 percent depletion

Phase 2 percent depletion
Effects of Pumping 35

pumping rates could be more finely adjusted and adjusted at 

finer time scales (daily or hourly for example) in response to 

changing conditions, such as drought conditions or operation 

of upstream reservoirs on the Sudbury River by the MWRA, 

to maximize pumping while minimizing effects on nearby 

surface-water resources. 

Under average conditions, pumping causes relatively 

small changes to simulated stream base flows (table 16, 

fig. 14A). Changes in base flow range from about -3 percent 

during months of greatest recharge and streamflow to about 

-17 percent at the end of summer when recharge and streamflow were at seasonal lows (fig. 14B). The maximum streamflow depletion under average conditions is about 9 percent 

when pumping is for 9 months instead of 12 months and about 

5 percent when pumping is for 6 months instead of 12 months. 

The maximum streamflow depletion takes place at different 

times of the year in response to alternate pumping schedules.

Percent changes to streamflow are more pronounced 

under drought conditions (fig. 14C). The largest percent 

streamflow depletions are caused by pumping for 12 months 

at maximum rates under dry climatic conditions (TransCalb411b), which causes streamflow during September to 

decrease by about 21 percent from 11.8 to 9.3 ft3/s. Under 

these same conditions, if wells are pumped for 9 months 

rather than 12 months (TransCalb412b), streamflow depletion 

is decreased by about 13 percent, from 11.8 to 10.4 ft3/s and 

if wells are pumped for 6 months (TransCalb413b) streamflow depletion is decreased by just 3.0 percent from 11.8 to 

11.4 ft3/s (table 16, fig. 14A). 

Scenarios TransCalb411b–13b simulate the return of 

10 percent of groundwater withdrawals from the Birch Road 

wells to Pod Meadow Pond, to maintain flow to and from the 

pond. Equivalent scenarios were not simulated in the Phase 1 

study. In each month of pumping (12 months in TransCalb411, 

9 months in TransCalb412, and 6 months in TransCalb413) 

10 percent of groundwater withdrawals are diverted into Pond 

Meadow Pond (stream segment 9). These diversions increase 

streamflow in the stream that traverses the Pod Meadow 

wetlands (stream segments 9–12,14–15) and in the Sudbury 

River downstream from the Pod Meadow wetlands (stream 

segments 14–16, 18, 21). Streamflow is increased during the 

months of the diversions and also during the summer months 

when there are no diversions or pumping. At the Sudbury 

River exit from the model the gains from diverting 10 percent 

of the withdrawals to Pod Meadow Pond are small (fig. 14C), 

but account for a large percentage of the outflow from Pod 

Meadow without pumping. 

Table 16. Simulated streamflow in the Sudbury River at the exit of the Phase 2 groundwater model.

[Highlighted cells indicate pumping at the Birch Road wells at 4.9 or 3.9 cubic feet per second]

Scenario Sudbury River flow (cubic feet per second)

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

Steady-state model, average recharge conditions

SsCalb101 128.8

SsCalb102 123.8

Transient model, average recharge conditions

TransCalb400 137.0 148.6 227.4 206.7 128.9 64.0 35.4 25.5 22.4 42.0 82.2 129.6

TransCalb401 135.2 147.1 226.6 206.3 128.7 63.8 35.3 25.5 22.3 42.0 82.2 129.6

TransCalb402 132.0 143.5 222.3 201.6 124.0 59.3 31.2 21.7 18.7 37.7 77.5 124.5

TransCalb403 132.8 144.0 222.7 201.9 124.2 59.4 32.2 23.4 21.0 39.6 78.9 125.5

TransCalb404 132.9 144.2 222.8 203.1 126.6 62.7 34.7 25.0 22.0 40.1 79.2 125.8

Transient model, dry (low recharge) conditions

TransCalb410a 83.8 97.4 127.4 133.3 81.7 45.6 22.0 13.4 11.8 18.0 49.5 82.7

TransCalb411a 78.6 92.2 122.2 128.1 76.8 40.9 18.0 10.3 9.3 13.5 44.7 77.5

TransCalb412a 79.5 92.8 122.7 128.4 77.0 41.1 18.8 11.4 10.4 15.6 46.2 78.6

TransCalb413a 79.7 93.0 122.8 129.7 79.3 44.3 21.3 12.9 11.4 16.2 46.6 78.9

Transient model, dry (low recharge) conditions, 10 percent return flow1

TransCalb411b 79.2 92.8 122.8 128.7 77.3 41.4 18.3 10.5 9.4 13.7 45.2 78.0

TransCalb412b 79.9 93.3 123.2 129.0 77.5 41.5 19.2 11.6 10.5 15.8 46.5 79.1

TransCalb413b 80.1 93.4 123.3 130.1 79.5 44.5 21.3 13.0 11.5 16.3 46.8 79.3

1

10 percent of groundwater withdrawals at the Birch Road wells are returned to Pod Meadow Pond.
36 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Figure 14. Simulated A, monthly streamflows and streamflow depletions in the Sudbury River at the model exit, with 

hypothetical groundwater withdrawals at the Birch Road wells, under B, average, and C, dry conditions. Simulated streamflow, in cubic feet per second

10

100

1,000

A. Monthly streamflow 

Percent change in simulated streamflow

B. Streamflow depletion under average base flow conditions 

0

-5

-10

-15

-20

C. Streamflow depletion under dry conditions

0

-5

-10

-15

-30

-20

-25

Average base flow conditions,

scenarios TransCalb400, 402, 403, 404

Dry conditions,

scenarios TransCalb410a, 411a, 412a, 413a

Scenarios TransCalb413a-b

Scenarios

TransCalb412a-b

Scenarios 

TransCalb411a-b

Grey lines indicate scenarios in 

which 10 percent of withdrawals

are returned to Pod Meadow Pond

Scenario TransCalb404

Scenario TransCalb403

Scenario TransCalb402

No pumping

6-month pumping

12-month pumping

9-month pumping

EXPLANATION

Jan. Feb. Mar. Apr. May June July Aug. Sept. Oct. Nov. Dec.

Month
References Cited 37

Summary and Conclusions

The Town of Framingham has proposed increasing water 

supplies by reactivating the Birch Road well site near the 

Sudbury River and the adjacent towns of Wayland, Sudbury, 

and Natick. Questions have been raised about whether or not 

groundwater withdrawals from the Birch Road supply wells 

will affect hydrology of the ecologically important Sudbury 

River Basin in east-central Massachusetts; specifically whether 

or not withdrawals will reduce flow in streams, such as the 

Sudbury River, or adversely affect nearby State and Federal 

conservation areas, such as the Great Meadows National Wildlife Refuge, and nearby wetlands and surface-water bodies, 

particularly Lake Cochituate. To address these questions, the 

U.S. Geological Survey carried out this study in cooperation 

with the Town of Framingham. A Phase 2 groundwater flow 

simulation model was developed, based on the existing Phase 

1 model, and used to improve understanding of groundwater 

surface-water interactions.

In the previously completed Phase 1 study (Eggleston 

and others, 2012) a numerical groundwater flow model was 

developed (with MODFLOW–NWT) that revealed a need for 

additional data and model refinements to better understand 

the complex aquifer system and the effects of groundwater 

withdrawals on surface-water bodies and wetlands in the study 

area. In this Phase 2 part of the study, the original groundwater 

flow model was revised to improve representation of groundwater and surface-water hydrology, stabilize the model, and 

reduce model error. The model was simplified by reducing the 

number of layers from 5 to 3 and adding the MODFLOW lake 

package (LAK) to simulate Lake Cochituate and Pod Meadow 

Pond and better represent interaction between the lakes and 

the aquifer. Model revisions improved stability and shortened 

run times, allowing use of automated parameter estimation 

software (PEST) to further refine the model hydraulic parameters and reduce simulation errors.

The calibrated Phase 2 groundwater flow model was 

applied to simulate hydrologic conditions under hypothetical 

pumping and climatic conditions. Model simulations indicate that under average base-flow conditions, the Birch Road 

wells have a small effect on flow in the Sudbury River during 

most months, even at the maximum proposed pumping rate of 

4.9 ft3/s (3.17 Mgal/d). Maximum percent simulated streamflow depletion in the Sudbury River at the outlet of the active 

model area caused by pumping takes place during seasonally 

dry periods and during severe drought conditions, decreasing 

streamflow by as much as 17 percent under average conditions and as much as 21 percent under dry conditions. By 

adjusting groundwater withdrawal rates on a monthly basis, 

reductions in simulated streamflow under sustained drought 

conditions, the most ecologically conservative assumptions, 

streamflow depletion changed from about 21 percent under 12 

months of maximum pumping to 13 percent under 9 months 

of pumping and to 3 percent under 6 months of pumping. The 

maximum streamflow depletion took place in September under 

the maximum simulated pumping rates and seasonally low 

recharge and streamflow, but changed to other months of the 

year according to the pumping schedule simulated. Simulations that return 10 percent of the Birch Road well withdrawals to Pod Meadow Pond indicate a modest reduction in the 

Sudbury River streamflow depletion, but provide a larger 

percentage increase to streamflow at the pond outlet and in the 

stream reach in the downstream wetland.

Model simulations indicate that groundwater withdrawals at the Birch Road site could be managed to substantially 

reduce adverse streamflow impacts caused by groundwater 

withdrawals during seasonally dry conditions and severe 

drought conditions. Results of the study indicate that well 

locations can have a large effect on the sustainable pumping 

rate and so should be chosen carefully in any potential future 

reactivation of the well field. The simulation model is a useful 

tool for evaluating alternative pumping rates and schedules not 

included in this analysis.

References Cited

Bicknell, B.R., Imhoff, J.C., Kittle, J.L., Jr., Jobes, T.H., and 

Donigian, A.S., Jr., 2000, Hydrological simulation programFortran, HSPF version 12 user’s manual: Mountain View, 

Calif., AQUA TERRA Consultants, EPA Contract No. 68-C98-010, 845 p.

Clapp, F. G., 1904, Relations of gravel deposits in the northern part of glacial Lake Charles, Massachusetts: Journal of 

Geology, v. 12, no. 3, p. 198–214. 

DeSimone, L.A., Walter, D.A., Eggleston, J.R., and Nimiroski, 

M.T., 2002, Simulation of ground-water flow and evaluation of water-management alternatives in the upper Charles 

River Basin, eastern Massachusetts: U.S. Geological 

Survey Water-Resources Investigations Report 02–4234, 

93 p., accessed June 27, 2014, at http://pubs.er.usgs.gov/

publication/wri024234.

DeSimone, L.A., 2004, Simulation of ground-water flow and 

evaluation of water-management alternatives in the Assabet 

River Basin, eastern Massachusetts: U.S. Geological Survey 

Scientific Investigations Report 2004–5114, 133 p. [Also 

available at http://pubs.usgs.gov/sir/2004/5114/.]

Doherty, J.E., 2010, PEST, Model-independent parameter 

estimation—User manual (5th ed., with slight additions): 

Brisbane, Australia, Watermark Numerical Computing, 

accessed June 27, 2014, at http://www.pesthomepage.org/

getfiles.php?file=pestman.pdf.

Doherty, J.E., and Hunt, R.J., 2010, Approaches to highly 

parameterized inversion—A guide to using PEST for 

groundwater-model calibration: U.S. Geological Survey 

Scientific Investigations Report 2010–5169, 59 p.
38 Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.

Doherty, J.E., Fienen, M.N., and Hunt, R.L., 2010, Approaches 

to highly parameterized inversion—Pilot-point theory, 

guidelines, and research directions: U.S. Geological Survey 

Scientific Investigations Report 2010–5168, 36 p.

Eggleston, J.R., Carlson, C.S., Fairchild, G.M., and Zarriello, 

P.J., 2012, Simulation of groundwater and surface-water 

interaction and effects of pumping in a complex glacial-sediment aquifer, east-central Massachusetts: U.S. Geological 

Survey Scientific Investigations Report 2012–5172, 48 p. 

[Also available at http://pubs.usgs.gov/sir/2012/5172/.]

Fortin, R.L, 1981, The groundwater resources of Wayland, 

Massachusetts: Wayland Department of Public Works, 65 p.

Gay, F., 1985, Estimated water and nutrient inflows and outflows, Lake Cochituate, eastern Massachusetts, 1977–79: 

U.S. Geological Survey Water Resources Investigations 

Report 84–4315, 59 p. [Also available at http://pubs.er.usgs.

gov/publication/wri844315/.]

Harbaugh, A.W., 2005, MODFLOW-2005, The U.S. Geological Survey modular ground-water model - the GroundWater Flow Process: U.S. Geological Survey Techniques 

and Methods, book 6, chap. A16, variously paged. 

Masterson, J.P., Carlson, C.S., and Walter, D.A., 2009, 

Hydrogeology and simulation of groundwater flow in the 

Plymouth-Carver-Kingston-Duxbury aquifer system, southeastern Massachusetts: U.S. Geological Survey Scientific 

Investigations Report 2009–5063, 110 p. [Also available at 

http://pubs.usgs.gov/sir/2009/5063/.]

Merritt, M.L., and Konikow, L.F., 2000, Documentation of a 

computer program to simulate lake-aquifer interaction using 

the MODFLOW ground-water flow model and the MOC3D 

solute-transport model: Water-Resources Investigations 

Report 00–4167, 146 p.

Niswonger, R.G., and Prudic, D.E., 2005, Documentation 

of the Streamflow-Routing (SFR2) Package to include 

unsaturated flow beneath streams—A modification to SFR1: 

U.S. Geological Survey Techniques and Methods, book 6, 

chap. A13, 50 p. [Also available at http://pubs.er.usgs.gov/

publication/tm6A37/.]

Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW-NWT, A Newton formulation for MODFLOW-2005: 

U.S. Geological Survey Techniques and Methods, book 6, 

chap. A37, 44 p. [Also available at http://pubs.er.usgs.gov/

publication/tm6A37/.]

SEA Consultants, Inc., 2008, Source final report, Birch Road 

well re-activation, Framingham, Massachusetts: Cambridge, 

Mass., SEA Consultants, Inc., 73 p.

SEA Consultants, Inc., 2009, Final environmental impact 

report and notice of project change, EEA No. 14197, for 

Birch Road well site reactivation and water treatment plant, 

Town of Framingham: Cambridge, Mass., SEA Consultants, 

Inc., 201 p.

U.S. Congress, 1999, Sudbury, Assabet, and Concord Wild and 

Scenic River Act: PUBLIC LAW 106–20—APR. 9, 1999, 

5 p. [Also available at https://www.congress.gov/106/plaws/

publ20/PLAW-106publ20.pdf.]

U.S. Department of Interior, 2009, Comments on MEPA 

EOEA No. 14197, Birch Road Well Field Redevelopment 

and Water Treatment Plant: Newton, Mass., Anthony R. 

Conte Regional Solicitor, October 20, 2009, 4 p.

U.S. Geological Survey, 2014, National Water Information 

System (NWIS): accessed December 3, 2014, at http://nwis.

waterdata.usgs.gov/usa/nwis/.

Zarriello, P.J., Parker, G.W., Armstrong, D.S., and Carlson, 

C.S., 2010, Effects of water use and land use on streamflow and aquatic habitat in the Sudbury and Assabet river 

basins, Massachusetts: U.S. Geological Survey Scientific 

Investigations Report 2010–5042, 160 p. [Also available at 

http://pubs.usgs.gov/sir/2010/5042/.]
For more information concerning this report, contact:

Director, New England Water Science Center

U.S. Geological Survey 

10 Bearfoot Road 

Northborough, MA 01532

dc_nweng@usgs.gov

or visit our Web site at:

http://newengland.water.usgs.gov

Publishing support by: 

The Pembroke and Rolla Publishing Service Centers
ISSN 2328-0328 (online)

http://dx.doi.org/10.3133/sir20155174

Eggleston and others—Groundwater and Surface-Water Interaction and Effects of Pumping, Complex Glacial-Sediment Aquifer, Phase 2, Mass.—SIR 2015–5174