Simulation of Groundwater and Surface-Water Interaction and Effects of Pumping in a Complex Glacial-Sediment Aquifer, East Central Massachusetts

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

Elevation, in feet above NGVD

West

East

Layer 2

Layer 1

Layer 4

Layer 5

Sudbury River

Lake Cochituate

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 Surface-

Water 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 MODFLOW-

NWT 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

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

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 (mi

) 2.590 square kilometer (km

)

Flow rate

cubic foot per day (ft

/d) 0.00001157 cubic foot per second (ft

/sec)

cubic foot per second (ft

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

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

/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-ll 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 simu-

lated 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 baseow 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 ow downstream in the

Sudbury River decreased at the same rate as pumping. How-

ever, surface water responded quickly to pumping stresses.

When pumping was simulated for 1 month and then stopped,

streamow 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 streamow depletions during

periods of low ow, which are of greatest concern to the eco-

logical integrity of the river. Results indicate that streamow

depletion during September, typically the month of lowest

ow, can be reduced by 29 percent by lowering the maxi-

mum pumping rates to near zero during September. Lowering

pumping rates for 3 months (July through September) reduces

streamow depletion during September by 79 percent as com-

pared to constant pumping. These results demonstrate that a

seasonal or streamow-based groundwater pumping schedule

can reduce the effects of pumping during periods of low ow.

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 streamows. 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 con-

nection 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 surface-

water features recently gained attention when the Town of

Framingham (g. 1) sought to reactivate production wells (the

Birch Road wells) that previously provided a local water sup-

ply 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

led 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 Ofce of the Solicitor for the U.S.

Department of the Interior that specically 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.

Sudbury River Sudbury River

Sudbury River

Dudley Pond

Lake Cochituate

Heard Pond

Pod Meadow

Pond

Great Meadows

National

Wildlife

Refuge

Great Meadows

National

Wildlife

Refuge

Bedrock islandBedrock island

WAYLAND

FRAMINGHAM

NATICK

SUDBURY

WAYLAND

FRAMINGHAM

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 wellsBirch 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

00'

30'

00'

Study Area

MASSACHUSETTS

BOSTON

50 KILOMETERS

50 MILES

Introduction 3

for migratory birds, sh, 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-ow periods, by develop-

ment pressures and by numerous groundwater and surface-

water withdrawals that affect streamows (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 pump-

ing. The purpose of this study by the U.S. Geological Survey

(USGS), made in cooperation with the Town of Framing-

ham, 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 natural-

resource 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) ground-

water 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 sched-

ules that could reduce the effects of additional groundwater

withdrawals on streams during low ows when the effects of

withdrawals are most pronounced.

Study Area

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

central Massachusetts (g. 1). The active groundwater model

area covers about 6.1 square miles (mi

) surrounding the

Birch Road wells in the towns of Framingham (1.7 mi

Wayland (4.0 mi

), Sudbury (0.4 mi

), and Natick (0.01 mi

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 (g. 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 mi

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 stratied 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 ne sands, which are the primary aquifer deposits.

Previous Investigations

The hydrology of the Lake Cochituate area was described

by Gay (1985). Surcial geology of the study area was rst

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 g. 1) was investigated by Friesz and

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

and Assabet River basins was developed by Zarriello and oth-

ers (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 water-

supply 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

ow and Lake Cochituate water levels under assumed con-

tribution or depletion of water from different surface-water

sources. Hypothetical pumping was simulated at a constant

rate of 6.65 cubic feet per second (ft

/s), with water coming

from Lake Cochituate, the Sudbury River, or combinations

of these sources to determine changes in lake level and river

ow. Simulation results indicated that reactivation of the

Birch Road wells could cause ows in the Sudbury River and

Cochituate Brook to be “substantially affected during periods

of low ow” (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 ow 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 stratication blanket a highly variable

bedrock surface and ll a deep bedrock valley that underlies

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

the regional stratigraphic framework, geotechnical sediment

classications, and geologic depositional processes respon-

sible for aquifer structure helps to appropriately represent the

aquifer system in a groundwater ow 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 ofce 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 dene elevations

and characteristics of hydrogeologic layers in the model.

Geophysical Data

Geophysical data were compiled from previous studies

and were collected during this study (g. 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 gure 2.

Passive seismic methods were used to measure depth

to bedrock at 32 sites (g. 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 anthro-

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

Zaf

, (1)

where

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

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 (g. 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 frequen-

cies (deeper depths to bedrock). Depths to bedrock at the

deeper control points in the study therefore were computed

from calibrated coefcients (a = 297.24 and b = 1.00) deter-

mined for Cape Cod (John Lane, U.S. Geological Survey,

written commun., 2011), where depths to bedrock are gen-

erally deeper and more closely matched depths than those

computed from locally derived coefcients. 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 coefcients and resonance frequencies

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

coefcients derived for Cape Cod (g. 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 (g. 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 prole 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) identied

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 (ft

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

Sudbury River

Pod Meadow

Pond

Dudley Pond

Lake Cochituate

Heard Pond

Bedrock islandBedrock 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 wellsBirch 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 (g. 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 conrm the bedrock topography in this area. Bedrock

topography likely determines groundwater ow 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 some-

what 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 weather-

ing and glacial movement. The upper bedrock in other parts of

Massachusetts has been shown to have an active groundwater

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

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

Simcox, 1994). Groundwater ow through the bedrock of the

study area was documented during construction of a water-

supply tunnel at depths of 200–500 ft below land surface.

Groundwater ow rates were weakly correlated with observed

lineaments (Mabee and others, 2002), and therefore linea-

ments were not represented in the model. In addition, ground-

water levels in the glacial sediments and streamow data do

not indicate that the hydraulic properties of the bedrock affect

shallow groundwater ow.

Glacial Sediment History

Overlying the bedrock in most of the study area are

stratied 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 (gs. 4, 5). Sedi-

ments 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 (g. 5), identied by areas

of stratied 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

ner 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

= 0.99

Resonance frequency

Estimated depth to bedrock

10 20 30

Resonance frequency, in Hertz

100

150

200

250

0 50 100 150 200 250

Estimated depth to bedrock, in feet

Depth to bedrock in boring, in feet

Local

Bedrock depths based on

horizontal to vertical

spectral ratio (HVSR) for

Cape Cod

Combined

100

150

200

250

EXPLANATION

Hydrogeology 7

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

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

Sudbury River

Dudley Pond

Lake Cochituate

Pod Meadow

Pond

Heard Pond

-100

-50

100

125

150

175

200

225

250

275

300

100

Bedrock islandBedrock island

EXPLANATION

Morphosequence 1

Undifferentiated fines

Undifferentiated coarse

Thick till

Fines at land surface

Fines at depth

Morphosequence 3

Morphosequence 2

Morpho-

sequence 1

Cross section (figure 6)

A A'

Bedrock elevation, in feet

Wayland production well

Borehole or seismic point used in

cross section and identifier

F1W-49F1W-49

Inactive

model area

Morphosequence 2

Morphosequence 3

Inactive

model area

Row 250, E-E’ (figure 11)

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

A''

A'''

Oxbow

Pod Meadow

Birch Road wellsBirch Road wells

Birch Road well

F1W-60

Seismic-DSW

WKB-51

WKW-117

WKW-27

Seismic-LA

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-21I

NS-21H

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

WKW-27

Seismic-LA

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-21I

NS-21H

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 ne-grained lacustrine deposits that are

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

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

extensive ne 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 (g. 6A). Exten-

sive ne 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 nes

consist of silty organic sediments, which may extend from the

land surface to bedrock. The ne deposits in the middle part

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

and slightly coarser (ne silts and silts) as indicated in cross-

section lines C–C′ and D–D′ (g. 6B). Contacts between sedi-

ments 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 difcult to explain. Permeable sand and gravel

deposits between the two sets of ponds should result in

groundwater ow 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 bot-

tom of the ponds. Sediment cores collected as part of an eutro-

phication 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 depos-

its 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 inows 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′,

gs. 6A and B) as evidenced by up to 12 ft of gyttja depos-

its 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 ground-

water are supplied by abundant precipitation, with average

annual precipitation in 2004–09 measuring about 50 inches

(in) at nearby Natick, Mass. (National Oceanic and Atmo-

spheric 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 northern-

most pond of Lake Cochituate to the south, Cochituate Brook

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

Pond (g. 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 north-

ernmost is shown in g. 7) connected by shallow, narrow

waterways that form a relatively long south-north trend-

ing lake. Total drainage area at the lake outlet is 17.5 mi

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 mi

surface

area and drains to the westward-owing Cochituate Brook

that connects to the northeastward-owing Sudbury River.

Streamow in Cochituate Brook (g. 7) was monitored by

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

from August 2010 through June 2012. Daily mean ow for

this entire period was 34.9 ft

/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 (g. 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 mi

and a total drain-

age area of 0.58 mi

measured at its outlet. Pod Meadow Pond

has a surface area of 0.01 mi

and a drainage area of 0.23 mi

Early topographic maps show Pod Meadow Pond as a wetland

that was likely modied by sand and gravel excavations in the

early to mid-20th century.

Outows from Pod Meadow Pond and Dudley Pond

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

provide data for groundwater model calibration. The quality

of streamow measurement data was generally considered

good to fair, but was poor (as dened by Kennedy, 1983) at

the lowest ows because of low stream depth and velocity.

Drainage areas at the outow measurement sites to Pod

Meadow Pond and Dudley Pond are 0.23 and 0.58 mi

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.

200

220

180

160

140

120

100

-20

Massachusetts Turnpike

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,0002,0000

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

Seismic-LA

Seismic-DSW

F1W-60

F1W-49

F1W-90

F1W-88

SEA-10

8-90

1-90

SEA-2

200

220

180

160

140

120

100

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

Boundaries between most contacts are uncertain

EXPLANATION

Fine to very fine sand or silt

Sand

Fines (silt and clay)

Sand and gravel

Till

Bedrock

Kettle depression

Vertical exaggeration X 20

Gyttja (organic muck)

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

200

180

160

140

120

100

-20

-40

-60

220

240

-80

200

180

160

140

120

100

-20

-40

-60

220

Dudley

Pond

Birch Road wells screened

interval 400 feet to south

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

Bedrock

Till

Sand and gravel

Sand

Bedrock

Till

Sand

Sand and gravel

Sand

Gyttja

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

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

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

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

130

140

150

160

Sudbury River

Heard Pond

Dudley Pond

Lake Cochituate

Pod Meadow

Pond

Bedrock islandBedrock 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

USGS partial-record streamgage

01098530

Dudley Pond outflow

WKW-119

F1W-84

MW-8

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 ft

/s, compared to the outow from Dudley Pond,

which ranged from 0 to 1.71 ft

/s. The relatively high and

consistent outow 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 outow is also evident from anecdotal

reports that the southwestern part of Pod Meadow Pond does

not freeze, which was conrmed during site visits in the winter

of 2010–11 (g. 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 ow upward, and the pond’s position just

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

ne-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 (g. 1).

The Sudbury River ows 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 mi

at its entrance to the study area and about

111 mi

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 ow

at Saxonville from October 1979 through September 2009 was

205 ft

/s with monthly mean ows ranging from 71 ft

/s in

September to 384 ft

/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 eleva-

tions of water bodies were used to develop a water-table map

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

water-table map indicates that regional groundwater ow 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, ground-

water ow in the shallow aquifer is towards the lake while

groundwater ow 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 low-

permeability lacustrine deposits in this area may cause locally

perched water table conditions, hydraulically separating the

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

levels from vertically paired wells in other parts of the study

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

shallow and deep parts of the aquifer, indicating that the

aquifer is unconned, 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.14

0.91

11-29-2011 0.12

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

Discharge affected by beaver dam that was impounding water; values

excluded from mean and standard deviation.

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 islandBedrock 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 gures 7 and 15; Elevation in North American Vertical Datum 1988; SS, used to calibrate steady-state model; TR, used to calibrate tran-

sient 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 les 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

eleva-

tion

(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 gures 7 and 15; Elevation in North American Vertical Datum 1988; SS, used to calibrate steady-state model; TR, used to calibrate tran-

sient 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 les 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

eleva-

tion

(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 nonpotable-

water 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

g. 7). These three wells have an average annual combined

withdrawal of 1.40 ft

/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 oper-

ated wells at the Birch Road site for municipal water supply

(g. 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 ft

/s (3.5 Mgal/d) but were likely

pumped at a rate closer to 4.9 ft

/s (SEA Consultants, Inc.,

1992). The wells were operated intermittently during the mid-

to late 1970s because of high iron and manganese concentra-

tions, 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, g. 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 rst large drinking-water

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

through 1951, when it was nally 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 numer-

ous groundwater and surface-water withdrawals that affect

streamow. Zarriello and others (2010) reported 6 production

wells upstream of the study area with total average annual

withdrawals of about 4.5 ft

/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 with-

drawal of about 5.3 ft

/s. The effects of these withdrawals are

included implicitly in the groundwater model through speci-

ed stream inows and are particularly evident during periods

of low ow.

The operation of surface-water reservoirs upstream of

the study area also affects streamow and can further decrease

low-ows 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 ow, can

cause large percentage changes in streamow 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 ow model of the aquifer based

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

developed to represent hydrogeologic conditions and simulate

groundwater ow and interaction with surface waters in the

study area. MODFLOW-NWT is a nite-difference ground-

water 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 under-

standing 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 surface-

water resources, and to identify further needs for data collec-

tion and model renements.

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 con-

structed 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 unconned 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 simula-

tions, can become wet again.

The model was initially calibrated to long-term steady-

state conditions that represent average annual conditions in

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

the model was modied to simulate transient, or unsteady,

conditions. Under steady-state simulations, stresses such as

recharge from precipitation, groundwater withdrawals, and

ow 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 (g. 9) of 9.0 mi

was spatially dis-

cretized 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 topog-

raphy near pumping wells could be accurately simulated. The

model boundary positions were chosen to correspond with

natural no-ow boundaries (no appreciable lateral ows into

or out of the model area) or to be sufciently far from the

main areas of interest so that the boundary would have little

effect on simulation results.

Vertically, the model was discretized into ve layers of

variable thickness that represent different aquifer sediments

as shown in gure 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 ne 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 ne

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 ow.

Layers 2 and 4 dene the principal aquifer in the study

area and in some areas are hydraulically separated by the ne

deposits represented in layer 3. The ne-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 ne 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 bed-

rock is close to the surface or where Lake Cochituate is pres-

ent. 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 nes, layers 2 and 4 were assigned a thickness equal to

one-half of the remainder of the thickness of the glacial strati-

ed deposits.

Groundwater Flow Model 19

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

Sudbury River

Pod Meadow

Pond

Dudley Pond

Lake Cochituate

Bedrock island

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

0.2

107.5

11.6

1.4

MV-1

Bedrock island

Grid cells 50 feet square

Pod Meadow Pond

outflow

Bedrock island

Grid cells 50 feet square

Simulated

Birch Road wells

Simulated

Birch Road wells

Pod Meadow Pond

Pod Meadow Pond outflow

Pod Meadow Pond

outflow

01098500

01098530

Dudley Pond outflow

Oxbow

Pod Meadow

100

120

140

160

180

200

220

240

260

280

100

120

140

160

180

200

220

240

260

280

300

320

340

360

Row 20

Column 20

Row 20

Column numbers

Model grid: 360 rows and 280 columns

Column 20

Row numbers

Row 250, E-E’ (figure 11)

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

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.

0.2

50,000

Simulated

streams

0.01

140

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

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

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 uxes 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 fth year. Simulation results

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

the rst 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-ow bound-

ary (g. 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

ows were assumed to be negligible through sediments under-

lying 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 (g. 4).

At the northeastern boundary, surface water gradients are low

over a wide wetland area, suggesting that groundwater gradi-

ents, and hence ow rates, are also low across the boundary.

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

that the likely limited groundwater ows through them would

have minimal effect on the model simulations.

Inows to Lake Cochituate were assumed to be fully

represented by the assigned stream inow at the southern end

of the lake. Lateral groundwater ow in the sediments beneath

Lake Cochituate across the southern model boundary is not

likely to be signicant because north-south surface-water gra-

dients are low in this area. However, there were few ground-

water level observations available to make a full assessment

of groundwater ow under the lake. Groundwater ow across

the northern boundary is likely small as indicated by the low

groundwater table gradient (g. 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 (g. 7), the Sudbury River acts as a hydraulic bound-

ary having no groundwater ow 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 no-

ow boundary.

The only exceptions to no-ow boundaries around the

perimeter of the active model area are inows and outows

to streams represented by the streamow-routing package

(SFR2; Niswonger and Prudic, 2005). Flows were specied

at the rst 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. Outow 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. Outow

by evapotranspiration was not explicitly simulated but was

incorporated into the net recharge values used in the model.

Hydraulic Conductivity

The rate of groundwater ow is determined by hydraulic

conductivity and water-level gradients. Hydraulic conductiv-

ity of aquifer material was assigned on the basis of sediment

lithology (table 4; gs. 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 conductiv-

ity varied somewhat but was generally uniform within a layer

(g. 10). The ratio of horizontal to vertical hydraulic conduc-

tivity (anisotropy) varied by layer from 5 to 20 (table 4). Verti-

cal 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 ne sand with silt; hence, these deposits do not

represent an aquitard (impermeable boundary to groundwater

ow) but rather a more restrictive layer to groundwater ow

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

considered slightly permeable for groundwater ow and were

assigned a hydraulic conductivity equal to that of ne depos-

its. 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 pres-

ent 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 coefcients, specic storage (Ss) and spe-

cic yield (Sy), dene the volume of water going in or out

of aquifer storage as groundwater levels rise and fall under

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

units of inverse feet (ft

-1

) is used to dene the volume of water

released per unit volume of aquifer per unit of groundwater

level change. When a model cell is partially saturated, Sy (per-

cent) is used and is dened as the volume of water that drains

by gravity per unit volume of aquifer. MODFLOW-NWT

automatically determines which storage coefcient to use from

the simulated water level and top and bottom elevations of

each cell. Uniform values were used to dene 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 modied during model calibration.

Streams

Streams are represented in layer 1 with the Streamow-

Routing Package (SFR2) by Niswonger and Prudic (2005).

SFR2 simulates ow between streams and the aquifer as

a head-dependent ux and tracks streamow as it moves

downstream. Modeled streamows represent baseow, the

groundwater component of streamow, and do not include

surface runoff or interow components of streamow. Simu-

lated streamow in the model is the combination of speci-

ed upstream inows 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 inltration into

the aquifer exceeds inow from upstream reaches.

Streams in the model area were dened by 21 segments

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

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

outows 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 (ne 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

Specic yield (fractional percent)

0.15 0.15 0.15 0.15 0.15

Specic 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.

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

200

100

-100

Vertical exaggeration X 10

Elevation in NAVD 88

Cochituate Brook

Sudbury River

Pod Meadow

Distance along model grid, in feet

Elevation, in feet

Fines

West East

Vertical exaggeration X 10

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

100

300

200

-200

-100

Sudbury River

Distance along model grid, in feet

Elevation, in feet

Lake Cochituate

Fines

Sand and gravel

50,000

140

Fines

B. Modeled layers and hydraulic conductivities, row 250

EXPLANATION

Hydraulic conductivity, in feet per day

Inactive model area

Inactive

Bedrock and till

Layer 3

A'''A''

E E'

Values as labeled

Fines

High conductivity zone for lake

140

Fines

Sand and gravel

Layer 3

Layer 2

Layer 4

Layer 5

0.2

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 (g. 9). Segments 5, 6, 8, 16, 18, and 21 dene the

Sudbury River, segments 1, 2, and 4 dene Lake Cochituate

and its outlet, segments 9–15 dene Pod Meadow Pond and

its outlet and wetland areas, and segments 19 and 20 dene

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

coefcient of 0.03.

Stream cells were assigned across the high conductiv-

ity zones representing Dudley Pond, Pod Meadow Pond, and

Lake Cochituate (g. 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.

Inows were specied to streams entering the active

model area on the basis of daily ow records from the

Saxonville streamgage (01098530) from January 1980 through

December 2010. For the steady-state model, inows of 107.5,

11.6, and 1.4, ft

/s were assigned to the Sudbury River, Lake

Cochituate, and the tributary to Cochituate Brook, respectively

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

ow observed at Saxonville (148 ft

/s) and were proportioned

by the drainage area to each boundary inow point. A small

inow of 0.2 ft

/s was assigned to Dudley Pond based on the

small drainage area to the pond outside of the active model

area. Streamow exits the model through the Sudbury River.

In the transient model simulations, stream inows were

the same for all scenarios and were assigned to represent

monthly low-ow conditions. The 75-percent daily ow

duration (ow value exceeded 75 percent of the time) was

determined for each month of the year from recorded daily

ows at the Saxonville streamgage from January 1980 through

December 2010. The 75-percent daily ow durations by

month at Saxonville were proportioned by drainage area,

as previously described, to specify inows at the model

boundary. No inow was assigned to Dudley Pond for

transient simulations. The 75-percent daily ow duration was

chosen to represent conservative monthly low-ow conditions

throughout the year in addition to the normal seasonal

variations in ows. For example, the Sudbury River at

Saxonville (01098530) streamgage had a 75-percent daily ow

duration for September of 18.0 ft

/s, whereas the 90-percent

daily ow duration for the entire period of record, a commonly

used low-ow stress indicator threshold, is 44 percent greater

(26.0 ft

/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 inow to

the aquifer and was simulated by using the MODFLOW-

NWT 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 oth-

ers, 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 inow to the

lake and the stream traversing the lake area.

For the transient model, recharge rates were varied

monthly to reect seasonal precipitation and evapotranspi-

ration variability over the year. Monthly recharge rates, as

percentages of annual recharge (table 6), were determined

from soil moisture budgets calculated with WATBUG soft-

ware (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 ground-

water 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. Long-

term 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

streamow. Adjustments of hydrologic parameter values and

boundary conditions were made within a range of values

determined from previous groundwater models of glacial-

sediment 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 t.

Steady-State Calibration

Steady-state calibration focused on streambed hydrau-

lic conductivity and vertical and horizontal aquifer hydrau-

lic conductivity to reduce simulation errors in pond levels,

groundwater levels, and streamows. 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 land-

surface 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 streamows 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 (g. 12),

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

added pumping stresses and decreased recharge cause addi-

tional 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 gure 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 Specied

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 inuences the

timing and magnitude of ows between streams and the aqui-

fer. Reducing streambed conductivity decreases ow 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 outows (segments 9 and 19, respec-

tively; g. 9), the average measured ows (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 xed at 5.0 for

layers 1 and 5 and calibrated for layers 2–4 (table 4). Verti-

cal 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

outows 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 outows. 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 streamow was evaluated for the stream

segments draining Dudley Pond and Pod Meadow Pond by

comparing simulated steady-state baseow in each stream

to observed average discharge measurements made from

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

74 percent higher than observed outow 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 outows 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 ux 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 ux -0.02 -0.05 -3.26

Change caused by pumping -1.22 -0.01 -1.58

Pond (1.22 and 0.7 ft

/s, respectively) and 16 percent lower

than observed outow from Dudley Pond (0.67 and 0.80 ft

/s,

respectively). As discussed previously, Pod Meadow Pond

outows are thought to be a result of focused groundwater

discharge and, therefore, were a priority in the model calibra-

tion, although these ows were difcult to simulate accurately.

Oversimulation of Pod Meadow Pond outows was likely

caused by poor model representation of local aquifer hydrau-

lic conductivity patterns around the Birch Road wells and

the pond. Outows from Dudley Pond are largely affected

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

outow stream going dry in August 2011 (table 1). Further-

more, because simulated outows from Dudley Pond reect

the simulation of baseow but not runoff and interow, the

undersimulation of measured Dudley Pond outows was not

considered unreasonable. Overall, the accuracy of simulated

baseow was considered acceptable for the purposes of this

study given the difference between baseow and total stream-

ow and the short period of available observations.

Steady-state simulated groundwater levels are generally

in agreement with the estimated groundwater table (g. 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 undersimu-

lated errors (g. 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 (g. 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 conductiv-

ity variations in the Birch Road well area. Local variations

in hydraulic conductivity could be examined in future model

modications 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

eld 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 (g. 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

reect 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 renements

to spatial distribution of hydraulic conductivity could improve

the model t.

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

Dry

Wayland production well

Active model cells:

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

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

Simulated water level, in feet

Observed water level, in feet

Transient Model Calibration

Transient model calibration used data for seasonal

groundwater levels and pumping-inuenced groundwater

levels collected during the 2006 aquifer test, but lake and

pond elevations were not used for the transient model calibra-

tion because they are strongly affected by surface runoff and

evapotranspiration, which were not explicitly simulated. Stor-

age coefcients 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 hydro-

geologic 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 (g. 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 dene simula-

tion 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 simula-

tion 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 con-

ditions; thereafter, a daily stress period was used (41 through

60) to represent the aquifer test. Model boundary conditions

(stream inows, recharge, pumping stresses) were assigned

to reect 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 (g. 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.

Pod Meadow

Pond

Dudley Pond

Heard Pond

Lake Cochituate

Sudbury River

Bedrock island

Bedrock

island

EXPLANATION

Simulated ground-

water 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

150

120

160

110

130

120

Groundwater Flow Model 31

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-123

WKW-124

F1W-84

F1W-91

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

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 gure 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 (g. 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 consid-

ered acceptable because simulated levels are monthly averages

and are expected to have less variation than observed instanta-

neous levels.

Simulation of the 2006 aquifer test produced groundwater

drawdowns that match observations relatively well—water

levels were oversimulated at four observation wells and under-

simulated at three observation wells (table 10 and g. 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 hav-

ing 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 ne and coarse sediment

deposits in this area. The simulated rate of drawdown also

reasonably matched the observed (g. 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

coefcients 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 (g. 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 simplica-

tion 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 conclu-

sions 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 uncer-

tainty 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 var-

ied hydraulic model parameters over likely value ranges

(table 11). Model parameters were modied 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 pump-

ing rate of 4.9 ft

/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 ux, from

the high-conductivity cells representing the lake to the sur-

rounding 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

) was dened as the number of days for 50 percent of

streamow depletion to occur following 1 month of pumping

from the Birch Road wells at a hypothetical rate of 4.9 ft

(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

ows continued to be simulated through December of year

5 (month 60). Streamow depletion in the Sudbury River

was calculated at the exit of the model by subtracting each

scenario’s ow from the corresponding ow 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

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

was interpolated to occur halfway through February and to

equal 14 days. Although some slight variation in t

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 rate

(cubic feet per second)

Streamflow response time

(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

Specic yield 1.58 -- -- 8.3 0 35.4

Specic storage 1.58 -- -- 8.3 16.6 7.8

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

Response time dened as time to reach 50 percent of streamow 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 conductance

(feet per day)

70 35 105 n33m2 n40 n50 tn51m tn52 tn62

Average aquifer hydraulic conductivity

(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

Specic yield (percent) 0.15 0.10 0.3 n33m2 -- -- tn51m tn58a tn68

Specic storage (ft

-1

) 10

-5

-6

-4

n33m2 -- -- tn51m tn59 tn69

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 inow from its

17.5 mi

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 rened and realistic simulations of streamows

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 ow of 0.5 ft

/s would

be added to simulated baseow. Additionally, in the transient

simulations no inow 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 ows do not change the main ndings of

this study (sources of water to the wells, timing of streamow

response to pumping, and the potential for pumping schedules

to ameliorate low ows). If the groundwater model is

combined with the surface-water model of Zarriello and

others (2010) to simulate full streamow, assignment of these

ows 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 ow across

the northern and southern boundaries of the model is a

potential limitation of the model, but few groundwater level

measurements are available to conrm the assumption. The

balance of evidence indicates there are likely low rates of

groundwater ow 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-ow conditions

at the boundaries, then pumping stress could potentially

draw water through these boundaries and the model would

oversimulate the effect of pumping on streamow. Therefore,

the no-ow boundaries are conservative, in terms of

simulating the maximum induced streamow.

Assignment of model parameter values, such as hydraulic

conductivity, Ss, and Sy, could be improved by using auto-

mated 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 ft

/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 ne 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 ne 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 coefcient for partially saturated aquifer

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

to 30 percent caused t

to increase from 8.3 to 35.4 days.

Decreasing Sy from 15 to 10 percent caused t

to decrease

from 8.3 to 0 days, meaning that t

was reached within

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

calculations contained streamow mass balance errors that

were large enough to affect the calculation of t

for run tn58a,

so the exact value of t

is uncertain. These mass balance

errors were caused by numerical artifacts that, while not large

relative to most ows in the model, were large relative to

the small differences in streamow used to calculate t

. The

timing of surface-water response is also sensitive to hydraulic

conductivity of sand and gravel (K

). Increasing average

model K

from 113 to 226 ft/d causes t

to decrease from 8.3

to 0 days. The effect of decreasing K

could not be assessed

because the model becomes numerically unstable unless K

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

However, decreasing K

would likely cause t

to increase

appreciably. The sensitivity of surface-water response time

to pumping reinforces the nding 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 congured, 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

ner 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. Rened

hydraulic conductivity assignment around the well site would

also likely reduce the oversimulation of baseow 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 sensitiv-

ity. The model is also sensitive to Sy, which adds uncertainty

to ndings 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 rst 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 baseow. 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 (scenar-

ios 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 reduc-

tion in total annual recharge of only 7.1 percent, but it occurs

in the 3-month window before the critical low-ow month

of September. The dry conditions, particularly in combina-

tion with the low stream inow values that are specied for

all transient simulations, create very low-ow ecologically

conservative conditions.

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

stream system, average annual streamow in the Sudbury

River at the Saxonville streamgage (01098530) is 205 ft

/s.

Hence, on average, the hypothetical 4.9 ft

/s withdrawal from

the Birch Road wells accounts for only a little over 2 percent

of average daily Sudbury River streamow. Concern arises

during low-ow conditions when water use and regulation rep-

resent an increasingly large proportion of the total streamow.

The potential effects of the additional stress on Sudbury River

low ows from Birch Road well withdrawals were addressed

by transient model simulations.

Pumping induces additional recharge to the aquifer sys-

tem from Lake Cochituate. Imposing a hypothetical constant

pumping rate of 4.9 ft

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

(scenario 2) induces additional recharge from Lake Cochituate

at a rate of 1.6 ft

/s, which is about equal to 32 percent of the

hypothetical pumping rate (table 9). Under constant pump-

ing, Dudley Pond receives slightly less groundwater discharge

from the aquifer, but more signicantly, groundwater dis-

charge 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 inows must equal outows; in

a transient model inows must equal outows plus or minus

changes in storage. A natural consequence of this requirement,

without a change to the boundary conditions, is that ground-

water 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

baseow reduction, whereas streams close to and downgradi-

ent from the pumping wells experience greater reduction. As

a percentage of the imposed pumping (4.9 ft

/s), the reduction

of baseow in streams (scenario 2) is illustrated by the width

of the stream lines shown in gure 18. Streamow depletions

in the Sudbury River, equal or nearly equal to the constant

4.9 ft

/s pumping rate, occur downstream of the oxbow. The

stream traversing Pod Meadow Pond (segment 9) experiences

the greatest percentage of streamow reduction (100 percent)

as streamow goes to zero under the steady-state imposed

pumping rate of 4.9 ft

/s; this reduction is not visible on

gure 18 because it falls below the display threshold value of

1.5 ft

/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-ow recharge

conditions. All the transient scenarios (scenarios 3–13) take a

conservative approach to streamow simulation by specify-

ing inow to stream segments at the model boundary as the

75-percent daily ow 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 inows and outows 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 inows to and outows from the model are 11.5 ft

/s;

83 percent of the inow 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. Outow 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 inows to and outows from

the model are 13.1 ft

/s; inow from recharge by precipitation

is the same as before, but decreases as a percentage of the total

inow (73 percent), whereas recharge from streams increases

(27 percent). Increased recharge from streams is mostly from

induced inltration of stream water from stream cells travers-

ing Lake Cochituate (table 5, g. 9). Most outow from the

model is still discharge to streams (52 percent) but decreases

from 10.1 to 6.8 ft

/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

outow from the model under average steady-state condi-

tions, 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.

Sudbury River

Sudbury River Sudbury River

Dudley

Pond

Lake Cochituate

Heard Pond

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

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 (g. 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 outow from the aquifer, although

in September pumping from Wayland wells accounts for

18 percent of total outow 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 (g. 19).

Evaluation of Pumping Strategies to Reduce

Low-Flow Stresses

The transient model was used to determine the tim-

ing of surface-water response to pumping and to evaluate

pumping strategies to reduce the effects of pumping during

low-ow 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 ft

/s for 1 month, and subsequent streamow depletions

are determined (g. 20). Associated streamow 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 baseow, or in other

-25

-20

-15

-10

-5

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 streamow depletion, decreases by about 40 percent.

Within 2 months, the reduction of baseow decreases by about

80 percent, and within 4 months, by about 90 percent (g. 20).

The fast response of surface water to pumping stresses

indicates the potential for pumping to be managed to reduce

streamow depletions during periods of low ow.

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 ows (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 ft

/s. Scenarios 6–8 and 11–13 simulated a reduced

pumping rate of 0.1 ft

/s for 1, 3, and 6 months of the year

and a pumping rate of 4.9 ft

/s in the remaining months. A

minimum pumping rate of 0.1 ft

/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 streamows. In actual operation, pumping

rates could be determined from streamow, 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

ows in the Sudbury River at the model exit under average

monthly climatic conditions and various pumping scenarios

100

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, g. 21). The greatest reduction

in streamow 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 baseow nearly the same as without

pumping. Simulated baseow under dry recharge conditions

is only slightly less than simulated baseow under average

recharge conditions.

Streamow depletion increases going downstream in

the Sudbury River as the effects of pumping accumulate

(g. 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

streamow depletion nearly reaches its full extent.

Reducing the pumping rate in September and in prior

months decreases streamow depletion compared to con-

stant pumping (g. 23). Under average recharge conditions,

reducing pumping for 1 month lowers September streamow

depletion from 4.4 to 3.1 ft

/s (scenarios 5 and 6), reducing

pumping for three months (scenario 7) lowers depletion to

0.8 ft

/s, and reducing pumping for 6 months (scenario 8)

drops depletion to 0.3 ft

/s. In terms of percentages, lowering

maximum pumping rates to near 0 during September reduces

September streamow depletion by 29 percent, and lowering

pumping rates for 3 months (July through September) reduces

streamow 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 streamow depletion from

4.2 to 3.0, 1.1, and 0.5 ft

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

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

Jan.

Feb. March April May June July Aug. Sept. Oct. Nov. Dec.

Constant pumping

(4.9 ft

/s)

Average Dry

Recharge condition

Period of reduced pumping (0.1 ft

/s)

Pumping at 4.9 ft

/s for periods of

12 months

(tn81)

(tn82)

(tn83)

(tn84)

(tn71)

(tn72)

(tn73)

(tn74)

Scenario name in parentheses;

/s, cubic feet per second

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 1 2 3

No pumping (scenario 4)

Pumping at 4.9 cubic feet per

second for periods of

6 months (scenario 8)

9 months (scenario 7)

11 months (scenario 6)

12 months (scenario 5)

Simulated streamflow, in cubic feet per second

River mile downstream from model 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 streamow than it does during average

recharge conditions (g. 23). Relative to streamow under no

pumping of the Birch Road wells, constant pumping of these

wells at a combined rate of 4.9 ft

/s decreases streamow

by 18.8 and 19.2 percent under average conditions and dry

recharge conditions, respectively (blue lines in g. 23).

When pumping is reduced to near zero during September

only, streamow depletion is 13.3 and 13.7 percent of total

streamow under average conditions and under dry conditions,

respectively (purple lines in g. 23). The streamow

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

streamow is based on assigned inow values of 75 percent

daily ow duration for each month, representing low-ow

conditions. Under average monthly ow conditions, the

percent streamow reductions would be substantially less.

Similar patterns of reduced streamow 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 con-

trols, particularly during periods of low ow when the stream

ecology is the most vulnerable to water-management opera-

tions. During summer months and periods of low precipita-

tion, the river is primarily fed by groundwater. Because of

the river’s close hydraulic connection with aquifers, ground-

water withdrawals can reduce streamow during low-ow

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 with-

drawals may further reduce streamows 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-ow periods. Goals of the

investigation included improving understanding of hydrology

in the study area, determining rates of ow 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-ow 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 uctuating

water table. The model is calibrated with geologic and hydro-

logic data compiled from prior studies and new data collected

during this study. Simulated groundwater levels and stream-

ows 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 ft

/s,

which is equal to about 32 percent of the simulated

4.9 ft

/s pumping rate.

• Groundwater withdrawals reduce ow in the Sudbury

River and tributary streams. The Sudbury River down-

stream 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, streamow 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

streamow depletions during periods of low ow by

altering pumping rates appropriately for seasonal or

anticipated ow conditions. Streamow depletion

during September, typically the month of lowest ow,

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 streamow 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 specic 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 Scientic 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 inows and out-

ows, 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 Milford-

Dedham 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 water-

supply 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/leviewer/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/leviewer/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/leviewer/

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/leviewer/

Rtn.aspx?rtn=3-0000629.

Boutt, D.F., Diggins, P., and Mabee, S.B., 2010, A eld study

(Massachusetts, USA) of the factors controlling the depth

of groundwater ow 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 north-

ern part of glacial Lake Charles, Massachusetts: Journal of

Geology, v. 12, no. 3, p. 198–214.

DeSimone, L.A., 2004, Simulation of ground-water ow and

evaluation of water-management alternatives in the Assabet

River Basin, eastern Massachusetts: U.S. Geological Survey

Scientic 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 ow and evalua-

tion 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 Ofce for Administration and

Finance: accessed November 22, 2010, at

http://www.mass.gov/mgis/p.htm.

Massachusetts Department of Conservation and Recre-

ation (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 ow in the

Plymouth-Carver-Kingston-Duxbury aquifer system, south-

eastern Massachusetts: U.S. Geological Survey Scientic

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 ow model and the MOC3D

solute-transport model: U.S. Geological Survey Water-

Resources 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, Surcial 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, Surcial 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, Massachu-

setts: 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, MOD-

FLOW-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

ow 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 ow process: U.S. Geologi-

cal 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 observa-

tion, 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, Micro-

tremor measurements used to map thickness of soft sedi-

ments: 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 investiga-

tions of the United States Geological Survey: chap. Al3,

book 3, Computation of continuous records of streamow,

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 Environmen-

tal 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, Cor-

relation of lineaments to groundwater inows 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 Streamow-Routing (SFR2) package to include unsatu-

rated ow beneath streams—A modication 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 nal report, Birch Road

well re-activation, Framingham, Massachusetts: SEA Con-

sultants, 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 investiga-

tion, tier classication, and phase II scope of work—Mass.

DEP RTN 3-27985: Sovereign Consulting, Inc., Manseld,

Mass., accessed May 2011, at http://public.dep.state.ma.us/

leviewer/Default.aspx?formdataid=0&documentid=51762.

Stone, J.R., and Stone, B.D., 2006, Surcial geologic map of

the Clinton-Concord-Grafton-Medeld 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/WellheadProtectionPlan6MB-

2June2011.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 streamow

and aquatic habitat in the Sudbury and Assabet river

basins, Massachusetts: U.S. Geological Survey Scientic

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

[email protected]

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