Decrease in stage

Increase in stage

Reservoir Surface

Reservoir Surface

Aquifer

Aquifer

Bedrock

Bedrock

Stage

Saturated

Unsaturated

Water table

Flow direction

Ground-Water Contributions to Reservoir Storage and

the Effect on Estimates of Firm Yield for Reservoirs

in Massachusetts 

Scientific Investigations Report 2006-5045

Prepared in cooperation with the

Massachusetts Department of Environmental Protection

U.S. Department of the Interior

U.S. Geological Survey
Cover: Diagram showing the response of the aquifer to a decrease and an increase in reservoir stage, respectively.
Ground-Water Contributions to Reservoir 

Storage and the Effect on Estimates of Firm 

Yield for Reservoirs in Massachusetts 

By Stacey A. Archfield and Carl S. Carlson

Prepared in cooperation with the Massachusetts Department 

of Environmental Protection

Scientific Investigations Report 2006–5045

U.S. Department of the Interior

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

Gale A. Norton, Secretary

U.S. Geological Survey

P. Patrick Leahy, Acting Director

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

For product and ordering information:

World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

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

natural hazards, and the environment:

World Wide Web: http://www.usgs.gov

Telephone: 1-888-ASK-USGS

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

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Archfield, S.A., and Carlson, C.S., 2006, Ground-water contributions to reservoir storage and the effect on estimates 

of firm yield for reservoirs in Massachusetts: U.S. Geological Survey Scientific Investigations Report 

2006–5045, 27 p.
v

Contents

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

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

Purpose and Scope ..............................................................................................................................4

The Existing Firm-Yield-Estimator Model for Reservoirs in Massachusetts...............................4

Determination of Ground-Water Contributions to Reservoir Storage...................................................5

Numerical-Model Simulations of Idealized Reservoir-Aquifer Geometries ...............................8

Analytical Simulations of Idealized Reservoir-Aquifer Geometries .............................................8

Effect of Reservoir Shape on the Flow Rate between the Reservoir and Aquifer.....................8

Application of the Ground-Water-Flow Equations to the Firm-Yield-Estimator Model ...................10

Modifications to the Ground-Water-Flow Equations....................................................................11

Numerical-Model Simulations of Ground-Water-Flow Rates and Reservoir 

Water Levels .................................................................................................................14

Modified Firm-Yield-Estimator Model Simulations of Ground-Water-Flow Rates and 

Reservoir Water Levels ...............................................................................................14

Comparison of Simulated Ground-Water-Flow Rates and Reservoir Water Levels.......14

Iterative Loop to Calculate the Change in Reservoir Stage.........................................................17

Effect of Uncertainty in the Input Variables to the Ground-Water-Flow Equations..........................20

Sensitivity of Firm Yield to the Reservoir Stage-Storage Relation..............................................20

Sensitivity of Firm Yield to Other Reservoir and Aquifer Input Variables ..................................21

Potential Effect of Ground-Water Contributions on Firm Yields for Reservoirs 

in Massachusetts...........................................................................................................................23

Summary........................................................................................................................................................26

Acknowledgments.......................................................................................................................................26

References Cited..........................................................................................................................................27
vi

Figures

1–3. Diagrams showing: 

1. Possible sources and losses of water for a hypothetical drinking-water 

 reservoir................................................................................................................................2

2. A, Reservoir water level at equilibrium with the water table in the surrounding 

aquifer; B, response of the water table to a decrease in the water level of the 

reservoir relative to the water table; or C, to an increase in the water level of the 

reservoir relative to the water table.................................................................................3

3. A, Surface-water-body shapes; B, with underlying model grid (shown with 

surface-water-body shape 10) for reservoir-aquifer systems; and C, reservoirstage fluctuations used to evaluate the effects of reservoir-aquifer geometry on 

the analytical solution to the ground-water flow equation derived by Rorabaugh 

(1964)......................................................................................................................................7

4. Graphs showing differences between the numerically and analytically simulated 

ground-water-flow rates, and timing of the flows for selected reservoir-aquifer 

geometries shown in figure 3 .....................................................................................................9

5. Flow chart of the modified firm-yield-estimator model........................................................12

6. Maps showing locations of four reservoir-aquifer systems in Massachusetts: A, 

Fitchburg Reservoir; B, Morse Reservoir; C, Bearhole Reservoir; and D, Millham 

Reservoir......................................................................................................................................13

7. Diagram of model grids used in the numerical simulations of four idealized reservoiraquifer systems patterned after the four reservoir-aquifer systems shown in figure 6. 

Reservoir A is patterned after Fitchburg Reservoir (fig. 6A); Reservoir B is patterned 

after Morse Reservoir (fig. 6B); Reservoir C is patterned after Bearhole Reservoir (fig. 

6C); Reservoir D is patterned after Millham Reservoir (fig. 6D)..........................................15

8–10. Graphs showing:

8. A, Comparison of simulated average ground-water flow rates; and B, reservoir 

water-level elevations for four idealized reservoir-aquifer systems during a 2-year 

period of average climate and streamflow conditions (April 1975 through March 

1977).....................................................................................................................................18

9. Sensitivity of firm yield to the closure criterion used in the iterative loop of the 

firm-yield simulations for four idealized reservoir-aquifer systems..........................20

10. Shape of the stage-storage relation for two drinking-water reservoirs: A,

Bearhole Reservoir; and B, Fitchburg Reservoir..........................................................21

11. Sensitivity of firm yield to reservoir and aquifer characteristics for four idealized 

reservoir-aquifer systems ................................................................................................22

12. Map showing locations of nine reservoirs that have shorelines in contact with sand 

and gravel aquifers.....................................................................................................................24

Tables

1. Reservoir and aquifer characteristics for four idealized reservoir-aquifer systems......16

2. Estimated reservoir and aquifer characteristics for nine reservoirs in Massachusetts 

whose reservoir shorelines are in contact with a sand and gravel aquifer. ...................25

3. Potential effect of ground-water contributions to reservoir storage on firm yield for 

nine reservoirs in Massachusetts. ..........................................................................................25
vii

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

square mile (mi2) 2.590 square kilometer (km2) 

Volume

gallon (gal) 3.785 liter (L) 

million gallons (Mgal) 3,785 cubic meter (m3)

cubic foot (ft3) 0.02832 cubic meter (m3

)

cubic mile (mi3) 4.168 cubic kilometer (km3) 

Flow rate

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

gallon per day (gal/d) 0.003785 cubic meter per day (m3/d)

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

Transmissivity*

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

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

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

ABBREVIATIONS

CHD Constant-Head Boundary

CN runoff-curve number

FYE Firm-Yield-Estimator

GIS Geographic Information System

MassDEP Massachusetts Department of Environmental Protection

QPPQ streamflow-probability-probability-streamflow

Conversion Factors, Horizontal Datum, and Abbreviations

Abstract

Potential ground-water contributions to reservoir storage 

were determined for nine reservoirs in Massachusetts that had 

shorelines in contact with sand and gravel aquifers. The effect 

of ground water on firm yield was not only substantial, but furthermore, the firm yield of a reservoir in contact with a sand 

and gravel aquifer was always greater when the ground-water 

contribution was included in the water balance. Increases 

in firm yield ranged from 2 to 113 percent, with a median 

increase in firm yield of 10 percent. Additionally, the increase 

in firm yield in two reservoirs was greater than 85 percent. 

This study identified a set of equations that are based 

on an analytical solution to the ground-water-flow equation 

for the case of one-dimensional flow in a finite-width aquifer 

bounded by a linear surface-water feature such as a stream. 

These equations, which require only five input variables, 

were incorporated into an existing firm-yield-estimator (FYE) 

model, and the potential effect of ground water on firm yield 

was evaluated. To apply the FYE model to a reservoir in 

Massachusetts, the model requires that the drainage area to 

the reservoir be clearly defined and that some surface water 

flows into the reservoir. For surface-water-body shapes having 

a more realistic representation of a reservoir shoreline than a 

stream, a comparison of ground-water-flow rates simulated 

by the ground-water equations with flow rates simulated by a 

two-dimensional, finite-difference ground-water-flow model 

indicate that the agreement between the simulated flow rates 

is within ±10 percent when the ratio of the distance from the 

reservoir shoreline to the aquifer boundary to the length of 

shoreline in contact with the aquifer is between values of 0.5 

and 3.5.

Idealized reservoir-aquifer systems were assumed to 

verify that the ground-water-flow equations were implemented 

correctly into the existing FYE model; however, the modified 

FYE model has not been validated through a comparison of 

simulated and observed data. A comparison of simulated and 

observed reservoir water levels would further define limitations to the applicability of the ground-water-flow equations 

to reservoirs in Massachusetts whose shorelines are in contact 

with a sand and gravel aquifer.

Introduction

The amount of water that is available for withdrawal 

from a drinking-water supply reservoir depends on the water 

balance of the reservoir. Reservoirs accumulate water from 

surface-water inflows and from precipitation on the reservoir 

surface; available storage is reduced by evaporation, releases 

to meet instream-flow requirements, and drinking-water 

withdrawals (fig. 1). Many of the terms in the water balance 

of a reservoir are derived from natural processes that cannot be altered. Consequently, to prevent reservoir failure, 

a withdrawal rate is chosen that maximizes the amount of 

water available for supply, but only allows for the complete 

depletion of available storage during one time interval over 

the period of simulation. The withdrawal rate that accomplishes this objective is termed the firm yield of the reservoir. 

Demands on a reservoir that are equal to or less than the 

firm yield ensure that the reservoir will always meet these 

demands, assuming the reservoir is affected by climatic and 

hydrologic conditions similar to those used in the evaluation. 

As demands for drinking-water supplies increase and available 

water is further reduced, calculation of the firm yield provides 

an upper boundary that can be used to evaluate the capacity 

of a reservoir to meet additional demands in excess of current 

withdrawals. Alternatively, the firm yield can be compared to 

current reservoir withdrawals to determine if these withdrawals could be sustained through a drought.

Water in available storage may be supplemented or 

reduced if a reservoir is in contact with an aquifer, due to the 

inflow of water from the aquifer to the reservoir or the outflow 

of water from the reservoir to the aquifer; however, because 

many reservoirs receive most of their water from a surfacewater source, the ground-water term in the water balance is 

assumed to be negligible and is ignored in the calculation of 

the firm yield. For a reservoir whose shoreline is in contact 

with an aquifer, as is the case for approximately one-quarter of 

the more than 70 drinking-water supply reservoirs in Massachusetts, the effect of ground-water contributions on firm yield 

is unknown because available models that are used to estimate 

firm yield were constructed without consideration for the 

ground-water contribution to reservoir storage. 

Ground-Water Contributions to Reservoir Storage and 

the Effect on Estimates of Firm Yield for Reservoirs in 

Massachusetts

By Stacey A. Archfield and Carl S. Carlson
Figure 1. Possible sources and losses of water for a hypothetical drinking-water reservoir.

The reservoir storage and corresponding stage will 

increase or decrease due to changes in the reservoir-water 

balance. Therefore, if a reservoir is in hydrologic connection with an aquifer, a change in reservoir stage will induce 

a change in the water level of the aquifer. For the simplest 

case, assume that the reservoir stage is in equilibrium with the 

water table in the adjacent aquifer (fig. 2A). If the reservoir 

stage decreases (fig. 2B), water will move from the aquifer to 

the reservoir, thereby decreasing the aquifer water table and 

increasing the reservoir stage until the reservoir stage and 

aquifer water table are again in equilibrium. If the reservoir 

stage increases (fig. 2C), water will move from the reservoir to 

