Documentation of the Streamflow-Routing (SFR2) 

Package to Include Unsaturated Flow Beneath 

Streams—A Modification to SFR1

Techniques and Methods 6-A13

Version 1.30, April 2010

A Product of the Ground-Water Resources Program

U.S. Department of the Interior

U.S. Geological Survey

Chapter 13 of

Section A, Ground Water, of

Book 6, Modeling Techniques

Fault

Stream

Intermittent stream

Intermittent stream

Bedrock

Streambed

Alluvium (unsaturated zone)

Alluvium (saturated zone)

Alluvium (saturated zone)

Alluvium (unsaturated zone)

Water table

Water table

A. Profile of an intermittent stream leaving the mountain front.

B. Cross section of stream separated from the water table by an unsaturated zone.
Documentation of the Streamflow-Routing 

(SFR2) Package to Include Unsaturated 

Flow Beneath Streams—A Modification 

to SFR1

By Richard G. Niswonger and David E. Prudic

Chapter 13 of

Section A, Ground Water, of

Book 6, Modeling Techniques

Techniques and Methods 6-A13

Version 1.0, APRIL 20

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: 2005

Revised: April 2010

For sale by U.S. Geological Survey, Information Services

Box 25286, Denver Federal Center

Denver, CO 80225

For more information about the USGS and its products:

Telephone: 1-888-ASK-USGS

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

Any use of trade, product, or firm names in this publication 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:

Niswonger, R.G., and Prudic, D.E., 2005, Documentation of the Streamflow-Routing (SFR2) Package to 

include unsaturated flow beneath streams—A modification to SFR1: U.S. Geological Survey Techniques and 

Methods 6-A13, 50 p.
iii

Preface

Performance of this computer program has been tested in models of hypothetical groundwater flow systems; however, future applications of the programs could reveal errors that 

were not detected in the test simulations. Users are requested to notify the U.S. Geological 

Survey (USGS) if errors are found in the documentation report or in the computer program. 

Correspondence regarding the report or program should be sent to:

U.S. Geological Survey, WRD

333 West Nye Lane, Room 203

Carson City, Nevada 89706 

Attention: Richard G. Niswonger

Although the computer program has been used by the USGS, no warranty, expressed or implied, 

is made by the USGS or the United States Government as to the accuracy and functionality of 

the program and related program material. Nor shall the fact of distribution constitute any such 

warranty, and no responsibility is assumed by the USGS in connection therewith. 

The revised Streamflow-Routing (SFR2) Package, MODFLOW-2000, and other ground-water 

programs are available for downloading from the USGS at the following World Wide Web 

(WWW) address: http://water.usgs.gov/software/ground_water.html.
iv

Contents

Preface ................................................................................................... iii

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

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

Purpose and Scope ................................................................................. 3

Acknowledgments ................................................................................. 3

Theory and Conceptualization of Unsaturated-Zone Flow .......................................... 3

Method ............................................................................................. 3

Applicability of Assumptions ..................................................................... 4

Surface-Water Flow and Infiltration ............................................................ 6

Coupling Unsaturated-Zone Flow in SFR2 to MODFLOW....................................... 7

Steady-State Simulations with SFR2 ............................................................ 8

Assumptions and Limitations ..................................................................... 8

Results and Discussion ................................................................................. 8

Model Comparisons .............................................................................. 8

Hypothetical Test Simulation ..................................................................... 12

Implementation of SFR2 with MODFLOW-2000 ...................................................... 15

Summary ................................................................................................... 15

References Cited.......................................................................................... 17

Appendix 1: Data Input Instructions for the Streamflow-Routing (SFR2) Package with 

Unsaturated Flow Beneath Streams ............................................................... 18

Appendix 2: Selected Input Data and Printed Results for the Hypothetical Test Simulation ... 34
v

Tables

Table 1. Hydraulic properties assigned in first test simulation assuming a column 

filled with sandy clay loam ............................................................... 9

Table 2. Input variables for second test simulation assuming vertical seepage through 

unsaturated zone that is restricted by a semipervious streambed .................. 10

Table 3. Comparison of calculated fluxes through a semipervious streambed 

resulting from four different stream stages from kinematic-wave 

approximation in the Streamflow-Routing (SFR2) Package with those from 

the Variably-Saturated Two-Dimensional Flow and Transport (VS2DT) Model ...... 10

Table 4. Unsaturated zone input variables used in example problem........................... 12

Figures

Figure 1. Conceptualization of a disconnected stream/aquifer system ........................ 2

Figure 2. Graph showing a wetting front through a uniform column of unsaturated 

material affected by a decrease in surface flux and results from a 

kinematic-wave approximation of the wetting front using leading and 

trailing waves .............................................................................. 5

Figure 3. Diagram showing discretization of unsaturated zone beneath a stream in a 

single MODFLOW cell ..................................................................... 7

Figure 4. Graphs showing comparison of results from the kinematic-wave 

approximation used in the Streamflow-Routing (SFR2) Package with the 

Variably Saturated Two-Dimensional Flow and Transport (VS2DT) Model ......... 9

Figure 5. Graph showing comparison of results from the kinematic-wave 

approximation in the Streamflow-Routing (SFR2) Package with the Variably 

Saturated Two-Dimensional Flow and Transport (VS2DT; Healy, 1990) Model 

for seepage through a one-dimensional column filled with a silt loam ............... 11

Figure 6. Diagram showing conceptualization of hypothetical test simulation ............... 13

Figure 7. Graphs showing effects of variations in thickness of a semipervious streambed... 14
vi

Conversion Factors

Multiply By To obtain

cubic meter (m3) 35.3147 cubic foot

kilometer (km) 0.621371 mile

meter (m) 3.28084 foot

meter per second (m/s) 3.281 foot per second
Abstract

Many streams in the United States, especially those 

in semiarid regions, have reaches that are hydraulically 

disconnected from underlying aquifers. Ground-water 

withdrawals have decreased water levels in valley aquifers 

beneath streams, increasing the occurrence of disconnected 

streams and aquifers. The U.S. Geological Survey modular 

ground-water model (MODFLOW-2000) can be used to model 

these interactions using the Streamflow-Routing (SFR1) 

Package. However, the approach does not consider unsaturated 

flow between streams and aquifers and may not give realistic 

results in areas with significantly deep unsaturated zones. 

This documentation describes a method for extending the 

capabilities of MODFLOW-2000 by incorporating the ability 

to simulate unsaturated flow beneath streams. 

A kinematic-wave approximation to Richards’ equation 

was solved by the method of characteristics to simulate 

unsaturated flow beneath streams in SFR1. This new 

package, called SFR2, includes all the capabilities of SFR1 

and is designed to be used with MODFLOW-2000. Unlike 

SFR1, seepage loss from the stream may be restricted by the 

hydraulic conductivity of the unsaturated zone. Unsaturated 

flow is simulated independently of saturated flow within each 

model cell corresponding to a stream reach whenever the 

water table (head in MODFLOW) is below the elevation of 

the streambed. The relation between unsaturated hydraulic 

conductivity and water content is defined by the BrooksCorey function. Unsaturated flow variables specified in SFR2 

include saturated and initial water contents; saturated vertical 

hydraulic conductivity; and the Brooks-Corey exponent. 

These variables are defined independently for each stream 

reach. Unsaturated flow in SFR2 was compared to the U.S. 

Geological Survey’s Variably Saturated Two-Dimensional 

Flow and Transport (VS2DT) Model for two test simulations. 

For both test simulations, results of the two models were in 

good agreement with respect to the magnitude and downward 

progression of a wetting front through an unsaturated 

column. A third hypothetical simulation is presented that 

includes interaction between a stream and aquifer separated 

by an unsaturated zone. This simulation is included to 

demonstrate the utility of unsaturated flow in SFR2 with 

MODFLOW-2000. This report includes a description of the 

data input requirements for simulating unsaturated flow in 

SFR2.

Introduction

For several decades, researchers have been modeling 

surface and subsurface flow as a continuum, for the purpose 

of considering all exchanges of water between the land 

surface and the underlying ground water. The U.S. Geological 

Survey (USGS) modular ground-water model, MODFLOW 

(McDonald and Harbaugh, 1988; and Harbaugh and others, 

2000), has been used to model surface- and ground-water 

interaction over a wide range of spatial scales. The application 

of MODFLOW to simulate the interaction of surface water 

with ground water ignores the time delay for water to flow 

from surface-water bodies at land surface to the water table. 

Where the depth to ground water is less than a few meters 

below the streambed, the time for water to flow through the 

unsaturated zone can be considered negligible. However, 

where the depth to ground water is tens to hundreds of meters 

below a streambed (fig. 1), the storage and flow of water 

Documentation of the Streamflow-Routing (SFR2) Package 

to Include Unsaturated Flow Beneath Streams—A 

Modification to SFR1

By Richard G. Niswonger and David E. Prudic
within the unsaturated zone may be important when assessing 

changes in the water resources of those areas. 

Because of the importance of water in arid regions, and, 

in particular, because of increased interest in the interaction of 

surface water and ground water, the capability of simulating 

unsaturated flow beneath streams has been added to the 

Streamflow-Routing (SFR1) Package (Prudic and others, 

2004) of MODFLOW-2000. This new package is called 

SFR2. Unsaturated flow beneath streams is simulated using a 

kinematic-wave approximation to Richards’ equation that is 

solved by the method of characteristics (Smith, 1983; Smith 

and Hebbert, 1983). The kinematic-wave approximation 

to Richards’ equation neglects diffusive forces and flow 

is assumed to occur in the vertical-downward direction 

through homogeneous and isotropic sediment. The method 

of characteristics is used to reduce the one-dimensional 

partial-differential equation resulting from the kinematicwave approximation to an ordinary differential equation that 

is solved by analytical integration. The solution tracks the 

movement of waves through the variably saturated zone and 

precludes the need to discretize the region between the water 

table and land surface. The Lagrangian formulation of the 

Figure 1. Conceptualization of a disconnected stream/aquifer system.

Fault

Intermittent stream

Streambed

Alluvium (saturated zone)

Alluvium (unsaturated zone)

Water table

Not to Scale 

Stream

Intermittent stream

Bedrock

Alluvium (unsaturated zone)

Alluvium (saturated zone)

Water table

Not to Scale 

A. Profile of an intermittent stream leaving the mountain front.

B. Cross section of stream separated from the water table by an unsaturated zone.

2 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
model allows seepage and unsaturated flow to occur over 

short time steps, while still maintaining suitably long time 

steps for ground-water flow in MODFLOW. This approach 

allows the simulation of flow and storage in the unsaturated 

zone, and for the simulation of a time delay for recharge while 

still maintaining the applicability of MODFLOW to basinscale problems. Assumptions made in this model may be an 

over-simplification for some natural systems where negative 

pressure gradients are important, such as in fine-grained or 

heterogeneous sediment.

There is another significant difference between 

simulating disconnected stream/aquifer interaction using 

SFR2 as compared to the instantaneous recharge approach 

using MODFLOW. The instantaneous approach simulates the 

saturation of the vadose zone after the onset of streamflow 

from the bottom up, where the water table must fill a volume 

of the pore space equal to the specific yield. SFR2 fills the 

pore space in the vadose zone from the top down, and because 

flow is assumed vertical in the vadose zone, the saturated 

region beneath the stream will be relatively narrow. This is 

in contrast to saturating the vadose zone from the bottom 

up where a significant proportion of ground-water flow will 

be horizontal and a much greater amount of pore space will 

need to be filled in order to saturate the vadose zone beneath 

the stream. Thus, intermittent streams that are hydraulically 

disconnected from the water table by several meters will be 

much more likely to reconnect to the water table using SFR2 

as compared to the instantaneous approach in MODFLOW. 

Reconnection between a stream and the water table can reduce 

the vertical hydraulic gradient and decrease the amount of 

seepage loss from a stream.

Purpose and Scope

The purpose of this report is to describe a method for 

simulating gravity-induced unsaturated flow beneath streams 

using MODFLOW. A description of the equations used in the 

kinematic-wave approximation for simulation of unsaturated 

flow beneath streams using MODFLOW was presented by 

Niswonger and Prudic (2004). This report summarizes that 

description, describes the assumptions and limitations of 

the method, presents the results from three test simulations, 

describes how the kinematic-wave approximation was added 

to SFR1 and MODFLOW-2000 (Harbaugh and others, 2000), 

provides a description of the data input for SFR2, and presents 

an example of data input and selected model output for the 

third test simulation. Details associated with streamflow 

routing are the same as in SFR1 (Prudic and others, 2004). 

Users of SFR2 should refer to the SFR1 documentation for 

information on these details.

Acknowledgments

The U.S. Geological Survey’s Ground-Water Resources 

Program funded this work. The original idea was developed 

through many conversations with Jon P. Fenske and Arlen 

D. Feldman of the Army Corps of Engineers Hydrologic 

Research Center, Davis, California and with Stanly A. Leake 

and Richard W. Healy of the U.S. Geological Survey. We 

appreciate comments by Paul M. Barlow, Arlen W. Harbaugh, 

and Eric D. Swain (U.S. Geological Survey) and Gordon 

Thrupp (S.S. Papadopulos and Assoc. Inc.).

Theory and Conceptualization of 

Unsaturated-Zone Flow 

Method

Vertical seepage through a homogeneous unsaturated 

zone can be approximated with kinematic waves (Smith, 

1983; Smith and Hebbert, 1983). This approximation is made 

by simplifying Richards’ equation. The one-dimensional 

Richards’ equation for vertical flow is:

∂

∂ = ∂

∂ = ∂ ∂

∂ − 















θ

∂

θ θ θ

t

q

z z

D

z

( ) ( ) K , (1)

where q is the volumetric water content; z is the elevation in 

the vertical direction (m); D( )q is the hydraulic diffusivity; 

K( )q is the unsaturated hydraulic conductivity; and t is time. 

Equation 1 is simplified such that vertical flux (q) is only 

driven by gravitational forces. Vertical flux is

q K = − ( )θ , (2)

where q is shown as positive upward. Substituting equation 2 

into equation 1 and eliminating the diffusive term yields:

∂

∂ +

∂

∂ = θ θt

K

z

( ) 0 . (3)

Theory and Conceptualization of Unsaturated-Zone Flow 3
Applicability of Assumptions

The method of characteristics solution to equation 3 is 

done by rewriting equation 3 so that q is taken partially in 

terms of t and z :

∂

∂ = − ∂∂

∂

∂

θ θ

θ

θ

t

K

z

( ) . (4)

Dividing both sides of equation 4 by ∂

∂

θ

z

gives the

characteristic solution for the kinematic-wave equation for 

unsaturated flow:

dz

dt

K = v d ∂

∂ = = ( ) ( ), θ

θ

θ θ 0 , (5)

where v( )q is the characteristic velocity, and q is a constant 

along a characteristic path. The characteristic velocity, which 

has positive and negative characteristics in the case of open 

channel flow, is restricted to the downward direction for the 

movement of a wetting front because the force of gravity can 

only cause the wetting front to move downward.

The derivative ∂

∂

K( )θ

θ is discontinuous over a sharp

wetting front that results from neglecting the diffusive term 

in equation 1 such that the spatial derivative in equation 3 is 

not defined in the absence of hydraulic diffusion. An analytic 

solution for the velocity may still be derived by considering 

the effects of diffusion and substituting an equivalent sharp 

wetting front of equal mass (Smith, 1983; Charbeneau, 1984). 

The solution to equation 1 for a wetting front that considers 

hydraulic diffusion can be found by integrating over a control 

volume containing a single wetting front (Charbeneau, 1984):

d

dt

dz K D

z z

z

z

z

θ θ θ θ+ −

∂

∂ ∫ ( ( ) ( ) ) =1

2

1

2

0 , (6)

where z1 and z2 are points above and below the wetting front 

at distances far enough such that ∂

∂ ≈θ

z

0 , respectively. 

D( )q can be neglected because ∂

∂ ≈θ

z

0 at z1 and z2 such that:

d

dt

dz K K

z

z

∫ θ θ θ + − = ( ) ( ) 2 1

1

2

0 , (7)

where K( ) θ1 and K(θ )2 are the values of K( )q at depths z1

and z2. Integrating over a profile containing a sharp front with 

equivalent mass gives:

θ θ θ

z

z

dz z z z z f f

1

2

∫ = − + − 1 1 2 2 ( ) ( ) , (8)

where zf is the depth of the sharp front. Combining equations 

7 and 8 gives:

dz

dt

u f K K = = s

( )− ( )

− ( , ) θ θ

θ θ

θ θ 1 2

1 2

1 2

, (9)

where us is the velocity of a sharp wetting front; q1 is the 

volumetric water content above a wetting front; and q2 is the 

volumetric water content below a wetting front.

Waves are generated by changes in infiltration resulting 

in waves with differing velocities. An increase in the surfaceinfiltration rate will cause a wetting front to form, which 

is represented by a lead wave. A decrease in the surface 

infiltration rate will cause internal drainage to occur, which 

is represented by a trailing wave. Attenuation of a lead wave 

occurs as a trailing wave of higher velocity overtakes it. When 

a trailing wave overtakes a lead wave, the lead wave water 

content becomes equal to the water content of the trailing 

