Prepared in cooperation with the Rhode Island Department of Health

Evaluating Prediction Uncertainty of Areas 

Contributing Recharge to Well Fields of Multiple 

Water Suppliers in the Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island

Scientific Investigations Report 2012–5114

U.S. Department of the Interior

U.S. Geological Survey
THIS PAGE INTENTIONALLY LEFT BLANK
Evaluating Prediction Uncertainty of Areas 

Contributing Recharge to Well Fields of 

Multiple Water Suppliers in the 

Hunt–Annaquatucket–Pettaquamscutt 

River Basins, Rhode Island

By Paul J. Friesz

Prepared in cooperation with the Rhode Island Department of Health

Scientific Investigations Report 2012–5114

U.S. Department of the Interior

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

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

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

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

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

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

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

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

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

U.S. Government.

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

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

Suggested citation:

Friesz, P.J., 2012, Evaluating prediction uncertainty of areas contributing recharge to well fields of multiple water 

suppliers in the Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island: U.S. Geological Survey Scientific 

Investigations Report 2012–5114, 46 p., at http://pubs.usgs.gov/sir/2012/5114/.
iii

Acknowledgments

The author especially thanks Clayton Commons and Richard Amirault (Rhode Island Department 

of Health, Office of Drinking Water Quality) and Ernie Panciera and Theodore Peters (Rhode 

Island Department of Environmental Management, Office of Water Resources) for supporting 

the study and for their suggestions. The author also thanks Jeffrey Starn (U.S. Geological Survey 

(USGS)) for providing a computer program to facilitate the Monte Carlo analysis. Finally, the 

author thanks Jack Eggleston, Martha Nielsen, and Rodney Sheets, all of the USGS, for their 

helpful reviews of the report.
THIS PAGE INTENTIONALLY LEFT BLANK
v

Contents

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

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

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

Description of Study Area and Previous Investigations ................................................................5

Overview of Original Models...............................................................................................................6

Simulation of Groundwater Flow.................................................................................................................9

Model Extent and Spatial Discretization...........................................................................................9

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

Hydraulic Properties...........................................................................................................................14

Groundwater Withdrawals and Return Flow .................................................................................15

Observation Data ................................................................................................................................15

Groundwater Levels ..................................................................................................................15

Base Flow....................................................................................................................................16

Calibration............................................................................................................................................17

Estimation of Model Parameters.............................................................................................17

Simulated Water Table and Water Budget............................................................................27

Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis ..........................27

Deterministic Areas Contributing Recharge ..................................................................................29

Probabilistic Areas Contributing Recharge....................................................................................35

Limitations of Model....................................................................................................................................42

Summary and Conclusions.........................................................................................................................42

References Cited..........................................................................................................................................44

Figures

1. Map showing location of study area, selected U.S. Geological Survey long-term 

network streamgages and observation wells, and a National Oceanic and 

Atmospheric Administration climatological station and tidal gage.....................................3

2. Map showing production wells, model extents, surficial geology, and selected 

drainage boundaries, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island..................................................................................................................................4

3. Map showing selected land uses, Hunt–Annaquatucket–Pettaquamscutt River 

Basins, Rhode Island....................................................................................................................8

4. Map showing section line, model-boundary types, groundwater-level observations, 

and continuous-record and partial-record streamgages, Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island..........................................................................11

5. Model cross section, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island................................................................................................................................13

6. Bar graph showing composite scaled sensitivities for model parameters, Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island............................................21

7. Scatterplots showing relation of A, weighted observation to weighted simulated 

equivalent, and of B, weighted residual to weighted simulated equivalent, Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island............................................23
vi

8. Map showing spatial distribution of weighted residuals and simulated water-table 

contours for calibrated, steady-state conditions, Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island..........................................................................24

9. Scatterplot showing normal probability of the weighted residuals, Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island............................................26

10. Map showing simulated deterministic areas contributing recharge to the 

production wells for 2004–08 average pumping rates, Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island..........................................................................30

11. Map showing simulated deterministic areas contributing recharge to the 

production wells for their maximum pumping rates, Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island..........................................................................31

12. Map showing simulated groundwater traveltimes to the production wells at 

their maximum pumping rates, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island................................................................................................................................34

13. Graphs showing model-parameter distributions A, before acceptance criteria and 

B, after acceptance criteria were applied to determine the probabilistic areas 

contributing recharge to the well centers, Hunt–Annaquatucket–Pettaquamscutt 

River Basins, Rhode Island .......................................................................................................36

14. Map showing simulated probabilistic areas contributing recharge to the well 

centers at their 2004–08 average pumping rates, Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island..........................................................................38

15. Maps showing simulated probabilistic areas contributing recharge to A, well 

center 1 and B, well centers 2, 3, and 4, at their maximum pumping rates, Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island............................................40

Tables

1. Withdrawal rates for production wells in the groundwater-flow model, Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island..............................................7

2. Definition of model parameters and statistics on parameter values, whether 

estimated or specified, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island................................................................................................................................10

3. Summary information on weighted residuals (observed minus simulated values) 

for groundwater-level observations, Hunt–Annaquatucket–Pettaquamscutt River 

Basins, Rhode Island..................................................................................................................16

4. Summary of the long-term mean annual streamflow and base-flow analysis for 

partial-record streamgages, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island................................................................................................................................18

5. Streamflow and drainage-area characteristics of selected continuous-record 

streamgages in Rhode Island and Connecticut.....................................................................20

6. Comparison of observed and simulated base flows, Hunt–Annaquatucket–

Pettaquamscutt Rvier Basins, Rhode Island..........................................................................21

7. Simulated steady-state average annual hydrologic budget for modeled area and for 

major units in the model, Hunt–Annaqutucket–Pettaquamscutt River Basins, 

Rhode Island................................................................................................................................28

8. Summary information concerning deterministic areas contributing recharge by 

well center, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island ...........32
vii

Conversion Factors, Datum, and Abbreviations

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

Area

square mile (mi2) 2.590 square kilometer (km2)

Flow rate

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

gallon per minute (gal/min) 0.06309 liter per second (L/s)

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

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

Hydraulic conductivity

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

Transmissivity*

foot squared per day (ft2/d) 0.09290 meter squared per day (m2/d)

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as follows:

 °C=(°F–32)/1.8

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29).

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

Altitude, as used in this report, refers to distance above the vertical datum.

*Transmissivity: The standard unit for transmissivity is cubic foot per day per square foot times 

foot of aquifer thickness [(ft3/d)/ft2]ft. In this report, the mathematically reduced form, foot 

squared per day (ft2/d), is used for convenience.

Mean annual streamflow and mean annual base flow: Mean annual is defined as the arithmetic 

mean of yearly flows.

Abbreviations

KCWA Kent County Water Authority

NKWD North Kingstown Water Department

RIDEM Rhode Island Department of Environmental Management

RIDOH Rhode Island Department of Health

RIEDC Rhode Island Economic Development Corporation

USGS U.S. Geological Survey
THIS PAGE INTENTIONALLY LEFT BLANK
Evaluating Prediction Uncertainty of Areas Contributing 

Recharge to Well Fields of Multiple Water Suppliers in 

the Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island

By Paul J. Friesz

Abstract

Three river basins in central Rhode Island—the Hunt 

River, the Annaquatucket River, and the Pettaquamscutt 

River—contain 15 production wells clustered in 4 pumping 

centers from which drinking water is withdrawn. These highcapacity production wells, operated by three water suppliers, 

are screened in coarse-grained deposits of glacial origin. The 

risk of contaminating water withdrawn by these well centers 

may be reduced if the areas contributing recharge to the well 

centers are delineated and these areas protected from land uses 

that may affect the water quality. The U.S. Geological Survey, 

in cooperation with the Rhode Island Department of Health, 

began an investigation in 2009 to improve the understanding 

of groundwater flow and delineate areas contributing recharge 

to the well centers as part of an effort to protect the source of 

water to these well centers. A groundwater-flow model was 

calibrated by inverse modeling using nonlinear regression to 

obtain the optimal set of parameter values, which provide a 

single, best representation of the area contributing recharge to 

a well center. Summary statistics from the calibrated model 

were used to evaluate the uncertainty associated with the 

predicted areas contributing recharge to the well centers. This 

uncertainty analysis was done so that the contributing areas to 

the well centers would not be underestimated, thereby leaving the well centers inadequately protected. The analysis led 

to contributing areas expressed as a probability distribution 

(probabilistic contributing areas) that differ from a single or 

deterministic contributing area.

Groundwater flow was simulated in the surficial deposits and the underlying bedrock in the 47-square-mile study 

area. Observations (165 groundwater levels and 7 base flows) 

provided sufficient information to estimate parameters representing recharge and horizontal hydraulic conductivity of the 

glacial deposits and hydraulic conductance of streambeds. 

The calibrated value for recharge to valley-fill deposits was 

27.3 inches per year (in/yr) and to upland till deposits was 

18.7 in/yr. Calibrated values for horizontal hydraulic conductivity of the valley-fill deposits ranged from 20 to 480 feet per 

day (ft/d) and of the upland till deposits was 16.2 ft/d. Calibrated values of streambed hydraulic conductance ranged from 

10,000 to 52,000 feet squared per day. Values of recharge and 

horizontal hydraulic conductivity of the valley-fill deposits 

were the most precisely estimated, whereas the horizontal 

hydraulic conductivity of till deposits was the least precisely 

estimated.

Simulated areas contributing recharge to the well centers 

on the basis of the calibrated model ranged from 0.19 to 

1.12 square miles (mi2

) and covered a total area of 2.79 mi2 

for average well center withdrawal rates during 2004–08 

(235 to 1,858 gallons per minute (gal/min)). Simulated areas 

contributing recharge for the maximum well center pumping 

capacities (800 to 8,500 gal/min) ranged from 0.37 to 3.53 mi2 

and covered a total area of 7.99 mi2 in the modeled area. 

Simulated areas contributing recharge extend upgradient of 

the well centers to upland till and to groundwater divides. 

Some areas contributing recharge include small, isolated areas 

remote from the well centers. Relatively short groundwater 

traveltimes from recharging locations to discharging wells 

indicated the wells are vulnerable to contamination from 

land-surface activities: median traveltimes ranged from 2.9 to 

5.0 years for the well centers, and 78 to 93 percent of the 

traveltimes were 10 years or less for the well centers. Land 

cover in the areas contributing recharge includes a substantial 

amount of urban land use for the two well centers in the Hunt 

River Basin, agriculture and sand and gravel mining uses for 

the well center in the Annaquatucket River Basin, and, for 

the well center in the Pettaquamscutt River Basin, land use is 

primarily undeveloped.

Model-prediction uncertainty was evaluated using a 

Monte Carlo analysis. The parameter variance–covariance 

matrix from nonlinear regression was used to create parameter 

sets that reflect the uncertainty of the parameter estimates and 

the correlation among parameters. The remaining parameters representing the glacial deposits (vertical anisotropy of 

valley-fill deposits and of till deposits, maximum groundwater 

evapotranspiration, and hydraulic conductance for headdependent cells representing a groundwater divide) that could 
2 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

not be estimated with nonlinear regression were incorporated 

into the variance–covariance matrix using prior information 

on parameters. Thus the uncertainty analysis was an outcome 

of calibrating the parameters to available observations and to 

information that the modeler provided. A water budget and 

model-fit statistical criteria were used to assess parameter sets 

so that prediction uncertainty was not overestimated. Because 

of the effects of parameter uncertainty, the size of the probabilistic contributing areas for each well center for both average 

and maximum pumping rates was larger than the size of the 

deterministic contributing areas for the well center. Thus, 

some areas not in the deterministic contributing area may actually be in the contributing area, including additional areas of 

urban and agricultural land use. Generally, areas closest to the 

well centers with short groundwater traveltimes are associated 

with higher probabilities, whereas areas distant from the well 

centers with long groundwater traveltimes are associated with 

lower probabilities. The deterministic contributing areas generally corresponded to areas associated with high probabilities 

(greater than 50 percent). Areas associated with low probabilities extended long distances along groundwater divides in the 

uplands remote from the well centers.

Introduction

Accurate delineation of areas contributing recharge 

to production wells is an essential component of Federal, 

State, and local strategies for the protection of drinkingwater supplies from contamination (U.S. Environmental 

Protection Agency, 1991). The area contributing recharge to 

a well is defined as the surface area where water recharges 

the groundwater and then flows toward and discharges to the 

well (Reilly and Pollock, 1993). At the State level, the Source 

Water Assessment Program of the Rhode Island Department 

of Health (RIDOH), Office of Drinking Water Quality, was 

established as the result of the 1996 Amendments to the 

Federal Safe Drinking Water Act. Since that time, RIDOH has 

assessed the susceptibility and risk of public-water supplies 

to contamination, and the agency has encouraged land-use 

planning within the areas contributing recharge to a 

production well.

The basins for three rivers in central Rhode Island—

the Hunt (H) River, the Annaquatucket (A) River, and the 

Pettaquamscutt (P) River—contain 15 production wells from 

which drinking water is being withdrawn (figs. 1 and 2). The 

town of North Kingstown Water Department (NKWD), Rhode 

Island Economic Development Corporation (RIEDC), and 

Kent County Water Authority (KCWA) withdraw drinking 

water from these large-capacity production wells, which 

supply a total average daily rate of 5.0 million gallons a day 

(Mgal/d) from unconfined sand and gravel aquifers of glacial 

origin. Groundwater withdrawals within this area, referred 

to as “the HAP” in this report, meet water-quality standards 

without need for treatment. Wells in the Hunt River Basin, 

however, are threatened by commercial, industrial, highway, 

and dense residential development (DeSimone, 1999; Hickey 

and Joubert, 2003). Nitrate concentrations are slightly elevated 

in some of these wells, indicating the presence of fertilizers 

or wastewater. In addition, the gasoline additive methyl 

tert-butyl ether (MTBE) has been detected in some of these 

wells (Richard Amirault, Rhode Island Department of Health, 

written commun., 2008). Wells in the Annaquatucket and 

Pettaquamscutt River Basins are in mostly rural areas, but 

30 to 40 percent of the currently protected areas could qualify 

to be converted to residential use. Land use upgradient of the 

wells in the Annaquatucket River Basin includes extensive 

agriculture and sand and gravel mining.

A U.S. Geological Survey (USGS) study in the HAP by 

Barlow and Dickerman (2001) designed a numerical model to 

simulate groundwater flow in the glacial valley-fill deposits; 

this model was updated by Barlow and Ostiguy (2007). One of 

the purposes of these studies was to determine areas contributing recharge to the production wells. Uplands bordering the 

valley-fill deposits were incorporated indirectly into the model 

by adding streamflow and groundwater discharge from these 

uplands at the edge of the model. Simulated areas contributing 

recharge to most production wells in the existing HAP models 

extended to this model boundary, indicating that uplands outside the model may be contributing water to the wells. Land 

cover in the uplands near some of the wells includes urban and 

agricultural land uses. If groundwater flow in the uplands was 

also to be simulated in a new model, the addition would help 

ensure that the entire area that contributes water to the production wells is included.

The original HAP models were calibrated to hydrologic 

data manually to provide a reasonable match between field 

(observed) and simulated groundwater levels and streamflows. 

Manual calibration, however, may not provide the optimal set 

of parameter values that give the best model fit to the observations. Results from a manually calibrated model also do not 

provide a means for quantitatively assessing model-prediction 

uncertainty concerning, in this case, the size, shape, and location of the area contributing recharge to a well. Because the 

simulated contributing areas to the HAP wells are generally 

narrow and long and, in some cases, do not cover the area 

above the well, the RIDOH and the Rhode Island Department 

of Environmental Management (RIDEM) are concerned about 

the practicality of applying these model results. Without an 

evaluation of the uncertainty associated with the model predictions, the contributing area to a well may be underestimated, 

thereby leaving the well inadequately protected.

The USGS, in cooperation with the RIDOH, began a 

2-year study in 2009 to increase understanding of groundwater 

flow in the HAP as part of an effort to protect the source of 

water to these 15 large-capacity production wells. The new 

study modified the original HAP models to incorporate the 

surrounding uplands directly into the model and calibrated the 

model to hydrologic data by inverse modeling using nonlinear 

regression. Nonlinear regression estimates the optimal set of 

model-parameter values. In contrast to parameter values that 
Introduction 3

Base from U.S. Geological Survey,

Rhode Island state plane projection, 1:1000,000, NAD 83

Sites west of Rhode Island are generalized (not to scale)

Pawcatuck

River

Block Island Sound

Queen

Chipuxet River

Usquepaug

River

Wood River

 River

0 10 MILES

0 10 KILOMETERS

42°

71°30' 71°45'

71°15'

41°30'

EXPLANATION

8452660 NOAA tidal gage and identifier

8452660

USGS streamgage and identifier (some sites off map)

USGS observation well and identifier

NOAA climatological station and identifier

41°45'

MASSACHUSETTS

CONNECTICUT

RHODE ISLAND

MASSACHUSETTS

01117500

01117000

01115098

NKW255

01117800

NKW255

01118000 01117468

01118300

01123000

374266

374266