the aquifer, thereby decreasing the reservoir stage and increasing the aquifer water table until the reservoir stage and aquifer 

water table are again in equilibrium. The aquifer response is 

further complicated when time-varying changes in reservoir 

stage are imposed on the system prior to achieving equilibrium. The resulting aquifer response will then be a function 

of both the new change in reservoir stage and the decaying 

effects of the previous changes in reservoir stage. 

When the sum of the inflows to the reservoir is less than 

the sum of the outflows, the reservoir stage will decrease 

because water in storage is being used to meet the withdrawals 

from the reservoir. A decrease in reservoir stage relative to the 

 Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
Equilibrium

Decrease in stage

Increase in stage

A.

B.

C.

Reservoir Surface

Reservoir Surface

Reservoir Surface

Aquifer

Aquifer

Aquifer

Bedrock

Bedrock

Bedrock

Bedrock

Bedrock

Bedrock

STAGE

SATURATED

UNSATURATED

WATER TABLE

FLOW DIRECTION

EXPLANATION

Figure 2. A, Reservoir water level at equilibrium with the water table in the 

surrounding aquifer; B, response of the water table to a decrease in the water level of 

the reservoir relative to the water table; or C, to an increase in the water level of the 

reservoir relative to the water table.

water table in the adjacent aquifer would supplement storage 

to the reservoir because water will be released from the aquifer 

to the reservoir (fig. 2B). As the reservoir is being drawn down 

and reservoir stage is decreasing, the resulting ground-water 

contribution from the aquifer to the reservoir could provide 

additional water that will prevent the reservoir from failing 

at a yield that would have otherwise caused the reservoir to 

fail, thereby increasing the firm yield. To evaluate the potential effect of ground water on firm yield, the U.S. Geological 

Survey (USGS), in cooperation with the Massachusetts 

Department of Environmental Protection (MassDEP), began 

a study to determine ground-water contributions to reservoir 

storage and to evaluate the effect of those contributions on 

estimates of firm yield.

Introduction 
Purpose and Scope

The purpose of this report is to present the results of 

the USGS/MassDEP study on the effect of ground-water 

contributions on estimates of firm yield. The report documents 

the methodology used to quantify the ground-water contributions to reservoir storage, explains how the methodology was 

applied to the existing firm-yield-estimator (FYE) model, and 

analyzes how the modified FYE model can be used to understand the effect of ground-water contributions on estimates of 

firm yield. 

This report identifies and evaluates the applicability of a 

set of equations that can be used to estimate the contribution of 

ground water to reservoir storage due to time-varying changes 

in reservoir stage for reservoir-aquifer systems for which a 

more detailed, two- or three-dimensional numerical groundwater-flow model is not available. The report discusses (1) the 

implementation of the set of equations into the existing version 

of the FYE model, (2) the sensitivity of an estimate of firm 

yield to changes in the input variables for four hypothetical 

reservoir-aquifer systems that represent the range of reservoirs 

and aquifers in Massachusetts, and (3) the potential differences 

in estimates of firm yield that result from the consideration 

of this ground-water contribution to reservoir storage for nine 

reservoir-aquifer systems in Massachusetts previously studied 

by Waldron and Archfield (2006) that have their shorelines in 

contact with sand and gravel aquifers. Although numericalmodel simulations were used to determine the applicability of 

the set of equations and the implementation of the equations, a 

comparison of the results from the modified FYE model with 

observed reservoir values is beyond the scope of this report. 

Furthermore, the firm-yield simulations in this report assume 

that all the water in usable reservoir storage is available for 

public supply, and the simulations do not include any consideration for the effect that this yield may have on streamflow 

downstream from the reservoir.

The Existing Firm-Yield-Estimator Model for 

Reservoirs in Massachusetts

To gain a better understanding of the margin between 

current demands to drinking-water reservoirs and their 

respective firm yields, MassDEP developed a guidance 

document that details a procedure to determine the firm 

yield of a reservoir in Massachusetts that receives water 

only from surface-water inputs (Massachusetts Department 

of Environmental Protection, 1996). The procedures 

in this document were used to develop the Firm-YieldEstimator (FYE) model (Massachusetts Department of 

Environmental Protection, 2000). Waldron and Archfield 

(2006) evaluated and modified the procedure used to estimate 

the reservoir water-balance terms and the firm yield.

Based on the MassDEP FYE guidance document (1996) 

and the modifications by Waldron and Archfield (2006), the 

firm yield is computed by solving the monthly water balance 

(eq. 1) of a reservoir during concurrent records of historical 

precipitation and streamflow values that include the drought of 

the 1960s:

Si= Si-1 + AwiQri + Ari(Pi – Ei) – αi ni Qy – Qri – Qsi – Qoi , (1)

where

i is the month;

Si is the water in usable storage at the end of the 

month, in million gallons;

Si-1 is the amount of water in usable storage at 

the end of the previous month, in million 

gallons;

Awi is the contributing drainage area, in miles

 squared;1

Qri is the streamflow per mile of reservoir 

drainage, in miles;1

Ari is the area of the reservoir surface, in miles 

squared;1

Pi is the precipitation, in miles;1

Ei is the evaporation from the reservoir, in 

miles;1

αi is the peak-use factor, dimensionless;

ni is the number of days in the month;

Qy is the yield, in million gallons per day;

Qri is the required release, in million gallons per 

month;

Qsi is the uncontrolled spill, in million gallons per

 month; and

Qoi is the withdrawal from the reservoir by other

 users, in million gallons per month.

1Precipitation, streamflow, and evaporation are in units of length and, when 

multiplied by Ari, become volumes (in cubic miles) that are converted to million gallons.

Precipitation data were obtained from the nearest 

climate station; evaporation was estimated using multivariateregression relations developed by Fennessey and Vogel (1996). 

Streamflow was determined by the QPPQ- (streamflowprobability-probability-streamflow) transform method, which 

is a method to determine streamflow at an ungaged site by 

relating the flow-duration curve at a gaged location with the 

flow-duration curve and basin characteristics at the ungaged 

location (see Fennessey (1994) for methodology; Waldron and 

Archfield (2006) provide an evaluation of and modifications to 

this methodology). 

To determine the firm yield for a reservoir or system 

of reservoirs, equation 1 is solved for each month of the 

historical record in chronological order. The yield, Qy, is 

initially set equal to zero and is incrementally increased. At 

each successive increase in yield, the monthly water balance 

is evaluated for all months of the simulation period. The yield 

is continually increased until a value is obtained that causes 

the reservoir to fail. When operating at this yield, the reservoir 

 Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
will not be able to meet the withdrawals in all months of the 

historical record; however, the previous yield value was able 

to meet the withdrawals for all months of the historical record. 

This yield—the previous value of Qy that had not resulted in 

reservoir failure—is considered the firm yield of the reservoir. 

Although the firm-yield model does not explicitly 

account for the contribution of ground water to reservoir 

storage, the QPPQ-transform method regression equations 

(Fennessey, 1994) used to determine the surface-water contribution to the reservoir include the soil-retention number, 

which is a measure of the amount of precipitation retained by 

the soils in the basin. The soil-retention number is determined 

by intersecting the drainage-basin boundary for the reservoir 

with maps of soil and land use. From this intersection, the 

percentage of each soil and land-use type in the drainage basin 

of the reservoir can be determined. Each combination of soil 

group and land cover is assigned a runoff-curve number (CN), 

which influences the expected runoff for a given amount of 

rainfall. A composite CN then is calculated by multiplying the 

fraction of the basin corresponding to a given CN by that number and adding the weighted CNs. Finally, the soil-retention 

number, in inches, is estimated from the CN by the use of the 

following equation (Natural Resources Conservation Service, 

1986):

Soil retention number = 1,000 10

CN - . (2)

Therefore, ambient ground-water flow to the reservoir that 

occurs as a function of soil retention is accounted for in the 

existing methodology to estimate the firm yield of reservoirs 

in Massachusetts; however, the contribution of ground 

water due to time-varying changes in reservoir stage is not. 

Approximately one-quarter of 70 selected drinking-water 

reservoirs in Massachusetts have shorelines that are in contact 

with sand and gravel aquifers. 

Determination of Ground-Water 

Contributions to Reservoir Storage

An analytical model of hydraulic interaction between a 

surface-water reservoir and a water-table (unconfined) aquifer 

was developed to calculate time-varying rates of ground-water 

inflow to, and storage losses from, surface-water reservoirs. 

The model is based on an analytical solution to the groundwater-flow equation for the case of one-dimensional groundwater flow in a finite-width aquifer bounded by a linear 

surface-water feature such as a stream (Rorabaugh, 1964). In 

this application, the surface-water feature is assumed to be a 

reservoir. The analytical model is a simplification of complex 

field conditions that provides approximate flow rates between 

ground-water and surface-water reservoirs in the absence of 

more detailed, time-consuming, and expensive, two- and threedimensional numerical ground-water-flow models.

The analytical solution derived by Rorabaugh (1964) to 

calculate the rate of flow between an aquifer and adjoining 

linear surface-water body for a constant rate of change in 

surface-water stage (h) is based on several assumptions, 

including:

The aquifer is of uniform width in the direction 

perpendicular to the surface-water body and extending to 

an impermeable boundary (or ground-water divide);

Aquifer hydraulic properties are isotropic and 

homogeneous; and

The surface-water feature that forms a boundary to the 

aquifer is straight, fully penetrates the aquifer, and is in 

direct hydraulic connection with the aquifer.

The second assumption requires that the transmissivity of 

the aquifer remains constant with time, which implies that, 

for a water-table aquifer, changes in water-table elevation are 

considered negligible when compared to the initial saturated 

thickness of the aquifer.

Because reservoir stage changes over time, it is necessary 

to use superposition to calculate total flow rate between 

the reservoir and aquifer as a function of time; that is, as a 

function of time-varying reservoir-stage changes. The use 

of superposition is mathematically acceptable here because 

the applied ground-water-flow equation is linear. In the 

superposition approach, incremental changes in the rate of 

flow between the reservoir and aquifer that occur in response 

to time-varying reservoir-stage changes are summed to 

estimate the total flow rate between the reservoir and aquifer. 

The calculation of total flow rate, QTOTi, between the 

reservoir and aquifer at present time ti is:

 

 

Q Q TOTi i

io

PRESENT PAST

t

Q j

j

= + =

å

1

, (3)

where

QPRESENTi is the flow rate at time ti between the 

reservoir and aquifer per unit length of 

reservoir shoreline resulting from the 

present change in reservoir stage that 

begins at time tio and ends at time ti , in 

units of square length per time; and

QPASTj is the flow rate between the reservoir 

and aquifer per unit length of reservoir 

shoreline resulting from a past change j in 

reservoir stage that begins at time tjo and 

ends at time tj, in units of square length 

per time.

1.

2.

3.

Determination of Ground-Water Contributions to Reservoir Storage 
The analytical solution to calculate the flow rate QPRESENTi

 resulting from the present change in reservoir stage that begins 

at time ti0 and ends at time ti is (Rorabaugh, 1964, eq. 5)

 

Q C aS

n

e PRESENT i a

n i

n T t t i io

= - - S

æ

è

ç

ç

ç

ç -

=

-

¥

1 å 8 1

2 2

2 2

2 4

1 3 

 ( )

, ... ç

ç

çç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷÷

,

 

(4)

where

QPRESENT i is the flow rate at time ti resulting from the present change in reservoir stage, in units of length

 squared per unit time per unit length of reservoir shoreline in contact with the aquifer;

Ci is the present rise or fall of the surface-water stage (hi – hio) divided by the time (ti – tio) over which the rise or fall

 has occurred, in units of length per time;

