Hydrologic Models—HSPF and MODFLOW

The precipitation-runoff model HSPF and the groundwater-flow model MODFLOW were used to evaluate the effects 

of water-management strategies in the Pawtucket River Basin 

as described in the previous chapters. HSPF and MODFLOW 

are designed differently and serve different purposes. Each 

model simulates hydrologic processes through the calculation 

of a water budget, but the components of the simulated water 

budget and how the budgets are calculated differ greatly. 

In HSPF, water budgets are simulated on the basis of 

hydrologic response units (HRUs) that characterize the water 

fluxes and storage in pervious and impervious areas from 

direct inputs of meteorological time-series data. HRUs provide 

a distributed calculation of surface runoff, interflow, and 

groundwater flow to streams by processes that determine the 

fate of water through various losses and storage components. 

Flows from HRUs are typically directed to streams and routed 

by the kinematic-wave method to simulate streamflow. Userspecified variable values control HRU processes and streamflow routing. HSPF capabilities make it an appropriate tool for 

the continuous simulation of total streamflow. 

In MODFLOW, water budgets are simulated on the basis 

of specified boundary conditions, stresses, and hydrogeologic 

properties of the aquifer by numerical finite-difference methods to determine groundwater heads and fluxes (baseflow). 

Boundary conditions, stresses, or both can remain constant 

(steady state) or vary in time (transient) through user-specified 

values at the beginning of each time increment (referred to as 

a stress period). MODFLOW capabilities make it an appropriate tool for simulation of groundwater elevation and baseflow 

discharge to streams. 

The appropriate model for simulating hydrologic conditions and posing questions about how those conditions vary 

in response to change depends largely on the nature of the 

problem and the water-management decisions to be made. The 

purpose of this section is to compare HSPF and MODFLOW 

simulation results for specific scenarios to provide a better 

understanding of capabilities and limitations of both models. 

This in part relies on understanding some of the functional 

differences between the models and how these differences can 

affect simulations. 

Functional Differences between HSPF and 

MODFLOW

HSPF and MODFLOW represent hydrologic processes 

in different ways that are described in detail in the respective 

model documentation. Two pragmatic differences between the 

models are the temporal and spatial representations of hydrologic conditions or features, often referred to as the model 

discretization, which can have a large effect on the simulation 

of hydrologic processes. Temporal discretization affects the 

representation of time-varying hydrologic conditions, whereas 

spatial discretization affects the representation of the hydrologic characteristics of the basin or aquifer. Spatial discretization is determined largely by the extent of the model domain, 

and by the underlying hydrologic processes represented in 

each model.

HSPF simulations are typically made with an hourly time 

step (as in the Pawcatuck River Basin model) driven by timeseries data with the same time step. MODFLOW simulates 

steady-state or transient conditions that reflect long-term 

average or time-varying conditions, respectively. Typically, 

transient simulations use a monthly stress period to reflect 

seasonal fluctuations (like most stress periods specified in 

this study); however, shorter stress periods can be specified 

(such as the weekly stress period specified for the 2002 summer period). For each stress period, a unique set of stressors, 

boundary conditions, or both are specified; typically, this 

entails changes in recharge (boundary condition), groundwater 

withdrawals (stressor), or both, but changes are not limited to 

these properties. The averaging of stressors such as withdrawals over time can affect simulation results. 

Spatially, HSPF is considered by some to be a lumpedparameter model because the HRUs that define the model 

are amorphic—that is, they are not explicitly represented 

spatially. Rather, they are represented by the total area of each 

HRU that contributes to a reach. Others consider HSPF to 

be a spatially distributed model because each HRU uniquely 

represents hydrologic processes that are appropriately distributed throughout the basin. The combined areas of all HRUs 

to a reach determine the overall hydrologic response at a 

specified point (node) in the model. Hydrologic fluxes to any 

point in the model can be computed by the summation of HRU 

water budgets that compose the contributing area to that point. 

Part 5. HSPF and MODFLOW—Capabilities, Limitations, and 

Integration

By Phillip J. Zarriello
108 Part 5. HSPF and MODFLOW—Capabilities, Limitations, and Integration

It should be noted that there is little basis for evaluating the 

accuracy of computed flows from a specific HRU because the 

model-calibration points (streamflow-gaging stations) generally receive water from the aggregate of many different HRUs.

In MODFLOW, spatial discretization is governed by the 

size of the model grid and the necessity to keep the total number of cells in the grid within manageable limits in terms of 

computational speed and output size. The model-grid spacing 

determines the level of detail in the spatial representation of 

hydrologic properties of the aquifer. The hydrologic properties of each cell are typically regionalized and calibrated to 

point information distributed over the model domain, such as 

well logs, observation well heads, and streamflow at gaging 

stations. The calibration processes for MODFLOW and HSPF 

are similar in that individual components of the models are 

