Cover photo: Charles River, Medway, Massachusetts

Courtesy of: John R. Eggleston
U.S. Department of the Interior

U.S. Geological Survey

Simulation of Ground-Water Flow and 

Evaluation of Water-Management 

Alternatives in the Upper Charles 

River Basin, Eastern Massachusetts

By LESLIE A. DESIMONE, DONALD A. WALTER, JOHN R. EGGLESTON, and MARK T. NIMIROSKI 

Water-Resources Investigations Report 02-4234

Prepared in cooperation with the

MASSACHUSETTS DEPARTMENT OF ENVIRONMENTAL MANAGEMENT and the

MASSACHUSETTS DEPARTMENT OF ENVIRONMENTAL PROTECTION

Northborough, Massachusetts

2002 
U.S. DEPARTMENT OF THE INTERIOR

GALE A. NORTON, Secretary

U.S. GEOLOGICAL SURVEY

Charles G. Groat, Director

The use of trade or product names in this report is for identification 

purposes only and does not constitute endorsement by the 

U.S. Government.

For additional information write to: Copies of this report can be purchased from:

Chief, Massachusetts-Rhode Island District U.S. Geological Survey

U.S. Geological Survey Branch of Information Services

10 Bearfoot Road Box 25286

Northborough, MA 01532 Denver, CO 80225-0286

or visit our Web site at

http://ma.water.usgs.gov
Contents III

CONTENTS

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

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

Purpose and Scope ...................................................................................................................................................... 5

Description of the Study Area..................................................................................................................................... 5

Previous Studies .......................................................................................................................................................... 6

Acknowledgments....................................................................................................................................................... 6

Upper Charles Stream-Aquifer System................................................................................................................................. 7

Hydrogeology.............................................................................................................................................................. 7

Water Withdrawals and Return Flow .......................................................................................................................... 9

Streamflow and Water Levels...................................................................................................................................... 18

Water Balance.............................................................................................................................................................. 28

Simulation of Ground-Water Flow........................................................................................................................................ 31

Steady-State Numerical Models.................................................................................................................................. 31

Spatial Discretization ........................................................................................................................................ 33

Boundary Conditions......................................................................................................................................... 33

Stresses.............................................................................................................................................................. 35

Recharge.................................................................................................................................................. 35

Water Withdrawals and Discharges......................................................................................................... 40

Hydraulic Properties.......................................................................................................................................... 40

Model Calibration ............................................................................................................................................. 42

Recharge.................................................................................................................................................. 43

Hydraulic Properties................................................................................................................................ 44

Calibration Results .................................................................................................................................. 44

Transient Numerical Models ....................................................................................................................................... 52

Temporal Discretization and Initial Conditions ................................................................................................ 52

Boundary Conditions and Stresses.................................................................................................................... 52

Hydraulic Properties.......................................................................................................................................... 54

Model Calibration ............................................................................................................................................. 54

Model Limitations ....................................................................................................................................................... 60

Evaluation of Ground-Water-Management Alternatives....................................................................................................... 61

Simulation of Increased Water Withdrawals and Altered Recharge ........................................................................... 62

Simulation of Increased Withdrawals ............................................................................................................... 62

Simulation of Altered Recharge ....................................................................................................................... 73

Simulation-Optimization of Water Withdrawals and Stream Base Flow in the Populatic Pond Area........................ 75

Theory and Methods of the Simulation-Optimization Approach...................................................................... 77

Application in the Populatic Pond Area............................................................................................................ 78

Hydrologic Stresses and Responses ........................................................................................................ 78

Management Model Applications ........................................................................................................... 82

Increase Withdrawals while Maintaining Stream Base Flow above Specified

Minimum-Flow Requirements ........................................................................................ 82

Increase Withdrawals while Maintaining Stream Base Flow at Current Levels ........................... 85

Increase Stream Base Flow while Maintaining Existing and Proposed Withdrawals................... 85

Limitations of Management-Model Applications......................................................................... 86

Summary ............................................................................................................................................................................... 87

References ............................................................................................................................................................................. 90
IV Contents

FIGURES 

1–3. Maps showing: 

1. Location and surficial geology of the study area, upper Charles River Basin, eastern Massachusetts......... 3

2. Thickness of the stratified glacial deposits in the upper Charles River Basin .............................................. 8

3. Supply wells, surface-water sources, and wastewater discharges in the upper Charles Basin...................... 13

4. Diagram showing summary of water withdrawals, discharges, and transfers in the upper Charles River 

Basin, 1989–98 ...................................................................................................................................................... 14

5. Graph showing average annual and monthly municipal withdrawals and wastewater discharges in the 

upper Charles River Basin, 1989–98: (A) Average annual municipal withdrawals and wastewater 

discharges, and (B) Average monthly municipal withdrawals and wastewater discharges ................................... 15

6, 7. Maps showing: 

6. Water and sewer systems in the upper Charles River Basin.......................................................................... 16

7. Observation wells, pond-measurement sites, and streamflow-measurement sites in the upper

Charles River Basin....................................................................................................................................... 19

8. Streamflow in the Charles River at the U.S. Geological Survey gaging station at Medway (01103280) 

for the period November 1997 through December 2000 .............................................................................. 24

9, 10. Graphs showing: 

9. Ground-water levels measured in selected observation wells in the upper Charles River Basin, 

eastern 1999–2000......................................................................................................................................... 25

10. Daily precipitation and measured streamflow during base-flow conditions, 1999–2000, and 

estimated monthly average base flow at selected streamflow-measurement sites in the upper 

Charles River Basin: (A) Charles River at North Bellingham (01103206), (B) Hopping Brook 

near West Medway (01103217), (C) Charles River near Millis (01103305), and (D) Bogastow 

Brook below Great Black Swamp near Millis (01103393)........................................................................... 28

11. Schematic diagram showing hydrogeologic units and flow components of the water balance for the 

upper Charles River Basin ..................................................................................................................................... 30

12. Map showing west and east model areas, grid, and boundary of the active model areas for the simulation 

models of the upper Charles River Basin .............................................................................................................. 32

13. Schematic section showing active model cells and model layers for the simulation models of the upper 

Charles River Basin ............................................................................................................................................... 34

14–16. Maps showing: 

14. Upland recharge subbasins and simulated surface-water features for the simulation models of the 

upper Charles River Basin, eastern Massachusetts: (A) West model area, and (B) East model area ............ 37

15. Areal recharge rate from precipitation recharge and septic-system return flow for the simulation 

models of the upper Charles River Basin...................................................................................................... 39

16. Hydraulic conductivity zones for the simulation models of upper Charles River Basin: (A) West 

model area, and (B) East model area............................................................................................................. 45

17. Graph showing relation between observed and model-calculated water levels and stream base flow for 

average conditions, 1989–98, for steady-state simulation models of the upper Charles River Basin:

(A) Water levels, and (B) Stream base flow........................................................................................................... 48

18. Map showing model-calculated steady-state water table in the upper Charles River Basin ................................. 49

19–31. Graphs showing: 

19. Monthly average recharge rates for transient simulation models and comparison of monthly average 

rates for flows to the active model area from upland till and bedrock areas with estimated monthly 

base flow in two streams draining upland watersheds in the upper Charles River Basin: (A) Monthly 

average recharge rates for stratified glacial deposits and rates for flows to the active model area from 

upland till and bedrock areas, and (B) Monthly average rates for flow to the active model area from 

upland till and bedrock areas and estimated average monthly base flow in the Charles River above 

Cedar Swamp Pond near Milford (01103110) and in Dirty Meadow Brook at Hollis Street near 

Holliston (011033885) .................................................................................................................................. 54

20. Model-calculated and observed water-level fluctuations during the average annual cycle for 

selected observation wells and ponds in the upper Charles River Basin: (A) Observation wells,

and (B) Ponds ................................................................................................................................................ 56
Contents V

21. Model-calculated and observed average monthly base flow in the Charles River, Bogastow Brook, 

and tributaries, upper Charles River Basin: (A) Charles River at Factory Pond at South Milford,

(B) Charles River at North Bellingham, (C) Charles River at West Medway, (D) Charles River near 

Millis, (E) Bogastow Brook below Great Black Swamp near Millis, (F) Hopping Brook near West 

Medway, (G) Mine Brook near Franklin, (H) Chicken Brook below Milk Pond near West Medway,

(I) Mill River below Bush Pond near City Mills, and (J) Dopping Brook at Whitney Street, Holliston...... 58

22. Relation between observed and model-calculated average March and September base flow for 

measurement sites in the upper Charles River Basin: (A) March, and (B) September.................................. 60

23. Permitted and registered withdrawal rates under the Massachusetts Water Management Act and 

1989–98 average withdrawal rates for municipal and large non-municipal water systems in the 

upper Charles River Basin............................................................................................................................. 62

24. Average monthly withdrawals and discharges by town or treatment facility for 1989–98 and two 

hypothetical scenarios of increased withdrawals in the upper Charles River Basin ..................................... 64

25. Model-calculated monthly average base flow for March and September along the Charles River 

for 1989–98 average conditions and two scenarios of increased withdrawals in the upper Charles 

River Basin: (A) March, and (B) September ................................................................................................. 67

26. Change in model-calculated base flow from 1989–98 flows at selected sites along the Charles River 

for two scenarios of increased withdrawals in the upper Charles River Basin: (A) Charles River 

above the Milford Treatment Facility (MTF) at the town line between Milford and Hopedale,

(B) Charles River below the MTF at station 01103140, (C) Charles River between Beaver and Stall 

Brooks at station 011032056, (D) Charles River between Chicken and Shepards Brooks at

station 01103260, (E) Charles River between the Charles River Pollution Control District Treatment 

Facility and Mill River, and (F) Charles River below Mill River at station 01103305................................. 68

27. Model-calculated monthly average base flow for March and September along Mine Brook for 

1989–98 average conditions and two scenarios of increased withdrawals in the upper Charles River 

Basin: (A) March, and (B) September ........................................................................................................... 69

28. Model-calculated monthly average base flow for March and September along the Mill River for 

1989–98 average conditions and two scenarios of increased withdrawals in the upper Charles River 

Basin: (A) March, and (B) September ........................................................................................................... 71

29. Model-calculated average monthly water levels in Kingsbury Pond for 1989–98 conditions and in 

two hypothetical scenarios of increased water withdrawals in the upper Charles River Basin .................... 72

30. Model-calculated monthly average base flow for March and September along Mine Brook for 

1989–98 average conditions and three scenarios of altered recharge in the upper Charles River 

Basin: (A) March, and (B) September ........................................................................................................... 74

31. Relation of increases in model-calculated monthly average base flow in Mine Brook and added 

recharge in the Mine Brook Subbasin, upper Charles River Basin............................................................... 75

32. Map showing aquifer area near Populatic Pond where optimization methods were applied in the upper 

Charles River Basin ............................................................................................................................................... 76

33–35. Graphs showing:

33. Relation of the time required for streamflow depletion to occur in response to pumping and the 

distance of the supply well from the nearest simulated stream for six supply wells in the

Populatic Pond area, upper Charles River Basin........................................................................................... 82

34. Monthly pumping rates for 1989–98 and for management-model applications for increased withdrawals 

or increased streamflow in the Charles River during the low-flow period in the Populatic Pond area of 

the upper Charles River Basin: (A) 1989–98, (B) Increased withdrawals with a minimum-streamflow 

requirement, (C) Increased withdrawals with no additional streamflow depletion, (D) Increased 

streamflow with existing and proposed withdrawals for two towns separately, and (E) Increased 

streamflow with existing and proposed withdrawals for two towns combined............................................. 84

35. Relation of minimum streamflow requirements and total summer pumping rates from managed 

wells as determined by optimal solutions to management models in the Populatic Pond area of 

the upper Charles River Basin....................................................................................................................... 85
VI Contents

TABLES

1. Hydraulic properties of stratified glacial deposits in the upper Charles River Basin, eastern Massachusetts, 

as determined by analysis of aquifer tests at public-supply wells .......................................................................... 10

2. Water withdrawals for existing and proposed municipal public-supply and large non-municipal sources 

and existing wastewater discharges in the upper Charles River Basin ................................................................... 11

3. Population on public water and sewer in eight towns in the upper Charles River Basin, 1990 and 2000.............. 17

4. Drainage-area characteristics of streamflow-measurement sites in the upper Charles River Basin....................... 20

5. Characteristics and water levels at observation wells and ponds in the upper Charles River Basin....................... 21

6. Drainage-area characteristics, streamflow, and estimated base flow at stream-gaging stations near the 

upper Charles River Basin ...................................................................................................................................... 23

7. Characteristics and water levels at long-term observation wells near the upper Charles River Basin ................... 23

8. Measured and estimated base flow at streamflow-measurement sites in the upper Charles River Basin............... 27

9. Estimated average annual water balance for the upper Charles River Basin, eastern Massachusetts, 1989–98 .... 29

10. Simulated water withdrawals and discharges in calibrated transient models and in scenarios 1 and 2 for 

existing and proposed municipal public-supply sources, large non-municipal sources, and wastewatertreatment facilities in the upper Charles River Basin ............................................................................................. 41

11. Model-calculated steady-state water levels and observed water levels at observation wells and ponds in the 

upper Charles River Basin ...................................................................................................................................... 47

12. Model-calculated steady-state stream base flow and observed stream base flow at measurement sites in the 

upper Charles River Basin ...................................................................................................................................... 50

13. Model-calculated steady-state hydrologic budget for the upper Charles River Basin............................................ 51

14. Changes in model-calculated average annual, March, and September stream base flow from 1989–98 

base flow for hypothetical scenarios of increased withdrawals and altered recharge in the upper Charles 

River Basin.............................................................................................................................................................. 65

15. Land-use characteristics, population density, and rooftop-drainage characteristics used to calculate rates 

of additional recharge in the Mine Brook aquifer area, upper Charles River Basin............................................... 73

16. Pumping wells, stream locations, and ponds in management models of the Populatic Pond area, upper 

Charles River Basin ................................................................................................................................................ 79

17. Hydrologic response coefficients for streamflow depletion at downstream observation points from 

pumping wells used to verify mass balance in the management models of the Populatic Pond area, 

upper Charles River Basin ...................................................................................................................................... 80

18. Hydrologic response coefficients for the Populatic Pond area of the upper Charles River Basin.......................... 81
Contents VII

CONVERSION FACTORS, VERTICAL DATUM, AND ABBREVIATIONS

CONVERSION FACTORS

VERTICAL DATUM

Sea level: In this report, “sea level” refers to the National Geodetic Vertical Datum of 1929 (NGVD of 

1929)—a geodetic datum derived from a general adjustment of the first-order level nets of the United States 

and Canada, formerly called Sea Level Datum of 1929.

ABBREVIATIONS

CRPCD Charles River Pollution Control District Treatment Facility

CRWA Charles River Watershed Association

DEM Digital elevation model

ET Evapotranspiration

GIS Geographic information system computer software

MADEM Massachusetts Department of Environmental Management

MADEP Massachusetts Department of Environmental Protection

MTF Milford Treatment Facility

MOVE.1 Maintenance of Variance Extension, type 1, method of correlation

No. Number

PART Computer program for streamflow partitioning (Rutledge, 1993)

PWSS Areas in the study area served by public water and septic systems

USGS U.S. Geological Survey

WMA Water Management Act

Multiply By To obtain

acre 0.4047 square kilometer

cubic foot per second (ft3/s) 0.02832 cubic meter per second

cubic foot per second per square mile [(ft3/sec)/mi2] 0.02832 cubic meter per second per square kilometer

foot (ft) 0.3048 meter

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

gallons per minute (gal/min) 0.0038 cubic meter per minute

inch (in.) 25.4 millimeter

inch per month (in/month) 25.4 millimeter

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

mile (mi) 1.609 kilometer

million gallons per day (Mgal/d) 0.04381 cubic meter per second

square foot per day (ft2/d) 0.0929 square meter per day

square mile (mi2) 2.590 square kilometer

Abstract 1

Simulation of Ground-Water 

Flow and Evaluation of WaterManagement Alternatives in the 

Upper Charles River Basin, 

Eastern Massachusetts

By Leslie A. DeSimone, Donald A. Walter, 

John R. Eggleston, and Mark T. Nimiroski

Abstract

Ground water is the primary source of 

drinking water for towns in the upper Charles 

River Basin, an area of 105 square miles in 

eastern Massachusetts that is undergoing rapid 

growth. The stratified-glacial aquifers in the basin 

are high yield, but also are thin, discontinuous,

and in close hydraulic connection with streams, 

ponds, and wetlands. Water withdrawals averaged 

10.1 million gallons per day in 1989–98 and are 

likely to increase in response to rapid growth. 

These withdrawals deplete streamflow and lower 

pond levels. A study was conducted to develop 

tools for evaluating water-management alternatives at the regional scale in the basin. Geologic 

and hydrologic data were compiled and collected 

to characterize the ground- and surface-water

systems. Numerical flow modeling techniques 

were applied to evaluate the effects of increased

withdrawals and altered recharge on groundwater levels, pond levels, and stream base flow. 

Simulation-optimization methods also were 

applied to test their efficacy for management of 

multiple water-supply and water-resource needs. 

Steady-state and transient ground-waterflow models were developed using the numerical 

modeling code MODFLOW-2000. The models 

were calibrated to 1989–98 average annual 

conditions of water withdrawals, water levels,

and stream base flow. Model recharge rates were 

varied spatially, by land use, surficial geology,

and septic-tank return flow. Recharge was

changed during model calibration by means of 

parameter-estimation techniques to better match 

the estimated average annual base flow; areaweighted rates averaged 22.5 inches per year for 

the basin. Water withdrawals accounted for about 

7 percent of total simulated flows through the 

stream-aquifer system and were about equal in 

magnitude to model-calculated rates of groundwater evapotranspiration from wetlands and ponds 

in aquifer areas. Water withdrawals as percentages 

of total flow varied spatially and temporally within 

an average year; maximum values were 12 to 13 

percent of total annual flow in some subbasins and 

of total monthly flow throughout the basin in 

summer and early fall.

Water-management alternatives were

evaluated by simulating hypothetical scenarios

of increased withdrawals and altered recharge

for average 1989–98 conditions with the flow 

models. Increased withdrawals to maximum Statepermitted levels would result in withdrawals of 

about 15 million gallons per day, or about 50

percent more than current withdrawals. Modelcalculated effects of these increased withdrawals 

included reductions in stream base flow that were 
2 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

greatest (as a percentage of total flow) in late 

summer and early fall. These reductions ranged 

from less than 5 percent to more than 60 percent

of model-calculated 1989–98 base flow along 

reaches of the Charles River and major tributaries 

during low-flow periods. Reductions in base

flow generally were comparable to upstream 

increases in withdrawals, but were slightly less 

than upstream withdrawals in areas where septicsystem return flow was simulated. Increased withdrawals also increased the proportion of wastewater in the Charles River downstream of treatment facilities. The wastewater component 

increased downstream from a treatment facility in 

Milford from 80 percent of September base flow 

under 1989–98 conditions to 90 percent of base 

flow, and from 18 to 27 percent of September

base flow downstream of a treatment facility in 

Medway. In another set of hypothetical scenarios, 

additional recharge equal to the transfer of water 

out of a typical subbasin by sewers was found to 

increase model-calculated base flows by about 12 

percent of model-calculated base flows. Addition 

of recharge equal to that available from artificial 

recharge of residential rooftop runoff had smaller 

effects, augmenting simulated September base 

flow by about 3 percent.

Simulation-optimization methods were 

applied to an area near Populatic Pond and the 

confluence of the Mill and Charles Rivers in

Franklin, Medway, and Norfolk. Water is withdrawn from six supply wells for three towns in this 

area. Management objectives in this analysis were 

to develop pumping schemes that (1) maximized 

water withdrawals while imposing specified constraints on streamflow depletion and (2) minimized streamflow depletion while meeting 

minimum requirements for water supply. Application of the optimization techniques indicated that 

hydrologic responses of pond levels and streamflow to pumping at different supply wells in the 

Populatic Pond area vary in time and duration. 

This variability suggests that water withdrawals 

could be managed to minimize the effects of 

increased withdrawals on streams and ponds. Simulation of several preliminary scenarios indicated 

the possibility that, with active management of 

water-supply sources, water withdrawals could be 

substantially increased from existing and proposed 

sources while base flow in the Charles and Mill 

Rivers was maintained above minimum average 

monthly flow requirements. Alternatively, base 

flow could be increased in the Charles River 

during low-flow periods while existing and proposed withdrawals were met. Finally, results from 

these scenarios indicated that collaborative management of water sources by towns could reduce 

the effect of withdrawals on stream base flow, and, 

by implication, would allow greater withdrawals 

from the Populatic Pond area without exceeding 

specified limits on streamflow depletion.

INTRODUCTION

Ground water is the primary source of drinking 

water for eight towns in the Upper Charles River Basin 

(fig. 1), an area of eastern Massachusetts that is undergoing rapid growth and development. More than 30 

public-supply wells withdraw water from sand and 

gravel aquifers that are high yield, but thin and discontinuous. The aquifers are in direct hydraulic connection 

with the Charles River and its tributaries, ponds, and 

wetlands, and often cross municipal boundaries. Along 

with one surface-water withdrawal and private sources, 

these wells serve a population of about 100,000, which, 

on average, has increased by about 15 percent, and

by more than 30 percent in a few towns, during the

past decade. Some of these towns also have had

large increases in commercial and industrial development in recent years. Continued rapid growth is projected, driven by an expanding metropolitan Boston 

area and the technology-based economic opportunities 

along Interstate 495 (Bouck, 1998; Massachusetts

Technology Collaborative, 1998; Metropolitan Area 

Planning Council, 2001). Towns, State agencies, business leaders, and residents are increasingly recognizing 

that water resources in the basin may fail to meet the 

needs of this continued development without adverse 

ecological and economic effects unless planning and 

management alternatives are investigated and implemented (Bouck, 1998; Massachusetts Technology

Collaborative, 1999).
3 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Echo 

Lake Lake 

Winthrop

Milford 

Pond

Louisa 

Lake

Populatic

Pond

Kingsbury

Pond

Box 

Pond

Lake 

Pearl

City Mills

Pond

Spring 

Pond

Bush

Pond

Lake 

Archer

Beaver Brook

Charles River

Charles River

Hopping B

S

rook

tall Brook

Charles River

M iscoe Brook

Bogastow

Broo k

Dirty Mead owBrook

Dopping

Brook

Jar Brook

Chicken Brook

Cress B r ook

Ea gle Brook

Miller Brook

Shepards Brook

Charles River

Mill River

Great Black Swamp

Mine Brook

Dix Brook

SHERBORN

MEDWAY

HOLLISTON

HOPKINTON

MILLIS

NORFOLK

ANP GAGE

BELLINGHAM

FRANKLIN

WRENTHAM

HOPEDALE HOPEDALE

01103280

MILFORD

BASIN BOUNDARY

TOWN BOUNDARY

TILL AND BEDROCK

STRATIFIED GLACIAL DEPOSITS

EXPLANATION

STREAM-GAGING STATION AND

 IDENTIFIER–ANP, American 

 National Power

01103280

0

0 1 2 

1 2 MILES

3 KILOMETERS

42o00'

42o

30'

73o

00' 72o

00'

71o

00'

70o

00'

41o

30'

CAPE COD 

BAY

ATLANTIC OCEAN

NANTUCKET

SOUND

AREA SHOWN

MAJOR BASIN BOUNDARIES

INDEX MAP SHOWING LOCATION OF UPPER CHARLES RIVER BASIN

MASSACHUSETTS

0 50 KILOMETERS

0 50 MILES

71o30'00"

71o

22'30"

42o15'00"

42o

07'30"

71o26'15"

42o

11'15"

42o

04'15"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 1. Location and surficial geology of the study area, upper Charles River Basin, eastern Massachusetts.
4 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Increased water withdrawals, combined with outof-basin or downstream transfers of wastewater and 

decreased natural recharge, which are associated with 

development, contribute to streamflow depletion and 

lowered water levels in ponds and wetlands in the 

upper Charles River Basin. Streamflow depletion, 

which is caused by reduced ground-water inflow, 

induced infiltration, in- and near-stream withdrawals, 

and pond-level changes, can adversely affect aquatic 

and riparian communities by reducing available water 

and altering habitat. Streamflow depletions also may 

have adverse effects on stream-water quality, because 

there is less ground-water discharge (base flow) to 

dilute contaminants from wastewater discharges and 

non-point pollution sources. These conditions and 

potential problems are typical of many suburban and 

rural communities in eastern Massachusetts and adjacent States.

In the upper Charles River Basin, streamflow 

depletion and altered pond levels already are problems 

in several areas (Bouck, 1998; Higgins, 2000;

Fiorentino and others, 2000). Reduced streamflows due 

to municipal withdrawals in the upper reaches of the 

mainstem Charles River in Milford (fig. 1) have 

resulted in poor habitat quality for aquatic life. Stratified glacial aquifers along Mine Brook in Franklin and 

Mill River in Norfolk are considered stressed by withdrawals from municipal supply wells; these withdrawals result in base-flow reductions in the rivers and 

lowered pond levels. The flow in Dopping Brook and 

Jar Brook in Holliston also may be reduced by future 

pumping. Nutrient enrichment in the Charles River and 

in-stream ponds occurs downstream of wastewatertreatment facilities in Milford and Medway. Lowered 

pond levels in Kingsbury Pond, a kettle-hole pond in 

Norfolk, have been a problem since the 1960s and have 

been associated with municipal pumping (Williams, 

1967). The stresses associated with existing watermanagement practices are likely to intensify, and new 

problem areas may develop, as a result of implementation of the multiple projects that currently are proposed. These projects include new supply wells in 

Franklin, Holliston, Norfolk, and Wrentham, a new 

surface-water withdrawal in Milford, and sewerinstallation projects in Bellingham and elsewhere. 

Nontraditional management alternatives, such as artificial recharge projects for wastewater and stormwater, 

also are proposed (DiBona and Eleria, 1999; Earth 

Tech, 1999; SEA Consultants, 1999). The location of 

some productive stratified glacial aquifers in the upper 

Charles River Basin along town boundaries makes 

implementation of some of these projects particularly 

contentious, leading to conflicts between towns over 

water rights and the adverse effects of withdrawals and 

wastewater discharges. 

A regional approach for evaluating waterresource development and management alternatives is 

needed in the upper Charles River Basin. Evaluation of 

proposed projects on a case-by-case basis does not adequately determine cumulative effects on the streamaquifer system, nor fully assess effects on downstream 

uses. In some cases, regional water-management practices may be reasonable alternatives for meeting water 

demands while minimizing adverse effects on water 

resources. A numerical ground-water-flow model is 

one tool that can be implemented at the basin-wide 

scale and used to assess the cumulative effects of largeor small-scale changes in the stream-aquifer system

on ground-water levels, pond levels, and streamflow 

depletion. When combined with optimization techniques for conjunctive management of surface- and 

ground-water systems, numerical simulation can be 

used to identify pumping strategies or other management schemes that simultaneously minimize adverse 

effects on surface waters and maximize ground-water 

yield (Barlow and Dickerman, 2001). 

Concern by State agencies and others over

the effects of population growth and development on 

water resources in the upper Charles River Basin and 

the need for a regional approach to evaluate water

management and development alternatives prompted

a study by the U.S. Geological Survey (USGS) in cooperation with the Massachusetts Departments

of Environmental Management (MADEM) and

Environmental Protection (MADEP) and the Charles 

River Watershed Association (Charles River Watershed 

Association). The objectives of this investigation were 

to define water-resource conditions within the sand and 

gravel aquifers and to evaluate the effects of selected 

ground-water development and return-flow alternatives 

on water resources in the basin through the development and application of ground-water-flow models. To 

ensure that the investigation adequately addressed 

issues of concern for water use and resources in the 

basin, representatives from towns, a water supplier, 

water dischargers, and a major water user in the study 
Introduction 5

area, along with MADEM and MADEP, participated in 

a Technical Advisory Committee that was coordinated 

by CRWA for the study. The water-use and management issues investigated in this study for the upper 

Charles River Basin are common to many other basins 

in eastern Massachusetts and adjacent States, where 

communities are striving to develop sustainable 

water-use practices.

Purpose and Scope

This report describes the development and application of numerical ground-water-flow models of the 

stream-aquifer system in the upper Charles River Basin 

in eastern Massachusetts. The study area of 105 mi2

encompasses most or all of Franklin, Holliston, Medway, and Milford, and substantial parts of Bellingham, 

Millis, Norfolk, and Wrentham (fig. 1). Steady-state 

and transient models were developed to evaluate 

regional effects of changes in the stream-aquifer 

system associated with increased water demands and 

development. Separate models were developed for 

aquifers in western and eastern parts of the study area. 

The models were developed and calibrated on the basis 

of hydrogeologic data collected during this and previous investigations. These data include lithologic information for aquifers; hydraulic properties of and 

recharge to aquifers; water levels measured in wells 

and ponds; streamflow measurements for the Charles 

River and its tributaries; and ground- and surface-water 

withdrawal rates and wastewater-discharge rates to the 

aquifers and the Charles River. Several of these data 

sets were compiled in cooperation with the CRWA. The 

models are representative of average withdrawal and 

hydrologic conditions in the stream-aquifer system for 

the 10-year period 1989–98. The USGS modular threedimensional finite-difference ground-water flow

modeling code, MODFLOW-2000 (McDonald and

Harbaugh, 1988; Harbaugh and others, 2000), was used 

to simulate the stream-aquifer system. Results of several hypothetical scenarios that incorporate changes in 

the stream-aquifer system, including increased pumping rates at existing withdrawals, new supply wells, and 

increased recharge, are described. The report also 

describes the application of simulation-optimization 

modeling in the area of Populatic Pond in Franklin, 

Medway, and Norfolk. The simulation-optimization 

model was used to investigate pumping strategies

for supply wells in this area that would minimize

streamflow depletion during critical low-flow periods. 

Description of the 

Study Area

The upper Charles River Basin encompasses an 

area of 105 mi2 in the Charles River Basin in eastern 

Massachusetts (fig. 1). The Charles River flows about 

80 mi from headwaters in Hopkinton to Boston Harbor 

and meanders extensively in its middle and lower 

reaches. The lower basin includes communities of metropolitan Boston and contains some of the most highly 

urbanized areas in Massachusetts. The upper basin, 

which is the focus of this study, is primarily suburban 

and semirural (Commonwealth of Massachusetts, 

2001). As defined for this study, the upper basin 

extends from the Charles River headwaters to approximately 1.5 river-miles downstream from the outlet of 

Populatic Pond; it also includes areas drained by 

Bogastow Brook, which joins the Charles River outside 

of the study area, and tributaries to Bogastow Brook 

upstream of Millis center (fig. 1). Most or substantial 

parts of Bellingham, Franklin, Holliston, Medway, 

Milford, Millis, Norfolk, and Wrentham are included in 

the study area. The boundary of the study area was 

chosen to include the areas of interest for water management and to coincide as much as possible with natural boundaries of aquifers, as determined from 

available data (Walker and others, 1975).

Topography in the upper Charles River Basin is 

gently rolling to hilly. Extensive low-lying and wetland 

areas are in the eastern parts and along stream valleys 

in Franklin and Holliston. The land surface generally 

slopes eastward, with elevations ranging from 550 ft 

above sea level in Milford and Hopkinton to about 120 

and 135 ft above sea level at the Charles River and 

Bogastow Brook outlets from the study area, respectively. The Charles River drops about 220 ft over about 

20 mi and is impounded by nine dams in the study

area (Myette and Simcox, 1992); dams also are 

common on the major tributaries. In addition to ponds 

and impoundments along streams, the basin contains 

several kettle ponds, including Kingsbury Pond in

Norfolk (fig. 1). Mean annual precipitation is about

46 in/yr, on the basis of data from six nearby climate 

stations (including West Medway where the average 
6 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

is 46.6 in/yr, 1957–2000; National Oceanic and

Atmospheric Administration, 2001) and is distributed 

uniformly throughout the year. 

Land use in the upper basin based on the most 

recent basin-wide data (land use, 1990; wetlands, 1998; 

MassGIS, 1997), was primarily forested (50 percent) 

and low- and medium-density residential (24 percent). 

Wetlands (8 percent), agriculture (5 percent), commercial and industrial uses (3.5 percent), and high-density 

residential use (3 percent) represented small fractions 

of land use in the upper basin in 1990. High-density 

residential, commercial, and industrial land uses in the 

study area are centered in Milford and Franklin and in 

the town centers of Medway, Millis, Holliston, and 

Bellingham. Milford is the most densely populated 

town, followed by Franklin and Medway. Since 1990, 

forested land use has decreased in area as residential, 

commercial, and industrial land-use areas have 

increased (for example, see SEA Consultants, 1999). 

Population has increased most rapidly between 1990 

and 2000 in Franklin (36 percent increase) and 

Medway (32 percent increase), and large population 

increases also have occurred in Wrentham (15 percent) 

and Holliston (14 percent) (U.S. Census Bureau data, 

from Commonwealth of Massachusetts, March 20, 

2001, and May 20, 2002). Total population in the upper 

Charles River Basin was about 100,000 in 2000. 

Previous Studies

Information on the hydrogeology and water 

resources of the upper Charles River Basin is available 

from many sources. Surficial and bedrock geology

of northern parts of the study area were mapped by

Volckmann (1975a, 1975b, 1977). Basic hydrologic 

data, including well and boring logs, water levels, and 

water quality, and the locations of high transmissivity 

zones, are available in Walker and others (1975, 1977). 

An analysis of aquifer yields on the basis of streamflow 

data was completed by Myette and Simcox (1992). 

Continuous-record streamflow data for the Charles 

River are available at a USGS stream-gaging station

in West Medway (USGS station number 01103280,

fig. 1), in operation since 1997, and at a streamgaging station operated by ENSR Corporation for 

American National Power, in Milford, since 1994 

(ENSR Corporation, 2000, and references therein). 

Streamflow data also are available for partial-record 

stations on the Charles River and its tributaries that 

were used for low-flow studies (Ries, 1993, 1994, and 

1999; Ries and Friesz, 2000). A recent assessment of 

water resources in the basin focused on water quality 

and habitat (Fiorintino and others, 2000). 

In addition to these basin-wide assessments and 

data-collection programs, numerous smaller-scale

studies have been conducted in the upper Charles River 

Basin, primarily by consultants to the towns. These 

include hydrogeologic investigations to locate watersupply sources, to determine wellhead-protection areas 

for public-supply wells, and to support specific projects 

such as wastewater discharges or large construction 

projects. Information available from these reports 

includes well and boring logs, hydrogeologic maps and 

sections, and results of aquifer tests and numerical

simulations. Reports that were used in this study 

include Amory Engineers (1989), Anderson-Nichols 

(1994, 1998), Coffin and Richardson (1977), D.L. 

Mahar (1986, 1989, 1993, 1995a, 1995b, 1996, 1997), 

Dufresne-Henry (1996, 1997), Earth Tech (1999),

Geophysical Applications (1995), Groundwater

Associates (1991), Haley and Ward (1994, 1995), 

Layne New England (1974), Metcalf & Eddy (1983a, 

1983b, 1997), SEA Consultants (1999), and Whitman 

and Howard (1991, 1992, 1996). 

Acknowledgments

The authors thank officials of Milford,

Bellingham, Franklin, Holliston, Medway, Milford, 

Medway, Norfolk, Wrentham, the Milford Water

Company, and the Milford and Charles River Pollution 

Control District Wastewater Treatment Facilities for 

their cooperation and assistance with data compilation 

and establishment of the hydrologic monitoring network. The assistance and support of the CRWA with 

data compilation, project coordination, and field work, 

particularly that of Kate Bowditch, Nigel Pickering, 

and Dudley Bonsal, also is greatly appreciated.

We specifically acknowledge the contributions of

Nigel Pickering, CRWA, in analysis of recharge and 

return-flow data and of Janet Stone, USGS, Hartford, 

Conn., in analysis of bedrock-elevation and aquiferlithologic data. Finally, we thank Mark Gerath, ENSR 

Consulting, Acton, Mass., for assistance with data 

compilation for the American National Power 

stream-gaging station in Milford. 
Upper Charles Stream-Aquifer System 7

UPPER CHARLES STREAMAQUIFER SYSTEM

Hydrogeology

Hydrogeologic units in the upper Charles River 

Basin include stratified glacial deposits, till, and bedrock (fig. 1). Stratified glacial deposits consist of sorted 

and layered sand, gravel, silt, and clay that was deposited in valleys and lowlands during the last glacial 

period (Volckmann, 1975a, 1975b). The major aquifers 

in the basin are contained in the stratified glacial deposits. Till is a poorly sorted, unstratified mixture of clay, 

silt, sand, gravel, cobbles, and boulders, deposited 

directly by the glacial ice. Till covers uplands in the 

study area and typically varies in thickness from 0 to 

about 15 ft (Volckmann, 1975a, 1975b); thicker till 

deposits underlie till hills (drumlins). Thin till layers 

also underlie stratified drift in some areas. The permeability of till is low but variable (0.6 to 2.7 ft/d for tills 

in New England; Melvin and others, 1992), and till 

generally is not used for water supply in the study area. 

Stratified glacial deposits and till are underlain by

bedrock that is mostly granite, diorite, and gneiss of 

Proterozoic and Lower Paleozoic age (Volckmann, 

1977). Bedrock outcrops are common in upland areas 

(Volckmann, 1975a, 1975b). Bedrock also is less permeable than stratified glacial deposits and yields water 

generally at low rates (for example, 10 gal/min or less); 

bedrock wells are used only for private supplies in the 

study area (Walker and others, 1975).

In addition to stratified glacial deposits, till, and 

bedrock, recent deposits in the study area include

alluvium, swamp deposits, and streambed deposits 

(Volckmann, 1975a, 1975b). Alluvium, mostly fine 

sand and silt, is along the Charles River in the eastern 

part of the study area. Swamp deposits of sand, silt, 

clay, and organic material overlie stratified glacial 

deposits and are areally extensive in low-lying areas, 

but generally form thin surficial layers. For groundwater-flow modeling purposes, alluvium and swamp 

