Prepared in cooperation with the Town of Framingham, Massachusetts

Simulation of Groundwater and Surface-Water 

Interaction and Effects of Pumping in a Complex 

Glacial-Sediment Aquifer, East Central Massachusetts

Land Surface

Bedrock Surface

100

300

200

-200

-100

0

Elevation, in feet above NGVD

West East

Layer 2

Layer 1

Layer 4

Layer 5

Sudbury River Lake Cochituate

Fines

Fines Sand and gravel

Fines

Bedrock and till

Layer 3

Sand and gravel

Inactive model

Vertical exaggeration X 10

Scientific Investigations Report 2012–5172

U.S. Department of the Interior

U.S. Geological Survey
Cover. Oblique view of the land topography and bedrock surfaces, and a cross section of the model layers between the surfaces.
Simulation of Groundwater and SurfaceWater Interaction and Effects of Pumping 

in a Complex Glacial-Sediment Aquifer, 

East Central Massachusetts

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

Prepared in cooperation with the Town of Framingham, Massachusetts

Scientific Investigations Report 2012–5172

U.S. Department of the Interior

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

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

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

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

To order this and other USGS information products, visit http://store.usgs.gov

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., 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/.)
iii

Acknowledgments

A stakeholder interest group assembled for the study provided feedback and guided selection 

of model scenarios that were examined by the U.S. Geological Survey. Members of the 

stakeholder interest group include representatives from the towns of Framingham and 

Wayland, Massachusetts Water Resources Authority (MWRA), Massachusetts Department 

of Environmental Protection and Department of Conservation and Recreation, Massachusetts 

Executive Office of Environmental Affairs, U.S. Fish and Wildlife Service, U.S. Environmental 

Protection Agency, National Park Service, Organization for the Assabet River, Water Supply 

Citizens Advisory Committee, and SuAsCo (Sudbury-Assabet-Concord) Watershed Community 

Council. A separate technical advisory committee provided scientific oversight and technical 

guidance to the study and included the following non-USGS members: David Ahlfeld (University 

of Massachusetts), Grant Garven (Tufts University), and Daniel Nvule (MWRA). The authors 

gratefully acknowledge the guidance from the stakeholders’ interest group and the technical 

advisory committee.

Rich Niswonger of the USGS provided valuable assistance in development of the MODFLOWNWT model, and Janet Stone and Byron Stone of the USGS substantially improved our 

understanding of glacial sedimentology at the site. The authors also extend special thanks to 

Peter Newton of Bristol Engineering Advisors, Inc., for his involvement in various aspects of the 

study ranging from sharing his knowledge of the study area to collecting field information used 

in the study.
THIS PAGE INTENTIONALLY LEFT BLANK
v

Contents

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

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

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

Purpose and Scope ..............................................................................................................................3

Study Area..............................................................................................................................................3

Previous Investigations........................................................................................................................3

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

Borehole Data........................................................................................................................................4

Geophysical Data..................................................................................................................................4

Bedrock ..................................................................................................................................................4

Glacial Sediment History .....................................................................................................................6

Water Resources ...........................................................................................................................................9

Surface-Water Features......................................................................................................................9

Groundwater Levels and Flow Paths...............................................................................................13

Water Use ............................................................................................................................................17

Groundwater Flow Model...........................................................................................................................18

Model Design.......................................................................................................................................18

Discretization..............................................................................................................................18

Model Boundaries .....................................................................................................................21

Hydraulic Conductivity..............................................................................................................21

Storage Coefficients..................................................................................................................22

Streams........................................................................................................................................22

Recharge.....................................................................................................................................25

Pumping.......................................................................................................................................25

Model Calibration................................................................................................................................25

Steady-State Calibration ..........................................................................................................25

Transient Model Calibration.....................................................................................................29

Model Uncertainty and Sensitivity Analysis ..................................................................................32

Model Limitations................................................................................................................................36

Simulated Aquifer and Streamflow Response........................................................................................37

Steady-State Simulations..................................................................................................................38

Transient Model Simulations ............................................................................................................38

Evaluation of Pumping Strategies to Reduce Low-Flow Stresses .............................................40

Summary and Conclusions.........................................................................................................................44

References Cited..........................................................................................................................................45
vi

Figures

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

2. Map showing well borings to bedrock, bedrock outcrops, passive seismic points, 

and seismic lines used to develop the hydrogeologic framework and bedrock 

surface topography, east central Massachusetts..................................................................5

3. Graphs showing passive seismic A, response resonance frequency at known 

bedrock depths and B, estimated and observed bedrock depths, east central 

Massachusetts..............................................................................................................................6

4. Map showing thickness of sediments above the bedrock surface in the study area, 

east central Massachusetts .......................................................................................................7

5. Map showing locations of glacial-sediment deposits in the study area, east central 

Massachusetts..............................................................................................................................8

6. Hydrogeologic cross sections running north-south (A–A′ and B–B′) and 

east-west (C–C′ and D–D′) through the study area, east central Massachusetts ..........10

7. Map showing surface-water features and interpolated groundwater table 

elevations in the study area, east central Massachusetts..................................................12

8. Photograph of Pod Meadow Pond showing open water along the southern shore in 

February 2011, Framingham, Massachusetts ........................................................................14

9. Map showing active model area and boundary conditions, east central 

Massachusetts............................................................................................................................19

10. Maps showing variability of hydraulic conductivity values of groundwater-model 

layers 1–4 used in the calibrated model, east central Massachusetts.............................20

11. Vertical cross sections of groundwater model layers and hydraulic conductivity 

values running A, north-south and B, east-west, east central Massachusetts ..............23

12. Maps showing model cells in layers 1–4 that dry in the calibrated steady-state 

model, east central Massachusetts ........................................................................................28

13. Graph showing observed and steady-state simulated groundwater levels, east 

central Massachusetts..............................................................................................................29

14. Map showing simulated steady-state groundwater levels and errors, east central 

Massachusetts............................................................................................................................30

15. Map showing wells used for transient model calibration of a glacial-sediment 

aquifer model in east central Massachusetts .......................................................................31

16. Simulated and observed groundwater hydrographs, under current conditions and 

no pumping at the Birch Road wells, east central Massachusetts ...................................33

17. Graph showing transient simulation of a 2006 aquifer test made with a daily time 

step compared to groundwater observations, east central Massachusetts ...................34

18. Map showing baseflow reduction (streamflow depletion) in stream-channel 

segments in response to simulated pumping at the Birch Road wells at a constant 

rate of 4.9 cubic feet per second (3.17 million gallons per day), east central 

Massachusetts............................................................................................................................39

19. Graph showing simulated transient monthly water budget under current average 

monthly conditions (no pumping of Birch Road wells), east central Massachusetts.....40

20. 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 (3.17 million gallons 

per day), east central Massachusetts ....................................................................................41

21. Graph showing simulated monthly streamflows in the Sudbury River at the model 

exit under five hypothetical pumping scenarios for the Birch Road wells, east 

central Massachusetts..............................................................................................................42
vii

22. Graph showing simulated streamflow along the Sudbury River for the month of 

September under average recharge rates and various pumping durations at the 

Birch Road wells, east central Massachusetts.....................................................................43

23. Graph showing percent reduction in Sudbury River simulated streamflow at the 

model exit in response to pumping under average and dry recharge conditions, 

east central Massachusetts .....................................................................................................43

Tables

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

Pond, east central Massachusetts..........................................................................................13

2. Borehole and groundwater observation wells used in the study, east central 

Massachusetts............................................................................................................................15

3. Average groundwater withdrawal rates from Wayland production wells, 1996–2000 

and 2002–2006, Wayland, Massachusetts..............................................................................17

4. Hydraulic conductivity of aquifer material and storage coefficients by layer in the 

calibrated model, east central Massachusetts.....................................................................22

5. Stream segments and properties in the groundwater flow model, east central 

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

6. Monthly recharge rates used in the transient groundwater flow model..........................25

7. MODFLOW-NWT input file variable values for NWT package for transient 

simulations of current conditions (scenario 4)......................................................................26

8. Hydrologic characteristics and observed and simulated water levels of Pod 

Meadow Pond, Dudley Pond, and Lake Cochituate, east central Massachusetts .........26

9. Simulated flows between Pod Meadow Pond, Dudley Pond, and Lake Cochituate 

and the aquifer, east central Massachusetts ........................................................................27

10. Observed and simulated water levels and drawdowns from aquifer test of the 

Birch Road wells, east central Massachusetts.....................................................................32

11. Model simulations used to determine sensitivity of induced recharge from 

Lake Cochituate and pumping response times to hydraulic parameter values, east 

central Massachusetts..............................................................................................................35

12. Changes in induced recharge from Lake Cochituate and streamflow response 

times resulting from changes in hydraulic parameter values, east central 

Massachusetts............................................................................................................................35

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

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

Massachusetts............................................................................................................................37

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

well withdrawals, east central Massachusetts ....................................................................38

15. Simulated monthly flow in the Sudbury River under different hypothetical pumping 

scenarios for the Birch Road wells, east central Massachusetts .....................................42
viii

Conversion Factors, Datum, and Abbreviations

Inch/Pound to SI

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)

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).

Abbreviations

EIR Environmental Impact Report

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

Ss Specific storage

Sy Specific yield

USFWS U.S. Fish and Wildlife Service

USEPA U.S. Environmental Protection Agency

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

and Effects of Pumping in a Complex Glacial-Sediment 

Aquifer, East Central Massachusetts

By John R. Eggleston, Carl S. Carlson, Gillian M. Fairchild, and Phillip J. Zarriello

Abstract

The effects of groundwater pumping on surface-water 

features were evaluated by use of a numerical groundwater 

model developed for a complex glacial-sediment aquifer 

in northeastern Framingham, Massachusetts, and parts of 

surrounding towns. The aquifer is composed of sand, gravel, 

silt, and clay glacial-fill sediments up to 270 feet thick 

over an irregular fractured bedrock surface. Surface-water 

bodies, including Cochituate Brook, the Sudbury River, Lake 

Cochituate, Dudley Pond, and adjoining wetlands, are in 

hydraulic connection with the aquifer and can be affected by 

groundwater withdrawals. 

Groundwater and surface-water interaction was simulated with MODFLOW-NWT under current conditions and 

a variety of hypothetical pumping conditions. Simulations 

of hypothetical pumping at reactivated water supply wells 

indicate that captured groundwater would decrease baseflow to 

the Sudbury River and induce recharge from Lake Cochituate. 

Under constant (steady-state) pumping, induced groundwater 

recharge from Lake Cochituate was equal to about 32 percent 

of the simulated pumping rate, and flow downstream in the 

Sudbury River decreased at the same rate as pumping. However, surface water responded quickly to pumping stresses. 

When pumping was simulated for 1 month and then stopped, 

streamflow depletions decreased by about 80 percent within 

2 months and by about 90 percent within about 4 months. The 

fast surface water response to groundwater pumping offers the 

potential to substantially reduce streamflow depletions during 

periods of low flow, which are of greatest concern to the ecological integrity of the river. Results indicate that streamflow 

depletion during September, typically the month of lowest 

flow, can be reduced by 29 percent by lowering the maximum pumping rates to near zero during September. Lowering 

pumping rates for 3 months (July through September) reduces 

streamflow depletion during September by 79 percent as compared to constant pumping. These results demonstrate that a 

seasonal or streamflow-based groundwater pumping schedule 

can reduce the effects of pumping during periods of low flow. 

Introduction

Glacial-sediment aquifers are an important source of 

water to communities in the northeastern United States. 

However, groundwater withdrawals from these aquifers are 

of growing concern because of their potential effects on 

surface-water resources, particularly streamflows. In the past, 

groundwater withdrawal limits were based on aquifer tests that 

determined the potential yield from wells with little regard to 

the effects of pumping on surface waters in hydraulic connection to the aquifer. More recent approaches to determining 

acceptable groundwater withdrawal rates include evaluating 

the long-term consequences of pumping on surface water 

resources and related ecosystems (Gleeson and others, 2011). 

The potential effects of groundwater withdrawals on surfacewater features recently gained attention when the Town of 

Framingham (fig. 1) sought to reactivate production wells (the 

Birch Road wells) that previously provided a local water supply from 1939 until about 1979. The wells were discontinued 

because of high iron and manganese concentrations and the 

availability of an alternative regional water supply. 

In 2009, in accordance with the Massachusetts 

Environmental Policy Act (MEPA), the Town of Framingham 

filed an Environmental Impact Report (EIR) to reactivate 

the wells (SEA Consultants, Inc., 2009). Letters of concern 

over the well reactivation were submitted during the MEPA 

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

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

Agency (USEPA), and the Office of the Solicitor for the U.S. 

Department of the Interior that specifically addresses Federal 

interests in the Sudbury River and the downstream Concord 

River corridors. Federal agency concerns were raised because 

the Sudbury River, which is near the proposed pumping wells, 

is designated by Congress as “Wild and Scenic,” requiring 

special resource protection. In addition, the river corridor 

includes the Great Meadows National Wildlife Refuge and 

the Minute Man National Historical Park. Of these, the 

Great Meadows National Wildlife Refuge is of particular 

concern because of its proximity to the Birch Road wells 

and its primary purpose to protect river and wetland habitats 
2 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

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

90

Sudbury River 

Sudbury River 

Cochituate Brook

Dudley Pond

Lake Cochituate

Heard Pond

Pod Meadow 

Pond

Great Meadows 

National

Wildlife 

Refuge

Bedrock island

WAYLAND

FRAMINGHAM

NATICK

SUDBURY

NATICK

SUDBURY

EXPLANATION

Wayland production well, 

 and identifier

Inactive model

 area

Wetland

Town boundary

Birch Road wells

Great Meadows 

 National Wildlife Refuge

Oxbow

Pod Meadow

Birch Road wells

MV-1

HH-2 HH-1

MV-1 0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET

0 250 500 METERS

42 o00'

42 o

30'

73 o

00'

71 o

Study Area 00'

MASSACHUSETTS BOSTON

0 50 KILOMETERS

0 50 MILES
Introduction 3

for migratory birds, fish, and other wildlife. In addition, 

environmental concerns were expressed in the MEPA review 

by state agencies, municipal agencies, and environmental 

interest groups, particularly about the effects of pumping on 

nearby Lake Cochituate, a state recreational resource.

The Sudbury River is considered stressed, particularly in 

its headwater reaches during low-flow periods, by development pressures and by numerous groundwater and surfacewater withdrawals that affect streamflows (Zarriello and 

others, 2010). Concern about potential impacts of additional 

groundwater withdrawals on surface-water resources and 

aquatic ecology has led to a need for better understanding 

of the hydrologic system in the area of the proposed pumping. The purpose of this study by the U.S. Geological Survey 

(USGS), made in cooperation with the Town of Framingham, is to contribute to that understanding by presenting a 

compilation of existing and new hydrogeologic data and a 

new groundwater simulation model that was used to assess 

interactions between groundwater and surface water. The 

purpose of the groundwater model is to simulate present and 

hypothetical pumping conditions and to assess the potential to 

manage groundwater pumping to minimize its effect on nearby 

surface-water features. Hydraulic connections between surface 

water and aquifers are a growing concern across the country, 

and the techniques used in this study contribute to balancing 

the competing demands of water-supply needs with naturalresource protection, particularly in glacial-sediment aquifers in 

the northeastern United States. 

Purpose and Scope

This report describes the development and calibration of 

a MODFLOW-NWT (Niswonger and others, 2011) groundwater model of the aquifer surrounding the Birch Road well 

site in Framingham, Massachusetts. The report also describes 

use of the model to simulate interaction between groundwater 

and surface-water features in the area under present conditions 

and scenarios of potential pumping from the Birch Road wells. 

The scenarios include alternative pumping rates and schedules that could reduce the effects of additional groundwater 

withdrawals on streams during low flows when the effects of 

withdrawals are most pronounced. 

Study Area

The study area is 16 miles (mi) west of Boston in east 

central Massachusetts (fig. 1). The active groundwater model 

area covers about 6.1 square miles (mi2) surrounding the 

Birch Road wells in the towns of Framingham (1.7 mi2), 

Wayland (4.0 mi2), Sudbury (0.4 mi2), and Natick (0.01 mi2).

The Sudbury River is the primary surface-water drainage 

feature, entering the study area from the southwest and 

discharging to the northeast. The Great Meadows National 

Wildlife Refuge (fig. 1) 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. Lake Cochituate is the largest 

water body in the study area and has a 17.5 mi2 watershed 

that is largely to the south and west of the study area. Lake 

Cochituate 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 meltwater deltaic deposits 

and proglacial lake deposits that range in texture from clay to 

coarse gravel and boulders. Most of these deposits are medium 

to fine sands, which are the primary aquifer deposits. 

Previous Investigations

The hydrology of the Lake Cochituate area was described 

by Gay (1985). Surficial geology of the study area was first 

described by Nelson (1974b, c) for the Framingham and 

Natick quadrangles and later updated by Stone and Stone 

(2006). Bedrock geology was described by Nelson (1974a, 

1975) and updated by Zen and others (1983). Pond-aquifer 

interaction for the South Pond part of Lake Cochituate (south 

of the area shown in fig. 1) was investigated by Friesz and 

Church (2001). A hydrologic watershed model of the Sudbury 

and Assabet River basins was developed by Zarriello and others (2010) to simulate effects of water use and land use. 

A variety of engineering consulting reports were 

published leading up to construction of the Massachusetts 

Water Resources Authority (MWRA) Metrowest watersupply tunnel that passes beneath the study area (Balsam 

Environmental Consultants 1986, 1987, and 1992). Seismic 

data were collected between 1989 and 1994 by GZA 

GeoEnvironmental, Inc. (1995). Aquifer characteristics, 

aquifer tests, and a one-layer groundwater model for the 

Birch Road well area were described by SEA Consultants, 

Inc., (1992, 2008) to evaluate reactivation of the Birch Road 

wells. Other observation wells and aquifer characteristics 

were described during investigations of nearby groundwater 

contamination (Sovereign Consulting, Inc., 2009; URS 

Corporation, 2003; IEP, Inc., 1983; Haley & Aldrich, 1996).

The Assabet-Sudbury River Basin study (Zarriello and 

others, 2010) includes end-member simulations of the effect 

of constant withdrawals from the Birch Road wells on river 

flow and Lake Cochituate water levels under assumed contribution or depletion of water from different surface-water 

sources. Hypothetical pumping was simulated at a constant 

rate of 6.65 cubic feet per second (ft3/s), with water coming 

from Lake Cochituate, the Sudbury River, or combinations 

of these sources to determine changes in lake level and river 

flow. Simulation results indicated that reactivation of the 

Birch Road wells could cause flows in the Sudbury River and 

Cochituate Brook to be “substantially affected during periods 

of low flow” (Zarriello and others, 2010). The authors note 

that “if pumping rates were varied, the effects on lake stage 
4 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

and flow could be mitigated if withdrawals were decreased at 

the appropriate times and by the appropriate amounts.” The 

current study looks at that possibility in greater detail.

Hydrogeology

Land surface elevations in the study area range from 

about 300 feet (ft) on the eastern margin to about 115 ft at 

the northeast outlet of the Sudbury River. Glacial-sediment 

deposits with complex stratification blanket a highly variable 

bedrock surface and fill a deep bedrock valley that underlies 

the north-south axis of Lake Cochituate. An understanding of 

the regional stratigraphic framework, geotechnical sediment 

classifications, and geologic depositional processes responsible for aquifer structure helps to appropriately represent the 

aquifer system in a groundwater flow model.