not explicitly calibrated; rather, the models are calibrated on 

the basis of the composite response of multiple components 

to point information such as streamflows at gaging stations or 

water levels at observation wells. 

In addition to the temporal and spatial limitations 

imposed by a specific model, a key factor to consider in the 

evaluation of simulation results is the hydrologic processes 

that are simulated. Alley and others (2002) articulated this 

consideration for groundwater models; however, the statement 

applies equally to other deterministic models: 

“Accuracy of model predictions is constrained by 

the correctness of the model (i.e., proper representation of relevant processes) and uncertainty in model 

parameters. The latter uncertainty is due to the 

limited accuracy with which parameter values can 

be measured and, more important, to the substantial 

heterogeneity inherent in aquifer characteristics. The 

inability to describe and represent this heterogeneity 

adequately is a fundamental problem in groundwater 

hydrology and will continue, even with improved 

models, to place limits on the reliability of model 

predictions. The links between spatial heterogeneity and model uncertainty also depend on the type 

of questions being asked. For example, reasonable 

estimation of head distributions in an aquifer may 

require only limited understanding of spatial heterogeneity. On the other hand, confidence in predictions 

of chemical concentrations at a specific location can 

be very sensitive to minor uncertainty in the spatial 

distribution of hydraulic properties, even for relatively homogeneous porous media.”

Although there are uncertainties associated with a given 

hydrologic process simulated by any model, fundamental 

differences in the process being considered can make the 

choice of the most appropriate model simple. For example, 

when the issue pertains to groundwater levels, MODFLOW 

would be the appropriate choice. When the issue pertains to 

streamflow, the choice of model can be more difficult because 

of inherent differences in the flow components simulated 

and the relation of the flow stressors to surface-water and 

groundwater interactions. These interactions are difficult to 

observe and measure, and are further complicated by many 

natural processes and human activities that commonly have 

been ignored in water-management considerations and policies 

(Winter and others, 1998). 

Comparison of Three Example HSPF and 

MODFLOW Results

In the Pawcatuck River Basin, HSPF and MODFLOW 

often produced different results for the same type of scenario, 

underscoring some of the fundamental differences between 

these models. Three examples were chosen to illustrate how 

(1) temporal discretization, herein referred to as time-step 

averaging, (2) spatial discretization, and (3) streamflowdepletion calculations can affect model results.

Example 1. Effects of a Pumped Well on 

Streamflow in Meadow Brook

The effects on streamflow of a hypothetical well (RIW–

481A) pumped at a constant rate of 1.55 ft3/s (1 Mgal/d) were 

simulated in lower Meadow Brook (MEAD2; fig. 5–1 and 

table A4–2). In MODFLOW, the pumped well was simulated 

at a distance of 200 ft from Meadow Brook under steady-state 

and transient conditions. Steady-state simulations represent 

average conditions from January 1, 2000, to September 30, 

2004; transient simulations represent monthly conditions 

over the same period, except for May through September 

2002, during which average weekly conditions are represented. Baseflow at the upstream boundary of MODFLOW in 

Meadow Brook was specified as the baseflow component of 

discharge calculated by HSPF at MEAD1 (fig. 5–1). Recharge 

values specified in MODFLOW for steady-state and transient 

simulations were derived from the flow into active-groundwater storage in HSPF as described in the recharge section in 

Appendix Part 3. 

MODFLOW simulated steady-state flow at the outlet of 

MEAD2 without RIW–481A pumping was about 6.8 ft3/s; 

with pumping, flow in the reach decreased by 1.36 ft3/s, or 

about 0.19 ft3/s less than the pumping rate of the well. The 

difference between streamflow loss and the withdrawal rate 

indicates that under steady-state conditions, simulated water 

to the well is captured from outside the surface-water subbasin 

divide; the captured water would otherwise likely discharge 

to the Pawcatuck River given the proximity of the well to the 

PAWC4 subbasin (fig. 5–1). 

HSPF simulated flow at the outlet of MEAD2 over the 

same period averaged about 14.1 ft3/s without RIW–481A 

pumping and about 12.5 ft3/s with the well pumping. The 

difference in average streamflow with and without pumping 

is equal to the withdrawal rate because, in the HSPF simulation, all water to the well is satisfied from flow that would 

otherwise discharge to Meadow Brook. Although the contribution of water to a pumped well from other reaches could be 