a is the perpendicular distance from the surface-water bank to the lateral boundary of the aquifer, in units of length;

S is the storage coefficient (or specific yield) of the aquifer, which is dimensionless; and

T is the transmissivity of the aquifer, in units of square length per time.

 The analytical solution to calculate QPASTj , the flow rate to the reservoir at present time ti resulting from a past change in

reservoir stage that began at time tjo and ended at time t

j

 is (Rorabaugh, 1964, eq. 7)

 

 

Q C aS

n

e PAST j n j

n T t t

a

i jo

= - - S

æ

è

ç

ç

ç

ç

ç

ç

- ( ) -

=

¥

1 å 8 1

2 2

2 2

2 4

1 3 



ç , ,...

ö

ø

÷

÷

÷

÷

÷

÷

÷

- -

æ

è

ç

ç

ç

ç

ç

çç

- ( ) -

=

¥

1 å 8 1

2 2

2 2

2 4

1 3 



n

e

n T t t

a

i j

S

n , ,...

ö

ø

÷

÷

÷

÷

÷

÷

÷

é

ë

ê

ê

ê

ê

ù

û

ú

ú

ú

ú

,

 

(5)

where

QPASTj is the flow rate at present time ti resulting from a past change j in reservoir stage ending at tj, in units of length

 squared per unit time per unit length of reservoir shoreline in contact with the aquifer;

Cj is the past rise or fall of the surface-water stage (hj – hjo) divided by the time (tj – tjo) over

 which the past rise or fall had occurred, in units of length per time;

j is the past time step;

a is the perpendicular distance from the surface-water bank to the lateral boundary of the aquifer, in units of length;

S is the storage coefficient (or specific yield) of the aquifer, which is dimensionless;

T is the transmissivity of the aquifer, in units of square length per time. 

The summation terms in equations 4 and 5 will approach zero when values of (T(t-to) / a2S) > 2.5; this can occur when t 

becomes relatively large (Rorabaugh, 1964). In the calculation of QPRESENTi (eq. 4), when the value of (T(ti-tio) / a2S) > 2.5, the 

flow between the aquifer and reservoir approaches the steady-state solution, CaS (Rorabaugh, 1964). In the calculation of QPASTj

(eq. 5), a past change Cj in reservoir stage will no longer contribute substantially to the calculation of QTOTi when 

(ti – tj) > (2.5a2S / T). For this case, the duration of time between the present time and the past change in reservoir stage 

becomes so large that the past change in reservoir stage has no substantial effect on the calculation of the flow rate between 

the reservoir and the aquifer at the present time. Calculation of QTOTi results in a flow rate per unit time per unit of reservoir 

shoreline in contact with the aquifer. Therefore, to determine the total volumetric flow rate entering the reservoir at time t

i

resulting from the present and past changes in reservoir stage, QTOTi is multiplied by the length of the reservoir shoreline in 

contact with the aquifer, L.

To maintain the assumption that flow between the reservoir and aquifer is horizontal and one-dimensional, the analytical 

method assumes that ground water is discharging along a stream channel with straight sides (Rorabaugh, 1964); however, few 

reservoirs have a shape that resembles a stream channel. To understand the effects of geometry on the one-dimensional solution 

to the ground-water-flow equation, ground-water-flow rates simulated by equation 3 were compared to ground-water-flow rates 

simulated by a two-dimensional, finite-difference ground-water-flow model for shapes of surface-water bodies that are a more 

realistic representation of a reservoir shoreline.

In the natural world, the map-view outlines of various ponds and lakes occur in a nearly infinite array of shapes; however, 

for comparative purposes, 12 simple, surface-water-body shapes were used to generalize the effect of reservoir-aquifer geometry 

on ground-water-flow rates (fig. 3A). The long-thin line (fig. 3A, shape no. 1) is the shape that is most representative of an idealized stream as described in the assumptions made in the derivation of the analytical solutions of Rorabaugh (1964); the remaining shapes represent various types of deviation from the idealized-stream shape (fig. 3A, shapes 2 to 12). 

 Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
A. SURFACE-WATER-BODY SHAPES

C. RESERVOIR STAGE USED IN THE NUMERICAL SIMULATIONS

1

3

6 7 9 10

4

2

5 8 11 12

0 2,000

0 500 1,000 METERS

4,000 FEET

B. MODEL GRID FOR RESERVOIR-AQUIFER SYSTEMS Rows

Columns

EXPLANATION 1

1 25 50 75 100 105

100

200

217

0 2,000

0 500 1,000 METERS

4,000 FEET

18.2 square 

miles

SURFACE-WATER-BODY SHAPE 10

MODEL GRID

BOUNDARY OF THE ACTIVE MODEL AREA

LENGTH OF RESERVOIR SHORELINE IN 

CONTACT WITH THE AQUIFER (L)

DISTANCE FROM THE RESERVOIR 

SHORELINE TO THE AQUIFER (a) 

L

a

EXPLANATION

9

44

45

46

47

48

49

50

0 60 120 180 240 300 360 420 480 540 600 660 720

TIME, IN DAYS

STAGE, IN FEET

SURFACE-WATER-BODY SHAPE 

AND IDENTIFIER—Sinusoidal 

time-varying hydraulic head 

fluctuation applied to area of shape 

using the Constant Head Boundary 

(CHD) Package in the numerical 

model simulation.

Figure 3. A, Surface-water-body shapes; B, with underlying model grid (shown with surface-water-body shape 10) for reservoiraquifer systems; and C, reservoir stage fluctuations used to evaluate the effects of reservoir-aquifer geometry on the analytical solution 

to the ground-water flow equation derived by Rorabaugh (1964).

Determination of Ground-Water Contributions to Reservoir Storage 
Numerical-Model Simulations of Idealized 

Reservoir-Aquifer Geometries

A finite-difference ground-water-flow model was used 

to calculate transient water budgets for idealized reservoiraquifer geometries. The finite-difference ground-water flow 

model used for these comparisons was MODFLOW-2000 

(Harbaugh and others, 2000). The idealized reservoir-aquifer 

geometries are different surface-water-body shapes (fig. 3A) 

incorporated into an underlying numerical model (fig. 3B) 

with identical hydrologic and spatial properties.

Each idealized reservoir-aquifer geometry (fig. 3B) was 

simulated using an identical underlying numerical-model grid 

covering an area of 18.2 square miles (fig. 3B). A transient, 

one-layer (50-ft thick) model consisting of 24 stress periods 

of 30 days in length (total simulated time of 720 days; two 

360-day years) was made to simulate the water budget for 

each of the 12 idealized-shape hydrologic systems (fig. 

3B) that include one of the shapes shown in figure 3A. The 

variable-spaced grid consists of 217 rows of increasing cellwidth (102.6 ft at center, increasing on the left to 396.8 ft 

and on the right to 408.8 ft) along the row and 105 columns 

of equal cell height (100.4 ft) down the column. A uniform 

value of 45 feet per day was specified for horizontal hydraulic 

conductivity. For each idealized-shape hydrologic system, 

the model was run under both confined and unconfined 

aquifer conditions. Values of storativity S were adjusted so 

that MODFLOW would calculate the same storage capacity 

and therefore produce similar water budgets when the model 

was run as a confined system compared to when the model 

run as an unconfined system. For unconfined conditions, a 

value for specific yield Sy is used by MODFLOW to calculate 

storage capacity. Here the value of Sy was set equal to 0.2 

(Sy = S = 0.2). For confined conditions, MODFLOW uses 

a value of specific storage Ss to calculate storage capacity. 

To determine a value for specific storage to yield a storage 

capacity similar to that for the unconfined case, the equation 

for a confined aquifer, S = bSs (Barlow and Moench, 1998; 

Fetter, 1994), where S equaled 0.2 and aquifer thickness b was 

50 ft. A Ss value of 0.004 was used for confined conditions. 

The computed flow rates between the reservoir and aquifer 

resulting from the unconfined simulations were used to 

evaluate the equations because there was no substantial 

difference between the unconfined and confined simulation 

results for the 12 idealized shapes.

The Constant Head Boundary (CHD) Package in 

MODFLOW-2000 (Harbaugh and others, 2000) was used to 

specify hydraulic head within the area corresponding to each 

different idealized surface-water-body shape (fig. 3A) for the 

length of the simulation. Numerical values that resemble a 

sinusoidal fluctuation when graphed against time were used 

to simulate the magnitude of the time-varying hydraulic head 

specified in the CHD input file (fig. 3C). The sinusoidal 

fluctuation, divided into segments corresponding to simulated 

stress periods 1 through 12 (year 1), had an initial value of 

47.0 ft, a high value of 49.5 ft, a low value of 44.5 ft, and a 

final value of 47.0 ft. The same head pattern was repeated for 

simulated stress periods 13 to 24 (year 2). No other stresses, 

such as pumping or recharge, were simulated. A value of 47.0 

ft was used elsewhere as starting head in the surrounding 

model cells.

Analytical Simulations of Idealized ReservoirAquifer Geometries

The idealized-shape hydrologic systems (fig. 3B) and 

2-year sinusoidal fluctuation in reservoir stage (fig. 3C) were 

used to analytically simulate flow rates between the reservoir 

and aquifer using the equations developed by Rorabaugh 

(1964). Although the numerical-model simulations track the 

flow of water between each cell over the entire simulated area 

to compute the total volumetric flow rate at the reservoiraquifer interface, the equations parameterize the geometry of 

the aquifer-reservoir geometries through the specification of 

L, the length of reservoir shoreline in contact with the aquifer 

and a, the distance from the reservoir shoreline to the aquifer 

boundary (fig. 3B). The values of L and a were estimated 

from the numerical-model grid (fig. 3B), where L was the 

perimeter of the idealized shape (because the entire shoreline 

of each surface-water body was in contact with the aquifer) 

and a was the average of several measured distances from the 

shoreline of the surface-water body to the aquifer boundary. 

The values of transmissivity T and storativity S used in the 

analytical simulations were equivalent to the values used in 

the numerical-model simulations: 2,250 feet squared per day 

(ft2/d) and 0.2, respectively. 

Effect of Reservoir Shape on the Flow Rate 

between the Reservoir and Aquifer

The timing and magnitude of the volumetric flow rates 

resulting from the numerical and analytical simulations 

of the reservoir-aquifer geometries were compared for 

the second year of the simulation period. The simulated 

volumetric flow rates were converted to units of million 

gallons per day and divided by L to obtain a volumetric 

flow rate in million gallons per day per foot of reservoir 

shoreline. Idealized shape 1 (fig. 3A)—the shape most 

resembling a stream channel—satisfies the assumptions 

used by Rorabaugh (1964) to derive the equations. For 

this case, the numerically and analytically simulated 

volumetric flow rates were nearly identical with respect to 

the timing and the magnitude of the flow rates (fig. 4). 

When the surface-water-body shape deviated from shape 

1 (fig. 3A), the difference between the timing of the numerically simulated and analytically simulated flow rates appear 

to be related to the ratio of a (the distance from the reservoir 

shoreline to the aquifer boundary) to L (the perimeter of the 

idealized shape) (fig. 4). As the ratio of a to L increased (L

smaller than a), the flow rates resulting from the numerical-

 Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
RATIO OF AQUIFER WIDTH TO THE LENGTH

OF RESERVOIR SHORELINE IN CONTACT 

WITH THE AQUIFER

181614121086420

DifferenceTime Lag

LAG (ANALYTICAL MINUS NUMERICAL) BETWEEN 

THE NUMERICAL-SIMULATION AND ANALYTICALSIMULATION RESPONSES, IN DAYS

PERCENT DIFFERENCE BETWEEN ANALYTICALLY 