Borehole Data

Borehole data were compiled from consultant reports 

(Balsam Environmental Consultants, 1986, 1987, and 1992; 

GZA GeoEnvironmental, Inc., 1995; Haley & Aldrich, 1996; 

SEA Consultants, Inc., 1992, 2008; Sovereign Consulting, 

Inc., 2009; URS Corporation, 2003; and Bristol Engineering 

Advisors, Inc., 2011), USGS reports (Gay, 1981, 1985), 

and well construction reports kept at the USGS office in 

Northborough, Mass. Sedimentary logs and corresponding 

well construction records were available for 162 boreholes 

in the study area. These logs helped to establish details of the 

glacial sediment history and were used to define elevations 

and characteristics of hydrogeologic layers in the model.

Geophysical Data

Geophysical data were compiled from previous studies 

and were collected during this study (fig. 2) to establish depth 

to bedrock. Seismic data were collected between 1989 and 

1994 by the MWRA as part of a study for construction of a 

water-supply tunnel that passes beneath the study area (GZA 

GeoEnvironmental, Inc., 1995). Depths to bedrock data were 

determined by seismic refraction methods along the seismic 

lines shown in figure 2. 

Passive seismic methods were used to measure depth 

to bedrock at 32 sites (fig. 2) where no borehole data were 

available and at 7 calibration sites where depths to bedrock 

were known. Passive seismic technology uses ambient ground 

vibrations caused by ocean waves, rainfall, wind, and anthropogenic activities to determine the thickness of unconsolidated 

sediments overlying bedrock (Ibs-von Seht and Wohlenberg, 

1999; Lane and others, 2008). A three-component seismometer 

was used in this study to record the resonance frequency from 

ground vibrations, and a spectral analysis was made to obtain 

resonance frequencies related to the sediment thickness using 

equation 1.

Z afr

b = 0 , (1)

where

Z is the depth to bedrock at a location, in feet; 

f

r0 is the fundamental resonance frequency, in 

hertz; and

a and b are constants determined by a nonlinear 

regression of data acquired at sites with 

known depths to bedrock.

Values for a (359.29) and b (-1.1979) were determined 

from the data collected at the seven control sites in the study 

area with known depths to bedrock (fig. 3A). The depths 

computed from equation 1 at the control sites generally better 

matched depths at higher resonance frequencies (shallower 

depths to bedrock) than at sites with lower resonance frequencies (deeper depths to bedrock). Depths to bedrock at the 

deeper control points in the study therefore were computed 

from calibrated coefficients (a = 297.24 and b = 1.00) determined for Cape Cod (John Lane, U.S. Geological Survey, 

written commun., 2011), where depths to bedrock are generally deeper and more closely matched depths than those 

computed from locally derived coefficients. Depths to bedrock 

were therefore determined at sites with resonance frequencies 

greater than about 2.5 hertz (depths less than about 120 ft) 

from locally derived coefficients and resonance frequencies 

less than about 2.5 hertz (depths greater than 120 ft) from 

coefficients derived for Cape Cod (fig. 3B). Of the 32 passive 

seismic sites without known bedrock depths, 19 had depths 

greater than 120 ft, and 14 had depths less than 120 ft.

Bedrock

The study area is underlain by crystalline bedrock 

that crops out in places along the study area boundaries, 

predominantly on the east side, and at one location in the 

middle of the study area (fig. 4). To represent the bedrock 

surface as accurately as possible, the borehole and geophysical 

data described above were combined to produce a gridded 

bedrock surface elevation map. The bedrock elevation map 

was derived from compiled bedrock elevation data including 

(1) 120 boreholes, (2) 103 points from the MWRA seismic 

data profile lines, (3) 32 passive seismic points collected for 

this study, and (4) 2,383 control points. The control points 

were locations of bedrock outcrop and points from hand drawn 

contours of the deep bedrock valley by Janet R. Stone (U.S. 

Geological Survey, written commun., May 2011) identified 

to better constrain the automated interpolation of the bedrock 

surface. Some control points beyond the boundaries of the 

study area were used to improve surface interpolation at the 

edges of the study area where data were more sparse. The 

bedrock grid was set equal to the size of the groundwater 

model cells (50 square feet (ft2)) so that each grid value 

represented the elevation of the center of a model cell.
Hydrogeology 5

Figure 2. Well borings to bedrock, bedrock outcrops, passive seismic points, and seismic lines used to develop the 

hydrogeologic framework and bedrock surface topography, east central Massachusetts.

90

Sudbury River Sudbury River 

Sudbury River 

Cochituate Brook

Pod Meadow 

 Pond

Dudley Pond

Lake Cochituate

Heard Pond

Bedrock island

EXPLANATION

Wayland production 

 well, and identifier

Bedrock data points

Inactive

model area

Borehole to bedrock

Passive seismic point

Coincident passive seismic 

 and bedrock borehole point

MWRA seismic 

 profile line

Well with lithologic log

Wetland

Area of bedrock outcrop

Birch Road wells

Oxbow

Pod Meadow

Birch Road wells

WKW-2

MV-1

HH-2 HH-1

MV-1

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET 

0 250 500 METERS 
6 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Elevation of the bedrock surface is highly variable with 

a deep bedrock valley following the north-south axis of Lake 

Cochituate. Bedrock surface elevations range from outcrops at 

about 300 ft on the east side of the study area to -118 ft at the 

bottom of the bedrock valley under Lake Cochituate (fig. 4). 

An isolated bedrock outcrop referred to as the “bedrock 

island” rises to an elevation of 165 ft just north of the Birch 

Road wells and drops steeply eastward to an elevation of 

about -40 ft. The bedrock valley generally continues to rise 

to the north of the bedrock island, but few data are available 

to confirm the bedrock topography in this area. Bedrock 

topography likely determines groundwater flow patterns in the 

study area because bedrock is much less permeable than the 

overlying glacial deposits.

Nelson (1974a, 1975) and Goldsmith (1991) present two 

somewhat different interpretations of bedrock lithology in 

the study area. Nelson (1974a, 1975) indicated that bedrock 

underlying the study area is composed of three units: (1) the 

Cherry Brook Formation west of Lake Cochituate and west 

and north of Dudley Pond, (2) the Westboro Quartzite east of 

Dudley Pond, and (3) the Dedham Granodiorite east of Lake 

Cochituate. Goldsmith (1991) presented bedrock as two units: 

(1) quartzite and (2) metavolcanic rocks overlapping somewhat with Nelson’s units. All the bedrock is thought to be 

relatively impermeable, although the upper bedrock surface is 

considered more fractured and permeable because of weathering and glacial movement. The upper bedrock in other parts of 

Massachusetts has been shown to have an active groundwater 

flow system that serves as a source of water to wells (Boutt 

and others, 2010; Mabee and others, 2002; and Hanson and 

Simcox, 1994). Groundwater flow through the bedrock of the 

study area was documented during construction of a watersupply tunnel at depths of 200–500 ft below land surface. 

Groundwater flow rates were weakly correlated with observed 

lineaments (Mabee and others, 2002), and therefore lineaments were not represented in the model. In addition, groundwater levels in the glacial sediments and streamflow data do 

not indicate that the hydraulic properties of the bedrock affect 

shallow groundwater flow.

Glacial Sediment History

Overlying the bedrock in most of the study area are 

stratified glacial deposits laid down in the last stages of 

glacial Lake Charles during the retreat of the Wisconsin ice 

sheet (Clapp, 1904). A thin layer (generally 1–10 ft thick) of 

low-permeability glacial till lies immediately over bedrock 

in most areas. Thicker till deposits form several drumlin hills 

on the east and west sides of the study area (figs. 4, 5). Sediments blanket till and bedrock over most of the study area and 

are interpreted as glacial meltwater deltaic deposits (Stone 

and Stone, 2006), similar to those recognized and mapped in 

many places in New England (Koteff and Pessl, 1981). As the 

glacier retreated to the north and northwest, it periodically 

paused and deposited gravel and sand at its terminus. Three 

glacial sediment morphosequences (fig. 5), identified by areas 

of stratified sediments contained between landforms and ice 

margins, have been documented in the study area (Stone and 

Stone, 2006; Janet Stone, U.S. Geological Survey, written 

commun., 2011). Generally, meltwater sediment deposits were 

finer grained where they settled in glacial Lake Charles away 

from the ice margin and coarser grained where they settled 

near the ice margin or mouth of the meltwater streams. 

Figure 3. Passive seismic A, response resonance frequency at 

known bedrock depths and B, estimated and observed bedrock 

depths, east central Massachusetts.

y =359.29x -1. 1979

R2 =0.99

0

A. Resonance frequency

B. Estimated depth to bedrock

10 20 30

Resonance frequency, in Hertz

0

50

100

150

200

250

0 50 100 150 200 250

Estimated depth to bedrock, in feet

Depth to bedrock in boring, in feet Depth to bedrock in boring, in feet

Local

Bedrock depths based on

horizontal to vertical

spectral ratio (HVSR) for

CapeCod

Combined

0

50

100

150

200

250

EXPLANATION
Hydrogeology 7

Figure 4. Thickness of sediments above the bedrock surface in the study area, east central Massachusetts.

90

Sudbury River 

Cochituate Brook

Sudbury

River

Bedrock island

EXPLANATION

Sediment thickness, 

in feet 

0–10

>10–50

>50–100

>100–150

>150–200

>200

Inactive

model area

Lakes and ponds

Oxbow

Pod Meadow

Birch Road wells

Wayland production well

Birch Road wells

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET 

0 250 500 METERS 
8 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 5. Locations of glacial-sediment deposits in the study area, east central Massachusetts.

90

Sudbury River 

Cochituate Brook

Sudbury

River

Dudley Pond

Lake Cochituate

Pod Meadow 

 Pond

Heard Pond

-100

-50

-50

0

0

0

25

50

50

50

50

75

75

75

75

75

75

100

100

125

125

125

125

150

150

150

175

175

175

200

200

225

225

250

250

275

300

100

100

100

Bedrock island

EXPLANATION

 Morphosequence 1 

Undifferentiated fines

Undifferentiated coarse

Thick till

Fines at land surface

Fines at depth

Morphosequence 3

 Morphosequence 2

Morphosequence 1

A A Cross section (figure 6) '

50 Bedrock elevation, in feet

Wayland production well

Borehole or seismic point used in 

cross section and identifier

F1W-49 F1W-49

Inactive

model area

 Morphosequence 2 

 Morphosequence 3 

Inactive

model area

Row 250, E-E’ (figure 11)

Column 125, A’’-A’’’ (figure 11)

E E'

A''