assigned in the HSPF model (in this example, 12 percent of 
Hydrologic Models—HSPF and MODFLOW 109

Figure 5–1. Hydrologic Simulation Program-FORTRAN (HSPF) and the modular groundwater-flow model (MODFLOW) 

representation of Meadow Brook and lower Wood River area, Pawcatuck River Basin, southwestern Rhode Island. 

(Surface-water sites described in table A2–4.)

RIW-481A

RIW-550

PAWC3

TOMA1

CANO1

WOOD5

PAWC4

CEDA1

MEAD1

MEAD2

COON1

WOOD6

PAWC5

WOOD4

0 1 2 MILES

0 1 2 KILOMETERS MEADOW

BROOK

WOOD

RIVER

PAWCATUCK

RIVER

I

*

*

0

Area of basin 

enlarged above

MEAD2

EXPLANATION

HSPF subbasin boundary

and reach name

Non-forested wetlands

Forested wetlands

Active groundwater-flow-model 

 area, cell size 100 by 100 feet

Groundwater-flow model grid

Simulated well and identifier

Streamflow-gaging station and identifier

RIW-481A

Partial-record streamflow-gaging station 

 and identifier

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

Rhode Island state plane projection, NAD83

Land-use data from RGIS, 1995

RIW-500

RIW 576

RIW-567

RIW-442

RIW-423

RIW-434

RIW-510

RIW-458

Diamond

Bog DIAMOND

BROOK

01118000

01118010

01117600

01117500

01117700

01118009

01118008

01118006

0111810

01117700
110 Part 5. HSPF and MODFLOW—Capabilities, Limitations, and Integration

the water to RIW–481A is likely obtained from PAWC4 under 

steady-state conditions), the amount of withdrawal that should 

be assigned to other reaches is not known without the information gained through MODFLOW simulations. 

MODFLOW simulations demonstrate that the 

contributing area to the pumped well is not bound by the 

surface-water divide. Under average 2000–04 conditions, 

HSPF simulated flow in Meadow Brook was about twice that 

of MODFLOW steady-state simulated flows, but because 

average streamflow was greater than the withdrawal rate 

of RIW–481A, the relative differences between models for 

simulated flows with and without pumping were minor. 

MODFLOW transient simulations indicate that the relative 

difference in flows simulated with and without pumping 

differs considerably between models at low flows. Transient 

MODFLOW simulations from January 1, 2000 through 

September 30, 2004 indicate that flow in MEAD2 ranged 

between zero and about 32 ft3

/s (fig. 5–2). During the normally 

wet spring months when the simulated streamflow was greater 

than the withdrawal rate, streamflow decreased in proportion 

to withdrawals from RIW–481A (1.55 ft3/s). During lowflow periods, particularly the summers of 2000, 2001, and 

2002, MODFLOW simulated streamflow stopped flowing 

at MEAD2 even without RIW–481A pumping. As a result, 

MODFLOW simulations indicate that withdrawals from RIW–

481A had little or no effect on streamflow because Meadow 

Brook had little or no simulated flow prior to pumping. Thus, 

MODFLOW-simulated changes in streamflow from pumping 

were limited to the periods when the brook had flow prior to 

the additional withdrawal; generally this resulted in extended 

periods of no flow because the simulated low flows prior 

to RIW–481A pumping were near or below the withdrawal 

rate during low-flow periods. This is particularly evident in 

simulated 2001 summer flows, when withdrawals composed 

an appreciable fraction of the simulated flows (fig. 5–2).

Daily HSPF simulations for the same period indicate that 

flow at MEAD2 ranged between about 2.0 and 100 ft3/s without RIW–481A pumping and between 0.3 and 100 ft3

/s with 

RIW–481A pumping (fig. 5–2). HSPF-simulated streamflow at 

the upstream streamflow-gaging station at MEAD 1 (fig. 5–1) 

is generally in good agreement with observed and estimated 

flows (Appendix Part 2; fig. A2–15) indicating that simulated flows in MEAD2 are better represented by HSPF than 

MODFLOW. The resulting interaction of pumping RIW–481A 

on streamflow, therefore, cannot be clearly defined by MODFLOW simulations, and the resulting decline in the water table 

around RIW–481A is likely oversimulated. 

Another consideration that could affect the interpretation of HSPF model results is changes in groundwater storage 

that are not directly simulated. Observed and HSPF-simulated 

flows at MEAD1 often fell below the 1.55 ft3/s pumping 

rate specified for RIW–481A. If flows in the reach near the 

0.1

1

10

100

1,000

DISCHARGE, IN CUBIC FEET PER SECOND

F A J A O D

2000

F A J A O D

2001

F A J A O D

2002

F A J A O D

2003

F A J A

2004

YEARS

Current

Pumped

Current

Pumped

MODFLOW

Daily mean Stress period

HSPF

EXPLANATION

Figure 5–2. Streamflow from January 1, 2000, through September 30, 2004, under 2000–04 (current) conditions and under current 

conditions with a well pumped at constant rate of 1 million gallons per day (pumped) simulated by the precipitation-runoff model 

Hydrologic Simulation Program-FORTRAN (HSPF) and the modular groundwater-flow model (MODFLOW) in lower Meadow Brook 

(MEAD2, RCHRES 48), Pawcatuck River Basin, southwestern Rhode Island. (Site location shown in figure 5–1 and described in 

table A2–4.) Months on x-axis are February, April, June, August, October, and December for years shown.
Hydrologic Models—HSPF and MODFLOW 111