wave. Consequently, this reduces the velocity, water content 

and flux of the lead wave. Conversely, when a lead wave 

overtakes a trailing wave or another lead wave of smaller 

magnitude, the overtaken wave(s) is (are) removed, and 

the water content and flux of the uppermost lead wave are 

maintained, resulting in rewetting.

In contrast to a wetting front that stays sharp due to 

gravity, trailing waves elongate with time due to gravity. 

Consequently, trailing waves must be either discretized into a 

series of incremental waves, or be represented by a function 

that describes internal drainage with time. In either case, 

internal drainage is determined based on the water content and 

relative location of points along a trailing wave. For example, 

points A and B shown in figure 2 represent two points along 

the trailing wave where the velocities may be determined. The 

velocity at points A and B is calculated from the characteristic 

solution (equation 5) and the ratio of their velocities is equal to 

the ratio of their depths:

z

z

u

u

dK d

dK d

B

A

B

A

B

A

= = =

=

( )

( )

( )/

( )/

θ

θ

θ θ

θ θ

θ θ

θ θ

, (10)

where qA and zA are the water content and depth of the 

trailing wave increment below land surface represented by 

point A in figure 2; and qB and zB are the water content 

and depth of trailing wave increment below land surface 

represented by point B in figure 2. 

4 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
In contrast to a lead wave that forms a discontinuity in the 

water content profile, the derivative, dK( )/ θ θ d 7, is continuous 

over a trailing front. Accordingly, the unsaturated hydraulic 

conductivity function described by Brooks and Corey (1966) 

can be used to evaluate dK d ( )/ θ θ7. The Brooks-Corey function 

for unsaturated hydraulic conductivity can be expressed as:

K Ks r

s r

( )θ θ θ

θ θ

ε

= −

−















 , (11)

where Ks is the saturated hydraulic conductivity; qr is the 

residual water content; qs is the saturated water content; and 

e is the Brooks-Corey exponent. 

The drainable water from the unsaturated zone, referred 

to as the specific yield or Sy, is equal to the porosity minus 

the specific retention (Sr), where Sr is the volume of water per 

unit volume of rock that is retained when the rock is drained 

by gravity (Meinzer, 1923, p. 28). Specific yield also is used 

by MODFLOW to calculate the amount of storage within an 

unconfined aquifer. Thus, continuity between the unsaturated 

zone and unconfined aquifers in MODFLOW is maintained 

through the specific yield by approximating Sr with qr. 

Residual water content generally is smaller than Sr, such that 

caution should be used when representing unsaturated flow 

variables in SFR2 with measured values. The value of qs-qr is 

represented by Sy for the remaining sections of this report. 

The velocity of the deepest point along a trailing wave is 

determined by taking the derivative of equation 11:

v

K

S S

s

y

r

y

=

 −



















−

ε θ θ ε 1

, (12)

where v and q are the velocity and water content at the 

deepest point along a trailing wave, respectively. A relation 

between the deepest point along a trailing wave and all other 

points along a trailing wave is determined by substituting the 

velocity defined by equation 12 into equation 10:

z z r

r

( )θ θ θ

θ θ 0

ε

= −

−

















−

0

1

, (13)

where z(θ) is the depth of a point on a trailing wave; q is the 

water content of a point on a trailing wave; and q0 and z0

are the water content and depth of the deepest point along a 

trailing wave, respectively.

Complications arise when routing waves resulting from 

rapid fluctuations in the infiltration rate. For example, a lead 

wave generated by a pulse of infiltration will decay as trailing 

waves overcome the lead wave. Equations 12 and 13 define the 

spreading of a trailing wave; however, they are not applicable 

subsequent to the intersection with a leading wave. Analytic 

relations can be derived for routing trailing waves during 

intersection with a lead wave but these equations are difficult 

to solve. A simpler approach is to discretize the trailing wave 

into steps and calculate velocities of each step based on a 

finite-difference approximation (Smith, 1983):

v

K K = ( ) ( ) θ θ θ − −

θ

∆

∆

, (14)

where ∆θ is the change in water content between two 

adjacent locations along a trailing wave.

As shown in figure 2, the finite-difference approximation 

of the true trailing wave does not conserve mass perfectly 

when a trailing wave intercepts a lead wave and vice versa. 

However, for problems involving rapidly varying infiltration 

rates, mass is conserved within an error of less than 0.05 

percent when a trailing wave is represented by 15 increments, 

and for most problems 10 increments are sufficient. The user 

specifies the number of increments used to represent a trailing 

wave with the SFR2 input parameter NSTRAIL.

Figure 2. A wetting front through a uniform column of unsaturated 

material affected by a decrease in surface flux and results from a 

kinematic-wave approximation of the wetting front using leading and 

trailing waves.

(Modified from Smith and Hebbert, 1983.)

Theory and Conceptualization of Unsaturated-Zone Flow 5

θi θt θ 1 t2

Trailing wave 

represented as a 

step function at 

time t2

Actual wetting front 

at time equals t1

Actual wetting front 

at time equals t2

Leading wave 

at time equals t2

 Leading wave 

at time equals t1

WATER CONTENT (θ)

DEPTH (z)

A

B
An algorithm was written in FORTRAN using equations 

9, 11, 12, 13, and 14 to simulate vertical unsaturated flow. 

This algorithm keeps track of the location, water content, 

velocity and flux of all waves in the unsaturated zone through 

time and the interaction of waves overtaking one another. 

Additionally, changes in the thickness of the unsaturated zone 

are considered with respect to waves reaching the water table 

and contributing to ground-water storage.

Boundary conditions used to represent infiltration from 

the land surface are formulated as an infiltration rate such that 

the specified flux at the surface is converted to water content 

in order to apply the characteristic solutions for unsaturated 

flow. A relation between water content and infiltration using 

equations 2 and 11 is:

θ θ

ε

= 



















 ( ) + < ≤ in

s

y r in s

q

K

S q K

1

0

/

, (15a)

θ θ= s K q s in < , (15b)

where q is the water content of a wave generated from 

infiltration, and qin is the vertical flux or infiltration rate 

(length per time). The water content is set to the saturated 

water content when the infiltration rate exceeds the saturated 

vertical hydraulic conductivity.

In gravity-dominated flow, the flux in the unsaturated 

zone is equal to the unsaturated vertical hydraulic conductivity 

and cannot increase above the value of the saturated hydraulic 

conductivity. Consequently, seepage across the streambed 

in SFR2 is limited by the underlying vertical hydraulic 

conductivity in the unsaturated zone. This differs from SFR1, 

in which seepage is only limited by streambed conductance 

and the head difference between the stream and aquifer, or 

when the stream is hydraulically disconnected from the water 

table, by the head difference between the stream and the 

bottom of the streambed. 

As with SFR1, numerical oscillations in SFR2 can occur 

when the streambed area is large relative to the corresponding 

model cell plan-view area. For these conditions, recharge from 

streambed seepage may be greater than storage and outflow 

from the model cell. To minimize oscillations, seepage in 

SFR2 is limited to the vertical hydraulic conductivity times 

the area of the model cell. When selecting the dimensions of 

model cells corresponding to stream reaches, oscillations are 

minimized when the plan-view area of a model cell is greater 

than the streambed area. This problem also may be alleviated 

by allowing seepage from a stream reach to be distributed to 

more than one MODFLOW cell. However, this capability is 

not presently available in SFR2.

Surface-Water Flow and Infiltration

Streamflow routing in SFR2 is done in the same manner 

as in SFR1 (Prudic and others, 2004). Streamflow is routed 

based on the continuity equation assuming steady, uniform 

flow such that the volumetric inflow is equal to the outflow 

minus all sources and sinks to the channel. 

In SFR2 (as in SFR1), flows are routed through a network 

of channels where flow is always in the same direction 

along channels, and where seepage (ground-water recharge 

or discharge) is constant for each MODFLOW time step. 

Ground-water recharge resulting from streambed seepage and 

flow through the unsaturated zone may be variable within 

a MODFLOW time step, because the time steps used in 

calculating flow through the unsaturated zone may be shorter 

than those used in the saturated zone. Thus, all water reaching 

the water table from the unsaturated zone is totaled over the 

MODFLOW time step. 

When the stream is separated from the water table by 

an unsaturated zone, the infiltration rate from streams is 

calculated using Darcy’s Law and assuming uniform flow 

between a given section of stream and a corresponding 

volume of aquifer. The seepage out of the stream when the 

flux does not exceed the vertical hydraulic conductivity of the 

unsaturated zone is computed using:

q

K

m = − h h s b ( ), (16)

where q is the infiltration rate; K is the saturated hydraulic 

conductivity of streambed sediments; hs is the hydraulic head 

in stream determined from top of streambed sediments; hb is 

the hydraulic head at the bottom of the streambed; and m is 

the thickness of streambed sediment as defined by the user.

The volume of water that seeps from a stream is 

calculated by multiplying the infiltration rate by the wetted 

area of the stream. The wetted area of the stream may be 

held constant or determined on the basis of the stream crosssectional dimensions, discharge, and stage. The relation 

between stage and discharge is calculated using Mannings 

equation. The option to simulate unsaturated flow between 

streams and aquifers may be used for a stream of constant 

wetted area and width or for an eight-point cross-section 

(fig. 3). The unsaturated zone may be discretized if a variable 

cross-section is specified (ICALC= 2). The unsaturated zone 

beneath a stream reach is discretized into ISUZN (as specified 

in the SFR2 input file) compartments. Water is routed 

independently within each of the unsaturated compartments 

(fig. 3). The option to simulate unsaturated flow beneath 

streams is not available when using ICALC = 0, 3, or 4. 

6 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Coupling Unsaturated-Zone Flow in SFR2 to 

MODFLOW

Water draining from the unsaturated zone is added to the 

ground-water flow equations in SFR2 when the water table 

is below the streambed. This approach differs from SFR1 in 

which seepage from the streambed instantly recharges the 

water table. Water leaving the unsaturated zone is subtracted 

from the right hand-side vector: 

AX B Q = − SFR2 , (17)

where A is a matrix containing the coefficients of the 

conductance equations (HCOF array) that are solved by 

MODFLOW-2000; X is a one-dimensional vector containing 

the ground-water heads that are solved by MODFLOW-2000; 

B is a one-dimensional vector containing all known terms in 

the conductance equations that are not multiplied by unknown 

head values (RHS array); and QSFR2 is the volumetric rate 

(volume per time) of recharge computed in a given model cell 

from SFR2. The layer, row, and column for each stream reach 

and underlying unsaturated zone are specified in the model 

input (Appendix 1). If the cell corresponding to the stream 

reach is inactive, the program will continue to route water 

in the unsaturated zone until it reaches the water table of the 

uppermost active cell. However, the hydraulic properties of 

the unsaturated zone (Ks , Sy , qs , e ) will remain constant even 

if the unsaturated zone spans more than one cell in the vertical 

direction. 

The quantity of water that recharges an active model 

cell from the unsaturated zone is dependent on the degree of 

fluctuation of the water table. If the water table is rising into 

an unsaturated zone that is near saturation, then the water 

table rise will be exaggerated because there is little empty 

pore space available to be filled by the rising water table. 

Alternately, if the water table is falling due to ground-water 

discharge or pumping at some other location, then it will take 

longer for a wetting front from the unsaturated zone to reach 

the water table because the water content is at residual in the 

interval over which the water table declined. 

Recharge resulting from a rising water table over a 

time step is calculated by first routing waves through the 

unsaturated zone using the base of the unsaturated zone as the 

water-table elevation at the end of the previous time step. Flow 

across the base is accumulated for each time step. The volume 

of water in the unsaturated zone through which the water 

table rose over a time step is added to the volume that flowed 

past the base. This sum is divided by the time step to obtain 

a volumetric recharge rate. Once the flow equation has been 

solved for the time step, a new base of the unsaturated zone is 

computed in each active cell.

Compartment 

width 

Compartment 

height 

Water table 

Streambed 

Stream 

Unsaturated zone 

compartment 

MODFLOW column or row width

(larger than maximum stream width)

Maximum stream width 

Figure 3. Discretization of unsaturated zone beneath a stream in a single MODFLOW cell. 

The unsaturated zone is discretized into seven compartments defined by eight points in this example. Dashed lines represent stream depth and infiltration 

when stream is at maximum width. 

Theory and Conceptualization of Unsaturated-Zone Flow 7
Recharge resulting from a declining water table over a 

time step is computed in a similar manner, except that the 

interval over which the water table declined during an iteration 

is set to the residual water content and waves are routed 

through the base of the unsaturated zone computed from the 

previous iteration. The volume of water that flowed across 

the base of the unsaturated zone is divided by the time step to 

obtain a volumetric recharge rate. 

A water balance for the entire unsaturated zone is tracked 

independently of the ground-water balance equations within 

MODFLOW. The mass balance for the unsaturated zone 

accounts for all water between the streambed and the water 

table. Water passing across the water table is accounted for in 

the ground-water budget. When the water table rises above the 

bottom of the streambed in a model cell, then the unsaturated 

zone mass balance ceases to be calculated for that cell, and the 

unsaturated zone storage and change in storage for that cell are 

zero. 

Steady-State Simulations with SFR2

A user may indicate that a stress period is steady state by 

using the character variable Ss/tr within the Discretization 

file (Harbaugh and others, 2000, p. 45). Steady-state 

simulations with SFR2 assume the infiltration rate through the 

streambed is equal to the ground-water recharge rate. SFR2 

will use the unsaturated flow conditions at the end of a steadystate stress period as initial conditions for subsequent transient 

stress periods. Steady-state water content profiles may be 

printed to a file for specified locations using the Gage Package 

(Appendix 1). 

Assumptions and Limitations

The kinematic-wave approximation does not account 

for the increased flux at the onset of wetting due to negative 

pressure gradients. Infiltration into initially dry material can 

be rapid due to large negative pressure gradients. However, 

gravity potential gradients dominate after the wetting front 

has advanced sufficiently far away from the surface, and 

infiltration slows to a value equivalent to the unsaturated 

hydraulic conductivity (Childs and Bybordi, 1969; Mein 

and Larson, 1973). Negative pressure gradients may also 

play a role in the flux of water during lateral and vertical 

redistribution in the unsaturated zone. Because negative 

pressure gradients are ignored, SFR2 likely under predicts 

the advancement of the wetting front at early times but over 

predicts during drying or redistribution.

The present form of the model does not simulate 

unsaturated flow through multiple layers of varying hydraulic 

properties. The model can simulate unsaturated flow through 

inactive cells to deeper active model cells. However, the 

hydraulic properties of an unsaturated column between the 

streambed and the uppermost active cell remain equal to their 

initial values during the simulation. 

If the vertical hydraulic conductivity for the 

unsaturated zone is specified within the Layer Property 

Flow (LPF) Package input file by setting ISFROPT= 2 

or 4, then the vertical hydraulic conductivity specified for 

the uppermost active cell beneath a stream reach will be 

used for the unsaturated zone for the entire simulation. If 

the user is representing confining beds using the vertical 

hydraulic conductivity defined in the LPF input file, then 

it is recommended that the user specify separate hydraulic 

conductivity values for the unsaturated zone within the SFR2 

input file by setting ISFROPT= 3 or 5. 

SFR2 can simulate streambed seepage in partially 

penetrating streams; however, due to the Dupuit approximation 

for simulating unconfined ground-water flow, errors may 

occur in the model. Anderson (2004) showed that when 

effective or calibrated streambed conductance is used to 

restrict seepage, errors in forecasted simulation results will be 

minimized. Seepage in partially penetrating streams also may 

be biased in regional models where a model cell area is much 

larger than the corresponding stream reach area. In this case, 

treating the streambed conductance as a calibration parameter 

improves model predictability. 

Results and Discussion

Model Comparisons 

Results from the kinematic-wave approximation in SFR2 

were compared with results from the USGS Variably Saturated 

Two-Dimensional Flow and Transport (VS2DT) Model 

(Healy, 1990) for two different problems. VS2DT solves 

Richards’ equation numerically in terms of pressure potential. 

Although VS2DT can simulate two-dimensional unsaturated 

flow, the kinematic-wave approximation is limited to onedimensional vertical flow. Thus, the two test simulations were 

limited only to vertical flow. The stream was represented by 

boundary conditions in both models and streamflow was not 

explicitly simulated.

The first test simulation consisted of simulating a 100-m 

tall column of materials having hydraulic properties typical 

of a sandy clay loam. The same hydraulic properties that 

relate unsaturated hydraulic conductivity to water content (K

and θ) and the capillary potential to water content (y and θ) 

were used in the simulation with SFR2 and VS2DT. Table 1 

lists the input variables used for the first test simulation. The 

simulations consisted of a variable-flux boundary at the top of 

the column and a zero-pressure head boundary at the bottom. 

8 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40

0 25 50 75 100

0

10

20

30

40

50

60

70

80

90

100

0

0.1

0.2

0.3

0.4

50 days

75 days

100 days

SURFACE FLUX, 

IN METER PER DAY 

DEPTH, IN METERS

WATER CONTENT, IN CUBIC METER PER CUBIC METER

TIME, IN DAYS 

Kinematic-wave approximation

Numerical solution of Richards equation