deposits overlying stratified glacial deposits (fig. 1) 

were combined with the stratified glacial deposits and 

considered to have the same hydraulic properties. 

Stratified glacial deposits occur along the 

Charles River and its major tributaries and cover about 

50 percent of the study area (fig. 1). The areal extent of 

stratified glacial deposits in the study area was determined from published and unpublished surficial geologic maps (Volckmann, 1975a, 1975b; B.D. Stone

and J.R. Stone, U.S. Geological Survey, written

commun., 2000). The thickness of the stratified glacial 

deposits was mapped by contouring the elevation of the 

underlying bedrock or till surface (J.R. Stone, U.S. 

Geological Survey, written commun., 2000). Data on 

depth to bedrock, till, or drilling refusal were obtained 

from about 800 well logs or borings from USGS files, 

from reports by private consultants cited previously, 

and from wells installed during this study. Thicknesses 

are greatest along a north-south trending preglacial 

bedrock valley in the eastern part of the study area

(fig. 2). Thicknesses of up to 140 ft are reached in this 

area near Great Black Swamp in Millis, along the Mill 

River in Norfolk and Franklin, and near Lake Pearl in 

Wrentham. Thick stratified glacial deposits also are 

along the Charles River near Box Pond in Bellingham. 

Aquifer thicknesses elsewhere are generally 70 ft or 

less (fig. 2). Stratified glacial aquifers in the western 

and eastern parts of the study area are physically discontinuous, but are connected by the Charles River, 

which flows across thin till and bedrock in Medway 

and Franklin (fig. 1). 

The distribution, thickness, and character of the 

stratified glacial deposits reflects their depositional history as well as the configuration of the preglacial bedrock surface. The stratified glacial deposits in the upper 

Charles River Basin were deposited as a series of sedimentary sequences during temporary standstills of the 

retreating ice sheet (Stone and others, 1998; J.R. Stone, 

U.S. Geological Survey, oral commun., 2001). These 

depositional sequences were formed in valleys and proglacial lakes controlled by local topography, ice distribution, and previously deposited glacial sediments 

(Stone and others, 1998). The sequences contain icecontact, fluvial, deltaic, and(or) lacustrine sediments. 

In each sequence, grain size decreases and sorting 

increases with increasing distance from the former ice 

margin. Vertical variations in texture resulted from prograding delta sequences (coarse over fine) or from the 

deposition of fine-grained or lacustrine sediments over 

collapsed coarse-grained sediments of an earlier depositional sequence. In the upper Charles River Basin, 

where the ice-sheet retreated northward, generally eastwest trending ice-margin positions resulted in depositional sequences that grade southwards from coarse to 

fine sediments (J.R. Stone, U.S. Geological Survey, 

written commun., 2001). In the southern part of the 

study area, where valleys were shallow and narrow, the 

stratified glacial deposits are predominantly collapsed 

coarse-grained deposits. Further northward, the finegrained deposits are better developed (J.R. Stone, U.S. 

Geological Survey, oral commun., 2001).
8 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

0 70 140 FEET

EXPLANATION

AQUIFER THICKNESS

TOWN BOUNDARY

BASIN BOUNDARY

0

0 1 2 

1 2 MILES

3 KILOMETERS SHERBORN

MEDWAY

HOLLISTON

HOPKINTON

MILLIS

NORFOLK

BELLINGHAM

FRANKLIN

WRENTHAM

MILFORD

HOPEDALE

71o30'00" 71o 42 22'30" o

15'

00"

42o

07'

30"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 2. Thickness of the stratified glacial aquifers in the upper Charles River Basin, eastern Massachusetts.
Upper Charles Stream-Aquifer System 9

Information on hydraulic properties of the stratified glacial deposits is available from aquifer tests at 

public-supply wells and geologic logs from wells and 

boreholes. Results of aquifer tests at supply-well sites 

in Bellingham, Franklin, Norfolk, Holliston, and

Milford indicated horizontal hydraulic conductivities 

that ranged from 68 to 247 ft/d (table 1). This range is 

typical of coarse-grained facies of stratified glacial 

deposits, valley-fill deposits in New England (Kontis 

and others, in press; Barlow and Dickerman, 2001; 

Stone and Dickerman, 2002).

Ground water in the stratified glacial deposits 

generally is unconfined in the study area. Confined 

conditions can exist locally, however, where finegrained facies of the glacial sequences described above 

are well developed. In the eastern part of the study area, 

silt and clay sediments underlie Populatic Pond and 

extend northwestward from the pond; silt and clay sediments also are along the Mill River near Kingsbury 

Pond in Norfolk (Groundwater Associates, Inc., 1991; 

D.L. Maher Co., 1996; Dufresne-Henry, Inc., 1996). 

These fine-grained sediments result in locally confined 

conditions in the underlying sand and gravel aquifer 

sediments. In the west model area, silt, clay, and peat 

sediments near land surface result in semi-confined 

conditions in the sand and gravel aquifer in wetland 

areas along the Charles River in northeast Bellingham. 

Ground water in the stratified glacial deposits 

(henceforth, “stratified glacial aquifers” or “aquifers”) 

generally flows from boundaries with till and bedrock 

uplands through the aquifers in the downstream direction. Locally, flow is directed towards streams; streams 

in the aquifer areas are predominantly gaining when 

unaffected by pumping and are in close hydraulic connection to the aquifers (Randall and others, 1988; 

Rosenshein, 1988). Stream reaches may be losing 

under natural conditions near upland boundaries, 

where abrupt changes occur in the permeability of

the material underlying the streams and water-table 

elevations may change rapidly over short distances. 

Infiltration of streamflow also may be induced by withdrawals from pumping wells adjacent to streams and 

ponds. 

Water Withdrawals and 

Return Flow

Water is withdrawn for public supply in the 

upper Charles River Basin from 33 wells or wellfields 

in the stratified glacial aquifers and from two locations 

on the Charles River in Milford (table 2 and fig. 3).

In addition, there are several private ground- and

surface-water withdrawals in the study area that exceed 

0.1 Mgal/d and are registered or permitted under the 

Massachusetts Water-Management Act (WMA); these 

include three golf courses and a power plant (table 2 

and fig. 3). Data on locations and monthly pumped volumes from 1989 to 1998 were obtained from town officials for municipal sources and from MADEP and 

MADEM for WMA sources. Average water withdrawals, discharges, and transfers during 1989 to 1998 

are summarized in figure 4. Several proposed new 

municipal ground-water wells were being considered 

for permits by the MADEP at the time of this study; the 

proposed wells are located in Franklin, Holliston,

Medway, Norfolk, and Wrentham (table 2). 

Total withdrawals for public supply from the 

stratified glacial aquifers and surface-water sources in 

the study area averaged 9.4 Mgal/d from 1989 to

1998 (table 2). Public-supply withdrawals increased 

about 15 percent during this time (fig. 5A). On average, 

sources were pumped at less than half their MADEPapproved maximum permitted withdrawal rates on

an annual basis (table 2). Withdrawals were greatest 

during summer months of June, July, and August;

these withdrawals were about 30 percent of annual

volumes (fig. 5B). Summer withdrawals, which 

increased by about 35 percent from 1989 to 1998, 

accounted for about half of the increased volume 

pumped during this period. A small amount of the 

water withdrawn for public supply is used consumptively. Consumptive use averages about 10 percent 

annually (Solley and others, 1993). Consumptive

use varies seasonally from nearly zero in winter 

months to about 30 percent in summer months (on

the basis of an analysis of monthly water withdrawals 

in the study area; N.B. Pickering, Charles River

Watershed Association, written commun., 2001), a typical pattern resulting from outside water use (for example, lawn watering). Withdrawals varied by town, 

reflecting differences in size and population density. 
10 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

1Confining unit present at site.

Table 1. Hydraulic properties of stratified glacial deposits in the upper Charles River Basin, eastern Massachusetts, as determined by an

public-supply wells

[Aquifer test well site: See table 2 for additional identification information; site locations shown on figure 3. ft, feet; ft/d, feet per day; ft2/d, square feet per day; gal/min

available]

Aquifer test well site

Predominant

grain size of

tested interval

Length 

of test 

(days)

Well

discharge

(gal/min)

Transmissivity (ft2/d) Saturated 

thickness

(ft)

Mean

hydraulic

conductivity

(ft/d)

Specific

yield 

Mean Minimum Maxi- mum

BL-12G ........................ Fine to coarse sand and 

gravel

5 384 6,020 4,950 16,400 26 231 10.001–0.5 Grou

FR-01G,-02G ............... Fine to coarse sand and 

gravel

2 -- 8,710 7,480 9,940 60 145 -- D.L. 

FR-07G ........................ Fine to coarse sand and 

gravel

2 -- 7,860 -- -- 61 129 -- D.L. 

FR-02P ......................... Sand and gravel -- -- 10,200 8,010 11,900 -- 231 10.0003–0.1 D.L. HL-01G ........................ Medium to coarse sand 

and gravel, silt

2 190 7,550 7,032 7,990 60 126 0.07–0.14 Whit

HL-02G ........................ Medium to coarse

sand and gravel

2 240 4,850 4,520 5,070 30 162 0.13–0.15 Whit

HL-05G ........................ Medium to coarse sand 

and gravel

2 440 12,800 11,900 13,800 52 247 0.19 Whit

HL-06G ........................ -- -- -- 11,230 -- -- 47 239 -- WhitMF-03G,-04G,-05G ..... Fine to coarse sand and 

gravel

2 50 2,450 2,260 2,650 36 68 10.01 D.L. 

NF-01G ........................ Fine to coarse sand and 

gravel

17 250 5,050 2,760 9,800 53 95 10.01–0.28 Dufr

NF-01P......................... Fine to coarse sand and 

gravel

23 -- 7,220 4,200 9,800 60 120 0.04 Dufr
Upper Charles Stream-Aquifer System 11

Table 2. Water withdrawals for existing and proposed municipal public-supply and large non-municipal sources and existing 

wastewater discharges in the upper Charles River Basin, eastern Massachusetts

[Source type: GW, ground water; SW, surface water. Well depth: All wells completed in the stratified glacial deposits unless otherwise noted. Maximum 

per-mitted withdrawal rate: Data from B.R. Bouck, Massachusetts Department of Environmental Protection, written commun., 2001 (municipal sources) 

and D. LeVangie, Massachusetts Department of Environmental Management, written commun., 1999 (non-municipal sources). Proposed withdrawals include 

sources currently submitted for approval by municipalities to the Massachusetts Department of Environmental Protection; Mgal/d, million gallons per day; 

--, not applicable or not known]

Well

identifier Town Source name Source type

Well 

depth

(feet)

Mean annual 

withdrawal or 

discharge rate, 

1989–98 

(Mgal/d)

Maximum

permitted

withdrawal

rate

(Mgal/d)

Municipal Public-Supply Withdrawals

BL-05G ............................. Bellingham Well No. 5 GW 32 0.20 0.29

BL-07G ............................. Bellingham Well No. 7 GW 50 .10 .42

BL-08G ............................. Bellingham Well No. 8 GW 48 .27 .68

BL-12G ............................. Bellingham Well No. 12 GW 62 .02 .55

FR-01G ............................. Franklin Well No. 1 GW 35 1.11 .47

FR-02G ............................. Franklin Well No. 2 GW 45 1.11 .72

FR-03G ............................. Franklin Well No. 3 GW 56 .30 .32

FR-04G ............................. Franklin Well No. 4 GW 73 .67 .92

FR-05G ............................. Franklin Well No. 5 GW 60 .17 .50

FR-06G ............................. Franklin Well No. 6 GW 60 .34 .53

FR-07G ............................. Franklin Well No. 7 GW 45 .28 .58

FR-08G ............................. Franklin Well No. 8 GW 51 .30 .26

FR-09G ............................. Franklin Well No. 9 GW 30 .26 .50

FR-10G ............................. Franklin Well No. 10 GW 30 .22 .50

FR-01P .............................. Franklin Proposed well No. 11 GW -- -- .72

FR-02P .............................. Franklin Proposed Populatic Pond well GW -- -- .47

HL-01G ............................. Holliston Well No. 1 Lake Winthrop GW 50 .05 .32

HL-02G ............................. Holliston Well No. 2 Maple St. GW 40 .06 .31

HL-04G ............................. Holliston Well No. 4 Washington St. GW 50 .14 .48

HL-05G ............................. Holliston Well No. 5 Central St. GW 60 .45 .71

HL-06G ............................. Holliston Well No. 6 Brook St. GW 60 .45 .86

HL-01P.............................. Holliston Proposed well No. 7 GW -- -- .86

MD-01G............................ Medway Well No. 1 Populatic St. GW 60 .41 .38

MD-02G............................ Medway Well No. 2 Oakland St. GW 69 .14 .59

MD-03G............................ Medway Well No. 3 Village St. GW 58 .33 .60

MD-01P............................. Medway Medway Proposed well GW -- -- .43

MF-01G............................. Milford Dilla St. wells No. 1 and 2 GW 32 .02 .68

MF-02G............................. Milford Clark Island wellfield GW 34 .43 .80

MF-03G,-04G,-05G .......... Milford Godfrey Brook wells GW 52 .36 .79

MF-01S ............................. Milford Charles River SW -- .41 21.57

MF-02S ............................. Milford Echo Lake SW -- 1.77 21.57

NF-01G ............................. Norfolk Well No. 1 Gold St. GW 49 .24 .43

NF-01P.............................. Norfolk Proposed Mill River well GW -- -- 1.08

WR-02G............................ Wrentham Well No. 2 GW 51 .32 .72

WR-03G............................ Wrentham Well No. 3 GW 61 .42 .95

WR-01P............................. Wrentham Proposed well GW -- -- 1.01
12 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

1Mean annual withdrawal rate determined as one-half of reported volumes for Franklin wells No. 1 and 2. 2Maximum permitted withdrawal rate is combined rate for Charles River and Echo Lake sources. 3Maximum permitted withdrawal rate includes Northeast Energy Association wells No. 1 through 5. 4Maximum permitted withdrawal rate includes Franklin Country Club well and reservoir.

Table 2. Water withdrawals for existing and proposed municipal public-supply and large non-municipal sources and existing 

wastewater discharges in the upper Charles River Basin, eastern Massachusetts—Continued

Well

identifier Town Source name Source type

Well 

depth

(feet)

Mean annual 

withdrawal or 

discharge rate, 

1989–98 

(Mgal/d)

Maximum

permitted

withdrawal

rate

(Mgal/d)

Large Non-Municipal Withdrawals

NEA-01G,-02G, -03G Bellingham Northeast Energy Association wells 

No. 1, 2 and 3

GW -- 0.53 30.66

NEA-04G, -05G Bellingham Northeast Energy Association wells 

No. 4 and 5

GW -- .01 3.66

GECC-01S Holliston Glen Ellen Country Club Bogastow Brook SW -- .04 .16

FCC-01G Franklin Franklin Country Club well GW -- .06 4.20

FCC-01S Franklin Franklin Country Club reservoir SW -- .04 4.20

MGCC-01G Holliston Maplegate Country Club well GW -- .06 .15

Wastewater Discharges

MTF Milford Milford Treatment Facility -- -- 3.56 --

CRPCD Medway Charles River Pollution Control District 

Treatment Facility

-- -- 3.92 --

Franklin and Milford each accounted for 30 percent of 

total annual withdrawals in the study area, and the six 

remaining towns each accounted for about 10 percent 

or less. Three new sources came on-line during this 

period, Bellingham well No. 12 (1997) and Franklin 

wells No. 9 and No. 10 (1990). 

Withdrawals by large nonmunicipal sources

that are registered or permitted under the WMA were 

variable as determined from the limited available

data. Water was withdrawn for golf-course irrigation 

(“Country Club” sources, table 2) from April to

October. Some of the water withdrawn for golf-course 

irrigation likely was returned to the aquifer as recharge. 

Withdrawals for Northeast Energy Association, a 

power-generating facility, were 100 percent consumptive (N.B. Pickering, Charles River Watershed

Association, written commun., 2001) and were evenly 

distributed throughout the year. Several other WMAregistered or permitted users that are not listed in

table 2 are in the study area, but location or pumping 

data were limited or unavailable for these sources; several may not be active. Total withdrawals from 1989 to 

1998 for the non-municipal sources estimated with 

available data averaged 0.74 Mgal/d (table 2). 

Water also is withdrawn from the stratified glacial and bedrock aquifers in the study area by domestic 

users and by small private suppliers. Water withdrawn 

from low-yielding wells for domestic supply generally 

is recharged to the aquifer through on-site septic systems, minus consumptive use. Consumptive use results 

in a small net (local) loss of water to the aquifer in the 

areas served by private supply. Comparison of waterdistribution and sewer lines also indicates that there

are minimal areas of private supply that are sewered, 

especially in the areas of stratified deposits (fig. 6). 

Water is transferred in and out of the study area by 

five municipal supply systems (fig. 4). Milford Water 

Company exported on average 0.13 Mgal/d from 1989 

to 1998, or about 4 percent of its total withdrawals, out 

of the Charles River Basin to Hopedale and Mendon. 

Some of this water is returned to the Charles River Basin 

as return flow to the Milford Water Treatment Facility 

(MTF). Bellingham, whose area is in both the Blackstone and Charles River Basins, distributes water 

throughout the town from supply wells in both basins. A 

net transfer of water into the Charles River Basin results, 

which equalled about 0.26 Mgal/d, or 17 percent of total 

withdrawals for Bellingham sources, in 1996 (D.F. 

DiMartino, Town of Bellingham, oral commun., 1999). 
Upper Charles Stream-Aquifer System 13

0

0 1 2 

1 2 MILES

3 KILOMETERS

BASIN BOUNDARY

TILL AND BEDROCK

STRATIFIED 

GLACIAL DEPOSITS

SUPPLY WELL OR SURFACE-WATER

 SOURCE AND IDENTIFIER—Identifier 

 ends in "G" for existing supply well, 

 in "P" for proposed supply well, and 

 in "S" for surface-water source. See table 2.

WASTEWATER-TREATMENT FACILITY

 DISCHARGE AND IDENTIFIER—MTF, Milford

 Treatment Facility; CRPCD, Charles River

 Pollution Control District Facility 

EXPLANATION

Echo 

Lake Lake 

Winthrop

City Mills

Pond

Milford 

Pond

Louisa 

Lake

Populatic

Pond

Kingsbury

Pond

Box 

Pond

Lake 

Pearl

Spring 

Pond

Bush

Pond

Lake 

Archer

HL-02G HL-01P

HL-06G

HL-01G GECC-01S

HL-05G

MD-01G

MD-02G

MD-03G

FR-08G

FR-02P NF-01G

NF-01P

FR-05G

FR-04G

WR-02G

WR-03G

WR-01P

FCC-01S

FR-09G

FR-10G FR-01P

FR-01G

FR-02G

FR-03G FR-06G

CRPCD

MTF

BL-05G

MTF

FR-07G

MGCC-01G

BL-08G

BL-07G

BL-05G BL-12G

NEA-01G,

 02G, 03G

NEA-04G,

 05G

MF-05G

MF-01S

MF-01G

MF-02S

HL-04G

MF-03G

MF-04G

FCC-01G

MF-02G

MD-01P

Beaver Brook

Charles River

Charles River

Hopping B

S

rook

tall Brook

Charles River

M iscoe Brook

Bogastow

Broo k

Dirty Mead owBrook

Dopping

Brook

Jar Brook

Chicken Brook

Cress Brook

Ea gle Brook

Miller Brook

Shepards Brook

Charles River

Mill River

Great Black Swamp

Mine Brook

Dix Brook

71o30'00" 71o

22'30" 42o

15'

00"

42o

07'

30"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 3. Supply wells, surface-water sources, and wastewater discharges in the upper Charles River Basin, eastern 

Massachusetts.
14 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

STUDY

AR

E

A

B

O

U

NDARY

BELLINGHAM–Supply 

wells in the Blackstone 

River Basin

MILLIS–Supply wells

in the lower Charles 

River Basin

NORFOLK–Septic-system

return flow in the lower 

Charles River Basin

WRENTHAM–Septic-system

return flow in the lower Charles 

and Taunton River Basins

Water users

10.1

Wastewatertreatment facilities

MTF CRPCD

Consumptive

use

1.8

MILFORD WATER COMPANY–

Water use in the 

Blackstone River Basin

0.13

3.1

Discharges

ANP 7.5

Echo

Lake

0.40 0.26

0.38

0.13 0.32

0.19

Infiltration

~ 2 

Infiltration

UPPER CHARLES RIVER

BASIN STUDY AREA

0.015

Septic systems

Golf-course 

irrigation 0.1

Withdrawals

10.1

Stratified glacial

aquifers

Surface water

0.2 0.175

3.1

~ 5

Figure 4. Summary of water withdrawals, discharges, and transfers in the upper Charles River Basin, eastern Massachusetts, 

1989–98, in millions of gallons per day. Dashed lines indicate transfers across the study-area boundary. (ANP, American 

National Power; CRPCD, Charles River Pollution Control District Treatment Facility; MTF, Milford Treatment Facility.)

Similarly, pumping from supply wells in Millis,

Norfolk, and Wrentham results in net transfers of water 

across the study-area boundary. These transfers were 

estimated from the proportions of total land area and 

pumped volumes within and outside of the study area 

for each town. The fact that all three supply wells for 

Millis are outside of the study area results in an estimated transfer of 0.40 Mgal/d to the study area. Pumping in Norfolk and Wrentham results in estimated net 

transfers of 0.13 and 0.38 Mgal/d, respectively, out

of the study area. These small transfers into and out

of the study area balance one another, such that

water use (consumptive plus non-consumptive) by 

municipal and large non-municipal users in the study 

area (10.1 Mgal/d) is about equal to total municipal and 

large non-municipal withdrawals (fig. 4 and table 2).

Within the study area, water is redistributed 

within towns by public-water systems. Populations 

served by municipal water supplies ranged by town 

from about 50 to more than 90 percent (table 3). In 

Milford, Franklin, Medway and Millis, sewer systems 

are extensive and serve about half or more of the 

towns’ populations (table 3 and fig. 6). Much of the 

water used for municipal supply in these towns is transported to treatment facilities that discharge to the 

Charles River, minus the small amount that is used 
Upper Charles Stream-Aquifer System 15

EXPLANATION

TOTAL MUNICIPAL WITHDRAWALS

WASTEWATER DISCHARGES FROM CRPCD

WASTEWATER DISCHARGES FROM MTF

A. ANNUAL B. MONTHLY

J FM MJ J A ASOND 1989

199019911992199319941995199619971998

15

10

5

0

 MILLIONS OF GALLONS 

PER DAY

Figure 5. Average annual and monthly municipal withdrawals and wastewater discharges in the upper 

Charles River Basin, eastern Massachusetts, 1989–98: (A) Average annual municipal withdrawals and 

wastewater discharges, and (B) Average monthly municipal withdrawals and wastewater discharges. 

(CRPCD, Charles River Pollution Control District Treatment Facility; MTF, Milford Treatment Facility.)

consumptively. Wastewater from Milford is treated and 

discharged at the MTF in Hopedale, Mass., and wastewater from Medway, Franklin, Millis, and Bellingham 

is treated and discharged at the Charles River Pollution 

Control District facility (CRPCD) in Medway (fig. 3). 

About half of the water transferred into the study area 

for use in Millis also is discharged from the CRPCD; 

this estimate is based on town-wide percentages of 

sewered populations. Bellingham, Holliston, Norfolk, 

and Wrentham currently have little or no sewering 

(table 3 and fig. 6). In these towns and in unsewered 

parts of the other towns, water withdrawn from the 

sand and gravel aquifers from municipal sources or 

imported from outside the study area is recharged as 

return flow through residential and commercial septic 

systems. 

Water also is redistributed in Milford by specific 

water-supply and wastewater-management practices. 

Water is withdrawn in winter months (January,

February, and March) from the Charles River above 

Milford Pond in excess of that used for immediate 

demand and is pumped into Echo Lake Reservoir. This 

transfer averaged 0.32 Mgal/d in 1998–2000 (H.C.

Papuga, Milford Water Company, written commun., 

2000). Wastewater from the MTF is diverted to a 

power-generation facility, American National Power 

(ANP), for use as cooling water (ENSR Corporation, 

2000). This diversion is allocated by regulatory permit 

in terms of flow in the Charles River, such that flows of 

3 ft3/s or greater must be maintained at the streamgaging station immediately downstream of the MTF 

discharge. Water use at the ANP power-generating 

facility is about 80 percent consumptive; the unused 

portion is returned to the MTF. The diversion began in 

1994 and has averaged 0.43 Mgal/d annually from 

1995 to 1998 (0.19 from 1989 to 1998); monthly averages range from less than 0.1 to 0.9 Mgal/d (Mark

Gerath, ENSR Corporation, written commun., 1998). 

Wastewater return flow from the CRPCD and 

MTF averaged 3.92 and 3.56 Mgal/d, respectively, 

from 1989 to 1998 (fig. 5C). Discharges at the MTF 

were variable but essentially unchanged during that 

period, whereas discharges at the CRPCD increased 

by about 40 percent. Increased discharge from the 

CRPCD likely resulted from the sewer extensions in 

the towns served by the facility (table 3). At both treatment facilities, flows were higher in winter and spring 

months, particularly March, April, and May, than at 

other times of the year (fig. 5D). These high flows 

likely result because of an increased infiltration of 

ground water into the sewer lines during periods of

seasonally high water levels. Discharges from the

treatment facilities represent substantial fractions of 

flow in the Charles River at both discharge locations. 
16 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

TILL AND BEDROCK

STRATIFIED GLACIAL DEPOSITS

SEWER LINE

WATER-DISTRIBUTION LINE

WASTEWATER-TREATMENT-FACILITY DISCHARGE

BASIN BOUNDARY

TOWN BOUNDARY

EXPLANATION

0

0 1 2 

1 2 MILES

3 KILOMETERS

CRPCD

MTF

71o30'00" 71o

22'30" 42o

15'

00"

42o

07'

30"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 6. Water and sewer systems in the upper Charles River Basin, eastern Massachusetts.
Upper Charles Stream-Aquifer System 17

Table 3. Population on public water and sewer in eight towns in the upper Charles River Basin, eastern Massachusetts, 1990 

and 2000

[Total population from U.S. Census Bureau data, obtained from Commonwealth of Massachusetts, 2001, 2002. Public-water and -sewer population data from 

N.B. Pickering, Charles River Watershed Association, written commun., 2000]

Town

Proportion

of town in

study area,

in percent

Total population Population on public

water, in percent

Population on public sewer,

in percent

1990 2000 1990 2000 1990 2000

Bellingham............. 52 14,877 15,314 96 96 5 29

Franklin .................. 91 22,095 29,560 76 80 51 61

Holliston................. 100 12,926 13,801 97 98 0 0

Medway.................. 100 9,931 12,448 79 73 37 45

Milford ................... 86 25,355 26,799 99 94 75 80

Millis ...................... 46 7,613 7,902 85 85 40 60

Norfolk................... 31 9,270 10,460 40 58 0 0

Wrentham............... 18 9,006 10,554 80 80 0 0

At the MTF, the annual average discharge volume corresponds to 5.5 ft3/s, which is about 20 percent of the 

average annual flow at the ANP stream-gaging station 

for water years1 1996 through 1998. The annual discharge at the CRPCD corresponds to 6.1 ft3/s, about 6 

percent of average annual flow at the Medway streamgaging station for water years 1997 and 1998. During 

low-flow periods, the relative magnitudes of the wastewater discharges compared to total streamflow are considerably higher than annual averages. 

Rates of return flow from septic systems were 

determined by CWRA as part of a recharge analysis for 

this study and for a CRWA watershed model of the 

Charles River Basin (N.B. Pickering, Charles River 

Watershed Association, written commun., 2001). Water 

and sewer lines (fig. 6) were used to identify areas 

receiving public water and served by on-site septic systems (PWSS). Return-flow rates were assigned to the 

PWSS areas, which corresponded to residential, commercial, and other urban land uses. Rates were based 

on estimates of population density (5.4 people per acre 

for high-density residential and 2.2 people per acre for 

other land uses) and per capita water use (80 gallons 

per person per day). This approach is similar to that 

used in a USGS watershed model of the Ipswich River 

Basin in eastern Massachusetts (Zarriello and Ries, 

2000). Septic-system return-flow rates thus calculated 

were from 2.1 to 3.2 in/yr for commercial, from 2.2 to 

3.2 in/yr for low-density residential, and from 5.8 to 

7.8 in/yr for high-density residential land uses. 

Summed across the entire study area, septic-system 

return flow was 3.1 Mgal/d. This volume includes 

some septic-system return flow of water that originated 

from outside the study area in Bellingham, Millis, 

Norfolk, and Wrentham. The return-flow volumes calculated in this way from population density and water 

uses were checked against measured withdrawals and 

wastewater discharges, and found to be consistent.

Ground-water infiltration to sewers is another 

likely transfer of water within the study area. However, 

it is difficult to quantify. Ground-water infiltration to 

sewers is generally considered to account for as much 

as 10 to 60 percent of a town’s wastewater flow to a 

treatment facility. As described previously, groundwater infiltration to sewers is greatest during times of 

seasonally high water levels; it also varies spatially and 

is affected by physical characteristics of the sewer 

infrastructure. One way of estimating ground-water 

infiltration to sewers in the study area is to compare 

wastewater discharge rates in winter and spring, when 

water levels are high, with discharge rates in summer 

and fall. Average monthly discharges from December 

to May averaged 2.3 Mgal/d greater than average 

monthly discharges from June to November for the 

treatment facilities in the study area. This difference

is about 30 percent of total wastewater discharges and 

is probably a reasonable estimate for ground-water 

infiltration to sewers in the study area. 

1A water year begins on October 1 of the year noted and ends 

on September 30 of the following year.
18 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Streamflow and Water 

Levels

Streamflow and water levels were measured 

during this study at 24 streamflow-measurement sites, 

47 observation wells, and 9 ponds (fig. 7). Measurement sites were distributed throughout the area of stratified glacial aquifers; some were clustered in areas of 

particular interest for purposes of water-resource

management (near Populatic and Kingsbury Ponds) or 

located at study-area boundaries. Streamflow and 

water-level measurement-site characteristics are listed 

in tables 4 and 5. Many of the streamflow sites were 

partial-record stations that had been measured in previous low-flow investigations (Ries, 1993, 1994, 1997, 

1999), and historical data were available for these sites 

(table 4). Observation wells were installed in 1998 and 

1999 for this study, except for two wells in Franklin 

(F2W76 and 78) and several wells in Norfolk 

(NNW108-112), which were existing wells at aquifer 

test sites. Measurements were made monthly at most 

sites from late summer 1999 through September 2000. 

At sites in the Dopping Brook aquifer area in east

Holliston, measurements were made from December 

1999 to November 2000. Streamflow measurements 

generally were made after several days of little or

no precipitation in order to represent base-flow

conditions. Streamflow and water-level data collected 

during this study are stored in the USGS NWIS

database.

Continuous streamflow data were available for 

two sites on the Charles River in the study area, but

the data were limited. The USGS maintains a streamgaging station in Medway (station number 01103280) 

and a stream-gaging station is operated by ENSR

Corporation for American National Power (ANP) at

the MTF in Hopedale (fig. 1). The ANP stream-gaging 

station is operated primarily for low-flow data, and

its stage-discharge relation is not well defined at

moderate and high flows (ENSR Corporation, 2000; 

R.S. Socolow, U.S. Geological Survey, oral commun., 

2001). Periods of record for both stream-gaging

stations are short, extending from October 1997 to

the present (2002) for the USGS stream-gaging station 

and from July 1994 to the present for the ANP 

stream-gaging station. 

Water year 2000, during which most of the 

streamflow and water-level measurements were made, 

was about average in terms of total annual precipitation, based on records from five NOAA weather stations in and near the study area (National Oceanic and 

Atmospheric Administration, 2001). Average streamflows during water year 2000 at several nearby streamgaging stations were within 20 percent of long-term 

average flows and were neither consistently higher nor 

lower than average (table 6). Ground-water levels at 

nearby observation wells were slightly higher than 

long-term averages (table 7). Thus, hydrologic conditions during the measurement period of this study were 

not exceptionally wet or dry. Streamflow in the Charles 

River at the Medway stream-gaging station (based

on continuous stream-gaging station data; fig. 8) averaged 103 ft3/s during water year 2000, with mean

daily streamflow in August 2000 equal to 24 ft3/s.

The summer of 1999, preceding water year 2000, was 

exceptionally dry; however, this period of extended low 

base-flow conditions and low ground-water levels was 

ended by two unusually large storms in September 

(including Tropical Storm Floyd), which increased 

stream base flow and ground-water levels in the study 

area to average conditions.

Ground-water levels in the stratified glacial aquifers generally follow topography and decrease from 

west to east and from upstream to downstream locations in the study area. Measured mean annual 

ground-water levels ranged from more than 250 ft 

above sea level in headwater basins in the western part 

of the study area to less than 140 ft above sea level in 

low-lying areas draining to Bogastow Brook and along 

the Charles River near the eastern boundaries of the 

study area (table 5). Ground-water levels typically 

reached maximum levels in late spring, and water 

levels declined through summer months; this pattern 

reflects seasonal patterns in recharge and evapotranspiration (fig. 9). Measured annual fluctuations in groundwater level ranged from more than 3 ft, in wells located 

on hillsides and near boundaries of the stratified glacial 

aquifers with till and bedrock uplands (for example, 

F2W72 and HTW48; fig. 9), to about 1 ft, in wells in 

low-lying areas near streams and ponds (for example, 

MWW51 and HTW 51; figs. 7 and 9). 

Measured streamflows and water levels during 

the study period were used to estimate average flows 

for the period 1989–98 (“current average conditions”). 

Base flow is the component of streamflow resulting 

from ground-water discharge and excludes direct

runoff; this flow is the component of streamflow

that is calculated in the ground-water-flow model. 
Upper Charles Stream-Aquifer System 19

BASIN BOUNDARY

TILL AND BEDROCK

STRATIFIED GLACIAL DEPOSITS

OBSERVATION WELL AND IDENTIFIER

POND-MEASUREMENT SITE

 AND IDENTIFIER

STREAMFLOW-MEASUREMENT SITE

 AND IDENTIFIER

EXPLANATION

0

0 1 2 

1 2 MILES

3 KILOMETERS

Echo 

Lake Lake 

Winthrop

City Mills

Pond

Milford 

Pond

Louisa 

Lake

Populatic

Pond

Kingsbury

Pond

Box 

Pond

Lake 

Pearl

Spring 

Pond

Bush

Pond

Lake 

Archer

Beaver Brook

Charles River

Charles River

Hopping B

S

rook

tall Brook

Charles River

M iscoe Brook

Bogastow

Broo k

Dirty Mead owBrook

Dopping

Brook

Jar Brook

Chicken Brook

Cress Brook

Ea gle Brook

Miller Brook

Shepards Brook

Charles River

Mill River

Great Black Swamp

Mine Brook

Dix Brook

71o30'00" 71o

22'30" 42o

15'

00"

42o

07'

30"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

HTW50

HTW49

MYW60

MNW20

MNW19

MYW61

MYW58 MYW59

NNW107

NNW108

NNW109

F2W76 NNW104

F2W77

NNW112

NNW103

NNW106

XUW68

XUW65

XUW64 XUW67

NNW105

F2W74

F2W73 F2W72

F2W68

F2W75 F2W69

F2W67

A6W59

A6W63

HVW40

A6W61

A6W53

A6W64

A6W57

A6W65

MWW51

HTW47

A6W62

Lake

Winthrop

Populatic Pond

Kingsbury Pond

City Mills Pond

Bush Pond

Lake Pearl

Box Pond

Milford Pond

01103389

011033885 01103381

01103393

01103300

01103305

01103295

01103292

01103280

01103110

01103120

01103140

011032053

011032056

01103217

01103240

01103225

01103234

01103260

011032515

MNW17

A6W60

A6W55

01103206

01103235

NNW110

NNW111

01103386

HTW48

HTW51

01103210

01103140

A6W52

A6W63

Box Pond

Old Mill Pond

Figure 7. Observation wells, pond-measurement sites, and streamflow-measurement sites in the upper Charles River Basin, 

eastern Massachusetts.
20 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

1Mean daily values at streamflow-gaging station for the period indicated.

Table 4. Drainage-area characteristics of streamflow-measurement sites in the upper Charles River Basin, eastern 

Massachusetts

[Site locations shown in figure 7. mi2, square miles; --, not determined]

Station

No. Station name

Years of

historical

data

Total

number of

measurements

Drainage-area characteristics

Area

(mi2)

Area of

sand and

gravel

(percent)

Mean

slope

(percent)

Mainstem Charles River

01103110 Charles River above Milford Pond near Milford none 15 3.3 3 3.4

01103120 Charles River at Milford 1967 14 8.3 17 3.2

01103140 Charles River at Factory Pond at South Milford none 14 12.8 20 3.1

011032053 Charles River at Bellingham none 15 17.9 28 2.8

011032056 Charles River near Bellingham none 14 19.3 31 2.7

01103206 Charles River at North Bellingham 1982–83 22 21.0 35 2.6

01103260 Charles River at West Medway 1968–71, 1982–83 31 60.0 38 2.4

01103280 Charles River at Medway 1997–2000 11,141 65.2 38 2.4

01103305 Charles River near Millis 1968, 1982–83,

1989–90

31 83.8 45 2.4

Tributaries to Charles River

01103210 Stall Brook at North Bellingham 1968–70, 1983 30 3.9 43 1.9

01103217 Hopping Brook near West Medway 1968–71, 1983, 

1989–90

30 9.9 40 2.4

01103225 Miscoe Brook at Washington Street near Franklin 1968–70, 1983 26 2.2 64 2.3

01103234 Mine Brook at Franklin 1983 11 8.6 50 2.5

01103235 Mine Brook near Franklin 1967, 1983, 1989 23 2.9 52 2.5

01103240 Mine Brook near Franklin 1968–70, 1983, 

1989–91

35 14.2 50 2.5