pumped well were less than the withdrawal rate, withdrawals in excess of the streamflow would be neglected without 

further modifications to the HSPF model to compensate for the 

loss in groundwater storage. These modifications can be made 

through the Special Action feature in HSPF that accounts 

for withdrawals in excess of streamflow to mimic depleted 

groundwater storage, which is then replenished when streamflow exceeds withdrawals (Zarriello and Ries, 2000). In this 

example, flows at MEAD2 did not fall below the withdrawal 

because of the additional contributing area between MEAD1 

and MEAD2. Had streamflow near the pumped well been less 

than the withdrawal rate, the effect of lost storage on subsequent streamflows would not be realized in HSPF simulations 

without adding a Special Action similar to the one developed 

by Zarriello and Ries (2000). 

In MODFLOW streamflow is tracked cell by cell; therefore, interaction between the stream and aquifer is closely tied 

to the model cells closest to the pumped well. Conceivably, 

this spatial accountability provides a better estimate of flow in 

the vicinity of the pumped well where the interaction between 

surface water and groundwater is most dependent on flow, 

but its success is contingent on the accuracy of MODFLOWsimulated streamflow. In this example, the MODFLOW 

simulated flows are not as representative of actual flows as are 

the HSPF simulated flows, which can lead to other problems 

in the interpretation of model results as described. Consideration should also be given the fact that MODFLOW simulates 

only the baseflow component of streamflow; thus, actual flows 

are underrepresented at times that could result in a different 

aquifer response to a withdrawal. It should also be noted that 

a reasonably good connection between the stream and the 

aquifer is assumed; if the stream and the aquifer are not well 

connected, a different set of conclusions could be reached.

Example 2. Effects of Withdrawals near 

Diamond Bog

The simulation of surface-water and groundwater 

interactions can be affected by the spatial discretization of the 

model. Simulation of a hypothetical well, RIW–550, in the 

Wood River subbasin at WOOD5 (fig. 5–1) near Diamond 

Bog, an important ecological resource, illustrates the effects 

of the model spatial discretization on simulation results. 

RIW–550 was simulated in MODFLOW about 500 ft from 

the Wood River and about an equal distance to Diamond 

Brook with a constant pumping rate of 1.55 ft3/s (1 Mgal/d). 

In MODFLOW, Diamond Brook is explicitly represented, 

whereas in HSPF, Diamond Brook was lumped into the 

representation of the WOOD5 subbasin. In both HSPF and 

MODFLOW simulations, the effect on the Wood River of 

pumping RIW–550 at 1.55 ft3/s was essentially the same, other 

than the fractional differences in flow associated with the flow 

components simulated by each model. However, the effects of 

the pumped well on Diamond Bog could be greatly influenced 

by the spatial discretization of MODFLOW.

In MODFLOW simulations, streamflow in Diamond 

Brook was critical in the calculation of the response of the 

aquifer to pumping as described in Part 3. In summary, the 

water-level change in Diamond Bog in response to pumping was highly dependent on the simulated flow in Diamond 

Brook. Withdrawals from RIW–550 were often comparable to 

or greater than the flow in the brook. As the simulated flow in 

the brook becomes less than the withdrawal rate, the drawdown from the pumping well increases and expands laterally, a 

response similar to that described for the transient simulations 

of a pumped well near Meadow Brook. Therefore, the accuracy of the water-level response in Diamond Bog to pumping 

is highly dependent on the accuracy of the simulated flow in 

Diamond Brook, which is unknown. Other factors can also 

affect the interaction between surface water and groundwater, 

such as the hydraulic connection between the brook and the 

aquifer and the specified elevation of the streambed; however, 

the accuracy of these model variable values is not well known 

for this area, a problem common in most groundwater models 

because of the difficulty in obtaining representative values 

over a wide area. 

In general, variable values are empirically estimated 

during the model calibration. The accuracy of these variables 

can be further compromised in MODFLOW as the variable 

values represent larger cell areas. If the cell sizes increase 

to where less than one cell separates Diamond Brook 

and Diamond Bog from Wood River, the resulting model 

would have insufficient detail (spatial discretization) to 

differentiate the effects of pumping on these features. Had 

a single groundwater-flow model been developed for the 

entire Pawcatuck River Basin, the spatial detail would likely 

have been inadequate to represent the influence of Diamond 

Brook on pumping near the bog. This is an example of why 

two separate groundwater-flow models were developed to 

allow for better spatial resolution of hydrologic properties 

in key areas of the basin. Simulations of surface-water and 

groundwater interactions are further complicated by temporal 

conditions that may provide additional sources of water to 

the pumped well, such as ponded water in the bog or adjacent 

wetlands, or localized time-varying recharge, and other factors 

not accounted for in the model. 