A'''

C

C'

D

D'

A'

A

B

B'

Oxbow Pod Meadow

Birch Road wells

Birch Road well

F1W-60

Seismic-DSW

WKB-51

WKW-117

Seismic-LA WKW-27

F1W-49

F1W-90

F1W-88

SEA-10

8-90 F1W-92

WKW-119

WKW-121

Seismic-DP-4

Seismic-LRCGR

WKX-4

WKW-118

WKW-81

T-22C T-23

WKW-29

WKW-35

WKW-95

1-90

SEA-2 WKW-54

WKW-53

SEA-11

NS-21H NS-21I

SEA-14

SEA-13

MW-5D

MW-2D

NS-21E

NS-21J

NS-21F

T-19

T-18

Seismic-SLN

Seismic-SSN

F1W-60

Seismic-DSW

WKB-51

WKW-117

Seismic-LA WKW-27

F1W-49

F1W-90

F1W-88

SEA-10

8-90 F1W-92

WKW-119

WKW-121

Seismic-DP-4

Seismic-LRCGR

WKX-4

WKW-118

WKW-81

T-22C T-23

WKW-29

WKW-35

WKW-95

1-90

SEA-2 WKW-54

WKW-53

SEA-11

NS-21H NS-21I

SEA-14

SEA-13

MW-5D

MW-2D

NS-21E

NS-21J

NS-21F

T-19

T-18

Seismic-SLN

Seismic-SSN

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET 

0 250 500 METERS 

Modified from Stone and Stone, 2006

Bedrock elevation in NAVD 88
Water Resources 9

The sand and gravel deltaic deposits are interspersed with 

lower permeability fine-grained lacustrine deposits that are 

generally a mix of fine sands, silts, and clays (shown in light 

blue on fig. 5). Along geologic cross sections A–A′ and B–B′, 

extensive fine deposits are generally present below 140 ft in 

elevation in the southernmost parts of the sections and are 

mostly overlain by more permeable coarse-grain deposits 

except where they underlie Lake Cochituate (fig. 6A). Extensive fine deposits also are found in the northern part of the 

study area to the north of cross-section lines A–A′ and B–B′. 

Few borehole logs are available from this area to characterize 

these sediments, but those that do exist indicate that the fines 

consist of silty organic sediments, which may extend from the 

land surface to bedrock. The fine deposits in the middle part 

of the study area near the Birch Road wells are less extensive 

and slightly coarser (fine silts and silts) as indicated in crosssection lines C–C′ and D–D′ (fig. 6B). Contacts between sediments were determined from well logs and from a theoretical 

understanding of the morphosequences in which they were 

deposited. However, the positions of boundaries are highly 

variable and poorly known in most locations. 

Beneath kettle depressions such as Dudley Pond, Lake 

Cochituate, and Pod Meadow Pond, the sediments collapsed 

as the ice blocks beneath them melted. The high water-surface 

elevation of Dudley Pond (153 ft) relative to the elevations 

of nearby Lake Cochituate (138.5 ft) and Pod Meadow Pond 

(125.8 ft) is difficult to explain. Permeable sand and gravel 

deposits between the two sets of ponds should result in 

groundwater flow causing the level of Dudley Pond to drop 

closer to the level of Lake Cochituate. The reason these ponds 

can maintain such a high hydraulic gradient may lie at the bottom of the ponds. Sediment cores collected as part of an eutrophication study of Dudley Pond by the Town of Wayland (IEP, 

Inc., 1983) indicated a layer of bottom muck sediments up to 

14 ft thick. This muck layer, referred to as gyttja, is partially 

decayed organic material that settles out of the water column 

though time and has a black gel-like consistency. These deposits are an impediment to seepage losses from the pond to the 

aquifer and a likely explanation for why Dudley Pond exists 

and why the surface level does not substantially drop during 

the late summer when inflows to the pond are typically small. 

Similar deposits are believed to underlie parts of northern 

Lake Cochituate (dark blue in sections B–B′, C–C′, and D–D′, 

figs. 6A and B) as evidenced by up to 12 ft of gyttja deposits determined from ground-penetrating radar surveys of the 

South Pond of Lake Cochituate (Friesz and Church, 2001).

Water Resources

Surface water generally is in hydraulic connection with 

groundwater in the study area. Both surface water and groundwater are supplied by abundant precipitation, with average 

annual precipitation in 2004–09 measuring about 50 inches 

(in) at nearby Natick, Mass. (National Oceanic and Atmospheric Administration (NOAA) station USC00195175). 

Surface-Water Features

The major surface-water features in the study area 

include the Sudbury River to the west and north, the northernmost pond of Lake Cochituate to the south, Cochituate Brook 

to the southwest, Dudley Pond to the east, and Pod Meadow 

Pond (fig. 7). The shallow and permeable aquifer system is 

generally in close hydraulic connection with the abundant 

surface-water features in the study area, but the connection 

may be locally constrained by gyttja deposits in lakes and 

ponds as previously described.

Lake Cochituate consists of four ponds (only the northernmost is shown in fig. 7) connected by shallow, narrow 

waterways that form a relatively long south-north trending lake. Total drainage area at the lake outlet is 17.5 mi2

. 

The lake is a series of kettle ponds formed following the 

last glacial retreat. Of these four connected ponds, only the 

northernmost pond—hereafter called “Lake Cochituate”—is 

in the study area. This part of the lake has a 0.31 mi2 surface 

area and drains to the westward-flowing Cochituate Brook 

that connects to the northeastward-flowing Sudbury River. 

Streamflow in Cochituate Brook (fig. 7) was monitored by 

the USGS (01098500) from October 1977 to June 1979 and 

from August 2010 through June 2012. Daily mean flow for 

this entire period was 34.9 ft3/s. Lake Cochituate stage has 

been continuously monitored by the USGS (01098499) since 

August 2010, during which time the lake level varied by 2.1 ft 

through December 2011 and was lowest during parts of the 

summer when levels dropped below the crest of the outlet 

spillway. Flow and stage data collected by the USGS are 

maintained in the National Water Information System (NWIS), 

which is available at http://waterdata.usgs.gov/nwis. 

Other ponds in the study area include Dudley Pond 

northeast of Lake Cochituate, Pod Meadow Pond north of 

Lake Cochituate, and Heard Pond in the north-central part of 

the study area (fig. 7). Heard Pond was not explicitly modeled 

because of its distance from the Birch Road wells and location 

on the opposite side of the Sudbury River. Dudley Pond, also a 

kettle pond, has a surface area of 0.14 mi2 and a total drainage area of 0.58 mi2

 measured at its outlet. Pod Meadow Pond 

has a surface area of 0.01 mi2 and a drainage area of 0.23 mi2. 

Early topographic maps show Pod Meadow Pond as a wetland 

that was likely modified by sand and gravel excavations in the 

early to mid-20th century. 

Outflows from Pod Meadow Pond and Dudley Pond 

(fig. 7) were measured monthly starting in March 2011 to 

provide data for groundwater model calibration. The quality 

of streamflow measurement data was generally considered 

good to fair, but was poor (as defined by Kennedy, 1983) at 

the lowest flows because of low stream depth and velocity. 

Drainage areas at the outflow measurement sites to Pod 

Meadow Pond and Dudley Pond are 0.23 and 0.58 mi2, 

respectively; however, discharge measurements (table 1) 

averaged nearly the same at the two sites except for the 

November 2011 measurements at Pod Meadow Pond outlet, 

which were affected by a beaver dam that caused water to 
10 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 6. Hydrogeologic cross sections running north-south (A–A′ and B–B′) and east-west (C–C′ and D–D′) through the study area, east central Massachusetts.

220200

80

180160

140120100

604020

-20

0

Massachusetts Turnpike 

Gyttja

Gyttja Lake Cochituate

Dudley Pond 

B' B

South North

WKB-51

WKW-117

WKW-27

WKW-118

WKX-4

T-22C

WKW-35

WKW-95,

projected

Bedrock

Sand and gravel

Sand

Fines

Sand Fines

Till

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

Distance, in feet

Elevation, in feet above NAVD 88

A' A

Cochituate 

Brook

Bedrock

Boundaries between units are mostly uncertain

Till

Lake Cochituate level

Sand and gravel

Fine to very fine sand or silt

Sand

Sand

Seismic-LA

Seismic-DSW

F1W-60

F1W-49

F1W-90

F1W-88

SEA-10

8-90

1-90

SEA-2

220200

80

180160140120100

604020

0

0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000

Distance, in feet

Elevation, in feet above NAVD 88

South North

Fines

Fine to very fine sand or silt

EXPLANATION

Boundaries between most contacts are uncertain

Fines (silt and clay)

Fine to very fine sand or siltSand

Sand and gravel

Till

BedrockKettle depression

Vertical exaggeration X 20

Vertical exaggeration X 20

Gyttja (organic muck)

90

Geologic interpretation assisted by Janet A. Stone, 2011

SECTION D–D’

SECTION C–C’

SECTION D–D’

SECTION C–C’
Water Resources 11

Figure 6. Hydrogeologic cross sections running north-south (A–A′ and B–B′) and east-west (C–C′ and D–D′) through 

the study area, east central Massachusetts.—Continued

Screened interval of SEA test wells,

800 feet to north

Artificial fill

0

200

80

180

160

140

120

100

60

40

20

-20

-40

-60

220

240

-80

200

80

180

160

140

120

100

60

40

20

0

-20

-40

-60

220

Dudley

Pond 

Birch Road wells screened 

interval 400 feet to south

Elevation, in feet above NAVD 88 Elevation, in feet above NAVD 88

Sudbury

River 

Pod Meadow 

Pond

C C'

West East

D D'

West East

T-18

T-19

MW-5D

MW-2D

NS-21F

NS-21H

Seismic-SSN

Seismic-SLN

NS-21I

SEA-11

SEA-10

F1W-92

WKW-119

WKW-121

Seismic-LRCGR

Seismic-DP4

NS-21J

SEA-13

SEA-14

WKW-54

T-22C

T-23

WKW-53

NS-21E

Lake Cochituate level

Sudbury

River 

Lake Cochituate

Dudley

Pond 

WKW-81

Boundaries between units are mostly uncertain

Boundaries between units are mostly uncertain

Bedrock

Till

Sand and gravel

Sand

Bedrock

Till

Sand

Sand and gravel

Sand

Gyttja

Gyttja

Gyttja

Gyttja

0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000

0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000

Distance, in feet

Distance, in feet

EXPLANATION

Fine to very fine

 sand or silt

Sand

Sand and gravel

Till

Bedrock

Gyttja (organic muck)

Kettle depression

Vertical exaggeration X 20

Vertical exaggeration X 20

Geologic interpretation assisted 

by Janet A. Stone, 2011

SECTION A–A’

SECTION B–B’

SECTION A–A’

SECTION B–B’
12 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 7. Surface-water features and interpolated groundwater table elevations in the study area, east central 

Massachusetts.

120

120

120

130

130

130

130

140

140

140

150

150

150

160

160

160

90

Sudbury River Sudbury River 

Cochituate Brook

Sudbury

River

Heard Pond

Dudley Pond

Lake Cochituate

Pod Meadow 

Pond

Bedrock island

EXPLANATION

Groundwater-table 

 elevation (approximate) 

General direction 

 of groundwater flow 

Wetland

120

Observation well used to draw water table,

 and identifier

Wayland production well, and identifier

USGS streamgage 

 and number

Inactive

 model area

Birch Road wells

Pod Meadow Pond outflow

01098500

01098500

USGS partial-record streamgage 

01098530

Dudley Pond outflow

WKW-119

F1W-84

MW-8

MW-1

MW-1

F1W-84

Oxbow

Pod Meadow

Birch Road 

wells

MV-1 HH-2

HH-1

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET

0 250 500 METERS

Groundwater elevation in NAVD 88
Water Resources 13

pool upstream of the measurement site. Discharge from Pod 

Meadow Pond also was relatively consistent, ranging from 

0.43 to 1.18 ft3/s, compared to the outflow from Dudley Pond, 

which ranged from 0 to 1.71 ft3/s. The relatively high and 

consistent outflow from Pod Meadow Pond, which has a small 

contributing area relative to Dudley Pond, suggests a larger 

groundwater discharge to Pod Meadow Pond than to Dudley 

Pond. This higher outflow is also evident from anecdotal 

reports that the southwestern part of Pod Meadow Pond does 

not freeze, which was confirmed during site visits in the winter 

of 2010–11 (fig. 8) when both Lake Cochituate and Dudley 

Pond were completely frozen and part of Pod Meadow Pond 

was not, indicating discharge of relatively warm groundwater. 

Factors that could contribute to the groundwater discharge 

at Pod Meadow Pond include its low topographic position, 

removal of surface material and lowered topography as a result 

of past mining activities, occurrence of high permeability sand 

and gravel, proximity to Lake Cochituate and related steep 

hydraulic gradient to the lake, increasing bedrock elevation 

forcing groundwater flow upward, and the pond’s position just 

north of the edge of the extensive layer of low-permeability 

fine-grained sediments.

Extensive wetlands exist adjacent to the Sudbury River 

and other streams in the northern part of the study area. Pod 

Meadow Pond drains north into the Pod Meadow wetland, 

which then drains into the oxbow on the Sudbury River. Just 

north of the oxbow, the Sudbury River meanders through 

extensive wetlands that are part of the Great Meadows 

National Wildlife Refuge (fig. 1).

The Sudbury River flows from the southwest toward the 

northeast through the study area and is the primary drainage 

feature that likely receives all groundwater discharge from 

the study area aquifer, either directly or indirectly from 

numerous tributaries. The drainage area of the Sudbury River 

is about 85 mi2 at its entrance to the study area and about 

111 mi2 at its exit from the study area. Groundwater levels 

indicate that there is little or no groundwater crossing beneath 

the river. Flow in the Sudbury River near the southwestern 

boundary of the study area has been continuously monitored at 

Saxonville (01098530) since October 1979. Mean daily flow 

at Saxonville from October 1979 through September 2009 was 

205 ft3/s with monthly mean flows ranging from 71 ft3/s in 

September to 384 ft3/s in April.

Groundwater Levels and Flow Paths

Groundwater level measurements collected in and near 

the study area for a variety of purposes were compiled for 

this study. Groundwater levels were obtained from NWIS or 

compiled from previous site investigations (Balsam, 1987, 

1992; SEA, 1992, 2008). Additional groundwater observations 

were also made during this study (Peter Newton, Bristol 

Engineering Advisors, Inc., written commun., 2011). The 

complete set of groundwater level observations (table 2) 

spans the period from 1931 through 2011. Water levels from 

individual observation wells span shorter periods, were 

obtained over a wide range of hydroclimatic conditions, and 

may not represent the long-term average. 

Water levels from observation wells and surface elevations of water bodies were used to develop a water-table map 

(fig. 7), which was manually interpolated and contoured. The 

water-table map indicates that regional groundwater flow is 

towards the Sudbury River. Groundwater discharges to the 

eastern boundary of Lake Cochituate from both the deep and 

shallow aquifer, but along the western lake boundary, groundwater flow in the shallow aquifer is towards the lake while 

groundwater flow in the deeper aquifer is away from the lake 

towards the Sudbury River. To the west of Lake Cochituate, 

water levels in shallow wells have been measured as much 

as 20 ft higher than water levels in deep wells. The lowpermeability lacustrine deposits in this area may cause locally 

perched water table conditions, hydraulically separating the 

upper and lower parts of the aquifer (fig. 6A). Groundwater 

levels from vertically paired wells in other parts of the study 

area (not shown in fig. 7) indicate little difference between 

shallow and deep parts of the aquifer, indicating that the 

aquifer is unconfined, and that deep and shallow levels of the 

aquifer are probably hydraulically well connected. 

Table 1. Discharge measurements at the outflows from Pod 

Meadow Pond and Dudley Pond, east central Massachusetts.

Date

Discharge 

(cubic feet per second)

Pod Meadow Pond Dudley Pond

03-23-2011 1.18 1.71

05-09-2011 0.79 1.21

06-14-2011 0.65 0.52

07-12-2011 0.43 0.08

08-01-2011 0.54 0.00

09-12-2011 0.74 0.61

10-07-2011 0.45 0.53

11-10-2011 0.141 0.912

11-29-2011 0.121 0.79

12-14-2011 0.58 1.03

01-09-2012 0.86 1.30

02-08-2012 0.68 0.92

03-08-2012 0.75 0.75

Mean 0.70 0.79

Standard deviation 0.21 0.52

1

Discharge affected by beaver dam that was impounding water; values 

excluded from mean and standard deviation.

2

Measurement made on 11-08-2011.
14 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 8. Pod Meadow Pond showing open water along the southern shore in February 2011, Framingham, 

Massachusetts.

North

East

Unfrozen water in shadow of 

north-facing slope where in 

winter water would be 

expected to freeze first and 

stay frozen longest.

Ice

Pod Meadow Pond

Pod Meadow 

Pond

Dudley 

Pond

Bedrock island

Direction photograph

was taken from,

February 11, 2011

Lake 

Cochituate

Area of 

open water
Water Resources 15

Table 2. Borehole and groundwater observation wells used in the study, east central Massachusetts.—Continued

[Well locations shown in figures 7 and 15; Elevation in North American Vertical Datum 1988; SS, used to calibrate steady-state model; TR, used to calibrate transient model; LI, used to determine lithology; WL, water level(s) used to develop model; WS, water-supply well; --, no data available. Data sources are: 1, SEA 

Consultants, Inc. (1992); 2, SEA Consultants, Inc., (2008); 3, Gay (1985); 4, Gay (1981); 5, Balsam (1987, volume I); 6, Balsam (1987, volume II); 7, Balsam 

(1992, volume I); 8, Balsam (1992, volume II); 9, Sovereign Consulting, Inc. (2009); 11, U.S. Geological Survey files in Northborough, Mass. (accessed 2011); 

13, National Water Information System (NWIS) online at http://waterdata.usgs.gov/nwis; 15, Wayland Wellhead Protection Committee and Bruce W. Young 

(2011)]

Well 

identifier

Data 

source Longitude Latitude

Land 

surface 

elevation

(feet)

Hole 

depth

(feet)

Bedrock 

elevation

(feet)

Screen elevation 

(feet)

Screen 

model 

layer

Water-level observations 

(feet)

Top Bottom Count Average Application

1-90 1 -71.3858 42.3301 136.0 63.0 -- 83.0 78.0 4 1 125.4 SS, LI

2-90 1 -71.3864 42.3305 135.5 79.0 56.5 82.5 77.5 4 1 125.5 SS, LI

3-90 1 -71.3872 42.3306 145.1 74.0 -- 94.1 89.1 4 4 127.3 SS, LI

5-90 1 -71.3878 42.3306 131.9 76.0 55.9 78.9 73.9 4 4 127.7 SS, LI

7-90 1 -71.3876 42.3324 145.6 48.0 97.6 102.7 97.7 4 4 120.0 SS, TR, LI

8-90 1 -71.3858 42.3293 172.3 56.0 116.3 88.3 83.3 4 4 127.2 SS, LI

F1W-41 3 -71.3864 42.3297 135.0 61.1 -- 90.0 72.0 4 1 128.9 SS, MD

F1W-42 3 -71.3875 42.3306 142.0 62.3 -- 76.6 61.6 4 1 128.2 SS, MD

F1W-43 3 -71.3886 42.3303 135.0 60.2 76.0 89.3 71.3 4 1 130.9 SS, MD

F1W-60 4 -71.3867 42.3144 138.0 80.1 57.8 62.8 57.8 4 1 137.3 SS, LI

F1W-64 11 -71.3961 42.3342 125.0 23.8 101.2 125.0 120.0 1 1 122.1 SS, LI

F1W-74 5 -71.3789 42.3144 140.0 71.0 69.0 95.0 90.0 4 1 138.2 SS, LI

F1W-84 3 -71.3856 42.3256 139.5 18.0 -- 127.5 124.5 2 16 137.5 SS, TR

F1W-85 3 -71.3847 42.3183 173.0 37.0 -- 149.0 146.0 2 19 149.9 WL, LI

F1W-87 3 -71.3858 42.3272 193.4 65.0 -- 134.4 131.4 2 49 149.6 WL

F1W-88 3 -71.3858 42.3272 193.3 163.0 30.3 104.8 101.8 4 11 127.6 SS, LI

F1W-89 3 -71.3858 42.3272 193.3 163.0 30.3 55.6 52.6 4 8 127.6 SS, TR, LI

F1W-90 3 -71.3856 42.3256 142.7 101.0 45.7 71.7 68.7 4 8 133.1 SS, LI

F1W-91 3 -71.3856 42.3256 142.8 101.0 45.8 47.8 44.8 4 8 134.0 SS, TR, LI

F1W-92 13 -71.3889 42.3142 191.2 203.0 -11.8 120.0 117.0 2 9 127.2 SS, TR, LI

F1W-93 13 -71.3833 42.3289 191.2 203.0 -- 73.3 70.3 4 10 127.2 SS, TR

F1W-94 4 -71.3833 42.3289 191.2 203.0 -11.8 20.5 17.5 4 9 127.1 SS, LI

HH-1 15 -- -- 125.9 42.0 -- -- -- 4 -- -- WS

HH-2 15 -- -- 125.7 47.0 -- -- -- 4 -- -- WS

MV-1 15 -- -- 124.5 80.0 -- -- -- 4 -- -- WS

MW-1 6 -71.3907 42.3305 176.5 60.5 -- 128.5 118.5 4 12 127.1 SS, TR, LI

MW-10 6 -71.3916 42.3320 162.4 49.0 -- 124.4 114.4 2 11 121.1 SS, LI

MW-11 6 -71.3916 42.3320 162.9 81.5 81.4 89.9 84.9 4 12 124.6 SS, LI

MW-12 6 -71.3901 42.3334 123.1 16.5 -- 119.1 109.1 1-2 9 118.3 SS, LI

MW-13 6 -71.3898 42.3310 171.6 54.0 -- 129.6 119.6 3-4 7 127.6 SS, TR, LI

MW-14 6 -71.3898 42.3309 171.4 60.0 -- 121.4 111.4 4 4 127.8 SS, LI

MW-14R 6 -71.3898 42.3309 171.4 60.0 -- 121.4 111.4 4 2 126.9 SS, LI

MW-15 8 -71.3893 42.3307 143.0 20.0 -- 133.0 123.0 2 2 127.0 SS, LI

MW-15D 8 -71.3892 42.3307 143.0 87.0 -- 75.0 65.0 4 2 127.3 SS, LI

MW-16 8 -71.3887 42.3312 142.8 26.5 -- 137.8 117.8 2 2 130.2 SS, LI

MW-2 6 -71.3918 42.3312 172.6 55.5 -- 127.6 117.6 4 12 126.5 SS, LI

MW-2B 6 -71.3918 42.3312 172.7 55.0 -- 129.7 119.7 3-4 11 126.9 SS, LI

MW-2D 6 -71.3917 42.3312 172.4 86.5 85.9 110.4 100.4 4 6 126.1 SS, LI

MW-3 6 -71.3920 42.3314 169.3 55.5 -- 125.3 115.3 4 12 128.2 SS, TR, LI

MW-4 6 -71.3915 42.3315 167.0 45.5 -- 133.0 123.0 3-4 12 124.4 SS, LI
16 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Table 2. Borehole and groundwater observation wells used in the study, east central Massachusetts.—Continued

[Well locations shown in figures 7 and 15; Elevation in North American Vertical Datum 1988; SS, used to calibrate steady-state model; TR, used to calibrate transient model; LI, used to determine lithology; WL, water level(s) used to develop model; WS, water-supply well; --, no data available. Data sources are: 1, SEA 

Consultants, Inc. (1992); 2, SEA Consultants, Inc., (2008); 3, Gay (1985); 4, Gay (1981); 5, Balsam (1987, volume I); 6, Balsam (1987, volume II); 7, Balsam 

(1992, volume I); 8, Balsam (1992, volume II); 9, Sovereign Consulting, Inc. (2009); 11, U.S. Geological Survey files in Northborough, Mass. (accessed 2011); 

13, National Water Information System (NWIS) online at http://waterdata.usgs.gov/nwis; 15, Wayland Wellhead Protection Committee and Bruce W. Young 

(2011)]

Well 

identifier

Data 

source Longitude Latitude

Land 

surface 

elevation

(feet)

Hole 

depth

(feet)

Bedrock 

elevation

(feet)

Screen elevation 

(feet)

Screen 

model 

layer

Water-level observations 

(feet)

Top Bottom Count Average Application

MW-5 6 -71.3913 42.3315 169.4 55.5 -- 125.4 115.4 4 12 125.7 SS, LI

MW-5D 8 -71.3912 42.3314 169.4 89.0 82.9 101.4 91.4 4 2 126.2 SS, LI

MW-6 6 -71.3917 42.3299 176.3 55.5 -- 132.3 122.3 4 2 131.9 SS, LI

MW-7 6 -71.3901 42.3317 173.1 60.5 -- 124.1 114.1 3-4 11 131.0 SS, LI

MW-8 6 -71.3919 42.3326 129.7 19.0 -- 121.7 111.7 1-2 12 120.6 SS, TR, LI

MW-8D 6 -71.3919 42.3326 129.6 56.0 73.6 84.6 74.6 4 6 121.7 SS, TR, LI

MW-9 6 -71.3910 42.3330 125.3 18.5 -- 121.3 111.3 1-2 12 117.9 SS, LI

MW-9D 6 -71.3910 42.3330 125.5 62.5 65.0 77.5 67.5 4 4 119.8 SS, LI

SEA-10 2 -71.3861 42.3282 180.6 72.0 108.6 114.6 109.6 4 3 127.9 SS, TR, LI

SEA-11 2 -71.3892 42.3279 174.2 77.0 -- 102.2 97.2 4 3 129.6 SS, TR, LI

SEA-12 2 -71.3893 42.3316 161.2 85.0 -- 81.2 76.2 4 3 121.8 SS, LI

SEA-13 2 -71.3894 42.3311 155.6 95.5 60.1 93.6 88.6 4 3 127.8 SS, LI

SEA-14 2 -71.3880 42.3314 144.8 39.5 105.3 110.3 105.3 4 3 121.9 SS, LI

SEA-15 2 -71.3888 42.3338 122.0 98.0 -- 58.0 53.0 4 3 117.6 SS, TR, LI

SEA-16 2 -71.3846 42.3265 186.2 123.5 -- 94.2 89.2 4 3 131.7 SS, LI

SEA-17 2 -71.3884 42.3305 139.0 60.0 -- 84.0 79.0 4 3 127.8 SS, LI

SEA-18 2 -71.3869 42.3291 149.0 59.0 90.0 95.0 90.0 4 3 127.7 SS, LI

SEA-2 2 -71.3859 42.3307 128.7 32.5 -- 101.2 96.2 4 3 125.6 SS, LI

SEA-3 2 -71.3844 42.3313 124.0 77.0 -- 65.0 60.0 4 3 122.4 SS, TR, LI

SEA-4 2 -71.3870 42.3300 143.3 80.0 -- 68.3 63.3 4 3 127.6 SS, TR, LI

SEA-7 2 -71.3877 42.3308 134.2 74.0 -- 69.2 64.2 4 3 127.9 SS, LI

SEA-8 2 -71.3863 42.3305 138.3 71.0 67.3 72.3 67.3 4 3 127.0 SS, LI

SEA-9 2 -71.3877 42.3301 139.0 112.0 -- 39.0 34.0 4 3 127.8 SS, LI

TW-1 2 -71.3870 42.3299 143.3 78.0 -- 84.3 69.3 4 -- -- WS, LI

TW-2 2 -71.3884 42.3305 139.0 61.0 -- 89.0 79.0 4 -- -- WS, LI

TW-3 2 -71.3872 42.3306 134.2 73.0 -- 79.2 64.2 4 -- -- WS, LI

TW-4 2 -71.3863 42.3305 138.3 67.2 -- 86.1 71.1 4 -- -- WS, LI

WKW-117 4 -71.3730 42.3181 138.4 32.0 -- 136.9 136.7 1 14 137.6 SS, LI

WKW-118 3 -71.3700 42.3225 181.0 38.0 -- 147.5 144.5 3 23 153.6 WL, LI

WKW-119 3 -71.3822 42.3286 175.1 65.0 -- 117.5 114.5 2 16 129.7 SS, TR, LI

WKW-120 3 -71.3714 42.3200 175.7 49.0 -- 158.7 155.7 -- 16 165.5 WL, LI

WKW-123 3 -71.3786 42.3281 150.2 100.0 -- 53.2 50.2 4 12 138.9 SS, TR, LI

WKW-124 3 -71.3786 42.3281 150.2 100.0 -- 128.9 125.9 2 12 144.2 TR, LI

WKW-2 13 -71.3681 42.3144 153.8 37.5 -- 122.8 120.8 -- -- -- LI

WKW-27 4 -71.3736 42.3194 155.0 76.5 -- 88.6 78.6 -- 1 149.0 WL, LI

WKW-30 13 -71.3883 42.3389 120.0 50.0 -- 80.0 75.0 4 1 118.8 SS, LI

WKW-52 13 -71.3894 42.3425 118.0 61.5 -- 65.0 60.0 4 1 117.0 SS, LI

WKW-53 11 -71.3819 42.3319 135.0 91.0 44.0 54.0 49.0 2 1 127.0 SS, LI

WKW-54 4 -71.3831 42.3314 130.0 82.0 48.0 58.0 53.0 4 1 128.0 SS, LI
Water Resources 17

Subsets of the groundwater-level data were used in model 

calibration. Groundwater level observations from 65 wells 

were used in calibrating the steady-state groundwater model. 

Twelve of these 65 wells have a single measurement, and 

the remaining have a median of 3 measurements spanning a 

median period of 30 days. Wells with more than 3 observations 

had a median difference between the minimum and maximum 

observations of 2.2 ft. Groundwater level observations from 

13 wells were used in calibrating the transient groundwater 

model. A set of groundwater level observations collected 

during a 2006 aquifer test (SEA Consultants, Inc., 2008) also 

was used for calibrating both the steady-state and transient 

models. Many observation wells were excluded from the 

model calibration because of uncertainty about well locations, 

land-surface elevations, well depths, or screened intervals. 

Water Use

A permit from the Massachusetts Department of 

Environmental Protection (MassDEP) is required for all 

public water-supply withdrawals and for large nonpotablewater withdrawals (industrial and agricultural, for example). 

The Town of Wayland has the only permitted groundwater 

withdrawals in the study area and pumps water from three 

production wells north of the Birch Road well site near the 

Sudbury River (MassDEP numbers 3315000–03G; –04G; 

–05G, labeled HH-1, HH-2, and MV-1, respectively, in 

fig. 7). These three wells have an average annual combined 

withdrawal of 1.40 ft3/s, equal to 0.90 million gallons per day 

(Mgal/d). Monthly average withdrawals peak in July with 

seasonal high withdrawals from May through September that 

are 9 to 34 percent greater than the annual average (table 3) 

based on records from 1996–2000 and 2002–2006.

From 1939 to 1979, the Town of Framingham operated wells at the Birch Road site for municipal water supply 

(fig. 7). The production wells consisted of a cluster of three 

wells about 1,400 ft north of Lake Cochituate and about 600 ft 

southeast of the Sudbury River. The wells had a combined 

pumping capacity of 5.3 ft3/s (3.5 Mgal/d) but were likely 

pumped at a rate closer to 4.9 ft3/s (SEA Consultants, Inc., 

1992). The wells were operated intermittently during the midto late 1970s because of high iron and manganese concentrations, and their use was eventually discontinued in 1979. The 

town currently obtains all of its water from MWRA from 

sources outside of the study area. Four large-diameter wells 

(Birch Road wells, fig. 7) that were installed for aquifer tests 

in 2005 may become the supply wells for future withdrawals 

and are treated as such in this study.

Lake Cochituate was the first large drinking-water 

supply for the City of Boston; the system operated from 1848 

through 1951, when it was finally abandoned because of 

declining quality and the availability of other water sources 

(http://www.mwra.state.ma.us/04water/html/hist2.htm). Use 

of Lake Cochituate had been in decline since the early 1900s 

as new water-supply sources came on line. In 1947, the lake 

was transferred to the State, and it is now managed by the 

Massachusetts Department of Conservation and Recreation 

(MassDCR), 2006). State and local municipalities maintain 

parks and boat access at several points along the lake for 

recreational purposes.

Wastewater discharge in the study area varies by town. 

The towns of Framingham and Natick export wastewater to 

the MRWA regional wastewater system, which discharges 

outside of the study area. The towns of Wayland and Sudbury 

are on private septic systems within the study area and return 

wastewater locally to the groundwater system. 

In the basin upstream of the study area are numerous groundwater and surface-water withdrawals that affect 

streamflow. Zarriello and others (2010) reported 6 production 

wells upstream of the study area with total average annual 

withdrawals of about 4.5 ft3/s from 1993 to 2003 in the Town 

of Natick. In the same publication, 24 production wells and 

7 surface-water withdrawals were reported to have operated 

from 1993 to 2003 with a combined annual average withdrawal of about 5.3 ft3

/s. The effects of these withdrawals are 

included implicitly in the groundwater model through specified stream inflows and are particularly evident during periods 

of low flow.

The operation of surface-water reservoirs upstream of 

the study area also affects streamflow and can further decrease 

low-flows during summer months. Three reservoirs in the 

Sudbury Reservoir system are the largest of these and are 

actively managed by the MWRA as an emergency supply. 

Three additional former supply reservoirs farther upstream 

in the basin are managed by the MassDCR. The operation of 

these reservoirs, particularly during periods of low flow, can 

cause large percentage changes in streamflow at the Sudbury 

River at Saxonville (01098530).

Table 3. Average groundwater withdrawal rates from 

Wayland production wells, 1996–2000 and 2002–2006, Wayland, 

Massachusetts. 

Month

Pumping rate 

(cubic feet per second)

HH-1 HH-2 MW-1 Total

Jan. 0.30 0.80 0.05 1.15

Feb. 0.46 0.60 0.09 1.15

Mar. 0.32 0.66 0.12 1.10

Apr. 0.41 0.71 0.16 1.29

May 0.56 0.78 0.18 1.52

June 0.67 0.91 0.25 1.83

July 0.72 0.87 0.29 1.87

Aug. 0.58 0.87 0.24 1.69

Sep. 0.52 0.96 0.19 1.66

Oct. 0.36 0.72 0.10 1.19

Nov. 0.40 0.55 0.16 1.11

Dec. 0.36 0.70 0.13 1.19

Average 0.47 0.76 0.16 1.40
18 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Groundwater Flow Model

A numerical groundwater flow model of the aquifer based 

on MODFLOW-NWT (Niswonger and others, 2011) was 

developed to represent hydrogeologic conditions and simulate 

groundwater flow and interaction with surface waters in the 

study area. MODFLOW-NWT is a finite-difference groundwater modeling software package that is the best available 

tool for simulating aquifers subject to model cells drying. The 

groundwater model was developed from our current understanding of the hydrogeology of the study area as part of a 

multiphase study to address questions about the effects, and 

the potential to mitigate the effects, of pumping on surfacewater resources, and to identify further needs for data collection and model refinements. 

Model Design

The numerical groundwater model was designed to 

simulate the complex aquifer system and the interaction of the 

aquifer with lakes and streams. The model was initially constructed by using MODFLOW-2000 software (Harbaugh and 

others, 2000) and was later converted to MODFLOW-NWT to 

prevent numerical instabilities caused by drying and rewetting 

of model cells. Cell drying is a common problem with models 

of shallow unconfined glacial-sediment aquifers (DeSimone 

and others, 2002; Masterson and others, 2009). Numerical 

instabilities can prevent the model from reaching a solution 

with acceptably small numerical errors. MODFLOW-NWT 

uses a numerical formulation that, unlike MODFLOW-2000, 

keeps aquifer cells active after they are dewatered. In this 

study, we refer to dewatered cells as “dry cells” with the 

understanding that they are still active and, in transient simulations, can become wet again.

The model was initially calibrated to long-term steadystate conditions that represent average annual conditions in 

the study area. After an acceptable steady-state calibration, 

the model was modified to simulate transient, or unsteady, 

conditions. Under steady-state simulations, stresses such as 

recharge from precipitation, groundwater withdrawals, and 

flow to or from water bodies remained constant, but under 

transient simulations, these stresses could change in time. The 

transient model was constructed to represent average monthly 

conditions, as represented by the observed groundwater levels 

spanning the period 1931 through 2011. 

Discretization

The total model area (fig. 9) of 9.0 mi2 was spatially discretized by a uniform grid of square cells, 50 ft on each side, 

spanning 360 rows (north-south direction) and 280 columns 

(east-west direction). This grid size was chosen so that water 

levels affected by the pumping and the steep bedrock topography near pumping wells could be accurately simulated. The 

model boundary positions were chosen to correspond with 

natural no-flow boundaries (no appreciable lateral flows into 

or out of the model area) or to be sufficiently far from the 

main areas of interest so that the boundary would have little 

effect on simulation results. 

Vertically, the model was discretized into five layers of 

variable thickness that represent different aquifer sediments 

as shown in figure 10. Layer 1 represents the top 10 ft of 

sediment below the land surface, except in a few areas 

noted below. Layers 2 through 4 have varying thicknesses 

representing aquifer sediments that collectively are as much 

as 261 ft deep. Layer 2 represents permeable sand and gravel 

deposits near the surface. Layer 3 represents fine lacustrine 

deposits separating the coarser, more permeable sediments in 

layers 2 and 4. Layer 4 represents permeable sand and gravel 

deposits deeper in the aquifer, including those below the fine 

deposits found at depth in the Lake Cochituate area. In areas 

where the hydrogeologic units that layers 2–4 represent are 

absent, a 0.5 ft thickness was assigned to these layers to meet 

model input requirements. At the margins of the active model 

area, the total thickness of layers 2–4 is as little as 1.5 ft. Layer 

5 (the bottom layer) represents the top 80 ft of bedrock and the 

mostly thin layer of till over the bedrock. The upper bedrock 

was included in the model because it is generally more 

fractured than the deeper bedrock and is considered a zone of 

active groundwater flow. 

Layers 2 and 4 define the principal aquifer in the study 

area and in some areas are hydraulically separated by the fine 

deposits represented in layer 3. The fine-grained component of 

layer 3 is most extensive in the southern part of the model area 

under and around Lake Cochituate, pinches out to the north of 

Lake Cochituate, and does not extend to the Birch Road wells. 

Layer 3 also represents extensive fine deposits in the northern 

part of the model area. 

Elevations of the land surface at the top of layer 1 were 

determined for each cell by interpolation of point elevation 

data derived from 1:5,000 orthophotos (MassGIS, 2003). 

Layer 1 extends to a depth of 10 ft, except in areas where bedrock is close to the surface or where Lake Cochituate is present. Layer 1 thins to 0.5 ft where bedrock is near the surface, 

mostly along the eastern edge of the active model area. Lake 

Cochituate, which is up to 65 ft deep, is represented in layer 1 

with the thickness determined from bathymetric data.

Layer 2 ranges from 0.5 to 125 ft thick and is thinnest 

where low-permeability sediments are at the surface in the 

northern and middle parts of the active model area, at the 

bedrock island, and in eastern and southwestern boundary 

cells where bedrock is near the land surface. Layer 3 ranges 

from 0.5 to 104 ft thick, with the thickest parts in the southern 

model area. In the northern model area, layer 3 is as much 

as 86 ft thick. Layer 4 ranges from 0.5 to 208 ft thick and 

is thickest in the bedrock valley in the southern part of the 

model. Where layer 3 is 0.5 ft thick, indicating an absence 

of fines, layers 2 and 4 were assigned a thickness equal to 

one-half of the remainder of the thickness of the glacial stratified deposits.
Groundwater Flow Model 19

Figure 9. Active model area and boundary conditions, east central Massachusetts. Sudbury River 

Sudbury River 

Sudbury River 

Pod Meadow 

 Pond

Dudley Pond

Lake Cochituate

Bedrock island

Birch Road wells Birch Road wells

MV-1

HH-2 HH-1

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET

0 250 500 METERS

No-flow

boundary

EXPLANATION

Inactive model area

Zone of high hydraulic

 conductivity representing 

 water

Simulated stream and 

 model segment number

USGS streamgage and 

 number

01098500

USGS partial-record site

Wayland production well, 

 and identifier 

Specified stream-inflow

 location and value, in 

 cubic feet per second

Birch Road wells

1

2

1.4

0.2

107.5

11.6 1.4

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

21

20

19

1

MV-1

Bedrock island

Grid cells 50 feet square

Pod Meadow Pond

Pod Meadow Pond 

outflow

Simulated

Birch Road wells

Simulated Pod Meadow Pond

Pod Meadow Pond outflow

Pod Meadow Pond 

outflow

01098500

01098530

Dudley Pond outflow

Oxbow

Pod Meadow

406080

100120140160180200220240260

280

40

60

80

100

120

140

160

180

200

220

240

260

280

300

320

340

360

Row 20

Column 20

Column numbers

Model grid: 360 rows and 280 columns Row numbers

Row 250, E-E’ (figure 11) Column 125, A’’-A’’’ (figure 11)

E E'

A''

A'''
20 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 10. Variability of hydraulic conductivity values of groundwater-model layers 1–4 used in the 