A. Flux-boundary condition assigned in simulation.

B. Distribution of vertical water content after 50, 75, and 100 days.

Figure 4. Comparison of results from the kinematic-wave approximation used in the Streamflow-Routing 

(SFR2) Package with the Variably Saturated Two-Dimensional Flow and Transport (VS2DT) Model (Healy, 

1990) .

Table 1. Hydraulic properties assigned in first test simulation assuming a column filled with sandy clay loam.

Saturated vertical

hydraulic conductivity

(meter per second)

Saturated water 

content

(cubic meter 

per meter)

Specific yield

(cubic meter 

per meter) 

Brooks-Corey 

exponent

(dimensionless)

Air entry 

pressure

 (meter)

3.4 ×10-6 0.35 0.27 3.5 -0.1

Results and Discussion 9
The initial water content was set equal to the difference 

between saturated water content and specific yield throughout 

the model. Figure 4 shows a comparison of the water 

content profile after 50, 75, and 100 days of simulation. The 

advancement of the wetting front in the simulation with SFR2 

is nearly the same as that with VS2DT (showing minimal 

difference at 50 and 75 days), whereas, during a drying period 

(75 to 100 days), the advancement of the wetting front is 

slightly greater in the simulation with SFR2 than that with 

VS2DT (see difference at 100 days). This difference suggests 

that negative pressure gradients resulting from infiltration 

into dry materials are attenuated within the first 50 days of 

infiltration and that during the last 25 days, capillary potentials 

slightly retard the advancement of the wetting front. 

The second test simulation was designed to assess 

the ability of SFR2 to simulate streambed infiltration into 

an unsaturated zone with variable stream depths. In this 

simulation, the saturated hydraulic conductivity of the 

streambed is less than that of the unsaturated zone. Table 2

lists the input variables used for the simulation. Results of 

fluxes calculated by SFR2 and VS2DT for four different 

stream depths are listed in table 3 and are in reasonably good 

agreement. The corresponding advancement of the wetting 

front after 11.6, 23.2, and 34.7 days of infiltration shows good 

agreement between the two models (fig. 5). 

Table 2. Input variables for second test simulation assuming vertical seepage through unsaturated zone that is restricted by 

a semipervious streambed.

[na, not applicable]

Saturated vertical 

hydraulic conductivity

(meter per second)

Saturated water 

content

(cubic meter

per meter)

Specific yield 

(cubic meter

per meter) 

Brooks-Corey 

exponent

(dimensionless)

Air entry 

pressure

 (meter)

Streambed 1 × 10-7 na na na na

Unsaturated zone 1 × 10-6 0.30 0.2 3.6 -0.20

Table 3. Comparison of calculated fluxes through a semipervious streambed resulting from four different stream stages 

from kinematic-wave approximation in the Streamflow-Routing (SFR2) Package with those from the Variably-Saturated TwoDimensional Flow and Transport (VS2DT) Model (Healy 1990). Fluxes correspond to the simulation time immediately following a 

change in the stream depth.

Stream depth (meter) 0.70 2.46 1.59 1.04

Flux (meter per second; SFR2) 9.54 × 10-7 1.83 × 10-6 1.40 × 10-6 1.12 × 10-6

Flux (meter per second; VS2DT) 9.44 × 10-7 1.78 × 10-6 1.37 × 10-6 1.11 × 10-6

10 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
0.05 0.10 0.15 0.20 0.25 0.30 0.35

WATER CONTENT, IN CUBIC METER PER CUBIC METER 

0

5

10

15

20

25

30

35

DEPTH BELOW TOP OF STREAMBED, IN METERS 

Kinematic-wave approximation 

Numerical solution of Richards’ equation

11.6 days

23.2 days

34.7 days

Figure 5. Comparison of results from the kinematic-wave approximation in the Streamflow-Routing (SFR2) Package with the 

Variably Saturated Two-Dimensional Flow and Transport (VS2DT; Healy, 1990) Model for seepage through a one-dimensional 

column filled with a silt loam. 

Variations in water content below the streambed (thickness of 1 meter) are the result of different calculated flux rates across the 

streambed caused by varying stream depth (from Niswonger and Prudic, 2004, fig. 4).

Results and Discussion 11
Hypothetical Test Simulation

A hypothetical aquifer was modeled assuming that an 

unsaturated zone separated the water table from a streambed 

(fig. 6). The hypothetical aquifer was divided into 7 cells in 

the y-direction (along columns), each 90-m long; 200 cells 

in the x-direction (along rows), each 100-m long; and one 

layer in the z-direction, 140–110-m thick (fig. 6A). A stream 

was placed along the middle row in the x-direction and had a 

maximum width at peak flow of 14 m and a depth of 6 m and 

was represented by an eight-point cross-section. A total of 200 

stream reaches (corresponding to the 200 aquifer cells in the 

x-direction) connected the individual cells. 

The initial water table was set at 104 m, resulting in an 

initial unsaturated zone that was 36 m and 6 m, respectively, 

below the streambed at the upstream and downstream ends of 

the model. No-flow boundaries were assigned along all sides 

of the model, except the downstream side where a generalhead boundary was used. The unsaturated zone for each cell 

with a stream was divided into 7 vertical sections beneath the 

stream (fig. 6B). Unsaturated flow was not considered in any 

other cells. 

Aquifer flow properties were considered homogeneous 

throughout the model, with a saturated horizontal hydraulic 

conductivity of 8 × 10-4 m/s, a saturated vertical hydraulic 

conductivity of 6 × 10-6 m/s, and a specific yield of 0.20. The 

streambed and unsaturated zone hydraulic properties are listed 

in table 4. Streamflow into the first reach was specified using 

monthly-averaged streamflows from the Cosumnes River in 

California. The input data files and selected output files are 

listed in Appendix 2.

Streambed thickness was varied between 0.5 and 1.5 m 

and the resulting effects on seepage, recharge, stream depth 

and water table beneath a selected model cell midway along 

the stream (stream reach 100) were analyzed. Because 

the streambed had a lower hydraulic conductivity than the 

unsaturated zone, seepage, recharge, and the water table 

elevation decreased with increased streambed thickness 

(fig. 7). Seepage was more sensitive to streambed thickness for 

greater stream depths (fig. 7C). 

This simulation suggests that seepage is sensitive to 

stream depth when the streambed hydraulic conductivity 

is lower than that of the materials in the unsaturated zone. 

However, stream depth becomes less important and the wetted 

area of the stream controls seepage when the streambed 

hydraulic conductivity is higher than that of the unsaturated 

zone. 

The example was limited to unsaturated flow beneath a 

disconnected stream whereby seepage occurred only from the 

stream into the streambed. SFR2 can also be used to simulate 

streams where an unsaturated zone routinely appears and 

disappears beneath a stream (Niswonger and Prudic, 2004), or 

to simulate ground-water discharge to streams as demonstrated 

with SFR1 (Prudic and others, 2004).

Table 4. Unsaturated zone input variables used in example problem.

[Abbreviations: m, meter; m/s, meter per second; m3/m3, cubic meter per cubic meter; na, not applicable]

Model input parameter Semipervious 

streambed

Unsaturated

 zone

Saturated vertical hydraulic conductivity, Ksv (m/s) 5 × 10-7 6 × 10-6

Saturated water content, θs (m3/m3) na .30

Specific yield, Sy (m3

/m3) na .20

Brooks-Corey exponent, ε na 3.5

12 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
A. Physical boundaries of modeled area.

B. Division of unsaturated zone into seven compartments beneath stream channel 

 (modified from Niswonger and Prudic, 2004).

Streambed

Stream

Water table

Stream 

0

0

140

104

0

0 14

6

Water table 

Unsaturated zone 

Unsaturated zone 

compartment

Saturated zone 

DISTANCE ALONG ROWS, IN KILOMETERS 

HORIZONTAL DISTANCE FROM LEFT EDGE OF STREAM, IN METERS

VERTICAL DISTANCE FROM 

THALWEG, IN METERS

DISTANCE 

ALONG COLUMNS, 

IN METERS 

VERTICAL DISTANCE FROM BASE 

OF MODEL, IN METERS 

20 

630 

0 

110

104

Figure 6. Conceptualization of hypothetical test simulation.

Results and Discussion 13
0

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.016

Streambed 1 meter thick

Streambed 1.5 meters thick

Streambed 0.5 meter thick

Water table when streambed 1 meter thick

Water table when streambed 1.5 meters thick

Water table when streambed 0.5 meter thick

100

104

108

112

116

120

124

128

132

0 300 600 900

TIME, IN DAYS

SEEPAGE, IN CUBIC 

METERS PER SECOND 

RECHARGE, IN CUBIC

METERS PER SECOND

1,200 1,500

Stream stage

Streambed

0

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.016

VERTICAL DISTANCE FROM BASE 

OF MODEL, IN METERS

A. Seepage rate across streambed.

B. Recharge rate at water table.

C. Water table for a stream reach midway along the stream.

Figure 7. Effects of variations in thickness of a semipervious streambed (from Niswonger 

and Prudic, 2004).

14 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Implementation of SFR2 with 

MODFLOW-2000

SFR2 consists of twenty subroutines that are used for the 

Ground-Water Flow (GWF1) Process of MODFLOW-2000 

(Harbaugh and others, 2000), three subroutines are used 

to read and write parameters and non-parameters, and one 

subroutine checks for incompatibilities with the Sensitivity 

Process in MODFLOW-2000. The unsaturated-flow 

computations were included within the individual modules of 

SFR2. Each subroutine or module incorporated in the main 

listing of MODFLOW-2000 is identified first by the GWF1 

Process, followed by SFR2 (version 2 of the StreamflowRouting Package), and then followed by the division of 

procedures used in MODFLOW-2000. Specifically, subroutine 

functions include: (1) memory allocation; (2) reading input 

data and initialization of variables at the beginning of the 

simulation; (3) setting values for the first stress period; (4) 

reading input data and initializing variables at the beginning 

of each stress period; (5) formulating characteristic and 

finite-difference equations; and (6) computing mass balance. 

These procedures occur within the following modules, 

GWF1SFR2ALP, GWF1SFR2RPP, GWF1SFR2RPS, 

GWF1SFR2AD, GWF1SFR2FM, and GWF1SFR2BD, 

respectively. 

Three additional modules, GWF1SFR2UHC, 

GWF1SFR2UZOT and GWF1SFR2SSET were added to the 

main program. The module GWF1SFR2UHC was added to 

define the saturated vertical hydraulic conductivity in the 

unsaturated zone, either as a variable that is read into SFR2 or 

obtained from the LPF Package. When the vertical hydraulic 

conductivity for the unsaturated zone beneath a stream reach is 

specified within the LPF Package, the hydraulic conductivity 

of the uppermost active cell at the beginning of the simulation 

will be used throughout the entire simulation. The module 

GWF1SFR2UZOT was added for the printing of an overall 

mass balance of the unsaturated zone beneath streams. The 

subroutine GWF1SFR2SSET was added to initialize wave 

characteristics after a steady-state simulation.

In addition to the modules from the main program, 

several new modules were added to SFR2 for routing 

kinematic waves through the unsaturated zone. These 

modules are named ROUTWAVESIT, ROUTWAVEST, 

UZFLOW, LEADWAVE, TRAILWAVE, CALC_UNSAT_

INFIL, UNSAT_MASS_BAL, and CHANNELAREA. 

Additionally, three GAGE Package modules (GWF1GAG5DF, 

GWF1GAG5I, and GWF1GAG5SO) were modified to print 

time series data of unsaturated zone infiltration, storage, and 

ground-water recharge, and to print time series of watercontent profiles in the unsaturated zone. 

Additionally, the MAIN program of MODFLOW-2000 

has been modified to designate “SFR” in position 44 in 

CUNIT. Whenever “SFR” is used at the beginning of a line in 

the NAME file, the program will use the specified unit number 

and file name to obtain information for SFR2 and will execute 

MODFLOW-2000 using SFR2.

Data input from SFR1 can be read into SFR2 without 

modification. Output from SFR2 can be printed: (1) to the 

main listing file; (2) to a separate listing file when requested 

in the SFR2 input; and (or) (3) to separate files for each 

specified stream reach that provides a time series of selected 

information using the Gage Package. Two new options were 

added to the Gage Package (Appendix 1) for printing time 

series data of the unsaturated zone beneath a selected reach. 

Option 6 will print time, stage, ground-water head, streambed 

seepage, change in unsaturated storage and recharge to ground 

water. Option 7 will print time and the unsaturated water 

content profiles for the selected stream reach. 

SFR2 has been programmed to apply leakage from a 

stream reach to a designated underlying cell. If the water level 

in the underlying cell becomes less than the bottom elevation 

of the cell, the cell is inactivated. When this happens, SFR2 

will begin searching vertically below the designated cell to 

find the uppermost active cell. If it finds an active cell in 

the vertical column, recharge at the water table is applied to 

that cell, however, the unsaturated zone hydraulic properties 

remain unchanged. If all cells in the vertical column beneath 

a stream reach are inactive, no interaction with ground water 

is allowed and flow in the stream reach is passed to the next 

reach. 

Summary

The U.S. Geological Survey modular ground-water model 

(MODFLOW-2000) has been used to model surface- and 

ground-water interaction over a wide range of spatial scales. 

Because of increased interest in the connection of surface 

water and ground water, the capability to simulate unsaturated 

flow beneath streams has been added to the StreamflowRouting (SFR2) Package of MODFLOW-2000. Unsaturated 

flow beneath streams is simulated using a kinematic-wave 

approximation and assumes that downward flow beneath 

streams is a result of gravity. SFR2 allows for the simulation 

of flow and storage in the unsaturated zone, and for the 

simulation of a time delay for recharge, while still maintaining 

the applicability of MODFLOW-2000 to basin-scale problems.

Summary 15
Unsaturated flow is simulated independently of saturated 

flow within each model cell that corresponds to a stream reach 

in which the water table is below the bottom of the streambed. 

The relation between the unsaturated hydraulic conductivity 

and water content in the unsaturated zone is defined on the 

basis of the Brooks-Corey function. Variables used to simulate 

unsaturated flow include saturated and initial water contents; 

saturated vertical hydraulic conductivity; and the BrooksCorey exponent. The values of variables can be different for 

each stream reach. 

In the simulation of unsaturated flow beneath streams, 

seepage flux (flow per unit area) across the streambed is 

limited to the vertical saturated hydraulic conductivity. This 

differs from the first version of the Streamflow-Routing 

(SFR1) Package in which the seepage losses are limited by the 

conductance term of the streambed and the head difference 

across the streambed. Additionally, when an unsaturated zone 

is simulated beneath a stream, recharge is not added to ground 

water unless water in the unsaturated zone drains to the water 

table. This approach also differs from SFR1, which assumes 

that water instantly recharges ground water.

The representative width of the stream may be held as 

a constant or determined on the basis of the stream cross 

section, discharge, and depth. When the water table rises, 

the quantity of water stored in the unsaturated zone above 

the residual water content is added to recharge. Thus, the 

water table can rise rapidly when the unsaturated zone is at 

or near saturation. A water budget for the unsaturated zone 

is calculated independently of the ground-water budget in 

MODFLOW-2000. 

Two test simulations were used to compare results from 

the kinematic-wave approximation in SFR2 with the U.S. 

Geological Survey’s Variably Saturated Two-Dimensional 

Flow and Transport (VS2DT) Model. The first test simulation 

checked the capability of SFR2 to simulate the advancement 

of a wetting front following a period of redistribution (drying), 

and the second test simulation checked the accuracy of 

representing a semipervious streambed. The simulations 

using SFR2 for both cases were in good agreement with the 

magnitude and downward progression of the wetting front 

with respect to the results from VS2DT. A third simulation 

assumed a stream was separated from ground water by an 

unsaturated zone. The simulation tested effects of different 

stream depths and streambed thicknesses on the timing of 

recharge to the aquifer. Results indicate that stream depth 

affects seepage loss when the hydraulic conductivity of the 

streambed is less than that in the unsaturated zone. However, 

stream depth becomes less important and the wetted area of 

the stream controls seepage loss when the streambed hydraulic 

conductivity is higher than that of the unsaturated zone.

The kinematic-wave approximation to unsaturated 

flow in SFR2 ignores the effects of negative potential 

gradients. Ignoring negative potential gradients during the 

onset of streamflow under predicts initial seepage across the 

streambed. Ignoring negative potential gradients after flow 

ceases over predicts the advancement of the wetting front. 

Seepage across a streambed also can exhibit a horizontal 

component of flow not accounted for in SFR2. For streambeds 

in which significant lateral flow occurs, SFR2 may over 

predict the advancement of the wetting front. Lastly, the model 

can simulate unsaturated flow through layers of cells but the 

hydraulic properties assigned to the unsaturated zone remain 

constant during the simulation. 

16 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
References Cited

Anderson, I.A., 2004, Modeling ground-water—surface-water 

interaction using the Dupuit approximation: Advances in 

Water Resources, v. 28, no. 4, p. 315-327.

Brooks, R.H., and Corey, A.T., 1966, Properties of porous 

media affecting fluid flow: American Society of Civil 

Engineers, Journal of Irrigation and Drainage, v. 101, 

p. 85-92.