011032515 Chicken Brook below Milk Pond near West Medway none 15 6.1 13 2.3

01103292 Eagle Brook near Wrentham 1967, 1983 23 7.5 68 2.7

01103295 Mill River below Bush Pond near City Mills none 15 9.4 72 2.8

01103300 Mill River near Norfolk 1968–71, 1983, 

1989–90

29 13.8 -- --

Bogastow Brook and Tributaries

01103381 Winthrop Canal at Lindon Street, Holliston none 11 2.7 11 1.9

01103386 Dopping Brook at Whitney Street, Holliston none 10 1.4 49 4.5

011033885 Dirty Meadow Brook at Hollis Street, Holliston none 11 2.1 12 1.8

01103389 Bogastow Brook at Central Street, Holliston none 11 12.5 38 2.7

01103393 Bogastow Brook below Great Black Swamp 

near Millis

1983 15 21.0 46 2.1
Upper Charles Stream-Aquifer System 21

Table 5. Characteristics and water levels at observation wells and ponds in the upper Charles River Basin, eastern 

Massachusetts—Continued

Well identifier

or pond name

Latitude

° ′ ′′

Longitude

° ′ ′′

Well depth

(feet below 

and surface)

Mean depth

to water 

(feet below

land surface)

Mean water-level elevation 

(feet above sea level)

Water year 2000 Estimated,

1989–98

Bellingham

A6W52.................... 42 06 15 071 29 50 17.0 7.88 229.76 229.70

A6W53.................... 42 05 35 071 28 48 57.8 40.15 212.70 212.69

A6W55.................... 42 06 31 071 27 52 17.2 6.55 208.49 208.47

A6W57.................... 42 05 46 071 27 48 15.9 11.14 220.94 220.82

A6W59.................... 42 06 09 071 28 59 28.9 8.62 215.56 215.51

A6W60.................... 42 07 20 071 27 36 30.5 13.13 205.78 205.69

A6W61.................... 42 05 54 071 29 11 51.2 35.79 219.05 218.99

A6W62.................... 42 07 47 071 28 14 24.4 12.73 210.97 210.96

A6W63.................... 42 05 36 071 28 48 29.2 10.30 220.26 220.27

A6W64.................... 42 05 20 071 28 14 26.4 9.88 224.13 224.09

A6W65.................... 42 05 50 071 27 40 22.4 6.84 203.28 203.16

Box Pond................. 42 05 41 071 29 05 -- -- 221.34 221.32

Franklin

F2W67..................... 42 06 44 071 26 06 55.4 20.17 178.90 179.01

F2W69..................... 42 05 23 071 24 24 32.4 21.51 257.99 258.08

F2W72..................... 42 03 28 071 24 30 27.3 20.33 286.57 286.51

F2W73..................... 42 03 26 071 25 38 45.2 29.99 254.31 254.38

F2W74..................... 42 02 34 071 25 45 37.5 11.02 265.48 265.50

F2W75..................... 42 05 32 071 26 03 29.0 19.96 188.51 188.51

F2W76..................... 42 07 54 071 22 58 -- 6.07 127.17 127.09

F2W77..................... 42 07 54 071 22 58 -- 5.21 126.96 126.84

Holliston

HTW47 ................... 42 11 22 071 28 09 -- 4.05 238.71 238.56

HTW48 ................... 42 13 42 071 26 09 20.0 8.88 205.61 205.24

HTW49 ................... 42 11 50 071 23 44 28.2 12.30 143.97 143.84

HTW50 ................... 42 12 22 071 25 13 12.0 6.71 167.16 167.07

HTW51 ................... 42 12 56 071 24 23 37.2 5.66 149.02 149.00

Lake Winthrop ........ 42 11 35 071 25 37 -- -- 1174.41 174.60

Hopedale

HVW40 ................... 42 06 59 071 30 22 21.8 12.42 238.32 238.36

Table 5. Characteristics and water levels at observation wells and ponds in the upper Charles River Basin, eastern 

Massachusetts

[Site locations shown in figure 7. All wells for which well depth is given are screened at bottom; screened interval equal to 5 feet. Latitude and longitude: In 

degrees, minutes, and seconds. -- not applicable or not known]
22 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

1Data for July, August, September, and November 2000 only.

2No data for January, April, and May 2000.

Medway

MNW17 .................. 42 08 45 071 25 19 26.8 9.49 183.07 183.02

MNW19 .................. 42 08 48 071 23 27 33.5 28.33 137.88 137.75

MNW20 .................. 42 09 57 071 23 44 32.8 11.57 140.86 140.71

Milford

MWW51 ................. 42 08 32 071 30 45 31.6 13.16 248.70 248.71

Milford Pond........... 42 08 39 071 30 47 -- -- 265.71 265.69

Millis

MYW58 .................. 42 09 50 071 21 42 33.5 12.62 144.75 144.71

MYW59 .................. 42 09 29 071 23 12 33.3 5.03 135.69 135.57

MYW60 .................. 42 10 59 071 23 35 21.9 7.01 138.51 138.48

MYW61 .................. 42 10 28 071 22 36 26.7 5.26 131.82 131.64

Norfolk

NNW103 ................. 42 07 19 071 22 03 29.9 5.79 130.24 130.25

NNW104 ................. 42 07 57 071 22 20 28.3 7.57 128.79 128.74

NNW105 ................. 42 05 48 071 21 12 24.1 9.63 173.17 173.31

NNW106 ................. 42 06 29 071 21 50 35.5 15.74 140.39 140.39

NNW107 ................. 42 08 02 071 22 41 -- 5.98 126.75 126.50

NNW108 ................. 42 07 52 071 22 36 -- 2.94 127.01 126.85

NNW109 ................. 42 07 18 071 22 25 -- 15.67 129.37 129.31

NNW110 ................. 42 06 48 071 22 10 -- 6.66 132.55 132.55

NNW111 ................. 42 06 46 071 22 06 -- 10.41 132.11 132.11

NNW112 ................. 42 07 31 071 22 52 -- 17.13 130.95 130.81

Bush Pond ............... 42 05 51 071 21 04 -- -- 173.77 173.84

City Mills Pond ....... 42 06 34 071 21 43 -- -- 148.87 148.83

Kingsbury Pond....... 42 07 17 071 22 13 -- -- 129.39 129.74

Populatic Pond ........ 42 07 44 071 22 42 -- -- 2126.98 127.27

Wrentham

XUW64 ................... 42 04 30 071 21 41 14.2 16.37 231.47 231.39

XUW65 ................... 42 04 15 071 21 31 43.4 25.60 205.73 205.65

XUW67 ................... 42 04 21 071 20 40 62.1 56.28 195.41 195.35

XUW68 ................... 42 03 56 071 20 37 40.2 22.20 204.20 204.20

Lake Pearl ............... 42 04 11 071 21 08 -- -- 2199.55 199.60

Old Mill Pond ......... 42 04 21 071 20 59 -- -- 194.43 194.33

Table 5. Characteristics and water levels at observation wells and ponds in the upper Charles River Basin, eastern 

Massachusetts—Continued

Well identifier

or pond name

Latitude

° ′ ′′

Longitude

° ′ ′′

Well depth

(feet below 

and surface)

Mean depth

to water 

(feet below

land surface)

Mean water-level elevation 

(feet above sea level)

Water year 2000 Estimated,

1989–98
Upper Charles Stream-Aquifer System 23

Table 6. Drainage-area characteristics, streamflow, and estimated base flow at stream-gaging stations near the upper Charles 

River Basin, eastern Massachusetts

[See Socolow and others (2001) for location information; ft3/s, cubic feet per second; in/yr, inches per year; mi2, square miles; --, not determined]

Station

No. Station name Period ofrecord

Drainage area characteristics

Area

(mi2)

Area of sand 

and gravel

(percent)

Mean slope 

(percent)

01096000 Squannacook River near West Groton, Mass. 1950–present 63.7 -- --

01097300 Nashoba Brook near Acton, Mass. 1964–present 12.2 61 2.3

01105730 Indian Head River at Hanover, Mass. 1967–present 30.3 68 2.3

01105600 Old Swamp River near South Weymouth, Mass. 1966–present 4.5 27 .5

01109000 Wading River near Norton, Mass. 1926–present 43.3 59 --

01111300 Nipmuc River near Harrisville, R.I. 1965–present 15.9 28 3.1

01175670 Sevenmile River near Spencer, Mass. 1961–present 8.7 13 5.4

Station

No.

Mean annual streamflow (ft3/s) Estimated mean annual base flow 

Period of 

record

(years)

1989–98 Water year

2000

Period of record 1989–98 Water year 2000

(in/yr) (ft3/s) (in/yr) (ft3/s) (in/yr) (ft3/s)

01096000 113 120 119 17.5 82.1 18.7 87.5 19.5 91.2

01097300 20.5 20.5 19.6 17.7 15.9 17.1 15.4 -- --

01105730 66.2 66.4 58.7 20.5 45.7 21.4 47.7 -- --

01105600 9.2 9.2 9.1 16.9 5.6 16.3 5.4 17.8 5.9

01109000 73.7 76.9 61.1 19.3 61.6 20.4 65.0 18.9 60.2

01111300 30.4 32.5 31.4 19.1 22.4 20.0 23.5 19.4 22.8

01175670 15.1 16.7 17.0 18.9 12.1 20.7 13.3 19.5 12.5

1Well open at bottom.

Table 7. Characteristics and water levels at long-term observation wells near the upper Charles River Basin, eastern 

Massachusetts

[Well location: See Socolow and others (2001) for additional location information]

Well

identifier Well location Period ofrecord

Well-screen

interval

(feet below

land surface)

Mean depth

to water

(feet below

land surface)

Mean water-level elevation 

(feet above sea level)

Period of 

record 1989–98 Water year2000

DVW10 ........ Dover, Mass. 1964–present 52–54 33.35 126.65 126.72 126.74

FXW3........... Foxborough, Mass. 1964–present 30–32 19.09 270.91 270.93 270.94

NNW27 ........ Norfolk, Mass. 1965–present 16–18 6.10 153.90 154.39 154.29

NXW54 ........ Northbridge, Mass. 1984–present 10–12 4.24 365.76 365.90 365.91

NSW21......... North Smithfield, R.I. 1947–present 116 8.04 230.64 230.87 231.24
24 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

JFM MJ J A ASOND

1997 1998

7,000

10,000

5,000

3,000

2,000

1,000

N D J FM MJ J A ASOND

1999

J FM MJ J A ASOND

2000

700

500

300

200

100

70

50

30

20

10

7

5

3

2

MEAN DAILY STREAMFLOW, IN CUBIC FEET PER SECOND

1

WATER YEAR 2000

AVERAGE FLOW, 1998–2000

AVERAGE FLOW, WATER YEAR 2000

Figure 8. Streamflow in the Charles River at the U.S. Geological Survey gaging station at Medway, Massachusetts 

(01103280) for the period November 1997 through December 2000.

Average base flows and water levels were estimated by 

(1) determining a functional relation between (a) flows 

and water levels at study measurement sites and

(b) same-day flows (mean daily) and water levels at 

nearby long-term continuous stream-gaging stations 

and observation wells; then, (2) determining current 

average conditions at the long-term sites; and, finally 

(3) using the functional relation to estimate current 

average conditions at the study measurement sites. 

Stream-gaging stations and observation wells used

in the analysis are listed in tables 6 and 7. Selected 

long-term stream-gaging stations were near the study 

area and largely unregulated. Their drainage areas

also spanned ranges of area and percent area covered 

by stratified glacial deposits similar to those for flowmeasurement sites (table 4). Selected long-term

observation wells also were nearby, screened in sand 

and gravel aquifers, and included a range of depths to 

water that was representative of the study area. 

Scatterplots of flows and water levels were used 

to determine if the relations of data at study measurement sites and at long-term sites were linear and to 

identify long-term sites that best correlated with each 

measurement site. For each measurement site, two to 

five long-term sites were identified. These included 

sites for which relations were log-linear and R2 values 

of correlations (on log-transformed data for flows) 

were greater than 0.8 (for a few sites with very good 

correlations overall, an R2 value of greater than about 

0.9 was used as the selection criterion). 

Equations relating flows and water levels at 

study and long-term sites were developed with the 

Maintenance of Variance Extension, Type 1 (MOVE.1) 

method (Hirsch, 1982). This method is applicable 

where relations are log-linear between same-day

discharges at partial-record sites and nearby streamgag-ing stations; it is commonly used to estimate lowflow statistics at partial-record sites (Bent, 1995, 1999; 

Ries and Friesz, 2000). The MOVE.1 equation is:

(1)

where

yi = streamflow or water level statistic of 

interest at the measurement site,

ym = mean instantaneous streamflow or water 

level at the measurement site,

sy = standard deviation of instantaneous 

streamflow or water level at the measurement

site,

yi ym

sy

sx

---- xi xm = + ⋅ ( ) –
25 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

284

289

285

286

287

288

F2W72

252

248

249

250

251

247

MWW51

206

211

207

208

209

210

A6W55

A SON D J FMAMJ J A SON

1999 2000

128

133

129

130

131

132

139

144

140

141

142

143

ASON D J FMAMJ J ASON

1999 2000

MNW20 NNW103

WATER LEVEL, IN FEET ABOVE SEA LEVEL

138

143

139

140

141

142

201

209

202

203

204

205

206

207

208 HTW48

147

151

148

149

150

152

HTW51 NNW106

A SON D J FMAMJ J A SON

1999 2000

Figure 9. Ground-water levels measured in selected observation wells in the upper Charles River Basin, eastern 

Massachusetts, 1999–2000.

For the continuous stream-gaging station on the 

Charles River at Medway (01103280), mean daily 

flows were used rather than mean instantaneous 

streamflow (ym in equation 1). In this study, the 

statistics of interest were mean annual or monthly base 

flow for streamflow sites and mean annual water level 

for observation wells. Base flow was chosen to 

represent streamflow because it is the base flow, or 

ground-water-discharge, component of streamflow that 

is simulated in a ground-water-flow model.

Mean annual and monthly base flow at continuous stream-gaging stations was estimated with an

automated hydrograph separation technique, PART 

(Rutledge, 1993, 1998). The computer program PART 

uses streamflow partitioning and is applied to a record 

of daily streamflows. Base flow is equated to 

xi = streamflow or water level statistic of interest 

at the long-term site,

xm = mean daily streamflow or water level at the 

long-term site, and

sx = standard deviation of mean daily streamflow 

or water level at the long-term site.
26 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

streamflows on days when antecedent and subsequent 

recession requirements are satisfied; these requirements describe the duration of substantial runoff and 

interflow after a streamflow peak. Linear interpolation 

is used to estimate base flow on days when flows do 

not meet these requirements. PART is applicable to 

basins of moderate size (about 1 to 500 mi2) with little 

or no flow regulation and where ground-water underflow at the stream measurement site is negligible. Mean 

annual base flow for continuous stream-gaging stations 

estimated with PART are listed in table 6.

Mean annual and monthly base flow for 1989 to 

1998 at continuous stream-gaging stations (table 6) 

was substituted into the MOVE.1 equations (equation 

1), along with means and standard deviations of daily 

and instantaneous streamflows for stream-gaging 

stations and measurement sites, respectively, to obtain 

estimated 1989–98 mean annual and monthly base flow 

at streamflow sites. For ground-water levels, mean 

annual water levels at long-term observation wells for 

1989 to 1998 (table 7) were used to estimate mean 

annual water levels at study observation wells (table 5). 

Two to five estimates were determined for each streamflow-measurement site and observation well; the 

number of estimates depended on the number of 

stream-gaging stations or observation wells that were 

identified as correlating well with each site or well, as 

described previously. Individual estimates from longterm sites (QBFi, in log units) were combined to generate a weighted-average base flow (QBF, in log units) or 

water level as (Ries and Friesz, 2000; K.G. Ries, U.S. 

Geological Survey, written commun., 2001):

(2)

where MSEBFi is the mean-square error (Helsel and 

Hirsch, 1992, p. 227) of the individual MOVE.1 

relations. Upper (CIhigh) and lower (CIlow) 90-percent 

confidence intervals for the weighted-average base 

flow estimates were determined as (Tasker and Driver, 

1988; K.G. Ries, U.S. Geological Survey, written 

commun., 2001): 

(3a)

and

(3b)

where MSEBF is the weighted-average MSE, in log 

units, and QBF, CIhigh, and CIlow are as defined 

previously but in arithmetic units (ft3/s). The weightedaverage MSE (MSEBF) is calculated from the meansquare error of the individual MOVE.1 relations 

(MSEBFi

) as (Ries and Friesz, 2000; K.G. Ries, U.S. 

Geological Survey, written commun., 2001):

(4)

Weighted-average base flow was retransformed back 

into arithmetic units after verification that the log 

retransformation bias was negligible (about 1 to 2 

percent as determined by the smearing estimator, 

Helsel and Hirsch, 1992). 

It should be noted that average base flow and 

water levels determined through the procedures 

described above do not take into account measurement 

errors at long-term or study-area sites, errors in the 

daily mean discharge records at the stream-gaging

stations, or errors in interpolated water levels at longterm observation wells. In addition, the average base 

flow and water-level estimates obtained for study

measurement sites with MOVE.1 relation with data 

from multiple long-term sites are not truly independent, 

because streamflow and water-level records at multiple 

long-term sites are correlated. Thus, the weighted average base flow and water-level estimates obtained with 

equations 2 and 4 may not be the best possible, and the 

actual confidence intervals are larger than those calculated with equations 3a and 3b (Ries and Friesz, 2000). 

Estimated mean annual base flow and water 

levels are given in tables 5 and 8. Estimated mean 

annual base flow at the continuous stream-gaging

station on the Charles River at Medway (01103280, 

91.0 ft3/s) was 88 percent of mean annual streamflow 

during water year 2000 (103 ft3/s). This flow value is 

consistent with estimates of the proportion of streamflow that is base flow at sites of stream-gaging stations 

in Massachusetts and Rhode Island that drain basins 

with both till and bedrock uplands and stratified glacial 

deposits (Bent, 1995, 1999). Estimated mean annual 

QBF

QBFi MSEBFi ( ) ⁄

i = 1

n

∑

1 MSEBFi ( ) ⁄

i = 1

n

∑

= ------------------------------------------------

CIlow QBF 10

1.645 MSEBF ( ) ⋅

⎝ ⎠ ⎛ ⎞ = ⁄

CIhigh QBF 10

1.645 MSEBF ( ) ⋅

⎝ ⎠ ⎛ ⎞ = ⋅

MSEBF

1

1 MSEBFi ( ) ⁄

i = 1

n

∑

= ----------------------------------------
Upper Charles Stream-Aquifer System 27

1No measurements in February or March 2000, due to ice conditions.

Table 8. Measured and estimated base flow at streamflow-measurement sites in the upper Charles River Basin, eastern 

Massachusetts

[Measured base flow: Streamflow measurements generally made after several days of little or no precipitation. August flows measured on August 28–30,

2000, unless otherwise indicated. --, not determined]

Station No.

Measured base flow, water year 2000 

(ft3/s)

Estimated base flow, 1989–98 (ft3/s)

Mean annual August

Mean annual,

water year 2000

August, water year 

2000 Flow

90-percent

confidence limits Flow

90-percent

confidence limits

Lower Upper Lower Upper 

Mainstem Charles River

01103110 2.4 0.1 2.0 1.2 3.2 0.3 0.2 0.5

01103120 4.8 .0 5.6 3.0 10.5 1.3 .7 2.5

01103140 18.7 1.2 15.8 10.6 23.7 5.6 3.7 8.3

011032053 21.4 3.4 20.1 14.7 27.5 6.4 4.6 8.7

011032056 23.3 5.5 24.2 17.2 34.0 9.2 6.5 12.9

01103206 33.0 5.6 28.8 21.4 38.9 10.4 7.7 14.1

01103260 169.8 15.4 86.9 62.4 121 25.8 18.5 35.9

01103280 -- -- 91.0 59.7 139 24.2 15.9 36.9

01103305 113 26.6 131 103 168 51.7 40.4 66.0

Tributaries to Charles River

01103210 5.5 0.8 7.6 0.8 71.7 1.1 0.1 10.1

01103217 14.1 .4 14.8 5.8 37.8 1.5 .6 3.8

01103225 3.3 .3 2.5 1.5 4.2 .5 .3 .9

01103234 -- -- 13.6 10.2 18.0 3.3 2.5 4.3

01103235 17.4 4.1 15.0 7.9 28.6 5.8 3.1 11.1

01103240 -- -- 24.4 14.9 40.0 7.7 4.7 12.6

011032515 10.2 .7 8.1 5.5 12.1 1.7 1.1 2.5

01103292 11.7 4.1 11.3 8.0 16.0 4.7 3.3 6.6

01103295 15.5 5.9 18.4 11.9 28.4 8.9 5.8 13.8

01103300 22.8 6.7 26.0 18.3 37.0 12.4 8.7 17.5

Bogastow Brook and Tributaries

01103381 -- 0.3 3.3 2.2 4.9 0.7 0.5 1.1

01103386 -- .1 1.0 .5 2.0 .1 .1 .3

011033885 -- .2 1.6 1.1 2.2 .3 .3 .4

01103389 -- 1.6 12.8 8.0 20.5 2.9 1.8 4.6

01103393 -- 1.6 30.2 20.9 43.6 6.2 4.3 9.0

base flow at other study flow-measurement sites ranged 

from 73 to 122 percent of measured mean annual base 

flow, based on monthly measurements during baseflow conditions in water year 2000; measured values 

all were within 90-percent confidence limits of estimated values. Comparison of estimated mean monthly 

base flow with individual streamflow measurements 

indicates that the monthly estimates match seasonal 

variations and the magnitudes of measured flows; 

exceptions are due to precipitation events or snowmelt, 

such as the February 2000 measurements at Charles 

River sites (fig. 10). Mean annual water levels estimated with MOVE.1 and data from long-term observation wells for 1989–98 (table 5) were nearly identical 

with measured water levels during water year 2000 (the 

median difference was less than 0.05 ft). 
28 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

C. 01103305 D. 01103393

50

0

150

50

100

FLOW, IN CUBIC FEET PER SECOND

JA SON D J FMAM J J A SON

2000

0

300

100

200

0

100

50

JA SON D J FMAM J J A SON

1999 1999 2000

0

100

B. 01103217

0

6

2

4

PRECIPITATION,

IN INCHES

0

6

2

4

MEASURED STREAMFLOW

ESTIMATED BASE FLOW

A. 01103206

Figure 10. Daily precipitation and measured streamflow during base-flow conditions, 1999–2000, and estimated monthly 

average base flow at selected streamflow-measurement sites in the upper Charles River Basin, eastern Massachusetts: 

(A) Charles River at North Bellingham (01103206), (B) Hopping Brook near West Medway (01103217), (C) Charles River 

near Millis (01103305), and (D) Bogastow Brook below Great Black Swamp near Millis (01103393). Precipitation data from 

the National Oceanic and Atmospheric Administration (2001).

Water Balance

A steady-state water balance for the upper 

Charles River Basin describes the stream-aquifer 

system and provides a conceptual framework for the 

ground-water flow models (table 9). In the steady-state 

water balance, inflows to and outflows from the combined stream-aquifer system are identified and quantified on an average annual basis. Stratified glacial 

aquifers and their associated streams as well as upland 

till and bedrock areas were included. The water balance 

described here includes the components simulated in

a ground-water-flow model; thus, the direct runoff 

component of streamflow is not included. 

Inflows include ground-water recharge from precipitation, ground-water recharge from septic-system 

return flow, and wastewater discharge to streams from 

treatment facilities (fig. 11). Outflows include stream 

outflows from the study area, withdrawals from supply 

wells and surface-water intakes, ground-water

evapotranspiration (ET), and infiltration into sewer 

lines (fig. 11). Water moves from till and bedrock areas 

in uplands to areas of stratified glacial aquifers in

valleys (fig. 1) through lateral ground-water inflow and 

streamflow from uplands (fig. 11). 

Ground-water recharge from precipitation is

the primary source of water. Precipitation recharge

(fig. 11) is the component of precipitation that 
Upper Charles Stream-Aquifer System 29

1Values correspond to recharge rates of 20–26 inches per year.

infiltrates the land surface and reaches the water table; 

it is equal to total precipitation minus direct surfacewater runoff and ET at or near land surface. Precipitation recharge rates vary temporally and spatially with 

climatic conditions, land-surface permeability and 

slope, vegetation, and soil saturation. Precipitation 

recharge rates, which can be estimated with a number 

of methods, were compiled from literature sources

for the water balance. Basin-wide precipitationrecharge rates have been estimated from long-term 

streamflow records to range from 23.8 to 25.5 in/yr

for basins in eastern Massachusetts and Rhode Island 

that contain about equal areas of stratified glacial 

deposits and till and bedrock uplands (Bent, 1995;

Barlow, 1997; Barlow and Dickerman, 2001); these 

values were about 52 to 54 percent of average annual 

precipitation. Rates from streamflow records in basins 

in southern New England with varying percentages of 

stratified glacial deposits and basin slopes ranged from 

17.5 to 28.1 in/yr (Bent, 1995; 1999). A range of

precipitation-recharge rates, 14 to 28 in/yr (reflecting a 

range of specific-yield values) was estimated from 

ground-water levels at a long-term observation well

in the Charles River Basin downstream of the study 

area (Myette and Simcox, 1992). Recent estimates of 

precipitation recharge to aquifers on western Cape 

Cod, where surface runoff is negligible, are about

26 in/yr, or 60 percent of average annual precipitation 

(Masterson and others, 1998). Precipitation-recharge 

rates to stratified glacial deposits, which may be higher 

than basin-wide rates in basins with a mix of till and 

stratified deposits, ranged from 19.3 to 29.4 in/yr at 

locations throughout New England and New York

(Barlow, 1997). For the upper Charles River Basin, 

recharge estimates of 20 to 26 in/yr (43 to 56 percent of 

average annual precipitation at West Medway), which 

are typical rates based on literature sources, correspond 

to inflow volumes of 100 to 130 Mgal/d to the study 

area (table 9).

Other inflows to the study area (fig. 11) are small 

in comparison with precipitation recharge (table 9). 

Septic-system return flow was estimated for 1989–98, 

at 3.1 Mgal/d, as described previously, or about 2 to 3 

percent of precipitation recharge. Wastewater return 

flow to streams from 1989–98 was 7.5 Mgal/d, or 

about 5 to 6 percent of precipitation recharge. Total 

inflows to the study area, thus, were estimated at 111 to 

141 Mgal/d (table 9).

The primary outflow from the study area is 

streamflow at the study-area boundary. Streamflow 

leaves the study area through the Charles River and 

Bogastow Brook in Norfolk and Millis (stations 

01103305 and 01103393, fig. 7). Streamflow at these 

locations includes direct runoff (stormflow) and base 

flow (ground-water discharge to streams, fig. 11); however, stormflow runoff, which is not simulated in 

ground-water-flow models, is not included in the water 

balance. Average annual base flow in the Charles River 

and Bogastow Brook at the study-area boundary was 

estimated at 104 Mgal/d (161 ft3/s, table 9). Confidence 

limits (90 percent) around these estimates yielded a 

range from 80 to 137 Mgal/d. Withdrawals for water 

supply are another outflow from the stream-aquifer 

system (fig. 11). Water-supply withdrawals were equal 

to 10.1 Mgal/d in 1989–98 (table 9). 

Ground-water ET, ground-water underflow, and 

infiltration to sewer lines are additional outflows from 

the stream-aquifer system that are difficult to quantify 

accurately. Ground-water ET, or ET directly from the 

water table (fig. 11), may occur in discharge areas such 

as along streams, in wetlands, or in other areas where 

the water table is close to land surface. On a basin-wide 

basis, ground-water ET is likely to be small relative to 

the total water balance because the areas where it can 

occur are small. Ground-water ET was estimated to be 

1 to 2 in/yr on the basis of streamflow records for about 

Table 9. Estimated average annual water balance for the 

upper Charles River Basin, eastern Massachusetts, 1989–98

[Numbers in parentheses correspond to 90-percent confidence intervals for 

base-flow estimates]

Water-balance component

Rate of flow

Million

gallons

per day

Cubic feet

per second

Inflow

Recharge from precipitation................ 1100–130 1155–201

Septic-tank return flow........................ 3.1 4.8

Wastewater discharge to streams......... 7.5 11.6

Total inflow...................................... 111–141 171–217

Outflow

Ground-water discharge to streams 

(stream base flow) ........................... 104 

(80–137)

161 

(124–212)

Water-supply withdrawals ................... 10.1 15.6

Ground-water evapotranspiration 

in wetlands....................................... 5.0 7.7

Total outflow.................................... 119 

1(95–152)

184

1(147–235)
30 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

A

B

D

1

2

4 C

C

F E

3

F

A

A

Water table

INFLOWS AND INTERNAL FLOWS

A. Precipitation recharge to uplands

B. Precipitation recharge to stratified glacial aquifers

C. Inflows from uplands to stratified glacial aquifers

D. Seepage from streams

E. Wastewater discharge

F. Septic-tank return flow

Till

Bedrock

Stratified glacial deposits

HYDROGEOLOGIC UNIT

Not to scale

OUTFLOWS

1. Discharge to streams

2. Ground-water evapotranspiration

3. Down-valley ground-water 

underflow (flow direction out of page)

4. Water withdrawals

Till and Bedrock Uplands Stratified Glacial Aquifers

EXPLANATION

Figure 11. Schematic diagram showing hydrogeologic units and flow components of the water balance for the upper 

Charles River Basin, eastern Massachusetts. (Modified from Randall and others, 1988, fig. 1.)

150 basins in the east-central United States; these estimates were calculated as the difference between average annual recharge and base-flow rates (Rutledge, 

1998). By means of a similar approach, ground-water 

ET was determined to be 2 to 3 in/yr for two basins

in southern Rhode Island. (Barlow and Dickerman, 

2001). One way of estimating ground-water ET from 

available data for the upper Charles River Basin is to 

use wetland areas, as mapped with the land-use data, 

and regional evaporation rates. Wetland areas represent 

areas where the water table is near land surface, and, 

thus, where ground-water ET may occur. Free-watersurface evaporation rates can be considered potential 

evaporation rates from a saturated surface and may be 

used to approximate actual ET rates, which are 

unknown (Farnsworth and others, 1982). Application 

of a growing-season free-water-surface evaporation 

rate of 21 in/yr (determined at a regional scale for eastcentral Massachusetts, Farnsworth and others, 1982)

to wetland areas in the basin yields a basin-wide 

ground-water ET estimate of about 1 in/yr. This estimate corresponds to an outflow rate of about 5 Mgal/d 

(table 9). This estimate, which is about 5 percent of 

streamflow and withdrawal outflows, indicates that 

ground-water ET is likely to be a small but substantial 

component of the water balance. Total estimated outflow from the stream-aquifer system, including the 

estimate of ground-water ET, was 119 Mgal/d, and 95 

to 152 Mgal/d represented the 90-percent confidence 

limits on the base-flow estimates (table 9).

Ground-water underflow (fig. 11) may occur as 

flow through the stratified glacial aquifers across studyarea boundaries. Ground-water underflow is unlikely to 

occur where Bogastow Brook exits the study area, 

because the stratified glacial deposits pinch out at the 

stream in this area. Near the Charles River in Norfolk, 
Simulation of Ground-Water Flow 31

however, stratified glacial deposits are thin but do not 

entirely pinch out at the study-area boundary, and 

ground-water underflow out of the study area may 

occur. Ground water also may flow in or out of the 

study area where the boundary is a surface-water divide 

that does not coincide with a ground-water divide. 

Generally, the assumption that ground-water and

surface-water divides coincide is valid in the upper 

Charles River Basin because of the rolling topography, 

shallow aquifers, and close hydraulic connection 

between surface and ground water; the assumption

also is supported by surface-water and land-surface 

elevations near the study-area boundary. Pumping 

wells may alter these conditions, however. In Millis, 

supply wells in stratified glacial deposits adjacent to 

the study area, along with the flat topography in this 

area, may result in some ground-water outflow across 

the study area boundary. In Norfolk, a supply well may 

draw water into the study area from stratified glacial 

deposits that extend across a surface-water divide east 

of the Charles River. The magnitude of these flows 

could not be quantified. However, they are likely to be 

small relative to the total water balance, given the magnitude of supply-wells withdrawals relative to total

outflows.

Infiltration of ground-water into sewer lines is a 

potential outflow in the water balance that is hard to 

quantify, but probably is small compared to the total 

water balance. Ground-water infiltration changes seasonally with water-table fluctuations and varies greatly 

with the age and integrity of sewer-system infrastructure. Ground-water infiltration represents some fraction 

of the total wastewater flows from sewage-treatment 

facilities, with estimates ranging from 10 to 60 percent. 

Few data were available for towns in the upper Charles 

River Basin, but ground-water infiltration and inflow of 

stormwater runoff into sewer lines averaged 26 percent 

of total wastewater flows for six towns in an adjacent 

basin (Earth Tech, 2001). Similarly, comparison of 

winter and summer wastewater-discharge rates in the 

basin suggested an infiltration rate of about 30 percent 

of total wastewater flows, as described previously.

At this rate, ground-water infiltration in the upper 

Charles River Basin would be about 2 Mgal/d, or

about 2 percent of total study-area outflows (table 9). 

Within the study area, lateral ground-water 

inflow and stream base flow from uplands (upland 

inflows; fig. 11) accounts for a large fraction of water 

inflows to the stratified glacial aquifers. Upland flows 

originate as precipitation recharge in till and bedrock 

areas. Recharge rates in these areas are poorly known, 

but may be as little as one-half of the rates in nearby 

stratified glacial aquifers (Kontis and others, in press). 

If a range of upland recharge rates from 50 to 100 percent of rates in stratified glacial aquifers is assumed, 

upland flows account for about 35 to 50 percent of total 

recharge inflows to stratified glacial aquifers. 

Overall, estimated water inflows—precipitation 

recharge and wastewater return flow—and outflows—

stream base flow at study-area boundaries, water withdrawals, ground-water ET, and ground-water infiltration to sewers—agree reasonably well (within 20 

percent; table 9), given the uncertainties in individual 

estimates. Ground-water outflows that were not 

directly quantified in the water balance, such as 

ground-water underflow, may account for some of the 

discrepancy between inflows and outflows. 

SIMULATION OF 

GROUND-WATER FLOW

Steady-State Numerical 

Models

Ground-water flow and head distribution across 

the upper Charles River Basin were simulated with the 

three-dimensional, finite-difference ground-water flow 

modeling code, MODFLOW-2000 (Harbaugh and

others, 2000). Ground-water flow and heads were only 

simulated in the stratified glacial aquifers. Upland till 

and bedrock areas, which constitute about 50 percent of 

the study area, were not simulated, although stream 

base flow and recharge from these areas were included. 

Active model areas are shown in figure 12. Physically 

discontinuous stratified glacial aquifers in west and 

east parts of the basin were simulated with separate 

models. The two models use the same discretization 

(cell size and orientation). Model-calculated flow in the 

Charles River at the exit from the west model is specified as inflow in the east model. Both steady-state and 

transient conditions were simulated. The steady-state 

models simulated average flow conditions during the 

10-year period, 1989–98. Areal recharge rates and 

hydraulic conductivity were estimated with inverse 

modeling (based on nonlinear regression) to obtain the 

best possible fit between observed and simulated head 

or flow values (Hill and others, 2000) and to calibrate 

the steady-state models. 
32 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

TILL OR BEDROCK

STRATIFIED GLACIAL DEPOSITS

MODEL GRID—Active model cells only 

shown except for first and last rows and columns 

BOUNDARY OF ACTIVE AREA OF MODEL

0

0 1 2 

1 2 MILES

3 KILOMETERS

BOUNDARY BETWEEN

 WEST AND EAST MODEL AREAS

EXPLANATION

BASIN BOUNDARY

150100150200250300325last

1

50

100

150

200

250

300

350

400

413

last

EAST AREA 

MODEL

WEST AREA 

MODEL

EAST AREA 

MODEL

WEST AREA 

MODEL

Figure 12. West and east model areas, grid, and boundary of active model areas for the simulation models of the 

upper Charles River Basin, eastern Massachusetts.
Simulation of Ground-Water Flow 33

Spatial Discretization

The model domain was discretized into a grid of 

420 rows and 325 columns of cells with uniform horizontal dimensions of 200 ft (fig. 12). West and east 

aquifer models are discretized with one model grid that 

contains both areas; in the west model, cells representing stratified glacial deposits of the east aquifer area are 

inactive, and opposite in the east model. The vertical 

discretization consists of two sloping layers of variable 

thickness (fig. 13). The bottom of the lower layer corresponds to the top of the bedrock and till surface, and 

was determined by intersecting the model grid with an 

elevation surface that was linearly interpolated from 

elevation contours of the bedrock and till surface. The 

bottom of the upper layer was set at 40 ft below a simulated water-table surface that was determined by means 

of uncalibrated, two-dimensional versions of the models. Where saturated thickness, based on the preliminary simulated water table, was less than 40 ft, only the 

upper layer is active, extending vertically to the bedrock surface. There were 11,197 active cells in the west 

model, (9,868 cells in layer 1 and 1,329 cells in layer 2) 

and 13,824 active cells in the east model (9,724 cells in 

layer 1 and 4,100 cells in layer 2). 

Horizontal boundaries of the active model were 

defined so as to exclude areas where the saturated 

thickness of the stratified glacial aquifers was thin

(fig. 12). Thinly saturated areas along the margins of 

the active model area (at the boundary with till and 

bedrock uplands) lead to numerical instabilities and the 

inability of the model to converge to a solution. Thinly 

saturated areas were identified using available data on 

aquifer geometry and preliminary model testing. Preliminary model testing also was used to identify and 

solve numerical problems with the models. Various 

two- and three-dimensional versions of an initial 

model, in which west and east areas were both active, 

were numerically unstable and did not yield steadystate solutions. A transient version was tested to

investigate a quasi-steady-state solution (defined in 

terms of negligible change in storage) as an alternative 

to a steady-state model, but also was unstable. Consequently, the initial model was simplified by simulating 