(latitude=41°40'08",

longitude=72°03'10")

CHW18

EXW238

EXW278

WGW206

WCW59

Study area

Hunt River

Annaquatucket

River

River Mattatuxet

River

Pettaquamscutt

Narragansett Bay

01117000

(latitude=41°28'29"

longitude=71°50'05")

Figure 1. Location of study area, selected U.S. Geological Survey (USGS) long-term network streamgages and observation 

wells, and a National Oceanic and Atmospheric Administration (NOAA) climatological station and tidal gage.
4 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

MattatuxetRiver

Chipuxet

River Basin

Frenchtown

71°33' 71°30' 71°27'

41°39'

41°36'

41°33' Well Center 3

Well Center 2

Well Center 4

Well Center 1

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Narragansett Bay

Coventry

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9 NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

Well Center 3

Well Center 2

Well Center 4

Well Center 1

NK4

SFH1

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook

Brook

Tributary

Hunt River

Pettaquamscutt

Sandhill Brook

River

Belleville Pond

Secret 

Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Stratified deposits

Morainal deposits

Till and bedrock outcrops

Water

Area where surface water drains 

southwest and groundwater drains east

Model boundary

Extent of original model (Barlow and Ostiguy, 2007)

Production well and identifier—drinking water

Wetlands Production well and identifier—other

EXPLANATION

Usquepaug Queen

River Basin

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83

0 1 2 KILOMETERS

0 1 2 MILES

Geology

Figure 2. Production wells, model extents, surficial geology, and selected drainage boundaries, Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Introduction 5

are calibrated manually, summary statistics from nonlinear 

regression can be used to provide a quantitative measure of 

the uncertainty of the predicted contributing area (Starn and 

others, 2010). The uncertainty that this study considers arises 

from the observation dataset used for this calibration and not 

from the model design. Results of the calibrated model (used 

for a deterministic contributing area) and the uncertainty 

analysis (used for a probabilistic contributing area) could be 

useful to water-resource managers when they assess the risk 

of contamination and implement land-use plans. Quantitative 

uncertainty analyses of the predicted contributing areas have 

been applied in only a few studies of small modeled areas. 

This study applies the uncertainty analysis in a large model 

with multiple water suppliers and multiple wells.

Purpose and Scope

This report documents the design and calibration of a 

numerical groundwater-flow model for the purpose of delineating the areas contributing recharge to 15 production wells 

supplying drinking water in the HAP River Basins. The 

NKWD operates a total of 11 production wells in the 3 river 

basins. In the Hunt River Basin, the RIEDC operates three 

production wells, and the KCWA operates one production 

well. The original HAP models were modified to simulate 

groundwater flow in the uplands and for use with inverse 

modeling. Changes to the lateral extent increased the modeled 

area from 18.0 to 47.0 square miles (mi2). The groundwaterflow model was calibrated by inverse modeling using nonlinear regression to estimate the optimal set of parameter values. 

Simulated areas contributing recharge to the wells based on 

the optimal set of parameter values (the deterministic contributing area) are shown on maps for selected pumping rates 

and for average, steady-state hydrologic conditions. Summary 

statistics from nonlinear regression were used to evaluate 

the uncertainty associated with the predicted areas contributing recharge to the well fields. Maps depict the results of 

the uncertainty analysis of the simulated area contributing 

recharge expressed as a probability distribution (the probabilistic contributing area).

Description of Study Area and Previous 

Investigations

The HAP study area is in central Rhode Island primarily in the towns of North Kingstown, East Greenwich, and 

Exeter (figs. 1 and 2). The study area is referred to as the HAP 

by this and most previous investigations because the major 

river basins are the Hunt, Annaquatucket, and Pettaquamscutt 

River Basins; the Pettaquamscutt is called the Mattatuxet in 

its headwaters. These three rivers, along with Cocumcossuc 

Brook and other small streams in the eastern part of the study 

area, discharge to Narragansett Bay. In the southwest, parts 

of the Usquepaug–Queen and the Chipuxet River Basins are 

also in the study area (fig. 2). Previous studies, most recently 

Dickerman and Barlow (1997), have shown that groundwater and surface-water divides do not coincide between the 

Annaquatucket River Basin and the Usquepaug–Queen and 

the Chipuxet River Basins. Surface water in these two basins 

flows southwest of the study area, eventually to the Pawcatuck 

River and Block Island Sound, but some of the groundwater in 

these basins flows eastward to the Annaquatucket River Basin. 

The total size of the HAP study area is 47.0 mi2.

Previous investigations by the USGS have analyzed 

the geology and hydrology of part or all of the study area. 

Surficial and bedrock geology have been mapped by Power 

(1957, 1959), Quinn (1952, 1963), Schafer (1961), Smith 

(1955, 1956), and Williams (1964). Reconnaissance studies 

of groundwater conditions, including measurements of 

groundwater levels in valley-fill deposits and in upland 

till, were done by Allen (1956), Allen and others (1959), 

Hahn (1959), and Johnson and Marks (1959). The first 

comprehensive investigation of the groundwater and surfacewater resources was done by Rosenshein and others (1968). 

As part of this study they completed contour maps of the 

bedrock surface and of the water table, saturated thickness, 

and transmissivity of the valley-fill deposits. Information on 

streamflows and streambed characteristics was also collected. 

A water-table map of the valley-fill deposits, based on waterlevel measurements made during October 1996 at nearaverage water-level conditions, was prepared by Dickerman 

and Barlow (1997). A numerical model was developed by 

Barlow and Dickerman (2001) to simulate groundwater flow 

in the valley-fill deposits. This same study used optimization 

techniques to evaluate alternative water-management 

strategies to reduce streamflow depletion and simulate the 

areas contributing recharge to the production wells for 1996 

average withdrawals. Barlow and Ostiguy (2007) modified 

this model to account for new hydrogeologic data and used 

the modified model to assess the effects of existing and 

proposed well withdrawals at a State fish hatchery in the 

Annaquatucket River Basin on groundwater and streamflow 

conditions, especially in nearby wetlands. In addition, the 

study simulated areas contributing recharge to the production 

wells for average withdrawals during 2003, and it included a 

NKWD proposed production well in the Annaquatucket River 

Basin. An overview of these 2001 and 2007 numerical models 

is described in the following section of this report.

In the HAP study area, glacial and post-glacial deposits 

overlie crystalline and metamorphosed sedimentary bedrock 

(fig. 2). These glacial deposits consist of till, stratified, and 

morainal deposits; most of these deposits were laid down 

during advance and retreat of the last (late Wisconsinan) 

ice sheet, which was retreating northward through southern 

Rhode Island between 17,000 to 18,000 radiocarbon years 

ago (Stone and others, 2005). A thin, discontinuous layer 

of till deposited directly on the bedrock by glacial ice is 

composed of a poorly sorted mixture of sediments ranging in 

size from clay to boulders. Stratified deposits consisting of 

well-sorted, layered sediments ranging in size from clay to 

gravel that were deposited by glacial meltwater overlie the 
6 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

till in the valleys, including an upland area of Frenchtown 

Brook Basin. Mixed sediments of till and stratified deposits 

are present in moraines, which formed at the ice margins 

(Schafer, 1961). The stratified deposits and the morainal 

deposits are referred to as valley-fill deposits in this report; 

they cover about 61 percent of the study area (fig. 2). Deposits 

in the valley range from a surface altitude of a few feet along 

the shoreline of Narragansett Bay to 250 feet (ft) along the 

upland-valley contact in the Usquepaug–Queen River Basin. 

Till deposits in the uplands rise to 480 ft in the headwaters 

of Frenchtown Brook. The production wells are screened 

in coarse-grained valley-fill deposits composed of sand and 

gravel. Transmissivity of the valley-fill deposits ranges from 

zero at the upland-valley contacts to about 40,000 feet squared 

per day (ft2/d) (Rosenshein and others, 1968); areas with the 

highest transmissivities (20,000 to 40,000 ft2/d) are in areas 

along most of the northeast-flowing Hunt River, including the 

production wells in the lower Hunt River Basin, and in the 

Annaquatucket River Basin in an area near Secret Lake and 

associated NKWD production wells. Post-glacial deposits of 

alluvium and peat locally overlie glacial deposits in floodplain and wetland areas (fig. 2).

The climate is humid and temperate with an average 

annual temperature of about 50°F and average annual precipitation of 48.7 inches (in.), according to records from 1941 

through 2008 for climatological station 374266 at Kingston, 

Rhode Island (National Oceanic and Atmospheric Administration, 2009) (fig. 1). The source of all water in the study area 

is ultimately from precipitation. Groundwater generally flows 

from topographical highs in the uplands toward streams and 

the valley-fill deposits. The groundwater system is recharged 

by direct infiltration of precipitation, stream leakage, and a 

small amount of wastewater from septic systems (Barlow and 

Dickerman, 2001). Groundwater is decreased by discharges 

to streams and surface-water bodies, by evapotranspiration 

from the water table, and by pumping at the production wells. 

The aquifer is generally in close hydraulic connection with 

the surface-water system (Barlow and Dickerman, 2001). A 

detailed description of the hydrogeology of the study area is 

available in Barlow and Dickerman (2001) and Rosenshein 

and others (1968).

Groundwater is withdrawn by 19 large-capacity production wells (fig. 2, table 1), 15 of which supply drinking water. 

The focus of this study is these 15 wells: 11 production wells 

of the NKWD (NK1–NK11), 3 wells of the RIEDC (3A, 9A, 

and 14A), and 1 well of the KCWA (KC1). The other four 

wells do not supply drinking water: three (wells SFH1–SFH3) 

provide water to a State fish hatchery in the Annaquatucket 

River Basin, and well IW is for an industrial site in the Hunt 

River Basin.

RIDOH and RIDEM delineate areas contributing 

recharge and determine susceptibility to contamination by 

well center. The 15 wells cluster in 4 well centers (fig. 2): well 

centers 1 and 2 (7 wells) are in the Hunt River Basin, well center 3 (5 wells) is in the Annaquatucket River Basin, and well 

center 4 (3 wells) is in the Pettaquamscutt River Basin.

Well center 1 consists of six wells (NK9, NK10, KC1, 

3A, 9A, and 14A) adjacent to the Hunt River in the lower part 

of the basin. No withdrawals were made from NK10 during 1996 or 2003, the years that the original models used to 

simulate areas contributing recharge. Well center 2 has a single 

well (NK6) near Scrabbletown Brook, a tributary to the Hunt 

River. Well center 3 consists of five wells (NK1, NK2, NK4, 

NK5, and NK11) in the Annaquatucket River Basin. Four of 

these wells are near Secret Lake or a tributary of the lake. The 

fifth well (NK11), upgradient of the others, began withdrawals 

in December 2008. Well center 4 consists of three wells (NK3, 

NK7, and NK8) in the Pettaquamscutt River Basin near Carr 

Pond and a tributary to the pond.

One measure for assessing well vulnerability to contamination that the RIDOH uses is the amount of high-intensity 

land use in the contributing area to a well or well center. 

These urban and agricultural land uses store, apply, or generate pollutants that have the potential to contaminate nearby 

water resources (Hickey and Joubert, 2003). Such uses include 

commercial, industrial, high- and medium-high density 

residential (one-fourth acre or less), medium-density residential (one-fourth to one acre), waste disposal, institutional, 

transportation, sand and gravel mining, and selected agriculture (cropland, orchards, and nurseries). The distribution of 

these urban and agricultural land uses (fig. 3) is important in 

assessing the potential sources of contaminants near the well 

centers. The remaining land cover not shown in figure 3 is not 

considered high-intensity land use, including forest, brushland, 

wetland, low density residential, and other types of agriculture 

(for example, pasture). Of the four well centers, well center 

1, which has the highest pumping rate, is in an area with the 

most urban land cover, whereas well center 4, with the lowest 

pumping rate, is in a mostly undeveloped area.

Overview of Original Models

Numerical models of the HAP aquifer were developed 

by Barlow and Dickerman (2001) and Barlow and Ostiguy 

(2007) to simulate groundwater flow in the surficial deposits 

in the valleys (fig. 2). A brief overview of the original steadystate models provides background information for the present 

study’s modifications that more accurately simulate groundwater flow in the uplands bordering the valleys and for the 

subsequent recalibration of the model with nonlinear regression. A detailed description of the original models is in Barlow 

and Dickerman (2001) and Barlow and Ostiguy (2007).

The original models simulated groundwater flow in the 

valley-fill deposits by a four-layered model with a uniformly 

spaced grid; each model cell was 200 by 200 ft. The model 

grid was oriented parallel to the northeast-trending valleys of 

the Hunt River and of Sandhill Brook and to the southeasttrending valleys of the Annaquatucket and Pettaquamscutt 

Rivers. In the vertical dimension, the model extended from 

the water table, defined by Dickerman and Barlow (1997), to 

the bedrock surface, defined by Rosenshein and others (1968) 
Introduction 7

Table 1. Withdrawal rates for production wells in the groundwater-flow model, Hunt–Annaquatucket–Pettaquamscutt River 

Basins, Rhode Island.

[gal/min, gallons per minute; ft3/s, cubic feet per second; KCWA, Kent County Water Authority; RIEDC, Rhode Island Economic Development Corporation; NKWD, North Kingstown Water Department; RIDEM, Rhode Island Department of Environmental Management; --, not applicable]

Water supplier 

well name

Water 

supplier

Well 

center

1996 withdrawal rate 2004–08 withdrawal rate 1Maximum withdrawal rate

(gal/min) (ft3/s) (gal/min) (ft3/s) (gal/min) (ft3/s)

Hunt River Basin

KC1 KCWA 1 226 0.50 563 1.25 1,200 2.67

3A RIEDC 1 155 0.35 135 0.30 1,100 2.45

9A RIEDC 1 113 0.25 142 0.32 1,100 2.45

14A RIEDC 1 272 0.61 171 0.38 1,100 2.45

NK6 NKWD 2 172 0.38 319 0.71 900 2.01

NK9 NKWD 1 819 1.82 470 1.05 2,000 4.46

NK10 NKWD 1 0 0 377 0.84 2,000 4.46

IW Industrial -- 174 0.39 2174 20.39 2174 20.39

Total for basin 1,931 4.30 2,351 5.24 9,574 21.34

Annaquatucket River Basin

NK1 NKWD 3 127 0.28 313 0.70 1,000 2.23

NK2 NKWD 3 95 0.21 95 0.21 700 1.56

NK4 NKWD 3 122 0.27 267 0.59 1,000 2.23

NK5 NKWD 3 343 0.76 237 0.53 750 1.67

NK11 NKWD 3 -- -- 3113 30.25 750 1.67

SFH1 RIDEM -- 438 0.98 4290 40.65 5290 50.65

SFH2 RIDEM -- 438 0.98 4400 40.89 5400 50.89

SFH3 RIDEM -- 0 0 4600 41.34 5600 51.34

Total for basin 1,563 3.48 2,315 5.16 5,490 12.24

Pettaquamscutt River Basin

NK3 NKWD 4 124 0.28 60 0.13 250 0.56

NK7 NKWD 4 13 0.03 95 0.21 325 0.72

NK8 NKWD 4 6 0.01 80 0.18 225 0.50

Total for basin 143 0.32 235 0.52 800 1.78

Total for all basins 3,637 8.10 4,901 10.92 15,864 35.36

1

Maximum-rated capacity of pump.

2

The 1996 average rate was used for Well IW.

3

The 2009 average rate was used for Well NK11.

4

Theodore Peters, Rhode Island Department of Environmental Management, written commun., 2010.

5

The 2004–08 average rate was used for Wells SFH1, SFH2, and SFH3.
8 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

Frenchtown

Well

Center 3

Well Center 2

Well Center 4

Well Center 1

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Mattatuxet

Narragansett Bay

Coventry

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9 NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

NK4

SFH1

Well

Center 3

Well Center 2

Well Center 4

Well Center 1

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well and identifier—drinking water

Production well and identifier—other

Selected land-use types 

High- to medium-high density residential (less than 1/4 acre)

Medium-density residential (1/4 to 1 acre)

Commercial

Industrial

Institutional (for example, hospital, schools)

Transportation

Waste disposal

Sand and gravel mining and quarries

Orchards and nurseries

Cropland

Geology

EXPLANATION

2003−04 Rhode Island Geographic Information System

land-use coverage

Base from U.S. Geological Survey, Rhode Island state plane projection

1:24,000, NAD 83

0 1 2 KILOMETERS

0 1 2 MILES

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook Brook

Tributary

Hunt River

Pettaquamscutt

Sandhill Brook

River

Belleville 

Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Fry

River

Figure 3. Selected land uses, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Simulation of Groundwater Flow 9

and modified by Barlow and Dickerman (2001). Thin surficial 

deposits, mostly along the valley edges (fig. 2), were excluded 

from the active area of the model because simulating thin 

deposits created numerical instabilities. Barlow and Ostiguy 

(2007) increased the lateral extent of the active area of the 

2001 model to include thin valley-fill deposits west and south 

of the State fish hatchery wells, and they resolved numerical 

instabilities by simulating all layers, including the top layer, 

with aquifer transmissivities that did not change with changing 

groundwater levels. This change to the lateral extent increased 

the size of the Barlow and Ostiguy (2007) model from 16.9 to 

18.0 mi2.

Hydrologic processes of recharge, evapotranspiration 

from the water table, well withdrawals, and the interaction 

between the aquifer and surface water (streams and ponds) 

were simulated in the original models. Recharge was from 

direct infiltration of precipitation and from wastewater in areas 

that received water supplies but were not sewered. Inflow 

from the uplands bordering the valley, either by streamflow 

or by groundwater, was applied to cells along the edge of the 

model. The southwest extent of the model was along a groundwater divide between the Annaquatucket River Basin and 

the Chipuxet and Usquepaug–Queen River Basins. Because 

NK11, a proposed well that would not begin withdrawals until 

late 2008, would be closer to this model boundary than were 

existing production wells (fig. 2), Barlow and Ostiguy (2007) 

modified that part of the groundwater divide in transmissive 

materials (valley-fill deposits) from a no-flow boundary to 

a general-head boundary. This type of boundary condition 

allows groundwater to flow across the simulated divide in 

response to changing stress conditions in the model.

The original steady-state models simulated long-term 

average hydrologic conditions represented by the 1941–96 

period. The models were manually calibrated to 22 groundwater levels and 5 streamflows. The groundwater levels had been 

measured on October 8, 1996, at near-average water-level 

conditions. Long-term average streamflow was calculated 

at one continuous-record streamgage and estimated at four 

partial-record streamgages.

Simulation of Groundwater Flow

For this study, groundwater levels and flows were 

simulated in the surficial deposits and the underlying 

bedrock in the HAP study area using a three-dimensional 

finite-difference numerical model code MODFLOW-2000 

(McDonald and Harbaugh, 1988; Harbaugh and others, 

2000). The original models were modified to (1) simulate 

groundwater levels and flows in the uplands (fig. 2) and 

(2) represent boundary conditions and hydraulic properties as 

parameters (table 2) for calibration by nonlinear regression 

and for evaluating model-prediction uncertainty. The 

groundwater-flow model was calibrated to long-term average, 

steady-state hydrologic conditions based on 172 groundwaterlevel and streamflow (base flow) observations.

Model Extent and Spatial Discretization

The geographic extent of the active area of the original models was expanded to include all of the thin surficial 

deposits in the valley and the predominately till uplands 

draining toward the valley (fig. 2). Topographical divides 

in the relatively low-permeability till where groundwater 

and surface-water divides are most likely to coincide were 

used for most of the lateral extent. In addition, the model 

was extended eastward to Narragansett Bay and one of its 

tributaries, Mill Creek—hydrologic features that serve as 

model boundaries (fig. 4). Thickness of surficial deposits near 

Narragansett Bay range from 0 to about 60 ft (Rosenshein 

and others (1968)). Modifying the southwest extent of the 

model along the groundwater divide in transmissive deposits 

separating the Annaquatucket River Basin from the Chipuxet 

and Usquepaug–Queen River Basins was beyond the scope 

of this study. These changes to the lateral extent of the model 

were made to (1) improve representation of the groundwater 

system, (2) move the edge of the model farther from pumping 

centers in most places, and (3) allow for additional streamflow observations in model calibration, or to more accurately 

simulate the drainage area of streamflow observations used in 

the original models. Information for constructing the model 

in these new areas was available from Rosenshein and others 

(1968) and from USGS surficial geological and topographical quadrangles. These changes to the 2007 model extent 

increased the active area from 18.0 to 47.0 mi2. The model 

grid was extended using the same cell size (200 by 200 ft) and 

orientation (43 degrees north of east) as the original models. 

The model grid consisted of 205 rows and 293 columns, and it 

included a total of 130,928 active cells.

Model layers were simulated by using a fixed 

transmissivity, including the top layer, to linearize and thereby 

simplify the numerical calculations and increase numerical 

stability. Barlow and Ostiguy (2007) changed the top layer of 

the model by Barlow and Dickerman (2001) from a variable 

transmissivity to a fixed transmissivity because of numerical 

instability caused by incorporating into the model thin surficial 

deposits near the State fish hatchery wells. Similar numerical 

instability would have occurred when simulating thin layers on 

sides of steeply sloping hills, such as the uplands. Calibration 

by inverse modeling when using nonlinear regression can 

also present convergence difficulties for the model. Another 

advantage of using a fixed transmissivity is that it increases 

the number of model simulations that converge for the 

analysis of the probabilistic contributing area. Simplifying 

the numerical calculations by using a fixed transmissivity 

is described by Hill (1998) and Hill and Tiedeman (2007). 

Barlow and Ostiguy (2007) reported that changes made to the 

model, including the fixed transmissivity for layer 1, did not 

result in substantial changes to simulated groundwater levels, 

streamflow, or hydrologic budgets.

In the original models, the valley-fill deposits were 

subdivided vertically into four model layers that extend from 

the water table to the bedrock surface. The lateral extent of 
10 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I. Table 2. Definition of model parameters and statistics on parameter values, whether estimated or specified, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island. [ft/d, feet per day; in/yr, inches per year; --, dimensionless or not applicable] Parameter name Parameter description Units Optimal or specified value 95-percent confidence interval (Lower value–Upper value) Coefficient of variation (dimensionless) Estimated by nonlinear regression K_MULT Multiplier of horizontal hydraulic conductivity of valley-fill deposits -- 0.81 0.67–1.00 0.10 K_TILL Horizontal hydraulic conductivity of till deposits ft/d 16.2 10.4–25.1 0.23 R_TILL Effective recharge rate on till deposits in/yr 18.7 13.1–24.3 0.15 R_VF Recharge rate on valley-fill deposits in/yr 27.3 22.2–32.4 0.09 SB_MULT Multiplier of streambed hydraulic conductance -- 2.6 2.0–3.5 0.15 Specified by prior information from the literature ETM Maximum evapotranspiration from the water table in valley-fill deposits in/yr 21 13.8–28.3 -- GHB_MULT Multiplier of hydraulic conductance -- 1 0.25–4.0 -- K_ROCK Horizontal hydraulic conductivity of bedrock ft/d 0.1 -- -- K_SW Horizontal hydraulic conductivity of surface water ft/d 50,000 -- -- KV_ROCK Ratio of horizontal to vertical hydraulic conductivity of bedrock -- 1 -- -- KV_TILL Ratio of horizontal to vertical hydraulic conductivity of till deposits -- 5 1–25 -- KV_VF Ratio of horizontal to vertical hydraulic conductivity of valley-fill deposits -- 5 1–25 -- R_SW Effective recharge rate on surface water in/yr 19.5 -- --
Simulation of Groundwater Flow 11

A A'

#

#

#

#

#

#

#

#

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

Frenchtown

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Narragansett Bay

Coventry

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9

NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

01116800

01116820

01116875

01116895

01116980

01117050

01117100

01117200

01117000

A

A'

01116980

01117000

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology Model cells

EXPLANATION

Line of cross section shown in figure 5

Groundwater-level observations

Group A

Group B

Group C

Streamgage—continuous record

Streamgage—partial-record

Stream cell

Constant-head cell

General-head cell

0 1 2 KILOMETERS

Base from U.S. Geological Survey, Rhode Island 0 1 2 MILES

state plane projection, 1:24,000, NAD 83

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook

Brook

Tributary

Hunt River

Pettaquamscutt

Sandhill Brook

River

Belleville Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Mattatuxet River

NK4

SFH1

Well Center 3

Well Center 2

Well

Center 4

Well Center 1

Figure 4. Section line, model-boundary types, groundwater-level observations, and continuous-record and partial-record 

streamgages, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
12 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

the active area decreased from the top layer to the bottom 

layer in conformity with the shape of the bedrock surface. 

For this study the model was relayered to provide a realistic 

representation of groundwater flow from the thin deposits in 

the uplands and near the valley-upland contact to the thicker 

valley deposits. A layer to represent bedrock was also added 

to allow for groundwater flow in bedrock areas, especially 

where surficial deposits are thin, such as beneath the uplands. 

Although the model layering was changed, the total thickness 

of the saturated valley-fill deposits was approximately the 

same in areas coincident to both models.

The model was relayered by first representing surficial 

deposits in one layer (the top layer), and bedrock in a second 

layer (the bottom layer). The top of the model grid, which is 

used to determine transmissivity of the top layer, was initially 

set at land-surface altitude, but then it was reset by running the 

two-layered model, temporarily adjusting hydraulic properties 

until the simulated water table compared favorably to available data (water-table map by Dickerman and Barlow (1997), 

and water-level and land-surface altitudes). This simulated 

water-table altitude was then used as the top of the model 

in subsequent simulations. The top layer was next subdivided into three layers (layers 1 to 3) based on surface-water 

features, lithology, and placements of the production well 

screens. Layer 1 is 5 ft thick to represent surface water in the 

ponds and to simulate shallow groundwater flow near surface 

water accurately. Layers 2 and 3 each represented 50 percent 

of the remaining deposits; the production wells are screened in 

these layers. Stratified, morainal, and till deposits are represented in all three layers. Shallow bedrock areas less than 7 ft 

from the top of the model were incorporated into surrounding 

surficial materials. The bottom layer (now layer 4) represented 

bedrock with a constant thickness of 200 ft throughout the 

model beneath the surficial deposits. An example of the model 

layering from the uplands of Hunt River Basin to Narragansett 

Bay is shown in figure 5.

Boundary Conditions

The model specified the same types of boundary conditions that the original models specified to represent sources of 

recharge and areas of discharge, except for a new boundary 

type to represent Narragansett Bay (fig. 4). Well withdrawals and wastewater return flow are discussed in a subsequent 

section. Recharge was defined on the basis of surficial geology and surface water (table 2). A single recharge parameter 

(R_VF) represented both the stratified deposits and morainal 

deposits; conceptually, because morainal deposits consist 

of permeable materials, recharge rates should be similar for 

both surficial materials. A recharge rate to ponds and lakes, 

defined by parameter R_SW, was not considered for parameter 

estimation but was instead specified for all model simulations as 19.5 inches per year (in/yr) by subtracting the average 

annual evaporation rate of a free-water surface from the average annual precipitation rate (Barlow and Dickerman, 2001). 

In new areas of the model, a third recharge parameter was 

defined for upland till (R_TILL).

Evapotranspiration from the water table in the valley-fill 

deposits was simulated with the evapotranspiration package of MODFLOW. Conceptually most evapotranspiration 

of groundwater in the study area occurs in wetlands and near 

surface-water features. A maximum evapotranspiration rate 

was simulated where the water table was at land surface, and 

this rate decreases linearly to zero where the water table was 

4 ft or more below land surface (Barlow and Dickerman, 

2001). Parameter ETM (table 2) was defined to represent this 

maximum evapotranspiration rate. Evapotranspiration from 

the water table in till deposits was not simulated directly in 

the model because simulated water levels may not be accurate 

enough to accurately simulate evapotranspiration. R_TILL, 

therefore, represents an effective recharge rate, which accounts 

for the effects of groundwater evapotranspiration, whereas 

the recharge parameter (R_VF) for stratified and for morainal 

deposits represents water that reaches the saturated zone, some 

of which may be lost through evapotranspiration. This same 

approach was used by DeSimone (2004) in a similar setting in 

eastern Massachusetts.

Stream-aquifer interactions were simulated as a headdependent flux boundary in layer 1 by using the stream-routing 

package (Prudic, 1989) developed for MODFLOW. The 

stream-routing package accounts for gains and losses of water 

in each stream cell, and it routes streamflow from upstream 

cells to downstream cells. Streams that flow into and out of 

ponds and lakes were simulated as flowing through these 

water bodies. Each stream-routing cell requires a conductance 

term that incorporates the geometry and the vertical hydraulic 

conductivity of the streambed. Field measurements of the 

vertical hydraulic conductivity of streambed sediments in the 

study area by Rosenshein and others (1968) ranged from 0.1 to 

15.2 feet per day (ft/d) for materials ranging from organicrich sediments to coarse sand and gravel. These values were 

used in the original models for different stream reaches based 

on field observations or proximity of the stream to known 

bed sediments. Most streambed conductance values in the 

original models for small streams generally were 4,000 ft2/d, 

for medium streams were 10,000 ft2/d, and for large streams 

were 20,000 ft2/d. These same conductances were used as a 

guide when assigning bed conductances to stream-routing 

cells in new areas of the model. Stream altitudes in new areas 

of the model were determined from, or interpolated between, 

topographical contours that intersected streams listed on 

USGS topographical quadrangles. The new model contained 

2,277 stream-routing cells. A dimensionless parameter 

(SB_MULT, table 2) was defined that multiplied the streambed 

conductance values. Thus a SB_MULT value of 1 would be 

equivalent to the streambed conductances assigned based on 

the original models.

A constant-head boundary was used to simulate 

Narragansett Bay and associated coves. The stage of the 

Bay and the coves was specified as 0.6 ft NGVD29 based 

on a long-term tidal gage in Narragansett Bay, Newport, R.I. 
Simulation of Groundwater Flow 13

Surface waterValley-fill depositsUpland till depositsBedrock

A A'

WEST EAST

Altitude, in feet relative to NGVD 29

VERTICAL SCALE GREATLY EXAGGERATED

EXPLANATION

400

Narragansett Bay

Modellayer1234

-300

0

100

-100-200

300200

0 2 KILOMETERS 10 1 2 MILES

Figure 5. Model cross section, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
14 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

(National Oceanic and Atmospheric Administration, 2010) 

(fig. 1). This altitude represents the stage of the Bay in 1996 

when the most accurate groundwater levels used in nonlinear 

regression were measured. The model contained 306 constanthead cells.

A general-head boundary, which represented the 

groundwater divide in transmissive materials between the 

Annaquatucket River Basin and the Chipuxet and Usquepaug–

Queen River Basins, required a hydraulic conductance 

and groundwater head for each boundary cell. The model 

contained a total of 131 general-head boundary cells in three 

layers along this divide. The hydraulic conductance value was 

calculated in the same manner as that in Barlow and Ostiguy 

(2007). This hydraulic conductance value was determined 

from aquifer transmissivity by Rosenshein and others (1968) 

and from model cell dimensions; values ranged from 250 

to 1,500 ft2/d. A dimensionless parameter (GHB_MULT, 

table 2) was defined that multiplied the hydraulic conductance 

values. The hydraulic conductivity (transmissivity divided 

by saturated thickness) in this area of the model was 50 ft/d. 

A GHB_MULT value of 1 is equivalent to these hydraulic 

conductances or to this hydraulic conductivity value. For 

groundwater head at each boundary cell, Barlow and Ostiguy 

(2007) initially assigned values based on water-table maps by 

Dickerman and Barlow (1997), Hahn (1959), and Johnson and 

Marks (1959) and then, during manual calibration, modified 

the groundwater head specified at each boundary cell so that 

simulated flow across the boundary was close to zero. After 

optimal model-parameter values were estimated by nonlinear 

regression for the present study, this boundary condition was 

activated and groundwater heads were again adjusted so that 

simulated flow across the divide was minimal. Although 

GHB_MULT was not calibrated using nonlinear regression, 

the uncertainty in the specified value was included in the 

prediction uncertainty analysis.

Hydraulic Properties

Hydraulic conductivity parameters were assigned on 

the basis of lithologic units (table 2). Horizontal hydraulic 

conductivity of valley-fill deposits in the original models was 

initially calculated as aquifer transmissivity divided by its 

saturated thickness. Transmissivity contours by Rosenshein 

and others (1968) and as modified by Barlow and Dickerman 

(2001) were drawn on the basis of aquifer-test results and 

lithologic logs. The aquifer-test results provided data on the 

relation between hydraulic conductivity and grain size, and 

the values determined from lithologic logs were based on this 

relation. These horizontal hydraulic values are 20 ft/d for very 

fine sand, 50 ft/d for fine sand, 100 ft/d for sand, 200 ft/d for 

sand and gravel, and 500 ft/d for gravel. Along with additional 

values determined by Dickerman (1984), these values have 

since been used in most Rhode Island groundwater studies. 

Saturated thickness was calculated by subtracting bedrocksurface altitude contours (Rosenshein and others, 1968, 

modified by Barlow and Dickerman, 2001) from water-table 

contours in Dickerman and Barlow (1997). However, Barlow 

and Dickerman (2001) adjusted some hydraulic conductivity 

values during manual calibration because of sparse data in 

some areas and because the transmissivity contours did not 

take into account the geologic framework of stagnation-zone 

retreat of glacial ice and of morphosequence deposition as 

described by Stone and others (2005) for other valley-fill 

settings. The final values from the manually calibrated model 

ranged from 25 to 587 ft/d, and they averaged approximately 

180 ft/d (Barlow and Ostiguy, 2007). This range included 

1,069 different values for horizontal hydraulic conductivity.

Model cells in new areas of the model between the edge 

of the 2007 model and the upland-valley contact that included 

valley-fill deposits were initially assigned horizontal hydraulic 

conductivity values based on values from the 2007 model for 

adjacent model cells. Stratified and morainal deposits in new, 

large areas of the model that are in the Narragansett Bay and 

Frenchtown Brook drainage areas were initially assigned a 

horizontal hydraulic conductivity value of 100 ft/d, which is 

representative of sand. Transmissivity contours by Rosenshein 

and others (1968) in the Narragansett Bay area were based on 

few data, and the upland area of Frenchtown Brook was not 

included in their analysis.

Three different parameterizations of the horizontal 

hydraulic conductivity of the valley-fill deposits were 

evaluated during initial parameter estimation. One of these 

schemes grouped the hydraulic conductivity values into 

three parameter zones, and a second into four parameter 

zones. The range of hydraulic conductivities in each 

parameter zone was guided by the relation between hydraulic 

conductivity and grain size mentioned earlier in this section. 

A third parameterization defined a dimensionless parameter, 

K_MULT, that multiplied the spatially varying hydraulic 

conductivity values. A K_MULT value of 1 would be equal 

to the hydraulic conductivity values in the original models 

and to those values initially assigned for new areas of valleyfill deposits. Model-fit statistics that can be used to compare 

models with alternative parameterizations, AIC and BIC 

(Hill and Tiedeman, 2007), indicated that the model using 

K_MULT had a slightly better model fit than the models using 

the other two parameterization schemes. Also, a comparison 

of weighted residuals indicated that the model using K_MULT 

more randomly distributed weighted residuals than did the 

other two schemes. Finally, the model that used K_MULT was 

more stable than the other two models during initial parameter 

estimation. Because these model-fit statistics and weighted 

residual comparisons indicated that the model with K_MULT 

was more accurate than the other two models, and because 

of its regression stability, the K_MULT model was used in 

subsequent model simulations.

Two lithologic units were not in the original models: 

till and bedrock. Parameter K_TILL represented horizontal hydraulic conductivity of till, and parameter K_ROCK 

represented horizontal hydraulic conductivity of bedrock. A 

parameter representing surface-water bodies, K_SW, which 
Simulation of Groundwater Flow 15

includes ponds, lakes, and Narragansett Bay and its coves, 

was specified for all model simulations. Model cells containing surface-water bodies were assigned a horizontal hydraulic 

conductivity value of 50,000 ft/d to simulate minimal resistance to flow and the corresponding flat gradient across these 

water bodies (Barlow and Dickerman, 2001).

The ratio of horizontal hydraulic conductivity to vertical 

hydraulic conductivity for the surficial deposits was represented by two parameters: KV_VF for valley-fill deposits 

and KV_TILL for upland till deposits. The ratio of horizontal 

hydraulic conductivity to vertical hydraulic conductivity for 

bedrock was defined by parameter KV_ROCK.

Groundwater Withdrawals and Return Flow

Average annual withdrawal rates for 1996 at the production wells were specified in the model using the rates listed 

in table 1. Barlow and Dickerman (2001) determined that the 

total withdrawal rates from each basin for 1996 were representative of total long-term average withdrawals from each basin 

during 1941–96 based on available records. Withdrawals by 

the State fish hatchery wells are discharged to the headwaters 

of the Annaquatucket River after use at the hatchery.

Wastewater return flow to groundwater, such as by septic 

systems, was specified in the model for all model simulations. Wastewater return flow in new areas of the model that 

received water supplies but are not sewered was determined 

in the same manner as in the original areas of the model. The 

amount of wastewater return flow was determined from information on the location and rates of water-supply deliveries 

available for 1996 (Barlow and Dickerman, 2001). The total 

amount of wastewater return flow applied to the model was 

2.1 cubic feet per second (ft3/s), which compared to 1.2 ft3/s in 

the original modeled area.

Observation Data

Parameter values in the groundwater-flow model were 

estimated using two types of observations—groundwater 

levels and base flow (groundwater discharge). Observations 

were weighted on the basis of methods described by Hill and 

Tiedeman (2007); the weights account for measurement error 

and for the difference in units between groundwater levels 

and base flows. More accurate observations are given larger 

weights than less accurate observations and thus have more 

influence in the regression and on the estimated parameter 

values. Different weights also allow for more observations 

to be included in the regression. The model was calibrated to 

172 observations located throughout the modeled area.

Groundwater Levels

Groundwater-level observations used in model calibration were selected from available groundwater levels compiled from previous USGS investigations. A subset (165 

groundwater levels) was chosen to provide a generally areal 

distribution of observations. Near the valley-upland contact, 

which are areas that have large differences in hydraulic properties, a closer spacing was used if observations were available. Also, groundwater-level observations were selected that 

were at, or near, average groundwater-level conditions, except 

in some areas of the model where either none or few groundwater levels were available at average conditions. Groundwater-level conditions were determined by long-term observation 

well NKW255 in the southern part of the study area, or for 

measurements that preceded those for this observation well, 

long-term observation well CHW18 in southern Rhode Island 

(fig. 1). Water levels have been measured monthly at NKW255 

since August 1954 and at CHW18 since October 1946. Both 

observation wells are screened in valley-fill deposits.

Groundwater-level observations were divided into three 

groups (A, B, and C), depending on their estimated accuracy 

(fig. 4). Group A included 33 groundwater levels measured 

October 8, 1996, near average groundwater-level conditions 

(Dickerman and Barlow, 1997). Group A observations were 

in valley-fill deposits screened in model layers 1, 2, and 3. A 

subset (22 groundwater levels) of this group was used in the 

manually calibrated models. Group B and Group C groundwater levels were originally measured to construct groundwater 

maps at a USGS quadrangle scale (Allen, 1956; Allen and 

others, 1959; Hahn, 1959; Johnson and Marks, 1959). Because 

most of these water levels were measured in shallow dug 

wells, they were included in model layer 1. Group B included 

83 groundwater levels measured August 1952 and periods 

of August–November 1954 at, or near, average groundwater 

conditions. These observations were located in valley-fill 

deposits. Group C included 49 groundwater levels measured in 

till or in a thin veneer of stratified deposits in the upland part 

of Frenchtown Brook Basin. Some of these groundwater levels 

were measured during the same near-average groundwaterlevel conditions of August–November 1954 that were used 

for Group B. Group C also included measurements made 

August 1949 and August 1951 at below average water-level 

conditions; most of these observations were in the upland part 

of Fry Brook and Mawney Brook Basins.

Observation weights are equal to the inverse of the sum 

of the individual variances for each type of measurement error 

(Hill and Tiedeman, 2007). The variance is the square of the 

standard deviation. Three types of measurement error were 

determined for each group: (1) water-level measurement error, 

(2) measuring-point altitude error, and (3) uncertainty in average conditions due to seasonal water-level fluctuations.

Water levels were measured with a steel tape or an 

electrical tape to within 0.01 ft for all three groups. To relate 

the water-level error to a normal distribution, the study used a 

95-percent confidence interval (1.96).

Measuring-point altitudes for Group A were surveyed to 

within 0.001 ft, and a 95-percent confidence interval (1.96) 

was also used. Measuring-point altitudes for Group B and 

Group C were based on land-surface altitude determined 

from USGS quadrangles. Land-surface contour intervals on 
16 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

these quadrangles are within 5 ft at the 90-percent confidence 

interval (1.65) (U.S. Geological Survey, 1980). In the uplands, 

where land-surface topography is steeper than in the valley 

and where the exact location of the well may not be known 

(thereby adding more uncertainty to the measuring-point altitude), an accuracy to within 10 ft was used for Group C.

Uncertainty in the groundwater-level measurements to 

represent average steady-state conditions because of seasonal water-level fluctuations was estimated for each group. 

One long-term observation well was used to assess average 

water-level conditions in the study area, a dynamic system 

with numerous hydrologic settings. Usually one water-level 

measurement was made at each observation well, and selected 

water-level observations used in nonlinear regression were 

made over several decades. The standard deviation of this 

measurement error was calculated as the range in the seasonal 

water-level variation divided by 4, which assumes that the 

observation represents the average and that the range represents the 95-percent confidence interval (Hill and Tiedeman, 

2007). Seasonal water-level fluctuations in valley-fill deposits 

typically vary less than in till deposits because of the larger 

storage capacity of valley-fill deposits. For 14 observation 

wells measured monthly from November 1995 to 1996 in the 

valley-fill deposits, water levels fluctuated about 2 to 4 ft, 

with the largest fluctuations at the highest water-level altitudes 

(Barlow and Dickerman, 2001). Hahn (1959) and Johnson and 

Marks (1959) reported seasonal fluctuations ranging from 3 to 

4 ft in the valley-fill deposits. Water levels at long-term well 

NKW255, screened in valley-fill deposits, typically fluctuated 

4 ft annually. For Group A and Group B, screened in valleyfill deposits, a range of 4 ft was assumed for all observations 

resulting in a standard deviation of 1 ft. For till deposits in the 

study area, Allen and others (1959) and Hahn (1959) reported 

seasonal fluctuations of 16 to 17 ft. Four USGS long-term 

network observation wells (EXW238, EXW278, WCW59, and 

WGW206), screened in till deposits in central Rhode Island 

near the study area (fig. 1), typically fluctuated annually from 

2 to 15 ft. For Group C, a range of 16 ft was assumed for all 

observations, or a standard deviation of 4 ft.

The standard deviation of the total error (square root 

of the total variance for each group) ranged from 1.0 ft for 

Group A to 7.3 ft for Group C. Corresponding observation weights ranged from 1.0 ft-2 for Group A to 0.02 ft-2 for 

Group C (table 3).

Base Flow

Long-term mean annual base flow calculated at eight 

partial-record sites and at the Hunt River streamgage 

(01117000) (fig. 4) were used for seven base-flow observations. In the original models, streamflows made at the partialrecord sites were related to streamflow at the same continuous 

streamgage by application of a graphical-correlation technique 

developed by Searcy (1959) to estimate mean annual streamflow at the partial-record site. For this study, the partial-record 

sites were related to multiple streamgages by a mathematical 

method termed Maintenance of Variance Extension, Type 1 

(MOVE.1) developed by Hirsch (1982) to calculate both mean 

annual streamflow and base flow. Base flow, or groundwater 

discharge, is a measure of effective recharge (recharge minus 

groundwater evapotranspiration).

Streamflow measurements made at the 8 partial-record 

sites were related to concurrent mean daily streamflow at several unregulated continuous streamgages in Rhode Island and 

Connecticut (fig. 1) to estimate the selected statistics. Hunt 

River streamgage (01117000) was not used in the analysis to 

determine selected streamflow statistics at the partial-record 

sites because of groundwater withdrawals upstream of the 

gage. In addition to instantaneous streamflow measurements 

at Annaquatucket River at Belleville (01117100), mean daily 

streamflow from continuous records (September 1961 through 

1964) were also available and used in the analysis. The drainage areas for the continuous streamgages have topographical, 

geologic, and climatic conditions similar to those of the drainage areas of the partial-record sites. Plots of log-transformed 

data were made to determine the quality and linearity of the 

relation between the partial-record sites and the continuous 

Table 3. Summary information on weighted residuals (observed minus simulated values) for groundwater-level observations, 

Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.

[--, not applicable]

Groundwater-level 

observations, by 

group (figure 4)

Number of 

observations

Weight

(feet-2)

Average

(feet)

Minimum

(feet)

Maximum

(feet)

Standard 

deviation

(feet)

Percentage within 

two standard errors 

of the regression

A 33 1 -0.97 -9.26 5.20 3.65 91

B 83 0.1 0.19 -8.28 10.45 3.02 92

C 49 0.02 0.11 -4.15 3.85 1.74 100

All 165 -- -0.07 -9.26 10.45 2.87 94
Simulation of Groundwater Flow 17

streamgages. Streamgages with the highest correlation coefficient (two to three gages) were used to estimate the selected 

streamflow statistics at the partial-record sites. The MOVE.1 

was used to provide an equation that related streamflow at the 

partial-record site to that at the continuous streamgage. Mean 

annual streamflow and mean annual base flow at the continuous streamgages for complete years of record were entered 

into the equation to estimate the corresponding mean annual 

streamflow and base flow at the partial-record sites. Base flow 

at the continuous streamgages was calculated by use of the 

automated hydrograph-separation technique PART (Rutledge, 

1998). The associated mean square error for each relation was 

used to combine the multiple estimates for each partial-record 

site into weighted-average estimates of the selected streamflow statistics to obtain the single best estimate. A detailed 

description of the MOVE.1 technique for streamflow analysis 

is described in Ries and Friesz (2000). Computer programs by 

Granato (2009) facilitated the streamflow analysis.

Information concerning the analysis for the partial-record 

sites is summarized in table 4 and that for the continuous 

streamgages in table 5. In addition to mean annual streamflow 

and mean annual base flow estimated for each partial-record 

site, table 4 lists the equivalent rate over the drainage area for 

six of the eight partial-record sites where the surface-water 

and groundwater divides are most likely to be similar. Mean 

annual base flow ranged from about 14.8 to 27.7 in/yr over the 

drainage areas, with a median of 22 in/yr. The lowest value 

was for Fry Brook (01116895), which drains predominately 

till deposits (98 percent). Base flow at all eight partial-record 

sites ranged from 62 to 92 percent of total streamflow with a 

median of 79 percent. The difference between total streamflow 

and base flow is overland runoff. The lowest percentage was 

for the till-dominated drainage area of Fry Brook, and the 

highest percentage was for Annaquatucket River (01117100). 

The Annaquatucket River Basin is covered by 90-percent 

valley-fill deposits, and its surface-water drainage area and its 

groundwater contributing area are not the same.

Mean annual base flow at six of the eight partial-record 

sites was used directly as observations in model calibration 

(table 6). Mean annual base flow at Scrabbletown Brook 

(01116800) and at Tributary to Hunt River (01116820) were 

not used as observations themselves but were used to compute 

an observation that represented a net gain in base flow along 

the main stem of the Hunt River. Mean annual base flow at the 

partial-record sites was estimated from multiple measurements 

made over several years and related to multiple long-term 

streamgages. For the purpose of weighting these base-flow 

observations in parameter estimation, these mean annual base 

flows were assumed to be accurate within 10 percent at the 

95-percent confidence interval.

The observation that represented a net gain in base 

flow along the main stem of the Hunt River (Hunt River 

observation, 14.1 ft3/s, table 6) was calculated by subtracting 

the total mean annual base flow estimated at the five partialrecord sites (gages 01116800–1116980) in the Hunt River 

Basin (23.9 ft3/s) (fig. 4 and table 4) from a mean annual base 

flow calculated from streamflow records at the Hunt River 

streamgage (01117000; 38.0 ft3/s) (fig. 4, table 5). Streamflow 

records at the Hunt River streamgage, although they may have 

been affected by groundwater withdrawals upstream of the 

streamgage, were used to determine mean annual base flow 

with the hydrograph-separation technique PART (Rutledge, 

1998). The Hunt River base-flow observation of 14.1 ft3/s is 

equivalent to 28.0 in/yr over 6.97 mi2 of mostly valley-fill 

deposits. The accuracy of the mean annual base flow at the 

Hunt River streamgage was assumed to be 5 percent; the total 

variance of the Hunt River observation was calculated by 

adding the variances of the five partial-record sites and the 

streamgage. The coefficient of variation of the Hunt River 

observation (standard deviation divided by the mean) from 

which the observation weight was calculated was 0.08 

(8 percent), which compares to a coefficient of variation of 

0.05 (5 percent) for each of the observations determined at the 

partial-record sites.

Calibration

The groundwater-flow model was calibrated with the 

Parameter Estimation (inverse modeling) process of MODFLOW-2000 (Hill and others, 2000; Hill and Tiedeman, 2007) 

using nonlinear regression that minimizes the differences, or 

residuals, between field (observed) and simulated water levels 

and base flows to obtain an optimal set of parameter values. 

The quality of this calibration was determined by analysis 

of the residuals and the accuracy of the estimated parameter 

values. Some parameters, however, may be insensitive to the 

available observations, or some parameters may be highly correlated with each other and therefore cannot be estimated by 

nonlinear regression. Values from the original models and the 

literature (prior information) were used to specify parameter 

values that could not be estimated by nonlinear regression.

Estimation of Model Parameters

Ten model parameters were evaluated with parameter 

estimation: two for recharge (R_VF, R_TILL), three for 

horizontal hydraulic conductivity (K_MULT, K_TILL, 

K_ROCK ), three for vertical anisotropy ( KV_VF, KV_TILL, 

KV_ROCK), one for streambed hydraulic conductance 

(SB_MULT), and one for the maximum evapotranspiration 

from the water table in valley-fill deposits (ETM) (table 2). 

As mentioned previously, parameter GHB_MULT was not 

considered for parameter estimation, but the uncertainty in 

the parameter value was included in the uncertainty analysis. 

Parameter sensitivities, shown by their composite scaled 

sensitivities in figure 6, indicate whether groundwater-level 

and base-flow observations provided sufficient information to 

permit an estimate of a given parameter (Hill and Tiedeman, 

2007). Parameters with higher sensitivities generally can be 

more precisely estimated than can parameters with lower 

sensitivities. Parameters with composite scaled sensitivities 
18 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I. Table 4. Summary of the long-term mean annual streamflow and base-flow analysis for partial-record streamgages, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.—Continued [USGS streamgage number: Locations shown on figure 4. mi2, square mile; ft3/s, cubic feet per second; in/yr, inches per year; Base-flow index is mean annual base flow divided by mean annual streamflow; --, not applicable] USGS stream- gage number Streamgage name Drain- age area (mi2

)

Area of 

glacial 

stratified 

deposits 

(percent)

Mean annual streamflow Mean annual base flow

Baseflow 

index 

(percent)

Continuous-record streamgage used in relation 90-percent confidence limits 90-percent confidence limits

Streamflow 

(ft3/s)

Lower 

(ft3/s)

Upper 

(ft3/s)

1Streamflow

(in/yr)

Base 

flow 

(ft3/s)

Lower 

(ft3/s)

Upper 

(ft3/s)

1Base 

flow 

(in/yr)

Number 

of measurements 

used in 

relation

USGS 

streamgage 

number

Streamgage name

Correlation 

coefficient

01116800 Scrabbletown Brook 

near Davisville, RI

1.32 24 2.72 1.92 3.86 27.9 2.03 1.43 2.88 20.9 0.75 16 01115098 Peeptoad Brook near North 

Scituate, RI

0.96

25 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.86

25 01123000 Little River near Hanover, CT 0.86

01116820 Tributary to Hunt 

River near Davisville, RI

1.26 20 3.15 2.15 4.61 33.9 2.57 1.76 3.76 27.7 0.82 17 01115098 Peeptoad Brook near North 

Scituate, RI

0.86

17 01117800 Wood River near Arcadia, RI 0.9026 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.85

01116875 Frenchtown Brook at 

State Route 2 near Davisville, RI

6.74 45 12.6 9.77 16.2 25.3 9.69 7.51 12.5 19.5 0.77 17 01115098 Peeptoad Brook near North 

Scituate, RI

0.94

17 01117800 Wood River near Arcadia, RI 0.9617 01123000 Little River near Hanover, CT 0.94

01116895 Fry Brook near East 

Greenwich, RI

3.11 2 5.47 3.03 9.87 23.8 3.39 1.88 6.12 14.8 0.62 16 01115098 Peeptoad Brook near North 

Scituate, RI

0.92

16 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.90

16 01123000 Little River near Hanover, CT 0.86

01116980 Sandhill Brook near 

East Greenwich, RI

3.53 296 7.92 4.44 14.1 30.4 6.19 3.47 11.0 23.8 0.78 17 01118000 Wood River at Hope Valley, 

RI

0.71

17 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.81
Simulation of Groundwater Flow 19

Table 4. Summary of the long-term mean annual streamflow and base-flow analysis for partial-record streamgages, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island.—Continued

[USGS streamgage number: Locations shown on figure 4. mi2, square mile; ft3/s, cubic feet per second; in/yr, inches per year; Base-flow index is mean annual base flow divided by mean annual streamflow; 

--, not applicable]

USGS 

streamgage 

number

Streamgage name

Drainage 

area

(mi2)

Area of 

glacial 

stratified 

deposits 

(percent)

Mean annual streamflow Mean annual base flow

Baseflow 

index 

(percent)

Continuous-record streamgage used in relation 90-percent confidence limits 90-percent confidence limits

Streamflow 

(ft3/s)

Lower 

(ft3/s)

Upper 

(ft3/s)

1Streamflow

(in/yr)

Base 

flow 

(ft3/s)

Lower 

(ft3/s)

Upper 

(ft3/s)

1Base 

flow 

(in/yr)

Number 

of measurements 

used in 

relation

USGS 

streamgage 

number

Streamgage name

Correlation 

coefficient

01117050 Cocumcussoc Brook 

near Wickford, RI

1.52 258 3.42 2.33 5.01 30.5 2.74 1.87 4.02 24.4 0.80 17 01115098 Peeptoad Brook near North 

Scituate, RI

0.88

17 01117800 Wood River near Arcadia, RI 0.8817 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.90

01117100 Annaquatucket River 

at Belleville, RI

6.44 290 20.5 16.0 26.3 -- 18.9 14.7 24.3 -- 0.92 1,154 01117500 Pawcatuck River at Wood 

River Junction, RI

0.71

1,154 01118000 Wood River at Hope Valley, 

RI

0.72

1,154 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.71

01117200 Mattatuxet River near 

Saunderstown, RI

4.81 262 11.1 7.07 17.4 -- 9.33 5.94 14.7 -- 0.84 17 01117468 Beaver River near Usquepaug, RI 0.71

24 01118300 Pendleton Hill Brook near 

Clarks Falls, CT

0.67

24 01123000 Little River near Hanover, CT 0.61

1Assumes drainage area is approximately the same as groundwater contributing area. Not applicable to Annaquatucket and Mattatuxet River basins.2Includes morainal deposits.
20 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I. Table 5. Streamflow and drainage-area characteristics of selected continuous-record streamgages in Rhode Island and Connecticut. [USGS streamgage number: Locations shown on figure 1. mi2, square mile; ft3/s, cubic feet per second; in/yr, inches per year. Base-flow index is mean annual base flow divided by mean annual streamflow] USGS streamgage number Streamgage name Drainage area (mi2

)

Area of glacial 

stratified 

deposits 

(percent)

Period of record 

analyzed

Mean 

annual 

streamflow 

(ft3/s)

Mean 

annual 

streamflow 

(in/yr)

Mean annual 

base flow 

PART 

(ft3/s)

Mean annual 

base flow 

PART 

(in/yr)

Baseflow 

index 

(percent)

01115098 Peeptoad Brook near North Scituate, RI 4.96 24 1995–07 10.4 28.5 7.79 21.3 7501117000 Hunt River near East Greenwich, RI 22.9 52 1941–07 47.3 28.1 38 22.5 8001117468 Beaver River near Usquepaug, RI 9.18 25 1975–2007 21.4 31.6 19.0 28.1 8901117500 Pawcatuck River at Wood River Junction, RI 99.3 47 1941–07 197.0 27.0 177.0 24.2 9001117800 Wood River near Arcadia, RI 35.2 23 1965–80; 1983–07 76.9 29.7 67.0 25.9 8701118000 Wood River at Hope Valley, RI 73.5 26 1942–07 156.7 29.0 133.2 24.6 8501118300 Pendleton Hill Brook near Clarks Falls, CT 3.99 8.4 1959–06 8.72 29.7 6.46 22.0 7401123000 Little River near Hanover, CT 30.0 18 1952–06 56.2 25.4 40.0 18.1 71
Simulation of Groundwater Flow 21

Table 6. Comparison of observed and simulated base flows, Hunt–Annaquatucket–Pettaquamscutt Rvier Basins, Rhode Island.

[USGS streamgage number: Locations shown on figure 4. ft3/s, cubic feet per second; --, not applicable]

USGS 

streamgage 

number

Streamgage or observation name

Observed 

base flow 

(ft3

/s)

Simulated 

base flow 

(ft3/s)

Observed minus 

simulated base flow 

(ft3/s)

Weighted 

residual

01116875 Frenchtown Brook at State Route 2 near Davisville, RI -9.69 -10.0 0.31 0.65

01116895 Fry Brook near East Greenwich, RI -3.39 -3.40 0.01 0.11

01116980 Sandhill Brook near East Greenwich, RI -6.19 -5.36 -0.83 -2.62

-- Hunt River net gain -14.1 -13.5 -0.60 -0.52

01117050 Cocumcussoc Brook near Wickford, RI -2.74 -3.47 0.73 5.24

01117100 Annaquatucket River at Belleville, RI -18.9 -9.39 -9.51 -9.86

01117200 Mattatuxet River near Saunderstown, RI -9.33 -8.69 -0.64 -1.35

Parameter name

R_VF

K_MULT

R_TILL

SB_MULT

K_TILL

ETM

K_ROCK

KV_VF

KV_ROCK

KV_TILL

Composite scaled sensitivity

5

4

3

2

1

0

Figure 6. Composite scaled sensitivities for model parameters (parameter information is provided on table 2), Hunt–

Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
22 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

that are about two orders of magnitude lower than that of the 

parameter with the highest value, or those with composite 

scaled sensitivities less than one, indicate that nonlinear 

regression may not be capable of estimating the parameter 

(Hill and Tiedeman, 2007).

Low sensitivities were associated with five parameters: 

KV_VF, KV_TILL, K_ROCK, KV_ROCK, and ETM (fig. 6). 

Including any of these five parameters in nonlinear regression resulted in the model not converging with reasonable 

values. Thus, these five parameters were either assigned values 

used in the original models or, if a parameter representing a 

geologic unit or hydrologic process was not in the original 

models, values from the literature. Both parameters representing the ratios of horizontal to vertical hydraulic conductivity 

of glacial deposits (KV_VF and KV_TILL) had low sensitivities. The original model value, 5 (table 2), for vertical anisotropy was used for KV_VF. This value was assigned based on 

stratified deposits of sand and gravel in Rhode Island and Cape 

Cod, Massachusetts (Dickerman and others, 1990; Masterson 

and Barlow, 1997; Barlow, 1997). Vertical anisotropy values 

in these referenced studies ranged from 1:1 to 50:1, with 

the highest values for fine sand and silt. Melvin and others 

(1992) summarized the hydraulic properties of till, which can 

be highly variable, from previous studies in southern New 

England. For till derived from crystalline bedrock, horizontal 

hydraulic conductivities ranged from 0.004 to 65 ft/d, and 

vertical hydraulic conductivities ranged from 0.013 to 96 ft/d. 

The present study used the same vertical anisotropic value of 5 

for KV_TILL that was also used for KV_VF.

Low sensitivities were also associated with bedrock 

parameters K_ROCK and KV_ROCK. Hydraulic conductivity 

of crystalline bedrock is generally low. Analysis of specificcapacity data from bedrock wells in eastern Connecticut 

indicated an average hydraulic conductivity of 0.5 ft/d 

(Randall and others, 1966). Lower values of 0.02 and 0.09 ft/d 

for crystalline bedrock in northern New Hampshire were 

determined through model calibration (Tiedeman and others, 

1997). Parameter K_ROCK was specified a value of 0.1 ft/d, 

and KV_ROCK was specified a ratio of 1 (table 2).

The final low-sensitivity parameter was for maximum 

evapotranspiration (ETM) from the simulated water table in 

valley-fill deposits. The original model value of 21.0 in/yr 

was used (table 2); this value is equal to the estimated mean 

growing-season rate of free-water-surface evaporation from 

shallow lakes (Farnsworth and others, 1982).

The five remaining parameters representing recharge 

(R_VF and R_TILL), horizontal hydraulic conductivity 

(K_MULT and K_TILL), and streambed hydraulic conductance (SB_MULT) were not highly correlated (parameter 

correlation coefficients greater than 0.95). Thus these five 

parameters were estimated by nonlinear regression.

The quality of model calibration can be determined 

by analysis of the weighted residuals, both numerically and 

graphically, and by the reasonableness of optimal parameter 

values and their associated confidence intervals; Hill and 

Tiedeman (2007) describe this analysis in detail. Weighted 

residuals should be randomly distributed and close to zero. 

The average weighted residual was -0.11 ft for all groundwater-level and all base-flow observations; it ranged from a minimum of -9.86 ft to a maximum of 10.45 ft. The sum of squared 

weighted residuals was 1,488 for the calibrated model. The 

calculated error variance (sum of squared weighted residuals 

divided by the difference between the number of observations 

and the number of parameters estimated by nonlinear regression) was 8.91, and the standard error of the regression (square 

root of the calculated error variance) was 2.98. Although these 

measures of the overall magnitude of the weighted residuals 

should, theoretically, equal 1, that is not commonly the case 

for groundwater models (Hill and Tiedeman, 2007).

Summary information concerning weighted residuals for 

groundwater-level observations is listed in table 3 for each 

head group category and for all heads combined. The average weighted head residual for all groundwater levels was 

close to zero (-0.07 ft), ranging from a minimum of -9.26 ft to 

a maximum of 10.45 ft. Ninety-four percent of all weighted 

head residuals were within two standard errors of the regression. Group A, the highest weighted head group in the regression, had the largest average weighted residual and the largest 

standard deviation, but 91 percent of these weighted residuals 

were within two standard errors of the regression.

A comparison of observed and simulated base flow 

is listed in table 6. Values for six of seven simulated base 

flows compared well with the observed values: five of seven 

weighted residuals are within one standard error of the 

regression, and six of seven are within two standard errors of 

the regression. For the observation of Annaquatucket River 

streamgage (01117100) at Belleville, R.I., however, simulated 

base flow is one-half the observed flow. This partial-record site 

is downstream of two large ponds, and its watershed does not 

coincide with its groundwater-contributing area, making it difficult for streamflow data to be well correlated with a similar 

long-term streamgage.

Weighted residuals for base flows and for groundwater 

levels are shown graphically in figure 7, which symbolizes the 

weighted head residuals by the category of their groundwaterlevel observation group. A comparison of weighted observed 

values and weighted simulated values (fig. 7A) indicated 

a good agreement; the correlation between them was 0.99. 

Figure 7B shows that weighted residuals are generally 

randomly distributed around zero for all weighted simulated 

values and that most weighted residuals are within two 

standard errors of the regression.

The spatial distribution of weighted head residuals is 

shown in figure 8. The spatial distribution of weighted baseflow residuals is not shown in the figure, but they are distributed randomly throughout the model. The spatial distribution 

of weighted head residuals is generally random in most areas 

of the model, except for clusters of weighted residuals that are 

either all positive or all negative in two areas.

The central part of the model along Sandhill Brook 

and the headwaters of Cocumcussoc Brook has a cluster of 

negative head residuals (simulated water level greater than 
Simulation of Groundwater Flow 23

Weighted simulated equivalent

-30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170

11.92

8.94

5.96

2.98

0

-2.98

-5.96

-8.94

-11.92

Weighted residual, (2.98 equals standard

error of the regression)

Groundwater-level observation

-30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

Weighted observation

180

170

160

150

140

130

120

110

100

90

80

70

60

50

40

30

20

10

0

-10

-20

-30

A

B

Group A

Group B

Group C

Streamflow observation

EXPLANATION

Figure 7. Relation of A, weighted observation to weighted simulated equivalent, and of B, weighted residual to weighted simulated 

equivalent, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
24 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

360

340

320

Frenchtown

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

0 1 2 KILOMETERS

0 1 2 MILES

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Narragansett Bay

Coventry

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9 NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

Weighted groundwater-level residual

-9.26 to less than -4.00

-4.00 to less than -2.00

-2.00 to less than 0.00

Greater than 0.00 to 2.00

Greater than 2.00 to 4.00

Greater than 4.00 to 10.45

Cook's D greater than critical value

20

20

40

40 60

10080

120 140 180160

200

220

240

20

40

40

60

60

60

60

60

80

100120

160

140

180

200

20

20

40

60

80

100

120

140

20

160

120

100

80

140160

180

200

220

220

240

260

280

300

320

340

360

380

400

400

400

260

260

280

300

300

320

340

360

380

40

EXPLANATION

400

380

300

280

260

160

140

120

140

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Simulated water table—

contour interval is 20 feet

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook

Brook

Tributary

Hunt River

Pettaquamscutt

Sandhill Brook

River

Belleville Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Mattatuxet River

Well Center 3

Well Center 2

Well Center 4

Well Center 1

NK4

SFH1

Well Center 3

Well Center 2

Well Center 4

Well Center 1

Figure 8. Spatial distribution of weighted residuals and simulated water-table contours for calibrated, steady-state conditions, 

Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Simulation of Groundwater Flow 25

observed water level), some of which are the largest in the 

model. Areas of the model between Sandhill Brook and 

Hunt River have thin surficial deposits, including some 

exposed till deposits and bedrock in the central part of the 

valley. The lower half of the area drained by Sandhill Brook 

contains high-density commercial and residential development 

where impervious surfaces may reduce recharge. The headwater area of Sandhill and Cocumcussoc Brooks consists of 

extensive wetlands. This complex surface-water feature is difficult to simulate accurately in a groundwater-flow model. In 

some models, wetlands in Rhode Island have been simulated 

with a high hydraulic conductivity to flatten the hydraulic gradient across the wetland (Dickerman and others, 1997; Friesz, 

2004; Masterson and others, 2007). Finally, stream reaches are 

more extensive in these wetlands based on a comparison of 

1:5,000 stream extent to the 1:24,000 used in the model; some 

of these stream reaches drain areas near the highest negative 

weighted residuals.

The southwest part of the model, roughly between 

Queens Fort Brook and Belleville Pond, has a cluster of 

positive head residuals (simulated water level less than 

observed water level). Although these observations are not in 

the immediate vicinity of the State fish hatchery wells, most 

were collected when pumping rates at these wells were less 

than that used in the calibrated model. These observations 

are also in a complex area of the model where loss from 

Queens Fort Brook flows toward the headwater reaches of 

Annaquatucket River.

These two areas indicate that the weighted residuals are 

not entirely randomly distributed and that there is some bias in 

model fit. Simplification of the aquifer and processes in these 

areas of the model may be the source of this model bias.

Insights into observations most important to the parameter estimates were calculated from influence statistics Cook’s 

D and DFBetas using the RESAN-2000 program (Hill and 

others, 2000). Observations with the most overall influence in 

the regression and on the resulting set of estimated parameters 

were determined from Cook’s D statistics. Twenty-three observations (6 base flows and 17 water levels), or 13 percent of the 

observations, had a Cook’s D value greater than the critical 

value of 0.023 (4 divided by the number of observations). Six 

of seven base-flow observations were considered influential 

to the set of estimated parameter values, thus indicating the 

importance of streamflow in this model calibration. In addition, two of the seven base-flow observations (Frenchtown 

Brook, 01116875, and Cocumcussoc Brook, 01117050) (fig. 4) 

were the only observations important to all of the individual 

estimated parameters, according to DFBetas statistics. 

Although the base-flow observations important to the set of 

estimated parameters are distributed throughout the model, 

head observations that dominate the regression are mostly in 

the southern half of the model (fig. 8). In addition, some of 

these head observations have relatively large residuals, both 

negative and positive. These observations with large residuals 

are in the headwaters of several streams previously mentioned 

and near the valley-upland contact in areas with steep hydraulic gradients.

Optimal values for the five recharge and hydraulic 

parameters (R_VF, R_TILL, K_TILL, K_MULT, SB_MULT) 

that were estimated using nonlinear regression are within 

a plausible range of values (table 2). Parameter R_VF, 

which represents precipitation recharge that reaches the 

water table, was 27.3 in/yr. A mean annual recharge rate of 

25.5 in/yr (Barlow and Ostiguy, 2007) was calculated by 

analysis of streamflow records of the Hunt River streamgage 

(01117000) during 1941–03 using the computer program 

RORA (Rutledge, 1993). This recharge rate represents an 

average over the entire Hunt River Basin, including areas 

of stratified deposits, till, wetlands, surface-water bodies, 

and a variety of land uses. Conceptually, recharge rates for 

stratified deposits is higher than this basin-wide average rate 

and lower for less permeable till. Effective recharge rate to 

till deposits, represented by R_TILL, was 18.7 in/yr, close to 

mean annual base flow of 22.5 in/yr estimated by PART from 

the streamflow records at the Hunt River streamgage (table 5). 

This effective recharge rate for till was, as expected, less than 

this basin-wide average rate. The optimal value for R_TILL 

was also within the range of effective recharge rates from 

till-dominated basins in southern New England, as calculated 

by computerized hydrograph-separation techniques from 

long-term streamflow records: effective recharge rates ranged 

from 16 to 24 in/yr when mean annual runoff ranged from 

27 to 31 in/yr (Bent, 1995, 1999; Friesz and Stone, 2007). The 

horizontal hydraulic conductivity of till (K_TILL), 16.2 ft/d, 

although within the range of plausible values summarized by 

Melvin and others (1992) and mentioned in the “Hydraulic 

Properties” section of this report, was higher than values 

of 1 to 10 ft/d that have been used in manually calibrated 

models (DeSimone, 2004; Masterson and others, 2007; Friesz 

and Stone, 2007) and values of 2 to 8 ft/d estimated from 

parameter estimation (Friesz, 2010) for till in Rhode Island 

and Massachusetts. Few data were available to define the 

contact between till and the bedrock surface in the uplands by 

Rosenshein and others (1968). Bedrock-surface altitude also 

varied substantially between data points in the uplands when 

assigning model-cell values. Optimal values for the remaining 

hydraulic parameters (0.81 for K_MULT and 2.6 for 

SB_MULT) were considered reasonable based on values 

reported by Rosenshein and others (1968) and because of the 

number of terms required to define streambed geometry.

The calibration results of the parameterized model and 

the original models are not directly comparable because 

of design changes and the difference in the number of 

observations and type of observations (base flow versus 

streamflow). However, the R_VF value of 27.3 in/yr in the 

parameterized model is relatively close to the manually 

calibrated value of 28.0 in/yr in the original models, as is the 

K_MULT of 81 percent of the original model values. 

SB_MULT, another parameter estimated by nonlinear 

regression and representing a process common to both 

models, was 2.6 times the values in the original models.
26 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

The uncertainty of the parameter estimate is indicated 

by the 95-percent linear confidence interval for each optimal 

value (table 2). For these linear confidence intervals to be 

valid, weighted residuals should be normally distributed and 

the model linear near the estimated optimal values (Hill and 

Tiedeman, 2007). If weighted residuals are independent and 

normally distributed, they plot on an approximately straight 

line on a normal probability graph (fig. 9). The correlation 

between weighted residuals and the normal order statistics 

for the calibrated model was 0.963. This value is near but 

less than the critical value for 172 observations, 0.984, at the 

5-percent significance level. The degree of model linearity can 

be quantified using the modified Beale’s measure, calculated 

with the BEALE-2000 program (Hill and others, 2000). The 

model is considered effectively linear if the modified Beale’s 

measure is less than 0.039 and nonlinear if it is greater than 

0.44. The modified Beale’s measure for the model was 0.46, 

indicating that the model is nonlinear. The confidence intervals 

listed in table 2 are thus approximate values.

The 95-percent confidence intervals for the parameter estimates are all within the ranges of reasonable values reported in the literature. A comparison of the relative 

precision of different parameter estimates can be made using 

the coefficient of variation (standard deviation of the estimated value divided by the optimal value) (table 2); a smaller 

coefficient of variance indicates a more precisely estimated 

value for the parameter. The coefficient of variations ranged 

from 0.09 to 0.23. Recharge parameter R_VF and horizontal 

hydraulic conductivity parameter K_MULT, with coefficient 

of variation values of 0.1 or less, were the most precisely 

estimated, whereas the horizontal hydraulic conductivity 

parameter K_TILL, with a coefficient of variation of 0.23, was 

the least precisely estimated. The order of the most to least 

precisely estimated parameter values follows the same order 

as that of the parameter sensitivities (fig. 6) because of the 

information provided by the observations in the regression.

The analysis of the weighted residuals and optimal 

parameter values indicated that the groundwater model is 

acceptable for the purposes of the study. Optimal parameter 

values are realistic, and their confidence intervals include reasonable values. Although the spatial distribution of weighted 

head residuals was not entirely random, model-fit statistics 

indicated that simulated values for groundwater levels and for 

base flow are generally close to observed values.

Weighted residual

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

Standard normally distributed number 3

2

1

0

-1

-2

-3

Groundwater-level observation

Group A

Group B

Group C

Streamflow observation

EXPLANATION

Figure 9. Normal probability of the weighted residuals, Hunt–Annaquatucket–Pettaquamscutt River Basins, 

Rhode Island.
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 27

Simulated Water Table and Water Budget

The altitude and configuration of the simulated water 

table for the calibrated model, shown in figure 8 at 20-ft 

contour intervals, are consistent with the conceptual model of 

groundwater flow in the study area and with regional groundwater maps by Dickerman and Barlow (1997) and Rosenshein 

and others (1962), as well as with groundwater maps on a 

quadrangle scale by Allen and others (1959), Hahn (1959), 

and Johnson and Marks (1959). Groundwater generally flows 

from topographically high areas and discharges to streams 

and surface-water bodies. In the uplands, the simulated water 

table approximately parallels the land surface, and simulated groundwater divides generally coincide with watershed 

divides. In some upland areas, however, especially in the 

Mattatuxet River Basin, simulated groundwater levels are 

above the land surface but the direction of groundwater flow is 

realistic. The water-table gradient is steepest in the till uplands 

and in valley-fill deposits near the contact and then it flattens 

in the more transmissive valley-fill areas.

The simulated groundwater budget for the calibrated 

modeled area, summarized in table 7, indicated that direct 

precipitation recharge accounts for most of the total inflow 

(77 percent). Streamflow loss accounts for most of the 

remaining inflow (21 percent); most of this streamflow 

loss occurs from tributaries downstream of the uplandvalley contact, which are near areas of abrupt changes in 

transmissivity, and at the downstream ends of ponds and lakes. 

Of the total inflow, 80 percent of the groundwater discharges 

to streams, and another 5 percent discharges directly to 

Narragansett Bay. Evapotranspiration from the water table 

and well withdrawals constitute about 8 percent each of the 

outflow. About one-half of the total well withdrawals are 

either returned to the modeled area by wastewater return flow 

or, in the case of the fish hatchery wells, to the headwaters of 

the Annaquatucket River; the remaining well withdrawals are 

exported from the modeled area.

In addition to the simulated groundwater budget for 

the whole model, water budgets for four major units of the 

modeled area (table 7) provide further insight. To simplify 

the water-budget calculations, the stratified deposits in the 

uplands west of Hunt River are included in the upland till 

category. Recharge is the dominant inflow to the uplands 

(87 percent). Groundwater discharge in the uplands is mostly 

to streams (52 percent), but a significant amount is to valleyfill deposits (38 percent). For the valley-fill deposits, recharge 

is also the main source of inflow (60 percent), although it is 

less of a percentage of the total inflow compared to that in 

the uplands. Upland sources, either directly or indirectly, also 

contribute a significant amount of inflow to valley-fill deposits 

(about one-third). These upland sources include shallow, 

lateral groundwater flow from till (15 percent); deep, vertical 

flow from bedrock (2 percent); and a large percentage of the 

streamflow loss (17 percent) to the valley-fill deposits.

Simulation of Areas Contributing 

Recharge and Prediction Uncertainty 

Analysis

Calibration of the HAP groundwater-flow model by 

inverse modeling using nonlinear regression provided an 

optimal set of parameter values. This optimal parameter set 

was estimated by minimizing the weighted residuals between 

the observation dataset (165 groundwater levels and 7 base 

flows) and simulated values. A predicted area contributing 

recharge to a production well in the HAP based on this optimal 

parameter set provides a single, most likely contributing area 

(deterministic contributing area). However, the parameter 

values were estimated with different levels of uncertainty; this 

uncertainty in the optimal values was based on the information 

that the observation dataset provided on the parameters. 

Parameter uncertainty and its associated effects on model 

predictions (spatial variability of the simulated contributing 

area to a well) can be evaluated by a stochastic Monte Carlo 

analysis. The parameter variance–covariance matrix from 

nonlinear regression can be used to create plausible parameter 

sets for the Monte Carlo analysis (Starn and others, 2000). 

The parameter variance–covariance matrix incorporates the 

uncertainty of the parameter estimates and the correlation 

among parameters from the calibrated model. The Monte 

Carlo analysis was done by replacing the parameter set in 

the calibrated model by a plausible parameter set multiple 

times. The probability of a particular location being in the 

contributing area to a production well was computed from 

these multiple model simulations.

A Monte Carlo analysis of this type whereby summary 

statistics from nonlinear regression were used to create parameter sets was described and applied by Starn and others (2000) 

for a well field screened in valley-fill deposits in a small modeled area. Lindsey (2005) applied the method for contributing 

areas to bedrock wells for those model parameters that could 

be estimated by the observations. The uncertainty analysis 

was also applied to wells in wetland and in coastal settings 

(Friesz, 2010). Starn and others (2010) compared the uncertainty of contributing areas to a well, depending on whether a 

model that used observations of only groundwater levels and 

streamflows or whether the model also included observations 

of atmospheric tracers. In both these referenced 2010 reports, 

prior information on model parameters was used to incorporate inestimable parameters into the uncertainty analysis, and 

a model-fit statistic was used to assess parameter sets so that 

prediction uncertainty would not be overestimated.

The present study delineated the areas contributing 

recharge for the 15 production wells that supply drinking 

water based on their 2004–08 average withdrawal rate and 

their maximum pumping rate (table 1). The total average 

annual withdrawal rate for the 15 production wells that supply drinking water increased from 2,587 gallons per minute 

(gal/min) (5.8 ft3/s) for 1996 to 3,437 gal/min (7.7 ft3/s) for 
28 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

Table 7. Simulated steady-state average annual hydrologic budget for modeled area and for major units in the model, 

Hunt–Annaqutucket–Pettaquamscutt River Basins, Rhode Island.

[ft3/s, cubic feet per second; Inflow and outflow may not equal and percent may not sum to 100 percent because of budget error and rounding]

Hydrologic budget component Flow rate 

(ft3

/s) Percent

Modeled area

Inflow

Recharge from precipitation 82.4 77

Recharge from wastewater return 2.1 2

Streamflow loss 22.5 21

Groundwater divides 0.0 0

Total inflow 107.0

Outflow

Streamflow 85.2 80

Narragansett Bay 5.6 5

Groundwater evapotranspiration 8.2 8

Well withdrawal 8.1 8

Groundwater divides 0.0 0

Total outflow 107.1

Upland deposits (mostly till)

Inflow

All recharge (precipitation and wastewater) 30.8 87

Streamflow loss 1.2 3

From valley-fill deposits 2.0 6

From bedrock 1.4 4

Total inflow 35.4

Outflow

Streamflow 18.4 52

Groundwater evapotranspiration 1.5 4

To valley-fill deposits 13.3 38

To bedrock 2.3 6

Total outflow 35.5

Valley-fill deposits

Inflow

All recharge (precipitation and wastewater) 52.9 60

Streamflow loss 15.3 17

Groundwater divides 0.0 0

From till 13.3 15

From ponds 5.3 6

From bedrock 1.3 2

Total inflow 88.1

Hydrologic budget component Flow rate 

(ft3

/s) Percent

Valley-fill deposits—Continued

Outflow

Streamflow 54.7 62

Narragansett Bay 5.6 6

Groundwater evapotranspiration 6.7 8

Well withdrawal 8.1 9

Groundwater divides 0.0 0

To till 2.0 2

To ponds 10.6 12

To bedrock 0.5 1

Total outflow 88.2

Freshwater ponds

Inflow

Recharge from precipitation 0.8 5

Streamflow loss 6.0 34

From valley-fill deposits 10.6 61

Total inflow 17.4

Outflow

Streamflow 12.0 69

To valley-fill deposits 5.3 31

Total outflow 17.3

Bedrock

Inflow

From valley-fill deposits 0.5 18

From till 2.3 82

Total inflow 2.8

Outflow

To valley-fill deposits 1.3 48

To till 1.4 52

Total outflow 2.7
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 29

the 2004–08 period, although individual wells may show an 

increase or decrease in withdrawal rate. Because there were no 

withdrawals from NK10 in 1996 or in 2003, areas contributing recharge to this well were not simulated in the original 

models. For the 2004–08 period, withdrawals from this well 

averaged 377 gal/min. For NK11, the average pumping rate 

listed represented the average rate for 2009 because this well 

did not begin pumping until December 2008.

The total maximum pumping rate for the production 

wells that supply drinking water, 14,400 gal/min (32.1 ft3/s), 

is 4.2 times greater than the 2004-08 average pumping rate of 

3,437 (7.7 ft3/s). The maximum pumping rate for the production wells represents the maximum-rated capacity of the well 

pumps. Maximum pumping rates used in the model simulation 

are not proposed, long-term (continuous) withdrawal rates. 

Instead, RIDOH considers the simulated areas contributing 

recharge at these maximum rates only when implementing 

land-use planning that protects the quality of the water that 

the production wells supply. Average withdrawals may change 

because of changes in water usage or changes in State policies. 

Areas contributing recharge for the maximum pumping rate 

therefore represent conservative, or larger, areas for land-use 

planning than if low pumping rates were to be used.

The rate of wastewater return was not changed from that 

used in the calibrated model for either of these pumping rate 

scenarios. Although 2004–08 average withdrawal rates were 

greater than those of the 1996 withdrawals, wastewater return 

was a small part of the total inflow to the groundwater system 

(table 7). Maximum pumping rates are much larger than the 

1996 withdrawal rates, but these maximum rates are hypothetical, and the distribution of the wastewater return in the 

model and the percent exported from the model are unknown. 

Thus the area contributing recharge to the production wells, 

especially at the maximum pumping rates, represents a conservative or a larger area than if wastewater return had been 

modified.

Deterministic Areas Contributing Recharge

Simulated deterministic areas contributing recharge 

and groundwater traveltimes to the production wells were 

determined on the basis of the calibrated steady-state model, 

for simulated pumping conditions, and by use of the particletracking program MODPATH (Pollock, 1994). The particletracking program calculates groundwater-flow paths and 

traveltimes on the basis of the head distribution computed by 

the groundwater-flow simulation. Areas contributing recharge 

were delineated by forward tracking of particles from recharge 

locations to the discharging wells. Consistent with the original 

models, particles were stopped at weak sinks, which remove 

only a part of the water that flows into the model cell. Barlow 

and Dickerman (2001) reported that using either option (of 

being stopped by or passing through a weak sink) minimally 

affected the simulated contributing areas to the majority of 

wells. Four particles for each model cell were used, thus each 

particle represents a surface area of 100 by 100 ft.

Shapes of the individual areas contributing recharge are 

strongly affected by nearby pumping wells and hydrogeologic 

features for both pumping rate scenarios (figs. 10 and 11). 

Most simulated areas contributing recharge for both average and maximum withdrawal rates extend upgradient from 

the wells to upland till and to groundwater and topographical 

divides, some of which serve as model boundaries. Some areas 

contributing recharge include small, isolated areas remote 

from the well. In addition, some of the contributing areas do 

not overlie a well. At the average pumping rate, four wells 

(NK5, NK7, NK8, and 3A) are not overlain by their contributing area, and two wells (NK3 and NK11) are overlain by 

the contributing area to a downgradient well (fig. 10). At the 

maximum pumping rate, although all wells are covered by a 

contributing area, three wells (NK3, NK10, and KC1) have 

the contributing area to a nearby well overlying it (fig. 11). 

Wells are screened in either layer 2 or layer 3, and recharge 

travels along paths above and around the screened interval of 

the well.

RIDOH and RIDEM assess the vulnerability of drinking 

water to contamination by well center. The size of the areas 

contributing recharge for the four well centers for the average 

pumping rate ranged from 0.19 mi2 for well center 4 to 

1.12 mi2 for well center 1 and covered a total area of 2.79 mi2

in the model (table 8). The size of the areas contributing 

recharge for the maximum pumping rate ranged from 0.37 mi2

for well center 4 to 3.53 mi2 for well center 3 and covered a 

total area of 7.99 mi2 in the model. Land cover in the areas 

contributing recharge to the well centers in the Hunt River 

Basin include substantial amount of urban land use (fig. 3). 

Land cover in the areas contributing recharge include agriculture and sand and gravel mining for the well center in the 

Annaquatucket River Basin, and, for the well center in the Pettaquamscutt River Basin, land cover is primarily undeveloped.

The size of the area contributing recharge to a well for a 

particular pumping rate is related to effective recharge rates 

from precipitation and, if applicable, from the quantity of 

water derived from other sources. These sources may include 

wastewater return flow and surface-water infiltration. Surfacewater infiltration may be from natural leakage, induced by 

pumping of a well, or both. For example, for the average 

pumping rate, well NK6, with the fourth largest withdrawal 

rate (319 gal/min) had the largest area contributing recharge 

(0.43 mi2), in part, because all of its water was derived from 

precipitation recharge and wastewater return flow. Three 

wells in well center 1 on the lower part of the Hunt River 

(NK9, NK10, and KC1), which have the three largest average 

pumping rates (377 to 563 gal/min), had smaller contributing 

areas than NK6 (0.24 to 0.39 mi2) because part of their water 

was derived from the Hunt River, and for KC1, part of its 

water was also derived from an unnamed tributary to the Hunt 

River that drains the uplands. Well 3A, which has an average 

pumping rate (135 gal/min or 0.30 ft3/s) greater than five other 

wells, had the smallest contributing area in the model. This 
30 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

Frenchtown

Well Center 3

Well Center 2

Well Center 4

Well Center 1

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Coventry

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9 NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

EXPLANATION

NK4

SFH1

Well Center 3

Well Center 2

Well Center 4

Well Center 1

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Well Center 1

NK9

NK10

3A

9A

14A

KC1

Well Center 2

NK6

Well Center 3

NK1

NK2

NK4

NK5

NK11

Well Center 4

NK3

NK7

NK8

Simulated deterministic areas contributing recharge

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83 0 1 2 KILOMETERS

0 1 2 MILES

Narragansett Bay

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook

Brook

Tributary

Hunt River

Pettaquamscutt

Sandhill Brook

River

Belleville 

Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Mattatuxet River

Figure 10. Simulated deterministic areas contributing recharge to the production wells for 2004–08 average 

pumping rates, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 31

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

Frenchtown

Well Center 3

Well Center 2

Well Center 4

Well Center 1

EXPLANATION

NK4

SFH1

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Narragansett Bay

Brook

Coventry

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9 NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

Well Center 3

Well Center 2

Well Center 4

Well Center 1

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83 0 1 2 KILOMETERS

0 1 2 MILES

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Well Center 1

NK9

NK10

3A

9A

14A

KC1

Well Center 2

NK6

Well Center 3

NK1

NK2

NK4

NK5

NK11

Well Center 4

NK3

NK7

NK8

Simulated deterministic areas contributing recharge

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook

Tributary

Pettaquamscutt

Sandhill Brook

River

Belleville 

Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Hunt River

Mattatuxet River

Figure 11. Simulated deterministic areas contributing recharge to the production wells for their maximum 

pumping rates, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
32 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I. Table 8. Summary information concerning deterministic areas contributing recharge by well center, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island. [gal/min, gallons per minute; mi2, square mile; yr, year] Well center Total pumping rate (gal/min) Size of area contributing recharge (mi2

)

Groundwater traveltimes

For 2004–08 average rate For maximum rate

2004–08 

average rate Maximum rate For 2004–08 average rate For maximum rate Median (yr) Percent 10 years or less Median (yr) Percent 10 years or less

1: (NK9, NK10, KC1, 3A, 9A, 14A) 1,858 8,500 1.12 3.03 4.7 78 5.0 782: (NK6) 319 900 0.43 1.06 4.8 90 4.6 873: (NK1, NK2, NK4, NK5, NK11) 1,025 4,200 1.05 3.53 5.3 70 4.3 784: (NK3, NK7, NK8) 235 800 0.19 0.37 4.1 90 2.9 93All well centers 3,437 14,400 2.79 7.99 4.9 77 4.5 80
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 33

well is near the confluence of Frenchtown Brook and the Hunt 

River where most of its water was derived from streamflow 

loss from Frenchtown Brook downstream of the upland-valley 

contact. The location of an area contributing recharge to a well 

can be strongly affected by nearby pumping and associated 

withdrawal rates. For example, NK9 and NK10 have the two 

largest maximum pumping rates (2,000 gal/min each). NK9 

had the largest contributing area in the model for the maximum rates, and its shape and sources of water were similar 

to its contributing area at the average rate. NK10, however, 

had the smallest contributing area in the model for maximum 

rates, even smaller than its contributing area at its lower 

pumping rate. Most of its pumped water was derived from the 

Hunt River (about 90 percent), and its remaining water was 

intercepted groundwater discharge in the valley-fill deposits 

southeast of the well. In contrast, at the average pumping rate, 

the well also intercepted precipitation recharge originating in 

the upland till northwest of the well.

At the maximum pumping rate, the area contributing 

recharge for the well centers expanded in all directions to 

capture enough water to balance the increased pumping rate, 

including downvalley from the wells. For two of the well 

centers (3 and 4) the contributing areas extended downvalley beneath ponds (Belleville and Carr), and additional areas 

extended beneath a lake (Secret Lake) and on the opposite side 

of these surface-water bodies to include small, isolated areas 

remote from the well centers (fig. 11). For NK2 this included 

an area in upland till and valley-fill deposits near the uplandvalley contact northeast of Belleville Pond. Particle tracks 

showed that recharging water originating in this till and in 

valley-fill deposits travels along deep groundwater-flow paths 

in the valley-fill deposits and, under pumping conditions, 

passes beneath Belleville Pond to NK2. Recharging water 

between this contributing area and Belleville Pond travels 

along shallow and intermediate-depth flow paths before it 

discharges to the pond.

The areas contributing recharge to the Annaquatucket 

River Basin wells (well center 3) extend westward from 

the wells to the edge of the model defined by topographical 

divides in upland till and a groundwater divide in transmissive 

materials represented by the general-head boundary (figs. 10 

and 11). In addition to recharge, wastewater return flow, and 

surface-water infiltration from the modeled area, the well center derived part of its water from induced groundwater from 

the Chipuxet River Basin. As previously mentioned in the 

“Boundary Conditions” section, simulated groundwater flow 

across this boundary was close to zero for the pumping rates 

used in the calibrated model (1996 average pumping rates). 

Although the average rates for 2004–08 for this well center 

increased by 338 gal/min from 687 to 1,025 gal/min, including 

withdrawals from the new well NK11 closest to the groundwater divide, the quantity of groundwater induced from the 

adjoining basin was minimal (about 1 percent each of NK11 

and NK5 withdrawals). At the maximum pumping rate for the 

well center (4,200 gal/min), the area contributing recharge not 

only expanded to capture enough water to balance the pumping rate, but also induced enough groundwater to supply 4 percent each of NK4 and NK5 and 16 percent of NK11 withdrawals. This is a total of 190 gal/min or 0.43 ft3

/s. An additional 

0.23 ft3/s of groundwater was induced by pumping along the 

easternmost general-head boundary cells, but this groundwater discharged to the headwaters of the Mattatuxet River. 

The induced groundwater withdrawn by the wells would be 

equivalent to 0.21 mi2 of land surface if the optimal value 

for precipitation recharge for valley-fill deposits (27.3 in/yr) 

was the only source of water and there was no groundwater 

evapotranspiration. This additional land-surface area in the 

Chipuxet River Basin, if adjacent to the divide, would include 

high-intensity land use that is mostly agricultural (cropland).

Simulated traveltime estimates from recharging locations to the production wells for the maximum pumping rate 

are shown in figure 12. Porosity values were specified in 

the model for MODPATH, but they affect only groundwater velocity and do not change the contributing areas to the 

wells. As in the original models, the valley-fill deposits were 

assigned a porosity of 0.35 based on values determined for 

similar deposits in southern Rhode Island (Allen and others, 

1963) and Cape Cod, Massachusetts (Garabedian and others, 

1991). Also, as in the original models, a porosity of 1 was 

specified for surface water. A porosity of 0.35 was specified 

for till deposits based on a range of values (0.22–0.50) determined for till deposits in southern Rhode Island (Allen and 

others, 1963) and southern New England (Melvin and others, 

1992). For bedrock, a porosity of 0.02 was assigned based on 

values for crystalline rock summarized in Meinzer (1923).

Traveltimes generally depend on where recharge enters 

the aquifer in relation to the production well. Water that 

recharges the aquifer near the wells has the shortest traveltimes and youngest water, whereas water originating in the 

till uplands, and for well center 1, also along the groundwater 

divide between the Hunt River and Sandhill Brook, had the 

longest traveltimes and the oldest water. Traveltimes ranged 

from less than 1 year to more than 350 years at each well 

center. Median traveltimes for the well centers ranged from 

2.9 to 5.0 years (table 8). A comparison of median traveltimes 

between maximum and average pumping rates indicated that 

increased pumping caused the percentage of young water 

to increase for three of the four well centers (table 8). Areas 

contributing water to the wells where traveltimes are 10 years 

or less ranged from 78 to 93 percent for the maximum pumping rate and from 70 to 90 percent for the average pumping 

rate for the four well centers. These percentages for traveltimes 10 years or less, and relatively short median traveltimes, 

indicate that the wells are vulnerable to contamination from 

activities on the land surface.
34 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

Frenchtown

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83

0 1 2 KILOMETERS

0 1 2 MILES

Well Center 3

Well Center 2

Well Center 4

Well Center 1

NK4

SFH1

Well Center 3

Well Center 2

Well Center 4

Well Center 1

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook Brook

Tributary

Pettaquamscutt

Sandhill Brook

River

Belleville 

Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Hunt River

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Mattatuxet

Narragansett Bay

Coventry

River

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK9 NK10

NK11

KC1

3A

9A

14A

IW

SFH1

SFH2

SFH3

Simulated traveltimes

Greater than 0 to 1 year

Greater than 1 to 3 years

Greater than 3 to 10 years

Greater than 10 to 25 years

Greater than 25 to 50 years

Greater than 50 years

EXPLANATION

Figure 12. Simulated groundwater traveltimes to the production wells at their maximum pumping rates, 

Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 35

Probabilistic Areas Contributing Recharge

A quantitative measure of the effects of parameter uncertainty on model predictions (the predicted contributing area) 

was done by a Monte Carlo analysis. A Monte Carlo analysis 

was used to obtain the probability of a recharge location being 

in the contributing area of a well center. The probability distribution is related to the information that the observation dataset 

provided on the estimated parameters, to prior information on 

specified parameters, and to the sensitivity of the simulated 

contributing area to the parameters. Hundreds of parameter 

sets generated from summary statistics of the calibrated 

model were used to run hundreds of model simulations in the 

Monte Carlo analysis. Because combinations of reasonable 

parameter values may result in unrealistic groundwater levels 

and streamflows, the parameter sets were evaluated using the 

pumping rates and associated observation dataset from the 

calibrated model. Those parameter sets that simulated realistic 

results were then used in Monte Carlo analyses for the two 

pumping scenarios.

Parameter values for the Monte Carlo analysis were created by the following equation (Starn and others, 2000, 2010):

b z = + σ μ (1)

where

b is a set of parameter values,

z is a vector of normally distributed random 

numbers with a mean of 0 and standard 

deviation of 1,

	 σ	 is the square root of the variance–covariance 

matrix calculated using Cholesky 

decomposition, and

μ is a vector of optimal parameter values.

Parameter values that could not be estimated by nonlinear 

regression and thus were not included in the parameter 

variance–covariance matrix of the calibrated model may still 

be important for model predictions (for this study the size, 

shape, and location of the area contributing recharge to the 

production wells). The remaining three glacial hydraulic 

parameters (KV_VF, KV_TILL, and GHB_MULT) and 

the parameter representing the maximum groundwater 

evapotranspiration in the valley-fill (ETM) were, therefore, 

incorporated into the parameter variance–covariance matrix 

using the specified parameter values. Parameters representing 

bedrock and surface water were not changed in this analysis. 

Incorporating these four parameters into the parameter 

variance–covariance matrix, however, caused large unrealistic 

uncertainties around the specified parameter values because 

the information that observations provided was insufficient. 

Prior information on these parameters from the literature was 

used to constrain this uncertainty.

A prior weight was used for KV_VF so that the lower 

limit of the 95-percent confidence was not less than 1. 

However, because the specified value for KV_VF was 5, 

which is relatively close to 1, the upper limit of the 95-percent 

confidence interval was constrained to 25, which may not 

include most plausible values based on study results in 

Rhode Island and Cape Cod, Massachusetts, mentioned in the 

section, “Estimation of Model Parameters.” This 95-percent 

confidence interval, however, does include most of the values 

for coarse-grained stratified deposits listed in these studies. A 

prior weight was applied to KV_TILL in the same manner as 

that for KV_VF, and it too resulted in a 95-percent confidence 

interval of 1 to 25. Prior weight for ETM was ±35 percent 

around the specified value of 21 in/yr, which is the growingseason rate of free-water-surface evaporation. The 95-percent 

confidence interval (13.8 to 28.3 in/yr) includes the annual 

rate of free-water-surface evaporation (28 in/yr; Farnsworth 

and others, 1982), a possible maximum value of groundwater 

evapotranspiration. A prior weight was used for parameter 

GHB_MULT so that the 95-percent confidence interval 

ranged from 0.25 to 4, which is equivalent to a hydraulic 

conductivity of 12 ft/d (close to a value representative of very 

fine sand) and 200 ft/d (sand and gravel) around the specified 

value of 50 ft/d (fine sand). Conceptually, a higher hydraulic 

conductivity value for the general-head boundary would allow 

pumping wells to induce more groundwater from adjacent 

basins than would a lower hydraulic conductivity value. The 

more groundwater induced from an adjacent basin, the smaller 

the contributing area to a well in the active model area when 

compared to a model simulation in which no groundwater or 

less water was induced across the groundwater divide.

The addition of these four inestimable parameters 

(KV_VF, KV_TILL, ETM and GHB_MULT) into the parameter variance–covariance matrix incorporated into the Monte 

Carlo analysis all parameter uncertainty potentially important 

for model predictions. Parameter uncertainties are from the 

observation dataset, but also from prior information on parameters that the modeler provided.

Parameter sets created by equation 1 are shown in 

figure 13A. The hydraulic parameters were log-transformed 

in the model. The parameter sets have a lognormal or normal 

distribution around the optimal or specified parameter value; 

the spread of these data indicates the certainty with which 

each parameter was estimated, considering the available 

observations and prior information. Parameters incorporated 

into the variance–covariance matrix using prior information 

(KV_VF, KV_TILL, GHB_MULT, and ETM) have the 

least certainty, which would be expected because a prior 

weight was used so as to incorporate most plausible values. 

Parameters sets, however, with glacial anisotropies less than 1 

(horizontal hydraulic conductivity less than vertical hydraulic 

conductivity) or greater than 50 were removed (conditioned) 

from the parameters sets because they were considered 

unrealistic for the aquifer as a whole.

For the Monte Carlo analysis, the model was first 

run with 500 parameter sets, and with the 1996 average 

pumping rates used in calibrating the model. (A Monte Carlo 

analysis with 600 parameter sets showed similar results for 

the probabilistic contributing area, indicating that 500 sets 
36 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

-3 -2 -1 210 3 4

0

1

2

3

4

0

1

2

3

4

0.00

0.04

0.08

0.12

0.16

0.00

0.04

0.08

0.12

0.16

GHB_MULT

K_MULT

K_TILL

KV_VF

KV_TILL

SB_MULT

ETM

R_VF

R_TILL

Probability density

Logarithm of parameter value, dimensionless, 

except for K_TILL, which is feet per day

Inches per year

A

B

-3 -2 -1 210 3 4

0 5 10 15 20 25 30 35 40

0 5 10 15 20 25 30 35 40

EXPLANATION EXPLANATION

Figure 13. Model-parameter distributions A, before acceptance criteria and B, after acceptance criteria were applied to determine 

the probabilistic areas contributing recharge to the well centers, (parameter information is provided on table 2), Hunt–Annaquatucket–

Pettaquamscutt River Basins, Rhode Island.
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 37

were sufficient to use in the analysis). The nine hydraulic, 

recharge, and groundwater evapotranspiration parameter 

values in each dataset replaced the corresponding parameter 

values in the calibrated model. The following three criteria 

for accepting a given parameter set were used: (1) the model 

converged, (2) the model mass balance was 2 percent or less, 

and (3) a model-fit statistic (calculated error variance) was 

less than a specified value. Of the 500 parameter sets run with 

MODFLOW, 425 sets (85 percent) converged, and of these, 

420 sets (84 percent) had a water budget of 2 percent or less, 

thus meeting two of the three criteria.

The third acceptance criterion, the model-fit statistic, 

was used so that model-prediction uncertainty would not be 

overestimated by using a parameter set that produced unrealistic groundwater levels or streamflows compared to that for the 

calibrated model. The value used for this criterion, however, 

can be model dependent and subjective. For example, a calculated error variance of 20 or less was used for a Monte Carlo 

analysis of areas contributing recharge to wells in a wetland 

setting, for which the calibrated model had a calculated error 

variance of 5.3 (Friesz, 2010); this criterion did not remove 

any datasets. In contrast, Starn and others (2010) used a 

standard error of regression of 20 (calculated error variance of 

400) for a comparable Monte Carlo analysis (prior information 

on selected parameters and the same types of observations) 

in a valley-fill setting, for which the calibrated model had a 

standard error of regression of 6.2 (calculated error variance 

of 38.4); this value for the criterion removed 26 percent of the 

datasets after the dataset met a water-budget criterion.

For this model application of the Monte Carlo analysis, 

the third criterion was a calculated error variance of 20. 

Most parameter sets had a calculated error variance of 15 

or less, and 4 sets had greater than 20. Common to these 

four parameter sets was the fact that the K_MULT and 

K_TILL both had values greater or less than their 95-percent 

confidence intervals. All four K_TILL values were less than 

its confidence interval, and three K_MULT values were 

greater, and one less, than its confidence interval. Although the 

individual parameter values were considered reasonable based 

on model calibration, the combination of these parameter 

values produced a poor model fit. These four parameter sets 

were removed from the analysis, and therefore, 416 sets 

(83 percent) of the 500 parameter sets run with MODFLOW 

fit the three acceptance criteria. The distribution of parameters 

after the unrealistic glacial anisotropies were removed, and 

after the acceptance criteria were applied, are shown in 

figure 13B. Except for the glacial anisotropies, the distribution 

of parameters, although slightly altered from the original 

parameter sets, indicated a generally lognormal or normal 

distribution. Using a higher, yet plausible, specified value of 

10 instead of 5 for the glacial anisotropies in the calibrated 

model would require less conditioning of these parameter sets 

and less altering of the distribution.

Monte Carlo analyses were then done by use of the 

parameter sets that fit the acceptance criteria for the 1996 

average withdrawal rates but using the 2004–08 average 

withdrawal rates and the maximum withdrawal rates. The criteria for the Monte Carlo analyses that used these 416 parameter sets and these pumping rates were that the water budget be 

2 percent or less for models that converged. For the 2004–08 

average withdrawal rate, 402 parameter sets (97 percent) fit 

these criteria and thus were run with the particle-tracking program. For the maximum withdrawal rates, 408 parameter sets 

(98 percent) fit these criteria. The probability that a recharge 

location would be in the area contributing recharge to the 

production wells was determined by dividing the number of 

times a particle at a given location was captured by a well by 

the total number of accepted particle-tracking simulations; this 

probability was expressed as a percentage.

The probabilistic areas contributing recharge to the 

well centers at each center’s 2004–08 average pumping 

rate are shown in figure 14, and those at the maximum 

pumping rate are shown in figure 15. For this analysis, the 

probabilistic distribution is not by well but by well center. 

Probabilistic contributing areas to individual wells may 

overlap, even though the deterministic contributing areas do 

not overlap under a steady-state simulation. The total size 

of the probabilistic contributing area for each well center 

for both average and maximum pumping rates was larger 

than the deterministic contributing areas for the well center 

because of the effects of parameter uncertainty. This indicated 

that some areas not in the deterministic contributing area, 

including additional areas of urban and agricultural land use, 

may actually be in the contributing area. Generally, areas 

closest to the well centers with short traveltimes are associated 

with higher probabilities, whereas areas distant from the 

well centers with long traveltimes are associated with lower 

probabilities. In most cases, areas with high probabilities 

(greater than 50 percent) generally coincided with the 

deterministic contributing areas for this model.

For well center 1, the deterministic contributing area 

for the average pumping rate (fig. 10) for the five closely 

spaced wells adjacent to the lower Hunt River extended 

southward from the wells in a relatively narrow width to the 

groundwater divide between the Hunt River and Sandhill 

Brook, and it extended northwestward from the wells to 

include most of a till hillslope between two tributary streams 

that drains toward the valley. The probabilistic contributing 

area (fig. 14) indicated that additional areas southward toward 

the groundwater divide and northwestward in the till uplands, 

mostly east of Fry Brook, may be in the contributing area. 

Additional recharge originating in upland areas, mostly along 

groundwater divides between Fry and Frenchtown Brooks 

and between Frenchtown Brook and small tributaries to 

Hunt River, including Scrabbletown Brook, may be in the 

contributing area to the well center. These are areas associated 

with low probabilities. Two separate Monte Carlo analyses 

were done for this well center to provide additional insights 

(results not shown): one to determine the probabilistic 

contributing area to only the five wells along the lowest part of 

the lower Hunt River, and one to well 3A near the confluence 

with Frenchtown Brook and the Hunt River. Areas of low 
38 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

71°33' 71°30' 71°27'

41°39'

41°36'

41°33'

Frenchtown

North Kingstown

West

Greenwich

East Greenwich

Warwick

West

Warwick

Exeter

Narragansett Bay

Brook

Coventry

NK4

NK5

NK7

NK8

3A

IW

SFH1

SFH2

SFH3

Well Center 1

Well Center 2

Well Center 3

Well Center 4

Probabilistic contributing areas to

well centers, in percent

Part of contributing area

to Well Center 4

Part of contributing area

to Well Center 1

Part of

contributing area

to Well Center 2

NK4

SFH1

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83

0 1 2 KILOMETERS

0 1 2 MILES

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Queens Fort Brook

Scrabbletown

Queen

Brook

River

Unnamed

Mawney Brook

Fry

Brook

Tributary

Hunt River

Pettaquamscutt

Sandhill Brook

River

Belleville 

Pond

Secret 

 Lake

Carr 

Pond Chipuxet River

Annaquatucket

Cocumcossuc Brook

Mill Creek

River

Mattatuxet River

NK6

NK9 NK10

KC1

9A

14A

Greater than 50 to 100

Greater than 10 to 25

Greater than 25 to 50

Greater than 0 to 10

EXPLANATION

NK1

NK2

NK11

NK3

Figure 14. Simulated probabilistic areas contributing recharge to the well centers at their 2004–08 average 

pumping rates, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 39

probabilities along both upland groundwater divides can be 

either to the well field along the lower Hunt River or to well 

3A. The probabilistic contributing area to well 3A did not 

include areas of high probabilities (greater than 50 percent), 

indicating that for most parameter sets used in the Monte 

Carlo analysis the well derived most of its water from stream 

infiltration. Areas of highest probabilities for this well (greater 

than 10 to 50 percent) are near the upland-valley contact and 

west of the well in the uplands; the deterministic contributing 

area for this well corresponded to these areas. In addition to 

areas west of the Hunt River, the probabilistic contributing 

area for well 3A also included an area of low probabilities 

south of the well near the groundwater divide between 

Hunt River and Sandhill Brook. As with the deterministic 

contributing area, the probabilistic contributing area did not 

overlie well 3A.

At the maximum pumping rate, the deterministic contributing area (fig. 11) to the five wells adjacent to the lower 

Hunt River expanded to include more valley-fill deposits 

south and east of the wells toward groundwater divides and to 

include more valley-fill deposits and till north and northwest 

of the wells toward topographical divides that serve as model 

boundaries. The deterministic contributing area to one of 

these five wells also included small areas in the till between 

Fry and Frenchtown Brooks. The deterministic contributing 

area for the maximum pumping rate for well 3A also extended 

upgradient on both sides of the Hunt River. The deterministic 

contributing areas generally corresponded to high probabilities 

(fig. 15A), including well 3A. Along margins of the probabilistic contributing area south and east of the well center and 

north and northwest of the five wells along the lower Hunt 

River, there is minimal spread in low probabilities because 

these areas are close to the wells and are constrained by 

hydrologic and model boundaries. Extensive areas associated 

with mostly low probabilities to both well 3A and to the wells 

along the lower Hunt River are along the groundwater divides 

in upland areas. The deterministic contributing areas in the 

uplands between Fry and Frenchtown Brooks are associated 

with mostly high probabilities for well 3A, but with probabilities mostly between 10 to 50 percent for the wells adjacent to 

the lower Hunt River.

For well center 2, the deterministic contributing area to 

NK6 for both average and maximum pumping rates extended 

southwestward from the well to the upland-valley contact and 

then northwestward in till uplands (figs. 10 and 11). The deterministic contributing area for the maximum pumping rate also 

included small, isolated areas west of the well adjacent to and 

beneath Scrabbletown Brook. The probabilistic contributing 

area to the well for both pumping rates had a large spread in 

low probabilities in the uplands compared to those for the valley-fill deposits (figs. 14 and 15B). The probabilistic contributing area for the average pumping rate included small, isolated 

areas west of the well near Scrabbletown Brook and adjacent 

uplands, areas that are associated with low probabilities. The 

probabilistic contributing areas for the maximum pumping rate 

included more areas west of the well adjacent to and beneath 

Scrabbletown Brook and toward the uplands; these areas are 

associated with a large range in probabilities. The deterministic contributing area in this area for the maximum pumping 

rate corresponded to an area of high probabilities (greater than 

50 percent) and to probabilities greater than 25 to 50 percent. 

For both pumping rates, these are additional areas that may 

be in the contributing area to the well that also have relatively 

short groundwater-flow paths and traveltimes.

For well center 3, the deterministic contributing area 

for the average pumping rate extended westward and 

southwestward from the wells to the edge of the model in 

relatively narrow bands, and it also underlies part of Secret 

Lake (fig. 10). The probabilistic contributing area indicated 

that additional areas west and southwest of the wells along the 

margins of the deterministic contributing area may actually be 

in the contributing area to the well (fig. 15B). This included 

some recharge near NK4 and NK5 that discharged to small 

tributaries to Secret Lake in the deterministic model. In 

contrast to the deterministic model, the Monte Carlo analysis 

indicated that NK5 may be overlain by the contributing area to 

the well center (probabilities of 2 percent and less). Additional 

areas beneath Secret Lake may also be in the contributing area 

to the well center.

The deterministic contributing area for the maximum 

pumping rate included all of the area between the well center 

and the west edge of the model, beneath most of Secret Lake 

and part of Belleville Pond, and small areas on the opposite 

side of Belleville Pond and upland till in Queens Fort Brook 

Basin (fig. 11). Along margins of the probabilistic contributing 

area (fig. 15B), there is minimal spread in low probabilities 

because of hydrologic and model boundaries; Belleville Pond 

and induced groundwater from Chipuxet River Basin can be 

major sources of water to the well center, depending on the 

values of a parameter set used in the Monte Carlo simulation. 

Withdrawals by the fish hatchery wells may also constrain the 

extent of the contributing area northwest of the well center. On 

the side of Belleville Pond opposite the well center, probabilities decrease in a radially outward pattern; most of the deterministic contributing area on this side of the pond coincided 

with probabilities greater than 10 to 50 percent. Additional 

areas in Queens Fort Brook that are not in the deterministic 

contributing area are in the probabilistic contributing area. 

These are areas of mostly low probabilities; the deterministic 

contributing area in the till uplands corresponded with probabilities of 25 percent and less.

For well center 4, the deterministic contributing area for 

both pumping rates extended upgradient of the well center to 

upland till on both sides of a small valley (figs. 10 and 11). 

The deterministic contributing area for the maximum pumping rate also included a small area downvalley beneath Carr 

Pond and on the opposite side of the pond in the uplands. The 

probabilistic contributing area for both pumping rates (figs. 14 

and 15B) indicated that additional areas upgradient of the well 

center in the valley and uplands may be in the contributing 

area, but the deterministic and probabilistic contributing areas 

also indicated that recharge discharges to headwater streams of 
40 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

71°33' 71°30' 71°27'

41°39'

41°36'

Frenchtown

North Kingstown

West Greenwich

East Greenwich

Warwick

West Warwick

Exeter

Coventry

3A

IW

Well Center 1

Probabilistic contributing areas to

well centers, in percent

NK10

IW

Hunt River

Mill Creek

Scrabbletown

Unnamed

Mawney Brook

Fry

Sandhill Brook

Brook

Brook

Tributary

Brook

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83

0 1 2 KILOMETERS

0 1 2 MILES

A

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Greater than 50 to 100

Greater than 10 to 25

Greater than 25 to 50

Greater than 0 to 10

EXPLANATION

NK9 NK10

KC1

9A

14A

NK6

Figure 15. Simulated probabilistic areas contributing recharge to A, well center 1 and B, well centers 2, 3, and 4, 

at their maximum pumping rates, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.
Simulation of Areas Contributing Recharge and Prediction Uncertainty Analysis 41

71°33' 71°30' 71°27'

41°36'

41°33'

North Kingstown

West Greenwich

East Greenwich

Exeter

Narragansett Bay

IW

SFH1

SFH2

SFH3

Well Center 2

Well Center 3

Well Center 4

Part of contributing area

to Well Center 4

Part of contributing area

to Well Center 3

Probabilistic contributing areas to

well centers, in percent

NK11

IW

Annaquatucket

Cocumcossuc Brook

Queens Fort Brook

Scrabbletown

Queen

Belleville 

Pond

Brook

River

Mattatuxet

Mill Creek

Chipuxet River

River

River

Secret 

 Lake

Carr 

Pond

Pettaquamscutt

River

Hunt River

Unnamed Frenchtown Sandhill Brook

Brook

Tributary

Base from U.S. Geological Survey, Rhode Island

state plane projection, 1:24,000, NAD 83

0 1 2 KILOMETERS

0 1 2 MILES

B

Valley-fill deposits

Till and bedrock outcrops

Water

Model boundary

Production well—drinking water

Production well—other

Geology

Greater than 50 to 100

Greater than 10 to 25

Greater than 25 to 50

Greater than 0 to 10

EXPLANATION

NK1

NK2

NK3

NK4

NK5

NK6

NK7

NK8

NK11

Figure 15. Simulated probabilistic areas contributing recharge to A, well center 1 and B, well centers 2, 3, and 4, 

at their maximum pumping rates, Hunt–Annaquatucket–Pettaquamscutt River Basins, Rhode Island.—Continued
42 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

this valley and uplands before it is available for infiltration as 

the stream flows near the well center. In addition to additional 

areas beneath Carr Pond for the maximum pumping rate, the 

probabilistic contributing area for both pumping rates included 

areas on the opposite side of the pond from the wells; these are 

areas with low probabilities. In the deterministic model for the 

average pumping rate, NK7 and NK8 were not covered by the 

contributing area for the well center but the probabilistic contributing area indicated, with probabilities greater than 25 to 

50 percent, that this area may be in the contributing area.

Limitations of Model

The finite-difference numerical model of the HAP study 

area is a regional-scale simulation of groundwater flow, water 

levels, and the interaction between groundwater and surface 

water. Simplification included consolidating parameters that 

represented hydraulic properties and boundary conditions into 

homogenous units and assigning these parameters to groups 

of model grid cells sized 200 by 200 ft. Further, recharge rates 

and their distribution do not take into account impermeable 

surfaces in urban areas, which are likely to reduce or redistribute recharge in these areas. For these reasons, the model may 

not be appropriate for simulating local-scale results, but, for 

evaluating aquifer vulnerability to contamination at a regional 

scale, the model is useful.

Groundwater flow in upland till and bedrock was greatly 

simplified in the model. Although hydraulic conductivity 

values and recharge rates can be highly variable in till deposits, only three parameters represented hydraulic properties 

and recharge rates in these unsorted materials. For example, 

recharge rates may range from near zero in low-permeability 

tills in areas of steep topography to values approaching those 

for water available for recharge in sandy tills on moderate 

slopes. Thus the parameter values represented an average 

value for these deposits. Groundwater flow in bedrock represented the bulk flow in the regional system instead of the flow 

through bedrock fractures. Thus the model accounts for the 

overall movement of groundwater through the upland tills and 

bedrock before it discharges to the valley-fill deposits where 

the production wells are located.

The areas contributing recharge to the well center in 

the Annaquatucket River Basin extended to the edge of the 

model in valley-fill deposits that a groundwater divide defined 

and that a general-head boundary represented. Although the 

quantity of groundwater induced across this boundary from the 

adjacent basin was simulated for the average and maximum 

pumping rates, it would be necessary to extend the model into 

the adjacent basin in order to delineate the actual area that 

contributes water to the well center. Additional geologic and 

hydrologic data would be needed in order to expand the model 

and to improve the understanding of the geohydrology in 

this area.

Uncertainty in the simulated areas contributing recharge 

to the wells was based on the observation dataset and not 

from model design. Additional groundwater-level and baseflow observations and other types of field observations may 

help reduce the uncertainty about the extent of the simulated 

contributing area by increasing the precision of the parameter value estimates. Additional observations may also help 

to increase the number of parameters that could be estimated 

using nonlinear regression, thereby decreasing the need for 

using prior information from the literature. The resulting 

uncertainty analysis would then be based solely on objective 

model-calibration data. Also, the seven base-flow observations 

used in the calibration represented the net gain in base flow 

over relatively large segments of stream reaches upstream of 

the observation. To gain a better understanding of groundwater 

and surface-water interactions and of streambed conductance 

near the pumping wells requires observations that represent a 

net gain or loss in base flow over shorter stream reaches near 

the wells.

The simulated groundwater traveltimes were based on the 

calibrated model and uniform porosity for each of the lithologic units. An uncertainty analysis of groundwater traveltimes based on the spatial variability and the plausible range 

in porosity was beyond the scope of this study. However, in 

a steady-state model, porosity does not affect the location or 

size of the simulated area that contributes recharge to a well. 

Finally, traveltimes do not take into account traveltime in the 

unsaturated zone between the land surface and the water table.

Summary and Conclusions

The U.S. Geological Survey (USGS), in cooperation with 

the Rhode Island Department of Health, Office of Drinking 

Water Quality, began a 2-year investigation in 2009 to 

increase understanding of groundwater flow and of the areas 

contributing recharge to 15 production wells in the Hunt, 

Annaquatucket, and Pettaquamscutt (HAP) Rivers Basins 

in central Rhode Island. These large capacity production 

wells are operated by the town of North Kingstown, the 

Rhode Island Economic Development Corporation, and 

the Kent County Water Authority. A total average daily 

rate of 5.0 million gallons per day is withdrawn from these 

production wells, which are screened in coarse-grained valleyfill deposits of glacial origin.

The area contributing recharge to a well is defined as the 

surface area where water recharges the groundwater and then 

flows toward and discharges to the well. Areas contributing 

recharge to the production wells were determined on the basis 

of a numerical steady-state groundwater-flow model representing long-term average hydrologic conditions. The study 

modified an existing model that simulated groundwater flow 

in the valley-fill deposits in order to (1) simulate flow in the 

adjoining uplands and (2) represent boundary conditions and 

hydraulic properties as parameters for calibration by inverse 
Summary and Conclusions 43

modeling using nonlinear regression and for evaluating modelprediction uncertainty. These changes increased the extent of 

the active modeled area from 18 to 47 square miles (mi2). The 

optimal parameter set from the calibrated model provided a 

single, best representation of the areas contributing recharge 

(deterministic contributing areas). The uncertainty analysis led 

to contributing areas expressed as a probability distribution 

(probabilistic contributing areas).

Groundwater flow in the HAP study area was simulated 

by a four-layer model representing surficial deposits and 

the underlying bedrock. The model was calibrated to 165 

groundwater-level observations and 7 streamflow (base flow) 

observations. Ten parameters representing boundary conditions and hydraulic properties were evaluated for nonlinear 

regression: two recharge parameters, three horizontal hydraulic-conductivity parameters, three vertical anisotropy parameters, a streambed hydraulic conductance parameter, and a 

maximum groundwater evapotranspiration parameter. Five of 

these parameters were estimated by nonlinear regression. The 

remaining five parameters, which represented vertical anisotropy of the surficial deposits, hydraulic properties of bedrock, 

and the maximum groundwater evapotranspiration from the 

valley-fill deposits, were specified on the basis of values previously reported in the literature because observations alone 

did not provide sufficient information on them. A hydraulic 

conductance parameter for head-dependent cells representing 

a groundwater divide was not considered for estimation, but 

it was included in model-prediction uncertainty. A model-fit 

statistic, the calculated error variance, which is a measure of 

the overall magnitude of the weighted residuals, was 8.91. 

Influence statistics indicated that 13 groundwater-level and 6 

base-flow observations had the most overall influence in the 

regression and in the optimal set of estimated parameters. Six 

of seven base-flow observations were influential to the set of 

estimated parameters, and in addition, two of the base-flow 

observations were the only observations important to all of the 

individually estimated parameters, thus indicating the importance of this type of observation in model calibration.

The estimated optimal parameter value for recharge to 

the valley-fill deposits was 27.3 inches per year (in/yr), and 

the estimated optimal parameter value for effective recharge 

to upland till deposits was 18.7 in/yr. The estimated optimal 

parameter value for horizontal hydraulic conductivity of 

upland till deposits was 16.2 feet per day (ft/d). Two estimated parameters defined multipliers, one of which applied 

to spatially varying horizontal hydraulic conductivity of the 

valley-fill deposits and the other to spatially varying streambed 

conductance. Calibrated horizontal hydraulic conductivity values of the valley-fill deposits ranged from 20 to 480 ft/d and 

calibrated streambed conductance values ranged from 10,000 

to 52,000 feet squared per day (ft2/d). Optimal parameter values representing processes in the valley-fill deposits (recharge 

and horizontal hydraulic conductivity) were the most precisely 

estimated, whereas the horizontal hydraulic conductivity of till 

was the least precisely estimated parameter value.

The simulated hydrologic budget indicated that recharge 

is the dominant inflow to the upland deposits. Groundwater 

in the uplands discharged mostly to streams (52 percent) and 

to the valley-fill deposits (38 percent). Most inflow to the 

valley-fill deposits was from recharge (60 percent), but upland 

sources contributed a significant amount (about one-third), 

either directly or indirectly. These upland sources of inflow to 

the valley-fill deposits include lateral flow from till (15 percent), upward flow from bedrock (2 percent), and a large 

percentage of the streamflow loss (17 percent).

The 15 production wells are clustered in four pumping 

centers: two centers in the Hunt River Basin contain seven 

of the wells, one center in the Annaquatucket River Basin 

contains five wells, and the fourth center in the Pettaquamscutt 

River Basin contains three wells. Areas contributing recharge 

to the well centers were simulated for each center’s 2004–08 

average withdrawal rates (ranging from 235 to 1,858 gallons per minute (gal/min)) and for each center’s maximum 

pumping capacities (ranging from 800 to 8,500 gal/min). 

Simulated areas contributing recharge extend upgradient of 

the well centers to upland till and to groundwater and topographical divides. Some areas contributing recharge include 

small, isolated areas remote from the well center including, 

for the maximum pumping rate scenario, on the opposite 

side of surface-water bodies from the well center. Particle 

tracks indicated that recharge originating in the upland till and 

valley-fill deposits travels along deep groundwater-flow paths 

in the valley-fill deposits and, under pumping conditions, 

passes beneath the surface-water body to the well center. At 

the average pumping rate, four of the wells are not overlain by 

the simulated contributing area because recharge travels along 

paths above and around the screen interval to the well. Simulated surface-water loss within some of the contributing areas 

affected their size. For the average pumping rates, the size of 

the areas contributing recharge to the four well centers ranged 

from 0.19 to 1.12 mi2 and covered a total area of 2.79 mi2, and 

for the maximum pumping rate, the size ranged from 0.37 to 

3.53 mi2 and covered a total area of 7.99 mi2. The well center 

in the Annaquatucket River Basin also derived part of its water 

from induced groundwater across a simulated groundwater 

divide with an adjoining basin; this induced groundwater was 

equivalent to recharge over an area of 0.21 mi2. Land cover in 

the areas contributing recharge includes a substantial amount 

of urban land use for the well centers in the Hunt River Basin, 

agriculture and sand and gravel mining for the well center in 

the Annaquatucket River Basin, and primarily undeveloped 

land for the well center in the Pettaquamscutt River Basin. 

Simulated groundwater traveltimes from recharge locations to 

production wells for the maximum pumping rate ranged from 

less than 1 year to greater than 350 years for each well center. 

Median traveltimes ranged from 2.9 to 5.0 years, and traveltimes of 10 years or less ranged from 78 to 93 percent for the 

well centers. These relatively short traveltimes indicated that 

the wells are vulnerable to contamination from activities on 

the land surface.
44 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

Parameter uncertainty and its associated effects on the 

simulated areas contributing recharge to the well centers were 

evaluated using a stochastic Monte Carlo analysis. Optimal 

parameter values and the parameter variance–covariance 

matrix from nonlinear regression were used to create parameter sets for the analysis. The parameter variance–covariance 

matrix preserves the uncertainty of the parameter estimates 

and the correlation among parameters from the calibrated 

model. The four inestimable parameters representing glacial 

hydraulic properties and maximum groundwater evapotranspiration were also incorporated into the parameter variance–

covariance matrix. Because observations did not provide 

enough information to constrain the uncertainty of these four 

parameters within realistic ranges around the specified values, 

prior information was required. Thus all parameters that may 

be important for model predictions were incorporated into the 

analysis. The uncertainty analysis was an outcome of calibrating the model to available observations, but it also depended 

on information provided by the modeler. Three acceptance 

criteria were used to assess parameter sets so that prediction 

uncertainty was not overestimated: the model converged, 

model mass balance was 2 percent or less, and the calculated 

error variance was 20 or less. Of 500 parameter sets using 

pumping values from the calibrated model, 416 fit the acceptance criteria. Four parameter sets that fit the first two criteria 

had a calculated error variance greater than 20. Common to 

these four sets was the fact that the two parameters representing horizontal hydraulic conductivity of the glacial deposits 

had values outside their 95-percent confidence interval. The 

2004–08 average pumping rate and the maximum pumping 

rate scenarios used these 416 parameter sets with the first two 

acceptance criteria.

The size of the probabilistic contributing areas for each 

well center for both average and maximum pumping rates 

was larger than the size of the deterministic contributing areas 

for the well center because of the effects of parameter uncertainty. Thus, some areas not in the deterministic contributing 

area, including additional areas with urban and agricultural 

land cover that has the potential to contaminate groundwater 

resources, may actually be in the contributing area. Generally, areas closest to the well centers with short groundwater 

traveltimes are associated with higher probabilities, whereas 

areas distant from the well centers with long groundwater 

traveltimes are associated with lower probabilities. In most 

cases, areas with high probabilities (greater than 50 percent) 

generally coincided with the deterministic contributing areas. 

In some cases, areas near the well centers where simulated 

streams and ponds intercepted recharge in the calibrated model 

were in the probabilistic contributing area, indicating that this 

recharge may instead go directly to a well. Three of the four 

wells not overlain by the deterministic contributing area for 

the average pumping rate were in the probabilistic contributing area that was associated with probabilities of 50 percent 

or less. Areas associated with low probabilities extended long 

distances along groundwater divides in the uplands remote 

from the well centers.

References Cited

Allen, W.B., 1956, Ground-water resources of the East Greenwich quadrangle, Rhode Island: Rhode Island Development 

Council Geological Bulletin 8, 56 p.

Allen, W.B., Hahn, G.W., and Tuttle, C.R., 1963, Geohydrological data for the upper Pawcatuck River Basin, Rhode 

Island: Rhode Island Water Resources Coordinating Board 

Geological Bulletin 13, 68 p., 1 pl., scale 1:24,000.

Allen, W.B., Johnson, K.E., and Mason, R.A., 1959, Groundwater map of the Crompton quadrangle, Rhode Island: 

Rhode Island Water Resources Coordinating Board 

GWM-3, 1 sheet, scale 1:24,000.

Barlow, P.M., 1997, Particle-tracking analysis of contributing 

areas of public-supply wells in simple and complex flow 

systems, Cape Cod, Massachusetts: U.S. Geological Survey 

Water-Supply Paper 2434, 66 p.

Barlow, P.M., and Dickerman, D.C., 2001, Numerical-simulation and conjunctive-management models of the Hunt–

Annaquatucket–Pettaquamscutt stream-aquifer system, 

Rhode Island: U.S. Geological Survey Professional Paper 

1636, 88 p.

Barlow, P.M., and Ostiguy, L.J., 2007, Simulation of hydrologic-system responses to ground-water withdrawals in 

the Hunt–Annaquatucket–Pettaquamscutt stream-aquifer 

system, Rhode Island: U.S. Geological Survey Open-File 

Report 2006–1226, 51 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 Water-Resources Investigations Report 95–4234, 

56 p.

Bent, G.C., 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, 68 p.

DeSimone, L.A., 1999, A vulnerability assessment of publicsupply wells in Rhode Island: U.S. Geological Survey 

Water-Resources Investigations Report 99–4160, 153 p.

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

evaluation of water-management alternatives in the Assabet 

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

Scientific Investigations Report 2004–5119, 133 p.

Dickerman, D.C., 1984, Aquifer tests in the stratified drift, 

Chipuxet River Basin, Rhode Island: U.S. Geological Survey Water-Resources Investigations Report 83–4231, 

39 p.
References Cited 45

Dickerman, D.C., and Barlow, P.M., 1997, Water-table 

conditions and stream-aquifer interaction in the Hunt–

Annaquatucket–Pettaquamscutt aquifer, central Rhode 

Island, October 7–9, 1996: U.S. Geological Survey WaterResources Investigations Report 97–4167, 1 sheet.

Dickerman, D.C., Kliever, J.D., and Stone, J.R., 1997, Hydrogeology, water quality, and simulation of ground-waterdevelopment alternatives in the Usquepaug–Queen groundwater reservoir, southern Rhode Island: U.S. Geological 

Survey Water-Resources Investigations Report 97–4126, 

48 p., 1 sheet.

Dickerman, D.C., Trench, E.C.T., and Russell, J.P., 1990, 

Hydrogeology, water quality, and ground-water development alternatives in the lower Wood River ground-water 

reservoir, Rhode Island: U.S. Geological Survey WaterResources Investigations Report 89–4031, 109 p.

Farnsworth, R.K., Thompson, E.S., and Peck, E.L., 1982, 

Evaporation atlas for the contiguous 48 United States: U.S. 

Department of Commerce, National Oceanic and Atmospheric Administration, Technical Report NWS 33, 

26 p., 4 pls.

Friesz, P.J., 2004, Delineation of areas contributing recharge 

to selected public-supply wells in glacial valley-fill and 

wetland settings, Rhode Island: U.S. Geological Survey 

Scientific Investigations Report 2004–5070, 57 p.

Friesz, P.J., 2010, Delineation and prediction uncertainty 

of areas contributing recharge to selected wells fields in 

wetland and coastal settings, southern Rhode Island: U.S. 

Geological Survey Scientific Investigations Report 2010–

5060, 69 p.

Friesz, P.J., and Stone, J.R., 2007, Simulation of ground-water 

flow and areas contributing recharge to production wells in 

contrasting glacial valley-fill settings, Rhode Island: U.S. 

Geological Survey Scientific Investigations Report 2007–

5133, 50 p.

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

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

sand and gravel, Cape Cod, Massachusetts: 2. Analysis of 

spatial moments for a nonreactive tracer: Water Resources 

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

Granato, G.E., 2009, Computer programs for obtaining and 

analyzing daily mean streamflow from the U.S. Geological 

Survey National Water Information System Web Site: 

U.S. Geological Survey Open-File Report 2008–1362, 

123 p. on CD–ROM, 5 appendixes.

Hahn, G.W., 1959, Ground-water map of the Slocum quadrangle, Rhode Island: Rhode Island Water Resources Coordinating Board GWM-2, 1 sheet, scale 1:24,000.

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.

Hickey, P., and Joubert, L., 2003, Kent County and North 

Kingstown Source Water Assessment: University of Rhode 

Island, Technical Report 03–06: Kingston, R.I., 86 p.

Hill, M.C., 1998, Methods and guidelines for effective model 

calibration: U.S. Geological Survey Water-Resources 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.

Hill, M.C., and Tiedeman, C.R., 2007, Effective groundwater model calibration—with analysis of data, sensitivities, 

predictions, and uncertainty: New Jersey, John Wiley and 

Sons, Inc., 455 p.

Hirsch, R.M., 1982, A comparison of four streamflow record 

extension techniques: Water Resources Research, v. 18, 

no. 4, p. 1081–1088.

Johnson, K.E., and Marks, L.Y., 1959, Ground-water map of 

the Wickford quadrangle, Rhode Island: Rhode Island Water 

Resources Coordinating Board GWM-1, 1 sheet, scale 

1:24,000.

Lindsey, B.D., 2005, Hydrogeology and simulation of source 

areas of water to production wells in a colluvium-mantled 

carbonate-bedrock aquifer near Shippensburg, Cumberland 

and Franklin Counties, Pennsylvania: U.S. Geological Survey Scientific Investigations Report 2005–5195, 49 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., Sorenson, J.R., Stone, J.R., Moran, S.B., and 

Hougham, A., 2007, Hydrogeology and simulated groundwater flow in the Salt Pond region of southern Rhode 

Island: U.S. Geological Survey Scientific Investigations 

Report 2006–5271, 56 p.

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

three-dimensional finite-difference ground-water flow 

model: U.S. Geological Survey Techniques of Water 

Resources Investigations, book 6, chap. A–1, 586 p.

Meinzer, O.E., 1923, The occurrence of ground water in the 

United States, with a discussion of principles: U.S. Geological Survey Water-Supply Paper 489, 321 p.
46 Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, R.I.

Melvin, R.L., de Lima, V., and Stone, B.D., 1992, The stratigraphy and hydraulic properties of tills in southern New 

England: U.S. Geological Survey Open-File Report 91–481, 

53 p.

National Oceanic and Atmospheric Administration, 2009, 

Annual climatological summary, accessed November 2009, 

at http://cdo.ncdc.noaa.gov/ancsum/ACS.

National Oceanic and Atmospheric Administration, 2010, 

Tides and currents–sea level trends, accessed January 2010, 

at http://tidesandcurrents.noaa.gov/sltrends/sltrends_station.

shtml?stnid=8452660.

Pollock, D.W., 1994, User’s guide for MODPATH/MODPATH_PLOT, version 3: A particle tracking post-processing 

package for MODFLOW, the U.S. Geological Survey 

finite-difference ground-water flow model: U.S. Geological 

Survey Open-File Report 94–464 [variously paged].

Power, W.R., Jr., 1957, Surficial geology of the Slocum 

quadrangle, Rhode Island: U.S. Geological Survey Geologic 

Quandrangle GQ–106, 1 sheet, scale 1:31,680.

Power, W.R., Jr., 1959, Bedrock geology of the Slocum 

quadrangle, Rhode Island: U.S. Geological Survey Geologic 

Quandrangle GQ–114, 1 sheet, scale 1:31,680.

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.

Quinn, A.W., 1952, Bedrock geology of the East Greenwich 

quadrangle, Rhode Island: U.S. Geological Survey Geologic 

Quandrangle GQ–17, 1 sheet, scale 1:31,680.

Quinn, A.W., 1963, Bedrock geology of the Crompton quadrangle, Rhode Island: U.S. Geological Survey Bulletin 

1158–B, 17 p.

Randall, A.D., Thomas, M.P., Thomas, C.E., Jr., and Baker, 

J.A., 1966, Water resources inventory of Connecticut; 

Part 1, Quinebaug River Basin: Connecticut Water 

Resources Bulletin, no. 8, 102 p.

Reilly, T.E., and Pollock, D.W., 1993, Factors affecting areas 

contributing recharge to wells in shallow aquifers: U.S. 

Geological Survey Water-Supply Paper 2412, 21 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–4135, 

81 p.

Rosenshein, J.S., Gonthier, J.B., and Allen, W.B., 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., 5 pls.

Rutledge, A.T., 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.

Schafer, J.P., 1961, Surficial geology of the Wickford quadrangle, Rhode Island: U.S. Geological Survey Geologic 

Quandrangle GQ–136, 1 sheet, scale 1:24,000.

Searcy, J.K., 1959, Flow-duration curves, manual of hydrology—part 2. low-flow techniques: U.S. Geological Survey 

Water-Supply Paper 1542–A, p. 1–33.

Smith, J.H., 1955, Surficial geology of the East Greenwich 

quadrangle, Rhode Island: U.S. Geological Survey Geologic 

Quandrangle GQ–62, 1 sheet, scale 1:31,680.

Smith, J.H., 1956, Surficial geology of the Crompton quadrangle, Rhode Island: U.S. Geological Survey Geologic 

Quandrangle GQ–94, 1 sheet, scale 1:31,680.

Starn, J.J., Bagtzoglou, A.C., and Robbins, G.A., 2010, Using 

atmospheric tracers to reduce uncertainty in groundwater 

recharge areas: Ground Water, v. 48, no. 6, p. 858–868.

Starn, J.J., Stone, J.R., and Mullaney, J.R., 2000, Delineation 

and analysis of uncertainty of contributing areas to wells 

at the Southbury Training School, Southbury, Connecticut: 

U.S. Geological Survey Water-Resources Investigations 

Report 00–4158, 53 p.

Stone, J.R., Schafer, J.P., London, E.H., DiGiacomo-Cohen, 

M.L., Lewis, R.S., and Thompson, W.B., 2005, Quaternary 

geologic map of Connecticut and Long Island Sound Basin: 

U.S. Geological Survey Scientific Investigations Map 2784, 

1:125,000 scale, with accompanying text, 72 p.

Tiedeman, C.R., Goode, D.J., and Hsieh, P.A., 1997, Numerical simulation of ground-water flow through glacial deposits 

and crystalline bedrock in the Mirror Lake area, Grafton 

County, New Hampshire: U.S. Geological Survey Professional Paper 1572, 50 p.

U.S. Environmental Protection Agency, 1991, Protecting the 

nation’s ground water: EPA’s strategy for the 1990’s: U.S. 

Environmental Protection Agency, Office of the Administrator, Report 21Z–1020, 84 p.

U.S. Geological Survey, 1980, Accuracy specifications for topographic mapping, in Technical instructions of the National 

Mapping Division: Reston, Va., chap. 1B4, p. 1–13.

Williams, R.B., 1964, Bedrock geology of the Wickford 

quadrangle, Rhode Island: U.S. Geological Survey Bulletin 

1158–C, 15 p.
Prepared by the Pembroke and West Trenton 

Publishing Service Centers.

For more information concerning this report, contact:

Director

U.S. Geological Survey

Massachusetts-Rhode Island Water Science Center

10 Bearfoot Road

Northborough, MA 01532

dc_ma@usgs.gov

or visit our Web site at:

http://ma.water.usgs.gov
Friesz—Evaluating Prediction Uncertainty of Areas Contributing Recharge to Well Fields of Multiple Water Suppliers, Rhode Island—SIR 2012–5114