AND NUMERICALLY SIMULATED MAXIMUM FLOW RATES 

RELATIVE TO THE NUMERICAL-SIMULATION RESULTS

360 420 480 540 600 660 720

360 420 480 540 600 660 720

0 1 2 3 4 5 6 7 8 Numerical Simulation

Analytical Formula

TIME, IN DAYS

TIME, IN DAYS

-10-15-20-25-30-35

-5

50

353025201510

360 420 480 540 600 660 720

0.0

-2.0-4.0-6.0-8.0

8.06.04.02.0

SEEPAGE RATE AT THE RESERVOIR-AQUIFER 

BOUNDARY, IN MILLION GALLONS PER DAY 

PER FOOT OF RESERVOIR SHORELINE X 10-5

0.0

-2.0-4.0-6.0-8.0

8.06.04.02.0

360 420 480 540 600 660 720

Numerical SimulationAnalytical Formula

Flow into aquifer

Flow into reservoir

Flow into aquifer

Flow into reservoir Numerical Simulation

Analytical Formula

Numerical SimulationAnalytical Formula

TIME, IN DAYS

SEEPAGE RATE AT THE RESERVOIR-AQUIFER 

BOUNDARY, IN MILLION GALLONS PER DAY 

PER FOOT OF RESERVOIR SHORELINE X 10-5

SURFACE-WATER-BODY SHAPE 1 SURFACE-WATER-BODY SHAPE 12

SURFACE-WATER-BODY SHAPE 2

SURFACE-WATER-BODY SHAPE 11

Flow into aquifer

Flow into reservoir

Flow into aquifer

Flow into reservoir

0.0

-2.0-4.0-6.0-8.0

8.06.04.02.0

SEEPAGE RATE AT THE RESERVOIR-AQUIFER 

BOUNDARY, IN MILLION GALLONS PER DAY 

PER FOOT OF RESERVOIR SHORELINE X 10-5

SEEPAGE RATE AT THE RESERVOIR-AQUIFER 

BOUNDARY, IN MILLION GALLONS PER DAY 

PER FOOT OF RESERVOIR SHORELINE X 10-5 0.0 -2.0 -4.0 -6.0 -8.0 8.0 6.0 4.0 2.0

SURFACE-WATER-BODY SHAPE 9

TIME, IN DAYS

Figure 4. Differences between the numerically and analytically simulated ground-water-flow rates, and timing of the flows for selected reservoir-aquifer geometries shown in figure 3.

Determination of Ground-Water Contributions to Reservoir Storage 
model simulations responded more slowly to the sinusoidal 

time-varying changes in reservoir stage than the resulting flow 

rates from the analytical simulations (fig. 4). 

With the exception of idealized shapes 2, 9, 11, and 

12, the differences in magnitude between the numerically 

simulated and analytically simulated volumetric flow rates 

were within ±10 percent (fig. 4, center graph). For values of 

a to L less than 1 (L larger than a), the analytically simulated 

flow rates were greater than the numerically simulated flow 

rates; for values of a to L greater than 1 (L smaller than 

a), the analytically simulated flow rates were less than the 

numerically simulated flow rates. Selected reservoir-aquifer 

geometries are shown in figure 4 to illustrate the differences 

between the flow rates determined from the numerical and 

analytical simulations. 

The four idealized shapes that had a greater than ±10-

percent difference between the two volumetric flow rates also 

were the reservoir-aquifer geometries with the most extreme 

values of a to L. Furthermore, the largest percent differences 

do not appear to be related to the shape of the reservoir 

shoreline but, rather, the values of a to L that resulted from 

the reservoir-aquifer geometry: idealized shape 2 is an ellipse, 

idealized shape 9 is a square, and idealized shapes 11 and 12 

are circles (fig. 3A and fig. 4). Reservoir-aquifer geometries 

with values of a to L between approximately 0.5 and 3.5 have 

simulated flow rates within ±10 percent (fig. 4).

Differences between the magnitude and timing of 

the flow rates resulting from the numerical and analytical 

solutions were expected because of the inherent differences 

between the numerical and analytical simulations. The 

equations (analytical simulations) assume a one-dimensional, 

horizontal flow, whereas the numerical simulations account for 

two-dimensional flow between the surface-water-body shapes 

and an aquifer. Additionally, the equations represent only the 

reservoir-aquifer geometry through the parameters a and L, 

and only the time-varying changes in reservoir stage determine 

the timing of the flow rates. Conversely, the numerical 

simulation accounts for the flow of water between each cell 

in the model grid (fig. 3B), and the timing and magnitude 

of the flow rate is not only a function of the time-varying 

changes in reservoir stage but also the integration of the 

fluxes between all cells in the model grid. Despite the simple 

geometric parameterization of the reservoir-aquifer interaction 

and the assumption of one-dimensional flow in the analytical 

simulations, the analytically simulated flow rates were 

reasonably comparable to the numerically simulated flow rates 

with respect both to the magnitude and timing of the flow rates 

for most of the idealized shapes. Furthermore, the shape of the 

reservoir did not appear to affect the agreement between the 

numerically simulated and analytically simulated flow rates; 

rather, the size of aquifer (represented by a) relative to the size 

of the reservoir (represented by L) appears to explain much 

of the variability in the agreement between the numerically 

simulated and analytically simulated flow rates for the various 

reservoir-aquifer geometries.

Application of the Ground-Water-Flow 

Equations to the Firm-Yield-Estimator 

Model 

Although the analytical equations developed by 

Rorabaugh (1964) were derived by assuming one-dimensional, 

horizontal flow between an aquifer and a stream, the equations 

provide a reasonable approximation of numerically simulated flow rates between a reservoir and aquifer (within ±10 

percent) when the ratio of a to L is between approximately 

0.5 and 2.3, regardless of the shape of the reservoir shoreline 

in contact with the aquifer (figs. 3A and 4). Additionally, the 

analytical equations require only five inputs to estimate the 

flow rate between the reservoir and aquifer: (1) a, the average 

distance from the reservoir shoreline to the aquifer boundary 

(fig. 3B); (2) L, the length of the reservoir shoreline in contact 

with the aquifer (fig. 3B); (3) T, aquifer transmissivity; (4) S, 

aquifer storativity; and (5) the relation between reservoir storage and reservoir stage, which is used to determine the timevarying changes in reservoir stage. Therefore, the analytical 

equations developed by Rorabaugh (1964) were incorporated 

into the existing FYE model. 

The modified FYE model was incorporated into the existing relational database structure developed by Waldron and 

Archfield (2006). The generalized framework for the modified 

FYE model is shown in figure 5. For the case where a reservoir is not in contact with an aquifer, the FYE model remains 

unchanged. When a reservoir is in contact with an aquifer, the 

five inputs must be estimated in addition to the data requirements of the original FYE model.

The reservoir water balance, including the ground-water 

contribution to reservoir storage due to time-varying changes 

in reservoir stage, can be expressed as

Si = Si-1 + Awi Qri + Ari (Pi – Ei) ± 

QGWi- αi ni Qy – Qri – Qsi – Qoi , (6) 

where

QGWi is the volume of water that moves between 

the reservoir and aquifer due to the present 

change in reservoir stage during the month 

and the time varying past changes in 

reservoir stage over previous months, in 

units of million gallons.

Although the change in reservoir stage over month i is 

assumed to be constant with time, the resulting flow rate is 

not (eqs. 4 and 5). The flow rate is a function of the difference 

between ti and tio, and time ti and tjo, which will be different for 

each point in time within the month. Therefore, the inclusion 

of the term QGWi in the water balance of a given reservoir 

implies that the resulting volumetric term QGWi is the integration of the flow rates between the reservoir and aquifer over all 

time t within the present month i:

10 Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
 

 Q Q GW TOT

t t

i dti

i io

=

-

ò

0

,

 (7)

where

ti is the number of days from the start of the simulation to the end of the present change in reservoir stage, and

tio is the number of days from the start of the simulation to the beginning of the present change in

 reservoir stage.

Equation 7 is an integration of equation 3 and requires the calculation of the ground-water flow resulting from all past timevarying changes as well as the present change in reservoir stage at all values of t throughout month i. A closed-form, analytical 

solution to equation 7 does not exist; hence, it is not possible to directly implement equation 7 into the FYE model. For this 

reason, an approximation to equations 3 through 5 and equation 7 was required to determine the flow rate during the month.

An additional challenge to the implementation of the equations arises from the feedback between reservoir storage and 

reservoir stage. To determine Ci and, in turn, QGWi, the reservoir stage at the end of the present month must be known; however, 

because the reservoir stage is estimated from the reservoir storage-stage relation, the stage in the reservoir is not known until the 

FYE model calculates the reservoir storage at the end of the present month. Yet, the reservoir storage at the end of the present 

month cannot be calculated until the values of the water-balance terms—including QGWi—are determined (eq. 5). An iterative 

loop was added to the modified FYE model to address this within-month feedback between reservoir storage, reservoir stage, 

and QGWi.

Modifications to the Ground-Water-Flow Equations

To approximate the equations that calculate QGWi, QPRESENTi and each QPASTj are computed at the middle of the month (t

approximately equal to 15 days). This flow rate is assumed to represent the average flow rate between the reservoir and aquifer 

during the present month i. The modified equations for QPRESENTi, QPASTj, and QGWi are 

 

Q C aS

n

e PRESENT i

n T t t

a S

i

i io

=- - -

æ -

è

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷÷

1 8 1

2 2

2

4

2 2

2





n=

¥

å

æ

è

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷÷

1 3, ,... ,

 

(8)

 

Q C aS

n

e PAST j

n T t t t

j

i i io

=- - -

-

æ -

è

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷÷

æ

è

ç

ç

1 8 1

2 2

2 2 2



 ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷÷

æ

è

ç

ç

ç

ç

ç

ç

ç

çç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷

÷÷

-

=

t

a S

jo

n

4 2

1 3, ,...

¥

å

æ

è

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷

÷÷

- - -

-

1 8 1

2 2

2 2





n

e

n T t

i

t t t i i j o - - æ

è

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷÷

æ

è

ç

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷÷

æ

è

ç

ç

ç

ç

ç

ç

ç 2 çç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷

÷÷

=

¥

å

æ

è

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ç

ö

ø

÷

÷

÷

÷

÷

÷

÷

÷

÷ 4 2

1 3

a S

n , ,... ÷

÷

÷

÷

÷

÷

÷÷

é

ë

ê

ê

ê

ê

ê

ê

ê

ê

ê

ù

û

ú

ú

ú

ú

ú

ú

ú

ú

ú

,

(9)

 and Q L Q t t GWi TOT i i = o ( )( )( ) - . (10)

To verify that the modified equations were applied correctly to the FYE model, numerical-model simulations of flow 

rates along with reservoir water levels were compared to FYE-simulated flow rates and reservoir water levels for four idealized hydrologic systems patterned after four reservoir-aquifer systems in Massachusetts (fig. 6): Reservoir A, patterned after 

Fitchburg Reservoir in Ashby; Reservoir B, patterned after Morse Reservoir in Leominster; Reservoir C, patterned after Bearhole Reservoir in West Springfield; and Reservoir D, patterned after Millham Reservoir in Marlborough. A 24-month (2-year) 

period—April 1975 through March 1977—that represented normal hydrologic and climatic conditions was used to simulate flow 

rates and water levels in the reservoirs. A 2-year simulation provided sufficient change in reservoir storage to compare the simulated ground-water-flow rates. Real reservoir-aquifer systems provided input variables that represent values in Massachusetts; 

however, because the simulations are an oversimplified representation of the actual reservoir-aquifer geometry, the numerically 