the study area with two, two-dimensional models for 

west and east aquifer areas, for which robust solutions 

were obtained. 

Boundary Conditions

Physically, the active model areas are bounded 

laterally by till and bedrock uplands or by areas of

thin stratified glacial deposits; they are bounded below 

by bedrock or till (figs. 12 and 13). In the numerical 

models, the lower boundaries of the active model areas 

are no-flow boundaries (fig. 13). Flow from adjacent 

uplands is simulated as specified flow along the lateral 

boundaries. No-flow boundaries are specified along 

ground-water drainage divides near downgradient 

boundaries of the model domain and in small areas 

near the headwaters of some tributary streams. In these 

areas, ground-water divides were determined from surface topography. The assumption that ground-water 

divides coincide with surface-water divides likely is 

valid in the upper Charles River Basin under natural 

conditions, because of the shallow aquifer depths and 

close connection between surface and ground water; 

however, there may be local deviation between surfaceand ground-water divides where pumping wells are 

close to divides. The study area boundary was chosen 

to minimize the use of surface-water divides as model 

boundaries, but this use could not be entirely avoided 

because of the interfingering geometry of stratified

glacial aquifers in the basin.

Streams are hydrologic boundaries within the 

models. Streams were simulated as head-dependent 

flow boundaries with a version of the Stream Routing 

Package (Prudic, 1989) that is compatible with

MODFLOW-2000. This package simulates hydraulic 

interaction between the aquifer and adjoining

streams and tracks the amount of water in each simulated stream. Simulated streams are divided into 

reaches, corresponding to individual model cells,

and segments, which are groups of reaches that are 

connected in downstream order (Prudic, 1989).

Water may flow in either direction at the boundary 

between the aquifer and streams. Flow, or leakage,

is calculated by multiplying the head difference 

between the stream stage and the underlying aquifer

by the specified streambed conductance. Stage in

the stream may be either specified or calculated.

In the upper Charles River Basin models, stream stage 

was specified. If the model-calculated head in the

aquifer is higher than the stream-stage elevation, 

ground-water discharges to the stream reach; when the 

model-calculated head in the aquifer falls below the 
34 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Areal recharge Stratified

Thick stratified

Stratified

ACTIVE MODEL CELLS—Thick stratified

glacial deposits

SIMULATED MODEL AREA

SIMULATED STREAM

SIMULATED POND

Till

Bedrock

NO-FLOW BOUNDARY TILL

THIN STRATIFIED GLACIAL DEPOSITS

BEDROCK

EXPLANATION

Till

Figure 13. Schematic section showing active model cells and model layers for the simulation models of upper Charles 

River Basin, Massachusetts.

stream-stage elevation, stream water leaks to the aquifer. The quantity of water in each stream reach is 

tracked and routed to downstream reaches. Simulated 

streams may go dry when stream leakage to the aquifer 

exceeds inflows from upstream reaches. 

Simulated streams in the models included the 

Charles River and Bogastow Brook and their major 

tributaries: Beaver Brook, Stall Brook, Hopping

Brook, Miscoe Brook, Dix Brook, Mine Brook, 

Chicken Brook, Shepards Brook, Eagle Brook, Mill 

River, Dirty Meadow Brook, Dopping Brook, and Jar 

Brook (figs. 1, 12A and B). Small headwater streams 

that drain substantial areas of upland also were 

included.

Stream-stage elevations used as input to the 

model were determined from DEM data (Elassal and 

Caruso, 1983) and measured elevations. The DEM

data provided elevations (based on 3-meter contour 

intervals) of stream surfaces at a grid resolution of 30 

meters. These data were linearly interpolated to generate surface elevations at a finer scale to more accurately 

represent elevations along stream channels. The resulting elevations for stream channels were manually 

edited to ensure that stream elevations decreased monotonically downstream throughout the active model 

areas and to ensure agreement with measured elevations at outlets from major impoundments along 

streams. A uniform water depth of 1 ft and streambed 

thickness of 5 ft were assumed to calculate streambed 

top and bottom elevations from stream-surface

elevations. 

Streambed conductances were determined initially from literature sources and assumed stream 

geometries. Conceptually, streambed conductance is 
Simulation of Ground-Water Flow 35

calculated from streambed hydraulic properties and 

geometry as (Prudic, 1989):

(5)

In practice, the values used to compute streambed 

conductance, particularly the hydraulic conductivity 

and thickness of the streambed, are poorly known. 

Vertical hydraulic conductivities of streambed 

sediments have been reported to range from 0.1 to

15 ft/d, with typical values of 1 to 2 ft/d, for streams 

draining stratified glacial deposits at locations in 

central Massachusetts, southern Rhode Island,

Long Island, New York, and northern New Jersey 

(Rosenshein, 1968; DeLima, 1991; Prince and others, 

1988; Dysart and Rheaume, 1999). In the upper 

Charles River Basin models, initial streambed 

conductances were calculated from assumed uniform 

stream widths and streambed thicknesses of 10 ft and

5 ft, respectively, stream lengths equal to the model cell 

size of 200 ft, and a streambed hydraulic conductivity 

of 5 ft/d. Streambed conductances were varied during 

model calibration, as described below.

Streamflow entering the active model area was 

specified at two locations where the Charles River 

enters the west and east models. Flow in the first reach 

of the segment representing the Charles River in the 

west model was specified at 2.3 ft3/s. Flow in the first 

reach of the segment representing the Charles River in 

the east model was specified equal to the simulated 

flow in the Charles River at the exit from the west 

model. Flow from the headwaters to the Charles River 

was specified to facilitate simulation of water withdrawals from Echo Lake Reservoir, which is located in 

this upland subbasin. 

Stresses

Recharge

Recharge to the models consisted of

infiltration from precipitation and wastewater return 

flow from septic systems. Recharge was applied to

the active model area as a spatially varying, specified 

flux to the uppermost active layer (fig. 13). In general, 

precipitation recharge varies spatially with land-surface 

permeability, which is a function of soil characteristics 

and land use. Recharge to the aquifer from septicsystem return flow occurs in areas where wastewater 

disposed of through septic systems originated as water 

withdrawn from the aquifer at a distant location, such 

as a public-supply well. Septic-system return flow in 

areas of co-located, private water supplies is assumed 

to result in no net inflow to the aquifer and was not 

explicitly represented in the models. 

The spatial distribution of recharge rates in the 

upper Charles River Basin was determined by CRWA 

from 1991 land-use data (1:25,000 scale; MassGIS, 

1997), the areal extent of stratified glacial deposits 

(Volckmann, 1975a, 1975b; B.D. Stone and J.R. Stone, 

U.S. Geological Survey, written commun., 2000), and 

digital data layers of water and sewer systems (fig. 6). 

Digital water- and sewer-system data was based on 

maps provided by each town. Twenty recharge categories were defined with seven types of land use (forest, 

open space, wetland, water, high-density residential, 

low-density residential, and commercial), two surficialgeology types to represent soil conditions (stratified 

glacial deposits and till or bedrock), and the presence 

or absence of septic-system return flow to develop a 

range of recharge rates for the study area. The recharge 

categories and associated rates were based on those 

used in a watershed model of the Ipswich River Basin 

in northeastern Massachusetts (Zarriello and Ries, 

2000). Recharge rates (excluding rates to wetlands, 

which were zero) ranged from 7.8 in/yr for commercial 

land-use areas to 28.9 in/yr for high-density residential 

land use in stratified glacial aquifers with return flow. 

where

CSB = streambed conductance, in ft2/d,

KSB = vertical hydraulic conductivity of the 

streambed, in ft/d,

LS = length of the stream, equal to the model cell 

dimension in the direction of streamflow, 

in ft,

WS = width of the stream, in ft, and

TSB = thickness of the streambed, in ft.

CSB

KSB LS WS ⋅ ⋅

TSB

= -------------------------------
36 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Both the area-weighted average rate for the study area, 

16.4 in/yr, and the rates corresponding to natural precipitation recharge (forested land cover), 10 or

17.5 in/yr for areas of till and bedrock or stratified glacial deposits (fig. 1), respectively, were low relative to 

recharge rates determined from literature sources and 

those used in the water balance. Thus, it was expected 

that recharge rates based on the CRWA analysis and 

Ipswich model would be increased during model calibration. To facilitate this process, recharge rates in 

model cells of the active model area were divided by 

the natural precipitation recharge rate for stratified glacial aquifers (17.5 in/yr); thus, a multiplier array was 

used to incorporate the spatial variability in recharge 

rate into the models. 

Flow from upland areas typically enters areas of 

stratified glacial aquifers through lateral ground-water 

inflow (unchanneled flow) or as flow in streams draining upland subbasins (channeled flow) (Randall and 

others, 1988). Inflows from upland areas through both 

these pathways were simulated as specified flow to the 

uppermost model layer at active model boundaries 

(except for flow from the headwaters of the Charles 

River, which was included as input the stream package). Upland areas were divided into subbasins to 

account for these inflows (fig. 14A and B). Subbasins 

were defined for each of the simulated streams on the 

basis of topography from points where the streams met 

active model boundaries. Subbasins were defined for 

upland areas between simulated streams on the basis of 

previously defined subbasin boundaries (MassGIS, 

1997) and aquifer geometry. Using the digital data 

layer of recharge polygons, an area-weighted recharge 

rate for each upland subbasin was calculated. This rate, 

applied to the entire subbasin area, was summed to 

apply as a specified flow to one or more model cells 

along the active model boundary (one cell for subbasins drained by streams; multiple cells for subbasins 

drained by lateral ground-water flow). This flow rate 

was added to the precipitation recharge rate for each 

boundary model cell, and converted to a recharge multiplier as described previously. The recharge data layer, 

adjusted to model-calibrated recharge rates (described 

below), is shown in figure 15.

The approach of using recharge multipliers to 

simulate boundary flows from the till and bedrock 

uplands to the active model area had advantages and 

disadvantages. It was used so that rates of inflow to the 

active model area from uplands, consisting of either 

lateral ground-water inflow or base flow in upland 

streams, could be varied, along with recharge rates to 

active model areas, during parameter estimation for 

model calibration. These inflows, which conceptually 

originate as precipitation recharge in uplands, were 

expected to be large components of the total hydrologic 

budgets for the stream-aquifer systems. Thus, it was 

important to vary the inflows from uplands when calibrating the model to obtain the best possible match 

between observed and model-calculated stream base 

flows. The approach required delineation of upland 

basins drained by channeled and unchanneled flow in 

uplands, as do other approaches for simulating 

recharge in uplands (for example, Kontis, 2001). Addition of base flow in streams draining uplands as a flow 

directly to the aquifer was based on the assumption that 

all base flow in these streams is lost to the aquifer 

immediately after the streams cross the upland/aquifer 

boundary. 

Ponds drained by streams and wetlands in the 

upper Charles models were conceptualized as areas of 

net ground-water discharge, where, on average, water 

levels were equal to or less than the surrounding water 

table. Saturated soils during most of the year and lowpermeability sediments likely result in no net recharge 

of water to aquifers from precipitation in most wetlands, under natural conditions. Similarly, precipitation 

onto surfaces of ponds drained by streams becomes 

direct stream (stormwater) runoff rather than contributing to ground-water recharge or constituting a component of the stream base flow that is simulated in the 

models. Water leaves these areas through ET, however. 

Thus, ponds drained by streams and wetlands were 

simulated as areas of net water loss from the streamaquifer system. Water loss from these areas in the 

active model area was simulated as a loss equal to the 

rate of free-water-surface evaporation in the growing 

season in east-central Massachusetts, or 21 in/yr (Farnsworth and others, 1982). Although vegetation in wetlands may result in actual ET rates that differ from freewater-surface evaporation rates, the direction and magnitude of these differences is uncertain; thus, freewater-surface evaporation rates are likely to be a good 

approximation of actual ET rates in wetlands (Mitsch 

and Gosselink, 1993). In uplands, ET from ponds 

drained by streams and wetlands was not simulated. 

However, wetland and pond areas in uplands are small 

relative to non-wet areas and ET from these areas was 

considered negligible. 
Simulation of Ground-Water Flow 37

0

0 1 2 

1 2 MILES

3 KILOMETERS

A. WEST MODEL AREA

UPLAND SUBBASIN DRAINED 

 BY CHANNELED INFLOW

UPLAND SUBBASIN

 DRAINED BY LATERAL 

 GROUND-WATER INFLOW

SIMULATED WETLANDS

SIMULATED PONDS

SIMULATED STREAMS

BOUNDARY OF THE ACTIVE MODEL AREA

BASIN BOUNDARY

BOUNDARY OF UPLAND RECHARGE SUBBASIN

SITE OF SPECIFIED INFLOW TO STREAM SEGMENT

EXPLANATION

SITE OF STREAMFLOW EXIT FROM ACTIVE MODEL AREA

71o30'00" 71o

22'30"

42o07'30"

42o

11'15"

42o

03'45"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 14. Upland recharge subbasins and simulated surface-water features for the simulation models of the upper 

Charles River Basin, eastern Massachusetts: (A) West model area, and (B) East model area.
38 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

0

0 1 2 

1 2 MILES

3 KILOMETERS

B. EAST MODEL AREA

71o22'30"

42o07'30"

42o

11'15"

42o

03'45"

71o

26'15"

From USGS and MassGIS data sources, Geographic Projection, 

Spheroid GRS 1980, Datum NAD 83 

UPLAND SUBBASIN DRAINED 

 BY CHANNELED INFLOW

UPLAND SUBBASIN

 DRAINED BY LATERAL 

 GROUND-WATER INFLOW

SIMULATED WETLANDS

SIMULATED PONDS

SIMULATED STREAMS

BOUNDARY OF THE ACTIVE MODEL AREA

BASIN BOUNDARY

BOUNDARY OF UPLAND RECHARGE SUBBASIN

SITE OF SPECIFIED INFLOW TO STREAM SEGMENT

EXPLANATION

SITE OF STREAMFLOW EXIT FROM ACTIVE MODEL AREA

Figure 14. Upland recharge subbasins and simulated surface-water features for the simulation models of the 

upper Charles River Basin, eastern Massachusetts: (A) West model area, and (B) East model area—Continued.
Simulation of Ground-Water Flow 39

0

10 to14

14 to 18

18 to 22

22 to 26

26 to 30

30 to 34

34 to 38

POND

AREAL RECHARGE

RATE, IN INCHES

PER YEAR

EXPLANATION

0

0 1 2 

1 2 MILES

3 KILOMETERS

BASIN BOUNDARY

71o30'00" 42o

15'

00"

42o

07'

30"

71o

22'30"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 15. Areal recharge rate from precipitation recharge and septic-system return flow for the simulation models of the 

upper Charles River Basin, eastern Massachusetts. Rates are based on the calibrated steady-state models.
40 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Ponds that are not drained by streams, typically 

kettle-hole ponds, have no outlet for surface-water 

runoff of precipitation. These areas were simulated as 

ground-water recharge areas and received simulated 

recharge of 22 in/yr in the active model area. This value 

is a representative recharge rate from literature sources 

cited previously and represents the estimated difference 

between precipitation and evaporation in the basin. 

Recharge rates for kettle ponds, and recharge rates representing ET from in-stream ponds and wetlands, were 

not varied in the calibration process.

Water Withdrawals and Discharges

Pumping wells, water withdrawals from surface 

water, and wastewater discharges to streams were simulated with the Well Package for MODFLOW-2000. 

Withdrawals from pumping wells were simulated as 

specified flows from the aquifer; surface-water withdrawals and discharges were simulated as specified 

flows from the model cells hydraulically connected to 

streams. Flow rates in the steady-state model (table 10) 

were equal to average annual withdrawal and discharge 

rates for 1989–98 for most sources (table 2). For irrigation (country club) withdrawals, flow rates were specified as 50 percent of the total withdrawal rates, to 

account for return flow (infiltration) of irrigation water. 

All simulated wells were in layer 1 for both models,

in accordance with available data about well-screen 

intervals (table 2). 

Several water withdrawals were outside the 

active model area. These were the town of Franklin 

supply well No. 9 and Franklin Country Club well and 

reservoir in the east model area and the Echo Reservoir 

in the west model area. The Franklin withdrawals were 

incorporated into the east model by reducing the inflow 

of water from the upland subbasin in which they were 

located by a volume equal to their combined annual 

average withdrawals (table 10). In the west model, 

inflow from the upland subbasin that includes the Echo 

Lake Reservoir was altered to exclude all flow that 

originated upstream of the reservoir outlet. This 

approach is based on the assumption that the reservoir, 

which stores water pumped to it in winter from a location further downstream on the Charles River, does

not have substantial stream releases during base-flow 

conditions. Water use from Echo Lake Reservoir is represented by the winter withdrawals for storage from 

the downstream Charles River location, which averaged 0.34 Mgal/d annually, based on available data 

(1998–2000).

Hydraulic Properties

Aquifer hydraulic properties required for the 

steady-state simulations were horizontal and vertical 

hydraulic conductivity. Aquifer-test data at publicsupply wells (table 1) and geologic logs from wells and 

boreholes were used to characterize hydraulic conductivities. The conceptual model of aquifer-sediment deposition also provided a framework for estimating 

hydraulic conductivity throughout the model area. A 

generalized approach was used to estimate the spatial 

distribution of hydraulic conductivity and values were 

adjusted during model calibration to match observed 

data.

Hydraulic conductivity was assigned in zones 

that were delineated based on aquifer and stream 

geometry and streamflow-measurement sites and were 

modified during model calibration (fig. 14A and B). 

Horizontal hydraulic conductivity (Kh) ranged from 

150 to 290 ft/d in the west model, and from 70 to

220 ft/d in the east model. Vertical hydraulic conductivity (Kv) was based on published ratios of vertical to 

horizontal conductivities, which generally range from 

1:3-1:5 for coarse-grained stratified glacial deposits to 

1:30-1:100 for fine-grained deposits (Dickerman and 

others, 1990; Masterson and Barlow, 1997, and references therein; Masterson and others, 1998; Stone and 

Dickerman, 2002). Vertical hydraulic conductivities 

ranged from 10 to 70 ft/d (Kv:Kh equal to 1:4 to 1:17) 

in the west model and from 5 to 30 ft/d (Kv:Kh equal to 

1:6 to 1:20) in the east model. Horizontal and vertical 

hydraulic conductivities were distributed in the same 

way in layers 1 and 2, except for cells in layer 1 that 

simulated ponds. Model areas used to simulate ponds 

were assigned horizontal and vertical hydraulic conductivity values of 10,000 ft/d. Simulating ponds as 

active model cells allowed pond levels to change with 

changing stresses in the aquifer. The large value of 

hydraulic conductivity effectively simulates the lack of 

resistance to flow through the ponds and results in little 

or no head drop across adjacent pond cells and in

realistic flow patterns in the aquifer surrounding the 

ponds.
Simulation of Ground-Water Flow 41

Table 10. Simulated water withdrawals and discharges in calibrated models and in scenarios 1 and 2 for existing and proposed 

municipal public-supply sources, large non-municipal sources, and wastewater-treatment facilities in the upper Charles River 

Basin, eastern Massachusetts—Continued

Identifier

Model location Simulated withdrawal or discharge rate (Mgal/d)

Layer Row Column

1989–98 Scenario 1 Scenario 2

Annual Summer Annual Summer Annual Summer

West Model

Charles River

MF-01S1 .......................... 1 163 61 0.84 0.12 0.84 0.12 0.84 0.12

MF-01G........................... 1

1

166

166

58

59

.008

.008

.009

.009

.34

.34

.34

.34

.34

.34

.34

.34

MF-02G........................... 1 179 56 .43 .46 .80 .80 .80 .80

MF-03G........................... 1 215 61 .12 .11 .26 .26 .26 .26

MF-04G........................... 1 217 62 .12 .11 .26 .26 .26 .26

MF-05G........................... 1 218 63 .12 .11 .26 .26 .26 .26

MTF................................. 1 240 66 (3.56) (2.72) (4.87) (3.75) (4.87) (3.75)

BL-12G............................ 1 278 130 .02 .04 .41 .35 .41 .38

Beaver Brook Subbasin

BL-05G............................ 1 278 98 .20 .20 .29 .29 .29 .29

NEA-01G, -02G, -03G .... 1 279 105 .53 .50 .66 .61 .66 .61

NEA-04G, -05G .............. 1 279 97 .010 .010 .014 .013 .014 .013

Stall Brook Subbasin

BL-07G............................ 1 234 123 .10 .15 .21 .37 .21 .37

BL-08G............................ 1 233 125 .27 .46 .44 .68 .44 .68

Hopping Brook Subbasin

HL-04G ........................... 1 137 117 .14 .24 .14 .24 .14 .24

HL-05G ........................... 1 98 215 .45 .47 .45 .47 .45 .47

Mine Brook Subbasin

FR-01G............................ 1 317 188 .11 .20 .24 .40 .19 .34

FR-02G............................ 1 319 189 .11 .20 .24 .40 .19 .34

FR-03G............................ 1 350 175 .30 .35 .32 .32 .32 .32

FR-06G............................ 1 341 182 .34 .39 .50 .53 .46 .53

FR-07G............................ 1 253 163 .28 .31 .48 .55 .41 .49

FR-10G............................ 1 376 202 .22 .25 .39 .48 .32 .41

FR-01P ............................ 1 384 161 .00 .00 .00 .00 .47 .47

MGCC-01G..................... 1 261 161 .030 .060 .075 .154 .075 .15

East Model

Charles River

FR-08G............................ 1 216 231 0.30 0.322 0.26 0.26 0.26 0.26

FR-02P ............................ 1 217 234 .00 .00 .00 .00 .72 .72

MD-01G .......................... 1 206 228 .41 .46 .41 .46 .37 .38

MD-03G .......................... 1 211 233 .33 .38 .34 .38 .37 .47

CRPCD............................ 1 212 239 (3.92) (3.52) (5.80) (5.85) (6.78) (6.33)

Table 10. Simulated water withdrawals and discharges in calibrated models and in scenarios 1 and 2 for existing and proposed 

municipal public-supply sources, large non-municipal sources, and wastewater-treatment facilities in the upper Charles River 

Basin, eastern Massachusetts

[Identifier: See table 2 for additional identification information; site locations shown in figure 3; identifiers ending in “G,” “S,” and “P” denote ground-water, 

surface-water, and proposed sources, respectively. Simulated withdrawal or discharge rate: Parentheses denote discharges. Average withdrawal or 

discharge rate for 1989–98: Rates used in calibrated steady-state and transient models. Average summer withdrawal or discharge rate: Average of 

monthly average June, July, and August rates. Mgal/d, million gallons per day; --, not applicable]
42 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

1Includes water transfers to Echo Lake Reservoir.

2Withdrawals simulated as reduced inflow from upland basin.

East Model—Continued

Bogastow Brook Subbasin

HL-01G ........................... 1 116 184 0.05 0.13 0.05 0.13 0.049 0.13

HL-02G ........................... 1 46 157 .06 .16 .056 .15 .057 .15

HL-06G ........................... 1 50 191 .45 .52 .45 .52 .45 .52

HL-01P ............................ 1 46 190 .00 .00 .00 .00 .35 .46

GECC-01S....................... 1 118 221 .02 .07 .08 .23 .08 .23

MD-02G .......................... 1 188 226 .14 .19 .14 .19 .16 .26

MD-01P........................... 1 165 221 .00 .00 .00 .00 .41 .43

Mill River Subbasin

FR-04G............................ 1 243 245 .67 .74 .92 .92 .85 .92

FR-05G............................ 1 255 249 .17 .20 .32 .40 .26 .34

FR-09G2 .......................... -- -- -- .26 -- .44 -- .37 --

NF-01G............................ 1 224 273 .24 .30 .38 .42 .24 .30

NF-01P ............................ 1 252 252 .00 .00 .00 .00 .26 .31

WR-02G .......................... 1 326 275 .32 .41 .40 .51 .32 .41

WR-03G .......................... 1 326 278 .42 .49 .52 .62 .42 .49

WR-01P........................... 1 351 274 .00 .00 .00 .00 .19 .23

FCC-01G, -01S2 .............. -- -- -- .05 -- .10 -- .10 --

Table 10. Simulated water withdrawals and discharges in calibrated models and in scenarios 1 and 2 for existing and proposed 

municipal public-supply sources, large non-municipal sources, and wastewater-treatment facilities in the upper Charles River 

Basin, eastern Massachusetts—Continued

Identifier

Model location Simulated withdrawal or discharge rate (Mgal/d)

Layer Row Column

1989–98 Scenario 1 Scenario 2

Annual Summer Annual Summer Annual Summer

Model Calibration

The steady-state models were calibrated by varying model input parameters—recharge, horizontal and 

vertical hydraulic conductivity, and streambed conductance—to obtain as close a match as possible between 

simulated and observed water levels and stream base 

flows. Calibration data consisted of stream base flow at 

23 measurement sites and water levels at 35 observation wells and 7 ponds. The models were calibrated to 

average annual conditions during 1989–98. This 10-

year period was chosen so that the calibrated models 

would reflect both current stresses and average hydrologic conditions. The observed values used in calibration consisted of estimated average annual stream base 

flow and water levels for the 1989–98 period.

A combination of inverse modeling and trialand-error methods was used to calibrate the steadystate models. In inverse modeling, model input parameters are solved for as unknowns, with automated mathematical methods (Hill, 1998) to determine parameter 

values that yield the best match between observed and 

simulated values. In trial-and-error methods, parameter 

values are adjusted manually until the match between 

observed and simulated values is considered reasonable. For the upper Charles models, inverse modeling 

was used to determine areal recharge rates, on the

basis of observed stream base flows. This approach

was based on the conceptual model of the study area 

and was supported by the water balance, in which all 

water enters the model area as recharge and exits as 

streamflow (except for minor outflows such as underflow, assumed to be negligible). The approach was
Simulation of Ground-Water Flow 43

considered appropriate for recharge because of the 

available streamflow data and the lack of long-term 

data from which recharge rates could be independently 

calculated. Trial-and-error methods were used to calibrate the models with respect to aquifer hydraulic properties. 

Recharge

Calibrated recharge rates were determined 

through parameter estimation with a version of 

UCODE, an inverse modeling code, that is incorporated into MODFLOW-2000 (Hill, 1998; Hill and

others, 2000). In UCODE, nonlinear regression is

used to determine parameter values that minimize 

values of a weighted least-squares (WLS) objective 

function. The WLS objective function quantifies the 

difference between observed and simulated values. The 

parameter-estimation process incorporates an analysis 

in which model sensitivities for the estimated parameter are determined at all specified observation points. 

The primary estimated parameter was the natural

precipitation-recharge rate for stratified glacial aquifers 

with forested land use. Cell-by-cell recharge rates used 

in the model were the product of this rate and a multiplier that incorporated the variations resulting from

differences in land use, surficial geology, and septicsystem return flow. Thus, adjustments to the natural 

precipitation-recharge rate made during the parameterestimation process also proportionately adjusted all 

spatially varying recharge rates specified in the model. 

Inflows to the active model area from uplands also 

were adjusted similarly, because these inflows were 

simulated through the recharge-multiplier values at 

active model boundaries. This approach allowed the 

automated adjustments during parameter estimation 

while preserving the spatial variations in recharge that 

reflected local land-surface permeability and return 

flow throughout the model area (fig. 15).

Recharge parameters were defined to correspond 

to natural precipitation-recharge rates in two aquifer 

zones of the east model. The northern zone included 

the aquifer associated with Dopping Brook, Jar Brook, 

and Bogastow Brook in Holliston, Millis, and Medway. 

The southern zone included the aquifer associated

with the Mill River and the Charles River in Wrentham, 

Franklin, Norfolk, and Medway. The two zones are 

separated by a ground-water divide and are drained 

by streams to different model exits. Observed values 

were stream base flow at four measurement sites in

the northern zone and five sites in the southern zone 

(sites shown in fig. 7). Uniform values of horizontal 

and vertical hydraulic conductivity of 250 ft/d and

50 ft/d, respectively, were used in the model during the 

estimation of the recharge parameter. For streambed 

conductance, a uniform value of 2,000 ft2/d was used, 

which corresponds to a streambed vertical hydraulic 

conductivity of 5 ft/d, as previously described. 

The rate of natural precipitation recharge in 

undeveloped sand and gravel areas used in the calibrated west and east models was 24.1 in/yr. Rates of 

natural precipitation recharge determined by parameter 

estimation were 23.9 in/yr and 24.3 in/yr, for northern 

and southern zones, respectively, of the east model. 

These rates (weighted by area) were averaged to yield 

24.1 in/yr for natural recharge precipitation for the east 

model. This rate yielded simulated stream base flows 

that closely matched observed flows in the east model 

(see below). A value of 24.1 in/yr also yielded a close 

match between simulated and observed stream base 

flow in the west model. Because of this close agreement and because natural precipitation recharge rates 

would not be expected to vary from east to west model 

areas, the calibrated value from the east model also was 

used in the west model. The natural precipitation 

recharge rate of 24.1 in/yr corresponded to spatially 

varying rates ranging from 10.7 in/yr for commercial 

land use areas to 38 in/yr for high-density residential 

land use in stratified glacial aquifers with return flow. 

Area-weighted rates across the model domain were 

21.5 and 23.0 in/yr for active areas of the west and east 

models, (excluding wetland areas and ponds drained by 

streams), respectively, and 21.9 and 24.0 in/yr for flows 

to the active model area from till and bedrock uplands 

of west and east model areas, respectively. The basinwide area-weighted rate was 22.5 in/yr. These rates are 

consistent with recharge rates determined for stratified 

glacial deposits and basins with a mix of stratified glacial aquifers and till and bedrock uplands in previous 

studies (Bent, 1995, 1999; Barlow, 1997; Barlow and 

Dickerman, 2001; Masterson and others, 1998; Myette 

and Simcox, 1992). 
44 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Hydraulic Properties

Hydraulic properties were adjusted manually

to obtain a close match between observed and simulated ground-water and pond levels. The degree of fit 

between observed and simulated values was measured 

with normal and weighted residual statistics, which

are calculated by MODFLOW-2000. Twenty welland pond-observation points were used in the west 

model and 30 points were used in the east model. 

Hydraulic properties adjusted during model calibration 

were horizontal and vertical hydraulic conductivity

and streambed conductance. 

Horizontal and vertical hydraulic conductivity 

values were adjusted in zones (fig. 16A and B), from 

initial values of 250 ft/d and 50 ft/d, respectively. Initial 

zones were modified to better match observed heads by 

adding separate zones in the east model for the Jar 

Brook aquifer in Holliston and for the eastern part of 

Millis. Calibrated horizontal hydraulic conductivity 

values (70 ft/d) were lower in these areas than in adjacent zones. Available lithologic data were sparse in 

these areas, but did indicate the presence of finegrained sediments, particularly in the Millis zone. As 

described previously, final calibrated values of horizontal hydraulic conductivity ranged from 150 to 290 ft/d 

in the west model and from 70 to 220 ft/d in the east 

model; calibrated values of vertical hydraulic conductivity ranged from 10 to 70 ft/d in the west model and 

from 5 to 30 ft/d in the east model (fig. 16A and B). 

Streambed conductances also were adjusted to 

improve the match between observed and simulated 

heads and flows. Streambed conductances were 

adjusted from initial values of 2,000 ft2/d in groups that 

corresponded to hydraulic conductivity zones and areas 

simulated as wetlands and ponds. In areas simulated as 

wetlands, calibrated streambed conductances generally 

were lower than in other areas, ranging from about 500 

to 2,000 ft2/d in both models; these values correspond 

to KSB values of about 1 to 5 ft/d, using reasonable 

assumptions about streambed thickness and stream 

width. In areas simulated as ponds, streambed conductances generally were higher than in other areas,

ranging from about 2,000 ft2/d to high values that

simulated no resistance to vertical flow. In areas not 

simulated as wetlands or ponds, streambed conductances generally ranged from about 1,000 to 7,000 ft2/d 

(KSB about 3 to 15 ft/d) in the west model and from 

about 5,000 to 15,000 ft2/d (KSB about 10 to 30 ft/d)

in the east model. Additional changes were made to 

streambed conductances to resolve numerical instabilities during development of transient versions of the 

model. Most changes were slight, affecting a few 

model cells only, except for changes made in several 

headwater streams in the southern part of the east 

model. In these areas, streambed conductances were 

lowered substantially, to about 30 ft2/d (KSB about

0.1 ft/d), which is still, however, within the range of 

typical values. 

Calibration Results

Calculated water levels for observation wells and 

ponds for the calibrated steady-state model are shown 

in table 11 and figure 17A. The mean absolute difference between observed and model-calculated water 

levels (mean absolute water-level residual) was 2.6 ft 

for the west model and 2.0 ft for the east model; both 

these values are about 1.4 percent of the total head 

(water-level) change across the simulated water table 

for the west (180 ft) and east (142 ft) models. In some 

cases, larger differences between observed and modelcalculated water levels may have resulted from the 

position of observation wells near the model boundaries (for example, well MNW17 in the west model; 

table 11), where model discretization or simulation of 

boundary stresses may have led to less accurate simulated water levels. The simulated water table for east 

and west model areas (fig. 18) is consistent with the 

conceptual model of flow in the stratified glacial aquifers. Water-table contours decrease in the downstream 

direction in tributary valleys and along the Charles 

River, and bend at large streams reflecting the effect of 

ground-water discharge. 

Calculated base flow at streamflow-measurement sites for the calibrated steady-state models is 

shown in table 12 and figure 17B. Observed and modelcalculated stream base flow at or near model exits

(stations 01103260 for the west model and stations 

01103305 and 01003393 for the east model) are in 

good agreement, differing by less than 4 percent. The 

mean difference between model calculated and 

observed base flow (mean absolute base-flow residual) 

was 1.6 ft3/s for the west model and 3.4 ft3/s for the 

east model, or 2.3 and 3.4 percent of the total range of 

observed flows in west and east model areas, respectively. At individual measurement sites, mean absolute 

base-flow residuals ranged from about 3 to 60 percent 

of observed values; large differences (greater than 30 

percent) generally were associated with sites on small 

streams with low flows (table 12). 
Simulation of Ground-Water Flow 45

HYDRAULIC CONDUCTIVITY ZONES—Upper 

 number is horizontal hydraulic 

 conductivity and lower number is vertical

 hydraulic conductivity, in feet per day

INACTIVE MODEL AREA

POND CELLS—Horizontal and vertical

hydraulic conductivity simulated as 10,000 feet per day

BOUNDARY OF THE ACTIVE MODEL AREA

BASIN BOUNDARY

260

 50

250

 40

250

50

290

70

270

 50

250

40

150

10

170

20

170

20

160

10

220

40

250

40

150

15

290

70

EXPLANATION

0

0 1 2 

1 2 MILES

3 KILOMETERS

A. WEST MODEL AREA

71o30'00" 71o

26'15"

42o07'30"

42o

11'15"

42o

03'45"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 16. Hydraulic conductivity zones for the simulation models of the upper Charles River Basin, eastern 

Massachusetts: (A) West model area, and (B) East model area.
46 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

HYDRAULIC CONDUCTIVITY ZONES—Upper 

 number is horizontal hydraulic 

 conductivity and lower number is vertical

 hydraulic conductivity, in feet per day

INACTIVE MODEL AREA

POND CELLS—Horizontal and vertical

hydraulic conductivity simulated as 10,000 feet per day

BOUNDARY OF THE ACTIVE MODEL AREA

100

10

290

70

100

5

130

20

70

5

130

20

70

5

185

20

180

30

220

30

EXPLANATION

B. EAST MODEL AREA

0

0 1 2 

1 2 MILES

3 KILOMETERS

71o26'15"

71o22'30"

42o07'30"

42o

11'15"

42o

03'45"

From USGS and MassGIS data sources, Geographic Projection, 

Spheroid GRS 1980, Datum NAD 83 

BASIN BOUNDARY

Figure 16. Hydraulic conductivity zones for the simulation models of the upper Charles River Basin, 

eastern Massachusetts: (A) West model area, and (B) East model area—Continued.
Simulation of Ground-Water Flow 47

Table 11. Model-calculated steady-state water levels and observed water levels at observation wells and ponds in the upper 

Charles River Basin, eastern Massachusetts—Continued

Well identifier or 

pond name

Model location Average annual water level 

Layer Row Column Model calculated

(feet above sea level)

Observed 

(feet above sea level)

Difference

(model calculated minus

observed, in feet) 

West Model

A6W52.................... 1 269 80 224.99 229.70 -4.71

A6W53.................... 1 289 103 214.62 212.69 +1.93

A6W55.................... 1 261 124 208.79 208.47 +.32

A6W59.................... 1 272 98 216.25 215.51 +.74

A6W60.................... 1 236 130 199.65 205.69 -6.04

A6W61.................... 1 280 94 218.85 218.99 -.14 

A6W62.................... 1 223 116 214.81 210.96 +3.85

A6W63.................... 1 285 94 220.17 220.27 -.10

A6W65.................... 1 282 129 205.85 203.16 +2.69

F2W67..................... 1 254 164 177.27 179.01 -1.74

F2W72..................... 1 354 200 287.07 286.51 +.56

F2W73..................... 1 355 175 248.55 254.38 -5.83

F2W74..................... 1 381 172 260.60 265.50 -4.90

F2W75..................... 1 291 165 189.99 188.51 +1.48

HTW47 ................... 1 114 118 240.04 238.56 +1.48

HVW40 ................... 1 247 67 234.81 238.36 -3.55

MWW51 ................. 1 200 59 250.38 248.71 +1.67

MNW17 .................. 1 193 182 173.44 183.02 -9.58

Box Pond................. 1 290 88 221.38 221.32 +.06

Milford Pond........... 1 176 57 265.44 265.69 -.25

East Model

F2W76..................... 1 219 235 129.96 127.09 +2.87

HTW48 ................... 1 44 163 204.34 205.24 -.90

HTW49 ................... 1 99 225 141.62 143.84 -2.22

HTW50 ................... 1 82 184 165.31 167.07 -1.76

HTW51 ................... 1 66 202 149.84 149.00 +.84

MNW19 .................. 1 192 224 141.15 137.75 +3.40