The HSPF model of the Pawcatuck River Basin was 

not discretized at a scale fine enough to directly evaluate the 

effects of pumping on Diamond Brook (the drainage area 

to Diamond Brook at its mouth is about 0.9 mi2 in a model 

constructed to simulate the entire 303-mi2 Pawcatuck River 

Basin), nor was the model capable of simulating groundwater 

levels. If the model had been developed at a scale sufficient to 

represent flow in Diamond Brook, then additional information 

would be required to determine the portion of the withdrawal 

that is captured groundwater discharge to, or induced infiltration from, the brook. Captured groundwater flow to or induced 

infiltration from Diamond Brook is complex, varies with time, 

and is best represented by a transient groundwater model. As 

previously noted, these simulations are highly dependent on 

the accuracy of simulated flow in Diamond Brook and on how 
112 Part 5. HSPF and MODFLOW—Capabilities, Limitations, and Integration

accurately the hydraulic connection between the brook and the 

aquifer in simulated. 

Example 3. Effects of Converting from SurfaceWater to Groundwater Withdrawals

Converting withdrawals from surface-water to groundwater sources can be beneficial to aquatic habitat because of the 

time lag and damping of peak demand on streamflow induced 

by withdrawing water from the aquifer instead of the stream. 

Simulating the change in streamflow caused by converting 

withdrawals from surface-water to groundwater sources can 

be greatly influenced by averaging irregular withdrawals over 

time and by the assumptions inherent in STRMDEPL, which 

is used in HSPF to compute the lag effects of a pumped well 

on streamflow depletion. This is illustrated by examining the 

effects of moving withdrawals from surface-water to 

groundwater sources in the lower Beaver River (BEAV3). 

In 2000–04, irrigation withdrawals were made directly 

from streams at seven locations in BEAV3 (fig. 5–3). 

MODFLOW and HSPF simulations examined the effects of 

moving five of these withdrawals to a series of groundwater 

wells located about 2,000 to 3,000 ft from the river to an 

area about midway between BEAV3 and QUEN7 (the lower 

Usquepaug-Queen River). Simulations included moving 

a surface-water withdrawal (AQU10A) from QUEN7 to 

the same area as the withdrawals moved from BEAV3. 

MODFLOW simulations also investigated the effects of 

several withdrawal scenarios that consolidated groundwater 

withdrawals, but results described in Part 3 indicate no 

appreciable differences in streamflow depletion of BEAV3. 

Consequently, effects of consolidating individual pumped 

wells are not addressed in this section.

HSPF simulations also included converting two downstream withdrawals (AB3B and AB5A) to groundwater 

withdrawals in BEAV3 (fig. 5–3). Wells AB3B and AB5A 

supply water to about 18 percent of the total irrigated area 

in the lower Beaver River subbasin. The general effects of 

time-step averaging and limitations of STRMDEPL, however, 

are not appreciably affected by differences in the withdrawals 

simulated with and without wells AB3B and AB5A in HSPF 

and MODFLOW, respectively. 

The simulated effects of an irregularly pumped well 

on streamflow depend on the temporal discretization of the 

model. The temporal discretization dictates the extent to which 

the withdrawal variations are averaged or smoothed in the 

simulation time step (HSPF) or stress period (MODFLOW). 

Irrigation withdrawals are particularly affected by time-step 

averaging because they are typically made only during part of 

the day.

In HSPF, time-step averaging is mostly an issue with 

postprocessing of model results because simulations are made 

with an hourly time step; however, averaging withdrawals input to the model could affect the simulated streamflow 

response. For example, when irrigation withdrawals that turn 

on and off during the day are averaged, the resulting data 

would damp peak withdrawals and cause withdrawals during 

periods of no pumping. This effect becomes more pronounced 

as the time period over which withdrawals are averaged 

increases or withdrawal variability increases, and is most 

pronounced as they increase simultaneously. The effect of 

time-step averaging in the postprocessing of HSPF-simulated 

streamflow is illustrated in BEAV3 in figures 5–4A and B. 

The simulated hourly hydrograph (fig. 5–4A) during July 

through early September 2002 show an intradaily oscillation 

that reflects the surface-water withdrawal pattern; the daily 

mean hydrograph (fig. 5–4B) for the same time interval does 

not show this oscillation because of the smoothing effect of 

averaging hourly values over the day. 

The effects of time-step averaging of variable withdrawals, such as irrigation withdrawals, become less relevant the 

farther the withdrawal is moved from the stream because of 

the inherent damping of streamflow depletion by the aquifer. 

It should also be noted that the extent to which time-step 

averaging affects model results is directly related to the rate 

of withdrawal in relation to the rate of streamflow. In general, the importance of time-step averaging decreases as the 

withdrawals compose a smaller fraction of total streamflow 

and increases as the withdrawals compose a larger fraction 

of streamflow. This effect can be seen in the HSPF-simulated 