and FYE-simulated flow rates and water levels were not intended to reproduce actual flow rates and water levels that occurred 

during the simulation period. 

Application of the Ground-Water-Flow Equations to the Firm-Yield-Estimator Model 11
Are the total 

number of 

failures > 1?

For each month i

of the simulation 

period

Is the reservoir in 

hydrologic connection 

with sand and gravel?

Output firm-yield 

(current yield 

value minus the 

increment value)

Increment

yield

Is this the last 

month of the 

simulation period?

Total the number 

of reservoir 

failures that 

occurred over the 

simulation period.

Solve for the storage in the 

reservoir Si at the end 

of month i using the known 

values in the water-balance 

equation.

Is the difference 

between Ci

n+1

 and Ci

n less

than the 

closure-criteria

value?

NO YES Next

month

NO YES

YES

NO

Add QGWi

n and Si-1 to obtain

new estimate of storage Si

n+1

in the reservoir at the end of

month i.

NO YES

Determine if the 

reservoir failed 

or overflowed.

If the reservoir overflowed:

1. Calculate overflow.

2. Set reservoir storage to 

 maximum capacity.

If the reservoir failed:

1. Activate failure flag.

2. Set reservoir storage to 

 zero.

Use the hypsographic curve 

to convert reservoir storages 

Si

n and Si-1 to reservoir 

stages hi

n

 and hi-1.

Calculate the quantity of water 

that has moved between the

reservoir and the aquifer,

QGWi

n

 resulting from the 

change in reservoir stage

over the month, Ci

n

.

Calculate

hi

n - hi-1 Ci

n

 = timei - timei-1

Does Ci

n

 = 0?

YES

Calculate new Ci

n+1

:

hi

n+1

 - hi-1 Ci

n+1

 = timei - timei-1

Use the hypsographic curve 

to convert reservoir 

storages Si

n+1

 and Si-1 to

reservoir stages 

hi

n+1 and hi-1.

NO

Set

n = 1

Set

n = n +1

EXPLANATION

reservoir storage at the end of the previous month

reservoir stage at the end of the previous month

number of estimates resulting from the iterative loop

nth estimate of the ground-water contribution resulting 

from the iterative loop

nth estimate of the change in reservoir stage over the 

month resulting from the iterative loop

nth estimate of reservoir stage resulting from the 

iterative loop

nth estimate of reservoir storage resulting from the 

iterative loop

n

hi

n

Ci

n

QGWi

n

Si

n

hi-1

Si-1

Figure 5. Flow chart of the modified firm-yield-estimator model.

12  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
Williams

Lake

Millham

Reservoir

MARLBOROUGH

NORTHBOROUGH

BERLIN

HUDSON

SOUTHBOROUGH

71°36' 71°34'48"

42°19'12"

42°21'36"

0 1 MILE

0 1 KILOMETER

Bearhole

Reservoir

Ashley

Pond

WEST SPRINGFIELD

WESTFIELD

HOLYOKE

0 1 MILE

0 1 KILOMETER

72°42' 72°39'

42°07'12"

42°10'12"

Haynes

Reservoir

Morse

Reservoir

LEOMINSTER

71°48' 71°46'

42°30'

42°32'

Ashby

Reservoir

Fitchburg

Reservoir

ASHBY

FITCHBURG

ASHBURNHAM

71°51'36" 71°49'12"

42°38'24"

42°39'36"

0 .5 MILES

0 .5 KILOMETERS

0 .5 MILES

0 .5 KILOMETERS

73° 71°

42°30'

41°30'

0 30 60 MILES

0 40 80 KILOMETERS

Base from U.S. Geological Survey digital data, 1:25,000

Lambert conformal conic projection. 1983 North American Datum

Massachusetts coordinate system mainland zone

EXPLANATION

SURFACE-WATER RESERVOIR

SAND AND GRAVEL AREA

TOWN BOUNDARY

SUBBASIN BOUNDARY

LENGTH OF RESERVOIR

SHORELINE IN CONTACT 

WITH THE AQUIFER 

A. B.

C.

D.

Figure 6. Locations of four reservoir-aquifer systems in Massachusetts: A, Fitchburg Reservoir; B, Morse Reservoir; C, 

Bearhole Reservoir; and D, Millham Reservoir.

Application of the Ground-Water-Flow Equations to the Firm-Yield-Estimator Model 13
Numerical-Model Simulations of Ground-WaterFlow Rates and Reservoir Water Levels

Reservoirs A through D were simulated by a transient 

one-layer model, and the layer type was specified as convertible between confined and unconfined. If a simulated head was 

above the top of the layer, it functioned as a confined layer; if 

the simulated head was below the top of the layer, it functioned as an unconfined layer. The reservoirs were simulated 

as a high transmissivity zone. Model grids were composed 

of cells 50 ft square and were developed to simulate closely 

the physical volumes of the reservoir and the sand and gravel 

aquifer of each respective system. The 50-ft square dimension 

also provided reasonable numerical solution times. 

All terms in the water-balance equation (eq. 6), with the 

exception of QGWi, were applied uniformly within each stress 

period. Monthly stress periods corresponded to the actual 

number of days per given month (28 to 31 days). To simplify 

the application of stresses to the model, the sum of the known 

inflows to the reservoir was applied as an injection well in the 

reservoir area, and the sum of the known outflows was applied 

as a pumping well in the reservoir area (fig. 7). 

To ensure the geometry did not affect the comparison 

between the simulations, the aquifer and reservoir were represented as rectangles (large streams) that extended the length of 

the aquifer, just as the idealized stream shape (fig. 3A, shape 

no. 1). The length of the reservoir shoreline was equal to the 

length-dimension of the aquifer, having length equal to the 

approximate length of the actual reservoir shoreline in contact 

with the aquifer (fig. 7). This distance was estimated from 

Geographic Information System (GIS) digital data layers of 

hydrography in Massachusetts at a 1:100,000 scale intersected 

with GIS digital data layers of surficial geology in Massachusetts at a 1:250,000 scale (Office of Geographic and Environmental Information (MassGIS), 2004a; Office of Geographic 

and Environmental Information (MassGIS), 2004b). The 

surface area and capacity of the reservoir were obtained from 

Waldron and Archfield (2006). To determine the width of the 

reservoir area (fig. 7), the surface area of the reservoir was 

divided by the length of the reservoir shoreline in contact with 

the aquifer (table 1). The capacity of the reservoir was divided 

by the surface area to obtain a depth for the reservoir (table 

1). Therefore, idealized reservoir dimensions were used in the 

simulations. 

The surface area of the aquifer also was estimated from 

GIS digital data layers of hydrography in Massachusetts at 

a 1:100,000 scale intersected with GIS digital data layers of 

surficial geology in Massachusetts at a 1:250,000 scale (Office 

of Geographic and Environmental Information (MassGIS), 

2004a; Office of Geographic and Environmental Information 

(MassGIS), 2004b). The width of the aquifer area was determined by dividing the surface area of the aquifer by the length 

of the reservoir shoreline in contact with the aquifer. The 

analytical equations assume that the reservoir fully penetrates 

the aquifer (Rorabaugh, 1964); therefore, the thickness of the 

aquifer was specified as equivalent to the depth of the reservoir in the simulations (table 1), although this assumption may 

be not true for all real reservoir-aquifer systems.

Hydrologic atlases were used to estimate an aquifer 

transmissivity for reservoirs B, C, and D (Brackley and 

Hanson, 1977). Because no other hydrologic information was 

available for the aquifers associated with reservoirs A, B, and 

C, specific yield values for the reservoir and aquifer were 

1 and 0.2, respectively. Due to the simulated geometry and 

inflow and outflow volumes applied to the reservoir, a high 

hydraulic conductivity value (5,000,000 ft/d) was applied to 

reservoirs A through D to simulate a flat water surface over 

the area of the reservoir. Reservoir D is located in an area for 

which a detailed ground-water-flow model has been developed 

and, therefore, recent estimates of the transmissivity and 

specific yield of the aquifer were available (DeSimone, 2004). 

Specific-yield values of 1 and 0.18 were specified for reservoir 

D and the aquifer, respectively (table 1). 

Modified Firm-Yield-Estimator Model Simulations 

of Ground-Water-Flow Rates and Reservoir 

Water Levels

Reservoirs A through D were simulated over the same 

2-year period (April 1975 through March 1977) as the numerical-model simulations. The relation between reservoir storage 

and reservoir stage was determined assuming a rectangularshaped reservoir; the FYE-model simulations represented 

the surface area of the reservoir as constant for all values of 

reservoir storage, and the relation between reservoir stage and 

reservoir storage as linear. The values of aquifer transmissivity 

T, aquifer storage coefficient S, aquifer width a, and length of 

reservoir shoreline in contact with the aquifer L were equivalent to the values used in the numerical-model simulations 

(table 1).

Comparison of Simulated Ground-Water-Flow 

Rates and Reservoir Water Levels

Average ground-water-flow rates and water levels 

resulting from the numerical and FYE-model simulations 

are compared in figure 8. The FYE simulation calculates the 

flow rate between the reservoir and aquifer at the middle of 

the month and assumes this flow rate represents the average 

flow rate over the entire month. The numerical simulation 

calculates the flow rate between the reservoir and aquifer at 14 

unequal time steps within each month. Therefore, the average 

flow rate between the reservoir and aquifer was computed for 

each month of the numerical simulation and compared to the 

flow rate simulated by the FYE model (fig. 8A). The water 

levels calculated by the FYE simulation are the levels that 

result from the stresses applied during the month, including 

the inflow or outflow of ground water that has taken place. 

Therefore, the resulting water level represents the water level 

14  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
1

1 25 50 75 100 122

2550

56

RESERVOIR D

1

1 25 37 1 13

2550

58

1

2539

injection well

Sand and gravel aquifer

pumping well

 RESERVOIR A RESERVOIR B

inflow

release

diversion 0

0 100 200 METERS

500 1,000 FEET

1

1 25 50 75 80

255075

100

125129

inflow

release

diversion

 RESERVOIR C 

injection well

pumping well

pumping well

injection well

injectionwellpumpingwell

EXPLANATION

Sand and gravel aquifer

Reservoir

Length of reservoirshoreline in contact withthe aquifer

columns

rows

columns

rows

columns

rows

columns

rows

Figure 7. Model grids used in the numerical simulations of four idealized reservoir-aquifer systems patterned after the four reservoir-aquifer systems shown in figure 6. Reservoir A is patterned after Fitchburg Reservoir (fig. 6

A); Reservoir B is patterned after Morse Reservoir (fig. 6

B); Reservoir C is patterned after Bearhole Reservoir (fig. 6

C); 

Reservoir D is patterned after Millham Reservoir (fig. 6

D).

Application of the Ground-Water-Flow Equations to the Firm-Yield-Estimator Model 15
Table 1. Reservoir and aquifer characteristics for four idealized reservoir-aquifer systems.

Idealized reservoir-aquifer system

Characteristics Reservoir A Reservoir B Reservoir C Reservoir D

Total number of rows in numerical-model grid 58 39 129 56Total number of columns in numerical-model grid 37 13 80 122Actual length of reservoir shoreline in contact with 

sand and gravel aquifer materials, in feet

2,900 1,950 6,450 2,800

Calculated actual volume of reservoir, in cubic feet 91,630,000 6,110,000 12,370,000 41,890,000Calculated actual volume of aquifer, in cubic feet 13,500,000 19,400,000 481,700,000 277,200,000Simulated length and width of reservoir, in feet 2,900 (58 cells), 

1,600 (32 cells)

1,950 (39 cells),

 150 (3 cells)