calibrated model, east central Massachusetts. Zones defined based on sediment morphosequences (figs. 5, 

6), borehole data, and geospatial interpolation.

70

0.2

50,000

10

Simulated

streams 

0.01

60

10

140

20

70

10

0.01

EXPLANATION

0.01

Hydraulic conductivity values, in feet per day

Layer 5 (not shown), constant value of 0.2 feet per day

0.2

10

20

60

70

140

50,000

Inactive

model area

Inactive

model area

Inactive

model area

Inactive

model area

C. Layer 3 D. Layer 4

A. Layer 1 B. Layer 2

N N

N N
Groundwater Flow Model 21

Time is discretized in transient simulations to represent 

month-to-month variability of the annual hydrologic cycle. 

Transient simulations of the annual cycle were repeated 

for 5 years to remove the effects of initial conditions and 

establish a dynamic equilibrium, in which simulated monthly 

water levels and fluxes are very nearly the same from year 

to year. Sixty monthly stress periods, representing 5 years of 

average monthly conditions, adequately established dynamic 

equilibrium conditions, as indicated by total storage changes 

of less than 1 percent in the fifth year. Simulation results 

were reported for the fifth year (months 49–60); results from 

the first 4 years (months 1–48) were used only to establish 

a dynamic equilibrium. The one exception to monthly stress 

periods is a simulation of a 20-day aquifer test, which was 

used for the transient model calibration described later in 

the report.

Model Boundaries

The active model area is surrounded by a no-flow boundary (fig. 9) that generally coincides with surface-water divides 

and areas of thin sediment. Bedrock is not near the surface 

at the southwestern and northeastern boundaries, where the 

Sudbury River enters and exits the model, or at the southern 

boundary where Lake Cochituate is located. Groundwater 

flows were assumed to be negligible through sediments underlying the Sudbury River at the model boundaries. At both the 

southwestern and northeastern boundaries, where the Sudbury 

River respectively enters and exits the model area, sediments 

are relatively thin and have relatively low permeability (fig. 4). 

At the northeastern boundary, surface water gradients are low 

over a wide wetland area, suggesting that groundwater gradients, and hence flow rates, are also low across the boundary. 

Further, these areas are far enough from the area of interest 

that the likely limited groundwater flows through them would 

have minimal effect on the model simulations.

Inflows to Lake Cochituate were assumed to be fully 

represented by the assigned stream inflow at the southern end 

of the lake. Lateral groundwater flow in the sediments beneath 

Lake Cochituate across the southern model boundary is not 

likely to be significant because north-south surface-water gradients are low in this area. However, there were few groundwater level observations available to make a full assessment 

of groundwater flow under the lake. Groundwater flow across 

the northern boundary is likely small as indicated by the low 

groundwater table gradient (fig. 7). Further, the potential for 

this boundary to affect simulation results near the Birch Road 

wells is minimal because of the large distance between the 

boundary and the wells. On the basis of observed groundwater 

levels (fig. 7), the Sudbury River acts as a hydraulic boundary having no groundwater flow beneath it and hydraulically 

isolating the pumping wells from areas to the west and north 

of the river. The bottom of layer 5, which is 80 ft beneath the 

top of bedrock and thin till, is considered an impermeable noflow boundary. 

The only exceptions to no-flow boundaries around the 

perimeter of the active model area are inflows and outflows 

to streams represented by the streamflow-routing package 

(SFR2; Niswonger and Prudic, 2005). Flows were specified 

at the first upstream cell representing the Sudbury River, 

Dudley Pond (to account for the small drainage area to the 

pond that is outside the active model area), and a small 

tributary to Cochituate Brook that has a small drainage area 

outside the model boundary. Recharge from precipitation was 

applied uniformly over the model surface through the recharge 

package (RCH) in MODFLOW-NWT to the highest active 

model cells. Outflow from the model includes discharge of 

the Sudbury River at the northeastern boundary, groundwater 

withdrawals by the Wayland production wells in the northern 

part of the model, and for some simulations, hypothetical 

groundwater withdrawals from the Birch Road wells. Outflow 

by evapotranspiration was not explicitly simulated but was 

incorporated into the net recharge values used in the model. 

Hydraulic Conductivity

The rate of groundwater flow is determined by hydraulic 

conductivity and water-level gradients. Hydraulic conductivity of aquifer material was assigned on the basis of sediment 

lithology (table 4; figs. 10, 11). Values were initially adopted 

from previous studies (DeSimone and others, 2002; Masterson 

and others, 2009) and then adjusted during model calibration. 

Each model layer can contain a variety of sediment types in 

this complex glacial-sediment aquifer, so hydraulic conductivity varied somewhat but was generally uniform within a layer 

(fig. 10). The ratio of horizontal to vertical hydraulic conductivity (anisotropy) varied by layer from 5 to 20 (table 4). Vertical anisotropy values were initially assigned based on previous 

studies (DeSimone and others, 2002; Masterson and others, 

2009) and then were adjusted during model calibration.

Sand and gravel are the most permeable sediments in the 

aquifer and were assigned the highest hydraulic conductivity 

values. Sediments are slightly coarser in layer 4 than in layer 

2 and were assigned a hydraulic conductivity value higher in 

layer 4 than layer 2. Fine sediment deposits represented by 

layer 3 contain a wide mix of sediment types but commonly 

consist of fine sand with silt; hence, these deposits do not 

represent an aquitard (impermeable boundary to groundwater 

flow) but rather a more restrictive layer to groundwater flow 

than the coarse-grain deposits. Bedrock and till in layer 5 were 

considered slightly permeable for groundwater flow and were 

assigned a hydraulic conductivity equal to that of fine deposits. The gyttja deposits are generally impermeable and were 

assigned a low hydraulic conductivity. Gyttja deposits were 

represented in the model by a low-conductivity layer below 

Dudley Pond on the basis of the coring study (IEP, Inc., 1983). 

However, a low-conductivity gyttja layer was not modeled 

below Lake Cochituate; although gyttja deposits may be present in some areas, no data were available to guide the areal 

distribution of the layer.
22 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Model cells representing Dudley Pond, Pod Meadow 

Pond, and Lake Cochituate (layer 1) were assigned a high 

hydraulic conductivity value of 50,000 ft/d. The high 

conductivity allows groundwater to move with little resistance 

between cells, causing groundwater levels to be nearly 

uniform within the cells, mimicking an open water body. 

This method of simulating surface water bodies has been 

analyzed in other studies (Hunt and others, 2003) and used in 

groundwater models of other northeastern glacial sediment 

aquifers (Masterson and others, 2009; DeSimone, 2004). In 