hourly streamflow under direct surface withdrawals 

(LT–Cdmd) during June 2002 when withdrawals composed 

only a small fraction of total streamflow relative to the later 

summer months when withdrawals composed a larger fraction 

of total streamflow (fig. 5–4A). 

In MODFLOW simulations, irrigation withdrawals 

were averaged over weekly or monthly stress periods. The 

smoothing effect of averaging irregular withdrawals over this 

period of time can diminish or mask the potential benefits 

of moving the withdrawal from a stream to a groundwater 

source, particularly when the withdrawal is large with respect 

to streamflow. During the 2002 irrigation season (June 

through September), estimated average monthly withdrawals 

in BEAV3 ranged from 0.46 to 2.01 ft3/s. During the same 

period, weekly, daily, and hourly withdrawals ranged from 

zero to a high of 2.05, 2.35, and 5.68 ft3/s, respectively 

(fig. 5–5). Averaging the hourly withdrawals over a day or 

over a month cut the peak hourly withdrawal rate by as much 

as 60 to 90 percent, respectively. In this particular example, 

the most pronounced effect of time-step averaging is when 

withdrawals occurred once during the week (June 11–18), 

yielding an average weekly rate about 94 and 86 percent less 

than the peak hourly and daily withdrawal rates, respectively. 

Note that these values include withdrawals from wells AB3B 

and AB5A, which were included in the HSPF simulations, 

but not in the MODFLOW simulations; exclusion of these 

withdrawals in the HSPF simulations would decrease the 

withdrawals from BEAV3 by about 20 percent from that 

shown. In MODFLOW, the simulated effect of direct 

withdrawals from the river was obtained by subtracting the 

average withdrawal from the simulated baseflow for a stress 

period. As a result, much of the short-term variation in the 
Hydrologic Models—HSPF and MODFLOW 113

outside of

model area

BEAVER

RIVER

USQUEPAUG

RIVER

AB5A

AB3B

AB2B

AB4A

AB2A

AB1A

PR-AB5A

PR-AB3B PR-AB2B

PR-AB4A

PR-AB3A

PR-AB2A

PR-AB1A

PAWC3

BEAV2

QUEN7

BEAV3

TANE1

QUEN6

PAWC2

AB3A

AQU10A

PR-AQU10A

BEAV3

EXPLANATION

HSPF subbasin boundary

and reach name

Non-forested wetlands

Forested wetlands

Groundwater-flow-model 

 grid cell 250 by 250 feet

Groundwater replacement well 

 and identifier

Streamflow-gaging station

 and identifier

0 1 2 MILES

0 1 2 KILOMETERS

AB2B Surface-water withdrawal and identifier

PR-AB2B

TANEY

BROOK

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

Rhode Island state plane projection, North American Datum of 1983

Land-use data from RGIS, 1995

71o38' 30’’

41o

30' 30’’

71o37' 

41o27' 30’’

I

*

*

0

Area of basin 

enlarged above

01117410

01117468

01117420

01117471

01117480

01117471

Partial-record streamflow-gaging 

 station and identifier

01117480

Figure 5–3. Hydrologic Simulation Program-FORTRAN (HSPF) and modular groundwater-flow model (MODFLOW) 

representation of the lower Beaver River area, Pawcatuck River Basin, southwestern Rhode Island. (Surface-water sites 

described in table A2–4 and groundwater sites described in table 4–1.)
114 Part 5. HSPF and MODFLOW—Capabilities, Limitations, and Integration

1

10

100

1,000

DISCHARGE, IN CUBIC FEET PER SECOND

F A J A O D

2000

F A J A O D

2001

F A J A O D

2002

F A J A O D

2003

F A J A

2004

SIMULATED STREAMFLOW 2000–04

Area 

enlarged 

above 

A

B

YEARS

LT-NoDmd (No withdrawals)

LT-Cdmd (Current withdrawals from surface-water)

LT-CDSWR (Current withdrawals from groundwater) 

0.1

10

100

10 20

MAY

31 10

JUNE

20 30 10

JULY

20 31 10

AUGUST

20 31 10

SEPTEMBER

20 30

1

DISCHARGE, IN CUBIC FEET PER SECOND

SIMULATED STREAMFLOW MAY–SEPTEMBER 2002

Surface-water 

withdrawals

Groundwater

withdrawals

MODFLOW HSPF

Weekly stress period Hourly time step

HSPF—Daily mean 

MODFLOW—Monthly stress period, except shaded area, which is weekly

Model and specified withdrawals scenario

Figure 5–4. Simulated streamflow from (A) May through September 2002 and (B) January 1, 2000, through September 30, 

2004, under current conditions with selected irrigation withdrawals converted from surface-water to groundwater sources 

with the precipitation-runoff model Hydrologic Simulation Program-FORTRAN (HSPF) and the modular groundwater-flow 

