A KINEMATIC DISTRIBUTED WATERSHED RAINFALL-RUNOFF MODEL:
ACCURACY AND PREDICTIVE UNCERTAINTY USING WEATHER RADAR
by
Brian Edward Skahill
B.S., University of Portland, 1990
M.S., Colorado State University, 1994
M.S., University of Colorado at Denver, 1999
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Civil Engineering
2000
2000 by Brian Edward Skahill
All rights reserved.
This thesis for the Doctor of Philosophy
degree by
Brian Edward Skahill
has been approved
by
o c
Date
Skahill, Brian Edward (Ph.D., Civil Engineering)
A Kinematic Distributed Watershed Rainfall-Runoff Model: Accuracy and Predictive
Uncertainty Using Weather Radar
Thesis directed by Professor Lynn E. Johnson
ABSTRACT
An event-based, kinematic, infiltration-excess, distributed rainfall-runoff
model was developed to acknowledge and account for the spatial variability and
uncertainty of several parameters relevant to storm surface runoff production and
surface flow. The model is compatible with raster Geographic Information Systems
(GIS) and spatially and temporally varied rainfall data. Principal raster data inputs
include hydrography data derived from a digital elevation model, soil texture, forest
type, forest density, lakes, reservoirs, land use and land cover, and impervious regions
within a basin. These GIS data sets are used to support computations related to
spatially-varied interception, accounting of input rainfall to lakes and reservoirs and
routing into reservoirs, spatially-varied infiltration, kinematic-wave overland flow-
routing, channel losses, and kinematic-wave channel flow routing. The main model
outputs include a volume summary, discharge hydrographs for each sub-basin as well
IV
as for the main basin outlet, and raster maps of various variables; such as, cumulative
infiltration and excess-rainfall.
Monte Carlo simulation and a likelihood measure are utilized to calibrate the
model; allowing for a range of possible system responses from the calibrated model.
The model structure supports a meaningful distributed model calibration, and allows
for the investigation of uncertainty in rainfall-runoff modeling. Using rain gauge
adjusted radar-rainfall estimates, the model was applied and evaluated to a limited
number of historical events for two watersheds within the Denver Urban Drainage
and Flood Control District (UDFCD) that contain mixed land use classifications. For
the events considered, the 95% uncertainty bounds obtained from the model envelop
almost all observed responses, not only at the main basin outlets, but also at a location
internal to one of the two basins, suggesting an acceptable model structure. While
based on a limited number of Monte Carlo simulations and considered events, Nash
and Sutcliffe efficiency score ranges, based on a comparison of observed and
simulated hydrographs, of -0.19-0.95/-0.75-0.81 were obtained from the calibrated
models for the two basins. The hydraulic conductivity of the soil was the single
model parameter most relevant to reducing runoff volume uncertainty. Uncertainties
in initial conditions and rainfall estimates were also shown to play a significant role in
model uncertainty.
v
This abstract accurately represents the content of the candidates thesis. I recommend
its publication.
Sign'd
vi
DEDICATION
This thesis is dedicated to Beth and Glen.
ACKNOWLEDGEMENTS
I would like to thank Dr. Lynn E. Johnson, my graduate thesis advisor, for
recognizing my ability to contribute and inviting me to participate. Thanks are also
due to Lynn for the guidance he provided throughout the duration of the project. I
would also like to thank Drs. James Heaney, Weldon Lodwick, Russell Qualls, and
Samuel Welch for serving on my committee and providing constructive criticism.
I would like to thank Michael Dixon from the Research Applications Program
of the National Center for Atmospheric Research for providing the processed radar
data and Kevin Stewart from the Denver Urban Drainage and Flood Control District
for providing the rain gauge and streamflow data. Thanks are also due to Bruce
Rindahl from HDR Engineering Inc. for providing rain gauge and streamflow data. I
would also like to thank Scott ODonnell from the Forecasting Systems Laboratory
for several helpful discussions.
Financial support for this work was provided by the Modernization Division
of the Forecasting Systems Laboratory, Boulder, CO. The College of Engineering at
the University of Colorado at Denver granted scholarship money. This financial
support is gratefully acknowledged.
I am most grateful to my family, for their enduring patience and support.
CONTENTS
Figures.....................................................................xiii
Tables......................................................................xvii
Chapter
1. Introduction................................................................1
1.1. Statement of Problem....................................................1
1.2. Objectives .............................................................4
1.3. Scope of Study .........................................................4
1.4. Organization of Thesis .................................................7
2. Literature Review...........................................................8
2.1. Advanced Distributed Rainfall-Runoff Models.............................8
2.1.1. Surface Runoff Modeling ............................................12
2.1.2. Existing Models ....................................................42
2.1.3. Geographic Information Systems and Hydrologic Modeling..............48
2.1.4. Calibration.........................................................51
2.1.5. Modeling Studies ...................................................60
2.2. Weather Radar .........................................................68
IX
2.2.1.
Weather Surveillance Radar-1988 Doppler Rainfall Algorithm
70
2.2.2. Radar-Rainfall Studies ...............................................75
2.2.3. Hydrologic Applications...............................................79
3. Model Description ...........................................................83
3.1. Introduction..............................................................83
3.2. Rainfall..................................................................85
3.3. Interception .............................................................87
3.4. Infiltration .............................................................87
3.5. Overland Flow ............................................................92
3.6. Channel Flow ............................................................101
3.7. Additional Model Components..............................................105
3.8. Code Verification........................................................107
3.8.1. Infiltration ........................................................107
3.8.2. Overland Flow .......................................................108
3.8.3. Connectivity of Overland Flow and Channel Flow .....................110
3.9. Model Testing ...........................................................113
3.9.1. The Buffalo Creek flash flood of July 12, 1996 ......................113
3.9.2. Organized versus Random Sampling from Cell to Cell ..............125
4. Model Accuracy and Predictive Uncertainty Using Rain Gauge Adjusted Radar-
Rainfall Estimates .........................................................128
4.1. Basin Descriptions ......................................................128
x
4.2. Historical Rainfall-Runoff Events .......................................130
4.3. Radar Rainfall Data......................................................132
4.3.1. Rain Gauge Adj ustment...............................................133
4.4. Model Database Development.............................................137
4.5. Model Initialization...................................................143
4.5.1. GLEAMS ..............................................................144
4.6. Model Application and Validation.......................................149
4.6.1. Model Results .......................................................152
4.6.2. Discussion of Model Results .........................................165
4.6.3. Evaluation of Distributed Model Output ..............................169
4.7. Model Sensitivity Analysis ............................................173
5. Summary and Conclusions ....................................................181
Appendix
A. List of Symbols ......................................................185
B. Contact Information for obtaining Model Source Code ..................189
C. TOPAZ Model Input File DNMCNT.INP for the Ralston Creek Basin at the
800m Grid Cell Scale..................................................190
D. TOPAZ Model Input File DNMCNT.INP for the Goldsmith Gulch Basin at the
500m Grid Cell Scale..................................................195
E. Source Code that supports the preparation of TOPAZ Model Input Files .200
xi
F. Source Code that supports the processing of TOPAZ Model Output Files.........211
G. Source Code that generates slope maps of the minimum and maximum slope
within each raster grid cell .................................................216
H. Source Code that generates raster and tabular data files for lake and reservoir IDs
within a region of interest ..................................................229
I. Input files prepared for the hydrology component of the GLEAMS model.........242
J. Source code that formats daily precipitation data into the appropriate format for
GLEAMS model application .....................................................259
K. Source code that formats daily temperature data into the appropriate format for
GLEAMS model application .....................................................274
References
289
I
I
i
FIGURES
2.1. Plot of em associated with the FTBS difference scheme for R=0.1 -0.9...29
2.2. Plot of error in propagation speed associated with the FTBS difference scheme
for R=0.1-0.9..........................................................30
2.3. Comparison of methods for linear advection experiment..................33
2.4. Comparison of one hour unit hydrographs obtained using two FD schemes ...44
2.5. Best model obtained from 250 Monte Carlo Simulations.................58
2.6. Worst model obtained from 250 Monte Carlo Simulations ...............58
2.7. Ninety-five percent uncertainty bounds obtained from the model for the June 1,
1991 storm on the Ralston Creek basin. Model rejection criteria set at 0.3 using
the Nash and Sutcliffe efficiency criterion as the likelihood measure..59
2.8. Ninety-five percent uncertainty bounds obtained from the model for the June 1,
1991 storm on the Ralston Creek basin. Model rejection criteria set at 0.8 using
the Nash and Sutcliffe efficiency criterion as the likelihood measure.....59
3.1. A sequence of four images of the rainfall for the Buffalo Creek flash flood event
of July 12, 1996 .........................................................87
3.2. A sequence of four images of cumulative for the Buffalo Creek flash flood
event of July 12, 1996 ...................................................92
xm
3.3. A sequence of four images of rainfall excess for the Buffalo Creek flash flood
event of July 12, 1996 ..................................................101
3.4. A channel grid cell with a rectangular cross section.....................103
3.5. Comparison of model results with the analytical solution for a one-dimensional
plane geometry...........................................................110
3.6. Derived channel network for the Buffalo Creek basin used to support model
testing..................................................................Ill
3.7. Comparison of model results with the dimensionless hydrograph solution of
Constantinides and Stephenson (1982) for a V-shaped basin................113
3.8. May 18, 1996 forest fire footprint within the Buffalo Creek basin .......116
3.9. Model Predictive uncertainty for the July 12, 1996 storm on the Buffalo Creek
watershed using 250 Monte Carlo simulations using original radar-rainfall
estimates......................................................................119
3.10. Model Predictive uncertainty for the July 12, 1996 storm on the Buffalo Creek
watershed using 250 Monte Carlo simulations using rain gauge adjusted radar-
rainfall estimates .........................................................119
3.11. Model discharge divided by the peak discharge for the given resolution.....124
3.12. Model discharge divided by the peak discharge obtained at the 100 m grid cell
scale.......................................................................124
3.13. Model predictive uncertainty for the July 28, 1997 event on Goldsmith Gulch.
Results based on 1,000 Monte Carlo simulations using (a) organized and (b)
xiv
random sampling. The Nash and Sutcliffe efficiency criterion was used for
the likelihood measure and the rejection criteria was set at 0.3............126
3.14. Model predictive uncertainty for the July 22, 1991 event on Ralston Creek
basin. Results based on 1,000 Monte Carlo simulations using (a) organized
and (b) random sampling. The Nash and Sutcliffe efficiency criterion was
used for the likelihood measure and the rejection criteria was set at 0.3.........127
4.1. Basin boundaries for the Ralston Creek and Goldsmith Gulch watersheds ....130
4.2. ALERT rain gauge locations within and surrounding the Ralston Creek ........135
4.3. ALERT rain gauge locations for the Goldsmith Gulch basin....................135
4.4. (a) DEM, (b) forest density, (c) soil texture classification, and (d) LULC data
for the Colorado Front Range region.........................................138
4.5. Raster map of the locations of lakes and reservoirs for a regional scale model
database developed for model application....................................142
4.6. Locations for climatological data provided in the GLEAMS model supplemental
database for the state of Colorado.........................................147
4.7. Locations for GLEAMS model application.....................................148
4.8. Ninety-five percent uncertainty bounds obtained from the model for the (a) June
1, 1991, (b) May 26, 1996, (c) August 4, 1997, (d) June 21, 1991, (e) July 22,
1991, (f) August 26, 1996, (g) September 18, 1996, and (h) July 30, 1997
storms on the Ralston Creek basin..............................................155
xv
4.9. Ninety-five percent uncertainty bounds obtained from the model for the (a) July
28,1997, (b) July 20, 1991, (c) July 12, 1996, (d) July 27, 1997, (e) August 18,
1995, and (f) July 30, 1997 storms on the Goldsmith Gulch basin .......158
4.10. Retained model ranges, and comparison with observed values, for (a) Time to
Peak, (b) Peak Discharge, and (c) Runoff Volume. (Ralston Creek Basin) ....160
4.11. Retained model ranges, and comparison with observed values, for (a) Time to
Peak, (b) Peak Discharge, and (c) Runoff Volume (Goldsmith Gulch) .....162
4.12. Best and worst retained simulations obtained from the calibrated model for the
(a,b) July 27, 1997, (c,d) August 18, 1995, and (e,f) July 30, 1997 validation
storms on the Goldsmith Gulch basin....................................164
4.13. Location internal to the Ralston Creek basin for which observed and computed
discharge data were compared ..........................................170
4.14. Ninety-five percent uncertainty bounds obtained from the model at the Simms
Street location for the (a) May 26, 1996, (b) August 4, 1997, (c) August 26,
1996, and (g) September 18, 1996 storms on the Ralston Creek basin ...........172
4.15. Sensitivity analysis results for (a) time to peak, (b) peak discharge, (c) volume
infiltrated, and (d) runoff volume ................................................180
xvi
TABLES
3.1. Ranges for Mannings roughness coefficient values based on land use .........99
3.2. Model results obtained to evaluate consistency across grid scale ...........123
4.1. Storm events evaluated on the Ralston Creek basin and radar locations......131
4.2. Storm events evaluated on the Goldsmith Gulch basin and radar locations ....132
4.3. Bias Correction Factors computed for each storm event.......................136
4.4. Computed Likelihood Measure values for each individual event. (Ralston Creek
Basin) .....................................................................160
4.5. Computed Likelihood Measure values for each individual event. (Goldsmith
Gulch Basin) ...............................................................162
4.6. Updated Computed Likelihood Measure values during model calibration.
(Ralston Creek Basin) ......................................................163
4.7. Updated Computed Likelihood Measure values during calibration. (Goldsmith
Gulch Basin) ...............................................................163
4.8. Nash and Sutcliffe efficiency scores obtained for the validation events on the
Goldsmith Gulch basin ......................................................163
xvii
1. Introduction
1.1. Statement of Problem
Forecasting and warning of flash floods is difficult for forecasters due to the
combination of intense rainfalls and brief time frames. In the United States, they are
responsible for more deaths than any other weather-related phenomenon. Due to the
expansion of residences and recreationalists into flash flood prone areas, localized
flash floods are accounting for a greater portion of flood damages and loss of life.
The modernization of the U.S. National Weather Service (NWS) has the goal
of providing improved forecasts and more reliable detection and prediction of severe
weather and flooding. Two major components of the modernization include the
installment of the weather radar system WSR-88D (Weather Surveillance Radar-1988
Doppler) and AWIPS-90 (Advanced Weather Interactive Processing System for the
1990s) workstations in NWS Weather Forecast Offices (WFOs).
As part of the NWS Next Generation Weather Radar program (NEXRAD),
nationwide delivery of over 160 WSR-88D radars was completed in 1997. The WSR-
88D radars provide almost complete coverage of the continental United States, and
each radar has a range of approximately 460 kilometers. The radar reflectivity data
product from the WSR-88D radars is capable of supporting the estimation of rainfall
rate with a spatial and temporal resolution on the order of one kilometer and six
minutes, respectively.
Data relevant to runoff generation and surface flow; such as, elevation, soils,
land use and land cover, and vegetative cover, are readily available today at multiple
scales, and Geographic Information Systems (GIS) technology provides a means for
data storage, handling, mapping, and evaluation. Also readily available are tools,
within or external to a GIS, that automatically extract hydrography data from a digital
elevation model. Computing power continues to increase at a rapid rate, and there are
numerous computational methods available to model subsurface and surface flow.
The radar-estimated rainfall data, data relevant to the partitioning of rainfall at
the land surface, and standard computational methods are all mutually compatible,
and have the potential to support regional scale flash flood forecasting operations. A
recent assessment of alternate flash flood warning procedures (Johnson 1998 and
2000) indicated the preference for the advanced distributed rainfall-runoff model
approach by NWS WFO forecasters. The advanced distributed model approach could
operate in the background, and output display of runoff as a grid across the basins and
region would be preferred as that format is compatible in format and interaction with
most of the other products on the AWIPS-90 workstation. Furthermore, collection of
surface runoff in the stream channels and routing of these contributing flows
downstream provides a single display graphic that can interface with flood impact
features typically located along the channels and at bridges. Research is required to
2
examine advanced distributed rainfall-runoff model accuracy while using radar
estimated rainfall data.
For a distributed model, the number of parameters depends upon whether or
not their spatial variability is considered. Current advanced distributed model
descriptions emphasize the models compatibility with GIS; thus, allowing for the
spatial variability of input, boundary conditions, and watershed characteristics
relevant to runoff generation and surface flow. Ironically, the purported advantage of
distributed models is also their clear weakness. They are highly parameterized, and,
as a result, often during the manual calibration process of an advanced distributed
rainfall-runoff model, for reasons of practicality, a single spatially uniform value is
arbitrarily assigned to several model parameters.
Field data is typically unavailable. Consequently, parameter values within
distributed models are often based on other physical information about the watershed.
For example, soil texture classification data are often used in modeling studies to
support the assignment of parameter values related to model infiltration. For a unique
soil texture class, these parameter values exhibit a high degree of uncertainty.
Needed, is an advanced distributed rainfall-runoff model structure that
1. accounts for the spatial variability and uncertainty of several parameters relevant
to storm surface runoff production and surface flow, and
2. also supports a meaningful distributed model calibration (Johnson 2000; Uccellini
1999).
3
1.2. Objectives
The objectives of this study are to:
1. Develop a complex distributed surface runoff model, compatible with raster GIS
and radar-rainfall data, that
A. accounts for uncertainties in rainfall-runoff modeling; and
B. supports a model calibration while allowing for the spatial variability of
model parameters;
2. Examine the accuracy and predictive uncertainty with which the model simulates
storm surface runoff for two basins within the Denver Urban Drainage and Flood
Control District; and
3. Approximate the relative contribution of uncertainties in rainfall, initial soil water
content, and model parameters to the models overall predictive uncertainty.
1.3. Scope of Study
In scope, this study is limited to:
1. convective rainfall, principally estimated using weather radar;
2. watersheds for which Horton overland flow is the dominant surface runoff
mechanism;
3. applications for which the kinematic wave model is valid; and
4. basin areas ranging from 15 km to 225 km .
4
Statistical procedures for storm climatology characterized instantaneous storm
properties (e.g., storm volume, precipitation area, precipitation flux), aggregate
properties for entire storm tracks (e.g., duration, speed and direction, precipitation
swath area and precipitation depth), and geographical distribution of storm and track
properties over the Colorado Front Range region (Dixon and Johnson 1998). The
storm tracking approach of Dixon and Johnson (1998) indicated that the mean storm
size was 32 km 90 percent of the storms are smaller than 150 km the mean velocity
of storm cells was 20 km/hr, and that 10 percent of the storms yield 75 percent of the
rainfall volume.
Radar rainfall imagery from the Mile High Radar, a NEXRAD prototype
radar, and the Denver WSR-88D radar (FTG) were obtained from the Research
Applications Program of the National Center for Atmospheric Research (NCAR),
where the Thunderstorm Identification, Tracking, Analysis and Nowcasting (TITAN)
program (Dixon and Weiner 1994) is run in real time. Evaluation of the model
described in this work required the availability of a high quality rainfall-runoff
database. The Denver Urban Drainage and Flood Control District (UDFCD)
maintains an Automated Local Evaluation in Real Time (ALERT) flood warning
system that primarily serves the Denver, Colorado metropolitan region. Two different
size basins located within the UDFCD, with different land use classifications, were
selected for the study based on the available radar, rain gauge, and streamflow data.
5
Michaud and Sorooshian (1994a) noted that it is commonly assumed that an
accurate rainfall-runoff model is one which, when excited with reasonably accurate
rainfall and initial soil water content, will produce simulated flows that closely match
the observed streamflow. They listed several factors that are responsible for the wide
variability of the accuracy of rainfall-runoff model simulations reported in the
literature. These include
1. whether calibration was performed,
2. whether results reflect the calibration events or split sample validation events,
3. the number and variety of events used to test the model,
4. differing opinions as to what is a good simulation and what is a poor simulation,
5. the reliability of the data against which model simulations are compared, and
6. whether results pertain to
A. real-time forecasting,
B. prediction of the flow resulting from a specific historical storm, or
C. prediction of the frequency distribution of, say, peak flows from a set of
historical storms. The model assessment in this study was carried out in ordinary
prediction mode (6.B.)
6
1.4 Organization of Thesis
Chapter two includes a literature review that focuses on surface runoff
modeling, distributed rainfall-runoff model calibration and predictive uncertainty, and
weather radar.
A comprehensive model description, including data requirements, parameter
estimation, model testing, and initial model application is provided in chapter three.
Chapter four focuses on model application, including database development,
initialization, calibration and validation, accuracy and predictive uncertainty, and a
model sensitivity analysis.
Chapter five contains conclusions and directions for further research.
A significant component of the research effort was the development of a
robust, general purpose, event-based, advanced distributed rainfall-runoff model.
Contact information for obtaining the model source code, written in ANSI C solely by
this author, is included in Appendix B.
7
2. Literature Review
The rapidly increasing power of computers, development of geographic
information systems (GIS), increasing availability of watershed databases in raster or
grid cell GIS format, weather radar rainfall estimates in raster format, and the
growing demand for more distributed predictions (Refsgaard and Abbott 1996) are a
few of the reasons why complex distributed rainfall-runoff modeling remains an
active applied research area. The literature review primarily focuses on surface runoff
modeling, complex distributed rainfall-runoff model calibration and predictive
uncertainty, and weather radar.
2.1. Advanced Distributed Rainfall-Runoff Models
Distributed-parameter rainfall-runoff models partition a basin into smaller
elements in the horizontal, and possibly, also the vertical, to account for the spatial
variability of processes, input, boundary conditions, and watershed characteristics
relevant to runoff generation and surface flow (Singh 1995). If at least one of the
hydrologic processes is based on a conservation law, or one of its approximations, in
one or more space dimensions, a distributed model is classified as advanced,
complex, or physically-based. Arguably, advanced distributed models are
conceptually more correct than traditional lumped methods, not only because of their
8
ability to incorporate large quantities of data describing model input and distributed
watershed properties, but also because of their sophisticated treatment of the flow
dynamics internal to a basin. In contrast to lumped conceptual models, which require
extensive records for the estimation of parameters, in theory, the parameters for
complex distributed models can be obtained from field measurements alone; allowing
for their application to ungauged basins, and for their use in evaluating alternative
land management strategies. Their ability to provide spatially distributed hydrographs
throughout a basin as well as detailed images of the spatial distribution of runoff
make them an attractive choice as a flash flood forecasting tool. Disadvantages of the
advanced distributed modeling approach include the significant time and
computational resources required for model development, data collection, database
development, parameter estimation, and model application.
Beven (1989) expressed concern regarding the practical application of the
current generation of physically-based models. It was emphasized that one must
differentiate between the two main objectives of simulation models: (1) gaining a
better understanding of hydrologic systems, and (2) prediction. It was noted that the
primary application of physically-based models to date has been devoted to the first
aim. Beven (1989) commented that several authors have not made a distinction
between these two objectives. Beven (1989) claimed that physically-based models are
conceptual models operating at the grid cell scale. For example, Beven (1989) asked,
What does a grid square average capillary potential mean? It certainly cannot be
9
measured in the field. Beven (1989) also noted that most of the distributed models
cannot predict overland flow over a part of a grid cell. Bevens (1989) comments
highlighted the potential for lumping of parameters and processes in physically-based
models.
Beven (1989) noted that physically-based models suffer the same
shortcomings as lumped conceptual models. For example, both modeling approaches
suffer from
1. model structure errors,
2. the inability to capture system heterogeneity,
3. the potential for a less than meaningful calibration and validation,
4. overparameterization, and
5. parameter interaction.
Commenting on a paper involving the application and calibration of the distributed
SHE model (Abbott et al. 1986), Beven (1989) noted that applying physical reasoning
while performing the calibration process did not appear to improve the calibration,
even when employed by expert users. Beven (1989) suggested, in the worst case, if
uncertainty is incorporated into the model predictions while performing what-if-
games, would one be able to distinguish statistically between the different scenarios?
Grayson et al. (1992b) reiterated several of the criticisms of physically-based
distributed models that were originally described by Beven (1989). They claimed that
the detailed modeling of the response of a basin to an input rainfall field is
10
transscientific; that is, it is a problem that cannot be answered by science. For
distributed models of the watershed, Grayson et al. (1992b) cautioned that values for
internal state variables, such as flow depths and velocities, may be an artifact of the
distributed model. Echoing the recommendation of Beven (1989), Grayson et al.
(1992b) urged predictions from distributed models to be associated with a realistic
assessment of the uncertainty in model predictions. Considering new directions, they
argued that complex distributed models should not be required to model the well-
behaved response of a basin. If distributed models are to be used for predictive
purposes, Grayson et al. (1992b) recommend that resources be allocated to field data
programs and the development of algorithms that apply at the model scale.
Woolhiser (1996) attempted to address the growing concern among some
researchers and practitioners regarding the physical basis and accuracy of advanced
distributed surface runoff computations. It was noted that the discretization of the
descriptive equations, the typical modeling approach, results in a lumping of
parameters at the chosen computational scale. Assuming no important mechanisms
have been neglected, Woolhiser (1996) suggested that the parameters must possess
some physical significance if the computational scale is comparable with the scales of
spatial variability of the system. Woolhiser (1996) further commented that as the
scale of the system being modeled increases and the computational scale becomes
large relative to the scale of variability of the system, variations in the parameter
fields will be missed or distorted and model parameters will lose some of their
physical significance.
Woolhiser (1996) admitted that it is difficult or impossible to calibrate a
model with numerous parameters; particularly when many interact. Woolhiser (1996)
also emphasized the importance of high resolution rainfall data in both space and time
to support hydrologic simulations not only at the large scale, but also at the small
scale. Woolhiser (1996) commented that poor spatial resolution of rainfall data would
contribute to difficulties in distributed model calibration. Woolhiser (1996)
commented that advanced distributed surface runoff models can be used in a
predictive sense for small basins. Woolhiser (1996) encouraged future research
involving carefully designed tests utilizing carefully checked hydrologic data.
2.1.1. Surface Runoff Modeling
Physically-based distributed rainfall-runoff models commonly utilize
physically measurable or remotely-sensed parameter values which describe the
spatially-varied character of the watershed topography, soils, vegetation, drainage
network, and rainfall. At a minimum, these parameter values serve as input to
numerical algorithms for the solution of the equations which govern the physics of
infiltration, overland flow, and channel flow for the modeling of the spatial and
temporal response of a watershed.
12
2.1.1.1. Surface Flow
Two components which can be found within almost all physically-based
distributed rainfall-runoff models account for the routing of excess rainfall over the
land surface and through an organized stream network.
2.1.1.1.1. Governing Equations For 1-Dimensional
Shallow-Water Flow
The conservation law for one-dimensional shallow-water flow over a plane of
small slope is given by
where x, t, h, q, g, ie, S0 and S/ represent the distance in the flow direction, time, flow
depth, flow rate per unit width, the acceleration due to gravity, excess rainfall rate,
slope of the plane, and the friction slope, respectively. Equation 2.1 can be rewritten
in a non-conservative form in terms of the same dependent variables by applying the
chain rule. The non-conservative form is
(2.1)
13
(2.2)
3\ .
dt~ 3x + 1
(ti\
v
A,=
l?-gh
-1 ^
-2 q
F,=
e
s*(s-s,l
in which A, is the negative of the Jacobian matrix of partial derivatives of f with
respect to the components of v. Equation 2.1 can also be rewritten in terms of the
flow depth and flow velocity, u, as
3 u
= A, + K
3t 3 x
( ; \
(h^ f-u -h') le
u = A 2 = F2 =
[uj l-g -uj \ nj
by carrying out the differentiations in Equation 2.1 and canceling terms using the
continuity equation. Liggett (1975), Yevjevich (1975), Cunge et al. (1980),
Constantinides and Stephenson (1982), and Chaudhry (1993), among others, list the
assumptions related to the development of Equations 2.1-2.3, commonly known as
the de Saint-Venant equations. Principal among these is the reduction of the vertical
momentum equation, via a scale analysis, to the equation for a hydrostatic pressure
distribution
14
p and p represent the pressure and the density of liquid water, respectively. An
additional significant assumption is that the effects of boundary friction and
turbulence are accounted for through empirical resistance laws for steady state flow,
such as the Chezy equation or Mannings equation. Derivations of the de Saint-
Venant equations can be found in numerous sources; Liggett (1975), Cunge et al.
(1980), and Chaudhry (1993), to name a few.
2.1.1.1.2. Kinematic Wave Model
Although several approximations are made to obtain the equations describing
unsteady free-surface flow, additional assumptions are often made to simplify the
solution procedure. The kinematic wave approximation to the continuity equation and
equation of motion is obtained by neglecting the local inertia, convective inertia,
spatial gradient of flow depth, and momentum-source terms in the equation of motion
so that we can take S0 = S/. In neglecting these terms it is assumed that the flow must
be gradually varied, the pressure forces due to backwater effects must be small, and
that the spatial gradient of the flow depth must be small in comparison to the land
surface slope. Based on empirical grounds, a functional relation between flow and
surface water depth is postulated: q = ahp This functional relation is substituted into
the equation of continuity to obtain a partial differential equation (PDE) with one
dependent variable. Lighthill and Whitham (1955) are credited with the description of
the kinematic wave.
2.1.1.1.3. Diffusive Wave Model
Despite being the most popular routing technique, it is difficult to satisfy all of
the simplifying assumptions used to derive the kinematic wave approximation in
many practical modeling situations. For example, watershed land surface slopes tend
to decrease and surface flow depths tend to increase as one proceeds downstream. In
addition, many watersheds contain topographic features which restrict flow capacity.
Similar to the kinematic wave, the diffusive wave approximation neglects the local
inertia, convective inertia, and momentum-source terms, yet retains the spatial
gradient of flow depth term in the equation of motion so that we can take
dh
Sf S The addition of a resistance law q = an which provides a
ox
functional relation between q and h, completes the model, with the terms a and /3
dependent upon whether the flow is laminar or turbulent. In general, for overland
flow, the relative magnitude of the local inertia, convective inertia, and spatial
gradient of flow depth terms in the equation of motion is (Ogden 1992; Richardson
1989)
16
d u d u d h
---> u--> g---
d t d x d x
(2.5)
As a result, Ogden (1992) noted that there would be no improvement in the accuracy
of runoff model predictions based on the diffusive wave approximation over the
kinematic wave formulation.
2.1.1.1.4. Criteria For Kinematic Wave Modeling
Woolhiser and Liggett (1967) reduced the shallow-water equations (Equations
2.1-2.3) describing one-dimensional, spatially-varied, unsteady flow over a plane of
small slope into non-dimensional form to facilitate a comparison of alternative cases.
The non-dimensional form of Equation 2.3 contains two parameters, F0 and k, and is
given by
dh, d(u,h,) ,
dtt dxt
du, du, 1 dh,
dt, U" dx, + dx.
= k
1--
.2 \
(2.6)
(2.7)
where k 2 is the kinematic flow number. The equations are obtained by
^0^0
introducing the normalizing quantities:
17
L0 = length of plane,
h0 = normal depth at x = L0 corresponding to a flow rate of ieL0,
u0 ieLJh0,
Q0 = ieL0 = normalizing discharge,
and subsequently defining the dimensionless variables
K = h/ha, u, = u/u, x. = x/L, t. = u0t/L0, Fa = ujjgh^, Qt = uh/Qo,
and the dimensionless excess rainfall rate (/e) = ijqn = 1. Woolhiser and Liggett
(1967) observed that at k oo, Equation 2.7 reduces to u\ = h, which implies that a
unique functional relationship exists between surface water depth and discharge; that
is, the general solution and the kinematic wave equation solution coincide.
Furthermore, based on numerical solutions of Equations 2.6 and 2.7, they
demonstrated that rising hydrographs for values of k larger than twenty and Froude
numbers greater than or equal to one half can be closely approximated by the
kinematic wave formulation. Most hydrologically significant cases were stated to
satisfy this criterion. For values of k less than ten, they recommended that the
complete set of equations be used.
18
Ponce et al. (1978) provided criteria for the applicability of the kinematic and
diffusive wave approximations to open channel flow. Morris and Woolhiser (1980)
expanded upon the work initiated by Woolhiser and Liggett (1967) and demonstrated
that for low values of the Froude number, the kinematic wave formulation is a
reasonable approximation provided F^k > 5. Daluz Vieira (1983) extended the work
initiated by Woolhiser and Liggett (1967) and Morris and Woolhiser (1980). Based
upon more than one hundred and fifty tests, Daluz Vieira (1983) provided graphical
information in the (F0,k) plane as to an appropriate approximation for two separate
lower boundary conditions (BC): a zero-depth-gradient BC and a critical-flow lower
BC. Pearson (1989) suggested a new criterion for using the kinematic wave
approximation that was consistent with previous criteria set by Woolhiser and Liggett
(1967), Morris and Woolhiser (1980), and Daluz Vieira (1983), yet an improvement,
specifically k>3 + 5/F*. The results of Woolhiser and Liggett (1967), Morris and
Woolhiser (1980), Daluz Vieira (1983), and Pearson (1989) indicate that the
kinematic routing method is a reasonable approximation for most hydrologic
applications of overland flow.
2.1.1.1.5. Kinematic Wave Model Applications
The relative simplicity of the kinematic wave formulation has resulted in its
widespread use to model the rainfall-runoff process for simple as well as complex
19
watershed geometries (Wooding 1965a,b; Brakensiek 1967; Kibler and Woolhiser
1970; Singh and Woolhiser 1976; Constantinides 1982; Constantinides and
Stephenson 1982; Goodrich 1990; Vieux et al. 1990; Goodrich et al. 1991; Michaud,
J.D. 1992; Michaud and Sorooshian 1994a). However, its application has not been
without controversy (Hromadka and DeVries 1988; Dawdy 1990; Goldman 1990;
Merkel 1990; Unkrich and Woolhiser 1990; Woolhiser and Goodrich 1990;
Hromadka and DeVries 1990; Ponce 1991; Goodrich 1992; Woolhiser 1992; Ponce
1992). Hromadka and DeVries (1988) demonstrated the grid dependence of the finite
difference scheme implemented within the HEC-1 kinematic wave channel routing
program (HEC 1979) by routing inflow hydrographs through a 25,000 ft (7,620 m)
long rectangular channel. Based on their results, Hromadka and DeVries (1988)
recommended a re-evaluation of the kinematic wave model for channel routing in
watershed models. The discussions that followed the paper by Hromadka and
DeVries (1988) criticized their results due to the unrepresentative examples that were
considered. Goldman (1990) noted that the kinematic wave model is a useful
approach for (urban) headwater basins. Woolhiser and Goodrich (1990) commented
that kinematic modeling does an excellent job for overland flow and for channels
with overland lateral inflow, provided adequate numerical models are used.
Ponce (1991) stated that error due to numerical dissipation and dispersion in
a numerical solution using the finite difference method is at the crux of the
controversy surrounding the kinematic wave model. Given the grid dependence of
20
calculated runoff hydrographs using finite difference approximations, it was noted to
be pointless to calibrate a kinematic wave model by varying a physical parameter
such as Mannings n and comparing observed data with simulated results. Ponce
(1991) recommended that the kinematic wave model be applied to basins less than 2.5
km2. Woolhiser (1992) observed that the controversy might be related more to the
choice and appropriate application of a finite difference method rather than the
kinematic wave model itself.
In a discussion following the paper by Ponce (1991), Goodrich (1992)
emphasized the statistical interpretation of Mannings n value for any application that
involves a geometric abstraction of a watershed into model elements. This is due to
the mixed representation of overland and concentrated flow within an overland flow
grid cell. Goodrich (1992) also stated that Ponce (1991) provided little supporting
evidence for the recommendation of limiting kinematic wave modeling applications
j
to basins less than 2.5 km Goodrich (1990) obtained good modeling results applying
a kinematic wave model for basins up to 6.3 km Michaud (1992) and Michaud and
Sorooshian (1994a) applied a kinematic wave model to the entire 150 km Walnut
Gulch watershed, an experimental watershed operated by the Agricultural Research
Service (ARS) of the United States Department of Agriculture (USDA). While model
results were not encouraging, approximately half of the difference between simulated
and observed flows was attributed to rainfall sampling errors resulting from
inadequate rain gauge densities (Michaud and Sorooshian 1994b). Goodrich (1992)
21
claimed that the errors related to rainfall estimation and excess rainfall determination
are greater than the error introduced due to application of the kinematic wave model,
provided the criteria for kinematic wave modeling are satisfied and a sufficient
computational grid density is used. It was noted that the problem is more what to
route than how to route (Cordova and Rodriguez-Iturbe 1983).
2.1.1.1.6. Kinematic Shock
Shock formation has been another issue of concern related to application of
the kinematic wave model (Borah et al. 1980; Croley, II and Hunt 1981; Ponce 1991;
Woolhiser 1992). Kinematic shocks will form whenever there is an abrupt increase in
hydraulic resistance or a decrease in the slope or lateral inflow. While not directly
discussing issues of numerical dissipation and dispersion, Croley II and Hunt (1981)
demonstrated that multiple peaks could be obtained while using the Lax-Wendroff
scheme, cited as a popularly used method for numerically approximating the
kinematic wave formulation in many watershed rainfall-runoff models. Ponce (1991)
outlined the necessary conditions for the development of a kinematic shock. These
include a kinematic wave, a low base-to-peak flow ratio, a hydraulically wide and
sufficiently long channel, and a high Froude number flow. Ponce (1991) stated that
these conditions are rarely satisfied. Woolhiser (1992) noted that kinematic shock
development rarely occurs as a result of wave steepening, but rather is a consequence
of the geometric representation of the watershed. It was further noted that the
22
damping character of the finite difference schemes is irrelevant because the shocks
occupy such a limited portion of the solution domain.
2.1.1.1.7. Numerical Dissipation and Dispersion
Advanced distributed rainfall-runoff models typically involve simultaneously
solving the partial differential equations describing unsteady free surface flow as
overland or open channel flow as well as, possibly, the equation describing the
process of vertical flow in the unsaturated zone. Normally, these equations are
discretized and solved using a variety of numerical methods; often times finite
difference methods. The introduction of finite differences for the approximation of
the partial derivative terms in the descriptive equations can induce error due to
numerical dissipation and/or dispersion. In particular, a numerical scheme may do a
poor job of resolving the form of a propagating wave and/or its speed of propagation.
The linearized flow equations compose a linear hyperbolic system of the form
du k d u
---= A--
d t d x
(h^ ' 0 -h0)
u = A -
,UJ ~8 0 ,
23
where h0 is a constant and A is a constant diagonalizable matrix (Liggett 1975;
Vreugdenhil 1994). Equation 2.8 can be reduced to an uncoupled system of equations
of the form
<7 U p <7 U
d t d x
(2.9)
where D = SAS' is a diagonal matrix containing the eigenvalues of A, S'1 is the
matrix of eigenvectors of A, and U = Su (Thomas 1995). Because of the decoupling
of the linearized flow equations, the linear, homogeneous, constant coefficient one
way wave equation can also serve as the model linear equation for the de Saint
Venant equations.
The following discussion is from Thomas (1995). The solution procedure for
solving linear partial differential equations often involves expressing the solution as a
superposition of Fourier modes of the form
l+M
(2.10)
where a> is the frequency of the wave, /? is the wave number, A = is the
r
wavelength, and A is the amplitude. To satisfy the equations, (5 and a> are related by
24
an equation co = co{/5), the dispersion relation. For example, if Ae
|(
is
substituted into the linear one way wave equation,
(P,+c(px= 0
(2.11)
cp is a solution if co = -cf3, where c is the constant speed of propagation. The solution
of a PDE is said to be dissipative if the Fourier modes do not grow with time and at
least one mode decays, while the solution of a PDE is said to be dispersive if Fourier
modes of differing wave numbers propagate at different speeds {co" 0). In
consideration of the dispersion relation for the linear one way wave equation and the
above definitions, it is apparent that it is neither dissipative nor dispersive; the wave
moves with velocity c with no decay of the amplitude.
The dissipative and dispersive properties of a finite difference scheme can be
determined in a similar fashion. The solution is assumed to be expressible as a
superposition of discrete Fourier modes of the form
where 0 < /? Ax < n. The general wave is then inserted into the finite difference
scheme to obtain a discrete dispersion relation, where to, in general, is complex;
(2.12)
25
co = a(p) + /^(p). If this form for co is inserted into the general discrete Fourier mode,
it takes the form
if) I *Ai- -^-InA/
(2.13)
The following observations can be made from Equation 2.13:
1. the difference scheme is dissipative if b > 0 for some /?,
2. the difference scheme is nondissipative if b = 0 for all /?,
-a
3. the difference scheme is dispersive if a 0 and is nonlinear in /3.
This method was applied to briefly review the dissipative and dispersive
qualities of the forward in time, backward in space, FTBS, difference scheme for
numerically approximating the linear one way wave equation,
... "+>_ ,n c^ (rn n )
V k _(P k . yP k *P k-\)
Ax
(2.14)
the model linear equation for the kinematic wave approximation and the complete de
Saint-Venant equations. The scheme is first order accurate and stable as an initial
cAt
value-scheme provided 0 < R =---< 1. As a result, 0 < R < 1 is a necessary
Ax
26
condition for the stability of the difference scheme for an initial-boundary-value
problem. The stability of the scheme considered as a scheme for an initial-value
problem is determined by computing the symbol of the scheme, p, and bounding the
magnitude of p by one. The symbol of the FTBS scheme is given by
p (<Ã‚Â£) = 1 R + R cos Ã‚Â£ -iR sin Ã‚Â£.
(2.15)
The symbol is determined by taking the discrete Fourier transform of the scheme and
placing the transformed scheme into the form
~+'
9 k
P
(fV;.
(2.16)
where
9(4)
I
Ã‚Â£=-oo
(2.17)
The dissipative and dispersive properties of a scheme can be obtained through the
symbol of the difference scheme, p ). If a discrete Fourier mode is inserted into a
difference scheme, the discrete dispersion relation
27
e
= e
iatd -b A/
e
p(pnx)
(2.18)
is obtained. As a result, the dissipative and dispersive properties of a scheme are
determined by evaluating and analyzing
respectively. Taking into consideration Equations 2.19 and 2.20, Figures 2.1 and 2.2
are plots which provide graphical information about the dissipative and dispersive
qualities of the scheme. Examining Figures 2.1 and 2.2, several observations can be
made:
1. the FTBS scheme is both dissipative and dispersive,
(2.19)
and
tan a At -
(2.20)
Rep(/?Ax)
28
2. the error due to dissipation and dispersion depends upon the ratios and
R =
cAt
Ax
3. the error due to dissipation and dispersion is improved by decreasing Ax,
cAt
4. the error due to dissipation and dispersion is improved by choosing R =--------as
Ax
close to the stability limit, 1, as possible,
5. there is no dispersive error for R = 0.5; however, the scheme is very dissipative
for this choice of R.
29
Figure 2.2. Plot of error in propagation speed associated with the FTBS difference
scheme for R=0.1-0.9.
Several investigators have employed this approach to describe how well the decay
and propagation of various Fourier modes of a difference scheme match with those of
the linear PDE or to explain the behavior of the approximations in terms of the decay
and propagation of the Fourier modes (Liggett and Cunge 1975; Abbott 1979; Ponce
et al. 1979; Cunge et al. 1980; Holden and Stephenson 1995). For example, Liggett
and Cunge (1975) performed such an analysis on the popular Preissmann scheme for
numerically approximating the linearized flow equations to illustrate how accuracy is
affected by choosing different values for the parameters At, Ax, and the weighting
parameter, 9. Liggett and Cunge (1975) supplemented their linear analysis by
30
numerically approximating the linearized flow equations using the Preissmann
scheme; considering various values for and 6.
An alternative approach to investigate the dissipative and dispersive properties
of a finite difference (FD) scheme is to find the modified partial differential equation
(Warming and Hyett 1974; Thomas 1995). This method is a truncated approximate
approach. The numerical solution of a PDE using a FD scheme will come closer to
satisfying a different PDE, called the modified PDE. Thomas (1995) developed the
modified PDE associated with the Lax-Wendroff scheme, Crank-Nicolson scheme,
and the leapfrog scheme, among others, for numerically approximating the linear one
way wave equation, the model linear equation for the kinematic wave approximation
and the complete de Saint-Venant equations. The modified PDEs for the Lax-
Wendroff scheme, Crank-Nicolson scheme, and the leapfrog scheme are shown in
Equations 2.21 2.23, respectively.
Ax
(2.21)
cAx
(2 + J?2)*=0
h, + chx +
12
(2.22)
31
cAx /
h,+ch,+ (\-R2)h^=0
(2.23)
Equations 2.21 2.23 include derivatives up to and including the fourth order
derivatives. Equation 2.21, the modified PDE for the Lax-Wendroff scheme, clearly
emphasizes that the dispersive term is dominant over the dissipative term. Comparing
Equation 2.21 with Equations 2.22 and 2.23, it is also apparent that the Lax-Wendroff
scheme has an advantage over both the Crank-Nicolson scheme and the leapfrog
scheme; it has some higher order dissipation available to damp the highly dispersive
modes; whereas, the other two dispersive schemes do not. This highlights what has
been commented by past authors, that having some numerical dissipation in your
solution is not altogether undesirable. Thomas (1995) noted that it might be required
if a neutrally stable scheme is applied to solve a nonlinear PDE. There is a relation
between dissipation and stability (Thomas 1995).
The linear one way wave equation with a wave celerity of 1 was numerically
approximated by Skahill and Johnson (1999c) using several FD schemes on the
interval [0,l] subject to periodic boundary conditions and an initial profile of 1 for
x e [0.4,0.6] and zero elsewhere. For the experiment, values of Ax = 0.01 and At =
0.008 were uniformly used for all of the methods. For the popular Preissmann
scheme, the weighting coefficient, 0, was taken to be 0.7, an average value
recommended when the method is used in distributed watershed models (Grayson et
32
al. 1995; Smith et al. 1995). This simple experiment was performed to illustrate the
dissipative and dispersive qualities of the methods. With the understanding that
4. event-based, complex distributed watershed rainfall-runoff models typically
assign an initial condition of zero excess flow depth for overland flow plane elements,
5. calculated negative excess flow depths cause such a model to crash, and
6. rainfall input to such models can be variable in space,
the simple numerical experiment was also performed to suggest why the Preissmann
scheme remains the popular method of choice for approximating the routing
equations within complex distributed rainfall-runoff models. The results from the
experiment are shown in Figure 2.3 for t = 5 (625 time steps).
-0.4
ugure 2.3. Comparison of methods for linear advection experiment.
33
More than likely, the reason for the popularity of the four-point implicit Preissmann
scheme, used within many of the current distributed watershed rainfall-runoff models,
is its dissipative qualities for the weighting coefficient ranges frequently
recommended. Of notable interest, examining Figure 2.3, is that the FTBS scheme is
less dissipative than the popular Preissmann scheme for the average weighting
coefficient cited. From an implementation and maintenance perspective, the FTBS
scheme is clearly preferred over the popular Preissmann scheme. The explicit FTBS
scheme is trivial to implement; a simple one line statement. Furthermore, in contrast
to the popular implicit Preissmann scheme, it does not require iteration to obtain a
solution. Woolhiser (1975) noted that if convergence of the iterative scheme is slow,
the advantage of unconditional stability may be only apparent. Skahill and Johnson
(1999a) supplemented the linearized results from past investigators (Liggett and
Cunge 1975; Abbott 1979; Ponce et al. 1979; Cunge et al. 1980) and demonstrated in
a clear and practical way that grid dependence of FD solutions is not exclusive to the
kinematic wave model.
2.1.1.2. Runoff Generation
The work of Michaud (1992) underscored the importance of the comment that
the problem is more what to route than how to route (Cordova and Rodriguez-
34
Iturbe 1983). For the validation storms considered in the study by Michaud (1992), on
average, only 11% of the rain resulted in runoff at the main basin outlet.
2.1.1.2.1. Richards Equation
The increasing availability of inexpensive yet powerful computational
resources has made Richards equation an attractive choice for estimating storm
surface runoff.
The partial differential equation (Richards equation) used to describe transient
isothermal one-dimensional vertical water movement in non-swelling unsaturated soil
is obtained by combining the equation of continuity with the Darcy flux law. This
equation may be written in three standard forms, with either the soil water pressure
head (relative to the atmosphere), y/, or the volumetric water content, 9, as the
dependent variable. The introduction of the soil water retention function, which
provides a functional relationship between 9 and y/, allows one to convert from one
form of the equation to another. The three forms of the unsaturated flow equation
include the ^-based form, the Abased form and the mixed form. These
equations are written as
(2.24)
35
(2.25)
dO dK{0)
dt dz
d 0 dK(y/)
dt + dz
(2.26)
In Equations 2.24 2.26, c{y/) dO/dy/ is the specific water capacity, K is the
hydraulic conductivity of the soil, D{0) = K(0)/C(0) is the soil water diffusivity, z is
the depth oriented positively downward, t is time, and the soil is assumed to be
isotropic.
The highly nonlinear 0- y/ and K y/(ox K 9) relationships exclude the
possibility for analytical solutions except for special cases. As a result, numerical
methods are typically employed to solve the unsaturated flow equation. Discretization
often results in a nonlinear system of equations which are solved by applying a
linearization and/or iterative technique.
The ^-based form of the unsaturated flow equation is preferred for the reason
that it can be directly used for layered soils as well as for saturated and unsaturated
soils. Its main drawbacks include potentially poor mass balance errors and poor CPU
efficiency for very dry initial profiles. Celia et al. (1990) pointed out that the reason
/ \d y/ d 0
for poor mass balance resides in the time derivative term. While C t// and
F dt dt
36
are mathematically equivalent in the continuous sense, their discrete analogs are not
equivalent due to the highly nonlinear nature of the specific capacity term. While the
pressure based formulations may suffer from poor mass balance, the water content
formulation can be written in conservation form; and, as a result, its approximations
conserve mass well. In addition, since the terms K{0) and D(Q) are less nonlinear
than K[i(/), the 0-based form is more CPU efficient than the ^-based method (Hills
et al., 1989). The primary disadvantages of the #-based models include their inability
to model water flow in soils with saturated regions and their inability to be directly
applied to water flow in layered soils.
More recent developments for modeling transient water flow in unsaturated
soils include those presented by Hills et al. (1989), Celia et al. (1990), and Pan and
Wierenga (1995, 1997). Hills et al. (1989) developed a water content based algorithm
for modeling one-dimensional unsaturated water flow into layered soils. They
demonstrated that the #-based formulation is insensitive to dry initial conditions, in
contrast to ^z-based models. Their Abased model was from one to three orders of
magnitude faster than the yz-based models they considered when applied to
infdtration into very dry soils.
Celia et al. (1990) greatly improved the mass balance errors associated with
^z-based models by solving the mixed form of Richards equation treating y/ as the
dependent variable. The key to their algorithm, denoted as the ^-based modified
37
Picard method; which begins with a backward Euler in time approximation coupled
with the Picard iteration scheme, involves a truncated Taylor series expansion of
Qn+sm+{ with respect to y/ about the point y/n+l'm. The superscripts n and m represent
the time and iteration level, respectively.
Pan and Wierenga (1995) described a new approach for solving Richards
equation which combines the advantages of Abased and ^-based models. Their
method introduces a simple nonlinear transformation of y/, denoted by
P,=vl( \+Pv) as the dependent variable with the modified Picard method. They
applied their model to twelve test cases involving saturated, unsaturated, layered, and
uniform one-dimensional soil profiles, with pressure and flux type boundary
conditions. They demonstrated that their method provided improved results over
previously considered transformations and that it is a fast, numerically robust scheme
for solving Richards equation for variably saturated, heterogeneous media. In
addition, their method is easy to incorporate into already existing ^/-based codes. Pan
and Wierenga (1997) expanded their approach to solve problems in two dimensions.
The numerical solution of Equations 2.24 2.26 assumes that appropriate
constitutive relationships between #and ^and K and y/{ox K and 9) are available.
Singh (1997) summarized a large number of prediction models of soil hydraulic
properties. Sources are available which relate the unsaturated flow parameters of
prediction models of soil hydraulic properties to soil textural classification (Clapp and
38
Homberger 1978; Cosby et al. 1984; Stephens 1995; Leij et al., 1996; Carsel and
Parish 1988).
Skahill and Johnson (1999a) implemented the method described by Pan and
Wierenga (1995) for numerically solving Richards equation in one dimension as a
component of the advanced distributed rainfall-runoff model described in detail in the
following chapter. They evaluated its feasibility to support watershed scale surface
runoff computations. While attractive; for example, given its potential to incorporate
the ground water table for simulations on watersheds where the water table plays a
role in surface runoff production, it was not attractive from a computational
perspective. The numerical solution of Richards equation requires a fine spatial
discretization to maintain stability. With the exception of highly monitored locations,
uncertainty in unsaturated soil hydraulic properties does not warrant application of
Richards equation (Ogden and Saghafian 1997; Reed and Maidment 1998).
2.1.1.2.2. Green and Ampt Infiltration Model
The Green and Ampt (GA) infiltration model (Green and Ampt 1911) has
received a fair degree of attention in the engineering literature (Goldman et al. 1990;
James and Kim 1990; Julien et al. 1995; Task Committee on Hydrology Handbook
1996; Ogden and Saghafian 1997; Reed and Maidment 1998; Wang and Hjelmfelt
1998). The model assumes a pistonlike profile for volumetric water content. The
infiltration rate is governed by two parameters: the hydraulic conductivity behind the
39
wetting front and the wetting front capillary pressure head. Ogden and Saghafian
(1997) provided three reasons for the popularity of the GA model. First, it has
physical significance. Second, at / = 0 and t = oo, the GA equation correctly models
free vertical imbibition and saturated flow, respectively. Third, its parameters have
been related to soil texture class with reasonable success (McCuen et al. 1981; Rawls
et al. 1982; Rawls et al. 1983).
Reed and Maidment (1998) noted that the physical meaning of the Green and
Ampt wetting front capillary pressure head parameter has been the subject of debate
in the literature. Miyazaki (1993) commented that several investigations have been
conducted to verify the physical meaning of the Green and Ampt wetting front
capillary pressure head parameter. The development of Neuman (1976), also
presented in Reed and Maidment (1998) and Miyazaki (1993), illustrated the physical
meaning of the Green and Ampt wetting front capillary pressure head parameter. The
derivation of Neuman (1976) involved integrating the Darcy flux law from the land
surface to the assumed abrupt wetting front and observing that the hydraulic
conductivity behind the wetting front is constant. The following result was obtained,
(2.27)
40
which is a physically based definition for the Green and Ampt wetting front capillary
pressure head parameter, y/f. K, Ks, and y/, represent the unsaturated hydraulic
conductivity, saturated hydraulic conductivity, and the initial pressure head of the
soil, respectively. Chapter three includes additional references related to the GA
infiltration model. References related to other mechanisms relevant to runoff
generation and surface flow are also provided.
2.1.1.3. Watershed Time to Equilibrium
Surface runoff modeling is influenced by a watersheds time to equilibrium,
defined as the time at which information about the zero flow depth at the upstream
boundary arrives at the downstream boundary. A watersheds time to equilibrium is a
function of its physical and geometrical properties as well as storm characteristics.
Equilibrium flow conditions occur when the duration of runoff production exceeds
the time of equilibrium. At equilibrium, assuming a steady rate of runoff production,
the peak discharge is constant and it equals the contributing area times the rate of
runoff production. Furthermore, at equilibrium, surface runoff increases linearly with
the rate of runoff production, justifying the application of linear methods to model the
rainfall-runoff process. Michaud (1992) listed several factors that would be expected
to influence whether or not equilibrium conditions occur. These include watershed
size, geometry and steepness, storm duration, areal extent of the storm, and
41
contributing area. Molnar and Julien (2000) noted that grid cell size would influence
distributed surface runoff model computations if partial equilibrium conditions exist.
2.1.2. Existing Models
Advanced distributed watershed scale surface runoff models were described
by Wang and Hjelmfelt (1998), Orlandini and Rosso (1996), Garrote and Bras (1995),
Julien et al. (1995), Vieux and Gauer (1994), Grayson et al. (1992a), Goodrich et al.
(1991), James and Kim (1990), Woolhiser et al. (1990), Beven et al. (1987), Field and
Williams (1987), Abbott et al. (1986), Jayawardena and White (1979), Ross et al.
(1977), and Singh and Woolhiser (1976), among others.
Wang and Hjelmfelt (1998) described a DEM based distributed rainfall-runoff
model that was developed to route overland flows from flat agricultural watersheds.
The model numerically approximates the one-dimensional diffusive wave
approximation using the MacCormack scheme. Overland flow is routed as a
cascading system, starting from the most upstream cells and ending with the
downstream cell. Infiltration rates are computed using the Green and Ampt equation.
The model does not include channel flow. Skahill and Johnson (1999a) implemented
the MacCormack method to numerically approximate the conservative form of the
kinematic wave approximation for overland flow in the overland flow component of
the advanced distributed rainfall-runoff model described in detail in the following
chapter. A uniform excess rainfall rate of one inch per hour was uniformly applied for
42
one hour to the Buffalo Creek watershed, a basin with a drainage area of
approximately 130 square kilometers located in the Rocky Mountain foothills
approximately 35 miles southwest of Denver, Colorado. Channel flow was not
simulated and a uniform overland Mannings n value of 0.03 was used for the
computation. The results obtained from the model using the MacCormack scheme,
and also the nonlinear forward in time, backward in space scheme (FTBS), are shown
in Figure 2.4. The dispersive character of the MacCormack scheme, which reduces to
the popular Lax-Wendroff scheme when it is applied to a linear problem, is easily
observed in the one hour unit hydrograph. The second order accurate MacCormack
method failed to complete the simulation. While the method may be applicable for
very flat sloped watersheds, the results of Skahill and Johnson (1999a) indicated that
the method would not apply in highly variable and steeply sloped terrain.
43
600
Figure 2.4. Comparison of one hour unit hydrographs obtained using two FD
schemes.
Orlandini and Rosso (1996) developed a simple and computationally
inexpensive network routing model based on the Muskingum-Cunge method. Surface
runoff is routed downstream, cell by cell, from the uppermost cell to the basin outlet.
Excess rainfall is determined externally, and input to the routing model at a one-hour
time step. Garrote and Bras (1995) described a distributed rainfall-runoff model that
was developed with the intent of formulating a parameterization of the runoff
generation process more physically based than that of conceptual models. It was
designed to support real-time flood forecasting in medium and large basins. The
runoff generation component of the model supported the infiltration-excess runoff
44
mechanism as well as exfiltration of subsurface flow. The distributed convolution
equation was utilized to determine the hydrograph at the main basin outlet.
Julien et al. (1995) described a two-dimensional physically-based surface
runoff model compatible with raster weather radar rainfall estimates that was
expressly designed for the purpose of fully utilizing raster GIS. Excess rainfall,
calculated using the Green and Ampt equation, is routed as overland flow using a
two-dimensional diffusive wave explicit method that does not route water in both the
x and y directions simultaneously, but rather sequentially in each direction (Ogden
1992). Excess rainfall is not routed as overland flow on the flood plain portion of
raster grid cells that contain channel segments, defined as drainage cells. Rather,
once a predefined retention depth is satisfied for a drainage cell, all incoming
overland flow from adjacent raster grid cells is passed directly to the channel. The
algorithmic approach for channel flow routing within the model is analogous to
overland flow; however, only one space dimension is considered.
Vieux and Gauer (1994) described a storm surface runoff model that was
internally integrated within a raster GIS. The model estimates infiltration rates for
each raster grid cell or sub-basin using the Green and Ampt equation. The one-
dimensional kinematic wave approximation is used to model overland flow and
channel flows. It is numerically approximated using the finite element method.
Grayson et al. (1992a) described an event-based runoff model that simulates
both the infiltration-excess and saturation-excess runoff mechanisms. The model uses
45
a contour-based method of terrain analysis. The Preissmann scheme in combination
with Newtons method is used to solve the kinematic wave approximation for the
routing of surface flow. Goodrich et al. (1991) described a two-dimensional
kinematic wave surface runoff model that operates on a triangular irregular network
DEM, and utilizes the finite element method. Further research was required to
incorporate infiltration into the model.
James and Kim (1990) developed a distributed watershed model to simulate a
single event storm utilizing remotely sensed data. The Green and Ampt equation was
used to calculate excess rainfall. Overland flow was modeled by numerically solving
the diffusive wave approximation of the two-dimensional de St. Venant equations
using an alternating direction implicit scheme. Channel flow routing was modeled by
numerically solving the one-dimensional de St. Venant equations using a four-point
implicit scheme.
Woolhiser et al. (1990) described an event-based, kinematic, infiltration-
excess, distributed rainfall-runoff model. The model routes excess rainfall as overland
or open channel flow using the popular implicit Preissmann scheme in combination
with Newtons method. Within the model, overland flow elements completely cover
the basin. As a result, rainfall is not directly applied to channels. Infiltration rates are
computed using the model described by Smith and Parlange (1978).
Beven et al. (1987) described a physically-based rainfall-runoff model that
solves the two-dimensional ^/-based form of Richards equation to model subsurface
46
flow. The kinematic wave approximation, numerically solved using the Preissmann
scheme in combination with iteration, is used to model surface flow.
Abbott et al. (1986) describe the mathematical and computational structure of
the European Hydrological System (SHE), an advanced, physically-based,
distributed, catchment model. Their model includes the primary components of the
land phase of the hydrologic cycle snowmelt, interception, evapotranspiration,
overland and channel flow routing, and unsaturated and saturated subsurface flow.
They implemented their model in a modular format, allowing for the flexibility not
only to easily update already existing components, but also to add new components.
The diffusive wave approximation of the two-dimensional de St. Venant equations
which describe the physics of gradually-varied overland flow is numerically solved
using an explicit finite difference scheme described by Preissmann and Zaoui (1979).
The approach for channel flow routing within the model is analogous to overland
flow; however, only one space dimension is considered and they use an implicit finite
difference described by Preissmann and Zaoui (1979). The one-dimensional y/-
based form of Richards equation is used to model unsaturated subsurface flow. The
equation is numerically solved using an implicit finite difference scheme.
All of the advanced distributed models cited above must be calibrated
manually. In addition, none were designed to provide a range of system responses
from the calibrated model. Melching (1995) commented that estimation of the
47
reliability of model output is important to hydrologic engineering and water-resources
planning.
2.1.3. Geographic Information Systems
and Hydrologic Modeling
A few of the more current advanced distributed models cited above were
expressly designed to capitalize on relatively recent advances in GIS technology and
the increasing availability of watershed databases in raster or grid cell GIS format
(Vieux and Gauer 1994; Julien et al. 1995). While discussing the application of GIS
in watershed modeling, Singh (1995) noted that the watershed and hydrology are
inherently spatial. The physical characteristics of the watershed; such as soils, land
use, and topography vary spatially. Furthermore, distributed watershed models can
directly use these data. A GIS serves as a spatial data management and processing
environment; providing the means for mapping and evaluation at multiple scales.
Johnson (1989) demonstrated the range of application of GIS functions for
hydrologic modeling with the MAPHYD package. The work by Johnson (1989)
showed that compatible data sets on terrain and radar-rainfall data could be fully
integrated into a watershed modeling system. Michaud (1992) commented on the
possibility that, in the future, advanced distributed rainfall-runoff models would be
applied to large numbers of watersheds on the basis of digitized data using GIS
technology. DeVantier and Feldman (1993) presented an overview of hydrologic
48
modeling with GIS. The different types of terrain data and hydrologic models were
analyzed. In support of flash flood operations, an assessment by Johnson (1998 and
2000) noted that GIS tools provide a means for integrating into a forecaster-usable
interface three dimensions that forecasters must address in recognizing and
forecasting flash flood threats. The recent work by DeBarry et al. (1999) provided an
excellent summary of GIS and distributed watershed modeling; discussing GIS data
types, map projections, data commonly required for hydrologic analysis using GIS
techniques, available data sources, limitations of available data, and the integration of
watershed hydrological analysis software and GIS techniques.
2.1.3.1. Digital Elevation Models and Hydrologic Modeling
Topography plays a critical role in hydrologic modeling. Landscape
topography can be digitized into an array of elevation values called a Digital
Elevation Model (DEM). DEMs are commonly stored in one of three data structures:
grids, contours, or triangulated irregular networks (TINs). Moore et al. (1991)
reviewed the advantages and disadvantages of each of these three data structures.
DeBarry et al. (1999) noted that automated extraction of surface drainage, channel
networks, drainage divides and other hydrography data from DEMs has established
itself in GIS over the past decade. This was attributed to the increasing availability
of DEMs and software tools that derive this data from DEMs (Jenson and Domingue,
1988; Mark, 1988; Moore at al., 1991; Martz and Garbrecht, 1992 and 1993). It was
49
further noted that automated techniques are faster and provide more precise and
reproducible measurements than traditional manual techniques applied to topographic
maps.
The advanced distributed rainfall-runoff model described in detail in the
following chapter relies on the topographic parameterization model TOPAZ
(Garbrecht and Martz 1999). This is primarily due to the tabular output files on
watershed drainage topology the TOPAZ model provides; which can be easily
modified and used as input data to support hydrologic modeling (DeBarry et al.
1999). The DEM processing in TOPAZ is based on the D-8 method and the critical
source area (CSA) concept. The D-8 method (Fairfield and Leymarie, 1991)
identifies the steepest downslope flow path between each cell of a raster DEM and its
eight neighbors, and defines this path as the only flow path leaving the raster cell. The
CSA concept defines the channels draining the landscape as those raster cells that
have an upstream drainage area greater than a threshold drainage area. The CSA
value defines a minimum drainage area below which a permanent channel is
maintained (Mark, 1984; Martz and Garbrecht, 1992). The CSA concept controls the
watershed segmentation and all resulting spatial and topologic drainage network and
subcatchment characteristics.
DeBarry et al. (1999) commented that the D-8 method recently has been
criticized on the grounds that it permits flow only in one direction away from a raster
grid cell. This approach fails to adequately represent divergent flow over convex
50
slopes (Freeman, 1991; Quinn and others, 1991; Costa-Cabral and Burges, 1994) and
can lead to a bias in flow path orientation (Fairfield and Leymarie, 1991). Although
the multiple flow direction algorithm seems to give superior results in the head-water
region of a source channel, DeBarry et al. (1999) noted that a single flow-direction
algorithm is superior in zones of convergent flow and along well-defined valleys
(Freeman, 1991; Quinn at al., 1991). For overland flow analysis on hillslopes the
multiple flow approach may be more appropriate. If the primary objective is the
delineation of the drainage network for large drainage areas with well-developed
channels, the use of a single flow-direction algorithm was stated by DeBarry et al.
(1999) to seem more appropriate.
2.1.4. Calibration
Complex distributed models are typically calibrated against a historical record
of observed events to compensate for uncertainties in data, model parameters, and
model structure (Binley et al. 1991). Ironically, the capability to incorporate
enormous quantities of spatially-varied data is the clear weakness for distributed-
parameter models. The large number of model parameters that may vary from cell to
cell, parameter interaction, and relatively significant computational requirements
make the calibration of a complex distributed model a difficult task. This was clearly
illustrated in the work by Michaud (1992) and Michaud and Sorooshian (1994a). For
example, Michaud (1992) and Michaud and Sorooshian (1994a) were not confident
51
that the optimal parameter set was identified for their application of an advanced
distributed rainfall-runoff model; even after subdividing the 150 km basin/study site
into eight regions for the calibration of the hydraulic conductivity for planes.
The work by Michaud (1992) and Michaud and Sorooshian (1994a) also
reinforced the observation made by Beven (1989) that physical reasoning may be of
little help in calibrating a physically-based distributed model. Binley et al. (1991)
stated two problems related to relying upon field measurements alone to calibrate a
physically-based model. The level of discretization required to maintain an accurate
numerical approximation with the descriptive partial differential equation(s) may
result in far too many parameters to determine by experiment or the field test may
very well not exist. While optimistic about the future for distributed rainfall-runoff
models, Bras (1999) noted that their calibration remains a challenge.
Several studies have suggested or observed that it is unlikely there is a unique
optimal parameter set for a specific complex distributed model structure (Binley et al.
1991; Beven and Binley 1992; Beven 1995; Beven et al. 1995; Freer et al. 1996;
Franks et al. 1998). Rather, there may be many parameter sets within a model
structure that equally represent the rainfall-runoff process in terms of some
subjectively chosen performance measure. Given the lack of a rigorous basis for
differentiating between the numerous acceptable parameter sets, which may lie in
different regions in the parameter space, Beven and Binley (1992) developed the
52
Generalized Likely Uncertainty Estimation (GLUE) procedure as a method for model
calibration and uncertainty estimation.
2.1.4.1. GLUE
The GLUE methodology of Beven and Binley (1992) does not acknowledge a
unique optimal parameter set. Instead, Monte Carlo simulation, a likelihood measure
or set of likelihood measures, and an acceptance/rejection criteria are used to
generate, evaluate, and accept or reject parameter sets as simulators of a basins
response. Beven and Binley (1992) noted that the term likelihood is used in a general
sense, and not in the restricted sense of maximum likelihood theory. The likelihood
measure is associated with a particular set of parameter values within a given model
structure. Parameter interaction and uncertainties in the input and boundary data are
implicitly reflected in the likelihood measure.
The results of the GLUE methodology of Beven and Binley (1992) are
dependent upon the choice of the likelihood measure and the rejection criteria. The
rejected models are assumed nonbehavioral and their likelihood measures are set to
zero; thereby removing them from the subsequent analysis. The likelihood measures
for the accepted parameter sets are utilized to estimate the uncertainty of the models
predictions. Beven and Binley (1992) commented that all simulations with a
likelihood measure significantly greater than zero are retained for consideration.
Freer et al. (1996) stated 0.3 to be a low model efficiency rejection criterion. They
53
commented that this would ensure wide enough uncertainty bounds to envelop most
of the observed responses during the calibration periods. The likelihood measures
associated with the retained models can be rescaled. Beven and Binley (1992) noted
that the traditional calibration search for a unique optimal parameter set is an extreme
case of the GLUE procedure in which the optimal model has a likelihood of one and
all other likelihood measures are set to zero.
Beven and Binley (1992) listed some example likelihood measures. Using the
notation of Beven and Binley (1992), two commonly used likelihood measures are
(2.28)
(2.29)
where Z,(<9,|f), erf, and erf represent the likelihood measure of the ith model
conditioned on the observations, error variance of the /th model, and observed
variance for the period under consideration, respectively. The likelihood measure
defined in Equation 2.28 is equivalent to the Nash and Sutcliffe (1970) efficiency
criterion when n= 1. Freer et al. (1996) and Franks et al. (1998) noted that higher
values of n accentuate the weight given to the better simulations; and this approach is
best suited to finding optimal parameter sets. Freer et al. (1996) noted that the
likelihood measure defined in Equation 2.29 has the characteristic that in applying
54
Bayes equation for updating likelihood measures, the error variance for each period
of data is given equal weight. Beven et al. (1995) noted that an important part of the
GLUE procedure is that it allows likelihoods to be defined according to different
types of calibration data.
Beven (1998) discussed four alternative methods for updating likelihood
measures as additional data becomes available
lM\y)=L(e,)L(e,\Y) (2.30)
LP(e,\y)=K(e,)+i(e,\r) (2.3 u
I(,|F)=maxji(9,),l(e,|y)} (2.32)
LF(e,\Y)=mm{lXe\L(e,\Y)} (2.33)
Lp(Gl | f) is the posterior likelihood for the /th model, L0(6,) is the prior
likelihood, and L(0l \ Y) is the likelihood calculated for the current data set Y. The use
of Bayes equation, Equation 2.30, has received the most attention. Equation 2.32,
defined as a fuzzy union, emphasizes the best performance of each model over all the
measures considered, while Equation 2.33, defined as a fuzzy intersection,
emphasizes the worst performance.
Beven and Binley (1992), Beven et al. (1995), Melching (1995), and Freer et
al. (1996) outline the requirements for the GLUE procedure. Applications of the
55
GLUE methodology of Beven and Binley (1992) to rainfall-runoff modeling include
the work of Beven and Binley (1992), Beven et al. (1995), Freer et al. (1996), and
Franks et al. (1998). Computational demands are high for the GLUE procedure.
However, computational constraints have also limited the success of manually
calibrating a physically-based distributed rainfall-runoff model (Michaud 1992;
Michaud and Sorooshian 1994a). Michaud (1992) expressed the need for an
automated calibration technique that would ease the application of an advanced
distributed rainfall-runoff model.
Beven et al. (1995) noted that a narrowing of the uncertainty bounds indicates
an improvement in the model calibration; whereas, a widening of the uncertainty
bounds indicates a reduction in the predictive capabilities of the model. Melching
(1995) identified the main shortcoming of the GLUE procedure to be the criteria for
acceptance and rejection of the models. If the acceptance criteria is set too high, then
there will only be a small number of accepted parameter sets, and a narrow,
unrealistic range for system responses. On the other hand, if the acceptance criteria is
set too low, then almost all of the parameter sets will be accepted, the predictive
range will be wide, and the model will be considered to be of little value.
For purposes of illustration, the GLUE methodology of Beven and Binley
(1992) was used to calibrate the advanced distributed rainfall-runoff model described
in detail in the following chapter for the June 1, 1991 event that occurred on the
-y
Ralston Creek basin, a 225 km watershed located within the UDFCD. The results
56
provided in Figures 2.5 2.8 illustrate the best likelihood measure, worst
likelihood measure, and the model predictive range for two different rejection criteria.
The Nash and Sutcliffe (1970) efficiency criterion was used as the likelihood
measure. It was computed for each model by comparing observed discharges with
computed discharges at 250 different time steps. The results are based on 250 Monte
Carlo simulations.
The attraction of the GLUE procedure is its simplicity and ability to be
applied to models of arbitrary complexity. Refsgaard et al. (1996) commented that the
GLUE procedure is the most comprehensive approach for assessing the uncertainty of
model predictions, and encouraged further development and application.
57
90
sb oo <> O rn
CN
time (MDT)
Figure 2.5. Best model obtained from 250 Monte Carlo Simulations.
Figure 2.6. Worst model obtained from 250 Monte Carlo Simulations.
58
Figure 2.7. Ninety-five percent uncertainty bounds obtained from the model for the
June 1, 1991 storm on the Ralston Creek basin. Model rejection criteria set at 0.3
using the Nash and Sutcliffe efficiency criterion as the likelihood measure.
Figure 2.8. Ninety-five percent uncertainty bounds obtained from the model for the
June 1, 1991 storm on the Ralston Creek basin. Model rejection criteria set at 0.8
using the Nash and Sutcliffe efficiency criterion as the likelihood measure.
59
2.1.5. Modeling Studies
Michaud and Sorooshian (1994a) addressed some of the capabilities and
limitations of the complex distributed model approach. The study had two major
objectives:
i. Compare a simple distributed model with a complex distributed model under
hydrologic and data conditions that are relevant to semiarid watersheds which
have been instrumented for flash flood forecasting, and
ii. Compare simulations produced with and without calibration.
The study site was the entire 150 km Walnut Gulch Experimental Watershed
operated by the ARS of the USDA. It was primarily chosen due to the high-quality,
long-term rainfall, runoff, topographic, and soils data available for the basin; thus
providing the basis for a meaningful validation of the calibrated model. Horton
overland flow was assumed to be the dominant surface runoff mechanism at the
Walnut Gulch watershed. That is, lateral movement of water in the subsurface and
saturation of the soil surface due to rising water tables was assumed to either not
occur or to be relatively unimportant. Thirty storms were used for model calibration
and validation. The six largest runoff events during the period 1974-1976 were used
for manual calibration. Twenty of the remaining twenty-four storms were those which
caused the twenty largest runoff events during the period 1957-1973 and 1977. The
four remaining storms were selected because the amount of runoff was unusually
60
small for the amount of rain. The study focused on simulation accuracy under data
constraints typical of ALERT (Automated Local Evaluation in Real Time)
watersheds. As a result, only 8 of the approximate 93 Walnut Gulch recording rain
gauges were used for the study; resulting in a distance of about 4.5 km between
gauges; assuming they are equally spaced. Radar-estimated rainfall data was not
considered in the study because it was not available. In addition, Michaud and
Sorooshian (1994a) comment that raingauges were more relevant to the existing
ALERT technology. The complex distributed model used in the study was a research
version of the KINEROS model (Goodrich 1990). KINEROS (Woolhiser et al.,
1990), an acronym for KINematic runoff and EROSion model, is an event-based,
infiltration-excess, kinematic, distributed rainfall-runoff model compatible with
spatially and temporally varied rainfall; developed specifically for semiarid regions.
The Soil Conservation Service (SCS) model was chosen for the simple model in the
study because it is widely used, has minimal data requirements, and can be applied to
ungauged basins. Both a spatially lumped and spatially distributed version of the SCS
model was considered in the study. For the distributed SCS model, rainfall input,
antecedent moisture conditions, and curve numbers were spatially distributed. In
addition, channel routing was accounted for through translation of the subarea
hydrographs. Before any routing was performed, channel losses were estimated. The
KINEROS model consisted of approximately 60 plane elements and 41 channel
elements. The watershed subdivision was performed manually. For the uncalibrated
61
KINEROS model, published guidelines and field observations were used to estimate
roughness values. Infiltration parameters were estimated based on soil textural
classification data from a county soil survey. For the uncalibrated SCS model,
Thiessen polygons were used to subdivide the basin, with one subarea per rain gauge.
Standard SCS techniques were used to estimate the curve numbers and unit
hydrographs. The hydrograph translation velocity was estimated using Mannings
equation and each channels loss parameter was estimated to be the product of the
average channel width (obtained from field surveys), channel saturated hydraulic
conductivity (estimated from the literature), and average flow duration. Manual
calibration of the KINEROS model was performed from six storms over a three-year
period. The calibration was based on a comparison of the simulated and observed
flows. For simplicity, the basin was subdivided into 8 areas to facilitate the
calibration of the hydraulic conductivity. In addition, five other parameters were
optimized:
i. Hydraulic conductivity for the main channels,
ii. n for the main channels,
iii. n for the lower order channels,
iv. n for the overland flow planes, and
v. the coefficient of variation of the hydraulic conductivity for the planes.
The model CREAMS (Knisel 1980) was used to determine initial soil moisture
profiles at each raingauge, utilizing daily rainfall data from the first of the year to the
62
start of each storm. Parameter estimates for the model CREAMS were based on
published guidelines. Initial soil profiles for the model subareas were obtained by
interpolating the values obtained at the eight gauges. The primary outcome of the
manual calibration of the KINEROS model was a correction for the time to peak. The
calibration did not provide for an improvement for volumes or the peak discharges.
The authors were not confident that the optimal parameter set was identified due to
computational and data constraints, and significant parameter interaction. The
calibration procedure for the SCS model was automated utilizing a search strategy
and an objective function involving squared errors in peak flow, volume, and time to
peak. The curve number and the time of concentration were calibrated for the lumped
SCS model, while the channel loss parameter, the hydrograph translation velocity, a
curve number multiplier, a time of concentration multiplier, and the difference
between curve numbers for dry and wet antecedent moisture conditions were
calibrated for the distributed SCS model. The automated calibration procedure for the
two SCS model configurations converged rapidly to a set of parameter values that
resulted in a local minimum for the objective function. The models were validated for
the 24 storms by comparing simulated and observed peak flow, runoff volume, time
to peak, and the ratio of peak flow to volume. The major observations and
conclusions from the comparison of the two modeling approaches are summarized
below:
63
i. The application of KINEROS on the Walnut Gulch Experimental Watershed
under ALERT-like data constraints resulted in poor simulations of peak
discharge and runoff volume. If Mannings n values were calibrated on
historic data, the model provided good simulations of time to peak.
ii. When no calibration was performed with the KINEROS model; that is,
parameter values were based on field surveys or values provided in the
literature; the simulations did a better job of simulating peak discharge;
however, they were biased for time to peak (they were consistently too slow).
iii. The calibrated spatially-distributed SCS model was slightly more accurate
than the KINEROS model. When no calibration was performed; the
KINEROS model was more accurate than the distributed SCS model. The
spatially lumped SCS model did not perform well at all, primarily due to the
requirement of averaging the rainfall over the entire basin.
iv. In comparison to the simple model, the complex spatially-distributed
KINEROS model required one order of magnitude user time for
implementation and two orders of magnitude more computer time; with this
additional effort providing no increase in accuracy.
v. The complex distributed KINEROS model was not recommended for ALERT
forecasting since the simpler alternative appeared to be just as accurate.
vi. Of the various factors which influence the accuracy of a simulation; for
example, input data, auxiliary conditions, model assumptions, parameter
64
estimates, model spatial and temporal resolution; the investigators noted that
the estimation of rainfall and infiltration were the two primary sources
contributing to the poor model simulations obtained from the KINEROS
model.
Despite the disappointing performance of the complex distributed model, the authors
stated that it would be shortsighted to discontinue research related to these types of
models, and that there is ample opportunity for continued work in this area; such as,
identification of dominant processes. Michaud (1992) encouraged the consideration
of a similar study for thunderstorm-generated floods on the east slope of the Rocky
Mountains.
Michaud and Sorooshian (1994b) investigated the effect of rainfall-sampling
y
errors on complex distributed rainfall-runoff simulations for the entire 150 km
Walnut Gulch Experimental Watershed. The study was motivated by the desire to
determine whether the poor results obtained from the complex distributed KINEROS
model on Walnut Gulch, under data constraints typical of ALERT watersheds, were
due to inadequate raingage densities. Three rainfall sampling schemes were
considered in the study:
1) a raingage density of approximately one gage per 20 km (a raingage density
typical of ALERT networks and the density used in the study of Michaud and
Sorooshian (1994a)),
'y
2) a raingage density of approximately one gage per 2.5 km and
65
3) a rainfall data spatial resolution of 4km x 4km.
Simulations obtained from the first rainfall sampling scheme resulted in reduced
runoff that, on the average, represented 58% of the observed peak flow.
Approximately half of the difference between simulated and observed flows was
attributed to rainfall sampling errors resulting from inadequate raingage densities.
2.1.5.1. Predictive Uncertainty
Melching (1995) described sources of uncertainty and methods for estimating
uncertainty, measures of output reliability, reliability analysis methods, a comparison
of these methods, and also reviewed how one can identify the key sources of
uncertainty using reliability analysis. There are four types of uncertainty that a
hydrologic modeler must contend with:
i. data,
ii. model parameters,
iii. model structure, and
iv. natural randomness.
Melching (1995) commented that at the time of writing the article, there had
not been a study which performed a complete reliability analysis considering
uncertainties in the rainfall data and all model parameters. Melching (1995)
commented that the CDF and PDF for various output variables of interest; such as,
peak discharge, runoff volume, and time to peak, can be used to assess model
66
reliability. For example, it was commented that an acceptability criterion could be
that 90% of the possible output values be within 30% of the predicted value.
Reliability analysis methods described included Monte Carlo simulation (MCS) and
Latin Hypercube simulation (LHS), among others. It was commented that MCS may
be the only method which can approximate the CDF and PDF for an output variable
of interest if complex system relationships exist. In addition, it was commented that
when the computing power is available, there can be no strong argument against the
use of MCS. Melching (1995) commented that reliability analysis can be used to
determine the key sources of uncertainty affecting the reliability of the models
output. As a result, reliability analysis can be used to guide data collection activities.
Binley et al. (1991) investigated the predictive uncertainty of the Institute of
Hydrology Distributed Model, IHDM, using two reliability methods, the method of
Rosenblueth and Monte Carlo simulation, on a 3.9 km basin in central Wales. Based
on a sensitivity analysis, four model parameters, saturated hydraulic conductivity,
saturated moisture content, initial capillary potential of the soil, and an overland flow
roughness coefficient, were calibrated by comparing observed and simulated
discharge hydrographs. The parameters were assumed to be uniform in space.
Statistics for the model parameters, subsequently used in the estimation of the
models reliability, were obtained from the calibration exercise; particularly, mean
and standard deviations. Their results suggested that the uncertainty bounds for
physically-based models can be quite wide, even when the parameters have been
67
constrained by the calibration process. It was also shown that the method of
Rosenblueth provides for a reasonable first estimate of the limits of predictive
uncertainty based on a small number of simulations.
Additional studies include those of Garen and Burges (1981), Goldman et al.
(1990), Melching et al. (1991), Beven and Binley (1992), Freer et al. (1996), and
Franks et al. (1998).
2.2. Weather Radar
The modernization of the U.S. National Weather Service (NWS) has the goal
of providing improved forecasts and more reliable detection and prediction of severe
weather and flooding. Two major components of the modernization include the
installment of the weather radar system WSR-88D (Doppler Weather Surveillance
Radar 1988) and AWIPS-90 (Advanced Weather Interactive Processing System for
the 1990s) workstations in NWS Weather Forecast Offices (WFOs).
Several authors have commented on the new opportunities the NWS
modernization program will provide for operational hydrology. Hudlow (1988)
described the hydrologic forecasting program of the National Weather Service (NWS)
in the United States, focusing on current and future forecasting technology. The
principal goal of the research and development effort of the NWS Hydrologic
Program is to produce and implement improved forecasting technology at the River
Forecast Centers (RFC) and Warning and Forecast Offices (WFO) nationwide.
68
Hudlow (1988) provided a schematic of a future integrated hydrometeorological
forecast system that would integrate information and provide feedback from models
running at the national, major river basin, and local WFO levels. Processes
incorporated in event-oriented flash flood models that would run at a WFO included a
local Quantitative Precipitation Forecast (QPF) model, precipitation, catchment input
variables, a watershed model, channel inflow, and a channel routing model.
Model/data linkages, soil moisture information, and atmospheric variables and
products would come from the RFCs and national models. Hudlow (1988)
emphasized that the NWS must improve quantitatively its ability to predict short scale
hydrologic events, such as flash floods.
A recent assessment (Johnson, 1998 and 2000) of four alternate flash flood
warning procedures using radar-estimated rainfall data was motivated by the
observation that a significant component of the modernized forecasting environment
is new capabilities for hydrologic tracking and forecasting. The study emphasized that
Weather Forecast Office (WFO) forecasters must address three fundamental
dimensions in recognizing and forecasting flash flood threats, and issuing appropriate
watch/waming messages. These dimensions are the driving hydrometeorological
(hydromet) phenomenon, the land surface response, and the flood impact features.
The land surface response dimension includes watersheds and the stream network,
soils and vegetation, and land use/urbanization. Impact features include rivers, towns,
developed areas, bridges and river crossings, and flood plain developments. Johnson
69
(1998 and 2000) noted that Geographic Information Systems (GIS) tools for spatial
data management and processing provide a means for integrating the three
dimensions into a forecaster-usable interface. The study by Johnson (1998 and 2000)
emphasized that hydrometeorological data alone is not enough by itself; these data
products must be coupled with the land surface response and flood impact feature
dimensions to obtain usable flash flood forecast products. It was concluded that for
the Colorado Front Range, where storms have a mean size of thirty square kilometers,
the level of spatial detail for assessing flood runoff and potential impacts is currently
inadequate. The study concluded with the recommendation that an advanced
distributed approach should be applied for the forecast region; coupled with a method
for assessing the reliability of the models predictions.
2.2.1. Weather Surveillance Radar-1988 Doppler
Rainfall Algorithm
Fulton et al. (1998) described the WSR-88D rainfall algorithm. The article
updated previous papers describing the original Precipitation Processing System
(PPS) design. It was specific to the WSR-88D software release, Build 9. The
authors provided a comprehensive description of the five principal algorithmic
components of the PPS, the two external support functions, current system
limitations, an overview of stage II and stage III processing, and opportunities for
future research.
70
The Next Generation Weather Radar (NEXRAD) program is a federal
triagency program of the National Weather Service (NWS; Department of
Commerce), Federal Aviation Administration (Department of Transportation), and
Air Force Air Weather Service and Naval Oceanography Command (Department of
Defense). The program has resulted in the installment of over 160 S-band (10 cm)
Weather Surveillance Radar-1988 Doppler (WSR-88D) radars across the United
States. The NEXRAD program has been a major component of the modernization of
the NWS.
The PPS consists of five components and two external support functions. The
five subalgorithms include preprocessing, rate, accumulation, adjustment, and
products, while the two external support functions, precipitation detection and rain
gauge data acquisition, operate independently of the PPS. The two external support
functions provide additional input data to the main PPS algorithm.
The two main functions of the precipitation detection function (PDF) are to
determine the scanning mode of the radar antenna and the processing mode of the
PPS. The PDF will automatically switch the scanning mode of the radar antenna from
clear air to precipitation mode whenever non-ground-clutter reflectivity echoes
exceed predefined intensity and areal coverage thresholds for the four lowest
elevation sweeps. In clear air mode the volume scan interval is ten minutes;
whereas, in precipitation mode, the volume scan interval is five or six minutes.
Different threshold values are used by the PDF to determine the processing mode of
71
the PPS. These values are typically lower than those used by the PDF to determine
the scanning mode of the radar antenna. They relate to the lower bound on resolvable
rainfall in the PPS algorithm. If the threshold values are exceeded, the PDF instructs
the PPS algorithm to begin accumulating rainfall regardless of whether the scanning
mode is clear air or precipitation. The end of a rainfall event is defined to occur
when there is no detectable precipitation by the PDF over a one-hour period. At this
time, the storm total precipitation product is set to zero.
The rain gauge data acquisition function (RGDAF) is the second PPS support
function within the radar product generator (RPG) computer. The RGDAF receives
real-time rain gauge data from an external gauge data support (GDS) computer and
places them in a database within the RPG for use by the PPS. The RGDAF is inactive
until the PDF detects rain. When rain is detected by the PDF, the RGDAF sends a
message to the GDS, requesting real time gauge reports as they become available at
the GDS. When rain is no longer detected by the PDF, the RGDAF sends a message
to the GDS requesting it to stop sending rain gauge reports.
The preprocessing component of the PPS assembles reflectivity measurements
from each volume scan into a fixed polar grid with resolution of one degree in
azimuth and one kilometer in range out to a maximum range of two hundred and
thirty kilometers. The reflectivity values to be used for each polar grid point are based
on the four lowest elevation angles. The chosen tilt angle depends on which angle is
closest to an optimum altitude, unless it is blocked more than 50% by terrain. The
72
bottom of the chosen beam must also clear the ground by at least 150 meters. The
lowest elevation angle that satisfies these criteria is used.
Biscan maximization (BM) is a procedure in the preprocessing component of
the PPS that assigns the maximum reflectivity from one of the two lowest elevation
angles to the polar grid bin. This was implemented to account for underestimation due
to beam blockage at the lowest elevation angle. Several studies have shown that for
the range of approximately 50-150 kilometers, BM may overestimate due to
intersection with the brightband. It is now performed only for ranges beyond 180
kilometers; where the second elevation angle is above the bright band.
Another part of the preprocessing component of the PPS includes correction
for occultation. No correction is made for bins where the blockage of the beam
exceeds 60% because the next higher tilt angle level is used instead to construct the
hybrid scan. The preprocessing component also includes quality control for isolated
reflectivity data that are large in magnitude but small in area.
The second main component of the PPS is the rain-rate algorithm. It consists
of four parts: (1) a conversion from reflectivity to rain rate, (2) correction for hail, (3)
a time continuity test, and (4) range degradation correction. This algorithm executes
every volume scan after completion of the preprocessing algorithm. It utilizes a
standard Z-R power law relationship. It was noted that the Z-R parameter settings are
designed to be adjustable at each site. The correction for hail sub-algorithm imposes a
maximum value on the reflectivity, which reflects an assumed maximum expected
73
instantaneous, rainfall-rate. This is performed to prevent overestimation from hail.
The time continuity test is a quality control step that checks whether to accept or
reject the current volume scan based on a comparison with the previous volume scan,
predetermined threshold levels, and expected normal precipitation development or
decay. A correction of the calculated rainfall-rate for range is currently not in
operation, but it is expected that a simple equation will be used to correct for range as
long term data become available.
The accumulation component of the PPS algorithm consists of three parts: (1)
integrating rate scans, (2) accounting for missing periods, and (3) removal of hourly
outliers. Storm total rainfall accumulation is computed by integrating consecutive rate
scans.
The fourth component of the PPS algorithm utilizes real-time rain gauge data
to adjust the radar-rainfall estimates; however, this component is currently not being
used because the communication system between the gauges and the WSR-88D has
not been completed. The adjustment will occur on an hourly basis. A mean field bias
will be computed for the entire radar scanning domain, based on a gauge-radar
multiplicative bias. This fourth component consists of four parts: (1) assembly of
hourly gauge accumulations, (2) assembly of hourly radar accumulations, (3) quality
control of gauge-radar pairs, and (4) a Kalman filter bias estimation.
The final processing step of the PPS is the products algorithm. The formatted
products generated at the RPG are sent to the principal user processor (PUP), an
74
external workstation, for display and manipulation. The products generated include
graphics products, digital products, and alphanumeric products.
The paper by Fulton et al. (1998) included a broad review of the current
deficiencies of the PPS algorithm. These include: (1) parameter optimization, (2)
bright band and snow, (3) range degradation, (4) reflectivity calibration and clutter
suppression, (5) anomalous propagation, (6) local gauge-radar biases, (7) attenuation,
and (8) dual polarization.
2.2.2. Radar-Rainfall Studies
Austin (1987) investigated several physical factors (the measurement of
equivalent radar reflectivity, uncertainty in Z-R relations, effects of variations in
precipitation with height, hail and vertical air motions in convective cells, anomalous
propagation, and discrepancies which result from sampling) which influence the
relationship between radar reflectivity measurements and surface rainfall rates, not
only from a theoretical perspective, but also based on a detailed comparison of
raingauge and radar records. Results were tabulated based on the physical factor, the
magnitude of its effect, and situations where it is significant. Methods were suggested
to compensate for the various effects. For intense convective storms, it was noted that
raindrops are larger than normal, and that measured reflectivity may be increased 1 -2
dBZ. A Z-R relation of the form Z = 400/?'3 was suggested to compensate for the
larger than average raindrops. The enhancement of radar reflectivity, by as much as
75
10 dBZ or more, due to the presence of hailstones was concluded to be the factor that
has the greatest effect on radar measurements of surface rainfall in intense convective
storms. Imposing a maximum expected rainfall rate was suggested to compensate for
this effect.
Rasmussen et al. (1989) implemented the NEXRAD precipitation processing
system at the Program for Regional Observing and Forecasting Services (PROFS) in
1988 and ran the algorithm in real time using data from the NCAR CP-2 10 cm
Doppler radar. For raingauge amounts greater than 0.1 inches, the raingauge and
radar estimates agreed to within a factor of two approximately 71% of the time. The
authors explained that a single Z-R relation is inadequate for the entire day and entire
surveillance area for the Colorado Front Range region. They further noted that the
area wide bias that is recomputed once per hour to be applied to all radar rain
estimates for the subsequent hour lags the evolution of the precipitation regime. The
authors believed the raingauge density of the PROFS mesonet was inadequate to
compute a good bias. The authors utilized the Z-R relation Z = 500/?'3. This relation
was used instead of the NEXRAD recommended Z-R relation because it produced
more reliable estimates for light rain. To account for the high frequency of occurrence
of hail for the Colorado Front Range region, the authors imposed a hail cap on the
radar reflectivity to compensate for hail induced reflectivity enhancement. Rasmussen
et al. (1989) recommended the consideration of differential reflectivity measurements
76
to compensate for the effects of hail, bright band, and the space-time variability of the
relation between reflectivity and rainfall-rate.
Smith and Lipschutz (1990) evaluated the accuracy of precipitation estimates
produced by the NEXRAD precipitation processing system by comparing the
precipitation estimates to rainfall measured by the raingauges in the Program for
Regional Observing and Forecasting Services (PROFS) mesoscale surface
observation network (Mesonet). Their evaluation included both stratiform and
convective events, as well as events that included snowfall. The hourly computed
radar-raingauge bias correction factor was restricted to lie between 0.5 and 2. To
account for the enhancement of reflectivity due to the presence of hail, a maximum
radar reflectivity of 50 dBZ was imposed. The authors noted that the Z-R relation
Z = 500/?13 produced the best agreement between the radar and the raingauges in
Colorado during 1988, so this Z-R relation was used throughout the year in 1989. The
radar and raingauge totals agreed within a factor of two for less than two-thirds of the
total number of storms studied. The bias adjustment was shown to have little effect on
the accuracy of the radar estimates of rainfall. However, the authors noted that a
denser raingauge network might have made the bias adjustment more significant.
Fulton (1999) examined the performance of the WSR-88D PPS for the
Buffalo Creek flash flood event of July 12, 1996. In particular, the study focused on
the sensitivity of the algorithms rain-rate threshold parameter and the performance of
the gauge-radar adjustment subalgorithm. On July 12, 1996, a thunderstorm event led
77
to intense flooding, erosion and property damage, and two deaths in the town of
Buffalo Creek, located approximately 35 miles southwest of Denver, Colorado. For
the study, 145 rain gauges were used to compare with radar-derived rainfall estimates
from the Denver WSR-88D radar. However, only one rain gauge sampled the flood
producing storm. The PPS significantly over-estimated rainfall relative to rain gauges
by about 60%. However, the radar estimate over the town of Buffalo Creek was
within 6% of a citizens gauge observation. A rain-rate threshold of 75 mm/h
(equivalent to 51 dBZ) was in use by the PPS on the day of the event. Fulton (1999)
concluded that radar-estimated rainfall estimates would still have been overestimated
even if the rain-rate threshold were reduced to 54 mm/h (49 dBZ). The default WSR-
88D Z-R relation was used in the study, and it was noted that this relation may not
have been optimal for the environment.
Fulton (1999) found the one-hour basin-average radar-rainfall never reaches
more than one-third of the one-hour Flash Flood Guidance (1992). Johnson (2000)
noted that basin averaging reduces input rainfall amounts unless the basin scale is the
same as the storm; this was the case in Buffalo Creek. Flash flood guidance is the
average rain needed over an area for a specified duration to initiate flooding on small
streams in an area or at a specific location immediately downstream of the area. It is a
product used by forecasters at WFOs, together with estimates and forecasts of
rainfall, to support the issuing of flash flood watches and warnings. Sweeney (1992)
proposed steps to modernize the areal FFG, stressing the need for more uniformity in
78
the determination of threshold runoff across the River Forecast Centers (RFCs).
Threshold runoff is the amount of runoff needed over an area to initiate flooding.
Current moisture conditions and threshold runoff are required to compute the FFG. If
the computed gauge-radar bias estimates had been applied to the radar estimates, the
basin average rainfall would have been even smaller for this case. Fulton (1999)
suggested that comparisons of FFG and basin-maximum radar rainfall would have
been more appropriate to evaluate the flash flood threat due to the small space and
time scales of the storm.
2.2.3. Hydrologic Applications
Various studies have focused on the influence of the spatial and temporal
sampling resolution of precipitation on the computed outflow (Milly and Eagleson
1988; Kouwen and Garland 1989; Ogden 1992; Ogden and Julien 1993, 1994;
Michaud and Sorooshian 1994b; Faures et al. 1995; Winchell et al. 1998). The basic
conclusions have been that runoff generation is highly sensitive to the spatial and
temporal scale of the rainfall input.
Winchell et al. (1998) were motivated by the observation that
1. there has only been one paper addressing how watershed computed outflows are
affected by changing the parameters of the reflectivity to rainfall transformation,
and
79
2. that there has been a bias towards the use of infiltration-excess models to
investigate the effects of rainfall sampling resolution on rainfall-runoff
simulations.
They addressed three unresolved questions related to the application of radar-
estimated rainfall data for rainfall-runoff modeling:
1. How significantly will errors in the precipitation data, due to the transformation of
reflectivity to rainfall, affect runoff simulations?
2. How significantly will the aggregation of the radar product in time and space
affect runoff simulations?, and
3. Do differently modeled surface-runoff mechanisms respond differently to these
sources of precipitation uncertainty?
They considered three scenarios for determining the values of a and b in the power
law model used to transform reflectivity into rainfall rate:
1. They used the default parameter values suggested by the NEXRAD algorithm
developers.
2. The parameters were optimized for the 7 storms that were considered for an area
surrounding the study site: the Timber Creek watershed; a 102 km watershed
located in north central Texas.
3. The parameter values were calibrated for each storm event separately,
acknowledging that these parameters vary from storm to storm.
80
In addition to these three methods for determining the parameters a and b, a
bias correction factor based on the ratio of the precipitation measured by the rain
gages to the precipitation estimated by the radar was either applied or not applied;
resulting in 6 different methods for developing radar-estimated precipitation for each
storm. The authors state that these six methods represent the range in precipitation
data obtainable by different processing techniques. The investigators also studied the
effect of spatial and temporal aggregation of the original 1 km 6 minute data on storm
surface runoff generation; also considering spatial scales of 2 km, 4 km, 8 km, and 16
km. Time resolutions, other than the original 6 min. included 24 min., 42 min., and 60
min. Therefore, for each storm, 20 different precipitation data sets were available.
The results they obtained for the sensitivity of runoff volume generation to
reflectivity-rainfall transformation uncertainty clearly illustrated that significant
errors can result. It was shown that the methods which incorporated storm dependent
Z-R parameter calibration produced better runoff simulations than the other methods.
In addition, the methods which included a bias correction factor resulted in better
runoff simulations. The investigators calculated the infiltration-excess response and
saturation-excess response, separately for the six precipitation data sets for the seven
storms; to determine if these two types of runoff responded differently to varying
precipitation input. They were also interested in studying the relationship between
error in runoff and error in rainfall. It was assumed that the true rainfall was
produced from the parameters a and b using the storm dependent approach including
81
bias correction. True runoff was defined to be the runoff produced from the true
rainfall. Their results indicated that for saturation-excess runoff, the relationship
between rainfall errors and runoff errors is fairly linear. On the other hand, for the
infiltration excess case, the relationship was shown to be highly nonlinear; with small
rainfall errors (on the order of 10% resulting in runoff errors on the order of 170%).
The study also included an investigation into the effect of precipitation sampling
resolution in space and time and how this affects saturation-excess and infiltration-
excess runoff production. Their results indicated there is virtually no sensitivity to the
temporal resolution for saturation-excess runoff; whereas, as the spatial scales were
changed, there were differences, but there was no consistent trend, and the changes
were not very large; on the order of 10-15%. For the case of infiltration-excess, the
trend was for decreasing amounts of runoff as the aggregation in time and space is
increased. The results for the sensitivity of the infiltration-excess case was consistent
with the results obtained from previous authors.
82
3. Model Description
3.1. Introduction
This chapter describes a complex distributed rainfall-runoff model that
expressly includes Monte Carlo simulation as an approach for model calibration and
uncertainty estimation. The model does also allow for manual calibration. It is event-
based, assumes the infiltration-excess runoff mechanism, and the kinematic wave
approximation of the equations describing one-dimensional unsteady free surface
flow as overland or open channel flow is used to route excess rainfall over the land
surface and through an organized stream network. The model is compatible with
raster GIS and spatially and temporally varied rainfall data. The objectives for
developing the model include the need for an advanced distributed rainfall-runoff
model structure that: (1) accommodates some of the issues that have been raised
regarding the calibration of such a model type; and (2) can be used to
investigate/evaluate uncertainty in rainfall-runoff modeling (Uccellini 1999).
Previous advanced distributed rainfall-runoff model applications of the GLUE
methodology of Beven and Binley (1992) have limited the number of model
parameters and negated their spatial variability for practical reasons. The model
described in this chapter has the potential to retain the spatial variability of model
parameters from cell to cell. An approach is employed to incorporate slope
83