Dudley Pond, where gyttja deposits are well documented (IEP, 

Inc., 1983), model cells beside and beneath the pond in layers 

1 and 2 were assigned a low hydraulic conductivity value 

of 0.01 ft/d to maintain a constant water level that matched 

observed conditions.

Storage Coefficients

The storage coefficients, specific storage (Ss) and specific yield (Sy), define the volume of water going in or out 

of aquifer storage as groundwater levels rise and fall under 

transient flow conditions. When a model cell is saturated, Ss in 

units of inverse feet (ft-1) is used to define the volume of water 

released per unit volume of aquifer per unit of groundwater 

level change. When a model cell is partially saturated, Sy (percent) is used and is defined as the volume of water that drains 

by gravity per unit volume of aquifer. MODFLOW-NWT 

automatically determines which storage coefficient to use from 

the simulated water level and top and bottom elevations of 

each cell. Uniform values were used to define Ss (0.00001 ft-1) 

and Sy (15 percent) over all model layers (table 4). The values 

were initially assigned based on aquifer test results (SEA, 

1992) and later modified during model calibration. 

Streams

Streams are represented in layer 1 with the StreamflowRouting Package (SFR2) by Niswonger and Prudic (2005). 

SFR2 simulates flow between streams and the aquifer as 

a head-dependent flux and tracks streamflow as it moves 

downstream. Modeled streamflows represent baseflow, the 

groundwater component of streamflow, and do not include 

surface runoff or interflow components of streamflow. Simulated streamflow in the model is the combination of specified upstream inflows and groundwater gains or losses to the 

stream computed by the model. Streams lose water to the 

aquifer when groundwater levels are below the stream level 

and gain water when groundwater levels are above the stream 

level. A stream reach can become dry if stream infiltration into 

the aquifer exceeds inflow from upstream reaches.

Streams in the model area were defined by 21 segments 

comprising 1,219 cells in layer 1 that represent the Sudbury 

River, Cochituate Brook, and minor tributaries, as well as 

outflows from Dudley and Pod Meadow Ponds and Lake 

Table 4. Hydraulic conductivity of aquifer material and storage coefficients by layer in the calibrated model, east 

central Massachusetts.

[--, no value assigned; bold values indicate the dominant type in the layer]

Layer 1 Layer 2 Layer 3 Layer 4 Layer 5

Horizontal hydraulic conductivity 

(feet per day)

Sand and gravel 70 60 70 140 --

Fines (fine sand, silt, and clay) 10 10 10 20 --

Till at surface 0.2 -- -- -- --

Bedrock at surface 10 10 10 20 --

Bedrock and till at depth -- -- -- -- 0.2

Pond muck (gyttja) 0.01 0.01 -- -- --

Lake and pond cells 50,000 -- -- -- --

Ratio of horizontal to vertical hydraulic conductivity

Anisotropy 5 10 20 10 5

Storage coefficients

Specific yield (fractional percent) 0.15 0.15 0.15 0.15 0.15

Specific storage (1/feet) 0.00001 0.00001 0.00001 0.00001 0.00001
Groundwater Flow Model 23 Figure 11. Vertical cross sections of groundwater model layers and hydraulic conductivity values running A, north-south and B, east-west, east central Massachusetts.

0

0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000 11,000 12,000 13,000 14,000 15,000 16,000 17,000 18,000

South North

A. Modeled layers and hydraulic conductivities, column 125

Layer 1

Layer 2

Layer 1

Area of Birch Road wells

Approximate length of cross section A-A’ in Figure 6

Bedrock island

200100

-100

10

Vertical exaggeration X 10

Elevation in NAVD 88

Cochituate Brook Sudbury River

Pod Meadow

Distance along model grid, in feet

Elevation, in feet 

10 10 20

70

70 Fines

20

Fines

Fines Fines

West East

Vertical exaggeration X 10

0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 9,000 10,000 11,000 12,000 13,000 14,000

Layer 2

Layer 1

Layer 4

Layer 5

300200100

-100-200

0

Sudbury River

Distance along model grid, in feet

Elevation, in feet 

Lake Cochituate

Fines

Fines

Sand and gravel

10

10

10 20 70 70 50,000

140

Fines

B. Modeled layers and hydraulic conductivities, row 250

EXPLANATION

Hydraulic conductivity, in feet per day

Inactive model area

Inactive 

Inactive 

Inactive 

Bedrock and till

60

Layer 3