MNW20 .................. 1 157 217 137.76 140.71 -2.95

MYW58 .................. 1 160 263 140.28 144.71 -4.43

MYW59 .................. 1 171 229 138.09 135.57 +2.49

MYW60 .................. 1 125 221 138.46 138.48 -.02

MYW61 .................. 1 141 143 126.09 131.64 -5.55

NNW103 ................. 1 237 256 129.46 130.25 -.79

NNW104 ................. 1 218 256 125.10 128.74 -3.64

NNW105 ................. 1 283 275 173.68 173.31 +.37

NNW106 ................. 1 262 260 148.96 140.39 +8.57

Table 11. Model-calculated steady-state water levels and observed water levels at observation wells and ponds in the upper 

Charles River Basin, eastern Massachusetts

[Site locations shown in figure 7. Observed water level: Estimated for 1989–98 from measurements made during 1999–2000, as described in text]
48 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

East model—Continued

NNW107 ................. 1 215 241 130.49 126.50 +3.99

NNW108 ................. 1 220 243 128.77 126.85 +1.92

NNW109 ................. 1 237 247 130.20 129.31 +.89

NNW110 ................. 1 252 253 134.50 132.55 -1.95

NNW112 ................. 1 231 237 131.27 130.81 +.46

XUW65 ................... 1 330 268 205.02 205.65 -.63

XUW67 ................... 1 237 250 196.93 195.35 +1.58

XUW68 ................... 1 339 288 200.58 204.20 -3.62

Bush Pond ............... 1 282 276 173.72 173.84 -.12

City Mills Pond ....... 1 263 264 149.06 148.83 +.23

Kingsbury Pond....... 1 237 250 130.18 129.74 +.44

Lake Pearl ............... 1 340 278 200.06 199.60 -.46

Lake Winthrop ........ 1 116 180 176.14 174.60 +1.54

Old Mill Pond ......... 1 330 278 194.59 194.33 +.26

Populatic Pond ........ 1 219 239 128.68 127.27 +1.41

Table 11. Model-calculated steady-state water levels and observed water levels at observation wells and ponds in the upper 

Charles River Basin, eastern Massachusetts—Continued

Well identifier or 

pond name

Model location Average annual water level 

Layer Row Column Model calculated

(feet above sea level)

Observed 

(feet above sea level)

Difference

(model calculated minus

observed, in feet) 

300

250

200

100

150

MODEL-CALCULATED WATER-LEVEL 

ELEVATION, IN FEET ABOVE SEA LEVEL

WEST MODEL

EAST MODEL

LINE OF EQUALITY

EXPLANATION

100 150 200 250 300

OBSERVED WATER LEVEL ELEVATION,

IN FEET ABOVE SEA LEVEL

MODEL-CALCULATED STREAM 

BASE FLOW, IN CUBIC FEET PER SECOND 20

40

60

80

100

120

140

0

0 20 40 60 80 100 120 140

OBSERVED STREAM BASE FLOW,

IN CUBIC FEET PER SECOND

A. WATER LEVELS B. STREAM BASE FLOW

Figure 17. Relation between observed and model-calculated water levels and stream base flow for average 

conditions, 1989–98, for steady-state simulation models of the upper Charles River Basin, eastern Massachusetts: 

(A) Water levels, and (B) Stream base flow. Observed values are estimates for 1989–98 from measurements made in 

1999–2000 as described in the text. Line is line of equality between measured and simulated water level.
Simulation of Ground-Water Flow 49

0

0 1 2 

1 2 MILES

3 KILOMETERS

BASIN BOUNDARY

TILL AND BEDROCK—West and east

 model areas shown

STRATIFIED 

GLACIAL DEPOSITS

SUPPLY WELL OR SURFACE-WATER SOURCE 

WASTEWATER-TREATMENT- 

 FACILITY DISCHARGE

MODEL-CALCULATED STEADY-STATE

 WATER-TABLE CONTOUR—Shows

 altitude of water table in feet above

 sea level. Contour interval 10 feet.

 Datum is sea level.

EXPLANATION

230

W E

270

250

250

250

270

250

250

230

290270

250 310 230

270 290 210

270

210 190

270

250

250

260

200

210

210 190 170

180

230

270 260

250

210

170 190

190

190

210

230 230 270 250

270

210

290 270 210

190

190

210

200 210 210

230 250

230

300

260

280

250

270

250

250

270

290

250

260

310

300

270

270

250 170

230210 210

190 180

190 170 190

170

150

170 150

150

190

180

170

160

150

130

130

130

150

150

130

150

130

150

170

170

210

230 210

250

230

330

310

230

280

220

200 220

220

280

300

320

220

240

260

280

260

260

280

240

240

200

180

220

240

260

140

160

140

160

180

190

220

240

200

140

140

140

160

160

180

200

240

220

260

220

260

240

260

240

220

180

200

280

240

200

160

140

160

320

290

180

280

200

220

WEST 

MODEL 

AREA

EAST 

MODEL 

AREA

200

200

180

180

200

190

71o30'00" 71o

22'30" 42o

15'

00"

42o

07'

30"

From USGS and MassGIS data sources, Geographic Projection, Spheroid GRS 1980, Datum NAD 83 

Figure 18. Model-calculated steady-state water table in the upper Charles River Basin, eastern Massachusetts.
50 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Table 12. Model-calculated steady-state stream base flow and observed stream base flow at measurement sites in the upper 

Charles River Basin, eastern Massachusetts

[Site locations shown in figure 7. Observed base flow: Estimated for 1989–98 from measurements made during 1999–2000, as described in text]

Station

No.

Model location Average annual stream base flow (cubic feet per second)

Layer Row Column Model

calculated

Observed

Difference (model 

calculated minus 

Flow observed)

90-percent confidence limits

Lower Upper 

West Model

Charles River

01103110 1 161 61 1.8 2.0 1.2 3.2 -0.2

01103120 1 197 58 6.9 5.6 3.0 10.5 +1.3

01103140 1 271 76 18.4 15.8 10.6 23.7 +2.6

011032056 1 282 131 27.7 24.2 14.7 27.5 +3.5

01103206 1 239 140 30.5 28.8 17.2 34.0 +1.7

01103260 1 210 186 87.1 86.9 21.4 38.9 +.2

Tributaries to Charles River

01103210 1 233 137 4.1 7.6 0.8 71.7 -3.5

01103217 1 200 136 15.3 14.8 5.8 37.8 +.5

01103225 1 375 167 2.7 2.5 1.5 4.2 +.2

01103234 1 305 192 11.3 13.6 10.2 18.0 -2.3

01103235 1 302 175 12.8 15.0 7.9 28.6 -2.2

01103240 1 232 170 19.0 24.4 14.9 40.0 -5.4

011032515 1 182 176 9.5 8.1 5.5 12.1 +1.4

East Model

Charles River

01103280 1 204 226 95.8 91.0 59.7 139 +4.5

01103305 1 216 263 126 131 103 168 -5.0

Tributaries to Charles River

01103292 1 326 281 6.8 11.3 8.0 16.0 -4.8

01103295 1 279 275 11.9 18.4 11.9 28.4 -6.5

01103300 1 236 258 19.4 26.0 18.3 37.0 -6.6

Bogastow Brook and Tributaries

01103381 1 85 176 3.8 3.3 2.2 4.9 +.5

01103386 1 66 201 1.3 1.0 .5 2.0 +.3

011033885 1 79 232 2.4 1.6 1.1 2.2 +.8

01103389 1 103 223 18.2 12.8 8.0 20.5 +5.4

01103393 1 132 243 29.1 30.2 20.9 43.6 -1.1

The steady-state, average annual hydrologic 

budget calculated with the calibrated model is shown in 

table 13. Total inflows were 70 Mgal/d to the west 

model area and 49 Mgal/d to the east model area. In 

both models, inflows from uplands were large components, representing about 70 percent of inflows to the 

west model and about 50 percent of inflows to the east 

model. These proportions are larger than estimates (30 

to 50 percent) made in the water balance for the study 

area, because upland inflows in the model budget 

include precipitation and return-flow recharge to thin 

areas of stratified glacial deposits that are excluded 

from the active model area. Leakage to the aquifer from 

streams was a small component (less than 10 percent) 

of both model budgets. Stream leakage resulted near 

the model boundary downstream of the specified 

streamflow in the west model; stream leakage occurred 

elsewhere in the models near pumping wells (induced 

infiltration) or downstream of constrictions in aquifer 

geometry. Inflows from wastewater discharge to 
Simulation of Ground-Water Flow 51

1Simulated as enhanced recharge at model boundaries.

2Does not include specified stream inflows to the Charles River, equal to 2.3 ft3/s for the west model, 97.1 ft3/s for the east model, and 2.3 ft3/s for the 

entire basin. These inflows are not part of the model budget calculations.

Table 13. Model-calculated steady-state hydrologic budget for the upper Charles River Basin, eastern Massachusetts

[ft3/s, cubic foot per second]

Hydrologic budget component

Rate of flow

West model East model Entire upper Charles 

Basin

Cubic feet 

per second

Million

gallons

per day

Cubic feet 

per second

Millions 

gallons

per day

Cubic feet 

per second

Million 

gallons

per day

Inflow

Recharge from precipitation and septic-tank 

return flow .............................................................. 19.7 12.7 25.4 16.4 45.1 29.1

Lateral ground-water inflow1 ..................................... 56.1 36.3 22.9 14.8 79.0 51.1

Streamflow from uplands1, 2 ...................................... 22.3 14.4 17.8 11.5 40.1 25.9

Stream leakage to aquifer........................................... 4.7 3.0 3.7 2.4 8.4 5.4

Wastewater discharge to streams................................ 5.5 3.6 6.1 3.9 11.6 7.5

Total inflow2 ........................................................... 108.3 70.0 75.8 49.0 184.1 119.0

Outflow

Ground-water discharge to streams............................ 97.1 62.8 64.5 41.7 161.6 104.5

Evapotranspiration from wetlands and ponds ............ 4.4 3.3 5.1 3.3 9.5 6.1

Water-supply withdrawal............................................ 6.7 4.3 6.2 4.0 12.9 8.3

Total outflow2 ......................................................... 108.2 69.9 75.8 49.0 184.0 118.9

Budget error (inflow-outflow)....................................... 0.1 0.0 0.0 0.0 0.1 0.1

streams also were small components of the model

budgets. Total net inflow to west and east models,

119 Mgal/d, was consistent with the range of values 

determined for the study area with the water balance 

(111 to 141 Mgal/d). 

Calculated outflows from the models consisted 

primarily of discharge to simulated streams; these outflows were 90 and 85 percent of total outflows from 

west and east models, respectively. Although small

relative to total outflows, model-calculated ET from

wetlands and ponds was similar in magnitude (though 

much less well quantified) to withdrawals for water 

supply, indicating the potential importance of this flux. 

Stream outflows that are calculated in the model budget 

represent ground-water discharge to streams within the 

active area. The calculated stream outflows do not 

include inflows specified at model boundaries in the 

stream package; specified inflows are routed through 

the model unaltered unless water in stream cells leaks 

into the aquifer, in which case the inflows are available 

for aquifer recharge. Thus, the simulated stream base 

flow at the exit of the west model was equal to the 

model-calculated flow from ground-water discharge to 

streams (97.1 ft3/s) minus stream leakage to the aquifer 

(4.7 ft3/s) plus the specified inflow for the Charles 

River in Milford (2.3 ft3/s), or 94.7 ft3/s. This value 

was specified as the inflow to the east model for the 

Charles River in Medway. Similarly, the simulated 

stream base flow at the exits of the study area (Charles 

River and Bogastow Brook in Norfolk and Millis)

was equal to the specified inflow to the east model,

plus the model-calculated flow from the ground-water 

discharge to streams for the east model (64.5 ft3/s), 

minus the model-calculated stream leakage to the

aquifer in the east model (3.7 ft3/s), or 155.5 ft3/s 

(100.5 Mgal/d). The simulated stream base flow at 

model exits plus other model outflows of evapotranspiration and water-supply-withdrawals yields total net 

outflows that are consistent with the values determined 

for the study-area water balance (119 Mgal/d with

a 90-percent confidence-interval range of 95 to

152 Mgal/d on the basis of base-flow estimates). 
52 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Withdrawals for water supply account for about 

7 percent of model-calculated average annual inflows 

or outflows on a basin-wide basis (table 13). Within 

individual tributary basins, however, water-supply 

withdrawals constituted more or less than this basinwide average. For example, water-supply withdrawals 

accounted for 13 percent of average annual inflows in 

the headwaters Charles River area in Milford (upstream 

of station 01103140), about 10 percent of inflows in the 

Stall Brook, Mine Brook, and Mill River Subbasins 

(areas upstream of stations 01103210, 01103240, and 

01103300, respectively), about 8 percent of inflows in 

the Bogastow and Dopping Brook Subbasins (upstream 

of station 01103989), and about 1.5 percent of inflows 

in the Hopping Brook Subbasin (upstream of station 

01103217). 

Transient Numerical 

Models

Transient models were developed to simulate the 

variations in hydrologic conditions within an average 

annual cycle. The transient models are based on the 

steady-state models but incorporate time-varying 

hydraulic stresses and boundary conditions. The spatial 

discretization of the model grid, boundary conditions 

other than specified flows, and spatial variations in 

stresses and hydraulic conductivities are the same in 

transient and steady-state models. The transient models 

are designed to simulate dynamic equilibrium, or the 

condition in which there is no net change in storage 

over the annual cycle (Barlow and Dickerman, 2001). 

The transient models were calibrated by comparing 

simulated stream base flow and water levels to average 

monthly flows and levels estimated for the 1989–98 

period.

With the transient models, the low-flow periods 

of the annual cycle can be simulated. Low-flow

periods, typically the late summer months, often are of 

particular concern in the evaluation of the effects of 

water-management alternatives. During these periods, 

the effects of water withdrawals and other management 

practices on aquatic life and stream-water quality often 

are greatest, because their effects are combined with 

naturally low flows and ground-water levels. Water 

demands also typically are highest during summer 

months. Thus, seasonally varying stresses and fluxes 

within the annual cycle often must be considered when 

the effects of water-management alternatives are tested.

Temporal Discretization and 

Initial Conditions

The annual hydrologic cycle was divided into 12 

monthly time periods that varied in length from 28 to 

31 days. Within each month, hydraulic stresses and 

boundary flows were uniform. Multiple model stress 

periods (with the same hydraulic stresses and boundary 

flows) were used to simulate some months to improve 

model stability, particularly during months of low 

recharge. Fifteen stress periods per year were used in 

the west model and 19 to 23 stress periods per year 

were used in the east model. Stress periods were 

divided into 30 time steps, regardless of length. The 

time steps increased in length by a factor of 1.3 during 

each stress period, for example, from 4.6 x10-4 to

0.9 day for stress periods of 4 days and from 3.6 x10-3

to 7.2 days for stress periods of 31 days. 

Water-level elevations from the calibrated 

steady-state models were specified as the initial conditions for the transient simulations. The calibrated 

steady-state water levels represent average annual

conditions. The transient simulations also began with 

stresses and boundary flows representing conditions

in the month of January. Because January conditions 

differ from average annual conditions, changes in

calculated water levels and flows in early parts of the 

transient simulations resulted from this initial discrepancy as well as from time-varying stresses. Transient 

simulations were run for repeated 1-year cycles until 

the effects of the initial conditions were eliminated and 

there was no change in storage over a 1-year cycle. 

After five annual cycles for both west and east models, 

the difference between flow into and out of storage was 

0.2 and 0.1 percent for west and east models, respectively. These difference indicated that net changes in 

storage were negligible, compared to the total water 

budget. Five annual cycles consisted of 75 stress periods for the west model and 107 stress periods for the 

east model.

Boundary Conditions and 

Stresses

Boundary conditions in the transient models 

were similar to those used in the steady-state models. 

As in the steady-state models, stream base flow was 

specified in the stream package at locations where the 

Charles River entered the active area of west and east 

models. For the Charles River at the upstream boundary of the west model, monthly average base flow was 
Simulation of Ground-Water Flow 53

varied on the basis of the monthly distribution of 

upland inflows, described below, and the base flow 

estimated from measurements at streamflow site 

01103110 (fig. 7). Monthly average base flow in the 

Charles River at the upstream boundary of the east 

model was equal to the model-calculated monthly average base flow in the Charles River at the west model 

exit. Stresses in the transient models were of the same 

type as in the steady-state models, but varied monthly. 

Average monthly recharge rates were based on 

the average annual recharge rate from the calibrated 

steady-state models, analyses of the monthly distribution of annual recharge in the Hunt River Basin in 

southern Rhode Island (Barlow and Dickerman, 2001), 

and observed streamflows at two sites near active 

model boundaries. Temporal recharge patterns in the 

Hunt River Basin were considered representative for 

the upper Charles River Basin because climatic conditions, average annual recharge rates, and the areal

percentage of stratified glacial deposits in the two 

basins are similar. The monthly distribution of annual 

recharge in the Hunt River Basin was determined

by Barlow and Dickerman (2001) from the analysis

of long-term streamflow records with a base-flow 

recession-curve displacement method developed by 

Rorabaugh (1964) and automated by Rutledge (1993). 

This approach to estimating average monthly varying 

recharge rates has the advantage over water-balance 

methods applied at the monthly scale. Because it is 

implemented with a daily time step, the Rutledge 

(1993) method can yield greater-than-zero recharge 

rates in summer months, whereas comparison of 

monthly average potential ET and monthly average 

precipitation would yield monthly average recharge 

estimates equal to zero. Although the recession-curve 

displacement method may under- or over-estimate 

recharge rates in months when ground-water ET, 

underflow, or water withdrawals (assumed negligible in 

this approach) are substantial outflows from the 

stream-aquifer system (Nicholson and Watt, 1997; R.S. 

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

2001), the method was applicable to the Hunt River 

Basin because these withdrawals were small relative

to streamflow (Barlow and Dickerman, 2001).

Average monthly recharge rates were calculated 

by varying the calibrated rate for natural precipitation 

recharge for stratified glacial aquifers proportionately 

through the year on the basis of the distribution of 

annual recharge among months in the Hunt River 

Basin. This approach made rates for all recharge categories uniformly proportional, because the cell-by-cell 

recharge rates applied in the model were the product of 

the natural precipitation rate to stratified glacial deposits and a multiplier array. Initially, the same distribution 

of annual recharge among months was used for 

recharge to interior cells of the active model area, 

which represents only precipitation and return-flow 

recharge to stratified glacial aquifers, and for enhanced 

recharge to cells along the active model boundary, 

which represents inflows from uplands as well as precipitation and return-flow recharge. This initial uniform 

monthly distribution of recharge was modified during 

model calibration. Final monthly recharge rates to interior cells (those representing stratified glacial deposits 

without inflow from uplands) and flow rates to cells 

receiving upland inflows are shown in figure 19A. The 

monthly rates for upland inflows and for recharge to 

stratified glacial deposits both correspond to the same 

average annual rate as used in the steady-state models. 

Flow rates to cells receiving upland inflows are more 

variable from month to month than recharge rates to 

interior aquifer cells. The variability in upland flow 

rates is comparable to the variability observed in estimated monthly base flow for streams draining upland 

basins in the upper Charles River Basin (stations 

01103110 and 011033885; fig. 19B).

Wetlands and ponds drained by streams were 

simulated in the steady-state models as areas of net 

specified water loss from the stream-aquifer system. 

The rate of loss was equal to the growing-season, freewater-surface evaporation rate, or 21 in/yr (Farnsworth 

and others, 1982). In the transient simulations, this loss 

was assumed to occur uniformly throughout the 

months of May to October, for an average monthly rate 

of 3.5 in/month. Ponds not drained by streams were 

simulated as net recharge areas (22 in/yr) in the steadystate models. This recharge was assumed to be distributed between summer and winter months by removing 

free-water-surface evaporation rates from precipitation 

during these periods (Farnsworth and others, 1982). 

Recharge to ponds not drained by streams was specified as 2.7 in/month from November to April and

0.9 in/month from May to October in the transient 

models.

Monthly withdrawals for water supply and discharges from wastewater-treatment facilities were set 

equal to average monthly volumes for 1989–98. Withdrawals at water supplies outside the active area of the 

east model, the town of Franklin supply well No. 9 and 

the Franklin Country Club well and reservoir (fig. 3), 
54 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

were not varied monthly. The average annual 

withdrawals were used for these sources. In steadystate and transient models, withdrawals at these sources 

were represented by reductions in the upland inflows, 

which were input to the models by use of a recharge 

multiplier. Echo Lake Reservoir (fig. 3), in the west 

model, also was outside of the active model area. Withdrawals at Echo Lake Reservoir were represented, as

in the steady-state model, by excluding the drainage 

area of Echo Lake Reservoir from the upland area contributing water to the active model area and by including winter withdrawals of water for transfer to Echo 

Lake Reservoir from the Charles River surface-water 

source (fig. 3). Water withdrawals for transfer from the 

Charles River source averaged 0.57 Mgal/d in January, 

1.30 Mgal/d in February, 1.78 Mgal/d in March,

1.56 Mgal/d in April, and 0 Mgal/d in other months, 

based on available data.

Hydraulic Properties

A uniform value of 0.28 for specific yield was 

used for stratified glacial deposits in most active model 

areas. This value is similar to values of specific yield 

determined for stratified glacial deposits on Cape Cod 

and in southern Rhode Island (Allen and others, 1963; 

Moench and others, 2000) and to values determined 

from laboratory studies (Morris and Johnson, 1967; 

Johnson, 1967; Nwankwor and others, 1984). For areas 

along the Mill River in the east model, where available 

data indicated locally confined conditions, specific 

yield values of 0.2 and 0.15 were used. These values 

were used to represent spatially averaged effects of the 

local confining units, which otherwise were not simulated. A specific yield of 1.0 (100 percent porosity) was 

specified for all simulated ponds. A uniform value for 

specific storage of 1.0 x 10-3 per foot was specified for 

stratified glacial deposits in most areas. This value corresponds to a storage coefficient of 1.0 x 10-2 for a 

10-ft saturated thickness of aquifer, which is within

the range of values determined for stratified glacial 

deposits in New England (Kontis and others, in press). 

Model Calibration

The transient models were calibrated to waterlevel fluctuations measured at observation wells

and ponds during the study (August 1999 through 

November 2000) and estimated average (1989–98) 

monthly base flow at streamflow-measurement sites 

(“observed values”). Trial-and-error adjustment of

storage parameters and, to some extent, the monthly 

distribution of annual recharge, was used to obtain reasonable agreement between simulated and observed 

values for both water levels and flows. Calibration of 

west and east models proceeded concurrently. Particular weight was given to the distribution of monthly 

streamflow at measurement sites near model exits on 

the Charles River (stations 01103260 and 01103305) 

and Bogastow Brook (01103393; fig. 7), because these 

sites integrate conditions throughout the model areas. 

Greater consideration also was given to simulated 

water-level fluctuations at observation wells near the 

horizontal center of aquifers than to simulated fluctuations at wells near model boundaries with inactive 

upland areas. 

TILL AND BEDROCK

AVERAGE MONTHLY 

RECHARGE, IN INCHES

PER MONTH

NORMALIZED AVERAGE

 MONTHLY FLOW, DIMENSIONSLESS

5

-0.5

0.0

0.5

1.0

1.5

0

1

2

3

4

-1.0

J FM MJ J A ASOND

STRATIFIED GLACIAL 

DEPOSITS

UPLAND RECHARGE

BASE FLOW AT 01103110

BASE FLOW AT 011033885

A.

B.

Figure 19. Monthly average recharge rates for transient 

simulation models and comparison of monthly average 

rates for flows to the active model area from upland till and 

bedrock areas with estimated monthly base flow in two 

streams draining upland watersheds in the upper Charles 

River Basin, eastern Massachusetts: (A) Monthly average 

recharge rates for stratified glacial deposits and rates for 

flows to the active model area from upland till and bedrock 

areas. Monthly average recharge rates for stratified glacial 

deposits are from Barlow and Dickerman (2001), and

(B) Monthly average rates for flows to the active model 

area from upland till and bedrock areas and estimated 

monthly average base flow in the Charles River above 

Cedar Swamp Pond near Milford (01103110) and Dirty 

Meadow Brook at Hollis Street near Holliston (011033885). 

Rates shown are normalized to average annual rates.
Simulation of Ground-Water Flow 55

Simulated and observed water-level fluctuations 

at selected observation wells and ponds are shown in 

figure 20. Fluctuations are shown relative to the average of model-calculated monthly water-level elevations 

for model-calculated values and relative to the average 

of water levels measured in water year 2000 for the 

observed values (table 5). In general, seasonal patterns 

of measured water-level fluctuations were reasonably 

well matched by simulated water levels. The timing of 

model-calculated and measured maximum and minimum water levels within the annual cycle was similar 

at most sites, although simulated water levels peaked 

1–2 months before measured water levels at some sites. 

Model-calculated water levels also tended to fluctuate 

more during the annual cycle than water levels measured during the study. The average difference between 

model-calculated and observed annual water-level

fluctuations was 0.9 ft (median 0.6 ft) for the 48 measurement sites (table 5) in the active model area. Large 

differences (greater than about 2 ft) generally occurred 

at observation wells adjacent to model boundaries (for 

example, F2W72, fig. 20A). Large differences also 

occurred at some ponds on the stream network (for 

example, Box Pond, fig. 20B) and may have resulted 

from the inability of the model to simulate streamflow 

hydraulics, such as pond storage. Differences between 

observed and simulated fluctuations also may have 

resulted from deviations of observed conditions during 

the measurement period from long-term average

conditions. 

Simulated and observed average monthly base 

flow at selected measurement sites on the Charles 

River, Bogastow Brook, and tributaries are shown in 

figure 19. Simulated and observed average monthly 

base flow agreed reasonably well at most sites.

Most model-calculated monthly base flow was within 

the range of values represented by 90-percent confidence intervals for observed base-flow estimates. Temporal patterns in base flow were simulated and the 

mean absolute base-flow residuals averaged 41 percent 

of observed flows overall. At sites on the Charles River, 

observed monthly flows were better matched, with the 

mean absolute base-flow residuals averaging 26 percent of observed values overall at eight measurement 

sites (table 12 and fig. 21A–D). The mean absolute 

base-flow residuals averaged 22 percent of observed 

values at each of the model exit sites on the Charles 

River (01103260 and 01103305) and 50 percent of 

observed values at the model exit site on Bogastow 

Brook (01103393), respectively. At most sites, including sites on the Charles River and some tributaries (for 

example, Chicken and Dopping Brook sites, fig. 21H

and J), model-calculated spring high base flow was 

larger than observed flows (fig. 22A; mean base-flow 

residual for March flows averaged 11 percent of 

observed values; mean absolute base-flow residuals 

averaged 29 percent of observed values). Modelcalculated late-summer low flows were less biased than 

high flows (mean base-flow residuals for September 

flows averaged 5 percent of observed values), but there 

was more scatter in observed and model-calculated 

values at the low flows (fig. 22B). Late-summer low 

flows occurred later in the annual cycle for simulated 

flows than for observed flows at most of these sites

(fig. 21). At a number of sites on tributaries draining 

proportionately large upland areas (for example,

Hopping Brook, 01103217, fig. 21F) or on tributaries 

in the Mine Brook aquifer area, spring flows were 

underestimated and summer flows overestimated.

The flashiness of streams in these areas may have been 

more strongly affected by lateral and channeled inflow 

from uplands than streams in other aquifer areas. 

Finally, discrepancies between model-calculated and 

observed base flow from the calibrated steady-state 

model resulted in large differences between modelcalculated and observed monthly average base flow for 

some sites in the transient models, including those in 

the Mill River aquifer area (for example, 01103295, 

fig. 21I). Generally, seasonal patterns in flows are

reasonably well matched at these sites, however.

During model calibration, aquifer storage properties and the distribution of annual recharge among 

months were varied to minimize the differences 

between simulated and measured water-level fluctuations and stream base flow. Specific yield was varied 

between values of 0.26 and 0.30 (except in zones where 

locally confined conditions were present). Although a 

specific yield of 0.30 improved the east model fit 

slightly with respect to water-level fluctuations, a value 

of 0.28 was retained for both models because it was 

considered more representative for aquifer sediments in 

the study area, on the basis of literature sources, and 

because of greater numerical instability in the west 

model with the higher value of specific yield. Specific 

storage also was varied, but had little effect on waterlevel fluctuations or stream base flow. The distribution 

of inflows from upland areas, simulated as enhanced 

recharge to cells along the active model boundary, was 

varied as described previously to increase the annual 

range of simulated base flow and better match the 

observed annual values. 
56 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

NNW108

5

4

3

2

0

1

-5

-4

-3

-2

-1

5

4

3

2

0

1

-5

-4

-3

-2

-1

5

4

3

2

0

1

-5

-4

-3

-2

-1

5

4

3

2

0

1

-5

-4

-3

-2

-1

A JJ SOND J FM MA ASOND A JJ SONDJ FM MA ASOND A JJ SOND JFM MA ASOND

WATER-LEVEL FLUCTUATION, IN FEET

F2W67

F2W72 HTW47

HVW40

HTW51

MNW20 MYW58

NNW108 NNW109 XUW67

A6W53 A6W60

MODELCALCULATED

OBSERVED

A.

1999 2000 1999 2000 1999 2000

Figure 20. Model-calculated and observed water-level fluctuations during the annual average cycle for selected 

observation wells and ponds in the upper Charles River Basin, eastern Massachusetts: (A) Observation wells, and 

(B) Ponds. Fluctuations are shown relative to the average of water levels measured in water year 2000 for 

the measured values (see table 5) and relative to the annual average of calculated monthly water levels for 

model-calculated values (see table 11).
Simulation of Ground-Water Flow 57

5

4

3

2

0

1

-5

-4

-3

-2

-1

A JJ SOND J FM MA ASOND

WATER-LEVEL FLUCTUATION, IN FEET

BOX POND POPULATIC POND

BUSH POND LAKE PEARL

A JJ SOND J FM MA ASOND

A JJ SOND J FM MA ASOND 5

4

3

2

0

1

-5

-4

-3

-2

-1 MODEL-CALCULATED

OBSERVED

KINGSBURY POND

B.

1999 2000 1999 2000

1999 2000

EXPLANATION

Figure 20. Model-calculated and observed water-level fluctuations during the annual average cycle for selected 

observation wells and ponds in the upper Charles River Basin, eastern Massachusetts: (A) Observation wells, and 

(B) Ponds. Fluctuations are shown relative to the average of water levels measured in water year 2000 for 

the measured values (see table 5) and relative to the annual average of calculated monthly water levels for 

model-calculated values (see table 11)—Continued.

Differences between simulated and observed 

water-level fluctuations and stream base flow may have 

resulted from various causes, including model calibration error, discretization effects, or inadequate simulation of aquifer geometry, storage properties, recharge or 

other hydrologic processes. Deviations between the 

annual cycle of water levels, as represented by modelcalculated, long-term monthly averages and as represented by the instantaneous measurements made during 

the study period also may have contributed to the differences between simulated and measured water-level 

fluctuations. The observed long-term monthly average 

streamflows, which were estimates from partial streamflow records, also may have contained multiple sources 

of error. Further modifications to aquifer-storage

properties, to specified rates of precipitation and returnflow recharge, or to specified rates of enhanced 

recharge representing inflows from uplands might 

have been made to improve the match between simulated and observed water-level fluctuations and streamflows. However, these modifications were judged to

be inappropriate given the limited availability of data 

for these variables and the expected sources of error 

inherent in the observations against which simulated 

water-level fluctuations and stream base flow were 

compared.

The average annual hydrologic budget for the 

entire upper Charles River Basin calculated with the 

transient models agreed well with the steady-state

budget. Inflow and outflow components, except for 

stream leakage to the aquifer, differed by 1 percent or 

less between transient and steady-state models; thus, the 

flows specified in the steady-state models were accurately represented in the transient models. Larger percentage differences between the steady-state and 

transient models in stream leakage, which constitutes 

about 5 percent of total flow, were probably insignificant. 
58 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

C. 01103305 D. 01103260

E. 01103393

J FM MJ J A ASOND

MODEL-CALCULATED

100

80

60

40

20

0

300

200

100

0

A. 01103140 B. 01103206

100

80

60

40

20

0

J FM MJ J A ASOND

STREAM BASE FLOW, IN CUBIC FEET PER SECOND

OBSERVED—Vertical bar

 shows the 90-percent 

 confidence interval

EXPLANATION

Figure 21. Model-calculated and observed average monthly base flow in the Charles River, Bogastow 

Brook, and tributaries, upper Charles River Basin, eastern Massachusetts: (A) Charles River at Factory 

Pond at South Milford (01103140), (B) Charles River at North Bellingham (01103206), (C) Charles 

River at West Medway (01103260), (D) Charles River near Millis (01103305), (E) Bogastow Brook 

below Great Black Swamp near Millis (01103393), (F) Hopping Brook near West Medway (01103217), 

(G) Mine Brook near Franklin (01103235), (H) Chicken Brook below Milk Pond near West Medway 

(011032515), (I) Mill River below Bush Pond near City Mills (01103295), and (J) Dopping Brook at 

Whitney Street, Holliston (01103386). Observed values are estimates for 1989–98 based on measured 

flows in 1999–2000, as described in text.
Simulation of Ground-Water Flow 59

J FM MJ J A ASOND

STREAM BASE FLOW, IN CUBIC FEET PER SECOND

J FM MJ J A ASOND

30

50

60

40

20

0

10

30

50

60

40

20

0

10

F. 01103217 G. 01103235

H. 011032515 I. 01103295

10

20

15

5

0

J. 01103386

62 76 106 103 60

MODEL-CALCULATED

OBSERVED—Vertical bar

 shows the 90-percent 

 confidence interval

EXPLANATION

Figure 21. Model-calculated and observed average monthly base flow in the Charles River, Bogastow 

Brook, and tributaries, upper Charles River Basin, eastern Massachusetts: (A) Charles River at Factory 

Pond at South Milford (01103140), (B) Charles River at North Bellingham (01103206), (C) Charles 

River at West Medway (01103260), (D) Charles River near Millis (01103305), (E) Bogastow Brook 

below Great Black Swamp near Millis (01103393), (F) Hopping Brook near West Medway (01103217), 

(G) Mine Brook near Franklin (01103235), (H) Chicken Brook below Milk Pond near West Medway 

(011032515), (I) Mill River below Bush Pond near City Mills (01103295), and (J) Dopping Brook at 

Whitney Street, Holliston (01103386). Observed values are estimates for 1989–98 based on measured 

flows in 1999–2000, as described in text—Continued.
60 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA MODEL -CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND

OBSERVED BASE FLOW,

IN CUBIC FEET PER SECOND

20

40

50

0

30

10

100

200

250

0

150

50

0 150 50 100 200 250 0 30 10 20 40 50 MODEL -CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND

OBSERVED BASE FLOW,

IN CUBIC FEET PER SECOND

A. MARCH B. SEPTEMBER

WEST MODEL

EAST MODEL

LINE OF EQUALITY

EXPLANATION

Figure 22. Relation between observed and model-calculated average March and September base flow for 

measurement sites in the upper Charles River Basin, eastern Massachusetts: (A) March, and (B) September. Observed 

values are estimates for 1989–98 from measurements made in 1999–2000 as described in the text. Line is line of 

equality between observed (estimated) and model-calculated values.

The average annual rate of inflow to and outflow from 

aquifer storage calculated with the transient models was 

21 Mgal/d, or about 18 percent of total flow through the 

stream-aquifer system. Water is added to aquifer storage 

from November through April and removed from storage 

from May through October. Hydrologic budgets calculated with the transient models also indicated that 

monthly water withdrawals, which averaged 7 percent of 

total flow on an annual average basis, range from about 4 

to 5 percent of model-calculated total monthly flow in 

winter and spring months (November to May) to about 

10 to 12 percent of model-calculated total monthly flow 

in summer and early fall months (June to October).

Model Limitations

The steady-state and transient flow models of the 

upper Charles River Basin stream-aquifer system provide a regional-scale simulation of ground-water flow in 

the stratified glacial aquifers in the study area. As with 

all mathematical models of natural systems, the simplifications and assumptions incorporated into the models 

result in limitations to their appropriate uses and to the 

interpretations that may be made of simulation results.

The upper Charles River Basin flow models simulate ground-water flow, water levels, and the interaction 

with surface-water features at the regional scale. Hydrologic processes and spatial variability in hydraulic properties and stresses are simplified and approximated to a 

degree consistent with this scale. The model calibration 

also represents the best fit to estimates and observations 

made throughout the upper Charles River study area. 

Thus, the agreement between simulated water levels or 

stream base flow in specific areas of the flow system 

may not be adequate to support local-scale model applications. The model discretization also may not be applicable for local-scale problems. In fact, model calibration 

generally should be taken into consideration whenever 

the effects of simulated changes in water-management 

practices or hydrologic conditions at specific observation 

points are evaluated. Because of the error associated 

with regional-scale calibration, simulation results for 

many applications may best be viewed as evidence of 

relative changes in water levels or base flow rather than 

as absolute changes at specific locations.

The effects of temporal and spatial discretization 

also impose limitations on model use. Hydrologic processes and hydraulic stresses were represented in the 

transient models as monthly averages. Simulation results 

are monthly average ground-water levels and flows. The 

models were not intended to be used to simulate changes 

that occur at finer time scales, for example, daily values, 

which may substantially exceed or fall below monthly 

average values. The spatial resolution of the simulation 

results was limited by the area of the 200x200-ft grid 

cell. Water withdrawals, discharges, and streamflow and 
Evaluation of Ground-Water-Management Alternatives 61

water-level observations were averaged within grid cells 

and their exact locations were approximated by the centers of the cells in which they occur. 

The upper Charles River Basin models described 

in this report simulate ground-water flow and water levels. Flows and water levels in surface-water features are 

simulated to the extent that they represent discharge 

areas or boundaries for the ground-water system. Thus, 

flows in simulated streams are base flow (ground-water 

discharges) and do not include the direct runoff component of streamflow. Moreover, the hydraulics of the

surface-water system, such as storage provided by 