model (MODFLOW) in lower Beaver River (BEAV3, RCHRES 43, 01117471), Pawcatuck River Basin, southwestern Rhode Island.

(Location shown in figure 5–3). Months on x-axis are February, April, June, August, October, and December for years shown. 

(Site location shown in figure 5–3 and described in table A2–4.)
Hydrologic Models—HSPF and MODFLOW 115

streamflow response to direct withdrawals is lost due to timestep averaging. 

For example, in June 2002 withdrawals from BEAV3 

included distinct periods of pumping and no pumping 

(fig. 5–5). The monthly average withdrawal in June 2002 

specified in MODFLOW was 1.07 ft3/s, which results in 

undersimulation of streamflow by about 1 ft3/s during days of 

no pumping and oversimulation by about 0.9 ft3/s on days with 

pumping. These values represent between 5 and 10 percent 

of the MODFLOW-simulated streamflow at the time. The 

amount by which streamflows are oversimulated or undersimulated depends on the length of the stress period and the 

withdrawals during that period. Additionally, withdrawals 

from groundwater are not as prone to the effects of time-step 

averaging as direct stream withdrawals because of the lag time 

and damping effects of withdrawals on streamflow depletion 

by the aquifer. It should be noted that although time-step averaging can cause inaccuracies in simulated flows, the models 

themselves are still valid because the mass balance of water 

within the relevant time interval is unaffected by time-step 

averaging (that is, there is no net gain or loss of water).

Converting withdrawals from surface- to groundwater 

sources in the BEAV3 subbasin also illustrates the limitations 

of STRMDEPL used in HSPF to simulate the effects of a 

pumped well on streamflow. In the HSPF simulation 

LT–CDSWR, five of the seven surface-water withdrawals were moved from BEAV3 to groundwater wells in the 

contributing area of the PAWC2 subbasin (fig. 5–3). As such, 

the effects of withdrawals from these wells on streamflow 

was removed entirely from BEAV3 and shifted to PAWC2. 

MODFLOW simulations, on the other hand, indicate that on 

average, only about 2 percent of the water came from outside 

the BEAV3 subbasin when surface-water withdrawals were 

converted to groundwater withdrawals in PAWC2; average 

streamflow from January 2000 through September 2004 in 

BEAV3 increased from 18.3 to 18.6 ft3/s. It should be noted 

that boundary conditions specified in MODFLOW did not 

force water to the pumped wells to come from a particular surface-water subbasin; rather, the source of water to the pumped 

well is determined on the basis of calculated groundwater-flow 

direction, which is affected by the specified hydraulic conductance and the computed hydraulic gradient between the aquifer 

and the stream. The effects of converting from surface-water 

to groundwater withdrawals on streamflow in BEAV3 may 

be diminished by withdrawals from AQU10A (moved from 

QUEN7), which was also moved to the same general area in 

PAWC2 as the five withdrawals moved from BEAV3. Furthermore, MODFLOW simulations during dry periods indicate 

that the source water to the six pumped wells comes from storage or areas outside the basin, or both; near the end of August 

2002, streamflow in BEAV3 increased from about 2.9 to 

4.1 ft3/s when withdrawals were converted from surface-water 

to groundwater sources, respectively. Thus, MODFLOW 

simulations indicate that HSPF oversimulated streamflow in 

BEAV3 when direct stream withdrawals were converted to 

groundwater withdrawals outside the drainage divide. 

Without prior knowledge of the source of water to a 

pumped well, allocation of stream depletion among different 

reaches in HSPF cannot be made. This is problematic when 

a pumped well alters the groundwater flow path in ways that 

MAY JUNE JULY AUGUST SEPTEMBER

WITHDRAWALS, IN CUBIC FEET PER SECOND

Weekly average

Monthly average

Daily average

Hourly

0

1

2

3

4

5

6

10 20 31 10 20 30 10 20 31 10 20 31 10 20 30

2002

Figure 5–5. Effects of time-step averaging of estimated irrigation withdrawals, May through September 2002, in the lower 

Beaver River (BEAV3, RCHRES 43, 01117471), Pawcatuck River Basin, southwestern Rhode Island. (Site location shown in 

figure 5–3 and described in table A2–4.)
116 Part 5. HSPF and MODFLOW—Capabilities, Limitations, and Integration

are not easily determined, such as in areas near the water-table 

divide or in areas where the water-table gradient is low or 

poorly defined. These conditions are further complicated by 

temporal changes in boundary conditions and largely unknown 

parameter values that control the interaction between 

groundwater and surface water.

Simple allocations of streamflow depletions on the basis 

of subbasin divides can be inaccurate and lead to incorrect conclusions. Allocations made on the basis of distance 

from neighboring streams alone can also be misleading; in 

this example, wells are about equal distance from BEAV3 

and QUEN7, but the MODFLOW simulations indicate the 

majority of the water withdrawn is water that would otherwise 