6,450 (129 cells),

100 (2 cells)

2,800 (56 cells),

800 (16 cells)

Simulated length and width of aquifer, in feet 2,900 (58 cells), 

250 (5 cells)

1,950 (39 cells),

500 (10 cells)

6,450 (129 cells),

3,900 (78 cells)

2,800 (56 cells),

5,300 (106 cells)

Simulated depth of aquifer and reservoir, in feet 19.75 20.89 19.18 18.70Simulated volume of aquifer, in cubic feet 14,318,750 20,367,750 482,472,900 277,508,000Transmissivity of sand and gravel aquifer, in feet 

squared per day

1,350 (Brackley and 

Hanson, 1977)

1,350 (Brackley and 

Hanson, 1977)

1,400 (Maevsky and 

Johnson, 1990)

1,800 (DeSimone, 2004)

Horizontal hydraulic conductivity, in feet per day reservoir = 5,000,000, 

aquifer = 68.36

reservoir = 5,000,000,

aquifer = 64.62

reservoir = 1E+8, 

aquifer = 74.86

reservoir = 5,000,000, 

aquifer = 96.25

Specific yield reservoir = 1, aquifer = 0.2 reservoir = 1, 

aquifer = 0.2

reservoir = 1, 

aquifer = 0.2

reservoir = 1,

aquifer = 0.18

Specific storage reservoir = 1, 

aquifer = 0.00001

reservoir = 1,

aquifer = 0.00001

reservoir = 1,

aquifer = 0.00001

reservoir = 1,

aquifer = 0.00001

16  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
in the reservoir at the end of the month and is compared to the 

water level from the final time step in the numerical simulation 

(fig. 8B).

The two simulated, average ground-water-flow rates and 

reservoir water levels show close agreement (fig. 8). The largest differences between the numerically and FYE-simulated 

average flow rates were less than 0.04 million gallons per day 

(Mgal/d) and are present at the peaks and troughs (fig. 8). The 

largest differences between the numerically and FYE-simulated water-level were less than 0.07 ft. For the majority of the 

simulation period, the differences between simulated water 

levels and simulated flow rates were less than 0.01 ft and 0.01 

Mgal/d, respectively (fig. 8). 

The comparisons between numerical and FYE-model 

simulations of ground-water-flow rates between the reservoir 

and aquifer and reservoir water levels appear to indicate that 

the modified Rorabaugh (1964) equations were incorporated 

into the FYE model appropriately. Furthermore, the results of 

the comparisons indicate the calculation of the total monthly 

contribution of ground water (assuming a constant, average 

ground-water-flow rate) compares well with numerically simulated, average ground-flow rates.

Iterative Loop to Calculate the Change in 

Reservoir Stage

The FYE model computes the reservoir storage at the 

end of each month in the simulation to determine if the usable 

storage in the reservoir is able to meet the demand to the 

reservoir. The tracking of reservoir storage at successive yields 

is ultimately how the firm yield of the reservoir is determined. 

To convert reservoir storage to reservoir stage, a bathymetric 

survey of the reservoir can be used to establish a continuous stage-storage relation from which one can interpolate a 

reservoir stage from a known reservoir storage. To enable the 

FYE model to track reservoir stage, the modified FYE model 

requires the entry of polynomial coefficients that provide a 

continuous function describing the stage-storage relation. 

At the beginning of the simulation, the FYE model 

assumes the reservoir storage and stage is at a maximum (Massachusetts Department of Environmental Protection, 1996). As the simulation steps forward in time, past 

changes in reservoir stage will contribute to the calculation of 

Si (eq. 9). Any past change in reservoir stage will have already 

been determined during the month when that change took 

place. The present ground-water contribution resulting from 

this past change in stage is only dependent on the duration of 

time that has passed between the past change and the present time. For this reason, the volume of water that exchanges 

between the reservoir and aquifer due to past changes in 

reservoir stage, ΣQPASTj, in present month i (eq. 9) can be 

determined outside of the iterative loop. Therefore, only Ci

and the resulting QPRESENTi terms in the ground-water-flow 

equations that determine QGWi are dependent upon the results 

of the iterative loop. 

To provide the initial estimate of the change in reservoir 

stage Ci

1, where the superscript denotes the estimate of 

Ci obtained from the nth pass (in this case, the first pass) 

through the iterative loop that has occurred during the present 

month i in the simulation, equation 6 is solved for Si

1 with 

consideration for the QPASTi terms in the equation for QGWi

(eq. 10) but not the value of QPRESENTi because it has yet to be 

determined. Through the storage-stage relation, the calculation 

of Si

1 is converted to a corresponding reservoir stage hi1. The 

value of hi

1 is then used to compute an initial estimate of Ci1

(fig. 5). Using the estimate of Ci

1, an estimate of QPRESENTi1

is obtained and used to re-solve the water-balance equation 

(eq. 6). Re-solving equation 6 will result in a new estimate of 

reservoir storage, Si

2, which will result in new estimates of hi2, 

Ci

2, and QGWi2. The iterative process repeats until the previous 

estimate of the stage change, Ci n-1, agrees with the new 

estimate of the stage change, Ci n, within a specified closurecriteria value (fig. 5). The closure criterion on the iterative 

loop is placed on Ci because it is the change in reservoir stage 

that controls the volume of water that exchanges between 

the reservoir and aquifer (eq. 10). When the closure criterion 

is met, the final estimate of Si is evaluated to determine if 

the usable reservoir capacity has been completely depleted 

and failure occurred. Then the FYE model steps forward to 

evaluate the next month.

The closure criterion determines the precision to 

which consecutive estimates of Ci n agree within the iterative loop of the FYE model. As the iterative loop converges 

towards a value of Ci, the closure criterion specifies when 

the values of Ci n and Ci n-1 are close enough to provide a 

reasonable estimate of the value of Ci , and in turn, a reasonable estimate of QPRESENTi

. As the value of the closure 

criterion decreases, successive estimates of Ci are required 

to be known to greater precision before the iterative loop 

will finish; however, as the closure criterion decreases, 

computational time will increase. The choice of too large a 

value for the closure criterion could result in large over- or 

underestimates of QPRESENTi for each month of the simulation, and, in turn, a large over- or underestimate of the firm 

yield. Alternatively, the choice of too small a value for the 

closure criterion will result in increased computational time.

To ensure firm-yield estimates are not affected by the 

choice of the closure-criterion value for the iterative loop, the 

difference in firm-yield estimates was compared for closure 

criterion values of 0.1 ft/d, 0.01 ft/d, 0.001 ft/d, and 0.0001 

ft/d. The firm yield at each closure criterion was computed 

for reservoirs A through D using the same 2-year period of 

stresses, reservoir and aquifer geometries, and reservoir and 

aquifer characteristics as was used to compare results from the 

numerical and FYE-model simulations (table 1). 

As the closure criterion decreased, the estimate of 

QPRESENTi was known to greater precision; however, as the 

closure criterion became more restrictive (less than 0.01 ft/d), 

the firm yield did not change within a value of 1,000 gallons 

(gal), even when the closure-criterion value was decreased 

further (fig. 9). Differences between firm yield estimates 

Application of the Ground-Water-Flow Equations to the Firm-Yield-Estimator Model 17
Figure 8. A, Comparison of simulated average ground-water flow rates; and 

B, reservoir water-level elevations for four idealized reservoir-aquifer systems during a 

2-year period of average climate and streamflow conditions (April 1975 through March 1977).

673.0672.5672.0671.5671.0670.5670.0669.5669.0668.5668.0667.5667.0666.5666.0665.5665.0

999.5

1007.51007.01006.51006.01005.51005.01004.51004.01003.51003.01002.51002.01001.51001.01000.51000.0

-0.02-0.04-0.06-0.08-0.10-0.12-0.14-0.16-0.18-0.20

0.200.180.160.140.120.100.080.060.040.020.00

-0.02-0.04-0.06-0.08-0.10-0.12-0.14-0.16-0.18-0.20

0.200.180.160.140.120.100.080.060.040.020.00

RESERVOIR A

RESERVOIR B

Firm-Yield-Estimator simulationNumerical simulation

Flow into reservoir

Flow into aquifer

Mar

'77

Apr

1975

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Mar

'77

Apr

1975

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Flow into reservoir

Flow into aquifer

Mar

'77

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Apr

'77

Mar

'77

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Apr

'77

SIMULATED RESERVOIR WATER-LEVEL 

ELEVATION, IN FEET

SIMULATED AVERAGE FLOW RATE, 

IN MILLION GALLONS PER DAY

SIMULATED AVERAGE FLOW RATE, 

IN MILLION GALLONS PER DAY

SIMULATED RESERVOIR WATER-LEVEL 

ELEVATION, IN FEET

Firm-Yield-Estimator simulationNumerical simulation

Firm-Yield-Estimator simulationNumerical simulation

Firm-Yield-Estimator simulationNumerical simulation

A. B.

18  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
Figure 8. A, Comparison of simulated average ground-water flow rates; and 

B, reservoir water-level elevations for four idealized reservoir-aquifer systems during a 

2-year period of average climate and streamflow conditions (April 1975 through March 1977).—Continued

236.5236.0235.5235.0234.5234.0233.5233.0232.5232.0231.5231.0230.5230.0229.5229.0228.5

165.1164.6164.1163.6163.1162.6162.1161.6161.1160.6160.1159.6159.1158.6158.1157.6157.1

-0.02-0.04-0.06-0.08-0.10-0.12-0.14-0.16-0.18-0.20

0.200.180.160.140.120.100.080.060.040.020.00

-0.02-0.04-0.06-0.08-0.10-0.12-0.14-0.16-0.18-0.20

0.200.180.160.140.120.100.080.060.040.020.00

Firm-Yield-Estimator simulationNumerical simulation

RESERVOIR C

Flow into reservoir

Flow into aquifer

RESERVOIR D

Mar

'77

Apr

1975

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Firm-Yield-Estimator simulationNumerical simulation

Flow into reservoir

Flow into aquifer

Mar

'77

Apr

1975

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Mar

'77

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Apr

'77

Mar

'77

May

'75

Jun

'75

Jul

'75

Aug

'75

Sep

'75

Oct

'75

Nov

'75

Dec

'75

Feb

'77

Jan

'77

Apr

'76

May

'76

Jun

'76

Jul

'76

Aug

'76

Sep

'76

Oct

'76

Nov

'76

Dec

'76

Mar

'76

Feb

'76

Jan

'76

Apr

'77

Firm-Yield-Estimator simulationNumerical simulation

Firm-Yield-Estimator simulationNumerical simulation

A. B.

WATER-LEVEL SIMULATED RESERVOIR 

ELEVATION, IN FEET

SIMULATED AVERAGE FLOW RATE, 

IN MILLION GALLONS PER DAY

SIMULATED AVERAGE FLOW RATE, 

IN MILLION GALLONS PER DAY

SIMULATED RESERVOIR WATER-LEVEL 

ELEVATION, IN FEET

Application of the Ground-Water-Flow Equations to the Firm-Yield-Estimator Model 19
were observed between closure-criteria values of 0.1 ft/d and 

0.01 ft/d as well as between closure-criteria values of 0.01 

ft/d and 0.001 ft/d; however, the difference between the firm 

yield determined from closure-criteria values of 0.001 ft/d and 

0.0001 ft/d were less than 0.001 Mgal/d for all idealized reservoir-aquifer systems (fig. 9). Therefore, a closure criterion of 

of 0.001 ft/d was used in the FYE model.

Effect of Uncertainty in the Input 

Variables to the Ground-Water-Flow 

Equations

Comparisons of ground-water-flow rates and reservoir 