impoundments and wetlands, were not simulated. The 

effects of these controls on surface-water flow were not 

included in the simulated stream base flows. Finally, the 

approach of representing stream stage by a fixed value 

representing average conditions may lead to some small 

inaccuracies in flow rates between aquifers and streams, 

particularly during periods of high flow.

Ground-water flow through till and bedrock are 

not directly simulated in the models. Ground water in 

fractured bedrock can have a widely variable area of 

recharge and natural discharge. Thus, although water 

withdrawals from bedrock aquifers may be simulated 

indirectly (as reductions in recharge), the effects of these 

withdrawals on the stream-aquifer system may not be 

appropriately addressed with the models. Assumptions 

about ground-water flow and recharge from till and bedrock to stratified glacial aquifers and streams also may 

place limitations on interpretation of simulation results. 

In the absence of any information about flow rates or 

pathways through bedrock, recharge in till and bedrock 

upland areas was routed directly to the lateral boundaries 

of the stratified glacial deposits, and no flow was 

assumed to occur between stratified glacial deposits and 

the underlying bedrock. Where flow through bedrock is 

substantial, the spatial or temporal distribution of simulated flows may not be adequately represented.

Finally, the models were calibrated to average 

stream base flow and water levels that were estimated 

from short-term monitoring data and correlations with 

long-term monitoring records from stations outside

the basin. Errors in estimated averages may have

contributed to the discrepancies between modelcalculated water levels and stream base flow and

actual conditions.

EVALUATION OF GROUNDWATER-MANAGEMENT

ALTERNATIVES

The ground-water-flow models were developed 

as a tool to evaluate the response of the upper Charles 

River stream-aquifer system to changes in watermanagement practices or hydrologic conditions. 

Increased water withdrawals, alternative pumping 

schedules for existing withdrawals, additional surfacewater discharges, land disposal of treated wastewater, 

sewering, or stormwater recharge are examples of 

water-management practices that may be tested with 

the models. Altered hydrologic conditions that could be 

simulated include drought conditions or conditions of 

altered recharge that may result, for example, from 

land-use changes. 

Two approaches were used to investigate alternative ground-water management practices and altered 

hydrologic conditions in the upper Charles River Basin. 

First, the flow models alone were used to determine the 

effects of increased withdrawals and altered recharge in 

several hypothetical scenarios. The scenarios were 

selected to represent possible future changes in water 

use in the basin, or to investigate the effects of watermanagement practices that could mitigate potential 

adverse effects of increased water withdrawals. Second, the flow models were used in conjunction with 

optimization techniques. This approach was used to 

quantify possible increases in withdrawals or streamflow that could be obtained by managing water withdrawals. The simulation-optimization approach is 

advantageous for many applications, because it determines the combination of practices to best (optimally) 

meet management goals for a particular set of constraints. It was applied in a sub-area of the model to 

demonstrate its use in addressing complex watermanagement questions in the valley-fill aquifer hydrogeologic setting. Results of both approaches are 

described in terms of changes in average monthly 

stream base flow and pond levels from the base flow 

and pond levels calculated with the calibrated transient 

models.
62 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Simulation of Increased 

Water Withdrawals and 

Altered Recharge

Hypothetical scenarios of increased water withdrawals and altered recharge in the upper Charles River 

Basin were identified for testing through consultation 

with the Technical Advisory Committee, MADEM, 

MADEP, CRWA, and town officials. The scenarios 

were (1) increasing water withdrawals from existing 

sources in the study area to levels currently permitted 

under the WMA; (2) increasing water withdrawals 

from increasing and proposed sources in the study area 

to levels currently permitted or proposed for permitting 

under the WMA; and (3) altering recharge in the Mine 

Brook aquifer area in Franklin to represent hypothetical changes in sewering or management of residential 

rooftop runoff.

Simulation of Increased 

Withdrawals

Two scenarios of increased withdrawals were 

tested with the flow models. In the first scenario (S1), 

water withdrawals at existing sources in the study area 

(table 2) were increased such that the total pumping 

rates for each town or non-municipal source were equal 

to the rates currently permitted under the WMA. In the 

second scenario (S2), water withdrawals in the study 

area were increased to the currently permitted rates 

under the WMA or, in cases where increases have been 

requested, to proposed levels. Withdrawals are made 

from both existing and proposed sources in S2. Proposed sources were identified as those in the process of 

being permitted by MADEP at the time of this study. In 

both S1 and S2, inflows from precipitation recharge 

remained the same and corresponded to average conditions during the model calibration period (1989–98).

WMA-permitted withdrawal rates for water systems in the upper Charles River Basin are shown in 

figure 23 (“WMA-permitted” in this report refers to 

withdrawals that are both registered and permitted 

under the WMA). Combined withdrawals from all 

sources in a water system within a river basin are permitted under the WMA as system-wide average annual 

rates (Duane LeVangie, Massachusetts Department of 

Environmental Protection, oral commun., 2000). Withdrawals from individual sources (“Zone II” approved 

rates) are separately permitted by MADEP (B.R. 

Bouck, Massachusetts Department of Environmental 

Protection, oral commun., 2000). For most towns and 

all non-municipal systems, system-wide average 

annual rates during 1989–98 were less than WMApermitted rates (fig. 23). System-wide average annual 

rates during 1989–98 met or slightly exceeded currently permitted WMA rates in Holliston and Medway; 

simulated withdrawals for sources in these towns were 

set to 1989–98 rates in S1 and to proposed rates in S2. 

Simulated withdrawals for sources in all other water 

systems were increased so that system-wide average 

annual rates equalled currently permitted rates for S1 

and equalled all currently permitted plus proposed rates 

for S2. For Norfolk, which is only partly in the study 

01 5 2 6 3 4

Bellingham

Franklin

Holliston

Medway

Milford

Norfolk

Wrentham

Franklin CC

Glen Ellen CC

Maplegate CC

NEA

CURRENT WMA-PERMITTED RATES

CURRENT PLUS PROPOSED WMAPERMITTED RATES

1989–98 WITHDRAWAL RATES

MILLIONS OF GALLONS

PER DAY

EXPLANATION

Figure 23. Permitted and registered withdrawal rates 

under the Massachusetts Water Management (WMA) 

and 1989–98 average withdrawal rates for municipal 

and large non-municipal water systems in the upper 

Charles River Basin, eastern Massachusetts. Rates are 

system-wide annual averages. (CC, country club; NEA, 

Northeast Energy Association.)
Evaluation of Ground-Water-Management Alternatives 63

area but entirely within the Charles River Basin, 

allowable increases under the WMA for S1 and S2 

were prorated. S1 and S2 withdrawals for Norfolk were 

increased to levels reflecting the contribution of 

sources in the study area to total Norfolk withdrawals 

during 1989–98. For Millis, with all its water-supply 

sources outside of the study area, neither water

withdrawals nor return flows were increased.

Simulated system-wide increases in average 

annual total withdrawal rates were distributed among 

months of the annual cycle on the basis of the monthly 

distribution of 1989–98 withdrawals for each system. 

For S1, increases were distributed among individual 

sources in each system in proportion to the historical 

(1989–98) contribution of each source to monthly total 

withdrawals in the system. Because withdrawals from a 

new source in Bellingham (well No. 12) began in mid1997, increased withdrawals were distributed among 

sources in that town on the basis of 1998 data. For S2, 

increased withdrawals were preferentially assigned to 

new sources. Average monthly withdrawals at municipal sources were not allowed to exceed their Zone II 

approved rates in either scenario. If distribution of 

system-wide increases according to the procedures 

described would have resulted in its Zone II permitted 

rate being exceeded at a source during any month, 

withdrawals at the source were set to the Zone II rate 

during the month and the excess withdrawal was distributed among the remaining sources. For Milford, the 

currently WMA-permitted withdrawals could not be 

met with existing sources without increasing surfacewater withdrawals above historical levels or exceeding 

Zone II permitted rates at ground-water sources. Consequently, withdrawals for Milford for S1 and S2

were set to Zone II permitted rates in all months for 

ground-water sources and to 1989–98 withdrawals for 

surface-water withdrawals. The monthly distribution

of simulated municipal withdrawals and discharges, 

totaled by town or facility, for S1 and S2 are shown

in figure 24; average withdrawal rates by source are 

listed in table 10. Overall, total simulated withdrawals 

increased in S1 by 3.2 Mgal/d (4.9 ft3/s) to 7.7 Mgal/d 

(11.9 ft3/s) in the west model and by 0.92 Mgal/d

(1.42 ft3/s) to 5.1 Mgal/d in the east model. In S2,

total simulated withdrawals increased by 2.9 Mgal/d 

(4.5 ft3/s) to 7.4 Mgal/d in the west model and

2.4 Mgal/d (3.7 ft3/s) to 6.6 Mgal/d in the east model. 

Withdrawals in the west model were slightly less in S2 

than in S1 because, for some towns, pumping at proposed wells located in the east model area in S2 

resulted in slightly lower pumping rates at existing 

wells in the west model area. Thus, total withdrawals 

simulated in S1 and S2 increased by 42 and 54 percent, 

respectively, over 1989–98 withdrawals.

Return flow from increased withdrawals for 

municipal water systems was routed to existing

wastewater-treatment plants or directly to the aquifer

as septic-system return flow, according to 1989–98 

wastewater-management practices in towns in the study 

area. Simulated additional return flow equaled the 

additional average monthly withdrawals per town, 

reduced by a factor representing consumptive use. 

Consumptive use, which averages about 10 percent 

annually, was varied from 0 percent in winter months 

to about 30 percent in June and July. For towns with 

sewers, additional return flow was added to simulated 

monthly discharges from the MTF (Milford) and the 

CRCPD treatment facility (Bellingham, Franklin, and 

Medway). For towns without sewers (Holliston, Norfolk, and Wrentham), additional septic-system return 

flow was added, as areal recharge. Recharge from additional septic-system return flow was applied in the 

same month as the additional withdrawal. For each 

town, the additional recharge was distributed between 

upland and active model areas in proportion to land 

areas. Additional return flow was assumed to be distributed throughout the total town areas for Norfolk 

and Wrentham. Additional monthly return-flow rates 

averaged 0.49 in/yr for Holliston (S2), 0.47 in/yr (S1) 

and 0.86 in/yr (S2) for Norfolk, and 0.14 in/yr (S1 and 

S2) for Wrentham. Return flow was not simulated for 

non-municipal sources. As in the calibrated models of 

current conditions, withdrawals for golf courses were 

reduced by 50 percent to account for return flow from 

irrigation use, and withdrawals for the Northeast 

Energy Association sources were considered to be 100 

percent consumptive.

Model-calculated stream base flow under scenarios of increased withdrawals (S1 and S2) is compared 

with model-calculated base flow for average 1989–98 

conditions in table 14 and figures 25–27. Comparisons 

are shown of average monthly flows for March and 

September. The magnitude of reductions in base flow 

relative to 1989–98 flows was variable throughout

the annual cycle, and was strongly affected by the 

timing of increased withdrawals at upstream sources. 
64 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

CRPCD TREATMENT 

FACILITY

MEDWAY

NORFOLK WRENTHAM HOLLISTON

FRANKLIN MILFORD

BELLINGHAM

3

2

1

0

6

4

2

0

5

3

1

6

4

2

0

5

3

1

3

2

1

0

3

2

1

0

3

2

1

0

3

2

1

0

9

4

2

0

6

3

1

7

5

8

9

4

2

0

6

3

1

7

5

8

SCENARIO 2

1989–98

SCENARIO 1

AVERAGE MONTHLY WITHDRAWAL OR DISCHARGE, IN MILLIONS OF GALLONS PER DAY

J FM MJ J A ASOND J FM MJ J A ASOND J FM MJ J A ASOND

MILFORD TREATMENT

FACILITY

EXPLANATION

Figure 24. Average monthly withdrawals and discharges by town or treatment facility for 1989–98 and two 

hypothetical scenarios of increased withdrawals in the upper Charles River Basin, eastern Massachusetts. 

(CRPCD, Charles River Pollution Control District.)
Evaluation of Ground-Water-Management Alternatives 65

Table 14. Changes in model-calculated average annual, March, and September stream base flow from 1989-98 base flow 

for hypothetical scenarios of increased withdrawals and altered recharge in the upper Charles River Basin, eastern 

Massachusetts—Continued

Streamflow site

Modelcalculated

1989–98

flow

(ft3/s)

Scenario

Change in model-calculated

stream base flow

Average 

annual

flow 

(ft3/s)

Average 

March

flow

(ft3/s)

Average

September

flow

(ft3/s)

West Model

Charles River

Charles River above the Milford Treatment Facility near the 

Milford/Hopedale town line ......................................................... 10.5 1

2

-2.1 (20)

-2.1 (20)

-2.3 (9)

-2.3 (9)

-0.7 (69)

-0.7 (69)

Charles River at 01103140 .............................................................. 18.6 1

2

-0.16 (0.9)

-0.16 (0.2)

+0.27 (0.7)

+0.27 (0.7)

+0.97 (18)

+0.97 (18)

Charles River at 011032056 ............................................................. 27.9 1

2

-0.55 (2)

-0.55 (2)

-0.11 (0.2)

-0.11 (0.2)

+0.60 (9)

+0.59 (9)

Charles River at 01103260 ............................................................... 87.5 1

2

-2.6 (3)

-2.7 (3)

-2.4 (1)

-2.7 (2)

-0.19 (5)

-0.19 (5)

Tributaries to Charles River

Beaver Brook ................................................................................... 4.7 1

2

-0.21 (5)

-0.21 (5)

-0.21 (3)

-0.21 (3)

-0.19 (10)

-0.19 (10)

Stall Brook at 01103120................................................................... 4.1 1

2

-0.36 (9)

-0.36 (9)

-0.31 (3)

-0.31 (3)

-0.34 (68)

-0.34 (68)

Hopping Brook................................................................................. 15.8 1

2

+0.01 (0)

+0.16 (1)

0.00 (0)

+0.16 (0.5)

+0.02 (0.8)

+0.11 (4)

Miscoe Brook at South Street above Mine Brook............................ 1.2 1

2

3a

3b

3c

0.00 (0)

-0.45 (37)

+0.11 (9)

+0.02 (1)

+0.03 (2)

0.00 (0)

-0.73 (34)

+0.15 (7)

+0.02 (1)

+0.04 (2)

0.00 (0)

-0.17 (36)

+0.08 (18)

+0.01 (3)

+0.02 (5)

Mine Brook at 01103225.................................................................. 2.7 1

2

3a

3b

3c

-0.01 (0.3)

-0.66 (25)

+0.23 (9)

+0.04 (1)

+0.06 (2)

-0.02 (0.4)

-1.02 (23)

+0.30 (7)

+0.05 (1)

+0.08 (2)

0.00 (0.2)

-0.34 (27)

+0.19 (15)

+0.03 (3)

+0.05 (4)

Mine Brook at 01103235.................................................................. 12.8 1

2

3a

3b

3c

-0.74 (6)

-1.15 (9)

+0.99 (8)

+0.16 (1)

+0.27 (2)

-0.92 (4)

-1.6 (7)

+1.2 (6)

+0.20 (0.9)

+0.34 (2)

-0.43 (7)

-0.64 (11)

+0.58 (10)

+0.10 (2)

+0.17 (3)

Mine Brook at 01103240.................................................................. 19.1 1

2

3a

3b

3c

-1.1 (6)

-1.4 (7)

+1.4 (7)

+0.23 (1)

+0.39 (2)

-1.3 (4)

-1.8 (5)

+1.7 (5)

+0.28 (0.8)

+0.47 (1)

-0.81 (10)

-0.92 (11)

+0.96 (12)

+0.16 (2)

+0.27 (3)

Table 14. Changes in model-calculated average annual, March, and September stream base flow from 1989-98 base flow 

for hypothetical scenarios of increased withdrawals and altered recharge in the upper Charles River Basin, eastern 

Massachusetts

[Streamflow site: See figure 7 for location of measurement sites; locations of sites along tributaries are immediately upstream of confluences with the Charles 

River unless otherwise indicated. Model-calculated 1989–98 flow: Average of model-calculated monthly flows from transient simulations. Change in modelcalculated base flow: Number in parentheses is percent of model-calculated 1989–98 average annual or monthly flow; ft3/s, cubic feet per second]
66 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

1Not determined because simulated stream is dry in September under 1989–98 conditions.

East Model

Charles River

Charles River at 01103280 ............................................................... 96.3 1

2

-2.5 (3)

-2.6 (3)

-2.3 (1)

-2.5 (1)

-1.0 (4)

-0.96 (4)

Charles River above Populatic Pond ................................................ 97.2 1

2

-1.9 (2)

-2.0 (2)

-1.7 (0.9)

-1.9 (1)

-0.32 (1)

-0.35 (1)

Charles River above Mill River........................................................ 104 1

2

+0.17 (1)

+0.53 (2)

+0.03 (.04)

+0.60 (0.3)

+1.8 (6)

+2.2 (8)

Charles River at 01103305 ............................................................... 127 1

2

-0.97 (0.2)

-0.46 (0.5)

-1.38 (0.6)

-0.50 (0.2)

+0.97 (3)

+1.42 (4)

Tributaries to Charles River

Mill River between 01103292 and 01103295 .................................. 11.6 1

2

-0.60 (5)

-0.60 (5)

-0.87 (5)

-0.87 (5)

-0.35 (6)

-0.35 (6)

Mill River at 01103300 ................................................................... 19.4 1

2

-0.98 (5)

-1.0 (5)

-1.3 (4)

-1.2 (3)

-0.69 (12)

-0.79 (14)

Mill River above Charles River........................................................ 22.2 1

2

-1.1 (5)

-0.99 (5)

-1.4 (3)

-1.1 (3)

-0.77 (12)

-0.79 (13)

Cress Brook above Mill River.......................................................... 1.5 1

2

-0.07 (4)

+0.07 (4)

-0.07 (2)

+0.07 (2)

10

10

Miller Brook above Mill River......................................................... 3.3 1

2

-0.23 (7)

-0.37 (11)

-0.22 (3)

-0.33 (4)

-0.22 (85)

-0.26 (100)

Bogastow Brook and Tributaries

Dopping Brook at 01103386 ............................................................ 1.4 1

2

0.00 (0)

-0.49 (36)

0.00 (0)

-0.93 (27)

0.00 (0.1)

-0.04 (26)

Great Black Swamp tributaries above Bogastow Brook .................. 7.7 1

2

+0.03 (0.4)

-0.59 (8)

+0.04 (0.3)

-0.65 (5)

+0.01 (0.2)

-0.4 (10)

Bogastow Brook above Great Black Swamp tributaries .................. 21.4 1

2

-0.10 (0.4)

-0.43 (2)

0.00 (0)

-0.8 (2)

-0.18 (3)

-0.06 (0.9)

Bogastow Brook at 01103393 .......................................................... 29.1 1

2

-0.07 (0.2)

-1.0 (4)

+0.03 (0.1)

-1.4 (3)

-0.18 (2)

-0.5 (4)

Table 14. Changes in model-calculated average annual, March, and September stream base flow from 1989-98 base flow 

for hypothetical scenarios of increased withdrawals and altered recharge in the upper Charles River Basin, eastern 

Massachusetts—Continued

Streamflow site

Modelcalculated

1989–98

flow

(ft3/s)

Scenario

Change in model-calculated

stream base flow

Average 

annual

flow 

(ft3/s)

Average 

March

flow

(ft3/s)

Average

September

flow

(ft3/s)
Evaluation of Ground-Water-Management Alternatives 67

260

240

220

20

200

180

160

140

120

60

100

80

40

0

0

20

40

30

10

0 2 4 6 8 10 12 14 16 18 20 22 Louisa Lake Tributary MF-03G,-04G,-05G MF-02G MF-01S MF-01G Godfrey Brook Beaver Brook Stall Brook BL-12G Hopping Brook MTF Mine Brook Chicken Brook Shepards Brook WEST MODEL EXIT EAST MODEL ENTRANCE MD-01G MD-03G FL-08G FL-02P Mill River CRPCD Louisa Lake Tributary MF-03G,-04G,-05G MF-01G MF-02G Godfrey Brook Beaver Brook Stall Brook BL-12G Hopping Brook MTF Mine Brook Chicken Brook Shepards Brook WEST MODEL EXIT EAST MODEL ENTRANCE MD-01G MD-03G FL-08G FL-02P Mill River CRPCD MF-01S

DISTANCE ALONG SIMULATED STREAM, IN MILES

MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND

1989–98 

SCENARIO 1

SCENARIO 2

FLOW SITE

A.

B.

Figure 25. Model-calculated monthly average base flow for March and September along the Charles River 

for 1989–98 average conditions and two scenarios of increased withdrawals in the upper Charles River 

Basin, eastern Massachusetts: (A) March, and (B) September. Flow sites shown are streamflowmeasurement sites (solid symbols; see figure 7) or streamflow-observation sites described in table 14 (open 

symbols).
68 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

3

-3

0

-6

JF N M A A M J J S O D JF N M A A M J J S O D

A. Above MTF at the town line between

Milford and Hopedale 

B. Below MTF, at station 01103140

3

-3

0

-6

C. Between Beaver and Stall Brooks,

at station 011032056

D. Between Chicken Brook and

Shepards Brooks, at station 01103260

E. Between CRPCD

and Mill River

F. Below Mill River, 

at station 01103305

3

-3

0

-6

SCENARIO 1

SCENARIO 2

CHANGE IN MODEL-CALCULATED BASE FLOW FROM 1989–98 FLOW,

IN CUBIC FEET PER SECOND

Figure 26. Change in model-calculated base flow from 1989–98 flows at selected sites along the Charles River 

for two scenarios of increased withdrawals in the upper Charles River Basin, eastern Massachusetts: (A) Charles 

River above the Milford Treatment Facility (MTF) at the town line between Milford and Hopedale, (B) Charles 

River below the MTF at station 01103140, (C) Charles River between Beaver and Stall Brooks at station 

011032056, (D) Charles River between Chicken and Shepards Brooks at station 01103260, (E) Charles River 

between the Charles River Pollution Control District Treatment Facility (CRPCD) and Mill River, and (F) Charles 

River below Mill River at station 01103305.

However, March and September typically were months 

of maximum and minimum monthly average flows, 

respectively, at most measurement sites. Thus, changes 

in base flow typically were proportionally largest or 

smallest relative to 1989–98 flows in these months. 

Along the Charles River, stream base-flow 

reductions from increased withdrawals were nearly 

balanced by flow augmentations from increased wastewater discharges in both scenarios overall. Results of 

S1 and S2 simulations were very similar (fig. 25). 

Along specific reaches of the river or during various 

months of the annual cycle, however, reductions or 

augmentations dominated (figs. 23 and 24). Generally, 

model-calculated base flows were reduced relative to 

model-calculated 1989–98 flows upstream of wastewater treatment facilities in most months and were

augmented relative to 1989–98 flows in low-flow 

months (July to October) downstream of treatment 

facilities (fig. 26). This results because withdrawals

at upstream sources affect stream baseflow over an 

extended period of time, rather than instantaneously. 

Relative to total average monthly flow for 1989–98, 

reductions in model-calculated base flows were greatest in September (table 14) or October. Base-flow 

reductions were associated with increased withdrawals 

(a) from sources along the Charles River, such as wells 

in Milford and Bellingham, and (b) from decreased 

inflows from tributaries, such as from Mine Brook. 
Evaluation of Ground-Water-Management Alternatives 69

0

0

20

40

30

10

5

13

50

45

35

25

15

0

10

5

15

1 2 4 6 8 10 12 14 3 5 9 11 7 FR-03G

FR-01P

Dix Brook

Miscoe Brook Mine Brook FR-06GFR-01G, -02G MGCC-01G

FR-07G

FR-03G

FR-01P

Dix Brook

Miscoe Brook Mine Brook FR-06G FR-01G, -02G MGCC-01G

FR-07G

DISTANCE ALONG SIMULATED STREAM, IN MILES

MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND

A.

B.

1989–98 

SCENARIO 1

SCENARIO 2

FLOW SITE

Figure 27. Model-calculated monthly average base flow for March and September along Mine Brook for 

1989–98 average conditions and two scenarios of increased withdrawals in the upper Charles River Basin, 

eastern Massachusetts: (A) March, and (B) September. Flow sites shown are streamflow-measurement 

sites (solid symbols; see fig. 7) or streamflow-observation sites described in table 14 (open symbols).
70 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Model-calculated September base flow was reduced 

relative to 1989–98 flows along three river reaches in 

S1 and S2: by about 70 percent (0.7 ft3/s) downstream 

of Milford sources, by about 2 to 4 percent (0.2 to 0.3 

ft3/s) between Stall and Mine Brooks, and by about 5 

percent (1.1 ft3/s) downstream of Mine Brook (fig. 25B

and table 14). In March, reductions in model-calculated base flow along these three reaches were larger in 

magnitude than in September but were smaller fractions (9, 1.5, and 1.5 percent, respectively) of total flow 

(fig. 25A and table 14). Model-calculated augmentations from wastewater discharges exceeded upstream 

depletions in September; as a result, flow in September 

increased downstream of treatment facilities. In March, 

model-calculated augmentations (less than 1 percent of 

total flow) were less than upstream reductions. The 

increased discharges from the treatment facilities 

resulted in wastewater constituting greater proportions 

of stream base flow during low-flow months. Treated 

wastewater constituted about 90 percent of simulated 

September flows immediately downstream from the 

MTF (S1 and S2) and about 27 (S1) and 32 (S2)

percent of simulated September flows downstream 

from the CRPCD treatment facility; corresponding 

values for simulated 1989–98 conditions were 80 percent for the MTF and 18 percent for the CRPCD treatment facility.

Along Mine Brook, simulated streamflows were 

reduced relative to 1989–98 flows from increased withdrawals at several Franklin municipal-supply wells and 

a golf course located along the stream and its tributaries (fig. 27). Return flow was routed out of the subbasin 

to the CRPCD treatment facility from these withdrawals. Thus, reductions in model-calculated stream base 

flow in Mine Brook were comparable to upstream 

increases in withdrawals from the aquifers (0.83 and 

1.02 Mgal/day, or 1.29 and 1.58 ft3/s, in S1 and S2, 

respectively) and averaging 1.1 ft3/s (S1) and 1.4 ft3/s 

(S2) over the annual cycle (table 14). Base-flow

reductions along specific reaches of Mine Brook varied 

depending on the proximity of the reach to individual 

supply wells and varied as proportions of total flow. 

Model-calculated 1989–98 September base flow in the 

aquifer along Mine Brook was reduced in S2 by more 

than 50 percent in its headwaters (Miscoe Brook), and 

in S1 and S2 by about 10 percent along downstream 

reaches (fig. 27). Results of S1 and S2 are dissimilar in 

upstream reaches of Mine Brook because of the 

increased withdrawals at the single proposed well at the 

headwaters of Miscoe Brook, which is included in S2 

but not in S1.

Similar results to those described for Mine 

Brook were obtained for other tributaries in areas 

where return flow was routed to wastewater-treatment 

facilities. Reductions in model-calculated stream

base flow along Stall Brook and Beaver Brook in

Bellingham and along a Great Black Swamp tributary 

in Medway were about equal to upstream withdrawals 

from the aquifers (tables 9 and 14). Along Stall Brook, 

reductions in simulated stream base flow were substantial and represented more than 60 percent of modelcalculated 1989–98 flows during low-flow months for 

S1 and S2 (table 14). 

Along tributaries in areas where return flow

was wholly or partly returned to the aquifer as added 

recharge, reductions in model-calculated base flow 

generally were less than increases in upstream

withdrawals from the aquifers. Reductions in modelcalculated base flow in tributaries were not eliminated 

even in areas where return flow was routed to the aquifer, primarily because of differences in the contributing 

areas to the streams and areas where return-flow 

recharge was applied (that is, intrabasin and interbasin 

transfers). For example, in S2, average annual withdrawals of 0.54 ft3/s from a new source in the aquifer 

along Dopping Brook (0.35 Mgal/d at HL-01P;

table 10) were distributed throughout the town of

Holliston. Thus, simulated average annual base flow in 

Dopping Brook was reduced by 0.49 ft3/s (measurement site 01103386; this reduction includes effects of 

an increased average annual withdrawal of 0.06 ft3/s at 

an upstream golf course), whereas simulated base flow 

was slightly increased in Hopping Brook (average 

annual increase of 0.16 ft3/s at the confluence with the 

Charles River; table 14) and in Chicken Brook, where 

withdrawals were unchanged but returnflow increased. 

Along the Mill River, reductions in modelcalculated stream base flow relative to 1989–98 flows 

averaged 1.1 ft3/s in S1 and 0.99 ft3/s in S2 above the 

Charles River (table 14). Results of S1 and S2 are similar, because increased withdrawals at the two proposed 

wells included in S2 largely are balanced by smaller 

additional withdrawals at existing wells in S2 than in 

S1 (fig. 28). Return flow from sources in the Mill 

River aquifer area either was routed to the CRPCD 

treatment facility (Franklin withdrawals) or to the aquifer as added return flow (Norfolk and Wrentham withdrawals), and average reductions in base flow over the 

annual cycle are about 75 percent of average increased 

withdrawals at the confluence with the Charles River. 
Evaluation of Ground-Water-Management Alternatives 71

0

20

40

30

10

5

50

45

35

25

15

0

8

4

10

MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND

6

2

0 1 2 4 6 8 10 35 9 7 WR-01P Miller Brook

Eagle Brook Mill River FR-05G

FR-04G

NF-01P

WR-02G, -03G

Cress Brook

WR-01P

Miller Brook

Eagle Brook Mill River FR-05G

FR-04G

NF-01P

WR-02G, -03G

Cress Brook

DISTANCE ALONG SIMULATED STREAM, IN MILES

A.

B.

1989–98 

SCENARIO 1

SCENARIO 2

FLOW SITE

Figure 28. Model-calculated monthly average base flow for March and September along the Mill River for 

1989–98 average conditions and two scenarios of increased withdrawals in the upper Charles River Basin, 

eastern Massachusetts: (A) March, and (B) September. Flow sites shown are streamflow-measurement 

sites (solid symbols; see fig. 7).
72 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Reductions in model-calculated base flow in the

September low-flow period for S1 and S2 are about 6 

percent of model-calculated 1989–98 flows downstream of Wrentham withdrawals and about 12 to 14 

percent of 1989–98 flows downstream of Norfolk and 

Franklin withdrawals (table 14 and fig. 28B). It should 

be noted that simulated base flow in the Mill River may 

be affected by aquifer geometry. Constrictions in the 

areal extent of the aquifer (fig. 13B) result in enhanced 

ground-water discharge to streams, some of which is 

returned to the aquifer as stream leakage downstream 

of the constrictions. This effect is indicated by the rapid 

increases in model-calculated September base flow 

from about stream miles 2.5 to 4.0, upstream of an 

aquifer constriction, and abrupt decrease in modelcalculated base flow where the aquifer broadens near 

flow-measurement site 01103295 (fig. 28B). The effect 

may be larger in model-calculated base flow than in 

actual flows because the constriction is more abruptly 

defined, relative to actual aquifer boundaries, in the 

active model area. 

Increased withdrawals in the Mill River aquifer 

area also resulted in changes in model-calculated base 

flow in small tributaries in the Populatic Pond area

(fig. 3 and table 14) and in water levels in Kingsbury 

Pond. Model-calculated September base flow in Miller 

Brook, near Franklin wells No. 4 and No. 5 (FR-04G 

and FR-05G) and the Norfolk proposed well (NF-01P), 

is reduced by 85 percent relative to 1989–98 flows in 

S1 and is entirely depleted in S2. In Cress Brook, 

model-calculated average annual flows are reduced by 

about 4 percent in S1 when withdrawals at the nearby 

existing Norfolk well are increased. In S2, when 

increased withdrawals for Norfolk are obtained from 

the proposed well along the Mill River, flows in Cress 

Brook are augmented by 4 percent because of increased 

septic-system return flow in this area. 

Water levels in Kingsbury Pond, which is near 

several existing wells and between two proposed wells 

in the Populatic Pond area, are reduced by about 0.8 ft 

on average over the annual cycle in both S1 and S2

(fig. 29). Kingsbury Pond is a kettle pond without surface water inflow or outflow. Water levels in the pond 

are representative of water levels in the aquifer and 

indicate a general lowering of the water table in 

response to increased withdrawals in this area. The 

response of Kingsbury Pond was similar in S1 and S2 

even though increases in withdrawals from wells in the 

aquifer near the pond (northern part of Mill River aquifer and Charles River aquifer, table 10) in S1 were 

nearly twice the increases in withdrawals from these 

wells in S2. This resulted because withdrawals from 

wells closest to the pond (particularly FR-04G and 

FR-05G) were increased less in S2 than in S1. 

In summary, the increased withdrawals in the 

upper Charles River Basin in scenarios 1 and 2 resulted 

in localized reductions in model-calculated stream base 

flow relative to model-calculated 1989–98 flows. 

Reductions ranged from less than 5 percent to more 

than 60 percent of simulated 1989–98 base flow during 

low-flow periods in the Charles River and its major 

tributaries. The largest reductions occurred in headwater areas or in small subbasins where withdrawals 

were large relative to base flow and where the withdrawn water was transferred out of the subbasin. Additional recharge from septic-system return flow from 

increased withdrawals, as simulated in S1 and S2, augmented model-calculated base flow in a few tributaries 

where withdrawals were not increased, and to some 

extent decreased the effects on model-calculated base 

flow in areas of increased withdrawals. However, the 

effects were not large. Overall, increased water withdrawals were transferred among subbasins but were 

134

132

131

130

133

126

128

129

JFM M A ASOND J J

127

SCENARIO 1

1989—98 OBSERVED

SCENARIO 2

 MODEL-CALCULATED WATER LEVEL, 

IN FEET ABOVE SEA LEVEL

EXPLANATION

Figure 29. Model-calculated average monthly water 

levels in Kingsbury Pond for 1989–98 conditions and in 

two hypothetical scenarios of increased water 

withdrawals in the upper Charles River Basin, eastern 

Massachusetts.
Evaluation of Ground-Water-Management Alternatives 73

returned nearly entirely to the stream and aquifer 

system within the model area. Thus, overall changes in 

simulated outflows from the study area were small, 

about 1 to 1.5 ft3/s on average over the annual cycle in 

S1 and S2, and resulted primarily (80 percent) from 

consumptive use. In S2, water (about 0.5 ft3/s on average) was transferred from the Bogastow Brook drainage area to the Charles River drainage area through 

withdrawals at a proposed well in Medway and routing 

of the water to the CRPCD treatment facility. Finally, 

increased withdrawals and transfer to treatment facilities potentially increased the effects of these discharges 

on water quality in the river by increasing the wastewater component of streamflow during low-flow periods.

Simulation of Altered 

Recharge

A third set of scenarios, in which recharge to the 

aquifer is increased through water-management alternatives (S3), was simulated for the aquifer along Mine 

Brook in Franklin. This area was selected to represent a 

typical area of the model domain where streamflow 

depletion currently occurs and is likely to increase with 

additional withdrawals. In the Mine Brook Subbasin, 

water is transferred out of the subbasin by sewering 

and discharge of treated wastewater to the Charles 

River at the CRPCD facility in Medway. First, recharge 

was increased by an amount equal to the out-ofsubbasin transfer currently resulting from sewering 

(S3a). About 50 percent of the Mine Brook Subbasin is 

sewered (fig. 6). Although it is unlikely that existing 

sewers would be removed, this scenario serves to quantify the effects of sewering on streamflow depletion in a 

typical aquifer area. Second, recharge was increased to 

simulate artificial recharge of residential rooftop runoff 

(S3b and c). Recharge of rooftop runoff is a water-management practice that often is proposed in densely populated areas to reduce stormwater flows. This scenario 

illustrates the relative magnitude of potential effects of 

this practice on water quantity in a less densely populated, suburban area. Recharge volumes associated 

with the management practices simulated in S3a and 

S3b were estimated by CRWA, as described below 

(N.B. Pickering, Charles River Watershed Association, 

written commun., 2002). In all three scenarios, water 

withdrawals, inflows from precipitation recharge, and 

1989–98 septic-system return flow remained the same 

as those used to simulate 1989–89 conditions in the

calibrated models.

Added recharge equivalent to water currently 

transferred out of the Mine Brook Subbasin through 

sewering was estimated at 1.47 in/yr. This value was 

based on an analysis of land-use data (1990), the areal 

extent of sewers, and population densities per land use 

category similar to that described previously for 1989–

98 septic-system return flow (N.B. Pickering, Charles 

River Watershed Association, written commun., 2002; 

table 15). Population values determined from 1990 data 

were adjusted using 2000 data to represent average 

1989–98 populations. Sewered land uses were grouped 

into categories of low-density residential (lots greater 

than one-half acre in area), medium-density residential 

(lot sizes from one-quarter to one-half acre), multifamily residential land use, and other developed land 

uses (table 15). Wastewater volumes of 140 gal/day

per person were estimated from 1989–98 flows to 

CRPCD from Franklin (N.B. Pickering, Charles River 

Watershed Association, written commun., 2002); this 

value includes ground-water infiltration into sewer 

lines. The additional recharge was added uniformly in 

all months of the annual cycle. 

1Includes commercial, industrial, transportation, urban open and public, and spectator-recreation land uses.

Table 15. Land-use characteristics, population density, and rooftop-drainage characteristics used to calculate rates of 

additional recharge in the Mine Brook aquifer area, upper Charles River Basin, eastern Massachusetts

[Data from N.B. Pickering, Charles River Watershed Association, written commun., 2002; --, not applicable]

Land use category

Area (acres) Population 

density

(people

per acre)

Number of 

rooftops

Effectiveness of rooftop 

impervious surface

(percent of total rooftop area)

Total Sewered Unsewered Low High

Low-density residential................. 760 240 520 2.2 832 7.5 12.5

Medium-density residential .......... 1,023 692 331 4.4 2,238 34 56

Multi-family residential ................ 134 87 47 15.2 340 45 75

Other developed1........................... 1,186 260 926 2.2 -- -- --