discharge to BEAV3. Although numerical groundwater-flowmodel simulations are best suited for this type of analysis, 

constructing an accurate model requires detailed field investigations, particularly at the local scale, to ensure that simulated 

streambed conductance, elevation, and aquifer heterogeneity 

adequately represent actual geohydrologic conditions.

Integrating HSPF and MODFLOW Models

Coupled groundwater and surface-water models help 

address some of the complexities of groundwater and surfacewater interactions and the problems associated with time 

discretization and simulation of flow components described 

above. The Integrated Hydrologic Model (IHM) formally 

couples HSPF and MODFLOW models (Ross and others, 

2003; Aly and others, 2003) and allows each model to work 

independently at its own time step or stress period. The results 

of each model are fed back and forth to each other at their 

specified time step or stress period to provide a more complete representation of hydrologic processes than cannot be 

achieved by the individual models alone. Although this can 

greatly improve the simulation of hydrologic processes and, 

thus, the potential model uses, the formally coupled model is 

still subject to the same spatial-discretization limitations as the 

individual models. The spatial limitations, as noted previously, 

can influence results and interpretation. 

Compounding the model limitations are uncertainties associated with the input variable values required by 

each model. Issues of model uncertainty lead to questions of 

equifinality—that is, alternative model structures and variable 

values can only be rejected (Beven and Binley, 1992; Beven, 

1993). Alternative model structures and variable values can be 

accepted that could potentially lead to conflicting results. The 

question remains as to whether fully integrated surface-water 

and groundwater models limit or broaden the range of possible 

acceptable model alternatives. Integrated models may limit the 

range of acceptable variable values because of the improved 

representation of hydrologic processes; however, the added 

complexity of an integrated model, coupled with an increased 

number of input variables, could also increase variable uncertainty and the range of model alternatives. 

IHM links HSPF Hydrologic Response Units (HRUs) 

with MODFLOW cells to pass information between their 

operational units. HSPF passes recharge from the active 

groundwater component (AGWO) of the HRUs to the appropriate MODFOW cell. MODFLOW, in turn, passes hydraulic 

head information to HSPF for the simulation of lower-zone 

storage processes and discharge of groundwater to stream 

reaches. In addition, IHM better accounts for head-dependent 

evapotranspiration losses in the unsaturated (vadose) zone. 

For these reasons IHM was tested early in the present study to 

determine its potential use.

Testing of HSPF and MODFLOW Integration in 

the Usquepaug-Queen River

The IHM model was tested in the 36-mi2 UsquepaugQueen River subbasin (RCHRES 1-20) (fig. A2–4) by 

using modified versions of the previously developed HSPF 

(Zarriello and Bent, 2004) and MODFLOW (Dickerman and 

others, 1977) models. The MODFLOW model was rebuilt 

using a 200- by 200-ft cell spacing that extended to the basin 

boundaries. Aquifer properties were similar to those used by 

Dickerman and others (1977), except upland tills are simulated 

in the modified groundwater-flow model. The HSPF model 

was modified to better represent head-dependent processes 

between MODFLOW and HRUs. This entailed duplicating 

the original 15 pervious-area HRUs (PERLNDs) in the model 

across the 20 model subbasins, resulting in 300 unique HRUs 

that were mapped to MODFLOW cells. Mapping HRUs to 

MODFLOW cells was done through GIS and the creation of 

an ACCESS database table that identifies the proportional area 

of each HRU and other hydrologic characteristics, such as 

field capacity, associated with each MODFLOW cell. For the 

IHM test basin, the database that maps HRUs to MODFLOW 

cells was about 70,000 records. In addition to these records, 

ACCESS tables are required to map MODFLOW cells to 

HSPF reaches and to pass MODFLOW-simulated baseflow 

back to HSPF for flow routing with non-baseflow components 

of flow simulated by HSPF. 

Once the IHM interface files were complete, the integrated model ran HSPF with an hourly time step and MODFLOW ran with a daily stress period. Although the IHM model 

initially ran, the model would inexplicably stop during the 

early part of the simulation period. This may have been caused 

by the instability of the MODFLOW simulation of upland tills, 

which is characterized by constant wetting and drying from 

variable recharge, a thin unconfined upper layer, and high 

relief. The problem was not investigated further because a 

fully coupled HSPF-MODFLOW model for the entire 

Pawcatuck Basin, which is about 10 times the size of the 

Usquepaug-Queen test subbasin, was determined to be impracticable for several reasons. These reasons include limitations 

of a workable MODFLOW cell size, the number of additional 

HRUs needed to represent head-dependent conditions throughout the basin, and the size of the database needed to link the 

models. It was also determined that the project objectives 

could be met by informally linking HSPF and MODFLOW, 

whereby HSPF would generate recharge for MODFLOW as 

described in the development of the MODFLOW models in 

Appendix Part 3.