water levels for idealized reservoir-aquifer systems indicated 

that the ground-water-flow equations were applied correctly to 

the modified FYE model. To apply the modified FYE model to 

real reservoir-aquifer systems in Massachusetts, one will need 

to estimate the input variables that are required to determine 

the ground-water contribution. Because the estimates of these 

input variables are likely to differ from the true values, the 

sensitivity of firm yield to the user-specified reservoir-aquifer characteristics (stage-storage relation, a, L, T, and S) was 

evaluated to quantify differences in firm yield. 

Sensitivity of Firm Yield to the Reservoir StageStorage Relation

Although the stage-storage relation for reservoirs A 

through D was simulated assuming a rectangular-shaped 

reservoir, the actual shape of a reservoir and the corresponding 

stage-storage relation is complex. Before the firm yield of a 

reservoir can be calculated, detailed bathymetric data are used 

to determine the stage-storage relation for the usable capacities 

of the reservoirs. For reservoirs A through D, a parsimonious, 

monotonically increasing polynomial was fit to the bathymetric data collected by Waldron and Archfield (2006) to provide 

a continuous relation between reservoir storage and stage. 

Using this relation, the firm yield was then calculated from the 

2-year simulation period (April 1975 through March 1977) 

for reservoirs A through D and compared to the firm yield 

determined from a stage-storage relation that was based on a 

rectangular-shaped reservoir.

Figure 9. Sensitivity of firm yield to the closure criterion used in the iterative loop of 

the firm-yield simulations for four idealized reservoir-aquifer systems.

-0.05

-0.04

-0.03

-0.02

-0.01

0

0.01

Reservoir A

Reservoir D

Reservoir B

Reservoir C DIFFERENCE IN FIR

M

YIELD, IN MILLION G

ALLONS PER DA

Y

Difference in firm yield between

closure criteria of 0.1 foot/day

and 0.01 foot/day.

Difference in firm yield between

closure criteria of 0.01 foot/day

and 0.001 foot/day.

Difference in firm yield between

closure criteria of 0.001 foot/day

and 0.0001 foot/day.

20 Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
During simulations of reservoirs C and D, the iterative 

loop did not converge when water levels approached the bottom of the reservoir. Near the bottom of the reservoir, large 

changes in reservoir stage induced an unrealistic ground-water 

flux that overwhelmed the values of the other terms in the 

water balance (fig. 10). This large overestimate in the groundwater flux prevented the iterative loop from converging or, 

in some cases, caused the iterative loop to converge at an 

unacceptably slow rate. Large changes in stage near the bottom of the reservoir are typical of reservoirs in Massachusetts 

because most drinking-water reservoirs were constructed on 

an existing stream. Therefore, the deepest part of the reservoir 

bottom has steeply sloping sides where the stream channel was 

located, and the upper part of the reservoir bottom is flat with 

gently sloping sides, representing the location of the former 

flood plain (fig. 10). Reservoir C was particularly sensitive to 

this numerical instability not only because of its steeply sloping stage-storage relation but also because the reservoir is in 

contact with a large aquifer, which resulted in large flow rates 

between the reservoir and aquifer even for small changes in 

reservoir storage. It should be noted that, because this numerical instability is caused by the strong curvature in the stagestorage relation (fig. 10), the instability is not apparent when 

the relation is linear, that is, when the reservoir has a rectangular shape.

Additional changes were made to the modified FYE 

model to resolve the numerical instability, and these modifications resulted in firm yields for reservoirs C and D. To 

test that these modifications did not affect the firm yield, the 

FYE model was rerun for reservoirs A and B; there was no 

difference between the firm yields calculated before and after 

the additional changes. Another resolution to the numerical 

instability would be for the user to fit a lower-order polynomial that does not account for the steep stage-storage relation near the bottom of the reservoir; however, it should be 

noted that the firm yield is sensitive to adjustments made in 

the stage-storage relation. To demonstrate this sensitivity, 

a comparison was made between the firm yield determined 

when the reservoir shape was represented as an ideal, straightsided reservoir square and the firm yield determined when the 

detailed bathymetric data were used to fit a stage-storage relation. This comparison showed that the firm yield was always 

greater when the more simplified reservoir shape was used in 

the simulation, with a maximum difference in firm yield of 

approximately 4 percent. 

Sensitivity of Firm Yield to Other Reservoir and 

Aquifer Input Variables

The sensitivity of firm yield to the reservoir and aquifer 

variables a, L, T, and S was quantified for reservoirs A through 

D (fig. 6), which represent a range of typical reservoir-aquifer 

systems, to understand the variability in firm yield that can be 

expected if estimates of the reservoir and aquifer characteristics contain uncertainty. For example, values of T and S may 

be difficult to estimate if detailed aquifer tests and well-bore 

Figure 10. Shape of the stage-storage relation for two drinking-water reservoirs in Massachusetts. As the reservoir drains, 

the slope increases sharply at the deepest part of the reservoir where the former stream channel is likely located. Two 

reservoirs are shown: A, Bearhole Reservoir, and B, Fitchburg Reservoir. 

140

145

150

155

160

165

0 20 40 60 80 100

STORAGE, IN MILLION GALLONS

RESERVOIR ELEVATION, IN FEET

Probable location of former stream channel

Probable location of former flood plain Probable location of former stream channel

Probable location of former flood plain

986

988

990

992

994

996

998

1,000

1,002

1,004

1,006

1,008

0 200 400 600 800

STORAGE, IN MILLION GALLONS

RESERVOIR ELEVATION, IN FEET

A. B.

Effect of Uncertainty in the Input Variables to the Ground-Water-Flow Equations  21
Figure 11. Sensitivity of firm yield to reservoir and aquifer characteristics for four idealized reservoir-aquifer systems.

-50

-40

-30

-20

-10

0

10

20

30

40

50

-100 -50 0 50 100

PERCENT CHANGE IN AQUIFER

CHARACTERISTIC

PERCENT CHANGE IN FIRM YIELD

RESERVOIR A RESERVOIR B

RESERVOIR C RESERVOIR D

Transmissivity

Distance from reservoir

shoreline to aquifer

boundary

Length of reservoir

shoreline in contact 

with aquifer

Storage coefficient

EXPLANATION

-50

-40

-30

-20

-10

0

10

20

30

40

50

-100 -50 0 50 100

PERCENT CHANGE IN AQUIFER

CHARACTERISTIC

-50

PERCENT CHANGE IN FIRM YIELD

-40

-30

-20

-10

0

10

20

30

40

50

-100 -50 0 50 100

PERCENT CHANGE IN AQUIFER

CHARACTERISTIC

PERCENT CHANGE IN FIRM YIELD

-50

-40

-30

-20

-10

0

10

20

30

40

50

-100 -50 0 50 100

PERCENT CHANGE IN AQUIFER

CHARACTERISTIC

PERCENT CHANGE IN FIRM YIELD

22  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
information are not available for the reservoir-aquifer system 

of interest. A sensitivity analysis can provide some information about how the firm yield is affected by these reservoir and 

aquifer characteristics.

For each study reservoir, the firm yield was determined 

from the entire historical record of inflows to and outflows 

from the reservoir by increasing and decreasing the values 

of a, L, T, and S (shown in table 1) by 40 and 80 percent for 

the simulation period. The percent change in firm yield was 

compared to the percent change in each reservoir-aquifer 

characteristic (fig. 11). For reservoirs in contact with a 

relatively small aquifer area (reservoir A and B), the firm yield 

is equally sensitive to changes in a, L, and S; however, this 

sensitivity does not change the firm yield by more than ±10 

percent(fig. 11). For reservoirs in contact with a relatively 

large aquifer area (reservoirs C and D), the firm yield is 

sensitive to changes in L, S, and T but not to changes in a

(fig. 11). The firm yield is most sensitive to changes in L, 

and the sensitivities of firm yield to changes in S and T are 

approximately equivalent. Reservoir C encountered numerical 

instability when S, L, and T were increased beyond their initial 

values; therefore, the change in firm yield due to increases in 

S, L, and T could not be determined (fig. 11).

Potential Effect of Ground-Water 

Contributions on Firm Yields for 

Reservoirs in Massachusetts

The modified FYE model assumes that both surface 

water and ground water contribute to reservoir storage. 

Therefore, if a reservoir has a clearly defined drainage area 

with some amount of surface-water inflow, the modified 

FYE model can be applied to reservoir-aquifer systems 

in Massachusetts; however, it is important to consider the 

applicability of the modified FYE model to a particular 

reservoir-aquifer system. For the case in which a reservoir is 

in contact with a large aquifer and the reservoir has a steep 

stage-storage relation, the modified FYE model may not result 

in a firm yield due to numerical instability. Furthermore, 

because it is not possible to measure actual ground-water-flow 

rates to compare with the ground-water-flow rates simulated 

by the modified FYE model, the results of the FYE model 

have not been validated by observed reservoir conditions. 

Despite these limitations, the modified FYE model offers 

a method to determine the firm yield of a reservoir whose 

shoreline is in contact with an aquifer in the absence of 

a more detailed, time-consuming, and expensive, two- or 

three-dimensional numerical ground-water-flow model. 

To determine the potential effects of ground-water 

contributions on firm yields, the firm yields of nine 

study reservoirs (fig. 12) that have reservoir shorelines 

in contact with sand and gravel (Waldron and Archfield, 

2006) were determined with and without the ground-water 

contributions. The reservoir and aquifer characteristics 

(table 2) and reservoir stage-storage relation were 

estimated for these nine reservoirs and the firm yields 

were then determined with and without consideration for 

the ground-water contribution to reservoir storage.

The firm yields were determined by using the entire 

historical record to solve the monthly water balance. Input 

variables needed to determine the ground-water contribution 

were estimated for each reservoir (table 2) in the following 

ways: L was estimated from GIS digital data layers of 

hydrography in Massachusetts at a 1:100,000 scale (Office 

of Geographic and Environmental Information (MassGIS), 

2004b); a was estimated from GIS digital data layers of 

surficial geology in Massachusetts at a 1:250,000 scale 

(Office of Geographic and Environmental Information 

(MassGIS), 2004a) and represents an average distance 

from the reservoir shoreline to the aquifer boundary; T

was estimated from hydrologic atlases or previous studies 

where a ground-water-flow model was developed (see table 

2 for references); a value of 0.2 was used for S; and the 

stage-storage relation was determined from the detailed 

bathymetric data from Waldron and Archfield (2006).

The often-used assumption in firm-yield estimations 

of negligible ground-water contributions for reservoirs in 

contact with aquifer sediments was discovered to be not a 

valid assumption for the reservoirs studied in Massachusetts. The firm yields for the study reservoirs consistently 

increased when the contribution of ground water determined 

by time-varying changes in reservoir stage was included 

in the reservoir water balance, with increases in firm yield 

ranging from 2 percent to 113 percent and a median increase 

in firm yield of 10 percent (table 3). Two reservoirs had an 

increase in firm yield greater than 85 percent. Therefore, 

potential changes in the firm yield due to the contribution of ground-water to reservoir storage are substantial; 

however, the amount by which the firm yield of reservoir 

increased was not found to be directly related to any one or 

combination of reservoir and aquifer characteristics. There 

may be cases for which the ground-water contribution to 

a reservoir whose shoreline is in contact with a sand and 

gravel aquifer is negligible; however, the firm yields of the 

study reservoirs were always affected by the inclusion of 

the ground-water contribution in the reservoir water balance 

(table 3). Therefore, without further study, it is not possible 

to establish a criterion as to when the ground-water contribution could be considered negligible for some reservoirs. 