Charbeneau, R.J., 1984, Kinematic models for soil moisture 

and solute transport: Water Resources Research, v. 20, no. 

6, p. 699-706.

Childs, E.C., and Bybordi, M., 1969, The vertical movement 

of water in stratified porous material, 1. Infiltration: Water 

Resources Research, v. 5, no. 2, p. 446-459.

Danskin, W.R., and Hanson, R.T., 2002, Appendix 2: 

Documentation and description of changes made to the 

streamflow-routing package in the MODFLOW groundwater flow model, in Hanson, R.T., Martin, Peter, and 

Koczot, K.M., Simulation of ground-water/surface-water 

flow in the Santa Clara-Calleguas Basin, Ventura County, 

California: U.S. Geological Survey Water-Resources 

Investigations Report 02-4136, p. 187-193. 

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—modularization concepts 

and the ground-water flow process: U.S. Geological Survey 

Open-File Report 00-92, 121 p.

Healy, R.W., 1990, Simulation of solute transport in variably 

saturated porous media with supplemental information on 

modification of the U.S. Geological Survey’s computer 

program VS2D: U.S. Geological Survey Water-Resources 

Investigations Report 90-4025, 125 p.

Konikow, L.F., Goode, D.J., and Hornberger, G.Z., 1996, A 

three-dimensional method-of-characteristics solute-transport 

model (MOC3D): U.S. Geological Survey Water-Resources 

Investigations Report 96-4267, 87 p. 

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground-water flow 

model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap. A1, 586 p.

Mein, R.G., and Larson, C.L., 1973, Modeling infiltration 

during a steady rain: Water Resources Research, v. 9, no. 2, 

p. 384-394.

Meinzer, O.E., 1923, Outline of ground-water hydrology, with 

definitions: U.S. Geological Survey Water-Supply Paper 

494, 71 p.

Merritt, M.L., and Konikow, L.F., 2000, Documentation of a 

computer program to simulate lake-aquifer interaction using 

the MODFLOW ground-water flow model and the MOC3D 

solute-transport model: U.S. Geological Survey WaterResources Investigations Report 00-4167, 146 p. 

Niswonger, R.G., and Prudic, D.E., 2004, Modeling variably 

saturated flow using kinematic waves in MODFLOW, 

in Hogan, J.F., Phillips, F.M., and Scanlon, B.R., eds., 

Groundwater Recharge in a Desert Environment: The 

Southwestern United States: American Geophysical Union, 

Washington, D.C., Water Science and Application Series, v. 

9, p. 101-112. 

Prudic, D.E., 1989, Documentation of a computer program to 

simulate stream-aquifer relations using a modular, finitedifference, ground-water flow model: U.S. Geological 

Survey Open-File Report 88-729, 113 p.

Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A new 

streamflow-routing (SFR1) package to simulate streamaquifer interaction with MODFLOW-2000: U.S. Geological 

Survey Open-File Report 2004-1042, p. 95. 

Smith, R.E., 1983, Approximate sediment water movement by 

kinematic characteristics: Soil Science Society of America 

Journal, v. 47, p. 3-8.

Smith, R.E., and Hebbert, R.H.B., 1983, Mathematical 

simulation of interdependent surface and subsurface 

hydrologic processes: Water Resources Research, v. 19, 

no. 4, p. 987-1001.

References Cited 17
Appendix 1: Data Input Instructions for the Streamflow-Routing (SFR2) Package 

with Unsaturated Flow Beneath Streams. 

MODFLOW Name File

The Streamflow-Routing Package is activated automatically by including a record in the MODFLOW name file using the 

file type (Ftype) “SFR” to indicate that relevant calculations are to be made in the model and to specify the related input data 

file. The user can optionally specify that stream gages and monitoring stations are to be represented at one or more locations 

along a stream channel by including a record in the MODFLOW name file using the file type (Ftype) “GAGE” that specifies 

the relevant input data file giving locations of gages. The modifications in SFR2 do not require any changes to the data input for 

SFR1. SFR2 is compatible with MODFLOW-2000 (Harbaugh and others, 2000), but not with earlier versions of MODFLOW.

Input Data Instructions

The modification of SFR2 to simulate unsaturated flow relies on the specific yield values as specified in the Layer Property 

Flow (LPF) Package, the Hydrogeologic-Unit Flow (HUF) Package, or the Block-Centered Flow (BCF) Package. When the 

option to use vertical hydraulic conductivity in the LPF Package is specified, the layer(s) that contain cells where unsaturated 

flow will be simulated must be specified as convertible. That is, the variable LAYTYP specified in LPF (or variable LTHUF in 

HUF) must not be equal to zero, otherwise the model will print an error and stop execution. 

Additional variables that must be specified to define hydraulic properties of the unsaturated zone are all included within the 

SFR2 input file. All values are entered in as free format. Data input for SFR1 works without modification if unsaturated flow is 

not simulated. Parameters can be used to define streambed hydraulic conductivity only when data input follows the SFR1 input 

structure (Prudic and others, 2004); parameters cannot be specified when simulating unsaturated flow beneath streams. The 

calculation of sensitivities for, or estimation of, parameters using the Sensitivity Process is supported by neither SFR1 nor SFR2. 

Additionally, the Ground-Water Transport Process is only available using the original SFR1 input structure and is not available 

when simulating unsaturated flow beneath streams.

For Each Simulation:

0. Data: [#Text]

Text A character variable (up to 199 characters) that is printed when the file is read. The “#” character must be in 

column 1, and, accordingly, the variable starts in column 2. Any characters can be included in Text.

Note 1: Item 0 can be repeated multiple times.

1. Data: NSTRM NSS NSFRPAR NPARSEG CONST DLEAK ISTCB1 ISTCB2 {ISFROPT} {NSTRAIL} 

{ISUZN} {NSFRSETS} {IRTFLG} {NUMTIM} {WEIGHT}{FLWTOL}

NSTRM An integer value that can be specified to be positive or negative. The absolute value of NSTRM is equal to 

the number of stream reaches (finite-difference cells) that are active during the simulation and the number of 

lines of data to be included in Item 2, described below. When NSTRM is specified to be a negative integer, it is 

also used as a flag for changing the format of the data input, for simulating unsaturated flow beneath streams, 

and (or) for simulating transient streamflow routing (for MODFLOW-2005 simulations only), depending on 

the values specified for variables ISFROPT and IRTFLG, as described below. When NSTRM is negative, 

NSFRPAR must be set to zero, which means that parameters cannot be specified.

NSS An integer value equal to the number of stream segments (consisting of one or more reaches) that are used 

to define the complete stream network. The value of NSS represents the number of segments that must be 

defined through a combination of parameters and variables in Item 4 or variables in Item 6.

NSFRPAR An integer value equal to the number of stream parameters (associated with one or more segments) to be 

defined. This variable must be zero when NSTRM is negative.

18 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
NPARSEG An integer value equal to (or exceeding) the number of stream-segment definitions associated with all 

parameters. This number can be more than the total number of segments (NSS) in the stream network because 

the same segment can be defined in multiple parameters, and because parameters can be time-varying. 

NPARSEG must equal or exceed the sum of NLST × N for all parameters, where N is the greater of 1 and 

NUMINST; that is, NPARSEG must equal or exceed the total number of repetitions of item 4b. This variable 

must be zero when NSTRM is negative.

CONST A real value (or conversion factor) used in calculating stream depth for stream reach. If stream depth is not 

calculated using Manning’s equation for any stream segment (that is, ICALC does not equal 1 or 2), then 

a value of zero can be entered. If Manning’s equation is used, a constant of 1.486 is used for flow units of 

cubic feet per second, and a constant of 1.0 is used for units of cubic meters per second. The constant must 

be multiplied by 86,400 when using time units of days in the simulation. An explanation of time units used in 

MODFLOW is given by Harbaugh and others (2000, p. 10).

DLEAK A real value equal to the tolerance level of stream depth used in computing leakage between each stream 

reach and active model cell. Value is in units of length. Usually a value of 0.0001 is sufficient when units of 

feet or meters are used in model.

ISTCB1 An integer value used as a flag for writing stream-aquifer leakage values. If ISTCB1 > 0, it is the unit 

number to which unformatted leakage between each stream reach and corresponding model cell will be 

saved to a file whenever the cell-by-cell budget has been specified in Output Control (see Harbaugh and 

others, 2000, pages 52-55). If ISTCB1 = 0, leakage values will not be printed or saved. If ISTCB1 < 0, all 

information on inflows and outflows from each reach; on stream depth, width, and streambed conductance; 

and on head difference and gradient across the streambed will be printed in the main listing file whenever a 

cell-by-cell budget has been specified in Output Control.

ISTCB2 An integer value used as a flag for writing to a separate formatted file all information on inflows and outflows 

from each reach; on stream depth, width, and streambed conductance; and on head difference and gradient 

across the streambed. If ISTCB2 > 0, then ISTCB2 also represents the unit number to which all information 

for each stream reach will be saved to a separate file when a cell-by-cell budget has been specified in Output 

Control. If ISTCB2 < 0, it is the unit number to which unformatted streamflow out of each reach will be 

saved to a file whenever the cell-by-cell budget has been specified in Output Control.

If NSTRM less than 0:

ISFROPT An integer value that defines the format of the input data and whether or not unsaturated flow is simulated 

 beneath streams. Values of ISFROPT are defined as follows:

 0 No vertical unsaturated flow beneath streams. Streambed elevations, stream slope, streambed thickness, and

 streambed hydraulic conductivity are read for each stress period using variables defined in Items 6b and 6c; 

 the optional variables in Item 2 are not used. 

 1 No vertical unsaturated flow beneath streams. Streambed elevation, stream slope, streambed thickness, and 

 streambed hydraulic conductivity are read for each reach only once at the beginning of the simulation using 

 optional variables defined in Item 2; Items 6b and 6c are used to define stream width and depth for ICALC =

 0 and stream width for ICALC = 1.

 2 Streambed and unsaturated-zone properties are read for each reach only once at the beginning of the

 simulation using optional variables defined in Item 2; Items 6b and 6c are used to define stream width and 

 depth for ICALC = 0 and stream width for ICALC = 1. When using the LPF Package, saturated vertical

 hydraulic conductivity for the unsaturated zone is the same as the vertical hydraulic conductivity of the 

 corresponding layer in LPF and input variable UHC is not read.

 3 Same as 2 except saturated vertical hydraulic conductivity for the unsaturated zone (input variable UHC) is

 read for each reach. 

 4 Streambed and unsaturated-zone properties are read for the beginning and end of each stream segment using

 variables defined in Items 6b and 6c; the optional variables in Item 2 are not used. Streambed properties can 

 vary each stress period. When using the LPF Package, saturated vertical hydraulic conductivity for the 

 unsaturated zone is the same as the vertical hydraulic conductivity of the corresponding layer in LPF 

 and input variable UHC1 is not read.

Appendix 1 19
 5 

20 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
One record for each stream reach:

2. Data: KRCH IRCH JRCH ISEG IREACH RCHLEN {STRTOP} {SLOPE} {STRTHICK} {STRHC1} 

 {THTS} {THTI} {EPS} {UHC}

KRCH An integer value equal to the layer number of the cell containing the stream reach.

IRCH An integer value equal to the row number of the cell containing the stream reach.

JRCH An integer value equal to the column number of the cell containing the stream reach.

ISEG An integer value equal to the number of stream segment in which this reach is located. Stream segments 

contain one or more reaches and are assumed to have uniform or linearly varying characteristics unless 

NSTRM < 0 and ISFROPT = 1, 2, or 3 (in which case streambed information is specified individually for 

each reach).

IREACH An integer value equal to the sequential number in a stream segment of this reach (where a reach corresponds 

to a single cell in the model). Numbering of reaches in a segment begins with 1 for the farthest upstream 

reach and continues in downstream order to the last reach of the segment. 

RCHLEN A real number equal to the length of channel of the stream reach within this model cell. The length of a 

stream reach can exceed the model cell dimensions because of the meandering nature of many streams. 

The length is used to calculate the streambed conductance for this reach. Also, the sum of the lengths of all 

stream reaches within a segment is used to calculate the average slope of the channel for the segment and 

subsequently other values, such as the elevation of the streambed and stream stage.

STRTOP A real number equal to the top elevation of the streambed. This variable is read when ISFROPT is 1, 2, or 3.

SLOPE A real number equal to the stream slope across the reach. This variable is read when ISFROPT is 1, 2, or 3. 

STRTHICK A real number equal to the thickness of the streambed. This variable is read when ISFROPT is 1, 2, or 3.

STRHC1 A real number equal to the hydraulic conductivity of the streambed. This variable is read when ISFROPT is 

1, 2, or 3.

THTS A real number equal to the saturated volumetric water content in the unsaturated zone. This variable is read 

when ISFROPT is 2 or 3. 

THTI A real number equal to the initial volumetric water content. THTI must be less than or equal to THTS and 

greater than or equal to THTS minus the specific yield defined in either LPF, BCF, or HUF. This variable is 

read when ISFROPT is 2 or 3. 

EPS A real number equal to the Brooks-Corey exponent used in the relation between water content and hydraulic 

conductivity within the unsaturated zone (Brooks and Corey, 1966). This variable is read when ISFROPT is 

2 or 3. 

UHC A real number equal to the vertical saturated hydraulic conductivity of the unsaturated zone. This variable 

is necessary when using BCF or HUF, whereas it is optional when using LPF. This variable is read when 

ISFROPT is 3. 

Note 5: Reach information is read in sequential order from upstream to downstream, first by segments, and then 

sequentially by reaches. If segments are not numbered sequentially in downstream order, then the inflow to a segment during 

the current MODFLOW iteration will be the outflow from an upstream segment calculated during the previous iteration. 

Lagging the inflow by one iteration does not change the solution for flows into and out of a segment after MODFLOW 

converges; however, this approach may require an additional iteration for the model to converge. Reaches must be listed and read 

sequentially because the order determines the connections of inflows and outflows within a stream segment. 

Note 6: The stream network is assumed to remain fixed geometrically over the duration of a simulation. The active part 

of the stream network, however, can be made to vary over time by making selected stream segments inactive for selected stress 

periods. This would be implemented by setting the streambed hydraulic conductivity, segment inflow, overland runoff, and direct 

precipitation to zero for the inactive segments in Items 4 or 6 for the specific stress periods when they are known to be inactive 

or dry.

Appendix 1 21
Note 7: If the model cell corresponding to a stream reach is inactive, the program will search for the uppermost active cell 

in the vertical column to apply the stream leakage. If there are no active cells or if the cell is a constant head, no interaction is 

allowed and flow in the reach is passed to the next reach. 

Note 8: When STRTOP, SLOPE, STRTHICK, and STRHC1 are specified for each reach, then HCOND1, THICKM1, 

ELEVUP, HCOND2, THICKM2, and ELEVDN are not read using Items 6b or 6c. 

Note 9: The residual water content for each cell is not specified by the user because it is calculated based on the specified 

saturated water content minus the specific yield of the active model cell corresponding to the stream reach. The calculation is 

performed internally to assure continuity between unsaturated and saturated zone storage. 

Note 10: Although unsaturated flow variables THTS, THTI, EPS, and UHC will not be used for reaches that are designated 

as ICALC = 0, 3, and 4 in Items 4b or 6a, values for these variables must be included for all reaches when ISFROPT = 2 or 3. 

Dummy values for these variables may be used for reaches that are designated as ICALC = 0, 3, and 4.

When NSFRPAR = 0, (including when NSTRM is negative), skip Items 3 and 4 and enter all stream segment data using Items 5 

and 6.

When NSFRPAR > 0, Items 3 and 4 are repeated NSFRPAR times:

3. Data: PARNAM PARTYP Parval NLST [INSTANCES NUMINST] 

PARNAM A set of characters used to name a parameter to be defined. This name can consist of 1 to 10 characters and 

is not case sensitive. That is, any combination of the same characters with different cases will be equivalent.

PARTYP Type of parameter to be defined. For the SFR2 Package, the only allowed parameter type is SFR, which 

defines values of the streambed hydraulic conductivity.

Parval A real number that is a parameter value that may be overridden by a value in the Sensitivity Process input 

file when ISENALL in that file is less than zero.

NLST An integer value that is the number of stream segments associated with this parameter (this value also 

indicates how many times Item 4 is repeated in the next block of input data).

INSTANCES An optional keyword that designates a parameter as time varying. The keyword is case-insensitive; that is, 

it may be entered in any combination of upper- and lower-case letters. If INSTANCES is present, it must be 

followed by a value for NUMINST. If INSTANCES is absent, PARNAM is non-time-varying and NUMINST

should not be present.

NUMINST An integer value that is the number of instances that are included in the definition of a time-varying 

parameter, where each instance is a sequence of Item 4 (Parts 4b through 4g) defining reaches and associated 

properties. If the keyword INSTANCES is present, NUMINST must be present and must be at least 1. If the 

keyword INSTANCES is absent, NUMINST should not be present.

When PARNAM is not time-varying, Part 4a is not read, and Parts 4b through 4g are read sequentially for each of NLST stream 

segments; and 

When PARNAM is time-varying, NUMINST instances are read. For each instance, Part 4a is read, and then Parts 4b through 4g 

are read sequentially for each of NLST stream segments (see notes 11 and 12): 

22 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
4a. Data: {INSTNAM}

INSTNAM The name of an instance associated with the parameter PARNAM specified in the corresponding Item 3 

(INSTNAM is read only if PARNAM is time-varying). The name can consist of 1 to 10 characters and is not 

case sensitive. That is, any combination of the same characters with different case will be equivalent. Names 

entered for INSTNAM must be unique for any given parameter, but names may be reused for instances 

associated with different parameters.

4b. Data: NSEG ICALC OUTSEG IUPSEG {IPRIOR} {NSTRPTS} FLOW RUNOFF ETSW PPTSW 

{ROUGHCH} {ROUGHBK} {CDPTH} {FDPTH} {AWDTH} {BWDTH}

NSEG An integer value of the stream segment for which information is given to identify inflow, outflow, and 

computation of stream depth.

ICALC An integer value used to indicate method used to calculate stream depth in this segment. 

	 • If ICALC ≤ 0, stream depth in each reach is specified at the beginning of a stress period and remains 

unchanged unless flow at the midpoint of a reach is zero, then depth is set to zero in that reach. 

	 • If ICALC > 0, stream depth is calculated and updated each iteration of the MODFLOW solver within a time 

step. 

	 • If ICALC = 1, stream depth is calculated using Manning’s equation and assuming a wide rectangular 

channel.

	 • If ICALC = 2, stream depth is calculated using Manning’s equation and assuming an eight-point channel 

cross section for each segment (which allows for the computation of a wetted perimeter and for changing 

hydraulic conductance of the streambed in relation to changes in flow). 

	 • If ICALC = 3, stream depth and width are calculated using a power function relating each to streamflow 

(Q) using equations 8 and 9 where DEPTH(y) = CDPTH × QFDPTH and WIDTH(w) = AWDTH × QBWDTH. 

	 • If ICALC = 4, stream depth and width are calculated using a table relating streamflow to depth and width 

(the table is defined in Part 4f). 

OUTSEG An integer value of the downstream stream segment that receives tributary inflow from the last downstream 

reach of this segment. If this segment (identified by NSEG) does not feed (or discharge into) another 

downstream (tributary) segment, then enter a value of “0” for this variable (Segment 7 in figure 3 is an 

example for which OUTSEG = 0). If the segment ends within the modeled grid and OUTSEG = 0, outflow 

from the segment is not routed anywhere and is no longer part of the stream network. One may wish to use 

this if all flow in the stream gets diverted into a lined canal or into a pipe. If the flow out of this segment 

discharges into a lake, set OUTSEG equal to the negative value of the lake identification number (where the 

minus sign is used as a flag to tell the model that flow enters a lake rather than a tributary stream segment). 

IUPSEG An integer value of the upstream segment from which water is diverted (or withdrawn) to supply inflow to 

this stream segment if this segment originates as a diversion from an upstream segment. If the source of a 

stream segment is discharge from a lake, set IUPSEG equal to the negative value of the lake identification 

number (where the minus sign is used as a flag to tell the model that streamflow into this segment is derived 

from a lake rather than a stream segment). If this stream segment (identified by NSEG) does not receive 

inflow as a diversion from an upstream segment, then set IUPSEG = 0.

IPRIOR An integer value that only is specified if IUPSEG > 0 (do not specify a value in this field if IUPSEG = 0 

or IUPSEG < 0). IPRIOR defines the prioritization system for diversion, such as when insufficient water 

is available to meet all diversion stipulations, and is used in conjunction with the value of FLOW (specified 

below). 

	 • When IPRIOR = 0, then if the specified diversion flow (FLOW) is greater than the flow available in the 

stream segment from which the diversion is made, the diversion is reduced to the amount available, which 

will leave no flow available for tributary flow into a downstream tributary of segment IUPSEG. 

Appendix 1 23
	 • When IPRIOR = -1, then if the specified diversion flow (FLOW) is greater than the flow available in the 

stream segment from which the diversion is made, no water is diverted from the stream. This approach 

assumes that once flow in the stream is sufficiently low, diversions from the stream cease, and is the 

“priority” algorithm that originally was programmed into the STR1 Package (Prudic, 1989). 

	 • When IPRIOR = -2, then the amount of the diversion is computed as a fraction of the available flow in 

segment IUPSEG; in this case, 0.0 ≤ FLOW ≤ 1.0.

	 • When IPRIOR = -3, then a diversion is made only if the streamflow leaving segment IUPSEG exceeds 

the value of FLOW. If this occurs, then the quantity of water diverted is the excess flow and the quantity that 

flows from the last reach of segment IUPSEG into its downstream tributary (OUTSEG) is equal to FLOW. 

This represents a flood-control type of diversion, as described by Danskin and Hanson (2002).

NSTRPTS An integer value specified only when ICALC = 4. It is used to dimension a table relating streamflow with 

stream depth and width as specified in Items 4e and 6e. NSTRPTS must be at least 2 but not more than 50. 

If the table exceeds 3 × 50 (for streamflow, stream depth, and width) values, then MAXPTS in the allocation 

subroutine GWF1SFR1ALP will need to be increased from 3 × 50 to 3 × (the desired maximum value). 

FLOW A real number that is the streamflow (in units of volume per time) entering or leaving the upstream end of a 

stream segment (that is, into the first reach). 

	 • If the stream is a headwater stream, FLOW defines the total inflow to the first reach of the segment. The 

value can be any number ≥ 0.

	 • If the stream is a tributary stream, FLOW defines additional specified inflow to or withdrawal from the 

first reach of the segment (that is, in addition to the discharge from the upstream segment of which this is 

a tributary). This additional flow does not interact with the ground-water system. For example, a positive 

number might be used to represent direct outflow into a stream from a sewage treatment plant, whereas 

a negative number might be used to represent pumpage directly from a stream into an intake pipe for a 

municipal water treatment plant. (Also see additional explanatory notes below.)

	 • If the stream is a diversionary stream, and the diversion is from another stream segment, FLOW defines the 

streamflow diverted from the last reach of stream segment IUPSEG into the first reach of this segment. The 

diversion is computed or adjusted according to the value of IPRIOR. 

	 • If the stream is a diversionary stream, and the diversion is from a lake, FLOW defines a fixed rate of 

discharge diverted from the lake into the first reach of this stream segment (unless the lake goes dry) and flow 

from the lake is not dependent on the value of ICALC. However, if FLOW = 0, then the lake outflow into the 

first reach of this segment will be calculated on the basis of lake stage relative to the top of the streambed for 

the first reach using one of the methods defined by ICALC. 

RUNOFF A real number that is the volumetric rate of the diffuse overland runoff that enters the stream segment (in 

units of volume per time). The specified rate is apportioned to each reach of the segment in direct relation to 

the fraction of the total length of the stream channel in the segment that is present in each reach. 

ETSW A real number that is the volumetric rate per unit area of water removed by evapotranspiration directly from 

the stream channel (in units of length per time). ETSW is defined as a positive value.

PPTSW A real number that is the volumetric rate per unit area of water added by precipitation directly on the stream 

channel (in units of length per time).

ROUGHCH A real number that is Manning’s roughness coefficient for the channel in all reaches in this segment. This 

variable is only specified if ICALC = 1 or 2.

ROUGHBK A real number that is Manning’s roughness coefficient for the overbank areas in all reaches in this segment. 

This variable is only specified if ICALC = 2.

24 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
CDPTH A real number that is the coefficient used in the equation: 

DEPTH = CDPTH × QFDPTH that relates stream depth in all reaches in this segment to streamflow. This variable 

is only specified if ICALC = 3.

FDPTH A real number that is the coefficient used in the equation:

DEPTH = CDPTH × QFDPTH that relates stream depth in all reaches in this segment to streamflow. This variable 

is only specified if ICALC = 3.

AWDTH A real number that is the coefficient used in the equation:

WIDTH = AWDTH × QBWDTH that relates stream width in all reaches in this segment to streamflow. This 

variable is only specified if ICALC = 3.

BWDTH A real number that is the coefficient used in the equation:

WIDTH = AWDTH × QBWDTH that relates stream width in all reaches in this segment to streamflow. This 

variable is only specified if ICALC = 3.

4c. Data: Hc1fact THICKM1 ELEVUP {WIDTH1} {DEPTH1} 

Hc1fact A real number that is a factor used to calculate hydraulic conductivity of the streambed at the upstream end of 

this segment from the parameter value (in units of length per time).

THICKM1 A real number that is the thickness of streambed material at the upstream end of this segment (in units of 

length). 

ELEVUP A real number that is the elevation of the top of the streambed at the upstream end of this segment (in units of 

length). 

WIDTH1 A real number that is the average width of the stream channel at the upstream end of this segment (in units of 

length). This variable is only specified if ICALC < 1. 

DEPTH1 A real number that is the average depth of water in the channel at the upstream end of this segment (units of 

length). This variable is only specified if ICALC = 0, in which case the stream stage in a reach is assumed to 

equal the elevation of the top of the streambed plus the depth of water. 

4d. Data: Hc2fact THICKM2 ELEVDN {WIDTH2} {DEPTH2}

Hc2fact A real number that is the factor used to calculate hydraulic conductivity of the streambed at the downstream 

end of this segment from the parameter value (units of length per time).

THICKM2 A real number that is the thickness of streambed material at the downstream end of this segment (units of 

length). 

ELEVDN A real number that is the elevation of the top of the streambed at the downstream end of this segment (units of 

length). 

WIDTH2 A real number that is the average width of the stream channel at the downstream end of this segment (units of 

length). This variable is only specified if ICALC < 1. 

DEPTH2 A real number that is the average depth of water in the channel at the downstream end of this segment (units 

of length). This variable is only specified if ICALC = 0, in which case the stream stage in a reach is assumed 

to equal the elevation of the top of the streambed plus the depth of water. 

Appendix 1 25
If ICALC = 2:

4e. Data: XCPT1 XCPT2 ... XCPT8 

 Data: ZCPT1 ZCPT2 ... ZCPT8

XCPTi A real number that is the distance relative to the left bank of the stream channel (when looking downstream) 

for the eight points (XCPT1 through XCPT8) used to describe the geometry of this segment of the stream 

channel. By definition, location XCPT1 represents the left edge of the channel cross section, and its value 

should be set equal to 0.0; values XCPT2 through XCPT8 should be equal to or greater than the previous 

distance.

ZCPTi A real number that is the height relative to the top of the lowest elevation of the streambed (thalweg). One 

value (ZCPT1 through ZCPT8) is needed for each of the eight horizontal distances defined by XCPTi. The 

location of the thalweg (set equal to 0.0) can be any location from XCPT2 through XCPT7. 

If ICALC = 4:

4f. Data: FLOWTAB(1) FLOWTAB(2) ... FLOWTAB(NSTRPTS)

Data: DPTHTAB(1) DPTHTAB(2) ... DPTHTAB(NSTRPTS)

Data: WDTHTAB(1) WDTHTAB(2) ... WDTHTAB(NSTRPTS) 

FLOWTAB A real number that is the streamflow (units of volume per time) related to a given depth and width. One value 

is needed for each streamflow that has a corresponding value of depth and width up to the total number of 

values used to define the table—FLOWTAB(1) through FLOWTAB(NSTRPTS). NSTRPTS is defined in Part 

4b.

DPTHTAB A real number that is the average depth (units of length) corresponding to a given flow. The number and 

order of values, DPTHTAB(1)through DPTHTAB(NSTRPTS)must coincide with the streamflow values 

FLOWTAB(1)through FLOWTAB(NSTRPTS).

WDTHTAB A real number that is the stream width (units of length) corresponding to a given flow. The number and 

order of values, WDTHTAB(1)through WDTHTAB(NSTRPTS), must coincide with the streamflow values 

FLOWTAB(1)through FLOWTAB(NSTRPTS).

If Ground-Water Transport (GWT) Process is active (Ftype “GWT” exists), then read Part 4g NSOL times (where NSOL is the 

number of solutes being simulated, see note 17):

4g. Data: {CONCQ(NSOL)} CONCRUN(NSOL) CONCPPT(NSOL)

CONCQ A real number that is the solute concentration associated with “FLOW” (Part 4b—specified flow entering the 

upstream end of the stream segment). Only specify this variable for stream segments where IUPSEG = 0.

CONCRUN A real number that is the solute concentration in the overland runoff entering the stream segment. 

CONCPPT A real number that is the solute concentration in precipitation that directly falls onto the stream surface.

NSOL An integer value that is the number of solutes that are being simulated. The version of the Ground-Water 

Flow (GWT) Process for MODFLOW-2000 (MODFLOW-GWT) documented by Konikow and others (1996) 

is limited to simulating one solute, and until a new version is released, NSOL must be 1.

Note 11: Item 4, Parts 4b through 4g must be completed sequentially for each of the NLST stream segments; that is, Parts 

4b through 4g for one stream segment must be entered before Parts 4b through 4g of the next stream segment. However, the data 

for stream segments need not be entered sequentially by stream segment number. For example, data for stream segment 2 can be 

entered before data for stream segment 1. 

26 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Note 12: When the defined parameter (PARNAM) is not time-varying, Item 4, Part 4a is omitted, and each stream segment 

controlled by parameter PARNAM is defined by one sequence of Parts 4b through 4g. When the defined parameter (PARNAM) 

is time-varying, Item 4, Parts 4a through 4g must be completed for each of NUMINST instances. For each instance, Part 4a 

is defined followed by Parts 4b through 4g for each of the NLST stream segments associated with that instance. For example, 

NUMINST = 2, NLST = 3, ICALC = 1 and the GWT Process is inactive, the sequence of records for Item 4 would be: 4a, 4b, 4c, 

4d, 4b, 4c, 4d, 4b, 4c, 4d, 4a, 4b, 4c, 4d, 4b, 4c, 4d, 4b, 4c, and 4d. 

Note 13: Record 4b will contain 8 to 13 variables; depending on the values of ICALC and IUPSEG. ICALC determines 

how stream depth is to be calculated; when ICALC is 1 or 2, depth is calculated using Manning’s equation, which, in turn, 

requires a channel roughness coefficient (ICALC = 1) or a channel and bank roughness coefficient (ICALC = 2). Similarly, 

Parts 4c and 4d will include 3 to 5 values. 

Note 14: A stream segment that receives inflow from upstream segments is allowed to have as many as ten upstream 

segments feeding it, as defined by the respective values of OUTSEG in Part 4b.

Note 15: Stream properties and stresses defined in Parts 4b are assumed constant and uniform within a single stream 

segment. Additionally, hydraulic conductivity, streambed thickness, elevation of top of streambed, stream width, and stream 

depth may vary smoothly and linearly within a single stream segment. For these variables, data values at the upstream end of the 

segment are described in Part 4c and data values at the downstream end of the segment are described in Part 4d. Values of these 

variables for individual reaches of a segment are estimated using linear interpolation. To make any variable the same throughout 

the segment, simply specify equal values in Parts 4c and 4d. The two elevations in Parts 4c and 4d are used in conjunction 

with the total length of the stream segment (calculated from RCHLEN given for each reach in Item 2) to compute the slope of 

the stream and the elevations for any intermediate reaches. The streambed thickness is subtracted from the top of streambed 

elevations to calculate the elevations of the bottom of the streambed (used in calculations of leakage).

Note 16: If Part 4e is included (for ICALC = 2), it is assumed that the cross-sectional geometry defined by these data is the 

same over the entire length of the segment. Similarly, if Part 4f is included (for ICALC = 4), it is assumed the tabulated relation 

between streamflow and stream depth and width is the same over the entire length of the segment. 

Note 17: If the Lake (LAK3) Package (Merritt and Konikow, 2000) is also implemented, then flow out of the lake into

a stream segment is dependent on the option used to compute stream depth (ICALC = 1, 2, 3, or 4). Constant discharge from a 

lake can be simulated no matter what value of ICALC is assigned to the stream segment emanating from the lake by assigning a 

positive value to FLOW in Part 4b.

Note 18: If a diversionary flow is large enough to warrant representation in the model, but is discharged into a pipeline, 

lined canal, or other structure or system that does not interact with the aquifer and the flow might exceed the available 

streamflow, then there is an alternative means to represent it. Instead of specifying a negative value of FLOW, we suggest 

representing the withdrawal by a single-reach diversionary stream segment, which would be located in the same model cell 

as the reach from the upstream segment (IUPSEG) from which the diversion is made; specifying the segment’s streambed 

hydraulic conductivity equal to 0 will preclude interaction with the aquifer and setting OUTSEG = 0 will remove the flow from 

the system. The diversion will then be subject to the constraints associated with the value of IPRIOR.

Note 19: The version of MODFLOW-GWT documented by Konikow and others (1996) is limited to simulating one solute, 

and until a new multi-species version is released, NSOL automatically defaults to a value of 1. 

Appendix 1 27
For Each Stress Period:

5. Data: ITMP IRDFLG IPTFLG {NP}

ITMP An integer value for reusing or reading stream segment data that can change each stress period. If ITMP = 0 

then all stream segment data are defined by Item 4 (NSFRPAR > 0; number of stream parameters is greater 

than 0). If ITMP > 0, then stream segment data are not defined in Item 4 and must be defined in Item 6 below 

for a number of segments equal to the value of ITMP. If ITMP < 0, then stream segment data not defined in 

Item 4 will be reused from the last stress period (Item 6 is not read for the current stress period). ITMP must 

be defined ≥ 0 for the first stress period of a simulation.

IRDFLG An integer value for printing input data specified for this stress period. If IRDFLG = 0, input data for this 

stress period will be printed. If IRDFLG > 0, then input data for this stress period will not be printed.

IPTFLG An integer value for printing streamflow-routing results during this stress period. If IPTFLG = 0, or 

whenever the variable ICBCFL or “Save Budget” is specified in Output Control, the results for specified time 

steps during this stress period will be printed. If IPTFLG > 0, then the results during this stress period will 

not be printed.

NP An integer value of the number of parameters used in the current stress period. The parameters being used are 

subsequently listed in Item 7 below. NP and Item 7 below are not read when NSFRPAR = 0. 

Note 20: In each stress period, the sum of ITMP plus the sum of all NLST values in Item 3 associated with the NP

parameters listed in Item 7 must equal or be less than (some stream segments may not be active during a stress period) the total 

number of stream segments in the stream network (NSS of Item 1). Stream segments defined by Items 3, 4, and 7 cannot be 

repeated using ITMP and Item 6.

If ITMP > 0:

6a. Data: NSEG ICALC OUTSEG IUPSEG {IPRIOR} {NSTRPTS} FLOW RUNOFF ETSW PPTSW 

{ROUGHCH} {ROUGHBK} {CDPTH} {FDPTH} {AWDTH} {BWDTH}

NSEG An integer value of the stream segment for which information is given to identify inflow, outflow, and 

computation of stream depth.

ICALC An integer value used to indicate method used to calculate stream depth in this segment. 

	 • If ICALC ≤ 0, stream depth in each reach is specified at the beginning of a stress period and remains 

unchanged unless flow at the midpoint of a reach is zero, then depth is set to zero in that reach. Unsaturated 

flow is not allowed.

	 • If ICALC > 0, stream depth is calculated and updated each iteration of the MODFLOW solver within a time 

step. 

	 • If ICALC = 1, stream depth is calculated using Manning’s equation and assuming a wide rectangular 

channel. Unsaturated flow is simulated when ISFROPT > 1 and the value of ICALC must remain 1 for the 

specified segment during the simulation.

	 • If ICALC = 2, stream depth is calculated using Manning’s equation and assuming an eight-point channel 

cross section for each segment (which allows for the computation of a wetted perimeter and for changing 

hydraulic conductance of the streambed in relation to changes in flow). Unsaturated flow is simulated when 

ISFROPT > 1 and the value of ICALC must remain 2 for the specified segment during the simulation.

	 • If ICALC = 3, stream depth and width are calculated using a power function relating each to streamflow 

(Q) using equations 8 and 9 where DEPTH(y) = CDPTH × QFDPTH and WIDTH(w) = AWDTH × QBWDTH. 

Unsaturated flow is not allowed.

	 • If ICALC = 4, stream depth and width are calculated using a table relating streamflow to depth and width 

(the table is defined in Part 4f). Unsaturated flow is not allowed.

28 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
OUTSEG An integer value of the downstream stream segment that receives tributary inflow from the last downstream 

reach of this segment. If this segment (identified by NSEG) does not feed (or discharge into) another 

downstream (tributary) segment, then enter a value of “0” for this variable (Segment 7 in figure 3 is an 

example for which OUTSEG = 0). If the segment ends within the modeled grid and OUTSEG = 0, outflow 

from the segment is not routed anywhere and is no longer part of the stream network. One may wish to use 

this if all flow in the stream gets diverted into a lined canal or into a pipe. If the flow out of this segment 

discharges into a lake, set OUTSEG equal to the negative value of the lake identification number (where the 

minus sign is used as a flag to tell the model that flow enters a lake rather than a tributary stream segment). 

IUPSEG An integer value of the upstream segment from which water is diverted (or withdrawn) to supply inflow to 

this stream segment if this segment originates as a diversion from an upstream segment. If the source of a 

stream segment is discharge from a lake, set IUPSEG equal to the negative value of the lake identification 

number (where the minus sign is used as a flag to tell the model that streamflow into this segment is derived 

from a lake rather than a stream segment). If this stream segment (identified by NSEG) does not receive 

inflow as a diversion from an upstream segment, then set IUPSEG = 0.

IPRIOR An integer value that only is specified if IUPSEG > 0 (do not specify a value in this field if IUPSEG = 0 

or IUPSEG < 0). IPRIOR defines the prioritization system for diversion, such as when insufficient water 

is available to meet all diversion stipulations, and is used in conjunction with the value of FLOW (specified 

below). 

	 • When IPRIOR = 0, then if the specified diversion flow (FLOW) is greater than the flow available in the 

stream segment from which the diversion is made, the diversion is reduced to the amount available, which 

will leave no flow available for tributary flow into a downstream tributary of segment IUPSEG. 

	 • When IPRIOR = -1, then if the specified diversion flow (FLOW) is greater than the flow available in the 

stream segment from which the diversion is made, no water is diverted from the stream. This approach 

assumes that once flow in the stream is sufficiently low, diversions from the stream cease, and is the 

“priority” algorithm that originally was programmed into the STR1 Package (Prudic, 1989). 

	 • When IPRIOR = -2, then the amount of the diversion is computed as a fraction of the available flow in 

segment IUPSEG; in this case, 0.0 ≤ FLOW ≤ 1.0.

	 • When IPRIOR = -3, then a diversion is made only if the streamflow leaving segment IUPSEG exceeds 

the value of FLOW. If this occurs, then the quantity of water diverted is the excess flow and the quantity that 

flows from the last reach of segment IUPSEG into its downstream tributary (OUTSEG) is equal to FLOW. 

This represents a flood-control type of diversion, as described by Danskin and Hanson (2002).

NSTRPTS An integer value specified only when ICALC = 4. It is used to dimension a table relating streamflow with 

stream depth and width as specified in Items 4e and 6e. NSTRPTS must be at least 2 but not more than 50. 

If the table exceeds 3 × 50 (for streamflow, stream depth, and width) values, then MAXPTS in the allocation 

subroutine GWF1SFR1ALP will need to be increased from 3 × 50 to 3 × (the desired maximum value). 

FLOW A real number that is the streamflow (in units of volume per time) entering or leaving the upstream end of a 

stream segment (that is, into the first reach). 

	 • If the stream is a headwater stream, FLOW defines the total inflow to the first reach of the segment. The 

value can be any number ≥ 0.

	 • If the stream is a tributary stream, FLOW defines additional specified inflow to or withdrawal from the 

first reach of the segment (that is, in addition to the discharge from the upstream segment of which this is 

a tributary). This additional flow does not interact with the ground-water system. For example, a positive 

number might be used to represent direct outflow into a stream from a sewage treatment plant, whereas 

a negative number might be used to represent pumpage directly from a stream into an intake pipe for a 

municipal water treatment plant. (Also see additional explanatory notes below.)

	 • If the stream is a diversionary stream, and the diversion is from another stream segment, FLOW defines the 

streamflow diverted from the last reach of stream segment IUPSEG into the first reach of this segment. The 

diversion is computed or adjusted according to the value of IPRIOR. 

	 • If the stream is a diversionary stream, and the diversion is from a lake, FLOW defines a fixed rate of 

discharge diverted from the lake into the first reach of this stream segment (unless the lake goes dry) and flow 

from the lake is not dependent on the value of ICALC. However, if FLOW = 0, then the lake outflow into the 

first reach of this segment will be calculated on the basis of lake stage relative to the top of the streambed for 

the first reach using one of the methods defined by ICALC. 

Appendix 1 29
RUNOFF A real number that is the volumetric rate of the diffuse overland runoff that enters the stream segment (in 

units of volume per time). The specified rate is apportioned to each reach of the segment in direct relation to 

the fraction of the total length of the stream channel in the segment that is present in each reach. 

ETSW A real number that is the volumetric rate per unit area of water removed by evapotranspiration directly from 

the stream channel (in units of length per time). ETSW is defined as a positive value.

PPTSW A real number that is the volumetric rate per unit area of water added by precipitation directly on the stream 

channel (in units of length per time).

ROUGHCH A real number that is Manning’s roughness coefficient for the channel in all reaches in this segment. This 

variable is only specified if ICALC = 1 or 2.

ROUGHBK A real number that is Manning’s roughness coefficient for the overbank areas in all reaches in this segment. 

This variable is only specified if ICALC = 2.

CDPTH A real number that is the coefficient used in the equation: 

DEPTH = CDPTH × QFDPTH that relates stream depth in all reaches in this segment to streamflow. This variable 

is only specified if ICALC = 3.

FDPTH A real number that is the coefficient used in the equation:

DEPTH = CDPTH × QFDPTH that relates stream depth in all reaches in this segment to streamflow. This variable 

is only specified if ICALC = 3.

AWDTH A real number that is the coefficient used in the equation:

WIDTH = AWDTH × QBWDTH that relates stream width in all reaches in this segment to streamflow. This 

variable is only specified if ICALC = 3.

BWDTH A real number that is the coefficient used in the equation:

WIDTH = AWDTH × QBWDTH that relates stream width in all reaches in this segment to streamflow. This 

variable is only specified if ICALC = 3.

6b. Data: {HCOND1} {THICKM1} {ELEVUP} {WIDTH1} {DEPTH1} 

{THTS1} {THTI1} {EPS1} {UHC1}

HCOND1 Hydraulic conductivity of the streambed at the upstream end of this segment (units of length per time). This 

variable is read for each stress period when NSTRM is positive or when NSTRM is negative and ISFROPT is 

0, 4, or 5.

THICKM1 Thickness of streambed material at the upstream end of this segment (in units of length). This variable is read 

each stress period for all segments when NSTRM is positive or when NSTRM is negative and ISFROPT is 0. 

When ISFROPT is 4 or 5, the variable is read each stress period for a segment when ICALC is 0, 3, or 4, and 

is read only the first stress period when ICALC is 1 or 2.

ELEVUP Elevation of the top of the streambed at the upstream end of this segment (in units of length). This variable is 

read each stress period for all segments when NSTRM is positive or when NSTRM is negative and ISFROPT is 

0. When ISFROPT is 4 or 5, the variable is read each stress period for a segment when ICALC is 0, 3, or 4, 

and is read only the first stress period when ICALC is 1 or 2.

WIDTH1 Average width of the stream channel at the upstream end of this segment (in units of length). This variable is 

read each stress period for all segments identified with an ICALC of 0 and is not dependent on ISFROPT. 

When ICALC is 1, the variable is read each stress period when NSTRM is positive or when NSTRM is negative 

and ISFROPT is 0 or 1, and is read only for the first stress period when ISFROPT is 2, 3, 4, or 5.

DEPTH1 Average depth of water in the channel at the upstream end of this segment (units of length). This variable is 

read only when ICALC is 0 and is not dependent on the value of NSTRM or ISFROPT. The stream stage in a 

reach is assumed to equal the elevation of the top of the streambed plus the depth of water.

THTS1 Saturated volumetric water content in the unsaturated zone beneath the upstream end of this segment. This 

variable is read for the first stress period when ICALC is 1 or 2 and ISFROPT is 4 or 5.

THTI1 Initial volumetric water content beneath the upstream end of this segment. THTI1 must be less than or equal 

to THTS and greater than or equal to the THTS minus the specific yield defined in either LPF, BCF, or HUF. 

This variable is read for the first stress period when ICALC is 1 or 2 and ISFROPT is 4 or 5. 

30 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
EPS1 Brooks-Corey exponent used in the relation between water content and hydraulic conductivity within the 

unsaturated zone beneath the upstream end of this segment. This variable is read for the first stress period 

when ICALC is 1 or 2 and ISFROPT is 4 or 5. 

UHC1 Vertical saturated hydraulic conductivity of the unsaturated zone beneath the upstream end of this segment. 

This variable is necessary when using BCF or HUF, whereas it is optional when using LPF. This variable is 

read only for the first stress period when ICALC is 1 or 2 and ISFROPT is 5. 

6c. Data: {HCOND2} {THICKM2} {ELEVDN} {WIDTH2} {DEPTH2}

{THTS2} {THTI2} {EPS2} {UHC2}

HCOND2 Hydraulic conductivity of the streambed at the downstream end of this segment (units of length per time). 

This variable is read for each stress period when NSTRM is positive or when NSTRM is negative and 

ISFROPT is 0, 4. or 5.

THICKM2 Thickness of streambed material at the downstream end of this segment (in units of length). This variable is 

read each stress period for all segments when NSTRM is positive or when NSTRM is negative and ISFROPT is 

0. When ISFROPT is 4 or 5, the variable is read each stress period for a segment when ICALC is 0, 3, or 4, 

and is read only the first stress period when ICALC is 1 or 2.

ELEVDN Elevation of the top of the streambed at the downstream end of this segment (in units of length). This variable 

is read each stress period for all segments when NSTRM is positive or when NSTRM is negative and ISFROPT

is 0. When ISFROPT is 4 or 5, the variable is read each stress period for a segment when ICALC is 0, 3, or 

4, and is read only the first stress period when ICALC is 1 or 2.

WIDTH2 Average width of the stream channel at the downstream end of this segment (in units of length). This variable 

is read each stress period for all segments identified with an ICALC of 0 and is not dependent on ISFROPT. 

When ICALC is 1, the variable is read each stress period when NSTRM is positive or when NSTRM is negative 

and ISFROPT is 0 or 1, and is read only for the first stress period when ISFROPT is 2, 3, 4, or 5.

DEPTH2 Average depth of water in the channel at the downstream end of this segment (units of length). This variable 

is read only when ICALC is 0 and is not dependent on the value of NSTRM or ISFROPT. The stream stage in 

a reach is assumed to equal the elevation of the top of the streambed plus the depth of water.

THTS2 Saturated volumetric water content in the unsaturated zone beneath the downstream end of this segment. This 

variable is read for the first stress period when ICALC is 1 or 2 and ISFROPT is 4 or 5.

THTI2 Initial volumetric water content beneath the downstream end of this segment. THTI2 must be less than or 

equal to THTS and greater than or equal to THTS minus the specific yield defined in either LPF, BCF. or 

HUF. This variable is read for the first stress period when ICALC is 1 or 2 and ISFROPT is 4 or 5.

EPS2 Brooks-Corey exponent used in the relation between water content and hydraulic conductivity within the 

unsaturated zone beneath the downstream end of this segment. This variable is read for the first stress period 

when ICALC is 1 or 2 and ISFROPT is 4 or 5.

UHC2 Vertical saturated hydraulic conductivity of the unsaturated zone beneath the downstream end of this segment. 

This variable is necessary when using BCF or HUF, whereas it is optional when using LPF. This variable is 

read only for the first stress period when ICALC is 1 or 2 and ISFROPT is 5.

If ICALC = 2:

6d. Data: XCPT1 XCPT2 ... XCPT8 

 Data: ZCPT1 ZCPT2 ... ZCPT8

See Item 4, Part 4e for variable definitions. These variables are read only for the first stress period if ISFROPT is 2, 3, 4, or 5 

(that is, unsaturated flow is active) because the channel geometry cannot change during the simulation when unsaturated flow is 

active. 

Appendix 1 31
If ICALC = 4:

6e. Data: FLOWTAB(1) FLOWTAB(2) ... FLOWTAB(NSTRPTS) 

 Data: DPTHTAB(1) DPTHTAB(2) ... DPTHTAB(NSTRPTS) 

 Data: WDTHTAB(1) WDTHTAB(2) ... WDTHTAB(NSTRPTS) 

See Item 4, Part 4f for variable definitions. 

If Ground-Water Transport (GWT) Process is active (Ftype “GWT” exists), then read Part 6f NSOL times (where NSOL is the 

number of solutes being simulated):

6f. Data: {CONCQ(NSOL)} CONCRUN(NSOL) CONCPPT(NSOL)

See Item 4, Part 4g for variable definitions.

Note 21: Item 6 must be completed ITMP times. The data need not be defined in sequential order by stream segment 

number. All active segments in the stream network must be defined for each stress period through a combination of ITMP and 

Item 6 and NP and Item 7.

Note 22: If ITMP ≤ 0, then Item 6 is excluded for this stress period. If ITMP < 0, then values for Item 6 from the previous 

stress period are reused. If ITMP = 0, then no Item 6 records are read, and all segments must be defined using parameters.

Note 23: If Item 6, Part 6f is included because solute transport is being simulated, then Part 6f should consist of one record 

(line) for each solute; each record must contain two or three values; and there must be as many records as the number of solutes 

being simulated (NSOL). The default value of NSOL is 1 and the presently documented version of MODFLOW-GWT is limited 

to simulating one solute. 

Note 24: All the explanatory notes applicable to Item 4 (except those related to parameters) are also relevant to Item 6. 

If NP > 0, then:

7. Data: Pname [Iname]

Pname The name of a parameter that is being used in the current stress period. Repeat Item 7 NP times (see Item 5).

Iname An instance name that is read only if Pname is a time-varying parameter. Multiple instances of the same 

time-varying parameter are not allowed in a stress period.

Stream Gaging (Monitoring) Station File (GAGE)

Cells of the model grid can be designated as “stream gaging station” locations. At each designated cell (or stream 

reach), the time, stage, streamflow out of that reach, the streambed seepage, unsaturated storage, change in unsaturated storage, 

and ground-water recharge will be written to a separate output file to facilitate model output evaluation and graphical post 

processing of the calculated data. Several options are available to print additional information. The input file for specifying 

gaging station locations is read if the file type (Ftype) “GAGE” is included in the MODFLOW name file. The output file will 

contain two header lines that provide relevant information (the text will be contained within quotes). 

32 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
For Each Simulation, if Gage Package is Used: 

1. Data: NUMGAGE 

NUMGAGE Number of gaging stations.

For Each Gaging Station:

2. Data: GAGESEG GAGERCH UNIT OUTTYPE

GAGESEG An integer value that is the stream segment number where gage is located.

GAGERCH An integer value that is the stream reach number where gage is located.

UNIT An integer value that is the unit number of output file for this gage.

OUTTYPE An integer value that is a flag for type of expanded listing desired in output file:

0 Use standard default listing of time, stage, outflow, and solute concentration.

1 Default values plus depth, width, and flow at midpoint.

2 Default values plus streambed conductance for the reach, head difference across streambed, and 

hydraulic gradient across streambed.

3 Default values plus solute load in stream (if GWT is active).

4 All of the above.

5 Used for diversions to provide a listing of time, stage, flow diverted, maximum assigned diversion rate, 

flow at end of upstream segment prior to diversion, solute concentration, and solute load.

6 Used for unsaturated flow routing to provide a listing of time, stage, ground-water head, streambed 

seepage, change in unsaturated storage, and recharge. 

7 Used for unsaturated flow routing to provide a listing of time and the unsaturated water content profile 

beneath the stream. Two profiles are printed. The first is the volume average water content of ISUZN

unsaturated cells beneath a stream reach. The second profile is the volume averaged water content 

beneath the low flow channel only. If ISUZN =1 then both water content profiles will be the same. 

Note 1: A unique unit number must be specified for each stream gaging station and the MODFLOW name file (see Harbaugh and 

others, 2000, p. 42-44) must specify the file for each unit. The file type in the name file must be “DATA”. 

Note 2: If the LAK3 Package is also active, a gaging station may be placed on a lake. In this case, only two required variables and 

one optional variable are read. The first variable should be the negative value of the lake number and the second is the unit number for 

the output. The original description of the Gage Package (Merritt and Konikow, 2000, p. 57) did not include additional print options 

listed for OUTTYPE. For a gaging station on a lake, a minus sign on the unit number is a flag indicating that OUTTYPE for a lake will 

be read. 

Note 3: Data Set 2 must include exactly NUMGAGE lines (records) of data. If NUMGAGE > 1, it is permissible to interleaf the Item 

2 lines for stream gaging stations with lines for lake gages. Data lines (records) within Item 2 can be listed in any arbitrary order. 

Note 4: Solute is not routed through the unsaturated zone. For this reason, solute routing between streams and ground water 

cannot be simulated when using the unsaturated flow option in SFR2. 

Note 5: Unsaturated flow is not simulated beneath lakes.

Appendix 2 33
Appendix 2: Selected Input Data and Printed Results for the Hypothetical Test 

Simulation. 

The hypothetical test simulation is designed to demonstrate the capability of version 2 of the Streamflow-Routing (SFR2) 

Package that includes unsaturated flow. The test simulation simulates recharge occurring from stream seepage in an alluvial 

basin that is separated from ground water by an unsaturated zone. The test simulation includes a sensitivity analysis of the 

streambed thickness on the timing and magnitude of recharge. A description of the basin-fill aquifer is presented in figure 6. 

Variables used to route flow in streams and flow through the aquifer and results from these simulations are presented in the 

section titled “Hypothetical Test Simulation” and in table 4 and figure 7. 

Selected sections of several key input and output data files are shown below—sometimes with annotations; gaps in the 

listings are indicated by an ellipsis. A complete set of the files is available for distribution over the Internet. Contents of some 

files are enclosed in a border and explanations are noted outside of the border; for other files, explanations are sometimes 

included as comments following a semi-colon on the line being explained. Information pertaining specifically to the SFR2 

Package is highlighted by shading.

34 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Listing of Selected Input Data Sets for the Hypothetical Test Simulation

Following (enclosed in a border) are the contents of the MODFLOW name file for the hypothetical test 

simulation; explanations are noted outside of border:

File name: Testsfr2.nam

LIST 12 testsfr2.lst → Output file for MODFLOW

BAS6 41 testsfr2.bs6 → Input file for Basic Package

LPF 42 testsfr2.lpf → Input file for Layer Property Flow Package

DIS 13 testsfr2.dis → Input file for Discretization Package

pcg 43 testsfr2.pcg → Input file for pcg solver Package

OC 44 testsfr2.oc → Input file for Output Control option 

SFR 45 testsfr2.sfr → Input file for Streamflow-Routing (SFR2) Package

GAGE 46 testsfr2.gag → Input file for Gage Package

global 11 testsfr2.glo → Output file for global output

DATA 15 testsfr2.ghb → Input file for general head boundary

DATA 83 testsfr2.gg1 → Seepage and recharge for stream reach 20

DATA 84 testsfr2.gg2 → Water content profiles for reach 20

DATA 85 testsfr2.gg3 → Seepage and recharge for stream reach 40

DATA 86 testsfr2.gg4 → Water content profiles for reach 40

DATA 87 testsfr2.gg5 → Seepage and recharge for stream reach 60

DATA 88 testsfr2.gg6 → Water content profiles for reach 60

DATA 89 testsfr2.gg7 → Seepage and recharge for stream reach 80

DATA 90 testsfr2.gg8 → Water content profiles for reach 80

 ↑ ↑ ↑

 1 2 3

1 Ftype (that is, the type of file)

2 Unit number

3 File name (name chosen to reflect contents of file)

Appendix 2 35
Following (enclosed in a border) are the contents of the Basic Package input file for the hypothetical test 

simulation; explanations are noted outside of border:

File name: testsfr2.bs6

# BAS6 Package input file for hypothetical test simulation ← 1

# Test of SFR2 Package ← 1

FREE Item 1: options ← 2

INTERNAL 1 (30I2) 1 Item 2: IBOUND, layer 1 ← 3

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

 1 1 1 1 1 1 1 1 1 1 

← 3

 -6999.00 HNOFLO ← 4

 CONSTANT 104.0 ← 5

1 Two header lines of comments. 

2 Options line.

3 IBOUND array.

4 Head value assigned to inactive cells.

5 Initial head information for transient simulation values are rounded to nearest tenth of meter.

36 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Following (enclosed in a border) are the contents of the Discretization Package input file for the hypothetical 

test simulation; explanations are noted outside of border:

File name: testsfr2.dis

# Discretization file input for hypothetical test simulation

# Units: seconds, meters ← 1

 1 7 100 50 1 2 Item 1: NLAY NROW NCOL NPER ITMUNI LENUNI ← 2

 1 Item 2: LAYCBD ← 3

CONSTANT 200. Item 3: DELR ← 4

 Item 4: DELC

INTERNAL 1.0 (FREE) 1 ← 4

90. ← 4

90. ← 4

90. ← 4

60. ← 4

90. ← 4

90. ← 4

90. ← 5

 125. 124. 124. 124. 124. 123. 123. 123. 123. 122. 

 122. 122. 122. 121. 121. 121. 121. 120. 120. 120. 

 120. 119. 119. 119. 119. 118. 118. 118. 118. 117.

 117. 117. 117. 116. 116. 116. 116. 115. 115. 115. 

 115. 114. 114. 114. 114. 113. 113. 113. 113. 112. 

 112. 112. 112. 111. 111. 111. 111. 110. 110. 110.

 110. 109. 109. 109. 109. 108. 108. 108. 108. 107. 

 107. 107. 107. 106. 106. 106. 106. 105. 105. 105. 

 105. 104. 104. 104. 104. 103. 103. 103. 103. 102.

 102. 102. 102. 101. 101. 101. 101. 100. 100. 100.

...

Previous 10 lines for model row 1 are repeated for model rows 2-7.

← 5

Constant 0.0 BOT layer 1 ← 6

2592000 30 1.1 TR Item 7: PERLEN NSTP TSMULT SS/TR ← 7

2592000 30 1.1 TR Item 7: PERLEN NSTP TSMULT SS/TR ← 7 ... ← 7

2592000 30 1.1 TR Item 7: PERLEN NSTP TSMULT SS/TR ← 7

1 Grid characterization data. 

2 Flags for quasi-3D simulation.

3 Row discretization data.

4 Column discretization data.

5 Top elevations for layer 1 read from following data.

6 Bottom elevations for layer 1 read from following data.

7 Time step and stress period information.

Appendix 2 37
Following (enclosed in a border) are the contents of the Layer-Property Flow Package input file for the hypothetical test 

simulation; explanations are noted outside of border:

File name: testsfr2.lpf

# LPF Package input file for hypothetical test simulation ← 1

# Test of SFR2 Package ← 1

0 -999 0 ← 2

1 ← 3

0 ← 3

1 ← 3

0 ← 4

0 ← 5

CONSTANT 8.0e-4 ← 6

CONSTANT 6.0e-6 ← 7

CONSTANT 0.0 ← 8

CONSTANT .2 ← 9

1 Comment lines. 

2 Flag for cell-by-cell flow terms, head in cell for wet-dry simulations, and flag related to wet-dry simulations.

3 Layer type, method for computing transmissivity, horizontal anisotropy.

4 Flag to determine vertical hydraulic conductivity.

5 Flag to determine if rewetting is active. 

6 Hydraulic conductivity data for layer 1 read from following data. 

7 Vertical hydraulic conductivity.

8 Specific storage for layer 1.

9 Specific yield for layer 1.

Following (enclosed in a border) are the contents of the General-Head Boundary Package input file for the hypothetical 

test simulation; explanations are noted outside of border:

File name: testsfr2.ghb

# GHB Package input file for hypothetical test simulation

# Test of SFR2 Package

7 -1 MXACTB IGHBCB ← 1

7 0 ITMP NP -- STRESS PERIOD 1 ← 2

1 1 100 80.2 400.1E-0

1 2 100 80.2 400.1E-0

1 3 100 80.2 400.1E-0

1 4 100 80.2 400.1E-0

1 5 100 80.2 400.1E-0

1 6 100 80.2 400.1E-0

1 7 100 80.2 400.1E-0

← 3

 -1 0 ITMP NP -- Stress Period 2

 -1 0 ITMP NP -- Stress Period 3

...

 -1 0 ITMP NP -- Stress Period 50

← 4

1 Maximum number of general-head boundary cells and unit number for storing cell-by-cell flow terms.

2 Read flags for stress period 1.

3 General-head boundary cells for stress period 1.

4 Read flags for stress period 2 (reusing information from stress period 1). 

Following (enclosed in a border) are the contents of the Preconditioned Conjugate-Gradient Package input file 

for the hypothetical test simulation; explanations are noted outside of border:

File name: testsfr2.pcg

 300 20 1 MXITER ITER1 NPCOND ← 1

.00001 .00001 .9 2 0 1 0.96 HCLOSE RCLOSE RELAX NBPOL IPRPCG MUTPCG DAMP ← 2

1 Maximum iterations, number of inner iterations, flag for matrix conditioning method.

2 Head change criterion, residual criterion for closure, relaxation parameter, upper bound on max eigenvalue NBPOL=2 when option is 

not used , printout interval for PCG, flag for printing convergence information, dampening factor.

38 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Following (enclosed in a border) are the contents of the Streamflow-Routing (SFR2) Package input file for the 

hypothetical test simulation; explanations are noted outside of border:

File name: Testsfr2.sfr

# SFR2 Package input file for hypothetical test simulation ← 1

# Test of SFR2 Package ← 1

-100 1 0 0 1.0 0.00001 -1 0 5 10 5 20 0 ← 2

 1 4 1 1 1 200.0 ← 3

 1 4 2 1 2 200.0 ← 3

 1 4 3 1 3 200.0 ← 3

 1 4 4 1 4 200.0 ← 3

 1 4 5 1 5 200.0 ← 3

 1 4 6 1 6 200.0 ← 3

 1 4 7 1 7 200.0 ← 3 ...

 1 4 94 1 94 200.0 ← 3

 1 4 95 1 95 200.0 ← 3

 1 4 96 1 96 200.0 ← 3

 1 4 97 1 97 200.0 ← 3

 1 4 98 1 98 200.0 ← 3

 1 4 99 1 99 200.0 ← 3

 1 4 100 1 100 200.0 ← 3

1 0 0 0 Item 5: stress period 1 ← 4

 1 2 0 0 .3 0.0 0.0 0.0 0.030 .04 Item 6a: ← 5

 .00000035 0.5 140. .3 .1 3.5 6.0e-6 Item 6b: ← 6

 .00000035 0.5 110. .3 .1 3.5 6.0e-6 Item 6c: ← 6

0. 2. 4. 6. 8. 10. 12. 14. ← 7

6.0 4.5 3.5 0. 0.3 3.5 4.5 6. ← 7

1 0 0 0 Item 5: stress period 2 ← 4

 1 2 0 0 2.533 0.0 0.0 0.0 0.030 .04 Item 6a: ← 5

 .00000035 Item 6b: ← 6

 .00000035 Item 6c: ← 6

1 0 0 0 Item 5: stress period 3 ← 4

 1 2 0 0 3.84 0.0 0.0 0.0 0.030 .04 Item 6a: ← 5

 .00000035 Item 6b: ← 6

 .00000035 Item 6c: ← 6 ...

1 0 0 0 Item 5: stress period 50 ← 4

 1 2 0 0 .9 0.0 0.0 0.0 0.030 .04 Item 6a: ← 5

 .00000035 Item 6b: ← 6

 .00000035 Item 6c: ← 6

1 Optional text line for comments. 

2 Basic stream network description (Input Item 1).

3 Stream reach descriptions, unsaturated flow variables(Input Item 2).

4 Read and print flags for stress period (Input Item 5).

5 General stream segment data (Input item 6a).

6 Stream segment data for upstream and downstream ends (Input Items 6b & 6c).

7 Stream segment data for 8-point cross section (Input Item 6d).

Appendix 2 39
Following (enclosed in a border) are the contents of the Output-Control input file for the hypothetical test simulation. 

Input was constructed in the optional style using words: 

File name: testsfr2.oc

# OC input file for hypothetical test simulation

# Test of SFR2 Package

HEAD PRINT FORMAT 14

DRAWDOWN PRINT FORMAT 0

PERIOD 1 STEP 10

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 1 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 5 STEP 30 

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 10 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 15 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 20 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 25 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 30 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 35 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 40 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 45 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

PERIOD 50 STEP 30

PRINT HEAD

PRINT BUDGET

SAVE BUDGET

40 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
Following are the contents of the Gage package input file for the hypothetical test simulation; explanations are 

noted as comments at the end of each record (Unit numbers are linked to a file name through the Ftype “DATA” in the 

MODFLOW name file):

File name: testsfr2.gag

8 NUMGAGE (number of gaging stations)

1 20 83 6 Stream Segment & Reach #; Output file for Gage 1; Output option

1 20 84 7 Stream Segment & Reach #; Output file for Gage 2; Output option

1 40 85 6 Stream Segment & Reach #; Output file for Gage 3; Output option

1 40 86 7 Stream Segment & Reach #; Output file for Gage 4; Output option 

1 60 87 6 Stream Segment & Reach #; Output file for Gage 5; Output option 

1 60 88 7 Stream Segment & Reach #; Output file for Gage 6; Output option 

1 80 89 6 Stream Segment & Reach #; Output file for Gage 7; Output option 

1 80 90 7 Stream Segment & Reach #; Output file for Gage 8; Output option

Appendix 2 41
Listing of Selected Output Files for the Test Simulation

Following are the contents of the main MODFLOW listing file of the hypothetical test 

simulation:

MODFLOW-2000

U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL

VERSION 1.15.01 04/05/2005

 This model run produced both GLOBAL and LIST files. This is the LIST file.

 # BAS6 Package input file for hypothetical test simulation 

 # Test of SFR2 Package 

 THE FREE FORMAT OPTION HAS BEEN SELECTED

 1 LAYERS 7 ROWS 100 COLUMNS

 50 STRESS PERIOD(S) IN SIMULATION

 BAS6 -- BASIC PACKAGE, VERSION 6, 1/11/2000 INPUT READ FROM UNIT 41

 5 ELEMENTS IN IR ARRAY ARE USED BY BAS

 GHB6 -- GHB PACKAGE, VERSION 6, 1/11/2000 INPUT READ FROM UNIT 15

 # GHB Package input file for hypothetical test simulation

 # Test of SFR2 Package

 No named parameters

 MAXIMUM OF 7 ACTIVE GHB CELLS AT ONE TIME

 CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0

 35 ELEMENTS IN RX ARRAY ARE USED BY GHB

42 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
 1618057 ELEMENTS OF RX ARRAY USED OUT OF 125000000

 501550 ELEMENTS OF IR ARRAY USED OUT OF 32000000

1

 # BAS6 Package input file for hypothetical test simulation 

 # Test of SFR2 Package...

...

SOLVING FOR HEAD 

 OUTPUT CONTROL FOR STRESS PERIOD 1 TIME STEP 10

 PRINT HEAD FOR ALL LAYERS

 PRINT BUDGET

 SAVE BUDGET

 HEAD DEP BOUNDS PERIOD 1 STEP 10

 BOUNDARY 1 LAYER 1 ROW 1 COL 100 RATE -0.300835 

 BOUNDARY 2 LAYER 1 ROW 2 COL 100 RATE -0.300730 

 BOUNDARY 3 LAYER 1 ROW 3 COL 100 RATE -0.297430 

 BOUNDARY 4 LAYER 1 ROW 4 COL 100 RATE -0.207811 

 BOUNDARY 5 LAYER 1 ROW 5 COL 100 RATE -0.297430 

Appendix 2 43

SFR2 -- STREAMFLOW ROUTING PACKAGE, VERSION 6, 04/27/2005

 INPUT READ FROM UNIT 45

# SFR2 Package input file for hypothetical test simulation

# Test of SFR2 Package

 NUMBER OF STREAM NODES IS 100

 NUMBER OF STREAM SEGMENTS IS 1

 NUMBER OF STREAM PARAMETERS IS 0

 NUMBER OF STREAM SEGMENTS DEFINED USING PARAMETERS IS 0

 MAXIMUM ERROR FOR STREAM LEAKAGE RATES IS 1.00E-05

 CONSTANT FOR MANNINGS EQUATION IS 1.0000E+00

 OPTION FOR UNSATURATED FLOW BENEATH STREAMBEDS IS ACTIVE 

 1618019 ELEMENTS IN RX ARRAY ARE USED BY SFR

 501513 ELEMENTS IN IR ARRAY ARE USED BY SFR

 3 ELEMENTS IN RX ARRAY ARE USED BY GAGE

 32 ELEMENTS IN IR ARRAY ARE USED BY GAGE
VOLUMETRIC BUDGET FOR UNSATURATED ZONE BENEATH STREAMS FOR TIME STEP 10 STRESS PERIOD 1

-------------------------------------------------------------------------------

--

CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T

 ------------------ ------------------------

STREAM LEAKAGE = 8553.5420 STREAM LEAKAGE = 3.4060E-02

CHANGE IN STORAGE = 8553.5176 CHANGE IN STORAGE = 3.4060E-02

RECHARGE TO GW = 0.0000 RECHARGE TO GW = 0.0000

 TOTAL IN = 8553.5420 TOTAL IN = 3.4060E-02

 TOTAL OUT = 8553.5176 TOTAL OUT = 3.4060E-02

 IN - OUT = 2.4414E-02 IN - OUT = 3.7253E-09

 PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00

...

HEAD IN LAYER 1 AT END OF TIME STEP 10 IN STRESS PERIOD 1

 ---------------------------------------------------------------------------

 1 2 3 4 5 6 7 8 9 10

 11 12 13 14 15 16 17 18 19 20

 21 22 23 24 25 26 27 28 29 30

 31 32 33 34 35 36 37 38 39 40

 41 42 43 44 45 46 47 48 49 50

 51 52 53 54 55 56 57 58 59 60

 61 62 63 64 65 66 67 68 69 70

 71 72 73 74 75 76 77 78 79 80

 81 82 83 84 85 86 87 88 89 90

 91 92 93 94 95 96 97 98 99 100

44 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1

 BOUNDARY 6 LAYER 1 ROW 6 COL 100 RATE -0.300730 

 BOUNDARY 7 LAYER 1 ROW 7 COL 100 RATE -0.300835 

...

STREAM LEAKAGE PERIOD 1 STEP 10

LAYER ROW COL. STREAM RCH. FLOW INTO STREAMBED FLOW OUT OF OVRLND. DIRECT

 SEG.NO. NO. STRM. RCH. LEAKAGE STRM. RCH. RUNOFF PRECIP 

 1 4 1 1 1 3.0000E-01 3.4738E-04 2.9965E-01 0.000E+00 0.000E+00

 1 4 2 1 2 2.9965E-01 3.4724E-04 2.9931E-01 0.000E+00 0.000E+00

 1 4 3 1 3 2.9931E-01 3.4710E-04 2.9896E-01 0.000E+00 0.000E+00

 1 4 4 1 4 2.9896E-01 3.4696E-04 2.9861E-01 0.000E+00 0.000E+00

STREAM STREAM STREAM STREAM STREAMBED FLOW TO CHNG. UNSAT.

ET HEAD DEPTH WIDTH CONDCTNC. WAT. TAB. STORAGE

0.000E+00 1.40279E+02 4.294E-01 2.326E+00 3.737E-04 0.000E+00 3.474E-04

0.000E+00 1.39979E+02 4.292E-01 2.326E+00 3.737E-04 0.000E+00 3.472E-04

0.000E+00 1.39679E+02 4.290E-01 2.326E+00 3.736E-04 0.000E+00 3.471E-04

0.000E+00 1.39379E+02 4.288E-01 2.326E+00 3.735E-04 0.000E+00 3.470E-04

...
 .........................................................................

 1 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 103.9 103.8 103.7 103.5 103.3 103.1 102.9 102.6

 102.3 102.0 101.6 101.0 100.4 99.3 96.9 92.6 85.5 75.1

...

7 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0 104.0

 104.0 104.0 103.9 103.8 103.7 103.5 103.3 103.1 102.9 102.6

 102.3 102.0 101.6 101.0 100.4 99.3 96.9 92.6 85.5 75.1

1

 VOLUMETRIC BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 10 IN STRESS PERIOD 1

 ------------------------------------------------------------------------------

 CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T

 ------------------ ------------------------

 IN: IN:

 --- ---

 STORAGE = 1302278.8750 STORAGE = 2.0056

 CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000

 HEAD DEP BOUNDS = 0.0000 HEAD DEP BOUNDS = 0.0000

 STREAM LEAKAGE = 0.0000 STREAM LEAKAGE = 0.0000

 TOTAL IN = 1302278.8750 TOTAL IN = 2.0056

 OUT: OUT:

 ---- ----

 STORAGE = 29.8633 STORAGE = 1.8395E-04

 CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000

 HEAD DEP BOUNDS = 1302252.2500 HEAD DEP BOUNDS = 2.0058

 STREAM LEAKAGE = 0.0000 STREAM LEAKAGE = 0.0000

 TOTAL OUT = 1302282.1250 TOTAL OUT = 2.0060

 IN - OUT = -3.2500 IN - OUT = -3.4690E-04

PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = -0.02

Appendix 2 45
 TIME SUMMARY AT END OF TIME STEP 10 IN STRESS PERIOD 1

 SECONDS MINUTES HOURS DAYS YEARS

 -----------------------------------------------------------

 TIME STEP LENGTH 37155. 619.25 10.321 0.43004 1.17738E-03

 STRESS PERIOD TIME 2.51132E+05 4185.5 69.759 2.9066 7.95791E-03

 TOTAL TIME 2.51132E+05 4185.5 69.759 2.9066 7.95791E-03

...

VOLUMETRIC BUDGET FOR UNSATURATED ZONE BENEATH STREAMS FOR TIME STEP 30 STRESS PERIOD 50

---------------------------------------------------------------------------------

CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T

 ------------------ ------------------------

STREAM LEAKAGE = 35678992.0000 STREAM LEAKAGE = 5.4863E-02

CHANGE IN STORAGE = 314356.0938 CHANGE IN STORAGE = 2.1510E-02

RECHARGE TO GW = 35373100.0000 RECHARGE TO GW = 3.3353E-02

 TOTAL IN = 35678992.0000 TOTAL IN = 5.4863E-02

 TOTAL OUT = 35687456.0000 TOTAL OUT = 5.4863E-02

 IN - OUT = -8464.0000 IN - OUT = 3.7253E-09

 PERCENT DISCREPANCY = -0.01 PERCENT DISCREPANCY = 0.00

HEAD IN LAYER 1 AT END OF TIME STEP 30 IN STRESS PERIOD 50

 ---------------------------------------------------------------------------

 1 2 3 4 5 6 7 8 9 10

 11 12 13 14 15 16 17 18 19 20

 21 22 23 24 25 26 27 28 29 30

 31 32 33 34 35 36 37 38 39 40

 41 42 43 44 45 46 47 48 49 50

 51 52 53 54 55 56 57 58 59 60

 61 62 63 64 65 66 67 68 69 70

 71 72 73 74 75 76 77 78 79 80

 81 82 83 84 85 86 87 88 89 90

 91 92 93 94 95 96 97 98 99 100

 .........................................................................

 1 115.2 115.2 115.2 115.2 115.2 115.1 115.1 115.0 115.0 114.9

 114.9 114.8 114.8 114.7 114.6 114.5 114.4 114.3 114.2 114.1

 114.0 113.8 113.7 113.6 113.4 113.3 113.1 112.9 112.8 112.6

 112.4 112.2 112.0 111.8 111.6 111.3 111.1 110.9 110.6 110.4

 110.1 109.8 109.5 109.2 109.0 108.6 108.3 108.0 107.7 107.3

 107.0 106.6 106.3 105.9 105.5 105.1 104.7 104.3 103.8 103.4

 102.9 102.5 102.0 101.5 101.0 100.5 100.0 99.5 98.9 98.4

 97.8 97.2 96.6 96.0 95.4 94.8 94.1 93.5 92.8 92.1

 91.4 90.7 89.9 89.2 88.4 87.6 86.9 86.0 85.2 84.4

 83.5 82.6 81.7 80.8 79.9 78.9 78.0 77.0 76.0 75.0

...

46 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1
7 115.2 115.2 115.2 115.2 115.2 115.1 115.1 115.0 115.0 114.9

 114.9 114.8 114.8 114.7 114.6 114.5 114.4 114.3 114.2 114.1

 114.0 113.8 113.7 113.6 113.4 113.3 113.1 112.9 112.8 112.6

 112.4 112.2 112.0 111.8 111.6 111.3 111.1 110.9 110.6 110.4

 110.1 109.8 109.5 109.2 109.0 108.6 108.3 108.0 107.7 107.3

 107.0 106.6 106.3 105.9 105.5 105.1 104.7 104.3 103.8 103.4

 102.9 102.5 102.0 101.5 101.0 100.5 100.0 99.5 98.9 98.4

 97.8 97.2 96.6 96.0 95.4 94.8 94.1 93.5 92.8 92.1

 91.4 90.7 89.9 89.2 88.4 87.6 86.9 86.0 85.2 84.4

 83.5 82.6 81.7 80.8 79.9 78.9 78.0 77.0 76.0 75.0

1

 VOLUMETRIC BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 30 IN STRESS PERIOD 50

 ------------------------------------------------------------------------------

 CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T

 ------------------ ------------------------

 IN: IN:

 --- ---

 STORAGE = 18276614.0000 STORAGE = 0.1495

 CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000

 HEAD DEP BOUNDS = 0.0000 HEAD DEP BOUNDS = 0.0000

 STREAM LEAKAGE = 35418412.0000 STREAM LEAKAGE = 3.3353E-02

 TOTAL IN = 53695024.0000 TOTAL IN = 0.1829

 OUT: OUT:

 ---- ----

 STORAGE = 17087914.0000 STORAGE = 0.0000

 CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000

 HEAD DEP BOUNDS = 36606024.0000 HEAD DEP BOUNDS = 0.1829

 STREAM LEAKAGE = 0.0000 STREAM LEAKAGE = 0.0000

 TOTAL OUT = 53693936.0000 TOTAL OUT = 0.1829

 IN - OUT = 1088.0000 IN - OUT = 2.5049E-05

 PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.01

 TIME SUMMARY AT END OF TIME STEP 30 IN STRESS PERIOD 50

 SECONDS MINUTES HOURS DAYS YEARS

 -----------------------------------------------------------

 TIME STEP LENGTH 2.49961E+05 4166.0 69.434 2.8931 7.92080E-03

 STRESS PERIOD TIME 2.59200E+06 43200. 720.00 30.000 8.21355E-02

 TOTAL TIME 1.29600E+08 2.16000E+06 36000. 1500.0 4.1068

Appendix 2 47
48 Report Title

Following are the abridged contents of Gaging Station 5 output file for the hypothetical test simulation: 

File name: testSFR2.gg5

“GAGE No. 5: K,I,J Coord. = 1, 4, 60; STREAM SEGMENT = 1; REACH = 60 “” DATA: 

 Time Stage Depth GW Head M-P Flow Leakage GW Rech. Chnge UZ Stor. 

Vol. UZ Stor.”

 1.5757E+04 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3928E-04 

5.3462E+00

 3.3091E+04 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.1227E+01

 5.2157E+04 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.7696E+01

 7.3130E+04 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

2.4812E+01

 9.6201E+04 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

3.2640E+01

 1.2158E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

4.1251E+01

 1.4949E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

5.0722E+01

 1.8020E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

6.1141E+01

 2.1398E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

7.2601E+01

 2.5113E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

8.5208E+01

 2.9200E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

9.9075E+01

 3.3696E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.1433E+02

 3.8641E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.3111E+02

 4.4081E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.4957E+02

 5.0065E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.6987E+02

 5.6647E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

1.9220E+02

 6.3888E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

2.1677E+02

 7.1852E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

2.4379E+02

 8.0613E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

2.7352E+02

 9.0251E+05 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

3.0621E+02

 1.0085E+06 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

3.4218E+02

 1.1251E+06 1.2257E+02 4.1743E-01 1.0400E+02 2.7957E-01 3.3929E-04 0.0000E+00 3.3929E-04 

3.8175E+02

... 
Heading 1 49

Following are the abridged contents of Gaging Station 6 output file for the hypothetical test simulation: 

File name: testSFR2.gg6

“GAGE No. 6: K,I,J Coord. = 1, 4, 60; STREAM SEGMENT = 1; REACH = 60 “

 “DATA: Time Depth Width Ave. Water Content Cell 1 Water Content”

...

2.5920E+06

 8.825E-01 1.2472E-01 1.8523E-01

 1.765E+00 1.2472E-01 1.8523E-01

 2.647E+00 1.2472E-01 1.8523E-01

 3.530E+00 1.2472E-01 1.8523E-01

 4.412E+00 1.2472E-01 1.8523E-01

 5.295E+00 1.2472E-01 1.8523E-01

 6.177E+00 1.2472E-01 1.8523E-01

 7.060E+00 1.2472E-01 1.8523E-01

 7.942E+00 1.2472E-01 1.8523E-01

 8.825E+00 1.2472E-01 1.8523E-01

 9.707E+00 1.1100E-01 1.3791E-01

 1.059E+01 1.0000E-01 1.0000E-01

 1.147E+01 1.0000E-01 1.0000E-01

 1.235E+01 1.0000E-01 1.0000E-01

 1.324E+01 1.0000E-01 1.0000E-01

 1.412E+01 1.0000E-01 1.0000E-01

 1.500E+01 1.0000E-01 1.0000E-01

 1.588E+01 1.0000E-01 1.0000E-01

 1.677E+01 1.0000E-01 1.0000E-01

 1.765E+01 1.0000E-01 1.0000E-01

...

1.2960E+08

 9.132E-01 1.2907E-01 1.9768E-01

 1.826E+00 1.2966E-01 1.9768E-01

 2.739E+00 1.3019E-01 1.9768E-01

 3.653E+00 1.3454E-01 1.9768E-01

 4.566E+00 1.3523E-01 1.9768E-01

 5.479E+00 1.3573E-01 1.9768E-01

 6.392E+00 1.3779E-01 1.9768E-01

 7.305E+00 1.4028E-01 1.9768E-01

 8.218E+00 1.4028E-01 1.9768E-01

 9.132E+00 1.4028E-01 1.9768E-01

 1.004E+01 1.4028E-01 1.9768E-01

 1.096E+01 1.4028E-01 1.9768E-01

 1.187E+01 1.4028E-01 1.9768E-01

 1.278E+01 1.4075E-01 1.9768E-01

 1.370E+01 1.4134E-01 1.9768E-01

 1.461E+01 1.4134E-01 1.9768E-01

 1.552E+01 1.4134E-01 1.9768E-01

 1.644E+01 1.4134E-01 1.9768E-01

 1.735E+01 1.4134E-01 1.9768E-01

 1.826E+01 1.4134E-01 1.9768E-01
50 Documentation of the Streamflow-Routing (SFR2) Package—A Modification to SFR1

This page left intentionally blank
Manuscript approved for publication, August 29, 2005

Prepared by the USGS, Tacoma Publishing Service Center

USGS Publishing staff

Bill Gibbs 

Linda Rogers

Angelia M. Thacker

For more information concerning the research in this report, contact the 

Director, Nevada Water Science Center

U.S. Geological Survey, 333 W. Nye Lane, Room 203

Carson City, NV 89706-0866

http://nevada.usgs.gov
TM 6-A13 Version 1.20 Niswonger and Prudic

Printed on recycled paper

Documentation of the Streamflow-Routing (SFR2) Package to Include 

Unsaturated Flow Beneath Streams—A Modification to SFR1