Undeveloped ................................. 5,038 0 5,038 0 -- -- --
74 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

0

0

20

40

30

10

5

13

50

45

35

25

15

10

5

0

15

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

DISTANCE ALONG SIMULATED STREAM, IN MILES

MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND MODEL-CALCULATED BASE FLOW, IN CUBIC FEET PER SECOND

1989–98 

SCENARIO 3a

SCENARIO 3b 

SCENARIO 3c 

FR-03G

FR-01P

Dix Brook

Miscoe Brook Mine Brook FR-06GFR-01G, -02G MGCC-01G

FR-07G

FR-03G

FR-01P

Dix Brook

Miscoe Brook Mine Brook FR-06G FR-01G, -02G MGCC-01G

FR-07G

A.

B.

FLOW SITE

Figure 30. Model-calculated monthly average base flow for March and September along Mine Brook for 

1989–98 average conditions and three scenarios of altered recharge in the upper Charles River Basin, 

eastern Massachusetts: (A) March, and (B) September. Flow sites shown are streamflow-measurement 

sites (solid symbols; see fig. 7) or streamflow-observation sites described in table 14 (open symbols).
Evaluation of Ground-Water-Management Alternatives 75

Added recharge resulting from artificial recharge 

of rooftop runoff was estimated as low and high values 

of 0.24 in/yr (3b) to 0.40 (3c) in/yr. These values were 

based on several assumptions (N.B. Pickering, Charles 

River Watershed Association, written commun., 2002). 

First, the number of residential rooftops in the Mine 

Brook aquifer area was estimated from land-use data, 

population densities per land-use category, and an

average number of people per residential household 

(table 15). The average number of people per household, 2.9, was determined from population data for 

Franklin and Millis. Multifamily residential structures 

were assumed to contain three households. Second, an 

average rooftop area of 2,000 ft2 was assumed, based 

on an analysis of orthophotographs (MassGIS, 1997; 

N.B. Pickering, Charles River Watershed Association, 

written commun., August 2002). Third, the fraction of 

water currently drained from residential rooftops such 

that it becomes surface-water runoff (effectiveness of 

rooftop impervious area) was assumed to range from 

7.5 percent, a low estimate for low-density residential 

land use, to 75 percent, a high estimate for multifamily 

residential land use. Finally, it was assumed that all of 

the water currently drained from residential rooftops 

could be captured and recharged to the aquifer. The 

added recharge was added uniformly in all months of 

the annual cycle. 

Model-calculated stream base flow in Mine 

Brook increased under the three scenarios of increased 

recharge (S3a, b, and c; table 14). Average annual 

increases in simulated base flow at measurement site 

01103240, the downstream end of the subbasin area 

where recharge was added, were 1.4 ft3/s for S3a 

(effects of sewers), 0.23 ft3/s for S3b (rooftop runoff 

recharge—low estimate), and 0.39 ft3/s for S3c (rooftop runoff recharge—high estimate). The increases 

were nearly equivalent to the volumes added in the 

drainage basin for the three scenarios: 1.54 ft3/s for 

S3a, 0.25 ft3/s for S3b, and 0.42 ft3/s for S3c, as 

expected under conditions of dynamic equilibrium. The 

streamflow augmentations increased in absolute magnitude downstream as drainage area and total flow 

increased (fig. 30). Model-calculated streamflow augmentations were about 12 percent of model-calculated 

1989–98 September base flow for the scenario of sewer 

effects and about 2 to 3 percent of September base flow 

for the scenarios of rooftop-runoff recharge (table 14). 

The effect of added recharge on model-calculated base 

flow was linear, but varied throughout the annual cycle 

(fig. 31). Recharge added in low-flow months went 

partly to replenish aquifer storage, such that increases 

in model-calculated base flow in May through October 

were less than (80 percent of) the added recharge rate. 

In high-flow months, from November through April, 

added water was released from storage and streamflow 

increased by amounts greater than (120 percent of) the 

added recharge. 

Simulation-Optimization of 

Water Withdrawals and 

Stream Base Flow in the 

Populatic Pond Area

Simulation-optimization methods were applied 

in a subarea of the study area near Populatic Pond and 

the confluence of the Charles and Mill Rivers (fig. 32). 

In the Populatic Pond area, water is withdrawn from 

the stratified glacial aquifer for public supply for

three towns—Franklin, Medway, and Norfolk. 

2.0

1.6

1.2

0.8

0.4

0.0

0.0 0.4 0.8 1.2 1.6 2.0

Mar.

Jan.

Nov.

May

July

Sept.

ADDED RECHARGE,

IN CUBIC FEET PER SECOND

INCREASE IN BASE FLOW,

IN CUBIC FEET PER SECOND

1.8

1.8

1.4

1.0 1.4

1.0

0.6

0.6

0.2

0.2

S3a

S3c

S3b

Figure 31. Relation of increases in model-calculated

monthly average base flow in Mine Brook (measurement site 

01103240) and added recharge in the Mine Brook Subbasin, 

upper Charles River Basin, eastern Massachusetts. Data 

points are for three model runs in which areal recharge is 

increased by 0.25, (S3b) 0.42 (S3c), and 1.54 (S3a) cubic

feet per second (0.24, 0.40, and 1.47 inches per year) to 

simulate recharge of rooftop runoff (0.25 and 0.42 cubic feet 

per second) and the effects of out-of-basin transfers by 

sewering (1.54 cubic feet per second).
76 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

MD-01G

MD-03G

FR-08G

FR-02P

NF-01G

FR-04G

NF-01P

FR-05G

EAST MODEL ACTIVE AREA 

CRPCD

EAST MODEL INACTIVE

 UPLAND AREA 

SIMULATED STREAM

TOWN BOUNDARY

C1

C2

C3

C4

M2

M1

K1

Kingsbury

Pond

Populatic

Pond

Charles RiverMill River

FR-05G MUNICIPAL SUPPLY WELL—See table 2 for

 identification information

C2 STREAMFLOW OBSERVATION POINT FOR

 MANAGEMENT MODELS—See table 15 for

 identification information

K1 POND-LEVEL OBSERVATION POINT FOR

 MANAGEMENT MODELS—See table 15 for

 identification information

0 0.5 1 MILE

0 0.5 1 KILOMETER

CRPCD WASTEWATER DISCHARGE—Charles River

 Pollution Control District Treatment Facility

MEDWAY

FRANKLIN

NORFOLK

N

EXPLANATION

Figure 32. Aquifer area near Populatic Pond where optimization methods were applied in the upper Charles 

River Basin, eastern Massachusetts.
Evaluation of Ground-Water-Management Alternatives 77

Water levels in ponds as well as streamflows are 

affected by these withdrawals. Simulation-optimization 

methods can provide information about how to minimize the effects of water withdrawals on the streamaquifer system, or, conversely, on how to maximize 

water withdrawals from an area within specific limitations (constraints) on hydrologic responses. The simulation-optimization approach was applied in the 

Populatic Pond area to test its efficacy for water management in a typical area of the upper Charles River 

Basin where multiple water-supply and water-resource 

needs coexist. Three management questions, identified 

through consultation with the Technical Advisory 

Committee, MADEM, MADEP, and CRWA, were 

posed: 

1. Can withdrawals at existing and proposed wells 

be increased, and by how much, while 

maintaining flow in the Charles and Mill Rivers 

above specified minimum-flow requirements?

2. Can existing withdrawal volumes be increased 

with no additional streamflow depletion?

3. Can streamflow depletion in the Charles River be 

reduced relative to current conditions while 

maintaining current water withdrawals?

Theory and Methods of the 

Simulation-Optimization Approach

The simulation-optimization approach relies on 

the numerical ground-water-flow model to simulate the 

hydrologic response of the stream-aquifer system to 

applied stresses such as water withdrawals. Optimization techniques are used to formulate and answer specific management questions about the applied stresses 

and hydrologic responses. The mathematical representation of these questions can be termed a management 

model that is formulated and solved separately from the 

flow model. Further information about the use of the 

simulation-optimization methods for ground-water 

resource management can be found in Ahlfeld and 

Mulligan (2000), and in a recent application (Barlow, 

1997; Barlow and Dickerman, 2001).

A management model has three components: 

objective function, decision variables, and constraints 

(Ahlfeld and Mulligan, 2000). (1) The objective function represents the goal of the management process. It 

typically is represented by a mathematical formula in 

which some quantity (typically a stress on the aquifer) 

is maximized or minimized. For example, it may be 

desired to maximize total pumping (QwT) at three 

supply wells (Maximize QwT = Qw1 + Qw2 + Qw3).

(2) The decision variables are the quantities that are to 

be determined in the solution of the management question. In the example given, pumping rates at the individual supply wells (Qwi) are the decision variables. 

(3) The constraints set limits on the values that decision 

variables may take in the solution to the question. Constraints for the example given could be that total withdrawals from all supply wells must exceed a specified 

value (Qw1 + Qw2 + Qw3 > 1 Mgal/d) or that pumping 

rates at individual wells cannot exceed established 

upper limits (Qw1 < 0.5 Mgal/d; Qw2 < 0.3 Mgal/d; and 

Qw3 < 0.2 Mgal/d). Constraints also may consist of 

limitations on the hydrologic response of the streamaquifer system, such as maximum allowed quantities of 

streamflow depletion or minimum ground-water levels 

at specified observation points. 

The management model is linked to the groundwater flow model through a matrix of response

coefficients (Gorelick and others, 1993; Ahlfeld and 

Mulligan, 2000). Response coefficients quantify the 

relation between a stress applied to the stream-aquifer 

system and the response of the system to the stress, in 

terms of streamflow depletion or altered ground-water 

levels. Response coefficients are calculated for streamflow depletion and changes in water levels at specific 

observation points and times as:

(6)

where 

Ri,j,tp,tr = response coefficient describing change at 

observation point j during month tr

caused by a change in pumping at well i

during month tp, dimensionless for

streamflow depletion and with units of 

ft/(ft3/d) for changes in water levels, 

Di,j,tp,tr = streamflow depletion or water-level 

change at observation point j during 

month tr caused by a change in pumping 

at well i during month tp, in ft3/d for 

streamflow depletion and in ft for waterlevel change, and

Qwi,tp = specified change in pumping at well i

during month tp, in ft3/d.

Ri j tp tr , , ,

Di j tp tr , , ,

Qwi tp ,

= --------------------
78 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Response coefficients are determined for specific 

locations of interest (observation points) in the active 

model area. A separate response coefficient represents 

the hydraulic response at each observation point to 

changes in the pumping rate at each pumping well 

during each stress period.

The response-matrix approach assumes a linear 

relationship between aquifer stresses and hydrologic 

responses. The assumption of linearity allows multiple 

responses to be added through superposition. Linearity 

is assumed both with the magnitude of the stress 

(streamflow and water-level changes are proportional 

to changes in pumping) and with the time of the

stress (streamflow and water levels change equally

in response to equal changes in pumping at all times

in the annual cycle). The assumption of linearity

may not be valid under some conditions; for example, 

when the proximity of a stressed well to a flow boundary causes non-linear head responses. In the upper 

Charles River Basin, aquifer geometry and the distribution of permeable materials is such that supply wells 

generally are located near streamflow boundaries. 

Thus, the assumption of linearity of aquifer response to 

imposed stresses was tested for management questions 

investigated in this study, as described below.

A linear programming computer software package, LINDO (Shrage, 1997), was used to solve management questions posed for the Populatic Pond area in 

this study. Solutions consist of optimal values for all 

the decision variables, which, for the questions posed, 

were average monthly pumping rates at each well 

included in the management models. 

Application in the 

Populatic Pond area

In the area near Populatic Pond (fig. 32), the 

aquifer is narrow, and supply wells for the three municipalities in many cases affect the same ground- and

surface-water resources. Issues of concern for water 

management in this area are (1) the multiple watersupply needs, and (2) the effects of the withdrawals on 

streamflows in the Charles and Mill Rivers and on 

water levels in Kingsbury Pond. Water-supply needs 

and hydrologic effects of withdrawals are of particular 

concern in summer months, when water use typically is 

greatest and recharge, streamflows, and water levels are 

low. Management models that formalize three management questions identified for testing were developed. 

The models were formulated similarly in terms of 

hydrologic stresses (decision variables) and responses, 

but differed in their objective functions and in the

constraints set on management model solutions.

Hydrologic Stresses and Responses

Hydrologic stresses included in the management 

models of the Populatic Pond area are water withdrawals at existing and proposed supply wells for the 

towns of Franklin and Norfolk (table 16). Monthly 

average pumping rates at these six wells constitute the 

decision variables in the models. Pumping rates for the 

Medway supply wells (MD-01G and MD-03G) were 

not included as decision variables. Changes in pumping 

rates at these two wells, which are immediately adjacent to the Charles River (fig. 32), would have little 

effect on the Mill River or Kingsbury Pond. For this 

reason, and to minimize the complexity of the management models, pumping rates at Medway wells were 

maintained at average 1989–98 levels. Similarly, 

wastewater discharge from the CRPCD treatment facility also was maintained at average 1989–98 levels. 

These simplifications were adopted for the hypothetical 

questions posed in this study and likely would be

modified in a more comprehensive application of the 

simulation-optimization approach in this area. 

Streamflows in the management models were 

represented by flows at six representative observation 

points along the Charles and Mill Rivers (table 16 and 

fig. 32). Observation points are model-cell locations 

where response coefficients were evaluated and where 

constraints were set in the management models. Observation points were selected near the confluences of

the rivers and upstream and downstream of managed 

water withdrawals. Four points were selected along the 

Charles River (C1–C4) and two were selected along

the Mill River (M1 and M2; fig. 32). Kingsbury Pond,

a kettle pond, also was represented by an observation 

point. Water levels in Populatic Pond, which is an 

extension of the Charles River, are controlled by stream 

stage, and, thus, an observation point in this pond

was not included in the management models. Observation points C1 and M1 approximately coincided

with measurement sites 01103280 and 01103300, 

respectively.

Response coefficients for streamflow depletion 

and changes in head were calculated for each observation point (table 16) by means of flow simulations with 

the east model and equation 6. A base simulation was 

run in which withdrawals were set to zero for all wells 

included in the management models. Then, four 
Evaluation of Ground-Water-Management Alternatives 79

1Approximate location of streamflow measurement site 01103280.

2Location of streamflow measurement site 01103300.

Table 16. Pumping wells, stream locations, and ponds in management models of the Populatic Pond area, upper Charles River 

Basin, eastern Massachusetts

[Identifier: See figure 32 for locations. Minimum flow: Based on a minimum-flow requirement of 0.21 cubic feet per second per square mile; --, not applicable or not used]

Identifier Name

Model location Role in management 

model Drainage

area

(square

miles)

Minimum

flow

(cubic

feet per

second) Layer Row Column Decision

variable

Observation

and (or)

constraint

location

Pumping wells

FR-04G ......... Franklin well No. 4 1 243 245 X -- -- --

FR-05G ......... Franklin well No. 5 1 255 249 X -- -- --

FR-08G ......... Franklin well No. 8 1 216 231 X -- -- --

FR-02P .......... Franklin proposed well 1 217 234 X -- -- --

NF-01G ......... Norfolk well No. 1 1 224 273 X -- -- --

NF-01P.......... Norfolk proposed well 1 252 252 X -- -- --

Stream locations

C1.................. Charles River at Medway1 1 205 225 -- X 65.2 13.7

C2.................. Charles River above Populatic Pond 1 216 236 -- X 66.2 13.9

C3.................. Charles River between Populatic Pond and 

Mill River

1 216 250 -- X 66.9 14.0

C4.................. Charles River above Mill River 1 218 261 -- X 67.6 14.2

M1 ................. Mill River near Kingsbury Pond2 1 236 258 -- X 13.8 2.90

M2 ................. Mill River above Charles River 1 220 262 -- X 16.0 3.36

Ponds

K1.................. Kingsbury Pond 1 237 250 -- X -- --

simulations were run for each well in which pumping 

was simulated at the well for a 1-month period. 

Response coefficients were calculated by dividing 

streamflow depletion (Di,j,tp,tr) by the pumping rate 

(Qwi,tp), as in equation 6. 

To test the assumption of linearity used in the 

response coefficient methodology, a high pumping rate 

(1.0 Mgal/d) and a lower pumping rate (0.5 Mgal/d) 

were simulated for each well during both January and 

July (a total of four simulations per well). For the 

assumption of linearity to be valid, response coefficients at an observation point generated from the four 

simulations for the same pumping well (R values with 

common i,j values; equation 6) should be equal or 

nearly so. The lower pumping rate of 0.5 Mgal/d equals 

the average State-permitted maximum pumping rate for 

the wells considered, and the high pumping rate of

1.0 Mgal/d exceeds the State-permitted maximum 

pumping rate for the wells, except for NF-01P which 

has a proposed limit of 1.1 Mgal/d. For the wells, 

times, and withdrawal rates tested, the response coefficients with common (i,j) values typically varied by less 

than 10 percent. The variation among response coefficients was the highest for observation points and times 

when the effect of the increased pumping was very 

small. This is likely due at least in part to numerical 

model errors because the variation in coefficients did 

not show any consistent correlation with either pumping rate or time of pumping. For coefficients that were 

greater than 0.20, the differences among coefficients 

were in all cases less than 20 percent. These typically 

slight disagreements among response coefficients indicate that the stream-aquifer system in the Populatic 

Pond area may be weakly nonlinear. Although a 

sequential linearization approach could be used to 

address the nonlinearities, such an approach was considered unnecessary based on previous studies (Barlow 

and Dickerman, 2001). 
80 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

A mass-balance test of the stream-aquifer system 

response in the Populatic Pond area also was performed by analysis of response coefficients. Stream 

cells are the only boundaries, other than pumping 

wells, where water can exit the model in the management-model area. Thus, streamflow depletions resulting from a pumping stress should equal the total 

pumping withdrawal, under steady-state conditions, at 

a stream location sufficiently downstream from the 

pumping well. Under transient conditions, the timing 

of pumping withdrawals and streamflow depletions at 

observation points may be offset, but mass balance 

should hold when total flows are summed over the 

annual cycle. This mass-balance relation means that the 

sum of response coefficients for each month of the 

annual cycle at a downstream observation point should 

equal 1.0 (Barlow, 1997). This test was verified for all 

wells upgradient of streamflow observation points

(fig. 32). The response coefficients for streamflow 

depletion in the Charles and Mill Rivers for these wells 

satisfied the mass-balance requirement to within 3

percent (table 17). 

Final response coefficients for the supply wells 

and observation points in the management models were 

determined by selecting values from one of the four 

simulation runs for each Ri,j pair. The selected values 

were close to middle of the range of values and best 

met the mass-balance relation. Final response

coefficients are given in table 18.

The variability in hydrologic responses to pumping stresses in the model area, in terms of streamflow 

depletion, is illustrated by the distribution of response 

coefficients among months of the annual cycle from

pumping at each supply well for 1 month at 1 ft3/s 

(0.65 Mgal/d) (table 18). For example, pumping at 

FR02P, which is immediately adjacent to the Charles 

River, results in rapid streamflow depletion in the 

Charles River but only small depletions in the Mill 

River. More than 90 percent of streamflow depletion 

from the Charles River from 1 month of pumping at 

FR-02P has occurred by the end of the first month after 

the pumping month (table 18). In contrast, pumping at 

FR-04G for 1 month results in streamflow depletions 

that extend for 5 subsequent months, primarily in the 

Mill River. This variability is summarized in figure 33, 

which shows the time required for 90 percent of the 

depletion that occurs in the Charles and Mill Rivers 

from pumping at each well for 1 month. The time 

required for depletions to occur increases with increasing distance of the pumping well from a simulated 

stream. This variability is advantageous for optimization analysis, in which pumping schedules at individual 

wells are managed to minimize streamflow depletion 

during critical time periods.

Response coefficients for the Kingsbury Pond 

observation point also indicate the variable effects of 

pumping at supply wells in the management-model 

area on pond-water levels. The largest effect resulted 

from pumping at FR-04G, as would be expected from 

its close proximity to the pond. This relation of pumping at FR-04G and water levels in Kingsbury Pond is 

consistent with the previous understanding of the 

ground-water flow system near the pond (Williams, 

1967; Bouck, 1998). Water levels in Kingsbury Pond 

also were affected, to a lesser degree than by pumping 

at FR-04G, by pumping at FR-05G and NF-01P. 

Response coefficients indicate that pumping at FR-04G 

at rate of 1 ft3/s (0.65 Mgal/d) for 1 month results in a 

decline in water levels at the pond of about 0.5 ft in 

each of the first 2 months following the beginning of 

pumping, with additional declines in subsequent 

months. Continuous pumping at FR-04G at this rate 

would result in a total quasi-steady-state decline in 

water levels of about 2 ft (monthly response coefficients are additive). These results are consistent with 

the results of the two scenarios of increased pumping 

described previously. It should be noted, also, that 

determination of the effects of pumping at supply wells 

on pond levels (as well as on streamflow) by analysis of 

response coefficients is subject to at least the same

limitations and assumptions as the numerical flow 

models, as wells as assumptions incorporated in the 

optimization approach. 

Table 17. Hydrologic response coefficients for streamflow 

depletion at downstream observation points from pumping 

wells used to verify mass balance in the management models 

of the Populatic Pond area, upper Charles River Basin, 

eastern Massachusetts

[See figure 32 and table 16 for observation-point locations]

Pumping

well

identifier

Sum of monthly response cofficients at

downstream observation points 

Mill River above 

Charles River

(C4)

Charles

River above 

Mill River

(M2)

Sum of

Mill and

Charles

Rivers

FR-04G ......... 0.79 0.24 1.03

FR-05G ......... .93 .05 .98

FR-08G ......... .03 .95 .98

FR-02P .......... .02 .98 .99

NF-01P.......... .96 .06 1.01
Evaluation of Ground-Water-Management Alternatives 81

Table 18. Hydrologic response coefficients for the Populatic Pond area of the upper Charles River Basin, eastern 

Massachusetts

[Response coefficents represent change in streamflow or water level per unit withdrawal (1 cubic foot per second) for 1 month of pumping at supply wells, 

dimensionless for changes in streamflow and in units of feet per cubic foot per second for changes in water level (1 cubic foot per second equals 86,400 cubic 

feet per day). Observation sites: C1-C4 and M1-M2 are streamflow sites, and K1 is a pond site; sites are described in table 16. Months: Month 1 is the 

month in which pumping occurs; --, no hydrologic response at observation point from pumping at supply well]

Supply 

well

Observation 

point

Months 

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

FR-04G ................. C1 -- -- -- -- -- -- -- -- -- -- -- --

C2 -- -- -- -- -- -- -- -- -- -- -- --

C3 0.007 0.030 0.045 0.037 0.022 0.016 0.010 0.006 -- -- -- --

C4 .007 .041 .052 .043 .031 .020 .015 .007 0.007 0.007 0.007 --

M1 .145 .192 .121 .073 .042 .025 .015 .010 .006 -- -- --

M2 .151 .220 .150 .098 .067 .040 .027 .018 .010 .007 -- --

K1 .504 .543 .362 .233 .155 .103 .065 .039 .026 .013 .013 --

FR-05G ................. C1 -- -- -- -- -- -- -- -- -- -- -- --

C2 -- -- -- -- -- -- -- -- -- -- -- --

C3 -- -- .007 .015 .007 .008 -- -- -- -- -- --

C4 -- -- .007 .015 .015 .009 .006 -- -- -- -- --

M1 .310 .307 .150 .073 .036 .019 .011 .007 -- -- -- --

M2 .295 .307 .159 .082 .042 .023 .013 .008 -- -- -- --

K1 .026 .103 .142 .116 .078 .065 -- -- -- -- -- --

FR-08G ................. C1 .015 .015 .015 -- -- -- -- -- -- -- -- --

C2 .401 .085 .023 .008 -- -- -- -- -- -- -- --

C3 .675 .192 .052 .017 .007 -- -- -- -- -- -- --

C4 .675 .195 .054 .018 .008 -- -- -- -- -- -- --

M1 -- .010 .006 -- -- -- -- -- -- -- -- --

M2 -- .012 .008 .006 .006 -- -- -- -- -- -- --

K1 -- .013 .013 -- -- -- -- -- -- -- -- --

FR-02P .................. C1 .009 .009 .009 -- -- -- -- -- -- -- -- --

C2 .319 .047 .013 .005 -- -- -- -- -- -- -- --

C3 .765 .147 .038 .012 .005 -- -- -- -- -- -- --

C4 .767 .151 .040 .014 .007 -- -- -- -- -- -- --

M1 -- -- -- -- -- -- -- -- -- -- -- --

M2 -- .005 .005 .005 -- -- -- -- -- -- -- --

K1 -- .015 .015 -- -- -- -- -- -- -- -- --

NF-01G ................. C1 -- -- -- -- -- -- -- -- -- -- -- --

C2 -- -- -- -- -- -- -- -- -- -- -- --

C3 -- -- -- -- -- -- -- -- -- -- -- --

C4 -- .015 .015 .015 .007 .005 -- -- -- -- -- --

M1 -- .006 -- -- -- -- -- -- -- -- -- --

M2 .111 .164 .112 .074 .050 .035 .025 .012 .010 .010 .013 0.007

K1 .006 .013 .019 .019 .019 .019 .013 .006 .006 -- -- --

NF-01P.................. C1 -- -- -- -- -- -- -- -- -- -- -- --

C2 -- -- -- -- -- -- -- -- -- -- -- --

C3 -- -- .015 .015 .007 .007 -- -- -- -- -- --

C4 -- -- .015 .015 .015 .007 .006 -- -- -- -- --

M1 .439 .214 .113 .062 .037 .021 .012 .009 .006 -- -- --

M2 .439 .221 .121 .071 .045 .025 .016 .010 .007 -- -- --

K1 .038 .115 .128 .103 .064 .051 .038 .013 .013 -- -- --
82 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Management Model Applications

Management models were developed and 

applied to address the three management questions 

described previously. For the first two questions,

the goals were to increase water withdrawals in the 

Populatic Pond area within specified constraints on 

streamflow depletion. For the third question, the goal 

was to reduce current streamflow depletion while

maintaining current total withdrawals. 

Increase Withdrawals while Maintaining 

Stream Base Flow above Specified 

Minimum-Flow Requirements

The first set of applications examined whether it 

would be possible to increase water withdrawals from 

existing and proposed sources in the Populatic Pond 

area, while maintaining streamflows in the Charles and 

Mill Rivers above specified minimum-flow requirements. Because water demand is greatest in June, July, 

and August, the objective function was formulated to 

maximize pumping during these months. The objective 

function for these applications is:

(7)

where NW is the total number of wells, NDt is total 

number of days in month t, and Qwi,t is the average 

pumping rate at well i in month t, as described 

previously. The value that is maximized is the sum of 

pumped volumes in June, July, and August from all 

wells in the management models (table 16). 

The value of the objective function was limited 

by constraints on streamflow depletion, pumping rates 

at wells, and, initially, the water level in Kingsbury 

Pond. Constraints were placed on streamflow depletion 

in the Charles River and Mill River, such that modelcalculated streamflows were not allowed to fall below 

the requirement of 0.21 (ft3/s)/mi2 at each of the 

streamflow observation sites listed in table 16. Imposition of the streamflow constraint at all streamflow 

observation points was designed to ensure that flows in 

the rivers were above the specified limit throughout the 

area of interest. A minimum-streamflow requirement 

of 0.21 (ft3/s)/mi2 currently is in use by MADEM for 

permitting activities that affect flows in the Charles 

River (V. Gartland, Massachusetts Department of

Environmental Management, written commun., 2002). 

It should be noted that minimum-flow requirements 

such as this one typically refer to total streamflow, 

rather than to the base-flow component of streamflow 

that is simulated by the ground-water flow models. 

However, during low-flow period such as summer 

months, most streamflow is base flow, so that use of 

minimum-flow requirements as constraints in the

management models was considered appropriate. 

Pumping rates at individual wells were constrained so as not to exceed MADEP-approved Zone II 

rates (table 2). The streamflow and Zone II constraints 

on pumping were imposed equally in all months of

the annual cycle. Pumping rates also were constrained 

so that combined pumping from all existing Franklin 

and Norfolk sources in each month of the annual cycle 

equaled or exceeded 1989–98 monthly average withdrawals, plus an added volume for each proposed well. 

10

0

1

2

3

4

5

6

7

8

9

0 2,500 500 1,000 1,500 2,000

FR-04G

FR-04G

FR-05G

FR-05G

NF-01G

NF-01G

NF-01P

NF-01P

FR-08G

FR-08G

FR-02P

FR-02P

CHARLES RIVER SITE AND IDENTIFIER

MILL RIVER SITE AND IDENTIFIER

DISTANCE TO NEAREST SIMULATED 

STREAM, IN FEET

TIME, IN MONTHS

PUMPING MONTH

FR-08G

FR-02P

EXPLANATION

Figure 33. Relation of the time required for streamflow 

depletion to occur in response to pumping for one month 

at supply wells and the distance of the supply well from 

the nearest simulated stream for six supply wells in the 

Populatic Pond area, upper Charles River Basin, 

eastern Massachusetts. Times shown are months, from 

the end of the pumping month, within which more than 

90 percent of the total streamflow depletion that results 

from the 1-month pumping stress in the nearest stream.

maximize NDtQwi t ,

t = 6

8

∑

i = 1

NW

∑
Evaluation of Ground-Water-Management Alternatives 83

The added volumes for proposed wells were equal

to one-half the estimated Zone II rate (table 2; B.R. 

Bouck, Massachusetts Department of Environmental 

Protection, written commun., 2002). These constraints 

were applied separately to Franklin and Norfolk 

sources, and ensured that total withdrawals for each 

town were at least as great as current withdrawals (as 

represented by 1989–98 averages). 

Initially, a constraint also was imposed on water 

levels in Kingsbury Pond. Water levels in Kingsbury 

Pond have declined in response to pumping since the 

1960s (Williams, 1967; Bouck, 1998). Current average 

water levels in Kingsbury Pond (table 5) were estimated at about 9 ft below pre-1960s levels of about

139 ft above mean sea level. Some septic-system

leachfields installed near the pond since these waterlevel declines have been below the 5-ft minimum depth 

to the water table required by State regulations (Title 

V) for septic systems (B.R. Bouck, Massachusetts 

Department of Environmental Management, oral

commun., 2002). Thus, a constraint was imposed to 

maintain water levels in Kingsbury Pond below

129.9 ft above sea level. This level was below the

average 1989–98 level. Inclusion of this constraint, 

however, resulted in no feasible solution to the management model. This outcome resulted because the 

high pumping rates required to keep water levels in 

Kingsbury Pond below the specified limit would have 

been greater than the maximum allowed Zone II pumping rates at wells. When constraints for maximum head 

in Kingsbury Pond were removed, a feasible solution 

was possible. 

Total withdrawals for June, July, and August

in the solution to this optimization problem averaged 

3.98 Mgal/d, more than twice the average withdrawals 

during 1989–98 in these months (1.55 Mgal/d). Average annual withdrawals were about doubled, from

1.38 Mgal/d to 2.7 Mgal/d. Most of the increases 

resulted from pumping at the proposed wells. However, 

pumping rates were increased at most wells in the management model in June, July, and August, such that all 

wells pumped at their maximum approved rates in 

these months (fig. 34A). Pumping rates were lower in 

months other than June, July, and August (although the 

total pumping for each town was still at or above 1989–

98 rates) because the other months were not part of the 

objective function. Rates also were lower in other 

months because lower rates in winter months help to 

maintain higher streamflow during the summer months 

when streamflow approaches the minimum flow rates.

These substantial increases in withdrawals 

were accompanied by additional streamflow depletion 

relative to model-calculated 1989–98 flows. Modelcalculated September streamflows were reduced by 

about 40 percent (2.4 to 3.0 ft3/s) at observation points 

in the Mill River, by about 5 percent (1.7 to 1.8 ft3/s) at 

observation points in the Charles River downstream of 

Populatic Pond, and by about 2 percent (0.6 ft3/s) at 

observation points in the Charles River upstream of 

Populatic Pond. Streamflow depletions in the Charles 

River downstream of Populatic Pond probably would 

be partly offset by increased discharges at the CRPCD 

treatment facility, which were not simulated. Modelcalculated stream base flows were not reduced to the 

constraint limit of 0.21 (ft3/s)/mi2 because the constraints on maximum pumping rates at individual wells 

were reached first. The increased withdrawals also 

were accompanied by a decline in model-calculated 

average annual water level in Kingsbury Pond of about 

0.3 ft.

Minimum streamflow requirements to maintain 

aquatic habitat are the subject of active research in 

Massachusetts (for example, Armstrong and others, 

2001). Alternative values could be considered applicable for streams in the upper Charles River Basin.

For example, the U.S. Fish and Wildlife Service 

(USFWS) recommends a minimum-flow requirement 

of 0.5 (ft3/s)/mi2 in summer months for New England 

streams (aquatic base flow or ABF; Armstrong and

others, 2001). Thus, alternative values to the

0.21 (ft3/s)/mi2 minimum-flow requirement may be 

considered in future watershed planning or permitting 

in the upper Charles River Basin. For this reason,

additional management models were formulated in 

which the effects of various minimum-streamflow 

requirements were investigated. Streamflowdepletion constraints were reformulated to impose

minimum-flow requirements of 0.35, 0.40, 0.43, and 

0.50 (ft3/s)/mi2. Other aspects of the management 

model, including the objective function, constraints on 

pumping, and the locations of streamflow constraints 

(observation points), remained the same. Each of the 

alternate minimum-flow requirements was implemented in a separate optimization formulation and 

solution. 
84 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

NF-01P

NF-01G

FR-02P

FR-08G

FR-05G

FR-04G

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

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

A. 1989–98

B. Increased withdrawals with minimum-flow

 requirement 

D. Increased streamflow with existing and proposed

 withdrawals for two towns separately

C. Increased withdrawals with no additional 

 streamflow depletion

E. Increased streamflow with existing and proposed

 withdrawals for two towns combined 

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

MONTH MONTH

SUPPLY WELLS

PUMPING RATE, IN MILLIONS OF GALLONS PER DAY

EXPLANATION

Figure 34. Monthly pumping rates for 1989–98 and for management-model applications for increased withdrawals 

or increased streamflow in the Charles River during the low-flow period in the Populatic Pond area of the upper 

Charles River Basin, eastern Massachusetts: (A) 1989–98, (B) Increased withdrawals with a minimum-streamflow 

requirement (0.21 cubic feet per second per square mile), (C) Increased withdrawals with no additional streamflow 

depletion, (D) Increased streamflow with existing and proposed withdrawals for two towns separately, and (E) Increased 

streamflow with existing and proposed withdrawals for two towns combined.
Evaluation of Ground-Water-Management Alternatives 85

The relation between total withdrawals in June, 

July, and August and minimum-streamflow requirements in the Populatic Pond area, as determined by the 

management models, is shown in figure 35. No reductions in total allowable summer withdrawals result

until the minimum-flow requirement is greater than 

0.39 (ft3/s)/mi2. For minimum-flow requirements of 

0.43 (ft3/s)/mi2 or greater, no solution could be found 

that allowed 1989–98 withdrawals in the area to be 

maintained. 

Increase Water Withdrawals while Maintaining

Stream Base Flow at Current Levels

A second set of model applications examined 

whether it would be possible to increase water withdrawals in the Populatic Pond area while maintaining 

current streamflows in the Charles and Mill Rivers. 

Current streamflows were represented by modelcalculated 1989–98 average monthly flows. Streamflow-depletion constraints were reformulated in these 

applications to allow no decrease in monthly average 

stream flows. The objective function was to maximize 

pumping in June, July, and August at existing and proposed sources. Constraints on pumping also were 

imposed initially such that combined total pumping for

each town (Franklin and Norfolk) was required to meet 

or exceed 1989–98 rates for each month with the addition of proposed sources as described previously.

However, no feasible solution was found with these 

constraints imposed. When these constraints were 

removed for the months of June, July, and August,

an optimal solution was found (fig. 34B). Total withdrawals for June, July, and August in the optimal solution averaged 1.95 Mgal/d, exceeding 1989–98 

withdrawals in these months by about 25 percent. Total 

average annual withdrawals were 2.19 Mgal/d, about 

60 percent greater than 1989–98 withdrawals. In this 

solution, combined withdrawals from Franklin sources 

were less than 1989–98 withdrawals in June and July, 

but greater in August (table 19). 

Most of the additional volume of water withdrawn in the optimal solution for this management 

model, relative to 1989–98 withdrawals, occurred at 

the proposed wells. However, pumping schedules at 

existing wells also differed from 1989–98 conditions to 

take advantage of the variability in hydrologic response 

to pumping stress at the wells. For example, FR-04G, 

which is farther from simulated streams than other 

wells in the management model (fig 30), was pumped 

at higher rates than in 1989–98 in June through

September and was not pumped at all from November 

through April. In contrast, pumping at FR-05G did not 

occur in May through July, but occurred at higher rates 

than in 1989–98 in August through December. 

Increase Stream Base Flow while Maintaining

Existing and Proposed Withdrawals

A third set of model applications was made to 

determine whether it would be possible to increase 

streamflows in the Populatic Pond area while maintaining withdrawals at 1989–98 levels, with added increments for proposed wells. This management model 

required that the objective function be formulated in 

terms of streamflows. The objective function for these 

applications is:

(8)

where QsdC4,8 is streamflow depletion in the Charles 

River at observation point C4 (Charles River above the 

Mill River) in August (eighth month). A more complex 

objective function could be formulated in which 

streamflow depletion is minimized at multiple 

observation points and at multiple times; however, that 

level of complexity was beyond the scope of this study. 

0.50

0.45

0.40

0.15

0.10

0.35

0.25

0.30

0.20

0.0

0.05

0 5 1234 6

TOTAL SUMMER PUMPING RATES, IN 

MILLIONS OF GALLONS PER DAY 

MINIMUM-STREAMFLOW REQUIREMENT, IN 

CUBIC FEET PER SECOND PER SQUARE MILE 

INFEASIBLE SOLUTION

ZONE II

RATES

EXCEEDED

Figure 35. Relation of minimum streamflow requirements 