Potential Effect of Ground-Water Contributions on Firm Yields for Reservoirs in Massachusetts  23
Figure 12. Location of nine reservoirs in Massachusetts that have shorelines in contact with sand and gravel aquifers.

RHODE ISLAND

CONNECTICUT

VERMONT

NEW HAMPSHIRE

NEW YORK

CAPE COD 

BAY

BOSTON BAY

BUZZARDS BAY

ATLANTIC

NANTUCKET SOUND

V ni eyard sound

OCEAN

EXPLANATION

TOWN BOUNDARY

Base from USGS and MassGIS

digital coverages, State Plane Coordinate System, NAD 83 

RESERVOIR AND IDENTIFIER

Danvers Water Department 1. Middleton PondFitchburg Water Department 2. Fitchburg Reservoir 3. Wachusett ReservoirLeominster Department of PublicWorks—Water Division 4. Distributing Reservoir 5. Goodfellow Reservoir 6. Morse ReservoirMarlborough Department ofPublic Works—Water and Sewer Division 7. Middleton Pond

Westfield Water Department 8. Granville ReservoirWest Springfield Water Department 9. Bearhole Reservoir

8

9

7

4

5

3 6

1

2

5

73o00'

42o30'

42o00'

41o30'

72o00'

71o00'

70o00'

50 MILES 0 10 20 30 40

10 20 30 40 50 KILOMETERS 0

24  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
Table 2. Estimated reservoir and aquifer characteristics for nine reservoirs in Massachusetts whose reservoir shorelines are in contact with a sand and gravel aquifer. 

[T, transmissivity; S, storage coefficient; a, distance from reservoir shoreline to aquifer boundary; L, length of reservoir shoreline in contact with the aquifer; DPW, Department of Public Works]

Water supplier Reservoir name L, in feet a, in feet T, in feet squared per day S, dimensionless

Danvers Water Department Middleton Pond 1,951.5 626.9 11,350 0.2Fitchburg Water Department Fitchburg Reservoir 2,900.0 271.6 11,350 .2

Wachusett Reservoir 10,918.6 904.0 12,000 .2

Leominster DPW–Water Division Distributing Reservoir 1,568.7 966.4 11,350 .2

Goodfellow Pond 121.0 1,233.5 11,350 .2Morse Reservoir 1,950.0 269.4 11,350 .2

Marlborough DPW, Water and Sewer Division Millham Reservoir 2,800.0 1,862.0 21,800 .18Westfield Water Department Granville Reservoir 9,219.5 463.4 3750 .2West Springfield Water Department Bearhole Reservoir 6,450.0 2,581.5 1,400 .2

1. Brackley and Hanson, 19772.DeSimone, 20043.Maevsky and Johnson, 1990

Table 3. Potential effect of ground-water contributions to reservoir storage on firm yield for nine reservoirs in Massachusetts.[Reservoir ID, identifier from figure 12; Mgal/d, million gallons per day; DPW, Department of Public Works]

Water supplier Reservoir name Reservoir ID

Reservoir firm yield 

without the groundwater contribution, in 

Mgal/d

Reservoir firm yield 

with the ground-water 

contribution, in Mgal/d

Percent increase 

in firm yield due 

to ground-water 

contributions

Danvers Water Department Middleton Pond 1 1.445 1.494 3Fitchburg Water Department Fitchburg Reservoir 2 1.341 1.367 2

Wachusett Reservoir 3 1.097 1.260 15

Leominster DPW

–Water Division Distributing Reservoir 4 .131 .279 113

Goodfellow Pond 5 .192 .199 4Morse Reservoir 6 .238 .259 9

Marlborough DPW, Water and Sewer Division Millham Reservoir 7 1.399 1.537 10Westfield Water Department Granville Reservoir 8 3.380 3.949 17West Springfield Water Department Bearhole Reservoir 9 1.360 2.525 86

Potential Effect of Ground-Water Contributions on Firm Yields for Reservoirs in Massachusetts  25
Summary

The amount of water that is available for withdrawal 

from a drinking-water supply reservoir depends on the water 

balance of the reservoir. Because many reservoirs receive 

most of their water from a surface-water source, the groundwater-inflow term in the water balance is commonly assumed 

to be negligible and is ignored in the calculation of the firm 

yield—the maximum yield that can be continually withdrawn 

from the reservoir. For a reservoir whose shoreline is in contact with an aquifer, however, as is the case for approximately 

one-quarter of the reservoirs in Massachusetts, the effect of 

ground-water contributions on firm yield can be substantial, 

but is unknown because available models that determine firm 

yield were constructed without consideration for the groundwater contribution to reservoir storage.

To evaluate the potential effect of ground water on firm 

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

Massachusetts Department of Environmental Protection, 

identified a set of equations that are based on an analytical 

solution to the ground-water-flow equation for the case of 

one-dimensional ground-water flow in a finite-width aquifer 

bounded by a linear surface-water feature such as a stream. 

These equations, which require only five input variables, were 

incorporated into the existing firm-yield-estimator (FYE) 

model, and the potential effect of ground water on firm yield 

was evaluated. 

To apply the FYE model to a reservoir-aquifer system 

in Massachusetts, the model requires that reservoir drainage 

areas be clearly defined and that some surface water flows into 

the reservoir. To determine the ground-water contribution to 

reservoir storage, the modified FYE model also requires the 

estimation of the reservoir storage-reservoir stage relation, the 

transmissivity (T) and storage coefficient (S) of the aquifer, an 

estimate of the length of reservoir shoreline that is in contact 

with the aquifer (L), and an estimate of the distance between 

the reservoir shoreline to the aquifer boundary (a). For surface-water-body shapes having a more realistic representation 

of a reservoir shoreline than a stream, a comparison of groundwater-flow rates simulated by the ground-water equations 

and ground-water-flow rates simulated by a two-dimensional, 

finite-difference ground-water-flow model indicate that the 

agreement between the simulated flow rates is within ±10 percent when the ratio of the distance from the reservoir shoreline 

to the aquifer boundary to the length of shoreline in contact 

with the aquifer is between values of 0.5 and 3.5.

The characteristic shape of the storage-stage relation 

for reservoirs in Massachusetts created numerical instability for some reservoirs when the water level in the reservoir 

approached the reservoir bottom. Steep slopes in the reservoir 

stage-storage relation resulted in ground-water-flow rates that 

were unrealistically high and prevented the calculation of firm 

yield. Reservoirs in contact with a large aquifer area are particularly susceptible to this problem. Fitting a monotonically 

increasing, parsimonious function to represent the reservoir 

stage-storage relation will minimize the possibly of numerical 

instability. Furthermore, differences in firm yield that result 

from oversimplification of the storage-stage relation are within 

5 percent.

For reservoirs in contact with a small aquifer area, the 

firm yield is equally sensitive to changes in a, L, and S; however, this sensitivity does not change the firm yield by more 

than ±10 percent. Therefore, uncertainty in the estimates of 

input variables a, L, T, and S for relatively small aquifer areas 

are unlikely to substantially affect the resulting firm-yield 

estimates. For reservoirs in contact with a large aquifer area, 

the firm yield is sensitive to changes in L, S, and T but not to 

changes in a. The firm yield is most sensitive to changes in 

L, and the sensitivities of firm yield to changes in S and T are 

approximately equivalent. 

The effect of ground-water contributions on firm yield 

was not only substantial, but furthermore, the firm yield of a 

reservoir in contact with a sand and gravel aquifer was always 

greater when the ground-water contribution to reservoir storage was included in the water balance. Increases in firm yield 

ranged from 2 to 113 percent, with a median increase in firm 

yield of 10 percent. Two reservoirs had an increase in firm 

yield that was greater than 85 percent. 

Idealized reservoir-aquifer systems also were assumed to 

verify that the ground-water-flow equations were implemented 

correctly into the existing FYE model; however, the modified 

FYE model has not been validated through a comparison of 

simulated and observed data. A comparison of simulated and 

observed reservoir data would further define limitations to the 

applicability of the ground-water-flow equations to reservoirs 

in Massachusetts whose shorelines are in contact with a sand 

and gravel aquifer.

Acknowledgments

The authors thank Paul Blain and the Massachusetts 

Department of Environmental Protection for their support of 

this project; USGS colleagues Stephen Garabedian, Paul Barlow, Marcus Waldron, and Forest Lyford for their invaluable 

technical assistance; and Philip Harte and John Masterson of 

the USGS and Thomas Lamonte of the Massachusetts Department of Environmental Protection for their report reviews.

26  Ground-Water Contributions to Reservoir Storage and the Effect on Estimates of Firm Yield
References Cited

Barlow, P.M., and Moench, A.F., 1998, Analytical solutions 

and computer programs for hydraulic interaction of streamaquifer systems: U.S. Geological Survey Open-File Report 

98-415A, 85 p.

Brackley, R.A., and Hansen, B.P., 1977, Water resources of the 

Nashua and Souhegan River Basins, Massachusetts: U.S. 

Geological Survey Hydrologic Investigations Atlas HA-276, 

2 sheets.

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

evaluation of water-management alternatives in the Assabet 

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

Scientific Investigations Report 2004-5114, 133 p.

Fennessey, N.M., 1994, A hydro-climatological model of daily 

streamflow for the northeast United States: Medford, Mass., 

Tufts University, Ph.D. dissertation, variously paged.

Fennessey, N.M., and Vogel, R.M., 1996, Regional models of 

potential evaporation and reference evapotranspiration for 

the Northeast USA: Journal of Hydrology, v. 184, p. 337-

354.

Fetter, C.W., 1994, Applied hydrogeology (3d ed.): Upper 

Saddle River, N.J., Prentice-Hall, 691 p.

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.

Maevsky, Anthony, and Johnson, D.G., 1990, Water resources 

of the Westfield and Farmington River Basins, Massachusetts: U.S. Geological Survey Hydrologic Investigations 

Atlas HA-716, 4 sheets.

Massachusetts Department of Environmental Protection, 1996, 

Estimating the firm yield of a surface water reservoir supply 

system in Massachusetts—A guidance document, Version 

1.0: Massachusetts Department of Environmental Protection, Office of Watershed Management, variously paged.

Office of Geographic and Environmental Information (MassGIS), Version 1.0, Software documentation: Prepared 

by Cambridge Environmental Inc. for the Massachusetts 

Department of Environmental Protection, Office of Watershed Management, variously paged.

Natural Resources Conservation Service, 1986, Urban hydrology for small watersheds: Natural Resources Conservation 

Service, Conservation Engineering Division, Technical 

Release 55, 164 p. (http://www.hydrocad.net/tr-55.htm)

Office of Geographic and Environmental Information (MassGIS), 2004a, Commonwealth of Massachusetts Executive 

Office of Environmental Affairs, 1:250,000-scale surficial 

geology digital data layer for sand and gravel deposits in 

Massachusetts; accessed on July 5, 2005, at http://www.

mass.gov/mgis/ftpsg.htm

Massachusetts Department of Environmental Protection, 

2004b, Commonwealth of Massachusetts Executive Office 

of Environmental Affairs, 1:100,000-scale hydrography 

digital data layer for surface-water bodies in Massachusetts; 

accessed on July 5, 2005, at http://www.mass.gov/mgis/

ftphd100_.htm

Rorabaugh, M.I., 1964, Estimating changes in bank storage 

and ground-water contribution to streamflow: International 

Association of Scientific Hydrology, Publication 63, p. 432441.

References Cited  27
For additional information write to:

Director, 

USGS Massachusetts–Rhode Island Water Science Center

10 Bearfoot Road

Northborough, MA 01532

or visit our Web site at

http://ma.water.usgs.gov