A''' A''

E E'

Values as labeledFinesHigh conductivity zone for lake

60

140

Fines

Sand and gravel

Layer 3

Layer 2

Layer 4

Layer 5

10

0.2

60 140 Sand and gravel

Bedrock and till

Layer 4

Layer 3

Layer 5

0.2
24 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Cochituate (fig. 9). Segments 5, 6, 8, 16, 18, and 21 define the 

Sudbury River, segments 1, 2, and 4 define Lake Cochituate 

and its outlet, segments 9–15 define Pod Meadow Pond and 

its outlet and wetland areas, and segments 19 and 20 define 

Dudley Pond and its outlet (table 5). Stream cell bottom 

elevations were determined from the digital elevation data 

of the land surface minus 1.0 ft to ensure that the channel 

bottoms were below land surface. Stream cell lengths are 

equal to the uniform model grid cell size of 50 ft, streambed 

thickness values were set to 1.0 ft, while stream widths and 

streambed hydraulic conductivity varied by stream (table 5). 

The hydraulic conductivity of the Sudbury River streambed 

downstream of the oxbow was increased by a factor of 10 

relative to hydraulic conductivity upstream of the oxbow 

to simulate the groundwater/surface-water exchange in the 

adjoining wetlands and the diffuse shallow stream network in 

the wetlands. Depths of water in the streams were calculated 

by the SFR2 package by using a uniform Manning’s roughness 

coefficient of 0.03. 

Stream cells were assigned across the high conductivity zones representing Dudley Pond, Pod Meadow Pond, and 

Lake Cochituate (fig. 9). These stream cells were assigned 

streambed elevations equal or close to the average water level 

of 138.5 ft for Lake Cochituate, 125.8 ft for Pod Meadow 

Pond, and 153.0 ft for Dudley Pond. Stream cells traversing 

the lake and ponds drain groundwater when it rises above the 

stream level but provide a source of recharge to the aquifer 

when stream levels exceed adjacent groundwater levels.

Inflows were specified to streams entering the active 

model area on the basis of daily flow records from the 

Saxonville streamgage (01098530) from January 1980 through 

December 2010. For the steady-state model, inflows of 107.5, 

11.6, and 1.4, ft3/s were assigned to the Sudbury River, Lake 

Cochituate, and the tributary to Cochituate Brook, respectively 

(fig. 9). These values total slightly less than the median daily 

flow observed at Saxonville (148 ft3/s) and were proportioned 

by the drainage area to each boundary inflow point. A small 

inflow of 0.2 ft3/s was assigned to Dudley Pond based on the 

small drainage area to the pond outside of the active model 

area. Streamflow exits the model through the Sudbury River. 

In the transient model simulations, stream inflows were 

the same for all scenarios and were assigned to represent 

monthly low-flow conditions. The 75-percent daily flow 

duration (flow value exceeded 75 percent of the time) was 

determined for each month of the year from recorded daily 

flows at the Saxonville streamgage from January 1980 through 

December 2010. The 75-percent daily flow durations by 

month at Saxonville were proportioned by drainage area, 

as previously described, to specify inflows at the model 

boundary. No inflow was assigned to Dudley Pond for 

transient simulations. The 75-percent daily flow duration was 

chosen to represent conservative monthly low-flow conditions 

throughout the year in addition to the normal seasonal 

variations in flows. For example, the Sudbury River at 

Saxonville (01098530) streamgage had a 75-percent daily flow 

duration for September of 18.0 ft3/s, whereas the 90-percent 

daily flow duration for the entire period of record, a commonly 

used low-flow stress indicator threshold, is 44 percent greater 

(26.0 ft3/s). 

Table 5. Stream segments and properties in the groundwater flow model, east central Massachusetts.

Stream 

segment 

number

Description Flows into 

segment

Streambed 

hydraulic 

conductivity 

(feet per day)

Width

(feet)

Steady-state 

flux

1 Lake Cochituate traverse stream 2 20 50 Losing

2–4 Cochituate Brook and tributary 6 20 15 Gaining

5, 6 Sudbury River upstream of oxbow 8 20 45 Gaining

7 Tributary 8 20 5 Gaining

8 Sudbury River at oxbow 16 20 45 Gaining

9 Pod Meadow Pond traverse 10 10 5 Gaining

10–15 Pod Meadow outlet, Pod Meadow swamp, and northern part of oxbow 16 20 5 Gaining

16, 18 Sudbury River downstream of oxbow and wetlands 21 200 45 Gaining

17 Heard Pond traverse and outlet 18 200 45 Gaining

19, 20 Dudley Pond outlet 21 2 5 Gaining

21 Sudbury River and wetlands Exits model 200 45 Gaining
Groundwater Flow Model 25

Recharge

Recharge from precipitation is the primary inflow to 

the aquifer and was simulated by using the MODFLOWNWT RCH package. Recharge in the steady-state model 

was assigned a spatially uniform value of 22 in/yr based on 

reported values for nearby watersheds (Zarriello and others, 2010; DeSimone, 2004; DeSimone and others, 2002). 

Recharge was assigned to all active cells except the high 

conductivity cells representing Lake Cochituate, as these cells 

were already saturated from the surface-water inflow to the 

lake and the stream traversing the lake area. 

For the transient model, recharge rates were varied 

monthly to reflect seasonal precipitation and evapotranspiration variability over the year. Monthly recharge rates, as 

percentages of annual recharge (table 6), were determined 

from soil moisture budgets calculated with WATBUG software (Wilmott, 1977) by using the Thornthwaite evaporation 

method from daily temperature and precipitation obtained 

from the Natick, Mass., weather station (NOAA-COOPID 

195175). The recharge rates are similar to values from groundwater studies in a neighboring basin by DeSimone (2004). 

Pumping

Groundwater withdrawal data that were included in 

all model simulations of the study area were from three 

production wells operated by the Town of Wayland. Longterm average and monthly average pumping rates (table 3) 

used in the steady-state and transient models were determined 

from reported 1996–2000 and 2002–2006 values. Hypothetical 

groundwater withdrawals from the Birch Road wells were 

included in the model for scenarios described later in 

the report.

Model Calibration

Hydraulic parameter values and boundary conditions 

were manually adjusted to improve the match between 

simulated and observed groundwater levels, lake levels, and 

streamflow. Adjustments of hydrologic parameter values and 

boundary conditions were made within a range of values 

determined from previous groundwater models of glacialsediment aquifers in Massachusetts. Adjustment of parameters 

in the NWT solver package also was required to achieve 

model convergence and reasonable run times (table 7). 

Calibration efforts were designed to determine the hydrologic 

controls on groundwater-surface water interaction and to 

determine areas of data need, rather than to obtain an optimum 

parameter fit. 

Steady-State Calibration

Steady-state calibration focused on streambed hydraulic conductivity and vertical and horizontal aquifer hydraulic conductivity to reduce simulation errors in pond levels, 

groundwater levels, and streamflows. Water levels in Dudley 

Pond, Pod Meadow Pond, and Lake Cochituate were used as 

calibration targets for the steady-state model (table 8). Lake 

and pond surface levels were determined from digital landsurface elevation data (MassGIS, 2003) so that they would be 

consistent with the land-surface elevations used in the model. 

These levels were given particular attention during calibration 

of streambed hydraulic conductivity for streams traversing the 

ponds. Measured streamflows coming out of Dudley and Pod 

Meadow Ponds (table 9) were also used in calibration of the 

steady-state model. 

For steady-state groundwater model calibration, the mean 

water-level observations from 65 wells (table 2) were used as 

calibration targets, and model error was calculated from the 

simulated minus the observed mean groundwater levels. The 

variability of groundwater level measurements over different 

seasons and climatic conditions is a source of uncertainty in 

developing a long-term steady state representation of water 

levels. The lack of systematic groundwater observations over 

time is a source of uncertainty that increases with increased 

groundwater level variability, fewer observations, and shorter 

periods of record. 

Model cell drying under steady pumping and climatic 

stresses presented a challenge and in some cases prevented the 

model from reaching a solution (converging). Convergence 

problems were addressed by changes to stress periods, time 

steps, and SFR and NWT input parameters. Under steady-state 

conditions, most cells in layer 1 dry except in low-lying areas; 

progressively fewer cells dry in lower layers 2–4 (fig. 12), 

and no cells dry in layer 5. In the transient model simulations, 

added pumping stresses and decreased recharge cause additional cells to dry, but cells can also rewet as appropriate. 

Table 6. Monthly recharge rates used in the transient groundwater flow model.

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

Inches 2.3 2.4 2.8 2.5 1.6 0.8 0.7 0.3 0.6 2.5 2.1 3.3 22.0

Percent of annual 10.6 11.1 12.9 11.4 7.1 3.8 3.0 1.3 2.8 11.4 9.7 14.9 100
26 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Table 8. Hydrologic characteristics and observed and simulated water levels of Pod Meadow Pond, Dudley Pond, and Lake 

Cochituate, east central Massachusetts.

[Observed and simulated levels in North American Vertical Datum 1988. Pond locations are shown in figure 1]

Pond Drainage area 

(square miles)

Hydraulic conductivity 

of bottom

Water level 

(feet)

Average 

observed

Average 

simulated

Difference, 

simulated-observed

Observed annual 

variation

Pod Meadow Pond 0.23 High 125.8 128.0 2.2 Unknown

Dudley Pond 0.58 Low 153.3 153.6 0.3 Unknown

Lake Cochituate 17.5 Mixed 138.5 139.5 1.0 ±2.1

Table 7. MODFLOW-NWT input file variable values for 

NWT package for transient simulations of current conditions 

(scenario 4).

Parameter Value

HEADTOL 0.01 (feet)

FLUXTOL 1,000 (feet per day)

MAXITEROUT 500

THICKFACT 0.00001

LINMETH 2

IPRNWT 1

OPTIONS Specified

IBOTAV 1

DBDTHETA 0.7

DBDKAPPA 0.00001

DBDGAMMA 0

MOMFACT 0.1

BACKFLAG 0

MAXBACKITER 20

BACKTOL 2

BACKREDUCE 0.6

IACL 2

NORDER 0

LEVEL 3

NORTH 3

IREDSYS 0

RRCTOLS 5.00E-04

IDROPTOL 0

EPSRN 5.00E-04

HCLOSEXMD 1.00E-03

MXITERXMD 150

Streambed vertical hydraulic conductivity influences the 

timing and magnitude of flows between streams and the aquifer. Reducing streambed conductivity decreases flow between 

the stream and the aquifer and increases the differences in 

water level between them. For each stream segment, long-term 

average groundwater levels in nearby observation wells were 

used to calibrate the streambed conductivity. For Pod Meadow 

Pond and Dudley Pond outflows (segments 9 and 19, respectively; fig. 9), the average measured flows (table 1) were also 

used to calibrate streambed and aquifer hydraulic conductivity 

values. Streambed conductivity (table 5) was the only stream 

property adjusted during calibration. 

Aquifer hydraulic conductivity values (table 4) were 

adjusted during the calibration process by changing values 

according to aquifer material and the layer multiplication 

factor. Vertical hydraulic conductivity was set by the ratio 

to horizontal hydraulic conductivity and was fixed at 5.0 for 

layers 1 and 5 and calibrated for layers 2–4 (table 4). Vertical anisotropy was largely determined by calibration to water 

levels in paired wells screened in the upper and lower parts of 

the aquifer.

Pond levels were slightly oversimulated in the 

steady-state model (table 8). When streambed and aquifer 

conductivities were calibrated, there was a tradeoff between 

simulation errors in pond levels, groundwater levels, and 

outflows from Dudley and Pod Meadow Ponds. Minimizing 

differences between simulated and observed groundwater 

levels was given priority over doing the same for pond 

levels and outflows. The oversimulation of pond levels was 

considered acceptable because the differences between the 

simulated and observed levels were well within the observed 

ranges of natural variability.

Simulated streamflow was evaluated for the stream 

segments draining Dudley Pond and Pod Meadow Pond by 

comparing simulated steady-state baseflow in each stream 

to observed average discharge measurements made from 

March 2011 to March 2012 (table 1). Simulated flow is about 

74 percent higher than observed outflow from Pod Meadow 
Groundwater Flow Model 27

Table 9. Simulated flows between Pod Meadow Pond, Dudley Pond, and Lake Cochituate and 

the aquifer, east central Massachusetts.

[Negative values are outflows from pond]

Pond water fluxes 

(cubic feet per second)

Pod Meadow Pond Dudley Pond Lake Cochituate

Steady-state current conditions

From groundwater 1.23 0.02 0.27

To groundwater -0.03 -0.06 -1.95

Net flux 1.20 -0.05 -1.68

Birch Road wells pumping at 4.9 cubic feet per second

From groundwater 0.93 0.01 0.17

To groundwater -0.95 -0.07 -3.43

Net flux -0.02 -0.05 -3.26

Change caused by pumping -1.22 -0.01 -1.58

Pond (1.22 and 0.7 ft3/s, respectively) and 16 percent lower 

than observed outflow from Dudley Pond (0.67 and 0.80 ft3/s, 

respectively). As discussed previously, Pod Meadow Pond 

outflows are thought to be a result of focused groundwater 

discharge and, therefore, were a priority in the model calibration, although these flows were difficult to simulate accurately. 

Oversimulation of Pod Meadow Pond outflows was likely 

caused by poor model representation of local aquifer hydraulic conductivity patterns around the Birch Road wells and 

the pond. Outflows from Dudley Pond are largely affected 

by climatic conditions, as can be seen from the Dudley Pond 

outflow stream going dry in August 2011 (table 1). Furthermore, because simulated outflows from Dudley Pond reflect 

the simulation of baseflow but not runoff and interflow, the 

undersimulation of measured Dudley Pond outflows was not 

considered unreasonable. Overall, the accuracy of simulated 

baseflow was considered acceptable for the purposes of this 

study given the difference between baseflow and total streamflow and the short period of available observations. 

Steady-state simulated groundwater levels are generally 

in agreement with the estimated groundwater table (fig. 14). 

Groundwater level simulations are acceptably accurate with a 

mean model error of -0.45 ft, average absolute error of 2.8 ft, 

and a balanced distribution of oversimulated and undersimulated errors (fig. 13). 

Not all model error is attributable to limitations of 

observed data. In some parts of the model area, steady-state 

simulated and observed groundwater levels did not agree as 

indicated by spatial patterns of simulation error around the 

Birch Road well area (fig. 14). The spatial patterns are likely 

caused by local hydraulic conductivity variations associated 

with the complex-sediment texture variability in this area. A 

previous single layer groundwater model of the Birch Road 

well area (SEA, 2008) included several zones of varying 

hydraulic conductivity near the Birch Road wells to reduce 

simulation errors. Hydraulic conductivity patterns assigned to 

the model in this current study were largely determined from 

sediment cores and glacial morphosequences, which did not 

indicate coherent patterns of local-scale hydraulic conductivity variations in the Birch Road well area. Local variations 

in hydraulic conductivity could be examined in future model 

modifications by using pilot point or regularization techniques 

to automate and optimize parameter value assignment. Efforts 

to remove local error patterns in the Birch Road well area by 

constructing and calibrating a detailed hydraulic conductivity 

field were not undertaken in this stage of the study because 

these errors are not expected to substantially affect answers to 

study questions.

The calibrated steady-state model also did not accurately 

reproduce high groundwater levels in some shallow 

observation wells close to Lake Cochituate that are seen in the 

map of observed water level measurements (fig. 7). Shallow 

wells F1W-85, -87, and WKW-27, -118, -120, -124 had 

high water levels; wells F1W-84 and F1W-88 did not. High 

observed water levels in some of these wells are thought to 

reflect locally perched water tables that could be reproduced 

in the model by lowering hydraulic conductivity values in 

cells in these areas, but doing so would come at the expense 

of increased model error at other locations. The borehole 

sediment logs for these areas did not justify changes to the 

model hydraulic conductivity values, but future refinements 

to spatial distribution of hydraulic conductivity could improve 

the model fit. 
28 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 12. Model cells in layers 1–4 that dry in the calibrated steady-state model, east central Massachusetts.

C. Layer 3 D. Layer 4

A. Layer 1 B. Layer 2

Simulated

streams 

Pod Meadow 

Pond

Dudley Pond

Lake 

Cochituate

Lake 

Cochituate

Lake 

Cochituate

Lake 

Cochituate

Pod Meadow 

Pond

Dudley Pond

Pod Meadow 

Pond

Dudley Pond

Pod Meadow 

Pond

Dudley Pond

EXPLANATION

Active model cells: Dry Wayland production well

Saturated

Inactive model area

Inactive

model

area 

Inactive

model

area 

Inactive

model

area 

Inactive

model

area 

Water bodies in layers 2–4 

shown for reference

N N

N N
Groundwater Flow Model 29

Figure 13. Observed and steady-state 

simulated groundwater levels, east 

central Massachusetts.

110

120

130

140

150

110 120 130 140 150

Simulatedwater level, in feet

Observed water level, in feet

Transient Model Calibration

Transient model calibration used data for seasonal 

groundwater levels and pumping-influenced groundwater 

levels collected during the 2006 aquifer test, but lake and 

pond elevations were not used for the transient model calibration because they are strongly affected by surface runoff and 

evapotranspiration, which were not explicitly simulated. Storage coefficients Ss and Sy were adjusted during calibration but 

were assigned uniform values for all layers (table 4). Uniform 

storage values were assigned because available information 

does not show consistent relations between Sy, which is the 

dominant storage parameter, and sediment texture or hydrogeologic unit.

The transient model was calibrated to observations of 

ambient groundwater levels made at 13 wells (table 2) that 

have a minimum of 6 water-level observations taken over 

a period of at least 150 days. The 13 wells (fig. 15) have a 

median of 10 observations per well and a median period of 

record of 346 days. The Head-Observation Package (HOB) 

by Hill and others (2000) was used to simulate groundwater 

levels at the times the observations were made, and differences 

between simulated and observed water levels define simulation error. Initial hydraulic head conditions for the transient 

simulation were based on steady-state simulated heads. Once 

the transient model was running, simulated hydraulic heads for 

the month of January under current conditions were assigned 

as initial conditions.

The transient model also was calibrated to groundwater 

level data collected during the April 26–May 15, 2006 aquifer 

test of the Birch Road wells (SEA, 2008). Transient simulation of the aquifer test (tn50m8) imposed the varying pumping 

rates of the four wells during the 20-day test. Similar to the 

other transient simulations, stress periods 1–40 in run tn50m8 

were used to establish only monthly dynamic equilibrium conditions; thereafter, a daily stress period was used (41 through 

60) to represent the aquifer test. Model boundary conditions 

(stream inflows, recharge, pumping stresses) were assigned 

to reflect known conditions during the dates of the test. For 

example, average April and May recharge rates were applied 

on the appropriate days. Groundwater drawdowns measured in 

seven wells (fig. 15) at the end of the test were used for model 

calibration (table 10). 

Simulated transient water levels with no pumping at 

the Birch Road wells (scenario 1 in table 13) have a mean 

simulation error of -0.21 ft and an absolute mean simulation 

error of 3.55 ft measured by the difference between simulated 

and observed groundwater levels at 13 observations wells. 
30 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Figure 14. Simulated steady-state groundwater levels and errors, east central Massachusetts.

90

Pod Meadow 

 Pond

Dudley Pond

Heard Pond

Lake Cochituate

Sudbury River 

Bedrock island Bedrock 

island

EXPLANATION

Simulated groundwater elevation,

in feet (NAVD 88)

Inactive 

model area

Wetland

120

Observation well used for 

steady-state model calibration; 

simulated minus observed 

groundwater elevation, in feet 

Wayland production well, 

and identifier

Birch Road wells

< -4

-4 to < -1

-1 to 1

> 1 to 4

> 4

Oxbow

Pod Meadow

Birch Road 

wells

Simulated

streams 

MV-1 HH-2 HH-1

HH-2

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET

0 250 500 METERS

Simulated groundwater levels in NAVD 88

140

130

140

140

150

150

150

120

160

160

160

110

130

130

130

130

120

120

120

Cochituate Brook
Groundwater Flow Model 31

90

Pod Meadow 

 Pond

Dudley Pond

Heard Pond

Lake Cochituate

Sudbury River 

Bedrock 

island

EXPLANATION

Inactive 

model area

Wetland

Model stream cells

Observation well, 

and identifier—

Wayland production 

well, and identifier 

Birch Road wells

Used for transient model 

 calibration

2008 aquifer test (table 10),

 MW-1 also used in transient

 model calibration 

Not used for calibration,

 but described in report

Oxbow

Pod Meadow

Simulated

streams 

MV-1

HH-2 HH-1

MV-1

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET

0 250 500 METERS

WKW-119

WKW-123WKW-124

F1W-84

F1W-91

F1W-91

F1W-85

F1W-85

F1W-87

WKW-27

WKW-120

WKW-2

WKW-118

F1W-89

F1W-93

MW-8

MW-8D

MW-1

SEA-11

SEA-9

SEA-10

SEA-3

SEA-15

7-90

MW-3

MW-13

F1W-92

SEA-11

Cochituate Brook

Sudbury

River

Figure 15. Wells used for transient model calibration of a glacial-sediment aquifer model in east central Massachusetts.
32 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Table 10. Observed and simulated water levels and drawdowns from aquifer test of the Birch Road wells, east central 

Massachusetts.

[Aquifer test done in 2006 by SEA Consultants, Inc. (from SEA Consultants, Inc., 2008). Observed and simulated levels in North American Vertical Datum 

1988. Well location shown in figure 15]

Observation 

well

Distance to nearest 

pumping well 

(feet)

Water level at start of test 

(feet)

Water level at end of test 

13 days later 

(feet)

Drawdown 

(feet)

Observed Simulated Observed Simulated Observed Simulated Error

7-90 693.6 120.3 123.6 119.0 120.9 1.3 2.7 -1.5

MW-1 630.4 128.3 126.6 122.4 124.2 5.9 2.4 3.5

SEA-10 667.1 127.9 132.6 121.8 130.3 6.1 2.3 3.8

SEA-11 936.7 129.5 130.4 126.6 127.6 2.9 2.7 0.2

SEA-15 1,193.8 117.7 121.9 117.3 120.6 0.4 1.3 -0.9

SEA-3 594.2 122.5 128.0 121.6 126.1 0.9 1.9 -1.0

SEA-9 194.6 127.6 128.5 117.2 115.6 10.5 13.0 -2.5

Means 4.0 3.8 0.2

Hydrographs of 4 of the 13 observation wells with the most 

frequent measurements, 2 near Lake Cochituate and 2 near the 

well site, indicate reasonable matches to seasonal variability 

by the transient simulations (fig. 16). The average range of 

groundwater level variation at a well is simulated to be 1.9 ft, 

while the observed average range is 2.6 ft. The difference 

between the average simulated and observed range is considered acceptable because simulated levels are monthly averages 

and are expected to have less variation than observed instantaneous levels. 

Simulation of the 2006 aquifer test produced groundwater 

drawdowns that match observations relatively well—water 

levels were oversimulated at four observation wells and undersimulated at three observation wells (table 10 and fig. 17). 

Simulated mean drawdown is 3.8 ft as compared to observed 

mean drawdown of 4.0 ft. Mean simulated drawdown error is 

0.2 ft, and the mean absolute drawdown error is 1.9 ft. There 

is no consistent correlation of simulated drawdown error with 

the distance from observation well to the nearest pumped well. 

However, there are spatial patterns in the errors, with wells to 

the south and southwest of the pumping wells generally having drawdowns simulated too low and wells to the north and 

northeast generally simulated too high. These patterns may be 

related to the complex structure of fine and coarse sediment 

deposits in this area. The simulated rate of drawdown also 

reasonably matched the observed (fig. 17).

These results indicate that the transient model can 

reproduce the aquifer test results reasonably well, but that 

local variability in hydraulic conductivity and perhaps storage 

coefficients could be improved in the model. Local variability 

patterns in hydraulic conductivity, particularly within 

about 1,500 ft to the south and west of Pod Meadow Pond, 

seem to be causing the spatial patterns of head errors in the 

steady-state simulations (fig. 14) and aquifer test simulations 

(table 10).

Model Uncertainty and Sensitivity Analysis

As with all models, uncertainties in the Framingham 

groundwater model are inherent because it is a simplification of a complex natural system. Knowledge of which model 

parameters contribute most to model uncertainty are useful 

for qualifying model results and the soundness of conclusions drawn from those results and for identifying future data 

needed to reduce uncertainties. With the research objectives in 

mind, two measures were chosen for assessing model uncertainty and parameter sensitivity—the rate of induced recharge 

from Lake Cochituate into the aquifer, and the timing of 

surface-water response to Birch Road well pumping. 

Sensitivity tests included 28 model runs that varied hydraulic model parameters over likely value ranges 

(table 11). Model parameters were modified individually, 

with calibrated values multiplied by factors of 0.5, 2, or 10, 

depending on their expected range. In some cases, the range of 

possible parameter values was limited by the model’s ability to 

converge to a solution, in which case the parameter range was 

narrowed to achieve convergence. For each scenario, changes 

in induced recharge rate and surface water response time were 

calculated (table 12).
Groundwater Flow Model 33

March June Sept. Dec. March June Sept.

125

129

1978 1979

133

137

141

145

F1W-84 WKW-119

Observed

Simulated

115

119

123

127

131

135

1985 1987 1989 1991 1993

Water level, in feet 

Year

MW-1 MW-8

Observed

Simulated

Water level, in feet 

Figure 16. Simulated and observed groundwater hydrographs, under current conditions and no 

pumping at the Birch Road wells, east central Massachusetts.
34 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

115

120

125

130

135

0 2 4 6 8 10 12 14

SEA-3

SEA-9

SEA-10

SEA-11

SEA-15

Observed Simulated

Time, in days since start of aquifer test

Water elevation, in feet

EXPLANATION

Observation well

Figure 17. Transient simulation of a 2006 aquifer test made with a daily time step compared to groundwater observations, 

east central Massachusetts. (Aquifer test by SEA Consultants, Inc., 2008; well locations shown in figure 15)

The sensitivity of induced recharge from Lake Cochituate 

was calculated by using the steady-state model with long-term 

average climatic conditions and a hypothetical constant pumping rate of 4.9 ft3

/s (3.17 Mgal/d) from the Birch Road wells 

(scenarios n40–n55). ZONEBUDGET (Harbaugh, 1990) was 

used to determine the rate of induced recharge, or flux, from 

the high-conductivity cells representing the lake to the surrounding aquifer cells. Variations of induced recharge caused 

by parameter value changes were calculated relative to the 

calibrated model (scenario n33m2). 

The sensitivity of surface-water response time was 

calculated by using transient model simulations with monthly 

average conditions (scenarios tn52–tn69). Response time 

(t

50) was defined as the number of days for 50 percent of 

streamflow depletion to occur following 1 month of pumping 

from the Birch Road wells at a hypothetical rate of 4.9 ft3/s 

(3.17 Mgal/d). The transient model was run for 5 years 

(60 months) under average monthly conditions, and the pumps 

were turned on for 1 month, January of year 4 (month 37). 

After the 1 month, the Birch Road wells were turned off and 

flows continued to be simulated through December of year 

5 (month 60). Streamflow depletion in the Sudbury River 

was calculated at the exit of the model by subtracting each 

scenario’s flow from the corresponding flow in the calibrated 

transient model with no Birch Road pumping (scenario 

tn51m). Because the transient simulations use a 1-month stress 

period, a linear interpolation was made between monthly 

responses to estimate the number of days to reach the t

50 target. 

For example, if 40 percent of total stream depletion was seen 

by the end of January, the month of pumping, and 60 percent 

of total depletion was seen by the end of February, then the 

t

50 was interpolated to occur halfway through February and to 

equal 14 days. Although some slight variation in t

50 could be 

expected if the pumps were run in another month than January 

or under different recharge conditions, the goal here was to 

assess sensitivity. Transient simulations were also run for 

scenarios in which pumps are run in other months and under 

varying recharge conditions.
Groundwater Flow Model 35

Table 12. Changes in induced recharge from Lake Cochituate and streamflow response times resulting from changes in 

hydraulic parameter values, east central Massachusetts.

[Model parameter changes shown in table 11; --, not tested]

Parameter values changed

Change in model results in response to parameter value changes

Induced recharge rate1

(cubic feet per second)

Streamflow response time2

(days)

Calibrated Lower Upper Calibrated Lower Upper

Streambed conductance 1.58 1.55 1.58 8.3 10.6 1.9

Aquifer hydraulic conductivity

Sand and gravel 1.58 -0.33 3.94 8.3 No solution 0

Fines 1.58 1.30 1.92 8.3 7.1 3.7

Dudley Pond gyttja 1.58 1.62 1.37 8.3 7.9 5.8

Bedrock and till at depth 1.58 1.57 1.63 8.3 10.9 0

Anisotropy 1.58 1.58 1.25 8.3 2.6 12.9

Storage coefficients

Specific yield 1.58 -- -- 8.3 0 35.4

Specific storage 1.58 -- -- 8.3 16.6 7.8

1

Change from no pumping; negative values indicate reduced flow from lake to aquifer.

2

Response time defined as time to reach 50 percent of streamflow depletion in the Sudbury River.

Table 11. Model simulations used to determine sensitivity of induced recharge from Lake Cochituate and pumping response times to 

hydraulic parameter values, east central Massachusetts.

[ft-1, inverse feet; --, not tested]

Parameters changed

Calibrated average model parameter 

values and range of values tested

Name of model scenario used to test parameter change

Induced recharge Streamflow response time

Calibrated Lower Upper Calibrated Lower Upper Calibrated Lower Upper

Average streambed conductance1

(feet per day)

70 35 105 n33m2 n40 n50 tn51m tn52 tn62

Average aquifer hydraulic conductivity1

(feet per day)

Sand and gravel 113 85 226 n33m2 n41 n51 tn51m tn53 tn63

Fines 12.5 6.25 25 n33m2 n42 n52 tn51m tn54 tn64

Dudley Pond gyttja 0.01 0.001 0.1 n33m2 n43 n53 tn51m tn55 tn65

Bedrock and till at depth 0.2 0.02 2 n33m2 n44 n54 tn51m tn56 tn66

Anisotropy 10 5 20 n33m2 n45 n55 tn51m tn57 tn67

Storage coefficients changes

Specific yield (percent) 0.15 0.10 0.3 n33m2 -- -- tn51m tn58a tn68

Specific storage (ft-1) 10-5 10-6 10-4 n33m2 -- -- tn51m tn59 tn69

1

Values vary in the model but averages were used to indicate relative change.
36 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

20-day aquifer test. Users should keep the following model 

limitations in mind when using the model to answer questions 

or design water resource management strategies. 

The model simulates Lake Cochituate’s exchange of 

lake water with the aquifer and the effects of groundwater 

pumping on this exchange through the representation of the 

lake as a high-conductivity zone. The model does not simulate 

other parts of Lake Cochituate’s water budget and so, by 

itself, the model cannot simulate lake levels. The model was 

purposefully constructed without using the LAK package 

(Merritt and Konikow, 2000) because Lake Cochituate’s 

hydrology is predominately controlled by inflow from its 

17.5 mi2 upstream watershed and a previous surface water 

model is available to explicitly simulate the lake (Zarriello, 

2010). In addition, the conductivity of bed sediments of Lake 

Cochituate is not well known, which adds uncertainty to the 

estimation of the exchange of water between the lake and the 

aquifer. Future work that combines the groundwater model 

with the surface-water model of Zarriello and others (2010) 

could allow groundwater effects to be included in more 

temporally refined and realistic simulations of streamflows 

and surface water levels. Unlike other ponds in the model 

that are assigned 22 in/yr of recharge, Lake Cochituate is 

assigned zero recharge. If 22 in/yr of recharge were assigned 

to Lake Cochituate, an annual average flow of 0.5 ft3/s would 

be added to simulated baseflow. Additionally, in the transient 

simulations no inflow is assigned to the stream traversing 

Dudley Pond to represent runoff from the pond’s drainage 

area outside the active model area. These small unassigned 

recharge and runoff flows do not change the main findings of 

this study (sources of water to the wells, timing of streamflow 

response to pumping, and the potential for pumping schedules 

to ameliorate low flows). If the groundwater model is 

combined with the surface-water model of Zarriello and 

others (2010) to simulate full streamflow, assignment of these 

flows should be revisited. Use of the LAK package within 

MODFLOW-NWT could help improve simulation of Lake 

Cochituate’s bottom sediments when more data are available 

to describe its hydraulic properties.

The assumption of no groundwater flow across 

the northern and southern boundaries of the model is a 

potential limitation of the model, but few groundwater level 

measurements are available to confirm the assumption. The 

balance of evidence indicates there are likely low rates of 

groundwater flow across these boundaries. This limitation is 

further tempered by the assignment of streams that provide 

an alternate route for water to move across these model 

boundaries. If the model misrepresents no-flow conditions 

at the boundaries, then pumping stress could potentially 

draw water through these boundaries and the model would 

oversimulate the effect of pumping on streamflow. Therefore, 

the no-flow boundaries are conservative, in terms of 

simulating the maximum induced streamflow. 

Assignment of model parameter values, such as hydraulic 

conductivity, Ss, and Sy, could be improved by using automated parameter estimation methods (Poeter and others, 2005; 

Simulation results indicate that induced recharge from 

Lake Cochituate is most sensitive to hydraulic conductivity of 

the sand and gravel materials (induced recharge increases from 

1.58 to 3.98 ft3/s when the average hydraulic conductivity 

of sand and gravel is increased by a factor of 2) and mildly 

sensitive to the hydraulic conductivity of fine sediments and 

anisotropy (table 12). These sensitivities are consistent with 

the hydrogeologic structure of the model in which layer 1 

cells representing Lake Cochituate are underlain by the sand 

and gravel of layer 2 and the fine sediment of layer 3. The 

sensitivity of lake recharge to the hydraulic conductivity of 

sand and gravel indicates that future data collection efforts 

should look more closely at lake bed sediments and the 

sediment units underlying the lake. Information on sediment 

textures and hydraulic conductivities of deposits underlying 

Lake Cochituate could improve model accuracy and lessen 

model uncertainty.

The timing of surface-water response is most sensitive 

to Sy, the storage coefficient for partially saturated aquifer 

cells. Increasing Sy from 15 percent, the calibrated value, 

to 30 percent caused t

50 to increase from 8.3 to 35.4 days. 

Decreasing Sy from 15 to 10 percent caused t

50 to decrease 

from 8.3 to 0 days, meaning that t

50 was reached within 

the same month as pumping (January). However, model 

calculations contained streamflow mass balance errors that 

were large enough to affect the calculation of t

50 for run tn58a, 

so the exact value of t

50 is uncertain. These mass balance 

errors were caused by numerical artifacts that, while not large 

relative to most flows in the model, were large relative to 

the small differences in streamflow used to calculate t

50. The 

timing of surface-water response is also sensitive to hydraulic 

conductivity of sand and gravel (Ksg). Increasing average 

model Ksg from 113 to 226 ft/d causes t

50 to decrease from 8.3 

to 0 days. The effect of decreasing Ksg could not be assessed 

because the model becomes numerically unstable unless Ksg

is close to or above the calibrated values (average 113 ft/d). 

However, decreasing Ksg would likely cause t

50 to increase 

appreciably. The sensitivity of surface-water response time 

to pumping reinforces the finding that better determination of 

the hydraulic conductivity of sand and gravel deposits could 

improve the model accuracy and lessen the model uncertainty. 

This test also indicates that better information on Sy could 

improve the model simulation of surface-water response times 

to pumping.

Model Limitations

The groundwater model is an imperfect, although still 

useful, representation of the study area aquifer. The model, as 

currently configured, is not suited for simulating groundwater 

variability at spatial scales of less than 50 ft or at time scales 

of less than 1 month. Variability at scales below these limits is 

averaged out in the model. However, there are some situations 

where it can be appropriate to modify the model for simulating 

finer spatial scales or time scales, as was done to simulate the 
Simulated Aquifer and Streamflow Response 37

Doherty and others, 2010), which could reduce or remove 

the spatial patterns in head simulation error around the Birch 

Road well site. This limitation should be kept in mind if a user 

wants to make detailed head predictions for this area. Refined 

hydraulic conductivity assignment around the well site would 

also likely reduce the oversimulation of baseflow from the Pod 

Meadow Pond drainage. 

Model sensitivity to bedrock hydraulic conductivity 

perhaps indicates that permeable bedrock is too thick in the 

model. Future modeling efforts should include adjustment of 

bedrock thickness to improve calibration and reduce sensitivity. The model is also sensitive to Sy, which adds uncertainty 

to findings about the timing of stream response to pumping. 

Simulated Aquifer and Streamflow 

Response

The calibrated steady-state and transient groundwater 

models were used to simulate pumping scenarios that address 

the study objectives (table 13). Water balances and basic 

Table 13. 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]

Scenario Hypothetical pumping of Birch Road wells 

(cubic feet per second)

Number Name Jan. Feb. Mar. Apr. May June July Aug. Sep. Oct. Nov. Dec.

Steady-state model scenarios

Average recharge conditions

1 n33m 0

2 n33m2 4.9

Transient model scenarios

Average recharge conditions

3 tn51m 4.9 0 0 0 0 0 0 0 0 0 0 0

4 tn50m7 0 0 0 0 0 0 0 0 0 0 0 0

5 tn71 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9

6 tn72 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 0.1 4.9 4.9 4.9

7 tn73 4.9 4.9 4.9 4.9 4.9 4.9 0.1 0.1 0.1 4.9 4.9 4.9

8 tn74 4.9 4.9 4.9 4.9 0.1 0.1 0.1 0.1 0.1 0.1 4.9 4.9

Dry (low recharge) conditions

9 tn85 0 0 0 0 0 0 0 0 0 0 0 0

10 tn81a 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9

11 tn82a 4.9 4.9 4.9 4.9 4.9 4.9 4.9 4.9 0.1 4.9 4.9 4.9

12 tn83a 4.9 4.9 4.9 4.9 4.9 4.9 0.1 0.1 0.1 4.9 4.9 4.9

13 tn84 4.9 4.9 4.9 4.9 0.1 0.1 0.1 0.1 0.1 0.1 4.9 4.9

features of the hydrologic system were first described from 

steady-state simulations of average conditions (scenarios 1 and 

2). Hypothetical pumping of the Birch Road wells was added 

to the steady-state model to examine the resulting long-term 

changes to the water balance and baseflow. Transient model 

simulations were used to examine the seasonal surface-water 

response to hypothetical pumping of the Birch Road wells 

under average (scenarios 3–8) and dry conditions (scenarios 9–13). Alternative pumping schedules were examined by 

using the transient model to determine how variable pumping 

schedules can reduce the effects of pumping on surface water. 

Dry conditions simulated in scenarios 9–13 are imposed 

by setting recharge to zero for the months of July, August, and 

September. Under monthly average conditions, recharge to the 

aquifer during these three months totals 1.6 in. This is a reduction in total annual recharge of only 7.1 percent, but it occurs 

in the 3-month window before the critical low-flow month 

of September. The dry conditions, particularly in combination with the low stream inflow values that are specified for 

all transient simulations, create very low-flow ecologically 

conservative conditions. 
38 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

stream system, average annual streamflow in the Sudbury 

River at the Saxonville streamgage (01098530) is 205 ft3/s. 

Hence, on average, the hypothetical 4.9 ft3/s withdrawal from 

the Birch Road wells accounts for only a little over 2 percent 

of average daily Sudbury River streamflow. Concern arises 

during low-flow conditions when water use and regulation represent an increasingly large proportion of the total streamflow. 

The potential effects of the additional stress on Sudbury River 

low flows from Birch Road well withdrawals were addressed 

by transient model simulations. 

Pumping induces additional recharge to the aquifer system from Lake Cochituate. Imposing a hypothetical constant 

pumping rate of 4.9 ft3/s (3.17 Mgal/d) at the Birch Road wells 

(scenario 2) induces additional recharge from Lake Cochituate 

at a rate of 1.6 ft3/s, which is about equal to 32 percent of the 

hypothetical pumping rate (table 9). Under constant pumping, Dudley Pond receives slightly less groundwater discharge 

from the aquifer, but more significantly, groundwater discharge to Pod Meadow Pond is almost entirely captured by the 

pumping wells. 

A fundamental principle of the steady-state groundwater 

model is the requirement that inflows must equal outflows; in 

a transient model inflows must equal outflows plus or minus 

changes in storage. A natural consequence of this requirement, 

without a change to the boundary conditions, is that groundwater discharge to streams must decrease by the amount of 

groundwater pumped from wells in steady-state simulations. 

The extent to which a particular stream reach is affected by 

withdrawals depends on its relative location to pumping wells 

and on neighboring hydrologic and boundary conditions. 

Streams upgradient of pumping wells experience little or no 

baseflow reduction, whereas streams close to and downgradient from the pumping wells experience greater reduction. As 

a percentage of the imposed pumping (4.9 ft3/s), the reduction 

of baseflow in streams (scenario 2) is illustrated by the width 

of the stream lines shown in figure 18. Streamflow depletions 

in the Sudbury River, equal or nearly equal to the constant 

4.9 ft3/s pumping rate, occur downstream of the oxbow. The 

stream traversing Pod Meadow Pond (segment 9) experiences 

the greatest percentage of streamflow reduction (100 percent) 

as streamflow goes to zero under the steady-state imposed 

pumping rate of 4.9 ft3/s; this reduction is not visible on 

figure 18 because it falls below the display threshold value of 

1.5 ft3/s.

Transient Model Simulations

Transient model simulations demonstrate the dynamic 

response of the groundwater and surface-water interaction to 

variable pumping rates under seasonal and low-flow recharge 

conditions. All the transient scenarios (scenarios 3–13) take a 

conservative approach to streamflow simulation by specifying inflow to stream segments at the model boundary as the 

75-percent daily flow duration for each month. All transient 

scenarios also simulate each of the Wayland wells pumping 

Table 14. Steady-state simulated water budgets with and 

without hypothetical Birch Road well withdrawals, east central 

Massachusetts.

Water balance

Current 

conditions

(scenario 1)

Birch Road wells 

pumping at 

4.9 cubic feet 

per second

(scenario 2)

Inflows 

(cubic feet per second)

Recharge from precipitation 9.5 9.5

Recharge from streams 2.0 3.6

Total 11.5 13.1

Outflows 

(cubic feet per second)

Pumping 1.4 6.3

Discharge to streams 10.1 6.8

Total 11.5 13.1

Mass balance error 0.0 0.0

Steady-State Simulations

The water balance of inflows and outflows under 

current and hypothetical conditions indicates the sources of 

water needed to satisfy withdrawals under average climatic 

conditions (table 14). Under current conditions (scenario 1), 

the total inflows to and outflows from the model are 11.5 ft3/s; 

83 percent of the inflow is from recharge by precipitation, 

and the remainder (17 percent) is from stream recharge. 

Stream recharge is mostly from river cells that traverse 

Lake Cochituate and lose water to the high-conductivity 

cells representing the lake. Outflow is mostly discharge to 

streams (88 percent), and the remainder (12 percent) is from 

withdrawal by the Wayland production wells. 

Under hypothetical steady-state pumping of the Birch 

Road wells (scenario 2), the total inflows to and outflows from 

the model are 13.1 ft3/s; inflow from recharge by precipitation 

is the same as before, but decreases as a percentage of the total 

inflow (73 percent), whereas recharge from streams increases 

(27 percent). Increased recharge from streams is mostly from 

induced infiltration of stream water from stream cells traversing Lake Cochituate (table 5, fig. 9). Most outflow from the 

model is still discharge to streams (52 percent) but decreases 

from 10.1 to 6.8 ft3/s compared to no pumping of the Birch 

Road wells in scenario 1. Withdrawals from the Wayland and 

the Birch Road wells account for 11 and 37 percent of the total 

outflow from the model under average steady-state conditions, respectively. It should be noted that although pumping 

intercepts groundwater that would otherwise discharge to the 
Simulated Aquifer and Streamflow Response 39

Figure 18. Baseflow reduction (streamflow depletion) in stream-channel segments in response to simulated pumping 

at the Birch Road wells at a constant rate of 4.9 cubic feet per second (3.17 million gallons per day), east central 

Massachusetts.

90

Sudbury River 

Sudbury River 

Sudbury River 

Dudley 

Pond

Lake Cochituate

Heard Pond

Bedrock island Bedrock island

Pod Meadow 

Pond

EXPLANATION

Streamflow depletion in 

cubic feet per second from 

steady-state pumping at 

4.9 cubic feet per second

1.5 to < 2.0

> 4.7 

2.0 to 2.5

Wetland

Wayland 

production well

Inactive

model area

Birch Road wells

Simulated stream cell

Oxbow

Pod Meadow

Birch Road wells

0 1,000

Base from U.S. Geological Survey and Massachusetts Geographic Information System data sources,

Massachusetts State Plane Coordinate System, Mainland Zone

2,000 FEET

0 250 500 METERS

Cochituate Brook

Sudbury

River
40 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

at the average reported monthly rate. When applicable to the 

scenario, pumping of the Birch Road wells was distributed 

among the four wells developed for the 2006 aquifer test by 

SEA Consultants, Inc. (2008).

Under monthly average conditions, with no pumping 

of Birch Road wells (scenario 4), water is lost from aquifer 

storage from May through September and replenished from 

October through April (fig. 19). A few stream segments, 

primarily the stream traversing Lake Cochituate, provide 

recharge to the aquifer throughout the year from a maximum 

of 29 percent in September to a minimum of 10 percent in 

December, when recharge from precipitation is at its normal 

annual lowest and highest levels, respectively. Discharge to 

streams is the primary outflow from the aquifer, although 

in September pumping from Wayland wells accounts for 

18 percent of total outflow from the aquifer when discharge 

from the aquifer to streams is at its seasonal low. Under 

average monthly recharge conditions, simulated discharge 

from the aquifer to streams is lowest during August and 

September and highest in April and May (fig. 19). 

Evaluation of Pumping Strategies to Reduce 

Low-Flow Stresses

The transient model was used to determine the timing of surface-water response to pumping and to evaluate 

pumping strategies to reduce the effects of pumping during 

low-flow conditions. The timing of surface-water response 

to pumping is characterized by the results from scenario 3, 

in which the Birch Road wells are simulated as pumping at 

4.9 ft3/s for 1 month, and subsequent streamflow depletions 

are determined (fig. 20). Associated streamflow depletions 

start in the same month as pumping and dissipate quickly after 

the pumping stops because of the high permeability of the 

aquifer and the proximity of streams to the pumping wells. As 

soon as pumping stops, the reduction of baseflow, or in other 

-25

-20

-15

-10

-5

0

5

10

15

20

25

Jan. Feb. March April May June July Aug. Sept. Oct. Nov. Dec.

Flow, in cubic feet per second

 Outflows Inflows

EXPLANATION

Inflows

From streams

From precipitation

Outflows

To existing well withdrawals

To streams

From storage (water table falling) To storage (water table rising)

Monthly groundwater-model flows

Note: As storage decreases (water table falls) it creates inflow to the active groundwater system, and 

as storage increases (water table rises), it creates an outflow from the active groundwater system. 

Inflows from precipitation are net, accounting for evapotranspiration losses.

Month

Figure 19. Simulated transient monthly water budget under current average monthly conditions (no pumping of Birch 

Road wells), east central Massachusetts.
Simulated Aquifer and Streamflow Response 41

words streamflow depletion, decreases by about 40 percent. 

Within 2 months, the reduction of baseflow decreases by about 

80 percent, and within 4 months, by about 90 percent (fig. 20). 

The fast response of surface water to pumping stresses 

indicates the potential for pumping to be managed to reduce 

streamflow depletions during periods of low flow.

Alternate pumping schedules and rates were simulated 

with the transient model (table 13) to evaluate possible 

strategies to reduce the impacts of pumping on surface 

water during average seasonal low flows (scenarios 4–8) 

and during more extreme dry conditions (scenarios 9–13). 

Scenarios 4 and 9 set baseline conditions for average and 

dry climatic conditions, respectively, in which there was no 

simulated pumping at the Birch Road wells. Scenarios 5 and 

10 simulated the Birch Road wells pumping at a constant rate 

of 4.9 ft3/s. Scenarios 6–8 and 11–13 simulated a reduced 

pumping rate of 0.1 ft3/s for 1, 3, and 6 months of the year 

and a pumping rate of 4.9 ft3/s in the remaining months. A 

minimum pumping rate of 0.1 ft3/s was used because the 

pumped water requires treatment, and this is the estimated 

minimum pumping rate required to keep a treatment facility 

operational without a major shutdown and startup (Peter 

Newton, Bristol Engineering Advisors, Inc., written commun., 

December 2011). 

Alternative pumping schedules were chosen as examples 

of how these rates could affect the usual seasonal cycle of 

low summertime streamflows. In actual operation, pumping 

rates could be determined from streamflow, weather forecasts, 

upstream reservoir operations, or other criteria, but for this 

analysis reductions in pumping the Birch Road wells were 

limited to 1, 3, and 6 months of the year. Simulated monthly 

flows in the Sudbury River at the model exit under average 

monthly climatic conditions and various pumping scenarios 

0

20

40

60

80

100

0

1

2

3

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

Percentage of total streamflow

Month

Percentage of total

Streamflow reduction, in cubic feet per second

Streamflow reduction

Figure 20. 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), east central Massachusetts.

(scenarios 5–8) show small changes between pumping and 

no-pumping scenarios during normal seasonal high recharge 

(December through May); however, changes are noticeable 

by September (table 15, fig. 21). The greatest reduction 

in streamflow occurs from July through September when 

pumping continues throughout most of the year (scenarios 

5 and 6). Reduced pumping for 3 and 6 months of the year 

(scenarios 7 and 8) keeps baseflow nearly the same as without 

pumping. Simulated baseflow under dry recharge conditions 

is only slightly less than simulated baseflow under average 

recharge conditions.

Streamflow depletion increases going downstream in 

the Sudbury River as the effects of pumping accumulate 

(fig. 22). At approximately 1.0 mi downstream from the model 

entrance, where drainage from the Pod Meadow area enters 

the Sudbury River, effects of pumping noticeably increase, and 

streamflow depletion nearly reaches its full extent.

Reducing the pumping rate in September and in prior 

months decreases streamflow depletion compared to constant pumping (fig. 23). Under average recharge conditions, 

reducing pumping for 1 month lowers September streamflow 

depletion from 4.4 to 3.1 ft3/s (scenarios 5 and 6), reducing 

pumping for three months (scenario 7) lowers depletion to 

0.8 ft3/s, and reducing pumping for 6 months (scenario 8) 

drops depletion to 0.3 ft3/s. In terms of percentages, lowering 

maximum pumping rates to near 0 during September reduces 

September streamflow depletion by 29 percent, and lowering 

pumping rates for 3 months (July through September) reduces 

streamflow depletion during September by 79 percent, as 

compared to constant pumping. Under dry conditions (low 

recharge), reducing year-round pumping by 1, 3, and 6 months 

(scenarios 11, 12, and 13) lowers streamflow depletion from 

4.2 to 3.0, 1.1, and 0.5 ft3/s, respectively. 
42 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Table 15. Simulated monthly flow in the Sudbury River under different hypothetical pumping scenarios for the Birch Road wells, east 

central Massachusetts.

[Scenarios 1, 4, and 9 represent simulations with no Birch Road pumping. Yellow shading highlights months of maximum pumping at a rate of 4.9 cubic feet 

per second]

Scenario Sudbury River flow 

(cubic feet per second)

Number Name Jan. Feb. Mar. Apr. May June July Aug. Sep. Oct. Nov. Dec.

Steady-state model scenarios

Average recharge conditions

1 n33m 128.8

2 n33m2 123.9

Transient model scenarios

Average recharge conditions

3 tn51m 134.6 147.4 226.7 205.9 128.6 64.4 36.2 26.3 23.1 42.2 82.1 128.9

4 tn50m7 136.8 148.6 227.4 206.2 128.8 64.5 36.3 26.4 23.2 42.3 82.1 129.0

5 tn71 131.7 143.5 222.2 201.1 123.9 59.8 31.7 21.9 18.8 37.3 77.1 123.8

6 tn72 132.0 143.7 222.4 201.3 124.0 59.9 31.8 21.9 20.1 38.4 77.8 124.3

7 tn73 132.3 143.9 222.5 201.4 124.0 59.9 33.5 25.0 22.3 39.5 78.5 124.7

8 tn74 132.4 144.0 222.6 203.4 127.1 63.5 35.7 25.9 22.9 39.8 78.8 124.9

Dry (low recharge) conditions

9 tn80 136.3 148.2 227.0 206.0 128.6 64.4 35.0 25.4 21.7 41.3 81.4 128.3

10 tn81 131.3 143.2 222.0 200.9 123.5 59.5 30.4 20.8 17.5 36.4 76.5 123.3

11 tn82 131.6 143.4 222.2 201.0 123.6 59.6 30.4 20.8 18.7 37.6 77.2 123.8

12 tn83 131.9 143.6 222.3 201.1 123.6 59.6 32.4 23.8 20.6 38.5 77.8 124.2

13 tn84 132.0 143.7 222.4 203.1 126.7 63.3 34.3 24.8 21.2 38.9 78.0 124.4

Figure 21. Simulated monthly streamflows in the Sudbury River at the model exit under five hypothetical pumping 

scenarios for the Birch Road wells, east central Massachusetts.

10

100

1,000

Simulated streamflow, in cubic feet per second

No pumping (scenario 4)

6 months (scenario 8)

9 months (scenario 7)

11 months (scenario 6)

12 months (scenario 5)

Pumping at 4.9 cubic feet per 

second for periods of

Jan. Feb. Mar. April May June July Aug. Sept. Oct. Nov. Dec.
Simulated Aquifer and Streamflow Response 43

-20

-18

-16

-14

-12

-10

-8

-6

-4

-2

0

Jan. Feb. March April May June July Aug. Sept. Oct. Nov. Dec.

Constant pumping (4.9 ft3/s) 

Average Dry

Recharge condition

Period of reduced pumping (0.1 ft3/s) 

Pumping at 4.9 ft3/s for periods of

12 months (tn81)

(tn82)

(tn83)

(tn84)

(tn71)

(tn72)

(tn73)

(tn74)

Scenario name in parentheses; ft3

/s, cubic feet per second 

1

10

100

1,000

Average

Dry 

Sudbury River 

discharge under 

average and dry 

recharge conditions

Percent reduction in streamflow

Streamflow, in cubic feet per second

11 months

9 months

6 months

Month

Figure 23. Percent reduction in Sudbury River simulated streamflow at the model exit in response to pumping under average 

and dry recharge conditions, east central Massachusetts. Streamflow depletion was determined relative to no pumping of the 

Birch Road wells.

0

10

20

30

0 1 2 3 4

No pumping (scenario 4)

Pumping at 4.9 cubic feet per 

second for periods of

6months (scenario 8)

9months (scenario 7)

11months (scenario 6)

12months (scenario 5)

Simulated streamflow, in cubic feet per second

River mile downstreamfrommodel boundary

Figure 22. Simulated streamflow along the Sudbury 

River for the month of September under average 

recharge rates and various pumping durations at the 

Birch Road wells, east central Massachusetts.
44 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

Under dry conditions, pumping causes a slightly larger 

percent decrease in streamflow than it does during average 

recharge conditions (fig. 23). Relative to streamflow under no 

pumping of the Birch Road wells, constant pumping of these 

wells at a combined rate of 4.9 ft3/s decreases streamflow 

by 18.8 and 19.2 percent under average conditions and dry 

recharge conditions, respectively (blue lines in fig. 23). 

When pumping is reduced to near zero during September 

only, streamflow depletion is 13.3 and 13.7 percent of total 

streamflow under average conditions and under dry conditions, 

respectively (purple lines in fig. 23). The streamflow 

reduction decreases further when pumping is reduced for 

3 and 6 months, decreasing to 3.9 and 4.8 percent and 1.3 

and 2.2 percent under average and dry recharge conditions, 

respectively. It should be noted that the percent reduction in 

streamflow is based on assigned inflow values of 75 percent 

daily flow duration for each month, representing low-flow 

conditions. Under average monthly flow conditions, the 

percent streamflow reductions would be substantially less. 

Similar patterns of reduced streamflow depletion could be 

achieved during other periods if reduced pumping rates were 

applied at other times.

Summary and Conclusions

The Sudbury River Basin in eastern Massachusetts is an 

ecologically important resource housing the Great Meadows 

National Wildlife Refuge, one of the eight National Wildlife 

Refuges in the State. The Sudbury River is considered stressed 

by urban development, water withdrawals, and reservoir controls, particularly during periods of low flow when the stream 

ecology is the most vulnerable to water-management operations. During summer months and periods of low precipitation, the river is primarily fed by groundwater. Because of 

the river’s close hydraulic connection with aquifers, groundwater withdrawals can reduce streamflow during low-flow 

periods. The Town of Framingham has proposed reactivating 

groundwater-supply wells along its northern border, near the 

Sudbury River and the adjacent towns of Wayland, Sudbury, 

and Natick. This proposal has raised concerns that these withdrawals may further reduce streamflows and adversely affect 

nearby State and Federal conservation areas and surface-water 

bodies, particularly Lake Cochituate and Great Meadows 

National Wildlife Refuge. 

In response to these concerns, the U.S. Geological 

Survey, in cooperation with the Town of Framingham, 

undertook this investigation to improve the understanding of 

the hydrogeology of the local aquifer system and the potential 

effects of the proposed pumping on nearby surface-water 

features. The study also examined whether groundwater 

pumping could be managed to minimize the effects of 

withdrawals during critical low-flow periods. Goals of the 

investigation included improving understanding of hydrology 

in the study area, determining rates of flow of water from 

Lake Cochituate to the Birch Road wells, evaluating potential 

effects of pumping on surface waters, and assessing the 

potential for managing pumping to reduce stresses on the 

hydrologic system. 

A numerical groundwater-flow model was developed for 

the study to simulate the hydrology of the glacial-sediment 

aquifer in northeastern Framingham and adjacent towns of 

Wayland, Sudbury, and Natick, by using MODFLOW-NWT. 

This model was chosen because it provides greater numerical 

stability than previous versions of MODFLOW for this type of 

hydrogeologic setting, where model cells are subject to drying 

because of their small saturated thickness and a fluctuating 

water table. The model is calibrated with geologic and hydrologic data compiled from prior studies and new data collected 

during this study. Simulated groundwater levels and streamflows have reasonably good agreement with observed values 

under various climatic and groundwater pumping stresses. 

Steady-state and transient simulations reveal details about 

the effects of proposed pumping on groundwater and surface 

water in the vicinity of the Birch Road wells: 

• Pumping the Birch Road wells captures groundwater 

from the surrounding aquifer and induces additional 

recharge from Lake Cochituate. Under constant 

(steady-state) pumping, the Birch Road wells induce 

recharge from Lake Cochituate at a rate of 1.6 ft3/s, 

which is equal to about 32 percent of the simulated 

4.9 ft3/s pumping rate.

• Groundwater withdrawals reduce flow in the Sudbury 

River and tributary streams. The Sudbury River downstream of the oxbow is depleted at a rate about equal 

to the rate of pumping of the Birch Road wells under 

steady-state conditions. 

• Streams respond quickly to changes in pumping. When 

the Birch Road wells are pumped for 1 month and 

then stopped, streamflow depletions decrease by about 

80 percent within 2 months and by about 90 percent 

within about 4 months.

• The fast response of surface water to pumping 

stresses provides the potential to substantially reduce 

streamflow depletions during periods of low flow by 

altering pumping rates appropriately for seasonal or 

anticipated flow conditions. Streamflow depletion 

during September, typically the month of lowest flow, 

could be reduced by 29 percent by lowering maximum 

pumping rates to near zero during September. 

Lowering pumping rates to near zero for 3 months 

(July through September) reduces streamflow depletion 

during September by 79 percent as compared to 

constant pumping. 

• Sensitivity analysis of surface-water response times 

and rates of induced recharge from Lake Cochituate 

suggest that model uncertainty could be reduced by 

better knowledge of the spatial distribution and values 

of specific yield and hydraulic conductivity of the sand 

and gravel sediments.
References Cited 45

Doherty, J.E., Hunt, R.J., and Tonkin, M.J., 2010, Approaches 

to highly parameterized inversion: A guide to using PEST 

for model-parameter and predictive-uncertainty analysis: 

U.S. Geological Survey Scientific Investigations Report 

2010–5211, 71 p. (Also available at http://pubs.usgs.gov/

sir/2010/5211.)

Friesz, P.J., and Church, P.E., 2001, Pond-aquifer interaction 

at South Pond of Lake Cochituate, Natick, Massachusetts: 

U.S. Geological Survey Water-Resources Investigations 

Report 01–4040, 42 p. (Also available at http://pubs.usgs.

gov/wri/wri014040/.)

Gay, F., 1981, Hydrologic data of the Lake Cochituate 

drainage basin, Framingham-Natick, Massachusetts: U.S. 

Geological Survey Open-File Report 82–342, 70 p. (Also 

available at http://pubs.er.usgs.gov/publication/ofr82342/.)

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/.)

Gleeson, T., Alley, W.M., Allen, D.M., Sophocleous, M.A., 

Zhou, Y., Taniguchi, M., and VanderSteen, J., 2011, 

Towards sustainable groundwater use—setting long-term 

goals, backcasting, and managing adaptively: Ground 

Water, v. 50, no. 1, p. 19–26, accessed September 1, 2012, 

at http://onlinelibrary.wiley.com/doi/10.1111/

j.1745-6584.2011.00825.x/full.

Goldsmith, Richard, 1991, Stratigraphy of the MilfordDedham Zone, eastern Massachusetts—An Avalonian 

Terrane, in Hatch, N.L., Jr., ed., 1991, The bedrock geology 

of Massachusetts: U.S. Geological Survey Professional 

Paper 1366–E, 62 p. (Also available at http://pubs.er.usgs.

gov/publication/pp1366EJ.)

GZA GeoEnvironmental, Inc., Jacobs Associates, and 

Sverdrup Civil, Inc., 1995, MWRA contracts 6054, 6055, 

and 6059, Geotechnical data report, Metrowest watersupply tunnel, construction packages 1, 2, and 3: v. III of 

IX, variously paged.

Haley & Aldrich, 1996, Response action outcome statement, 

1455 Concord Street, Framingham, Massachusetts: 

Cambridge, Mass., Haley & Aldrich, Inc., Mass DEP 

RTN-3-13162, 183 p., accessed September 1, 2012, at 

http://public.dep.state.ma.us/fileviewer/Rtn.

aspx?rtn=3-0013162.

Hanson, B.P., and Simcox, A.C., 1994, Yields of bedrock wells 

in Massachusetts: U.S. Geological Survey Water-Resources 

Investigations Report 93–4115, 43 p. (Also available at 

http://pubs.er.usgs.gov/publication/wri934115/.)

Model simulations made thus far have improved the 

understanding of groundwater/surface-water interactions in 

the study area. Additional data for calibrating the model could 

further improve the understanding of this complex system and, 

along with model simulations designed to achieve the most 

effective pumping strategies, could substantially reduce the 

effects of withdrawals on surface waters in the area.

References Cited

Balsam Environmental Consultants, Inc., 1986, Preliminary 

site investigation report, New England Sand and Gravel 

Company, Saxonville (Framingham), Massachusetts: 

Salem, N.H., Balsam Environmental Consultants, Inc., 

123 p., accessed September 1, 2012, at http://public.dep.

state.ma.us/fileviewer/Rtn.aspx?rtn=3-0000629.

Balsam Environmental Consultants, Inc., 1987, Volume I 

and II site investigation report, New England Sand and 

Gravel Company property Saxonville (Framingham), 

Massachusetts: Salem, N.H., Balsam Environmental 

Consultants, Inc., 115 p. and 215 p., accessed September 1, 

2012, at http://public.dep.state.ma.us/fileviewer/

Rtn.aspx?rtn=3-0000629.

Balsam Environmental Consultants, Inc., 1992, Additional 

hydrogeologic investigation/ground water modeling, New 

England Sand and Gravel site, Saxonville, Massachusetts, 

Volumes I and II: Salem, N.H., Balsam Environmental 

Consultants, Inc., 123 p., accessed September 1, 2012, 

at http://public.dep.state.ma.us/fileviewer/

Rtn.aspx?rtn=3-0000629.

Boutt, D.F., Diggins, P., and Mabee, S.B., 2010, A field study 

(Massachusetts, USA) of the factors controlling the depth 

of groundwater flow systems in crystalline fractured-rock 

terrain, Hydrogeology Journal, v. 18, no. 8, p. 1839–1854, 

DOI: 10.1007/s10040-010-0640-y.

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., 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/.)

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. (Also available at http://pubs.er.usgs.gov/publication/

wri024234.)
46 Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts

MassGIS, 2003, Digital orthophoto elevation, Commonwealth 

of Massachusetts Executive Office for Administration and 

Finance: accessed November 22, 2010, at 

http://www.mass.gov/mgis/p.htm.

Massachusetts Department of Conservation and Recreation (MassDCR), 2006, Public access management plan 

update—Quabbin Reservoir watershed system: 18 p., 

accessed March 22, 2012, at http://www.mass.gov/dcr/

watersupply/watershed/documents/quabbinaccessch1.pdf.

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: U.S. Geological Survey WaterResources Investigations Report 00–4167, 146 p. (Also 

available at http://pubs.er.usgs.gov/publication/wri004167.)

Nelson, A.E., 1974a, Bedrock geologic map of the 

Natick quadrangle, Middlesex and Norfolk Counties, 

Massachusetts: U.S. Geological Survey Geologic 

Quadrangle Map GQ–1208, scale 1:24,000, accessed 

September 1, 2012, at http://ngmdb.usgs.gov/Prodesc/

proddesc_10746.htm.

Nelson, A.E., 1974b, Surficial geologic map of the 

Framingham quadrangle, Middlesex and Worcester 

Counties, Massachusetts: U.S. Geological Survey Geologic 

Quadrangle Map GQ–1176, scale 1:24,000, accessed 

September 1, 2012, at http://ngmdb.usgs.gov/Prodesc/

proddesc_10691.htm.

Nelson, A.E., 1974c, Surficial geologic map of the 

Natick quadrangle, Middlesex and Norfolk Counties, 

Massachusetts: U.S. Geological Survey Geologic 

Quadrangle Map GQ–1151, scale 1:24,000, accessed 

September 1, 2012, at http://maps.ngmdb.us/dataviewer/.

Nelson, A.E., 1975, Bedrock geologic map of the Framingham 

quadrangle, Middlesex and Worcester Counties, Massachusetts: U.S. Geological Survey Geologic Quadrangle Map 

GQ–1274, scale1:24,000, accessed September 1, 2012, at 

http://maps.ngmdb.us/dataviewer/.

Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, MODFLOW-NWT, A Newton formulation for MODFLOW-2005: 

U.S. Geological Survey Techniques and Methods 6–A37, 

44 p. (Also available at http://pubs.er.usgs.gov/publication/

tm6A37/.)

Harbaugh, A.W., 1990, A computer program for calculating 

subregional water budgets using results from the U.S. 

Geological Survey modular three-dimensional ground-water 

flow model: U.S. Geological Survey Open-File Report 

90–392, 46 p. (Also available at http://pubs.er.usgs.gov/

publication/ofr90392/.)

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW-2000, the U.S. Geological Survey 

modular groundwater model—User guide to modularization 

concepts and the ground-water flow process: U.S. Geological Survey Open-File Report 00–92, 121 p. (Also available 

at http://pubs.er.usgs.gov/publication/ofr200092/.)

Hill, M.C., Banta, E.R., Harbaugh, A.W., and Anderman, 

E.R., 2000, MODFLOW-2000, the U.S. Geological Survey 

modular ground-water model—User guide to the observation, sensitivity, and parameter-estimation processes and 

three post-processing programs: U.S. Geological Survey 

Open-File Report 00–184, 210 p. (Also available at 

http://pubs.er.usgs.gov/publication/ofr00184/.)

Hunt, R.J., Haitjema, H.M., Krohelski, J.T., and Feinstein, 

D.T., 2003, Simulating ground water–lake interactions—

Approaches and insights: Ground Water, v. 41, no. 2, 

p. 227–237.

Ibs-von Seht, Malte, and Wohlenberg, Jtirgen, 1999, Microtremor measurements used to map thickness of soft sediments: Bulletin of the Seismological Society of America, 

v. 89, p. 250–259.

IEP, Inc., 1983, Diagnostic feasibility study Dudley Pond 

Wayland, Massachusetts, Town of Wayland Massachusetts: 

Surface Water Quality Study Committee, April 1983, 145 p.

Kennedy, E.J., 1983, Techniques of water-resources investigations of the United States Geological Survey: chap. Al3, 

book 3, Computation of continuous records of streamflow, 

p. 47. (Also available at http://pubs.er.usgs.gov/publication/

twri03A13/.)

Koteff, Carl, and Pessl, Fred, 1981, Systematic ice retreat in 

New England: U.S. Geological Survey Professional Paper 

1179, 20 p. (Also available at http://pubs.er.usgs.gov/

publication/pp1179/.)

Lane, J.W., White, E.A., Steele, G.V., and Cannia, J.C., 2008, 

Estimation of bedrock depth using the horizontal-to-vertical 

(H/V) ambient-noise seismic method, in Symposium on the 

Application of Geophysics to Engineering and Environmental Problems, April 2008, Philadelphia, Pa., Proceedings: 

Denver, Colo., Environmental and Engineering Geophysical 

Society, 13 p.

Mabee, S.B., Curry, P.J., and Hardcastle, K.C., 2002, Correlation of lineaments to groundwater inflows in a bedrock 

tunnel: Ground Water, v. 40, no. 1, p. 37–43.
References Cited 47

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 6–A13, 50 p. 

(http://pubs.er.usgs.gov/publication/tm6A13)

Poeter, E.P., Hill, M.C., Banta, E.R., Mehl, Steffen, and 

Christensen, Steen, 2005, UCODE_2005 and six other 

computer codes for universal sensitivity analysis, 

calibration, and uncertainty evaluation: U.S. Geological 

Survey Techniques and Methods 6–A11, 283 p.

SEA Consultants, Inc., 1992, Report on prolonged pumping 

test for the Concord/Sudbury river basin Cochituate 

well supply-water supply feasibility study, Framingham, 

Massachusetts: SEA Consultants, Inc., Cambridge, Mass., 

38 p.

SEA Consultants, Inc., 2008, Source final report, Birch Road 

well re-activation, Framingham, Massachusetts: SEA Consultants, Inc., Cambridge, Mass., 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: SEA Consultants, Inc., Cambridge, 

Mass., 201 p.

Sovereign Consulting, Inc., 2009, Phase I initial site investigation, tier classification, and phase II scope of work—Mass. 

DEP RTN 3-27985: Sovereign Consulting, Inc., Mansfield, 

Mass., accessed May 2011, at http://public.dep.state.ma.us/

fileviewer/Default.aspx?formdataid=0&documentid=51762.

Stone, J.R., and Stone, B.D., 2006, Surficial geologic map of 

the Clinton-Concord-Grafton-Medfield quadrangle area in 

east central Massachusetts: U.S. Geological Survey 

Open-File Report 2006–1260–A, 1 pl. (Also available at 

http://ngmdb.usgs.gov/Prodesc/proddesc_80659.htm.)

URS Corporation, 2003, Response action outcome 

statement, New England Sand and Gravel site, Saxonville, 

Massachusetts, RTN 3-0629, HQ AFCEE: URS 

Corporation, Portland, Maine, 62 p.

Wayland Wellhead Protection Committee, and Young, 

B.W., 2011, Wayland wellhead protection plan: Wayland, 

Mass., Wayland Wellhead Protection Committee, 86 p., 

accessed September 1, 2012, at http://www.wayland.ma.us/

Pages/WaylandMA_DPW/WellheadProtectionPlan6MB2June2011.pdf.

Wilmott, C.J., 1977, WATBUG–A Fortran IV algorithm for 

calculating the climatic water budget: Climatology, v. 30, 

p. 2–4.

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/.)

Zen, E-an, Goldsmith, Richard, Ratcliffe, N.M., Robinson, 

Peter, Stanley, R.S., Hatch, N.L., Shride, A.F., Weed, 

E.G.A., and Wones, D.R., 1983, Bedrock geologic map of 

Massachusetts: U.S. Geological Survey, scale 1:250,000, 

accessed September 1, 2012, at http://ngmdb.usgs.gov/

Prodesc/proddesc_16357.htm.
THIS PAGE INTENTIONALLY LEFT BLANK
Prepared by the Pembroke Publishing Service Center.

For more information concerning this report, contact:

Director

U.S. Geological Survey

Massachusetts-Rhode Island Water Science Center

10 Bearfoot Road

Northborough, MA 01532

dc_ma@usgs.gov

or visit our Web site at:

http://ma.water.usgs.gov
Printed on recycled paper

Eggleston and others—Simulation of Groundwater and Surface-Water Interaction in a Glacial-Sediment Aquifer, Massachusetts—SIR 2012–5172