and total summer pumping rates from managed wells, as 

determined by optimal solutions to management models in 

the Populatic Pond area of the upper Charles River Basin, 

eastern Massachusetts. Rates shown are sums of 

average daily pumping rates in June, July, and August for 

supply wells FR-04G, FR-05G, FR-08G, FR-02P, NF01G, and NF-01P (see tables 4 and 16).

minimizeQsdC4 8,
86 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Constraints on pumping and streamflow-depletion rates 

were the same as in the applications in which June, 

July, and August withdrawals were maximized and the 

minimum-flow requirement was set at 0.21 (ft3/s)/mi2.

August stream base flow in the Charles River 

above Mill River in the optimal solution of the objective function was 40.2 ft3/s. This value was about

1 ft3/s greater than model-calculated August base flow 

(39.3 ft3/s) for 1989–98 withdrawals at existing sources 

and withdrawals from proposed wells equal to one-half 

their estimated MADEP-approved rate as described 

previously (fig. 34C). In this solution, Franklin and 

Norfolk sources separately withdrew the same water 

volumes in all months of the year as in 1989–98, with 

added increments for proposed wells. Thus, Franklin 

and Norfolk wells were managed separately to achieve 

reductions in streamflow depletion. Franklin wells 

were managed by eliminating or reducing pumping

at wells adjacent to the Charles River (FR-08G

and FR-02P) in summer relative to 1989–98 rates, 

increasing pumping at these wells and at FR-05G in 

winter, and eliminating pumping at FR-04G in winter 

(table 20). Norfolk’s wells were managed by eliminating pumping in summer at NF-01G, which is the closer 

of the two wells to the Charles River, and increasing 

summer pumping at NF-01P, which is farther from the 

Charles River. The changes in pumping schedules also 

resulted in a decline in the average annual water level 

in Kingsbury Pond of about 1 ft.

An alternative formulation for this question

was tested to examine how collaborative management 

of water supplies could result in further reductions

in streamflow depletion. In this application, the constraint that Franklin and Norfolk sources separately 

withdrew the same water volumes as in 1989–98 (with 

increments for proposed wells) was replaced with a 

constraint that total withdrawals from all Franklin

and Norfolk sources jointly met this requirement.

This is equivalent to assuming that water-supply

systems for these towns were connected and regulated 

so that supplies could be shared. The optimal solution 

for this application resulted in August stream base

flow in the Charles River above Mill River of 45.8 ft3/s, 

an increase of about 14 percent relative to the optimal 

solution for separate management of water-supply

systems. The additional reductions in streamflow 

depletion resulted from decreased pumping at Franklin 

sources adjacent to the Charles River and increased 

pumping at NF-01G, in summer, and at NF-01P, 

throughout the year, relative to 1989–98 rates (fig. 34D

and E). Thus, collaborative management of water systems in the Populatic Pond area results in a better solution in terms of the model objective because the 

optimization analysis can make use of greater variability in hydrologic responses to pumping stresses to 

better manage individual withdrawals. 

Limitations of Management-Model

Applications

The management-model applications for the 

upper Charles River Basin described in this study are 

subject to limitations that should be considered when 

their results are interpreted. In all applications of management models in this study, allowable streamflow 

depletion is quantified relative to model-calculated, 

rather than observed, stream base flows. Thus, differences between observed and model-calculated flows 

(model calibration error) should be considered when 

the results of the management-model applications are 

interpreted. Theoretically, alternative formulations of 

the management models could specify streamflow 

depletion relative to observed rather than modelcalculated flows. Observed streamflows, however, are 

available only for a limited number of locations. In 

addition, simulation-optimization results would need to 

be interpreted in terms of differences between observed 

flow values, represented by estimated monthly averages for 1989-98, and model-calculated values (modelcalibration error). Model-calibration error, in general, 

resulted from multiple, poorly known sources of error 

and varied spatially and over the annual cycle. Thus, it 

could be difficult to formulate streamflow constraints 

that eliminate model calibration error consistently and 

accurately throughout the management model. This 

effort was beyond the scope of this study. 

The results of the management-model

applications that include streamflow constraints are

for average monthly base-flow conditions. They do not 

represent dry years or describe the variability in stream 

base flow within a month during the average year. In 

dry years, stream base flow is less than long-term average flows, represented by the 1989–98 time period in 

this study. In average and dry years, mean daily streamflows on some days in any month are less than the 

monthly average. Thus, mean daily streamflows occurring under pumping schedules determined by

management-model applications that include

minimum-flow constraints are likely to fall below the 

specified minimum flow constraints on some days, 
Summary 87

because the streamflow constraint is applied in the 

models as a limit on monthly average values. Finally, 

minimum-flow requirements typically are set in terms 

of total streamflow, which is easily measured, rather 

than base flow, which is calculated by the 

ground-water-flow models.

SUMMARY

Ground water is the primary source of drinking 

water for towns in the upper Charles River Basin, an 

area in eastern Massachusetts that is undergoing rapid 

growth. Stratified glacial aquifers in the basin are thin 

and discontinuous and are in close hydraulic connection with streams, ponds, and wetlands. Increased 

water withdrawals, combined with out-of-basin or 

downstream transfers of wastewater and decreased

natural recharge from changes in land use, have 

stressed water resources by contributing to streamflow 

depletion and lowering water levels in ponds and

wetlands. The stratified glacial aquifers extend across 

municipal boundaries, potentially leading to conflicts 

over water availability and downstream effects.

Numerical modeling tools were developed by the U.S. 

Geological Survey, in cooperation with the

Massachusetts Departments of Environmental Management and Environmental Protection, to evaluate the 

effects of water-management alternatives and to 

address the regional-scale needs for information on 

water resources in the basin.

The upper Charles River Basin is an area of

105 mi2 that encompasses most or substantial parts of 

Towns of Bellingham, Franklin, Holliston, Medway, 

Milford, Millis, Norfolk and Wrentham. Land use is 

primarily forested and residential. Population has 

increased by about 15 percent in the past decade, to 

about 125,000 in 2000. Water withdrawals from 33 

public-supply wells or wellfields, 2 surface-water

withdrawals from the Charles River, and several large 

non-municipal sources averaged 10.1 Mgal/d in 1989–

98 and are likely to increase. Six new public-supply 

wells were being permitted at the time of this study. 

About half of the population in the basin is served by 

public sewers, and the treated wastewater discharges

to the Charles River from facilities in Milford and 

Medway. 

Stratified glacial aquifers provide the primary 

source of water for municipal supply and large nonmunicipal users. Stratified glacial aquifers occur along 

the Charles River and major tributaries and are thin and 

irregular in shape. Areally, they cover about 50 percent 

of the study area. Stratified glacial aquifers consist of 

depositional sequences that include sand, gravel, silt, 

and clay. Maximum thicknesses of 120 to 130 ft are 

reached in a preglacial bedrock valley in the eastern 

part of the study area, but elsewhere aquifer thickness 

generally is 70 ft or less. Stratified glacial deposits in 

western and eastern parts of the study area are discontinuous, but are hydrologically connected by the 

Charles River. Ground water generally is unconfined, 

but local confining conditions are present in areas 

where silt and clay sediments are well developed.

Streamflow averaged 103 ft3/s in the Charles 

River at Medway during water year 2000. Average 

annual base flow, estimated for the period 1989–98 

from correlations of measurements made during the 

study with records at long-term stream-gaging stations, 

was 88 percent of average annual streamflow during 

water year 2000. Estimates were similarly made of 

average annual and monthly base flow at 24 sites where 

flow measurements were made during the study on the 

Charles River and its tributaries. Estimated mean 

annual base flow at these sites ranged from 73 to 122 

percent of measured annual base flow, on the basis of 

monthly measurements during base-flow conditions in 

water year 2000; measured values all were within the 

90-percent confidence intervals of estimated values. 

Measured ground-water levels at 50 observation 

wells ranged from 250 ft above sea level in headwater 

areas in the western part of the study area to less than 

140 ft above sea level in low-lying areas along streams 

near the eastern boundary. Measured annual fluctuations in ground-water levels varied with topography 

and position of the observation well with respect to 

aquifer boundaries. Annual fluctuations during the 

study ranged from 1 ft to more than 3 ft, and typically 

reached maximum elevations in late spring. Mean 

annual water levels for 1989–98, estimated by correlation with nearby long-term observation wells, were 

nearly identical with measured water levels during 

water year 2000.

Inflows to the upper Charles River Basin 

included ground-water recharge from precipitation, 

ground-water recharge from septic-system return flow, 

and wastewater discharge to streams from treatment 

facilities. These inflows were estimated at 100 to

130 Mgal/d for precipitation recharge, 3.1 Mgal/d for 

septic-system return flow, and 7.5 Mgal/d for wastewater discharges to streams. Outflows included stream 
88 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

base flow at downstream boundaries, water withdrawals, ground-water evapotranspiration, possible groundwater underflow at study area boundaries, and infiltration into sewer lines. Outflows through streamflow 

were estimated at 104 Mgal/d (90-percent confidence 

intervals for streamflow corresponding to 80 to

137 Mgal/d). Outflows for water withdrawals were 

10.1 Mgal/d, and outflows through ground-water 

evapotranspiration from wetlands were estimated at 

about 5 Mgal/d. Small transfers for water supply and 

wastewater discharge in and out of the study area balanced each other. The inflows and outflows considered 

in the water balance agreed within 20 percent. 

Steady-state and transient ground-water-flow 

models were developed for the upper Charles River 

Basin with MODFLOW-2000. Separate models were 

developed for east and west aquifer areas to improve 

numerical stability of the models. Stratified glacial 

aquifers only were included as active model areas. The 

models were calibrated to 1989–98 conditions of water 

withdrawals, water levels, and stream base flows. 

Recharge rates were simulated as varying spatially in 

active model areas and in inactive till and bedrock 

uplands; the rates varied by land use, surficial geology, 

and return flow. Recharge rates were calibrated to estimated 1989–98 average annual stream base flow with a 

parameter-estimation technique. Area-weighted rates 

were 21.5 and 23.0 in/yr for active areas of the west 

and east models, respectively, and 21.9 and 24.0 in/yr 

for flows to the active model area from till and bedrock 

uplands of the west and east models, respectively. 

Horizontal hydraulic conductivity ranged from 70 to 

290 ft/d and vertical hydraulic conductivity ranged 

from 5 to 50 ft/d in the calibrated models. In the transient models, specific yield varied from 0.15 to 0.28 

and specific storage was set at 1.0 x 10-3 in the

calibrated models. 

Total inflows to the calibrated steady-state 

models were 119 Mgal/d. In both west and east models, 

inflows from uplands, which were simulated as 

additional recharge at active model boundaries, were 

large components of total inflows. Inflows from 

uplands were about 70 percent of inflows to the west 

model and about 50 percent of inflows to the east 

model. Calculated outflows consisted primarily of discharge to simulated streams, which were 90 and 85 

percent of total outflows from the west and east models, respectively. Although small relative to total outflows, simulated ground-water evapotranspiration from 

wetlands was of similar magnitude (although much less 

well quantified) to water withdrawals; this result indicates the importance of this flux. Water withdrawals 

accounted for 7 percent of simulated total fluxes overall through the stream aquifer, and varied from about 

1.5 percent (Hopping Brook subbasin) to 13 percent 

(area near the headwaters of the Charles River) of total 

fluxes in subbasins of the study area. In the transient 

models, water withdrawals ranged from 4 to 5 percent 

of total monthly flows in winter and spring months to 

about 10 to 12 percent of total monthly flows in 

summer and early fall.

Two hypothetical scenarios of increased withdrawals in the upper Charles River Basin were simulated with the numerical flow models. Increasing 

withdrawals to levels allowed in the basin under the 

Massachusetts Water Management Act would result in 

withdrawals of about 15 Mgal/d, or about 50 percent 

more than current withdrawals. Simulated effects of 

these increased withdrawals include reductions in 

stream base flow that are greatest in proportion to total 

flow in late summer and early fall, but varied through 

the annual cycle as a result of the timing of increased 

withdrawals at upstream sources. Reductions in modelcalculated stream base flow ranged from less than 5 

percent to more than 60 percent of model-calculated 

1989–98 base flow during low-flow periods in the 

Charles River and major tributaries. 

Along the Charles River, reductions in stream 

base flow from increased ground-water withdrawals 

were nearly balanced by flow augmentations from 

increased wastewater discharges overall. Along specific reaches of the river, reductions or augmentations 

dominated. Generally, model-calculated base flow

was reduced relative to 1989–98 flows upstream of 

wastewater-treatment facilities in most months and was 

augmented downstream of treatment facilities in lowflow months. Model-calculated September base flow in 

the Charles River was reduced relative to 1989–98 

flows along three river reaches: by about 70 percent 

downstream of Milford sources, by about 2 to 4 percent 

between Stall and Mine Brooks, and by about 5 percent 

downstream of Mine Brook. The proportion of wastewater in the Charles River downstream of treatment 

facilities increased relative to 1989–98 conditions from 

80 to 90 percent of model-calculated September base 

flow in Milford and from 18 to 27 percent in Medway. 

Along tributaries such as Mine Brook, where 

return flow was routed to a wastewater discharge on the 

Charles River, reductions in model-calculated stream 

base flow was about equal to upstream withdrawals. 
Summary 89

Model-calculated 1989–98 September base flow in 

Mine Brook were reduced by more than 50 percent in 

headwaters (Miscoe Brook) and by about 10 percent 

along downstream reaches. Along tributaries in areas 

where return flow was wholly or partly returned to the 

ground-water flow system as added recharge, reductions in model-calculated base flow generally were less 

than increases in upstream withdrawals. Reductions in 

model-calculated base flow in tributaries were not 

eliminated even in areas where return flow was routed 

to the aquifer, primarily because of differences in the 

contributing areas to the streams and areas where 

return-flow recharge was applied. Along the Mill River, 

where return flow is routed to a treatment facility and 

to septic systems, average reductions in base flow

over the annual cycle are about 75 percent of average 

increased withdrawals. Reductions in model-calculated 

base flow for the Mill River during the September lowflow period were from about 12 to 14 percent of 1989–

98 flows. Increased withdrawals in the Mill River aquifer area also resulted in changes in water levels in 

Kingsbury Pond, which were reduced by about 0.8 ft 

on average over the annual cycle in scenarios of 

increased withdrawals. 

In a third set of hypothetical scenarios, additional recharge equal to the transfer of water out of a 

typical subbasin by sewers was simulated and was 

found to increase model-calculated base flows. Flow 

augmentations from this added recharge were about 12 

percent of model-calculated 1989–98 base flows, about 

equal to base flow reductions simulated in the scenarios 

of increased withdrawals. Addition of recharge equal to 

that available through artificial recharge of residential 

rooftop runoff had smaller effects, augmenting simulated September base flow by about 3 percent. The 

effect of added recharge on model-calculated base flow 

was linear, but varied throughout the annual cycle. 

Recharge added in low-flow months went partly to 

replenish aquifer storage. Thus, increases in modelcalculated base flow in low-flow months were less than 

(60 percent of) the added recharge rate, whereas in 

high-flow months, added water was released from storage and streamflow increased by amounts greater than 

(110 percent of) the added recharge.

Simulation-optimization methods were applied 

in the aquifer area near Populatic Pond and the confluence of the Mill and Charles Rivers to demonstrate the 

use of these methods in the basin. Water is withdrawn 

from six supply wells for the three towns in this area. 

Management models were developed to maximize 

water withdrawals within specified constraints on 

streamflow depletion and to minimize streamflow 

depletion within minimum requirements for water supply. Application of the simulation-optimization methods indicates that hydrologic responses to pumping 

from different supply wells vary in time and duration in 

the Populatic Pond area. This variability suggests that 

water withdrawals could be managed to minimize the 

effects of increased withdrawals on streams and ponds. 

Solutions of preliminary management models

for the Populatic Pond area suggested that water withdrawals could be substantially increased from existing 

and proposed sources while maintaining stream base 

flow in the Charles and Mill Rivers above minimum 

flow requirements with active management of withdrawals. Using a minimum-streamflow requirement of 

0.21 (ft3/s)/mi2, management model solutions indicated 

that annual average withdrawals from existing and proposed sources in the Populatic Pond area could be 

about doubled. The minimum streamflow constraints 

were applied on a monthly average basis, however,

so that mean daily flows below the specified levels 

would likely occur in average years. Managementmodel solutions also indicated that it was possible to 

increase or maintain existing withdrawals with higher 

minimum-streamflow requirements, but that 1989–98 

withdrawals could not be maintained for minimum 

flow requirements above 0.43 (ft3/s)/mi2. A second set 

of management models indicated that it might be possible to increase withdrawals by about 60 percent, while 

allowing no further reductions in stream base flow in 

the Charles and Mill Rivers, by managing withdrawals 

jointly between two towns in the Populatic Pond area. 

Withdrawals could be increased by about 25 percent if 

sources for the two towns are managed separately in 

these scenarios. Finally, a third set of management 

models indicated that base flow during the low-flow 

period in the Charles River might be increased without 

reducing withdrawals. The result was obtained by eliminating pumping in wells along the river in summer 

months, increasing pumping at these wells in the

winter, and increasing pumping at wells farther from 

the river in summer months. This solution again was 

improved by collaborative management of sources 

between towns. This analysis indicates that the simulation-optimization approach could be useful for

water management in areas of the upper Charles River 

Basin and elsewhere in Massachusetts where multiple 

water-supply and water-resource needs coexist.
90 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

REFERENCES

Ahlfeld, D.P., and Mulligan, A.E., 2000, Optimal 

management of flow in groundwater systems: San 

Diego, CA, Academic Press, 185 p.

Allen, W.B., Hahn, G.W., and Tuttle, C.R., 1963, 

Geohydrologic data for the upper Pawcatuck River 

Basin, Rhode Island: Rhode Island Water 

Resources Coordinating Board GWM-3, 1 sheet.

Amory Engineers, 1989, Test-well reconnaissance 

summary in Bellingham, Massachusetts: 

September 1999, variously paginated.

Anderson-Nichols & Company, Inc., 1994, Pumping 

test proposal BRP WS 18—Conduct pumping test 

70 gallons per minute or greater at the Point 

Property site test well No. 4-93 and No. 5-93, 

Millis, Massachusetts: July 1994, variously 

paginated.

_____1998, New source approval BRP WS 18—

Conduct pumping test 70 gallons per minute or 

greater at the Point Property site test well No. 4-93 

and No. 5-93, Millis, Massachusetts: August 1998, 

variously paginated.

Armstrong, D.S., Richards, T.A., and Parker, G.W., 

2001, Assessment of habitat, fish communities, 

and streamflow requirements for habitat 

protection, Ipswich River, Massachusetts, 1998–

98: U.S. Geological Survey Water-Resources 

Investigations Report 01-4161, 72 p.

Barlow, P.M., 1997, Dynamic models for conjunctive 

management of stream-aquifer systems of the 

glaciated northeast: Storrs, CT., Ph.D. Thesis, 

University of Connecticut, 256 p.

Barlow, P.M., and Dickerman, D.C., 2001, Numericalsimulation and conjunctive-management models 

of the Hunt–Annaquatucket–Pettaquamscutt 

stream-aquifer system, Rhode Island: U.S. 

Geological Survey Professional Paper 1636, 79 p.

Bent, G.C., 1995, Streamflow, ground-water recharge 

and discharge, and characteristics of surficial 

deposits in Buzzards Bay Basin, Southeastern 

Massachusetts: U.S. Geological Survey WaterResources Investigations Report 95-4234, 56 p.

_____1999, Streamflow, base flow, and ground-water 

recharge in the Housatonic River Basin, Western 

Massachusetts and parts of Eastern New York and 

Northwestern Connecticut: U.S. Geological 

Survey Water-Resources Investigations Report 98-

4232, 24 p.

Bouck, B.R., 1998, The upper Charles River 

dilemma—A growing basin, in National Water 

Well Association: Burlington, VT, Proceedings of 

the Northeast Focus Ground Water Conference, 

October 20–21, 1998, p. 94–107.

Coffin & Richardson, 1977, Report on pumping tests in 

the Godfrey Brook well field for the Milford Water 

Company: April 1977, variously paginated.

Commonwealth of Massachusetts, 2001, 

Commonwealth Communities, accessed March 20, 

2001, and May 20, 2002, at URL http:// 

www.state.ma.us/cc.

D.L. Maher Co., 1986, Hydrogeologcial report 

summarizing the results of the test well program 

undertaken within the Town of Franklin, Mass.: 

Spring 1986, variously paginated.

_____1989, Test well exploration adjacent to gravel 

packed well No. 1 during February 1989, Franklin, 

Massachusetts: March 1989, variously paginated.

_____1993, Town of Franklin, Massachusetts, 

Delineation of Zone II, existing well 3 and 

wellfield 6: July 1993, variously paginated.

_____1995a, Town of Franklin, Massachusetts, 

Delineation of Zone II, existing wells 1, 2, and 7: 

March 1995, variously paginated.

_____1995b, Results of a 48-hour prolonged pumping 

test on a new 12 inch diameter satellite well 

constructed within the Godfrey Brook wellfield for 

the Milford Water Company, Milford, Mass.: 

March 1995, variously paginated.

_____1996, Results of a prolonged pumping test, at 

site T.W. No. 1-86, Populatic St., Franklin, Mass., 

during October, 1994: May 1996, variously 

paginated.

_____1997, Zone II delineation for gravel packed well 

No. 2, Oakland Street, Medway Massachusetts: 

November 1997, variously paginated. 

DeLima, Virginia, 1991, Stream-aquifer relations and 

yield of stratified-drift aquifers in the Nashua 

River Basin, Massachusetts: U.S. Geological 

Survey Water-Resources Investigations Report 88-

4847, 47 p.

DiBona, Pam, and Eleria, Anna, 1999, Run of the river: 

Auburndate, MA, Streamer, Newsletter of the 

Charles River Watershed Association, Spring 

1999, v. 30, p. 6–7.
References 91

Dickerman, D.C., Trench, E.C.T., and Russell, J.P., 

1990, Hydrogeology, water quality, and groundwater development alternatives in the lower Wood 

River ground-water reservoir: U.S. Geological 

Survey Water-Resources Investigations Report 89-

4031, 109 p. 

Dufresne-Henry, Inc., 1996, Report on Mill River 

Pump Test, Norfolk Water Department, Norfolk, 

Massachusetts: February 1996, variously 

paginated.

_____1997, Report on Gold Street well zone II: 

January 1997, variously paginated.

Dysart, J.E., and Rheaume, S.J., 1999, Induced 

infiltration from the Rockaway River and water 

chemistry in a stratified-drift aquifer at Dover, 

New Jersey: U.S. Geological Survey WaterResources Investigations Report 96-4068, 112 p.

ENSR Corporation, 2000, American National Power, 

Milford, Massachusetts, Charles River Monitoring 

Program, 1999, Data report No. 10: June 2000, 

variously paginated. 

Earth Tech, Inc., 1999, Comprehensive wastewater 

management plan and draft environmental impact 

report phase II document, Holliston, 

Massachusetts: Concord, MA, September 1999, 

variously paginated.

_____2000, Comprehensive wastewater management 

plan and final environmental impact report, Phase 

III Document, Holliston, Massachusetts: Concord, 

MA, January 2000, variously paginated.

_____2001, Comprehensive wastewater management 

plan and environmental impact report, Phase I—

Needs analysis, Assabet Consortium: Concord, 

MA, May 2001, variously paginated. 

Elassal, A.A., and Caruso, V.M., 1983, Digital 

elevation models, in McEwen, R.B., Witmer, R.E., 

and Ramey, B.S., eds., USGS Digital Cartographic 

Data Standards: U.S. Geological Survey Circular 

895-B, 40 p.

Farnsworth, R.K., Thompson, E.L., and Peck, E.L., 

1982, Evaporation atlas for the contiguous 48 

United States: Washington, DC, National Oceanic 

and Atmospheric Administration, NOAA 

Technical Report NWS 33, 26 p.

Fiorentino, J.F., Kennedy, L.E., and Weinstein, M.J., 

2000, Charles River Watershed 1997/1998 water 

quality assessment report: Worcester, MA, 

Massachusetts Department of Environmental 

Protection, Division of Watershed Management, 

Report Number 72-AC-3, variously paginated.

Geophysical Applications, 1995, Seismic refraction 

survey, O’Brien and Hospital Brooks, regional 

storm water mitigation basin project, Milford, 

Massachusetts: June 1995, variously paginated.

Gorelick, S.M., Freeze, R.A., Donahue, David, and 

Keely, J.F., 1993, Groundwater contamination—

Optimal capture and containment: Boca Raton, 

FL, Lewis Publishers, 385 p.

Groundwater Associates, Inc., 1991, Evaluation of the 

8-inch diameter test well (site No. 12, Charles 

River Basin), Bellingham, Massachusetts: June 

1991, variously paginated.

Haley and Ward, Inc., 1994, Interim report on water 

quality at the Populatic Street water supply well 

No. 1—March 1994, Geophysical Applications, 

1995, Seismic refraction survey, O’Brien and 

Hospital Brooks, regional storm water mitigation 

basin project, Milford, Massachusetts: June 1995, 

variously paginated.

_____1995, Report on second phase of water quality 

investigation at the Populatic Street water supply 

well No. 1—September 1995, Geophysical 

Applications, 1995, Seismic refraction survey, 

O’Brien and Hospital Brooks, regional storm 

water mitigation basin project, Milford, 

Massachusetts: June 1995, variously paginated.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and 

McDonald, M.G., 2000, MODFLOW-2000, the 

U.S. Geological Survey modular ground-water 

model—User guide to modularization concepts 

and the ground-water flow process: U.S. 

Geological Survey Open-File Report 00-92, 121 p.

Helsel, D.R., and Hirsch, R.M., 1992, Statistical 

methods in water resources: Amsterdam, The 

Netherlands, Elsevier, 522 p.

Higgins, R., 2000, “H2O aplenty, but ...”: Boston 

Sunday Globe, May 14, 2000, p. 1.

Hill, M.C., 1998, Methods and guidelines for effective 

model calibration: U.S. Geological Survey WaterResources Investigations Report 98-4005, 90 p.

Hill, M.C., Banta, E.R., Harbaugh, A.W., and 

Anderman, E.R., 2000, MODFLOW-2000, the 

U.S. Geological Survey modular ground-water 

model—User guide to the observation, sensitivity, 

and parameter-estimation processes and three 

post-processing programs: U.S. Geological Survey 

Open-File Report 00-184, 209 p.

Hirsch, R.M., 1982, A comparison of four streamflow 

record extension techniques: Water Resources 

Research, v. 18, no. 4, p. 1081–1088.
92 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Johnson, A.I., 1967, Specific yield—Compilation of 

specific yields for various materials: U.S. 

Geological Survey Water-Supply Paper 1662-D, 

74 p.

Kontis, A.L., 2001, Computer program for simulation 

of variable recharge with the U.S. Geological 

Survey modular finite-difference ground-water 

flow model (MODFLOW): U.S. Geological 

Survey Open-File Report 00-173, 75 p.

Kontis, A.L., Randall, A.D., and Mazzaferro, D.L., in 

press, Regional hydrology and simulation of flow 

of stratified drift aquifers in the glaciated 

northeastern United States: U.S. Geological 

Survey Professional Paper 1415C, in press.

Layne New England Company, 1974, Town of 

Medway, Massachusetts, report of evaluation test 

well and seven day pumping test: August 1974, 

variously paginated.

Massachusetts Technology Collaborative, 1998, The I495 regional overview, November 1998, accessed 

March 13, 2001, at URL http://www.mtpc/org/ 

495solutions/495reports.htm.

_____1999, 495 water and sewer task force, January 

1999, accessed March 13, 2001, at URL http:// 

www.mtpc/org/495solutions/495reports.htm.

MassGIS, 1997, MassGIS datalayer descriptions and 

guide to user services: Boston, MA, 

Commonwealth of Massachusetts, Executive 

Office of Environmental Affairs, 131 p.

Masterson, J.P., and Barlow, P.M., 1997, Effects of 

simulated ground-water pumping and recharge on 

ground-water flow in Cape Cod, Martha’s 

Vineyard, and Nantucket Island Basins, 

Massachusetts: U.S. Geological Survey WaterSupply Paper 2447, 79 p.

Masterson, J.P., Walter, D.A., and LeBlanc, D.R., 1998, 

Delineation of contributing areas to selected 

public-supply wells, Western Cape Cod, 

Massachusetts: U.S. Geological Survey WaterResources Investigations Report 98-4237, 45 p.

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

modular three-dimensional finite-difference 

ground-water flow model: U.S. Geological 

Techniques of Water-Resources Investigations, 

book 6, chap. A1, 586 p.

Melvin, R.L., de Lima, Virginia, and Stone, B.D., 

1992, The stratigraphy and hydraulic properties of 

till in southern New England: U.S. Geological 

Survey Open-File Report 91-481, 53 p.

Metcalf & Eddy, 1983a, Report concerning Phase 2 of 

groundwater exploration program, to Mr. Harry 

Sabin, Chairman, Board of Water and Sewer 

Commissioners, Medway, Massachusetts: March 

1983, variously paginated.

_____1983b, Report on the completion Phase 2 of 

groundwater exploration program, to Mr. Gardner 

Rice, Chairman, Board of Water and Sewer 

Commissioners, Medway, Massachusetts: March 

1983, variously paginated.

_____1997, Final report to Milford Water Company, 

Milford, Massachusetts, Water master plan and 

analysis of reservoir yield: March 1997, variously 

paginated.

Metropolitan Area Planning Council, 2001, Regional 

comprehensive development strategy, accessed 

March 13, 2001, at URL http/www.mapc.org.

Mitsch, W.J., and Gosselink, J.G., 1993, Wetlands: 

New York, Van Nostrand Reinhold, (2d ed.), 

722 p.

Moench, A.G., Garabedian, S.P., and LeBlanc, D.R., 

2000, Estimation of hydraulic parameters from an 

unconfined aquifer test conducted in a glacial 

outwash deposit, Cape Cod, Massachusetts: U.S. 

Geological Survey Open-File Report 00-485, 

139 p.

Morris, D.A., and Johnson, A.I., 1967, Summary of 

hydrologic and physical properties of rock and 

soil materials, as analyzed by the hydrologic 

laboratory of the U.S. Geological Survey 1948–60: 

U.S. Geological Survey Water-Supply Paper 1839-

D, 42 p.

Myette, C.F., and Simcox, A.C., 1992, Water resources 

and aquifer yields in the Charles River Basin, 

Massachusetts: U.S. Geological Survey WaterResources Investigations Report 88-4173, 50 p.

National Oceanic and Atmospheric Administration, 

2001, National Virtual Data System, Climate data, 

Daily surface data, accessed May 12, 2001, at 

URL http://www.ncdc.noaa.gov/climate/ 

stationlocator.html.
References 93

Nicholson, R.S., and Watt, M.K., 1997, Simulation of 

ground-water flow in the unconfined aquifer 

system of Toms River, Metedeconk River, and 

Kettle Creek Basins, New Jersey: U.S. Geological 

Survey Water-Resources Investigations Report 97-

4066, 99 p.

Nwankwor, G.I., Cherry, J.A., and Gillham, R.W., 

1984, A comparative study of specific yield 

determinations for a shallow sand aquifer: Ground 

Water, v. 22, no. 6, p. 764–772.

Prince, K.R., Franke, O.L., and Reilly, T.E., 1988, 

Quantitative assessment of the shallow groundwater flow system associated with Connetquot 

Brook, Long Island, New York: U.S. Geological 

Survey Water-Supply Paper 2309, 28 p.

Prudic, D.E., 1989, Documentation of a computer 

program to simulate stream-aquifer relations using 

a modular, finite-difference, ground-water flow 

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

88-729, 113 p.

Randall, A.D., Francis, R.M., Frimpter, M.H., and 

Emery, J.M., 1988, Northeastern Appalachians, in

Back, William, Rosenshein, J.S., and Seaber, P.R., 

Hydrogeology: Boulder, CO, The Geological 

Society of America, The Geology of North 

America, v. O-2, p. 177–187.

Ries, K.G., III, 1993, Estimation of low-flow duration 

discharges in Massachusetts: U.S. Geological 

Survey Open-File Report 93-38, 59 p.

_____1994, Development and application of 

generalized-least-squares regression models to 

estimate low-flow duration discharges in 

Massachusetts: U.S. Geological Survey WaterResources Investigations Report 94-4155, 33 p.

_____1997, August median streamflows in 

Massachusetts: U.S. Geological Survey WaterResources Investigations Report 97-4190, 27 p.

_____1999, Streamflow measurements, basin 

characteristics, and streamflow statistics for lowflow partial record stations operated in 

Massachusetts from 1989 through 1996: U.S. 

Geological Survey Water-Resources Investigations 

Report 99-4006, 162 p.

Ries, K.G., III, and Friesz, P. J., 2000, Methods for 

estimating low-flow statistics for Massachusetts 

streams: U.S. Geological Survey Water-Resources 

Investigations Report 00-4136, 81 p.

Rorabaugh, M.I., 1964, Estimating changes in bank 

storage and ground-water contribution to 

streamflow: International Association of Scientific 

Hydrology, Publication 63, p. 352–363.

Rosenshein, J.S., 1968, Hydrologic characteristics and 

sustained yield of principal ground-water units, 

Potowomut–Wickford Area, Rhode Island: U.S. 

Geological Survey Water-Supply Paper 1775,

38 p.

_____1988, Alluvial valleys, in Back, William, 

Rosenshein, J.S., and Seaber, P.R., Hydrogeology: 

Boulder, CO, The Geological Society of America, 

The Geology of North America, v. O-2, p. 177–

187.

Rutledge, A.T., 1993, Computer programs for 

describing the recession of ground-water discharge 

and for estimating mean ground-water recharge 

and discharge from streamflow records: U.S. 

Geological Survey Water-Resources Investigations 

Report 93-4121, 45 p.

_____1998, Computer programs for describing the 

recession of ground-water discharge and for 

estimating mean ground-water recharge and 

discharge from streamflow records—Update: U.S. 

Geological Survey Water-Resources Investigations 

Report 98-4148, 43 p.

SEA Consultants, Inc., 1999, Town of Bellingham 

comprehensive water resources management plan/ 

Environmental impact report, Phase I 

Environmental Conditions, water supply, 

stormwater, & wastewater needs analysis: 

Cambridge, MA, variously paged.

Schrage, L.E., 1997, Optimization modeling with 

LINDO (5th ed.): Pacific Grove, CA, Brooks/Cole 

Publishing Company, 470 p.

Socolow. R.S., Whitley, J.S., Murino, Dominic, Jr., and 

Ramsbey, L.R., 2001, Water Resources Data 

Massachusetts and Rhode Island: U.S. Geological 

Survey Water-Data Report MA-RI-00-1, 459 p.
94 Simulation of Ground-Water Flow and Evaluation of Water-Management Alternatives in the Upper Charles River Basin, Eastern MA

Solley, W.B., Pierce, R.R., and Perlman, H.A., 1993, 

Estimated use of water in the United States in 

1990: U.S. Geological Survey Circular 1080, 76 p.

Stone, J.R., Shafer, J.P., London, E.H., DiGiacomoCohen, Mary, Lewis, R.S., and Thompson, W.B., 

1998, Quaternary geologic map of Connecticut 

and Long Island Sound Basin: U.S. Geological 

Survey Open-File Report 98-371, 69 p.

Stone, J.R., and Dickerman, D.C., 2002, Surficial 

materials map of the Big River Area, Central 

Rhode Island: U.S. Geological Survey WaterResources Investigations Report 01-4169, 1 pl.

Tasker, G.D., and Driver, N.E., 1988, Nationwide 

regression models for predicting urban runoff 

water quality at unmonitored sites: Water 

Resources Bulletin, v. 24, No. 5, p. 1090–1101.

Volckmann, R.P., 1975a, Surficial geologic map of the 

Holliston Quadrangle, Middlesex, Norfolk, and 

Worcester Counties, Massachusetts: U.S. 

Geological Quadrangle Map GQ-1217, 1 sheet.

_____1975b, Surficial geologic map of the Medfield 

Quadrangle, Middlesex, Norfolk, and Worcester 

Counties, Massachusetts: U.S. Geological 

Quadrangle Map GQ-1218, 1 sheet.

_____1977, Bedrock geologic map of the Holliston 

Quadrangle, Middlesex, Norfolk, and Worcester 

Counties, Massachusetts: U.S. Geological Survey 

Miscellaneous Investigations Map I-1-53, 1 sheet.

Walker, E.H., Caswell, W.W., and Wandle, S.W., Jr., 

1977, Hydrologic data of the Charles River Basin, 

Massachusetts: U.S. Geological Survey 

Hydrologic-Data Report No. 19, 53 p.

Walker, E.H., Wandle, S.W., Jr., and Caswell, W.W., 

1975, Hydrology and water resources of the 

Charles River Basin, Massachusetts: U.S. 

Geological Survey Hydrologic Atlas HA-554, 2 

sheets.

Whitman & Howard, Inc., 1991, Zone II delineations 

for the Godfrey Brook, Clark Island, and Dilla 

Street wellfields, Milford, Massachusetts: 

February 1991, variously paginated.

_____1992, Final report on delineation of Zones II and 

III for municipal wells 2, 3, and 4; water-resources 

evaluation; and aquifer-resource mapping, 

Wrentham, Massachusetts: November 1992, 

variously paginated.

_____1996, Numerical groundwater flow modeling 

and Zone II delineations for wells No. 1, 2, 3, 4, 5 

and 6, Holliston, Massachusetts: March 1996, 

variously paginated.

Williams, J.R., 1967, Drastic lowering of Kingsbury 

Pond, Norfolk, Massachusetts: U.S. Geological 

Survey Open-File Report, July 9, 1967, 33 p.

Zarriello, P.J., and Ries, K.G., III, 2000, A 

precipitation-runoff model for analysis of the 

effects of water withdrawals on streamflow, 

Ipswich River Basin, Massachusetts: U.S. 

Geological Survey Water-Resources Investigations 

Report 00-4029, 99 p.