Copyright
by
Kwabena Oduro Asante
2000
Approaches to Continental Scale River Flow Routing
by
Kwabena Oduro Asante, B.S., M.S.
Dissertation
Presented to the Faculty of the Graduate School of
The University of Texas at Austin
in Partial Fulfillment
of the Requirements
for the Degree of
Doctor of Philosophy
The University of Texas at Austin
May, 2000
Approaches to Continental Scale River Flow Routing
Approved by
Dissertation Committee:
David R. Maidment, Supervisor
James S. Famiglietti, Cosupervisor
Francisco Olivera
Randall J. Charbeneau
Daene McKinney
v
Acknowledgements
To my parents, Kofi and Akua, and siblings, Akosua, Kwaku and Afua
for their love and support through the many years of school,
To Dr David R. Maidment, my supervisor, for the interest he took
in developing me as a person, a hydrologist and a teacher
and for his infectious enthusiasm,
To Dr James S. Famiglietti, my cosupervisor, for his global perspective
and his flexibility in allowing me to pursue my goals,
To Dr Francisco Olivera, my teacher and friend, for his patient instruction
and his insistence that equations remain a part of every modeling decision,
To Dr Randall Charbeneau and Dr Daene McKinney
for their service as members of my dissertation committee,
To Dr Seann Reed, Dr Zichuan Ye, Ferdinand Hellweger and the other GIS
Hydro researchers whose prior work laid the foundation for mine,
To Cai Ximing, Mary Lear, Ty Lehman, Kathy Rose, the staff of CRWR
and the many fine people I met in Texas,
I am truly grateful and give thanks.
vi
I would also like to acknowledge the following institutions for the use of their
data and hydrologic models
? US Geological Survey EROS Data Center, Sioux Falls, SD
? National Center for Atmospheric Research, Boulder, CO
? Hydrologic Engineering Center of the US Army Corp of Engineers
? Center for Research in Water Resources, Univ. of Texas at Austin
vii
Approaches to Continental Scale River Flow Routing
Publication No._____________
Kwabena Oduro Asante, PhD
The University of Texas at Austin, 2000
Supervisor: David R. Maidment
Cosupervisor: James S. Famiglietti
Recent concerns about global climate change and a series of large-scale
hydrologic events such as the Mid-West flood of 1993 and the El Nino of 1997
have focused attention on the need to track the flow of water through the entire
hydrologic cycle. On the land surface, databases of routing parameters and
routing models are required to describe the movement of runoff generated by
Global Circulation Models (GCMs) and other soil water balance models over the
earth?s surface. In this study, a terrain analysis is performed using 30 arcsecond
Digital Elevation Model (DEM) data to develop a global database of terrain
derived routing parameters. A computationally efficient grid based routing model,
called a source to sink (STS) model, is implemented in this study. It routes flow
directly from the point of generation to the desired observation point. The STS
model also allows for easy interaction with models of other phases of the
viii
hydrologic cycle by incorporating the boundaries of their modeling units into the
definition of its own modeling units. A continental scale STS model is created and
parameterized for each continent using datasets derived from the terrain analysis.
A process is defined for determining additional velocity and dispersion
parameters from observed flow data. Another continental scale routing model is
developed using the watershed based approach of the Hydrologic Modeling
System (HMS). Hydrologic elements and routing parameters for the HMS model
are derived from the same terrain data used in parameterizing the STS model.
Basin responses from the two models are compared for various spatial and
temporal resolutions and parameter distributions to determine the implications of
their respective conceptual models. These comparisons show that basin responses
in the STS model are relatively independent of spatial and temporal scale while
the HMS model is scale dependent with regard to both spatial and temporal
resolution. Basin responses for a fine resolution HMS models were successfully
duplicated in a STS model for both the uniform and non-uniformly distributed
velocity and dispersion parameter case.
ix
Table of Contents
LIST OF TABLES XII
LIST OF FIGURES XIII
CHAPTER 11
Introduction .............................................................................................................1
1.1 Problem Statement....................................................................................1
1.2 Objectives .................................................................................................4
1.3 Research Outline.......................................................................................5
CHAPTER 27
Literature Review ....................................................................................................7
2.1 Conceptual Models of a River Basin........................................................7
2.1.1 The Watershed Based Approach ..................................................8
2.1.2 The Cell to Cell Routing Approach..............................................9
2.1.3 The Source to Sink Routing Approach.......................................10
2.2 Methods of Characterizing Flow ............................................................10
2.2.1 Translation With Incidental Dispersion......................................11
2.2.2 Translation With An Approximate Dispersion Process .............12
2.2.3 Translation With Fully Described Dispersion ............................13
2.3 Continental Scale Routing Models .........................................................17
2.4 Scaling Issues in Routing .......................................................................22
2.4.1 Scale Dependence of Routing Methods......................................23
2.4.2 Scale Dependence of Routing Parameters..................................25
CHAPTER 327
Data Development .................................................................................................27
3.1 Data Sources ...........................................................................................27
3.2 Terrain Database Development ..............................................................29
3.3 Preprocessing For the Source to Sink Model .........................................40
x
3.4 Preprocessing For the Hydrologic Modeling System.............................55
CHAPTER 464
Methodology..........................................................................................................64
4.1 The Source to Sink Modeling Approach ................................................64
4.1.1 STS Model Assumptions............................................................64
4.1.2 STS Model Components.............................................................67
4.1.3 STS Routing Methodology.........................................................70
4.1.4 Estimating STS Routing Parameters ..........................................78
4.2 The Hydrologic Modeling System Approach.........................................83
4.2.1 HMS Model Assumptions ..........................................................84
4.2.2 HMS Model Components...........................................................85
4.2.3 HMS Routing Methodology .......................................................88
4.2.4 Estimating HMS Routing Parameters.........................................93
CHAPTER 596
Model Applications ...............................................................................................96
5.1 The application basins ............................................................................96
5.1.1 The Congo Basin ........................................................................97
5.1.2 The Nile Basin ..........................................................................100
5.2 Source to sink model runs ....................................................................102
5.2.1 Longitudinal Decomposability of the STS Model....................103
5.2.2 Effect of Spatial Resolution of the STS Model ........................107
5.2.3 Effect of Spatial Distribution of STS Parameters.....................113
5.2.4 Effect of Temporal Resolution of STS Model..........................122
5.3 Hydrologic Modeling System Model Runs ..........................................124
5.3.1 Longitudinal Decomposability of the HMS Model..................127
5.3.2 Effect of Spatial Resolution of the HMS Model ......................132
5.3.3 Effect of Spatial Distribution of HMS Parameters...................136
5.3.4 Effect of Temporal Resolution of HMS Model........................141
xi
5.4 Comparison with Observed Flows .......................................................151
5.4.1 Comparing STS Simulated Flows with Observed Flows .........152
5.4.2 Comparing HMS Simulated Flows with Observed Flows .......158
CHAPTER 6 164
Conclusions and Recommendations....................................................................164
6.1 Conclusions ..........................................................................................164
6.1.1 Inferences from Global Terrain Analysis .................................166
6.1.2 Inferences from Global STS Model Implementation ...............167
6.1.3 Inferences from HMS Model Development .............................169
6.1.4 Inferences from Comparing STS and HMS Models ................170
6.2 Recommendations ................................................................................172
APPENDICES 174
Appendix 1: AML Code for Source to Sink Model Preprocessing.....................175
Appendix 2: Parameters of the T42 mesh ...........................................................179
Appendix 3: FORTRAN Code for Source to Sink Routing Model.....................183
Appendix 4: AML Code for Computing Equivalent Muskingum Parameters
from a Field of Velocities and Dispersion Coefficients .............................193
Appendix 5: 10 day flows in the Nile basin, 1953-57 .........................................195
REFERENCES 197
VITA 204
xii
List of Tables
Table 3-1: Projection Centers for Various Continents ..........................................31
Table 3-2: Sample Projection File for Africa ........................................................31
Table 3-3: Inland Catchments Identified ...............................................................34
Table 5-1: Effect of Spatial Resolution on Lagged and Dispersed Basin Response
.....................................................................................................................112
Table 5-2: Flow Routing Parameters used in the Nile Basin ..............................118
Table 5-3: Effect of STS Routing Interval on Congo Basin Response ...............123
Table 5-4: Parameters for the Congo Basin Delineation.....................................126
Table 5-5: Effect of Muskingum n on River Reach Response ............................130
Table 5-6: Range of Dispersion Coefficients of a Muskingum Routing Reach with
24 hour Routing Interval and Velocity of 0.3m/s........................................131
Table 5-7: Effect of HMS Enforced Minimum Subbasin Lag ............................143
Table 5-8: Effect of Min. Subbasin Lag Rule on Basin response, Lag Routing .146
Table 5-9: Effect of Subbasin Lag Rule on Basin response, Muskingum Routing
.....................................................................................................................148
Table 5-10: Routing Parameters Computed by Method of Moments .................155
xiii
List of Figures
Figure 3-1: Inland Catchments of the World.........................................................33
Figure 3-2: Flow Direction Arrows .......................................................................35
Figure 3-3: Flow Direction Grid of Africa ............................................................36
Figure 3-4 (a): Histogram of Flow Directions for North America ........................37
Figure 3-4 (b): Histogram of Flow Directions for Africa......................................37
Figure 3-5: Flow Accumulation Grid of Africa.....................................................39
Figure 3-6: Flow Length Grid of Africa................................................................40
Figure 3-7: Components of a Source-To-Sink Model ...........................................42
Figure 3-8: Sink Locations for the Continent of Africa ........................................44
Figure 3-9 (a): Drainage Basins of North America ...............................................46
Figure 3-9 (b): Drainage Basins of Africa.............................................................46
Figure 3-9 (c): Drainage Basins of South America ...............................................47
Figure 3-9 (d): Drainage Basins of Australia ........................................................47
Figure 3-9 (e): Drainage Basins of Europe............................................................48
Figure 3-9 (f): Drainage Basins of Asia ................................................................48
Figure 3-10: Linking Modeling Units for the African Continent..........................50
Figure 3-11 (a): Parameterized Sources of North America...................................52
Figure 3-11 (b): Parameterized Sources of Africa.................................................52
Figure 3-11 (c): Parameterized Sources of South America...................................53
Figure 3-11 (d): Parameterized Sources of Australia ............................................53
Figure 3-11 (e): Parameterized Sources of Europe ...............................................54
Figure 3-11 (f): Parameterized Sources in Asia ....................................................54
Figure 3-12: Delineated Rivers Coverage of Africa..............................................58
Figure 3-13: Delineated Watersheds of Africa......................................................58
Figure 3-14: Line Diagram of Hydrologic Element Connectivity ........................60
Figure 3-15: HMS Model Representation of Africa..............................................61
Figure 4-1: Diagram of STS Model Components .................................................68
Figure 4-2: Diagram of HMS Model Components................................................85
Figure 4-3: Elements of an HMS Basin Model .....................................................86
Figure 4-4: Soil Conservation Service Triangular Hydrograph ............................90
Figure 5-1: The Four Major Basins of Africa........................................................97
Figure 5-2: Area Map of the Congo Basin ............................................................98
Figure 5-3: Major Tributaries and Flow Gages of the Nile Basin.......................100
Figure 5-4: Source to Sink Model of the Congo Basin .......................................103
Figure 5-5: Setup for Evaluating Longitudinal Decomposability .......................104
Figure 5-6: Source to Sink versus Cell to Cell Routing ......................................105
Figure 5-7: Comparison of Diffusion Type STS and CTC Discharge ................106
Figure 5-8: Congo Basin STS Models at 30, 10 and 5 arcminute Resolutions ...108
xiv
Figure 5-9(a): Comparison of Congo Basin Responses for 30 and 10 arcminute
Resolution STS Model with Pure Lag Routing...........................................110
Figure 5-9(b): Comparison of Congo Basin Responses for 10 and 5 arcminute
Resolution STS Model with Pure Lag Routing...........................................110
Figure 5-10: Effect of Spatial Resolution on Dispersed STS Hydrograph..........111
Figure 5-11: Nile Basin (a) Parametric Zones and (b) Sources...........................115
Figure 5-12: Effect of Spatial Distribution of Velocity with Uniform Dispersion
Coefficient ...................................................................................................117
Figure 5-13: Effect of Spatial Distribution of Dispersion Coefficient with Uniform
Velocity .......................................................................................................119
Figure 5-14: Effect of Spatial Distribution of Dispersion Coefficient with Non-
Uniform Velocity.........................................................................................121
Figure 5-15: Effect of STS Routing Interval on Congo Basin Response............123
Figure 5-16: Sample HMS Schematic for Congo Basin......................................125
Figure 5-17: HMS Schematic for Longitudinal Decomposability Runs .............127
Figure 5-18: Changes in River Reach Response with Muskingum n..................129
Figure 5-19: Effect of Delineation Threshold on Congo Basin Response for HMS
Model with Muskingum Routing ................................................................133
Figure 5-20: Comparison of 1000 Cell Threshold HMS Model Response with the
STS Model Responses in the Congo Basin .................................................134
Figure 5-21: Comparison of 1000 Cell Threshold HMS Model Response with the
STS Model Responses in the Congo Basin .................................................135
Figure 5-22: Basin Response for Iterations of Muskingum x and n Parameters.139
Figure 5-23: Equivalent HMS and STS Basin Response for Distributed
Parameters ...................................................................................................140
Figure 5-24(a): Effect of Minimum Subbasin Lag for 24hr-Interval Lag Routing
.....................................................................................................................145
Figure 5-24(b): Effect of Min. Subbasin Lag for 12hr-Interval Lag Routing .....145
Figure 5-25(a): Effect of Min. Subbasin Lag for 24hr-Interval Muskingum
Routing ........................................................................................................147
Figure 5-25(b): Effect of Min. Subbasin Lag for 12hr-Interval Muskingum
Routing ........................................................................................................147
Figure 5-26(a): Effect of Routing Interval on Response with Lag Routing........149
Figure 5-26(b): Effect of Routing Interval on Response with Muskingum Routing
.....................................................................................................................149
Figure 5-27: Inflow and Discharge Hydrographs in the Sudd Marshes ..............153
Figure 5-28: Excess Input and Discharge Runoff after Baseflow Separation for (a)
Period 1, (b) Period 2 and (c) Period 3 ........................................................154
Figure 5-29: Routing Excess Runoff using Computed Velocities and Dispersion
Coefficients for (a) Period 1, (b) Period 2 and (c) Period 3 ........................156
Figure 5-30: Effect of Routing Input Flow using Mean Velocity and Dispersion
Coefficients..................................................................................................157
xv
Figure 5-31: Nile Basin (a) Gauge Locations and (b) HMS Model Components
.....................................................................................................................160
Figure 5-32: Observed and HMS Simulated Flows at Malakal...........................161
Figure 5-33: Observed and HMS Simulated Flows at Khartoum .......................162
Figure 5-34: Observed and HMS Simulated Flows at Aswan.............................163
1
CHAPTER 1
Introduction
1.1 PROBLEM STATEMENT
The desire to understand the movement of water over the earth?s surface
has long been a fascination. This desire has mainly arisen from the need to
evaluate the amount of water available at a particular location to meet local
demand as well as the risk of flooding due to excess water. Hydrologists have
been called upon to determine the quantity of flow along rivers at various points
in time and space. Until recently, the scope of these problems has been limited to
individual communities and consequently, most hydrologic studies have focused
on small watersheds. Routing methods have been developed based on these small
watersheds. Recent concerns about global climate change and a series of large-
scale hydrologic events such as the Midwest flood of 1993 and the El Nino of
1997 have focused attention on the need to track the flow of water through the
entire global hydrologic cycle. In the climate modeling community, in particular,
there is a need for models which route runoff generated by the land surface
component of Global Circulation Models to the ocean cells used in the ocean
modeling component. Approaches for routing water through large-scale systems
are slowly evolving in response to these modeling needs.
One evolving approach for large scale routing involves dividing the
earth?s surface into a series of cells using a uniform grid. A vertical soil water
2
balance is performed on each of these cells, and the resulting excess runoff is
routed from one cell to the next until it gets to an outlet or sink. Watershed
boundaries are not explicitly delineated. Instead, the flow direction assigned to
these coarse resolution cells dictates where the watershed divide should lie. The
flow length is simply estimated on the basis of the distance between the centers of
adjacent cells. This approach is not computationally efficient since routing is
performed throughout the watershed even at locations where discharge values are
not required. A trade off must be made between routing detail and the number of
modeling units. When used at the continental scale, the routing process in grid
based models is often reduced to very simple representations in order to increase
computational efficiency. A routing model which computes discharge only at
required locations allows a more detailed description of the routing process
without sacrificing the spatial resolution of the model.
Grid based approaches differ significantly from other hydrologic models
which use a watershed as the primary modeling unit. These watershed-based
models are built on the fact that water flows on or near the land surface for a
while before entering the river system. The land surface is therefore divided into
watersheds with river reaches serving to connect the watersheds to each other and
to the discharge location at the basin outlet. Therefore, in modeling river basins,
different routing methods are defined for the respective watershed and in-stream
phases of flow. Models of land-atmosphere interactions are used to generate input
runoff for the routing model. Since the atmosphere is a spatially continuous
system, it is more convenient to employ the uniform cells used in atmospheric
3
circulation models. The modeling units in the land-atmosphere interaction models
and the watershed-based routing models are hence not coincident. Although there
are ways of accounting for this discrepancy, it can be a real impediment when
dealing with runoff which is variable in both space and time. There is hence an
advantage to performing routing in a gridded environment when working with
time-varying input runoff. The ideal situation would be to work in a routing
environment which incorporates the boundaries of both the drainage basin and the
soil water balance modeling units.
The relationship between the grid based approach and the watershed
modeling approach is not clearly understood. The results of the two modeling
approaches cannot be compared because the boundaries of the modeling units do
not coincide. There is also a discrepancy in the spatial scale at which the two
modeling approaches have been applied. Grid based models were developed for
continental scale modeling while watershed based approaches have traditionally
been used for local and regional scale modeling. In order to compare the two
modeling approaches, the basin boundaries defined in one model should be
duplicated in the other. Applying the respective modeling approaches to basins of
comparable spatial scale also helps to bridge the discrepancy in spatial scale
between climate and watershed models. Current concerns about global and
continental scale hydrologic events provide the impetus for performing such
intercomparisons in large river basins.
The development of continental scale routing models requires detailed
information about the distribution of geographic features and parameters which
4
influence flow. Geographical Information Systems (GIS) have recently been used
in deriving vector representations of watersheds and rivers for hydrologic
modeling. GIS based approaches have also been used for deriving hydrologic
modeling parameters from gridded representations of the earth. GIS thus provides
an ideal environment for developing and parameterizing continental scale routing
models.
In summary, there is a need for cell-based routing models to support
continental scale hydrologic modeling. Existing models are limited by their
inability to accurately define basin boundaries, by inconsistencies between their
modeling units and those of associated models and by their computational
intensity. A new routing approach which seeks to resolve these inconsistencies in
modeling unit definition while only performing computations at required
locations is implemented in this study.
1.2 OBJECTIVES
A global database of hydrologic parameters is required to support surface
water routing around the earth. Such a database allows the modeler to derive
routing parameters for any modeling units defined in a study area. The first
objective of this study is to develop a GIS database to support large-scale surface
water routing globally.
The second objective is to implement a modeling framework which
incorporates basin boundaries in a grid based model while maintaining
computational efficiency by only performing routing at desired locations. The
5
efficacy of the source to sink modeling framework is demonstrated by using it to
route runoff generated by the land surface component of a Global Climate Model
at the continental scale.
A variety of topographic and climatic conditions are encountered when
routing is performed at the continental scale. The third objective of this study is to
examine of the robustness of the source to sink approach as compared to the
watershed based approach in continental scale applications.
1.3 RESEARCH OUTLINE
In meeting the first objective, existing datasets which could be used in
deriving hydrologic parameters globally were identified. These datasets were
processed to develop a database of routing parameters.
In meeting the second objective, a continental scale routing model using
the source to sink approach was developed for all the continents except Antarctica
which has no rivers. The model was used to route 10 years of daily runoff data
from a Global Climate Model to receiving ocean cells.
In meeting the third objective, routing models were developed for selected
drainage basins using the watershed based and grid based approaches. The models
were parameterized from the global database of hydrologic parameters thus
producing two independent and yet consistent hydrologic representations of the
same drainage basins. Various simulations were performed to determine
differences and similarities in the response of the modeling approaches to changes
in modeling parameters. To verify its suitability for continental scale routing, the
6
watershed based model was used to route flow in a large drainage basin, and the
results were compared with actual flows in the basin.
The remainder of this dissertation is divided as follows. Chapter 2
describes the various conceptual models and methods used in routing. It also
provides a review of prior developments in continental scale routing. In Chapter
3, sources of input data and the development of the global database of hydrologic
parameters are described. The procedure for developing the each of the
continental scale routing models is also outlined. The theoretical basis of the two
continental scale modeling approaches are discussed in Chapter 4 along with
methods of determining the required model parameters. This is following in
Chapter 5 by the presentation of results of model runs, and the discussion of the
results in Chapter 6. In Chapter 7, the findings of the study are summarized with
recommendations for further study.
7
CHAPTER 2
Literature Review
2.1 CONCEPTUAL MODELS OF A RIVER BASIN
In developing models of hydrologic systems, it is recognized that
processes that occur in nature cannot be reproduced exactly. Even if such a
representation of natural processes were possible, the resulting model would be so
complex that it would be impossible to solve mathematically. Some simplifying
assumptions are consequently made to reduce the complex system to a set of
simple descriptive processes (Dooge, 1973, pp.148). Only those processes
considered critical to the functioning of the system are described. In this chapter,
the conceptual models and methods used in routing are discussed. Existing large-
scale routing models are reviewed along with their limitations. Problems of scale
dependence affecting routing models are also discussed. To begin with, the
hydrologic system is reduced to a simple conceptual model. The control volume
proposed by Chow et al (1988, pp.7) provides a conceptual model that is
applicable to hydrologic systems in general. They define a hydrologic system as
?a structure or volume in space, surrounded by a boundary, that accepts
water and other inputs, operates on them internally, and produces them as
outputs?.
8
Hydrologists develop models by trying to relate the input into the control
volume with the output using a system response or transformation. The first task
in modeling must therefore be the definition of the control volume. There are
several approaches to defining a control volume for routing applications. The
most common of these can be grouped into three classes as described in sections
2.1.1 to 2.1.3 below.
2.1.1 The Watershed Based Approach
Hydrologists have for a long time envisaged the land surface component
of the hydrologic cycle as being made up of river reaches with associated
drainage areas. Reaches are defined as stretches of river with no major tributaries
in between, and the portion of land draining directly into each river reach is
defined as its associated watershed. River junctions are defined wherever two or
more river reaches meet, and diversions are defined wherever two or more
reaches form downstream of a single reach. Other natural features such as lakes
and springs and manmade features such as dams, reservoirs and canals may also
be included in the definition of the hydrologic system. In the systems approach,
each of these features is regarded as a control volume, and a separate system
response is defined to transform the inputs in each control volume into discharge
at its downstream end. Sequential links are established among the control volumes
such that flow can be passed from one control volume to the next until it arrives at
a terminal point such as the continental margin or an inland water body with no
outlet. The Hydrologic Modeling System (HMS) developed by the Hydrologic
9
Engineering Center of the US Army Corps of Engineers (Peters, 1998) and the
MIKE BASIN modeling system developed by the Danish Hydraulic Institute
(Storm and Kjelds, 1999) are examples of models built on this concept.
2.1.2 The Cell to Cell Routing Approach
A second approach to modeling flow over the land surface involves
discretizing the land surface into smaller segments, and routing flow from one
segment to the next until it arrives at a terminal point. The underlying assumption
in this approach is that the land surface within each of these segments is
homogenous. It is therefore not necessary to distinguish between overland and in-
stream flow processes. For convenience and computational efficiency, the land
segments used are often equaled sized cells and hence the name cell to cell
routing. The approach can however be used with other irregularly shaped land
segments. The control volume in a cell to cell routing approach is the land
segment. Inputs to each segment are transferred to the next downstream segment
using a unique response function. The response function is expected to capture the
essence of the various flow processes occurring within the control volume. Its
growing usage in surface flow modeling has been spurred by the need to couple
land surface models with models of the atmospheric and subsurface flow phases
of the hydrologic cycle. Watershed boundaries defined in the watershed modeling
approach (see section 2.1.1) are of little relevance in these other phases. Examples
of surface water routing models using this approach include those developed by
V?r?smarty et al (1989), Miller et al (1994), Sausen et al (1994) and Coe (1997).
10
2.1.3 The Source to Sink Routing Approach
The final approach to characterizing surface flow involves subdividing the
land surface into smaller segments and routing the flow from each segment
directly to a discharge location such as the continental margin or an inland
depression. In this approach, the control volume is drawn around the flow path
linking the segment to the discharge location. The system transformation
describes the change in distribution of the input when it arrives at the discharge
location, allowing for losses within the control volume. The flow path may be
composed of different flow regimes such as overland flow and in-stream flow.
Water may flow at different speeds through these various flow regimes. The
response function defined for each land segment is in effect a summary of all the
transformations that occurs along the flow path between the source of input and
the discharge location (also called the sink). The name source to sink routing is
used by Olivera et al (2000) to describe all models of this kind. Examples of
source to sink models include the Clark unit graph method (Clark, 1945), the
modified Clark method (Kull, 1998), the methods of Naden (1992,1993) and
Olivera and Maidment (1999).
2.2 METHODS OF CHARACTERIZING FLOW
After reviewing control volume definitions used in river basin modeling,
the next step is to define the flow routing processes and the nature of the system
transformations used to characterize these processes. As mentioned in section 2.1
of this dissertation, the task in hydrologic modeling is to discern the essential
11
characteristics of flow. With regard to transformation of flow through a control
volume, these may be summarized as translation, redistribution and translational
losses. Translation refers to the simple movement of water from the input location
to the discharge location. During translation, the shape of the input pulse or
impulse is altered, mainly as a result of three mechanisms. The shape may change
due to shearing forces between the flowing water and the channel surface with
which it comes in contact (Maidment, 1993, pp.14.14). The shape may also
change due to the damping effect of off-stream storage created by the constant
variation in channel dimensions, and by passage through wide water bodies such
as lakes and reservoirs (Henderson, 1966, p.355). These first two mechanisms are
usually combined into a single mechanism known as flow redistribution. The third
mechanism by which the shape of the input pulse may change is by loss of flow
through seepage and evaporation. This third mechanism, referred to as
translational loss, is usually treated separately since it implies changes in the mass
balance.
2.2.1 Translation With Incidental Dispersion
The simplest form of routing involves describing the movement of water
through a control volume using the continuity equation along with a mathematical
relationship between discharge and control volume storage. In these methods,
there is no parameter for explicitly controlling dispersion. Redistribution is
achieved as an incidental by-product of the routing process. These methods are
jointly referred to as level pool routing (Maidment, 1993, pp.10.6). The most
12
commonly used level pool routing method is the linear reservoir method (Zoch,
1934) in which discharge from the control volume is assumed to be directly
proportional to storage. Another common approach to level pool routing involves
the use of approximate numerical methods such as the Runge Kutta to determine
discharge from the control volume (Maidment, 1993, pp.10.6). In this method,
water is routed sequentially through several discontinuous segments of a river
reach. The Puls and Modified Puls methods (Peters, 1998) in which a relationship
is established between the storage and discharge from a river reach based on
observed data are also examples of level pool routing.
2.2.2 Translation With An Approximate Dispersion Process
The second class of routing methods are those that contain a parameter
which, though not descriptive of the dispersion process, can be used to control the
magnitude of dispersion experienced by an input passing through the control
volume. The two most common examples of this class of methods are the Cascade
of Linear Reservoirs attributed by Montes (1998) to Kalinin and Milyukov (1958)
and the Muskingum method developed by the US Army Corp of Engineers
(Henderson, 1966, p.362).
The Cascade of Linear Reservoirs approach takes advantage of the scale
dependence of linear reservoir models by routing input flow through a series of
linear reservoirs as it passes through the control volume. While the output of a
single linear reservoir is an exponential distribution, subsequent convolutions of
this distribution yield gamma distributions (Olivera, 1996). The response of a
13
cascade of reservoirs containing more than one reservoir is consequently a gamma
distribution and the number of reservoirs in the cascade can be used to modify the
magnitude of dispersion experienced by an input passing through the control
volume.
The Muskingum method approximates the storage effects in a river reach
with a trapezoidal control volume consisting of a prism with an overlying wedge.
The flow through the prism is assumed to be dependent on discharge only while
flow through the wedge is dependent on the difference between the inflow and
discharge (Henderson, 1966, p.362). A factor, the Muskingum x, is applied to the
wedge storage to account for the rate of attenuation of flow as it passes through
the control volume. The assumption made in applying the Muskingum factor is
that in cases where the rate of attenuation is low, the area of the wedge storage
will approximate that of a triangle, and thus x will assume a value of 0.5. In all
other cases, the area of the wedge is a concave parabola, and x will assume a
value between 0 and 0.5. The Muskingum x can thus be used to determine the
amount of dispersion experienced by flow passing through the control volume.
2.2.3 Translation With Fully Described Dispersion
The third class of routing methods trace the movement of flow through the
control volume with transformations which fully describe both the translation and
dispersion processes. These methods generally rely on the momentum equation
for unsteady, non-uniform flow (Chow et al, 1988, p.31) shown as Equation 2-1.
14
Equation 2-1 0S
x
h
g
x
VV
t
V
f =
?
?
?
?
?
?
+
?
?
+
?
?
+
?
?
where A is the cross sectional area of the channel
V is the flow velocity
S
f
is the friction slope
x
h
?
?
is the slope of piezometric head
g is the gravitational constant
x is the distance along the longitudinal axis of the channel
t is time
The momentum equation describes the change of form of an input
subjected to the forces of friction and gravity. Since these are the same forces
responsible for the deformation of the input as it travels through the control
volume, this equation in effect provides a full description of the dispersion
process. The momentum equation and the continuity equation (given as Equation
2-2 below) for flow routing are together referred to as the St. Venant equations
after the 19
th
century French mathematician who derived them in 1871.
Equation 2-2 0
)(
q
t
A
x
AV
=?
?
?
+
?
?
where
x
AV
?
? )(
is the rate of change of flow rate with distance
t
A
?
?
is the rate of change of area of flow with time
15
q is lateral inflow or outflow along the channel
There is no analytic solution to these equations in their complete form.
Some simplifying assumptions are therefore made to enable the equations to be
solved. These solutions are therefore limited in applicability to conditions under
which the simplifying assumptions can be shown to hold true. The assumptions
and their resulting solutions can be classified into three categories, namely
dynamic, kinematic and diffusion wave solutions.
In dynamic wave routing, the full St. Venant equations are used but
assumptions are made about the condition of flow at the boundaries. These
boundary flow assumptions provide a starting point for computations which then
proceed using one of several numerical integration techniques. The most common
of these techniques are described by Montes (1998, pp.470). They include the
method of characteristics (Massau, 1889, Richardson, 1910 and Lin, 1952),
explicit methods (Stoker, 1953) and implicit methods (Preissmann, 1961, Cunge,
1961, 1966, Lai, 1966). Many of these numerical integration methods can only be
applied over small space and time steps over which the magnitude of error
introduced by the approximation of the integral is small compared to the
magnitude of the flow. These numerical stability limitations have been explored
in the works of Von Neumann (1950) and Courant and Freidriechs (1948).
Dynamic wave routing is generally unsuitable for continental scale applications
because the numerical stability limitations make the methods too computationally
intensive to be practicable at that scale.
16
Kinematic and diffusion wave routing methods on the other hand make
some assumptions about the condition of the flow parameters. These assumptions
allow one or more of the terms in the momentum equation to be dropped thus
making it possible to find an exact analytical solution to the St. Venant equations.
The development of the kinematic wave approach is attributed by Montes (1998,
p.547) to the works of Kleitz (1877), Graeff (1875), Forchheimer (1930), and
Lighthill and Whitham (1954). It assumes that within the control volume, velocity
is constant in time and space, and the water surface is parallel to the channel bed
(Chow et al, 1988, p.281). The momentum equation is consequently reduced to
two terms expressing the equality of bed slope and friction slope. The kinematic
wave approach cannot be used, on its own, for routing flows at the continental
scale because the assumption of a constant water depth does not allow the flood
wave to attenuate as it travels through the control volume. Peak attenuation is
considered to be an important property in continental scale routing. Kinematic
wave routing has however found numerous applications in characterizing
overland flow (Iwagaki, 1995, Henderson and Wooding, 1964, Woolhiser and
Liggett, 1967).
The diffusion wave routing approach differs from the kinematic wave
approach in that it does not make the assumption of a constant water depth within
the control volume. The main parametric assumptions made in this approach are
the spatial and temporal invariance of flow velocity. The method was first
proposed by Hayami (1951) who demonstrated that by making this simplification,
an exact solution of the St. Venant equations could be obtained. Hayami
17
subdivided the input into a series of impulses to which the diffusion wave impulse
response function could be applied. The discharge due to each of these impulses
was then aggregated at the outlet.
The assumption of time invariant parameters (velocity and dispersion
coefficient) within the control volume is one obvious limitation of the diffusion
wave method. In real river systems, large flows are translated faster than smaller
flows. A solution to this problem was provided by the multiple linearization
technique of McQuiver and Keefer (1974) in which the input is further subdivided
vertically with different response functions being applied to different magnitudes
of flow. Subsequent researchers have also addressed the problem of spatial
invariance of parameters by defining ways to estimate a single flow path
parameter for a path consisting of several segments each with its own parameters
(Olivera, 1996). The flow path parameters thus computed result in a discharge
distribution with properties similar to that resulting from sequential routing
through multiple segments.
2.3 CONTINENTAL SCALE ROUTING MODELS
V?r?smarty et al (1989) developed a coupled water balance and water
transport model for continental scale routing in South America. In their model, the
land surface is divided into a 0.5 by 0.5 degree grid, and flow is routed
horizontally from cell to cell. Flow through the system is modeled using the
continuity equation and a linear reservoir storage function. However, the model
attempts to subdivide the storage into separate flood plain and in-stream
18
components. This introduces a non-linearity in the model which cannot be
resolved analytically. A Runge-Kutta integration technique is used to solve for the
two storage components. Flow velocity in this model is computed using the mean
annual discharge and the empirical velocity functions proposed by Leopold
(1953). The main limitation of this model is that it contains several empirical
relationships, the validity of which cannot be verified.
Miller et al (1994) proposed to allow for spatially varying storage
characteristics of the cell to cell model by making the transfer coefficient in the
V?r?smarty model a function of excess storage instead of total storage. In
essence, this model implies that in-cell storage is caused by the inability of the
cell to transfer water to its downstream neighbor. A sill depth is defined for each
cell as the depth of the river channel, and this depth is used to determine a
minimum cell storage which must be met before flow is transmitted from a cell to
a neighboring cell. Flow in excess of this minimum storage is transmitted to the
next downstream cell by a linear storage function as in the V?r?smarty model.
The model routes flow globally at a resolution of 2.0 by 2.5 degrees. The flow
direction for cells in this model is determined by visual inspection. Flow velocity
is determined by an empirical function relating slope within each grid cell to a
reference slope and a reference velocity. Again, there is no physical basis for the
storage and velocity relationships. The definition of a minimum cell storage is
inconsistent with the mechanisms of storage in river systems defined in section
2.2 of this dissertation. Besides, no method is suggested for computing the depth
of river channels in each grid cell globally.
19
Coe (1997) proposed to improve Miller?s (1994) cell to cell approach by
using lakes and other natural depressions in the land surface to determine the cell
storage. In this model, a 2 by 2 degree grid is used to divide the land surface into
routing cells. As in Miller?s model, discharge is assumed to occur only when the
minimum cell storage is exceeded. Consequently, the discharge rate is a function
of the difference between the amount of water in the cell and its minimum storage
requirements. The minimum storage of each 2 by 2 degree cell is determined by
filling the sinks in a finer resolution grid, in this case a 5 by 5 minute grid. Using
this approach, regions with no internal drainage are assigned a minimum storage
of zero. Sinks in DEMs may result from natural depressions or from errors in the
DEM. The raw DEM cannot therefore be used for determining storage without
eliminating the erroneous pits. No method is suggested for distinguishing between
erroneous pits and true natural depressions in this application. It is consequently
unclear whether the cell storage determined from the 5 by 5 minute grids in this
application are representative of actual storage capacities in the land surface.
An alternative method of determining discharge in the cell to cell
environment is proposed by Sausen et al (1994). In this model, cells on the land
surface are assumed to have no storage. All the water available in each cell is
allowed to flow in any one of four orthogonal directions. The quantity of water
flowing in each direction is determined by an empirical function involving a
translation velocity and the land surface slope. As with other empirical functions,
there is no physical basis for its use. The validity of the relationship is further
thrown into question by the generation of negative discharges during the study.
20
Kite et al (1994) used a spatially distributed version of the SLURP
(Simple Lumped Reservoir Parametric) model to route runoff generated by a
GCM. In this approach, watersheds are divided into grouped response units
(GRU) using a 10 km by 10 km sized mesh. Within each GRU, flow is routed to
an outlet using a velocity calculated for each land use category. The velocity is
determined from Manning?s equation by assuming a hydraulic radius of 1. This
assumption is meant to imply that the width of the river is much larger than its
depth. However, a review of the equation for hydraulic radius would indicate that
as the width of a channel increases the hydraulic radius does in fact tend towards
the depth of flow, not 1. Manning?s roughness coefficient, n, is assigned to each
cell based on its land cover category allowing velocity to be computed. The
velocity is then used to compute isochrones of travel time as in the Clark
hydrograph method (Clark, 1945). Flow at the watershed outlet is determined by
aggregating flows from the different isochrones. Flow is routed from GRU to
GRU until it arrives at a sink. For this routing a nonlinear, empirical storage
function is used to compute storage and outflow. As with other empirical
functions, there is no physical basis for its use, and its associated parameters must
be determined by calibration.
A more physically based solution is proposed by Naden (1992,1993) who
applied a network response function to the Severn and Thames river basins. The
response function is derived by convolving the linear solution of the diffusion
wave model with the network width function. A network width function is
analogous to an area distance diagram and is defined as the number of river
21
channels within a given distance of the basin outlet divided by the total number of
channels in the basin. In this study, the basins are divided into grids cell (40 by 40
kilometers wide) which only serve to define areas of uniform runoff generation.
Flow routing is performed through the stream network. The main difficulty with
this approach is that it is not easy to discern the relationship between the network
width function and the distance-time diagram created by the diffusion equation.
The latter is a response function describing the longitudinal distribution of input
runoff when it is subjected to both an advective and a dispersive process. The
network width function on the other hand is a measure of the distribution of flow
channels within a basin. No mathematical or physical justification is given for the
convolution. In addition, it is difficult to justify the use of a network width
function when distance area diagrams and spatially distributed flow length grids
can be computed using Geographic Information Systems.
Lohmann et al (1996, 1998) also used a diffusion wave model coupled
with a land surface parameterization scheme to route GCM runoff. In this model,
the land surface is seen conceptually as consisting of watersheds and streams,
each of which has a response function. Hence even when the land surface is
subdivided with a grid, each cell must still be assigned two response functions to
account for the respective watershed and stream segments of flow within the grid
cell. The main limitation of this model lies in the conceptual model. The grid cell
response function accounts for the transmission of flow from within the grid cell
to the outlet of the cell. To use the same response to account for flow entering the
cell from a neighboring cell would not be accurate since such flow enters the river
22
system without entering the watershed. The response function for routing flows
generated within the grid cell is determined by deconvolution from observed flow
data. Even if flow data were available to define grid cell response functions
globally, it would still be difficult to justify the assumption of a grid cell as
consisting of rectangular watersheds which fit exactly into the cell.
Hagemann and Dumenil (1998) apply a cascade of linear reservoirs to
flow routing a cell to cell model. In their model, a grid cell is viewed as being
equivalent to a watershed of equal resolution. Flow processes observed in a 2500
km
2
watershed are consequently applied to 0.5 degree resolution grid cell. Within
each cell, three separate flow paths are defined. Runoff generated within the grid
cell may be routed as overland flow using a cascade of linear reservoir routing or
as baseflow using a linear reservoir. In addition, flow entering the cell from
upstream cell is routed using another cascade of linear reservoirs. Flows from
these three routing processes are aggregated at the outlet of the cell and passed on
to the next cell. The main limitations of this model are that the number of
reservoirs in a cascade is not a physically based parameter and can only be
obtained by calibration. Given that there are three different responses, three
different velocities must also be defined for each grid cell. The conceptual model
of a grid cell being equivalent to a watershed is also difficult to justify.
2.4 SCALING ISSUES IN ROUTING
The field of hydrologic engineering evolved out of the need to address
local flooding problems. Most hydrologic methods were thus developed with
23
local scale analysis in mind. Today, the scale of hydrologic analysis is shifting to
larger scale problems and concerns. The applicability of methods of analysis
developed at the local scale to large-scale problems must consequently be
addressed. Problems of scale dependence also arise when measurements made at
a point are extrapolated to another scale. In section 2.4.1, the problem of scale
dependence of routing methods is addressed. This is followed in section 2.4.2 by a
brief discussion of current research in scale dependence of hydrologic parameters.
2.4.1 Scale Dependence of Routing Methods
Routing methods are transformations which are used in conjunction with
conceptual models to represent flow processes. If the routing methods are scale
independent, then routing parameters determined for a given model at one
resolution should be applicable at different resolutions. Methods can become
scale dependent if assumptions about the state of parameters become invalid or
the solution of the transformation equations becomes numerically unstable.
An example of scale dependence due to assumptions about the state of
routing parameters is found in storage routing methods in which linear
relationships are used to link control volume storage to inflow and outflow is
assumed. These relationships imply that within the control volume, velocity is
spatially and temporally invariant (Montes, 1998, pp.582). A similar assumption
is made in the simplification of the momentum equation in kinematic wave
routing. Over certain routing intervals, the magnitudes of the first two terms of the
momentum equation are small even for the variable velocity case (Miller and
24
Cunge, 1975). The magnitude of these terms may however become significant as
the routing interval and reach length increases, making the methods scale
dependent. Another illustration of this type of scale dependence is found in the
linear reservoir method. Inflows to the control volume are instantaneously
distributed throughout the reach thus ensuring uniform storage throughout the
reach. Discharge from the reach is then computed on the basis of the control
volume storage. At an hourly time interval, it is easy to envisage a 200m long
reservoir in which all input is instantaneously distributed throughout the reservoir.
It is not as realistic to envisage input being distributed instantaneously over a
1000 km reservoir in one hour simply because such a time period is too short for
such a redistribution to take place. A linear reservoir model is thus made scale
dependent by the underlying assumption of instantaneous inflow redistribution.
Scale dependence in routing may also be brought about by a numerically
unstable solution. Simplifying assumptions are made in solving the mathematical
equations describing flow. These assumptions result in a deviation of the estimate
from the real value. A solution is said to be numerically unstable one when the
magnitude of the deviation of the estimate from the real value is significant
compared with the magnitude of the real value itself. In the Muskingum method
for example, discharge from a control volume is predicted on the basis of
discharge during the preceding time period, inflow in the preceding time period
and inflow in the current time period. The simple integration technique used in
arriving at the solution was shown by Cunge to result in the elimination of second
order derivatives thus introducing an error in the solution (Montes, 1998, pp.572).
25
This error (similar in form to the diffusion term in the diffusion wave equation) is
responsible for the observed attenuation of the hydrograph when the Muskingum
method is used for routing flow. The magnitude of this error varies in proportion
to the square of the length of the reach. Applying the method sequentially over a
reach will therefore not result in the same discharge as a single application across
the entire reach. In effect, the Muskingum method is scale dependent because the
approximation made in deriving the method results in an error whose magnitude
is a function of length. As shown in the work of Cunge (1969) and Koussis
(1978), the magnitude of this error can in fact be computed and taken into
account. However, unless this correction is made, the Muskingum equation results
in a scale dependent model.
In continental scale routing, the use of scale independent methods allows
parameters determined from finer resolutions models to be incorporated into
coarser resolution models. Scale independent methods also serve to minimize the
propagation of systematic errors when routing is performed over large distances.
2.4.2 Scale Dependence of Routing Parameters
Another approach to the problem of scale dependence involves the use of
statistical methods to translate parameters measured at one location to another.
This approach, commonly known as scaling theory, has produced numerous
relationships many of which are discussed in Gupta et al (1986) and Sposito
(1999). The approach is based on the premise that there is an intrinsic relationship
between hydrologic parameters and the topographic conditions at the point of
26
measurement. Statistical methods are used to deduce this relationship from
measurements of the parameter and a corresponding topographic index. The
relationship is then expanded to cover all similar areas irrespective of spatial
resolution. The applicability of scaling theory is not addressed in this research
because the range of possible solutions is infinite when an intrinsic relationship is
deduced from two sets of data. The presence of measurement errors in limits the
reliability of the solution to the data set from which it is deduced. The addition of
more data sets simply increases the degrees of freedom in the solution.
27
CHAPTER 3
Data Development
3.1 DATA SOURCES
The Digital Chart of the World distributed by Environmental Systems
Research Institute (ESRI, 1993) includes digital maps of the world?s rivers. The
database was derived from operational navigation charts of Australia, Canada, the
United Kingdom and the United States. The rivers in this database were digitized
from 1:1,000,000 scale maps. The DCW represents the most complete coverage
of the rivers of the world currently available in digital format. A coarser
resolution coverage is also available in the Arc World dataset (ESRI, 1992) also
compiled and distributed by ESRI. The river coverage was digitized from
1:3,000,000 scale maps. Neither of these river coverages incorporates network
connectivity.
The River Nile in North-Eastern Africa is one of the most studied rivers in
the world. Flow gauge dating as far back as 620 AD have been collected in the
basin (Hurst and Phillips, 1938). A major effort to systematize and publish these
datasets was undertaken by a team led by Hurst. The results of these efforts are
published in a series of volumes entitled ?The Nile Basin? (Hurst and Phillips,
1931-). Flow data from these volumes as well as those reported in Shahin (1985)
are used in deriving flow parameters in this study.
28
30 arc-second Digital Elevation Models (DEM) produced by the EROS
Data Center of the US Geological Survey (USGS, 1996) are used to derive
hydrologic characteristics for both the source-to-sink and the HMS models. The
DEMs which were developed as part of the GTOPO30 mapping project are
currently the finest resolution global terrain representation available. The DEMs
can be downloaded from the USGS site, over the Internet at
?http://edcdaac.usgs.gov/gtopo30/gtopo30.html?
Runoff generated from a Global Circulation Model (GCM) is also
available. The Community Climate Model (CCM3) created by the National
Center for Atmospheric Research (NCAR) (Kiehl et al, 1996) was used to
simulate daily runoff generation for a period of 10 years (Branstetter and
Famiglietti, 1998). Runoff was generated over T42 grid cells with a resolution of
approximately 2.8 by 2.8 degrees. While the resulting data may not be
representative of actual values on the ground surface, they provide spatially
varying input that can be used for modeling. The use of GCM runoff also allows
routing to be performed in an environment similar to that encountered by global
climatologists and other users of continental scale routing models.
River discharge data is required for the calibration and verification of the
routing models. The Global River Discharge Database contains monthly flow
time-series from 949 runoff stations around the world (V?r?smarty et al, 1996).
The database was assembled by reviewing discharge of rivers from published data
29
sources. These data provide a basis for determining flow routing parameters such
as velocity and dispersion coefficient. The data are available over the Internet at
?http://www.rivdis.sr.unh.edu/?
3.2 TERRAIN DATABASE DEVELOPMENT
As stated in the section 1.2 of this dissertation, the development of a
database to support routing at the global scale is of prime importance. The utility
of a modeling system is determined not only by its ability to represent a
hydrologic system but also the ease with which such a representation can be
attained. At the global scale, data availability and computational efficiency
determine the ease of representation of a hydrologic system. A terrain analysis is
undertaken to develop a global database of terrain-derived parameters required for
continental scale hydrologic modeling. The database is used to parameterize two
continental scale routing models as described in sections 3.4 and 3.5 of this
dissertation. Drainage basins and modeling units are delineated from Digital
Elevation Models (DEMs) in accordance with the conceptual models defined in
those applications. DEMs are the most detailed description of surface topography
available in digital form. They are gridded representations of the land surface with
the ground surface elevation of each corresponding piece of land stored in a grid
cell. The determination of hydrologic modeling parameters from DEMs relies on
the simple premise that water flows downhill, and in the direction of steepest
decent. The parameters that influence the response of an area of land to input
runoff can hence be determined from surface elevation.
30
The first objective in GTOPO30 DEM processing is to ensure that all land
surface elevations are referenced to a single horizontal datum. In meeting this
objective, ocean cells and other water bodies are first removed from the DEM.
The USGS sets a datum of 65536 m for all cell located within water bodies and
for land surface cells with elevations below sea level (USGS, 1996). In addition,
ocean cells are assigned dummy elevation values of ?9999. The GTOPO30 DEM
was prepared for analysis by assigning null values to the ocean cells. Land surface
cells with elevation values greater than 32768 were restored to their true negative
elevations by resetting their datum to zero. The resulting land surface grid was
raised by 1000m to ensure a positive land surface elevation everywhere. Next, a
map projection was used to transform the land surface grid from the World
Geodetic Survey ellipsoid of 1984 (WGS84) used to represent the earth in the
original DEM onto a flat surface. The map projection selected for this study is the
Lambert Azimuthal Equal-Area projection used by the USGS in a similar
delineation (USGS, 1997). In this projection, different central longitudes and
latitudes are defined for each continent. All other projection parameters are kept
constant. A projection center which lies close to the centroid of each continent is
chosen in order to minimize areal distortion. The projection centers for the
various continents are given in Table 3-1 below.
31
Table 3-1: Projection Centers for Various Continents
Continent Central Longitude Central Latitude
Africa 2000 500
Asia 10000 4500
Australia 13500 -1500
Europe 2000 5500
NorthAmerica -10000 4500
SouthAmerica -6000 -1500
A sample projection file is also shown in Table 3-2 below. The parameters
of the output projection include the radius of the sphere of reference (6378137m),
the longitude of the center of the projection (20.0
o
), the latitude of the center of
the projection (5.0
o
), the false easting (0.0) and the false northing (0.0).
Table 3-2: Sample Projection File for Africa
INPUT
PROJECTION GEOGRAPHIC
UNITS DD
PARAMETERS
OUTPUT
PROJECTION LAMBERT_AZIMUTH
UNITS METERS
PARAMETERS
6378137.00000
2000
500
0.0
0.0
END
32
The second objective in DEM processing involves correcting erroneous
cell values to allow for the derivation of hydrologically correct flow paths. In
DEM processing, a pit is a cell or group of cells completely surrounded by higher
elevation cells. The delineation programs used in the analysis of these data cannot
define flow direction for the cells at the bottom of such pits. While naturally
occurring pits exist on the earth? surface, many of the pits in the DEM are the
byproduct of the DEM production process. Such artificial pits were located and
filled while retaining natural ones. This was done by first filling all pits in the
DEM. The filled areas were then isolated by subtracting the elevations of cells in
the unfilled DEM from those of the filled one. Natural sinks were identified by
comparing the shape of the filled area with those of naturally occurring lakes and
depression identified by visual inspection of areas maps. Sources of maps used
include the Times Atlas of the World (Bartholomew etal., 1975), the UNESCO
Atlas of the World Water Balance (Korzoun et al, 1978), Watersheds of the
World (Revenga etal,1998) and the Perry-Casta?eda Library (PCL) Map
Collection, available over the Internet at
?http://www.lib.utexas.edu/Libs/PCL/Map_collection/Map_collection.html?
After distinguishing between natural and artificial pits, a single null value
cell was placed at the lowest point of the natural pits in the unfilled DEM (Olivera
et al, 2000). When the resulting DEM was subsequently filled, natural pits were
preserved while artificial ones were eliminated. Cells in the catchment of the
unfilled natural pits (called inland catchments) are thus allowed to drain internally
33
in the direction of the null value cell. All other cells drain towards on the
continental margin. Figure 3-1 also shows the locations of the major inland
catchments identified around the world.
Figure 3-1: Inland Catchments of the World
Table 3-3 lists the names of the natural depressions identified for each
continent along with the location (in geographic coordinates) of the corresponding
null value cell. Due to the large spatial extent of many of these features, the
geographical locations given in the table are only approximate.
34
Table 3-3: Inland Catchments Identified
Continent Inland Catchments Longitude Latitude
Africa Lake Chad 14.45 21.09
Lake Eyasi 35.13 -3.59
El Wahat Kharga 31.36 22.87
Lake Turkana 36.43 2.91
El Djouf -5.83 21.09
Lake Rukwa 32.06 -7.83
Lake Baringo 36.04 1.75
Okavango Delta 25.84 -20.82
Asia Aral Sea 59.84 44.86
Lake Balkash 76.07 46.69
Oz Issyk Kul 77.39 42.41
Lop Nor 90.38 39.95
Kuruk Gol 92.11 42.70
Turfan Depression 89.30 42.68
Australia Lake Eyre 140.04 -28.55
Lake Frome 140.24 -29.47
Europe Caspian Sea 50.11 43.02
Dead Sea 37.45 31.33
North America The Great Salt Lake -118.02 39.80
Salton Sea -115.83 33.31
Saline Valley -112.54 41.11
Death Valley -117.83 36.70
South America Lake Poopo -67.49 -17.87
Salar de Arizaro -67.69 -24.85
Salar de Atacama -68.27 -23.53
Salina del Antofelia -68.03 -26.25
La Carachi Pampa -67.45 -26.54
35
The resulting DEMs were then used in the determination of flow direction
and other hydrologic parameters of interest. The GIS software used in this
analysis, Arc Info, defines a single flow direction for each cell based on the
direction of steepest descent from a given cell to one of the eight neighboring
cells. The underlying assumption in this analysis is that each grid cell discharges
all its flow to a single neighboring cell. Flow direction is denoted by a value
assigned to each grid cell. Figure 3-2 shows the cell values assigned for the
respective directions of flow.
Figure 3-2: Flow Direction Arrows
The flow direction grid is an important component of the database because
it is a required input for the derivation of other database components such as
watershed coverages, river networks and flow length grids. It must consequently
be checked for errors. One simple test performed for this data is the computation
of histograms of cell values. The presence of cell values other than those shown in
Figure 3-2 denotes the presence of sinks which must be removed by running the
Fill function on the corrected DEM again. Figure 3-3 shows the flow direction
grid derived for the continent of Africa.
36
Figure 3-3: Flow Direction Grid of Africa
37
Figure 3-4 (a): Histogram of Flow Directions for North America
Figure 3-4 (b): Histogram of Flow Directions for Africa
38
Figures 3-4 (a) and (b) show the histograms of the flow direction grid of
North America and Africa. The histograms indicate that the flow direction
algorithm tends to favor the cardinal directions (north, south, east and west) over
the diagonal directions (northeast, northwest, southeast and southwest). For the
entire dataset (of the whole world) 65.8% of grids cells had flow in a cardinal
direction as compared to 34.2% diagonals. This indicates that the flow direction
algorithm used in the model is biased in favor of flow through the cardinal
directions. No cell values other than those corresponding to the respective flow
directions in the eight direction pour point model shown in Figure 3-2 were
identified in the histograms of the various continents. This indicates that the flow
direction grids computed for the respective continents do not contain any sinks.
A flow accumulation grid is another basic database component derived
from the flow direction grid. Flow accumulation is measure of the number of
upstream cells that drain into a given grid cell. Since each cell in a flow direction
grid represents an equal area of land, the value of flow accumulation in any cell
can be converted to upstream drainage area by multiplying by the area of the cell.
If a weight grid containing a spatially distributed variable such as runoff is
supplied, performing flow accumulation results in a grid of the total quantity of
the variable upstream of each cell, a quantity that can be used for steady state
analysis. Figure 3-5 shows the flow accumulation grid of Africa derived in this
study. The darker areas in the figure are areas of high flow accumulation.
39
Figure 3-5: Flow Accumulation Grid of Africa
A key determinant of the response of a river basin is the relative distance
of various parts of basin from the basin outlet. This distance, often referred to as
the flow length, can also be determined from the flow direction grid. Arc Info
contains a function for determining flow length from a given cell to its most
upstream or downstream cell. The flow length function also includes a weight
option which can be used for aggregating quantities along the flow path to or from
a given cell. Figure 3-6 shows the flow length grid derived for Africa.
40
Figure 3-6: Flow Length Grid of Africa
3.3 PREPROCESSING FOR THE SOURCE TO SINK MODEL
The development of the database of hydrologic parameters was motivated
in part by need to route runoff generated by the land surface component of Global
Circulation Models to the ocean cells used in the ocean modeling component. In
this section, a preprocessing procedure which allows a surface water routing
model to interact with both a soil water balance model and an ocean model is
41
described. The procedure was applied to each of the earth?s continents to justify
its label as a continental scale model. The Land Surface Model, LSM (Bonan,
1996) of the Climate Community Model, CCM3 (Kiehl et al, 1998), developed by
the National Center for Atmospheric Research was used for runoff generation.
The ocean modeling component of CCM3 was not used directly in this study but
discharge outlets were defined in a manner consistent with the modeling units
used in that component. The goal in this analysis was to demonstrate a procedure
for implementing continental scale routing. Hence the parameters and inflows
used in this section were not expected to match any observed flows. The routing
procedure adopted in the study is described in chapter 4, and a copy of the routine
code is included in appendix 3 of this dissertation. Testing of the routing model
with various scenarios of routing parameters is described in chapters 5. The
procedure defined in this study can be used with other continental scale runoff
generation and discharge receiving models after minor adjustments for input data
format.
In this application, the source to sink approach (described in section 2.1.3)
is preferable to the watershed and cell to cell approaches (described in sections
2.1.1 and 2.1.2 respectively) because of its ability to route flow directly to
required locations. Continental scale source to sink models were consequently
developed for each continent excluding Antarctica which has no rivers. The
processing required to create these models is described in this section.
Preprocessing programs developed (by the present author and Francisco Olivera)
42
in the Arc Macro Language (AML) of the GIS software, Arc Info, are included in
Appendix 1 of this dissertation.
Figure 3-7: Components of a Source-To-Sink Model
The conceptual model for source to sink routing consists of three
components namely a source which receives input, a sink which serves as a
discharge location and a flow path which links the source to the sink. Figure 3-7
shows a source to sink model consisting of 8 sources linked to a single sink by
43
flow paths. In this study, sources, sinks and flow paths are defined to allow for
easy interaction with the modeling units defined in the respective land surface and
ocean modeling components of the CCM3. Each source is linked is linked to a
single LSM modeling unit from which it receives input runoff and to a single
Ocean model unit to which it discharges its flow. In other words, there must be a
one to one relation between sources and sinks and between sources and runoff
generating units. The process of establishing these one to one relations is
described below.
A line coverage consisting of the continental margin and the outline of
inland catchments was first created from filled DEM (described in section 3.2).
The coverage was created by identifying all grid cells attached to at least one null
value cell, putting them into a separate grid and vectorizing them. A second
coverage was created by generating a vector mesh of comparable resolution to the
Ocean Model units. For this study, a 3 x 3 degree mesh was used to represent the
ocean modeling units. The two coverages were intersected to break up the
continental margin into a series of small segments each wholly contained within a
single ocean modeling unit. The individual segments were then converted back
into raster format using their unique segment numbers as value fields. The grid
cells thus created were used as sinks for the division of the land surface into
drainage basins. Figure 3-8 shows the sinks defined for the continent of Africa.
Each mesh spacing corresponds to an ocean modeling unit. The gray boxes along
the continental margin are ocean modeling units which receive their inflow
directly from the surface water routing model. The sections of the continental
44
margin within each receiving unit serve as sinks. Inland catchments are similarly
represented by gray boxes located inside the continental margin. All other land
surfaces are represented by a light shade of gray.
Figure 3-8: Sink Locations for the Continent of Africa
The boundaries of the drainage basins serve to demarcate areas within
which sources associated with a particular sink can be defined. In this study, grid
45
cells identified as sinks were in conjunction with the 1 by 1 kilometer flow
direction grids (described in section 3.2) to delineate the drainage basins. The
resulting drainage basins delineated for the various continents are shown in
Figures 3-9 (a) to (f). Each of the resulting basins bears the id-number of the
ocean or inland sink to which its drainage flows. All discharges generated within
it can hence be linked to the correct ocean modeling unit. The continental divide
between Europe and Asia was established by delineating all basins draining into
the inland catchment formed by the Caspian Sea as part of Europe. All basins to
the east of this were delineated as part of Asia. This results in the definition of
Europe that is much larger than dictated by the political boundaries. It however
reduces the computer requirements for modeling the Asian continent which is the
world?s largest with an area of over 50 million km
2
. This is far more surface area
than that of any other continent and would have strained the processing capacity
of the computers available for this study. By assigning the Caspian Sea drainage
basin to Europe, the model of Asia was reduced to a more manageable size with
approximately 38 million grid cells.
46
Figure 3-9 (a): Drainage Basins of North America
Figure 3-9 (b): Drainage Basins of Africa
47
Figure 3-9 (c): Drainage Basins of South America
Figure 3-9 (d): Drainage Basins of Australia
48
Figure 3-9 (e): Drainage Basins of Europe
Figure 3-9 (f): Drainage Basins of Asia
49
With the link between the drainage basins and ocean models established,
the next step in preprocessing involves the establishment of a one to one relation
between the land surface and the LSM modeling in which runoff is generated. For
this purpose, the drainage basins must be divided into smaller areas such that each
small area fits exactly into a single LSM modeling unit. Runoff in the LSM can be
generated at various resolutions. The runoff used in this particular study was
generated on a Gaussian grid with a horizontal resolution of T42 (approximate 2.8
by 2.8 degrees of latitude and longitude). T42 stands for spectral (spherical
harmonic basis functions) triangulation 42. It is a mesh defined in geographic
coordinates with uniform mesh spacing (in degrees) along the longitudinal
direction but slightly non-uniform mesh spacing in the latitudinal direction
(Hortal, 1991). The non-uniformity in the latitudinal direction is to account for
differences between the spherical or geocentric coordinate system used in
atmospheric modeling and the ellipsoidal or geodetic coordinate system used in
land surface modeling. Reed and Maidment (1995) describe the differences
between the two coordinate systems and outline a procedure for accounting for
these differences in hydrologic modeling. The slightly non-uniform T42 mesh
translates to a uniform mesh in spherical coordinates thus ensuring a one to one
relation between the land surface and atmospheric modeling units. The parameters
of the T42 mesh are provided in Appendix 2 of this dissertation. During each
routing time step, the LSM generates a single value of runoff in (mm/s) for each
50
T42 unit. The drainage basin coverage of each continent was intersected with the
T42 mesh to create a new polygon coverage. Each polygon in the new coverage is
associated with a single LSM modeling unit from which it receives its input
runoff and a single ocean model unit to which it discharges its output flow. Figure
3-10 shows the distribution of these polygons on the Africa continent. The
intersection of LSM modeling units (represented by a white mesh on the land
surface) and the drainage basin boundaries (represented by a thick black line)
results in a model with a few sources defined in each drainage basin.
Figure 3-10: Linking Modeling Units for the African Continent
51
A response function can be defined to describe the transportation of the
runoff from each of the sources in Figure 3-10 to its associated ocean cell. Such a
response would imply that flows from all areas within each source travel through
a single flow path and arrive at the basin outlet at the same time. A single one of
these sources covers an area of about 98000 km
2
. This is equivalent to an area
roughly 310 km wide and 310 km long. The assumption of a single flow path for
the entire polygon is consequently not realistic, and the resulting model would not
capture the features of the discharge hydrograph effectively. Smaller source areas
are required to allow for better representation of the spatial variability of flow
length. A 0.5 by 0.5 degree mesh was intersected with the polygon coverage
formed by the intersection of the T42 mesh and the basin coverage to create
smaller sources for the source to sink model. The use of a 0.5-degree mesh
ensured a maximum source area of approximately 3600 km
2
which corresponds to
an area about 60 by 60 km wide at the equator. While this choice of resolution
was made in line with the expected resolution of future GCMs, a finer or coarser
resolution could just as easily have been selected. The implications of the choice
of spatial resolution are explored in section 5.2.2 of this dissertation. The sources
defined for each of the continents are shown in Figures 3-11 (a) to (f) below.
52
Figure 3-11 (a): Parameterized Sources of North America
Figure 3-11 (b): Parameterized Sources of Africa
53
Figure 3-11 (c): Parameterized Sources of South America
Figure 3-11 (d): Parameterized Sources of Australia
54
Figure 3-11 (e): Parameterized Sources of Europe
Figure 3-11 (f): Parameterized Sources in Asia
55
3.4 PREPROCESSING FOR THE HYDROLOGIC MODELING SYSTEM
An important aspect of this research is the validation of the routing
process. The typical approach to model validation is to run the model with a set of
input runoff values for which there is a corresponding set of observed discharges.
The discharges produced by the model can then be compared with the observed
values. This approach would not suitable for validating a routing model which is
only concerned with the horizontal transport component of the hydrologic cycle.
It would be impossible to determine which errors to attribute to the runoff
generation scheme used to produce input data and which to attribute to the routing
scheme. A more representative approach would be to compare the discharges
from the routing model to those predicted by another routing model given the
same input data and corresponding routing parameters. The choice of the second
routing model is obviously an important factor in such a comparison since it
serves as a baseline. For this study, the Hydrologic Modeling System (HMS)
model (Peters, 1998) developed by the Hydrologic Engineering Center of the US
Army was selected because of its widespread use and acceptance in the
engineering community. The HMS model has the additional advantage that it
takes into account spatially distributed parameters which can be derived from the
same digital elevation models (DEMs). Since both the HMS and the STS model
can be parameterized from the same DEM, a true comparison can be made of
their respective routing methods.
56
As with the STS model, an HMS model of the continental of Africa was
developed to justify its label as a continental scale model. The preprocessing
required to define hydrologic elements and the element connectivity in HMS is
described below. The GIS preprocessor (Olivera and Maidment, 1999) developed
at the Center for Research in Water Resources of the University of Texas at
Austin, CRWR-Prepro, is used in the delineation. The filled DEM, flow direction
and flow accumulation grids described in section 3.2 of this research were used in
the delineation. The first objective in the HMS delineation was the definition of
all the natural hydrologic elements namely streams, outlets and watersheds.
Streams can derived from DEMs by specifying a minimum drainage area required
to constitute a stream. From previous experience working with 30 arc-second
DEMs, a threshold area of 1000 km
2
or higher results in a reasonable
representation of the river networks. A 10000 km
2
threshold was selected for the
continental scale model runs. A finer threshold of 1000 km
2
was however
specified for a secondary model of the Congo basin created for more detailed
analysis. For the 30 arc-second DEMs, the threshold value also translates to the
value of flow accumulation since each cell has an area of 1 km
2
. Hence all cells
with flow accumulation values of 10000 or higher were labeled as being stream
cells.
The stream cells were then grouped into reaches or links and each link was
assigned a unique label to differentiate it from other links. 1791 links were
identified in the Congo basin during this delineation. From the definition of flow
accumulation, the outlet cell in each link is also the cell with the highest flow
57
accumulation. CRWR-Prepro contains a function which identifies the outlet cell
for each link based on the flow accumulation value of all stream cells. The
preprocessor also contains an option which allows the user to specify additional
outlets though this option was not exercised in this delineation.
Outlets serve an important role in delineation because they provide
convenient starting points for watersheds. By dividing a drainage area into
watersheds an HMS modeler can determine over which areas to apply runoff
transformation approaches for the overland portion of flow and where to
commence in-stream routing approaches such as the Muskingum method. A
watershed grid was delineated from the outlet points and flow direction grid, and
a vector representation of the entire basin was obtained by converting the
respective watershed and stream grids to vector coverages. In vectorizing the
streams, a minimum stream length of 710 m was specified to prevent the
occurrence of single-cell streams. This value was determined as the approximate
length of a single 1 km by 1 km cell with flow through one of its corners. A
higher threshold may be specified as long as it is understood that the removal of a
given length of stream is justifiable only if routing flow over the specified length
does not result in any significant difference between the input and output. The
vectorized stream and watershed coverages delineated in Africa at a 10,000 km
2
threshold are shown in Figures 3-12 and 3-13, respectively.
58
Figure 3-12: Delineated Rivers Coverage of Africa
Figure 3-13: Delineated Watersheds of Africa
59
The next step in the preprocessing is the calculation of hydrologic
parameters. In this step, CRWR-Prepro uses the previously derived flow
accumulation grids, as well as flow length grids and other routing parameters to
compute the numerical characteristic of each hydrologic element. For watershed
runoff transformation, the initial and constant losses were set at zero for all the
watersheds, and a single lag time was used to translate flow generated within the
watershed to its outlet. The routing interval was specified as one day in line with
the frequency of runoff generation, and a flow velocity of 0.3 m/s was used.
These specifications imply that when input runoff is applied to the watershed, an
initial portion of the input (expressed in mm) is lost through evaporation and
infiltration. The remainder of the runoff is transported to the outlet by pure
translation with a constant loss rate in mm/hr. The travel time to the outlet is
determined as the larger of 0.6 times the lag time along the longest flow path in
the watershed and 3.5 days. The lag time along the longest flow path is in turn
determined by finding the cell with the longest downstream flow length,
calculating the distance from that cell to the outlet and dividing the result by the
watershed velocity.
Stream routing parameters were also determined for the Congo river
reaches. For this determination, the Muskingum method was specified with a
Muskingum x of 0.3 and an in-stream flow velocity of 0.3 m/s. This specification
results in reaches with a lag time greater than the routing interval are routed using
the Muskingum method while the rest were routed by pure translation. The
Muskingum K and n parameters for each element are computed based on the flow
60
velocity and Muskingum x provided by the user. A full discussion of the
equations and routing methods used in the preprocessing and flow routing
processes is given in section 4.2.3. They are only mentioned here to highlight the
capabilities of the preprocessor in calculating routing parameters for each model
element. The HMS model created of Africa consists of 1551 watersheds, 701 river
reaches, 197 sinks and 693 junctions. Figure 3-14 shows the line diagram of
hydrologic elements derived for Africa using CRWR-Prepro.
Figure 3-14: Line Diagram of Hydrologic Element Connectivity
The final step in preprocessing involves the writing of a basin file which
can be interpreted by HMS. CRWR-Prepro performs this transcription by
61
converting the parameterized watershed and river network coverages into a line
diagram with the connectivity between the elements outlined (Hellweger and
Maidment, 1996). It then creates a text file listing each hydrologic element, its
attributes and the element to which it is connected at its downstream end. The
resulting file was imported into HMS and used to generate a basin model such as
the one shown in Figure 3-15 below.
Figure 3-15: HMS Model Representation of Africa
A number of mistakes were encountered when the basin file was first
imported into HMS. These mistakes manifested themselves either during the
62
import stage or when the imported basin model is used to perform a simulation.
The type of errors encountered most frequently were those resulting for
erroneously created sources, and those resulting from parameters in excess of the
permitted range. Erroneous sources arise due to the presence of sliver polygons in
the watershed coverage. Sliver polygons are small polygons that are erroneously
created when the topology of a coverage is rebuilt. They are often not visible at a
glance but they are identified by CRWR-Prepro as sources because they appear to
be small closed polygons which are not attached to an outlet. The presence of
such sliver polygons is easy to detect when HMS is run because they result in
sources with no associated inflow. These sources can be removed by deleting
them from the basin file along with the reaches that link them to the next
downstream element, usually a junction.
The problem of parameters falling outside their permissible range was also
encountered when the HMS basin file was used for flow routing. HMS has a
specified maximum Muskingum K value of 150 hours. An error message is
generated when HMS encounters a reach with a K value in excess of this value
during routing. To correct this error, such reaches were each subdivided into 2 or
more equal reaches. The total lag time of the original reach was subdivided
among the new reaches, and new values of Muskingum n were also computed for
each new reach. Network connectivity was restored by linking the new reaches
sequentially to each other and connecting the first and last reaches in the sequence
to the original upstream and downstream elements, respectively. To preserve the
labeling structure used in the original delineation, new reaches were labeled by
63
introducing an additional digital after the last digit of the id-number of the old
reach. These changes were made by manually editing the basin file.
The creation and parameterization of a HMS model of Africa served to
illustrate that a continental scale watershed based routing model could be created
and parameterized based on available DEM data. Having previously created a
continental scale source to sink model based on the same DEM data, two
independent and yet consistent representations of the land surface were obtained.
Key features of basin responses computed from the two models can thus be
compared. The methods used for computing hydrographs in the respective models
are described in chapter 4, and in chapter 5, various simulations are performed to
determine the effect of various modeling parameters on the basin responses.
64
CHAPTER 4
Methodology
4.1 THE SOURCE TO SINK MODELING APPROACH
A framework for implementing the source to sink model at the continental
scale is described in this section. The assumptions made in the development of the
framework, the model components and the routing methodology are described
along with methods of estimating routing parameters. In these discussions, the
concept of a hydrologic system is adopted to explain the functioning of the
framework. A similar description is presented of the watershed based modeling
system, HMS, in section 4.2 of this chapter.
4.1.1 STS Model Assumptions
A hydrologic system can be classified by the nature of its transformation
as linear or non-linear. The system is said to be linear if the transformation meets
the criteria for the application of superpositioning otherwise it is non-linear.
Dooge (1973, pp.18) gave the mathematically interpretation of linearity as
follows:
?If a transformation X
1
(t) results in Y
1
(t) and X
2
(t) results in Y
2
(t) then
for a linear system, X
1
(t) + X
2
(t) results in Y
1
(t) + Y
2
(t)?
65
An assumption of linearity is made in the Source to Sink routing model. In
other words, a decomposition of the inputs at the source into a series of smaller
segments results in an identical output at the sink as would result from a routing
the lumped input. Olivera (1996) outlines three major features that characterize
linear hydrologic systems. These include longitudinal decomposability which
allows rivers to be split into smaller segments and modeled separately, areal
decomposability which allows basins to be subdivided and modeled as a series of
separate watersheds, and superpositioning of subarea responses which allows the
results of simulations in different watersheds to be aggregated
Dooge (1973) presents a rather complete review of linear systems theory
and its application in hydrology. In this review, he demonstrates that the
assumption of linearity results in some distinct advantages in hydrologic
modeling. A linear transformation allows the differential equations describing a
hydrologic process to be reduced to a set of algebraic equations which can be
solved by serial substitution. The reduction may be achieved using numerical
approximations or mathematical transformations. The assumption of a linear
system also allows for the modeling of processes that have not been described
with mathematical equations. If, for a set of system inputs, a corresponding set of
outputs is available, then the mathematical transformation can be derived. The
inputs and outputs can be broken down into a series of segments of equal
duration. The system transformation can be identified by matrix operations or
optimization techniques as the transformation that maps the inputs onto the
outputs. This type of transformation identification is only possible if all elements
66
of a system are linear. The introduction into a river network of a controlled
reservoir, for example, would render the whole network non-linear, making the
application of linear systems invalid.
The assumption of linearity has been widely used in hydrology and it
underlies many routing methods such as the unit hydrograph, the Muskingum
method and the linear solutions of the St. Venant equations. However, not all
routing methods used in hydrology assume linearity. The form of the dynamic
wave equation used in overland flow routing for example is non-linear. Also, the
use of a storage function in which storage is a function of inflow or outflow raised
to a power other than one results in a non-linear system. Hence linearity must be a
stated assumption if superpositioning is to be applied in the source-to-sink routing
method.
The second important assumption made in the source-to-sink approach is
that transformation parameters relating system input to output are time invariant.
In other words, a given system input would result in the same output irrespective
of the time of application of the input. For a source to sink model, this implies that
the physical location of the flow path linking a source to a sink does not vary with
time. In addition, the parameters at each location along that flow path remain
constant throughout the routing period. According to Dooge (1973, pp.18), time
invariance can be expressed mathematically as follows:
?If a transformation X(t) results in Y(t) then
for a time invariant system, X(t+r) results in Y(t+r)
where t and t+r are points in time?
67
This assumption is essential to the definition of a system transformation or
instantaneous unit hydrograph which can be applied for subsequent events. If a
system transformation changes with time then the outputs due to two events
occurring within a few minutes of each other cannot be superimposed because
they undergo different transformations. Hydrologic systems rely on the time
invariance of parameters to predict the output of future events. If a system
transformation is time variant then no relationships can be deduced from past
events, and consequently, no predictions can be made for the output of future
events based on what happened in the past. Implicit and explicit methods used to
solve the full St. Venant equations are examples of time variant solutions.
4.1.2 STS Model Components
In source-to-sink modeling, a river basin is viewed as a series of sources
each defined and parameterized not in relation to neighboring cells but in relation
to an observation point or sink located further downstream. Sources interact with
sinks through a response function. The model is not fully spatially distributed in
the sense that the location of water is not tracked throughout the system. Rather,
water is provided as input at the source and only measured as it flows into the
sink. The flow path linking each source to its sink serves to trace the control
volume within which flow is routed. If the location of a sink changes, a new
response function must be defined for each source since the definition of the flow
path is also altered. Figure 4-1 provides an illustration of the major components
68
defined for the STS model implemented for routing Land Surface Model (LSM)
runoff to the continental margin in this study.
Figure 4-1: Diagram of STS Model Components
The source file contains a listing of all the sources within a study area
along with characteristics of their associated flow paths. The parameters listed for
each source include a unique identification number (source-id), its area (source
area), the identification number of its associated sink (sink-id), the length of flow
path from source to sink (flow length), the effective velocity of flow from source
to sink (velocity), the effective dispersion coefficient along the flow path
(dispersion coefficient), the effective first order loss coefficient along the flow
path (loss coefficient) and a key field linking each source to a unique modeling
unit in the runoff generation model (runoff key field). The sink file contains a
69
listing of the all the sinks within the study area by the same sink-id used in the
source file. The sinks listed may include both inland catchments and river basin
outlets located at the continental margin. The parameter file contains simulation
specific information including the routing interval (in seconds), the number of
input events, the number of sinks and the number of sources. The routing code for
the STS model was developed in the FORTRAN programming language.
Separate versions of the routing code are maintained for pure translation and
dispersed model runs. The details of the routing process are described in section
4.1.3. The flow files consist of the input runoff files (Branstetter and Famiglietti,
1998) generated by the Land Surface Model (LSM) and the output files
containing routed flows generated by the routing code. The LSM runoff is
reported for each T42 cell in millimeters per second. A single runoff file is
generated during each routing interval. In this study, the runoff files were named
in terms of the routing interval during which they were generated so that
?h0001.txt?, for example, is the runoff file generated in the first routing interval.
For computing basin response, a single runoff file with a uniform runoff input was
provided. The output files generated by the routing code in this study include a
file containing the routed flow time series (in cubic meters per second) for each
sink, a file containing the total input runoff (in cubic meters) for each drainage
basin for each time period, and a file containing the response function along the
flow path (in terms of fraction of inflow) for each source. The actual routing
procedure is described in section 4.1.3 followed by methods of deriving routing
parameters in sections 4.1.4.
70
4.1.3 STS Routing Methodology
Routing is performed in a space first, time second order. Hence given a set
of spatially distributed inputs, the Fortran program routes each input for all the
sources before going to the next input. The results of routing individual input sets
are superimposed at the sinks. The values of discharge at the outlet are not correct
until all the sets of inputs have been routed. Input runoff for each routing time
step is generated as an ASCII file. The link between records in the file and
sources in the basin file of the source to sink model is established through a key
field that is transferred from the T42 mesh to the source polygon coverage when
the two are intersected during preprocessing as described section 3.3 of this
dissertation. For this study, the key field was created from the X and Y coordinate
indices (listed in Appendix 2 of this dissertation) using equation <number
below>. The origin of coordinate indices (1,1) is defined at the southwestern
corner of the T42 mesh at (0
o
, -87.864
o
) in geographic coordinates.
Equation 4-1 )129(1000 XYK ?+=
where K is the value of the key field
Y is the Y coordinate index with values ranging from 1 to 64
X is the X coordinate index with values ranging from 1 to 128
Flow in the STS model is routed directly from a source to its associated
sink along a flow path. A response function is defined for each flow path which
71
describes the shape of the hydrograph at the sink given an instantaneous input at
the source. The response function includes a measure of the travel time, storage
effects and losses due to evaporation and infiltration along the flow path. No
additional flow is allowed to enter the flow path between the source and the sink
since a separate source is defined for each geographic location. Even though, in
reality, flows from the various sources may pass though the same channel, flows
from each source can be treated as separate entities because of the assumption of
linearity. Travel time is computed by dividing the downstream flow length by the
translational velocity. The term translational velocity is synonymous with wave
celerity and is used to differentiate it from flow velocity at a given cross section.
If different velocity zones exist, their combined effect along the flow path can be
determined by expressing velocity in terms of flow time through each segment
and adding up to get the total flow time for the entire reach. Velocities along a
flow path cannot be averaged since velocity is a vector quantity without
accounting for length and direction of the velocity vector. On the other hand, flow
time is a scalar quantity and can be summed up irrespective of direction. Olivera
and Maidment (1996, 1999) define a GIS based approach for performing the
summation of flow time in a spatially distributed field using a weighted flow
length scheme. For a flow path composed of n segments or grid cells, the total
flow time through the flow path is given by
Equation 4-2
?
=
=
n
i i
i
v
l
t
1
72
where l
i
is the length of each flow path segment or grid cell
v
i
is the flow velocity through the segment or grid cell
If a weight grid of
i
v
1
is supplied, the weighted flow length function results in a
spatially distributed grid of flow time from each cell to the outlet.
Storage effects that describe the deformation of an input pulse as it
travels downstream along a channel must similarly be taken into account. The
storage effect has been attributed to the shear forces between the sides of the
channel and the moving water (Fischer et al, 1979). These forces result in a
velocity profile across the stream section such that particles closest to the sides of
the channel are transported at a much lower velocity than those closer to the
center of the stream. The result is that input water is spread out longitudinally
along the river channel as it travels downstream. Additional spreading out results
from the damping effect of irregular channel cross sections (Henderson, 1966).
Such irregularities result in the storage of water along the banks thus delaying the
arrival of some portions of the input at the downstream location. Several
approaches for dealing with storage are documented in the literature. Conceptual
models such as the linear reservoir and the cascade of linear reservoirs (described
in the section 2.2 of this dissertation) have been used for this purpose. Linear
reservoir models describe the storage effects between two interconnected
elements. In the source-to-sink method, the sink is remote from and usually
unconnected to the source. Hence a linear reservoir cannot be applied directly.
The Lag and Route method uses pure translation to account for advection to a
remote outlet followed by a linear reservoir or cascade of linear reservoirs
73
(Maidment, 1993, pp.26.10) to account for storage effects. Both the cascade of
linear reservoirs and the lag and route methods are classified as advection
methods with approximate dispersion parameters as discussed in section 2.2.2 of
this dissertation. While these approaches may yield reliable results in large scale
routing (Dooge, 1973, pp.257), there is only an approximate physical basis for
their use. Methods that derive directly from the St. Venant equations are
preferable because they incorporate a full description of the dispersion process.
Maintaining the assumption of linearity and time invariance made earlier,
the momentum equation
Equation 4-3 0=
?
?
?
?
?
?
+
?
?
+
?
?
+
?
?
fS
x
h
g
x
V
V
t
V
can be reduced to
Equation 4-4 0=+
?
?
fS
x
h
by assuming a constant velocity.
Equation 4-5 0=
?
?
=
?
?
x
V
t
V
The impulse response function resulting from this simplification was
shown by Hayami in 1951 (Montes, 1998) to take the following form
Equation 4-6
()
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
Dt
xVt
PDt
x
txh
4
exp
4
),(
2
3
74
where
D is the dispersion coefficient in m
2
/s
V is the flow translational velocity in m/s
x is the distance from the source in meters
t is the time after input in seconds
P is the numerical constant ? = 3.14159
This equation is referred to as the passage times distribution. It can be
convolved with the input runoff at a source to generate flow at the sink. The flow
resulting from numerous sources can be superimposed at an associated sink such
that for any time, t, the total discharge, Q, at the sink is given by
Equation 4-7
?
=
=
n
i
i
i txhIQ
1
),(*
where n is the number of sources associated with the sink
iI is the input runoff in modeling unit i at time t
),( txh is the response of element i at the sink and
* is the convolution integral
Assuming a first order decay process for translational losses, an additional
term can be introduced into the exponent of the response function (Maidment, ed.,
1993 , pp.14.18) such that
75
Equation 4-8
()
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
= Kt
Dt
xVt
PDt
x
txh
4
exp
4
),(
2
3
where K is a positive, first order decay constant with units of (second)
-1
The exponential term in the first passage time distribution ensures that the
response asymptotically approaches zero without ever touching it. This results in
a long time series the latter terms of which are close to but greater than zero. The
small terms cannot be ignored because they would result in a mass balance
problem which may become significant for long term modeling applications. In
order to apply the first passage times distribution in a mass conservative manner,
it must be applied over a limited range of time steps with normalized coefficients
of the impulse response functions which sum up to one for unit input. Fischer et al
(1979) suggests that for simulations involving the first passage times distribution,
the significant width of the dispersed cloud can be assumed to be 4 standard
deviations. The width of the impulse response function, W, is consequently given
by Equation 4-9 below.
Equation 4-9 DtW 32=
where D is the dispersion coefficient in m
2
/s
t is the time after input in seconds
In this application, the width of the response function is used to estimate
the number of time intervals over which an input is distributed when it arrives at
76
the outlet. By ensuring that the total area under the response over this width is
equal to 1, the first passage times distribution is applied in a mass conservative
fashion without taking into account the endless tail that results from the
exponential function asymptotically approaching zero.
The assumption of a constant velocity, V, and dispersion coefficient, D,
used to simplify the momentum equation does not limit the utility of the model in
areas where D and V are not constant. Olivera (1996, 1999) demonstrates the use
of the flow length function in estimating a single dispersion parameter in place of
a spatially distributed field of coefficients. In his approach, the dispersion
coefficient in the first passage times distribution is replaced with the flow path
Peclet number such that the response function becomes:
Equation 4-10
?
?
?
?
?
?
?
?
?
?
=
ii
i
ii
i
tt
tt
ttPt
tu
/)/(4
)]/(1[
exp
/)/(2
1
)(
2
where )(tu
j
is the response at the outlet of flow path i
i
t is total flow time along flow path I
P is the numerical constant ? = 3.14159
i
? is the Peclet number along flow path i with j segments and is
computed as follows:
Equation 4-11
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=?
j
j
j
j
j
j
j
i
x
v
D
x
v
3
2
1
77
The use of the Peclet number however introduces another parameter
whose magnitude does not lend itself to direct interpretation. From McCuen
(1998), the total variance along a flow path can be determined as the sum of
variance along the component segments. This property derives from the fact that
variance is a scalar quantity. For a given flow path, i, made up n segments (j = 1
to n), the net variance is given by Equation 4-12.
Equation 4-12
?
=
?
?
?
?
?
?
?
?
=
n
j
j
j
j
i
L
v
D
iance
1
3
2var
In this study, a single value of the respective diffusion coefficients and
flow velocities are computed for each flow path based on the scalar nature of the
parameters. The value of single flow velocity that results in the same travel time
as the non-uniform segment velocity along a given flow path i is computed from
Equation 4-13 as follows:
Equation 4-13
?
=
?
?
?
?
?
?
?
?
=
n
j j
j
ii
v
L
Lv
1
where Li is the length of the single reach
This is in effect a division of the flow path length by the total flow time.
Similarly, a single value of dispersion coefficient which results in the same flow
path variance can be computed for each flow path from Equation 4-14.
Equation 4-14
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=
?
=
3
1
3
i
i
n
j
j
j
j
i
v
L
L
v
D
D
78
Computing the flow path parameters in this way allows the first passage
times distribution to be used in its original form thus eliminating the need to
introduce another parameter such as the Peclet number.
4.1.4 Estimating STS Routing Parameters
In addition to the distance of the source from the sink, the three main
parameters that must be estimated for a source to sink model are the diffusion
coefficient, D, the flow velocity, V, and the loss coefficient, K. The physical
meaning of these parameters is described in this section. This is followed by a
description of the method adopted for estimating the three parameters from
observed flow data.
The dispersion coefficient is a measure of the amount of stretching that an
input undergoes. In many applications in water quality modeling, the dispersion
coefficient is derived from channel properties (Fischer, 1967). The dispersion
coefficient used in the derivation of the diffusion wave method is a compound D
which incorporates the effects of both off-stream storage and spreading due to
friction between flowing water and channel surface, as shown in Equation 4-15.
Equation 4-15
WF
DDD +=
where D is the compound dispersion coefficient in m
2
/s
D
F
is the dispersion coefficient due to friction in m
2
/s
D
W
is the dispersion coefficient due to off-stream storage in m
2
/s
79
This compound dispersion coefficient cannot therefore be derived from
channel section properties measured a specific location as this would exclude the
cumulative effect of off-stream storage (Henderson, 1969). It must therefore be
derived from flow measurements. Dispersion coefficients compiled by Fischer et
al (1979) and Seo and Cheong (1998) may be used as first estimates. Empirical
functions such as Equation 4-16 suggested by Fischer (1967) for computing
longitudinal dispersion in natural streams may also be used as a first estimate if
river cross section data is available.
Equation 4-16
gRSd
wV
D
F
22
01.0
=
where D
F
is the dispersion coefficient due to friction in m
2
/s
V is the average stream velocity in m/s
w is the width of the stream in m
d is the gravitational acceleration in m/s
2
R is the average depth of the stream in m
g is the hydraulic radius in m
S is the friction slope expressed as a decimal
However, these equations require the use of a hydraulic radius which in turn
requires the modeler to make an assumption about the quantity of flow. Also, they
only take into account the dispersion due to friction.
The diffusion wave method describes the propagation of a flood wave
downstream. The flow translation velocity used in the method is thus a measure of
wave celerity not in-stream flow velocity at a location. It is a measure of how
80
quickly a flood wave is propagated through the flow channel. The wave celerity
for a wide rectangular channel is 1.5 or 1.667 times as high as in-stream flow
velocity depending on whether one adopts Chezy?s equation (Henderson, 1969) or
Manning?s equation (Chow et al, 1988) as a starting point. As with the diffusion
coefficient, flow translation velocity cannot be determined from equations which
include hydraulic radius as these require an assumption about the quantity of
flow. Approaches involving the replacement of the hydraulic radius with a steady
parameter such as drainage area (Maidment et al, 1996) hold promise. However,
in the opinion of this author, there is insufficient verification of the approach with
velocity data to justify their use globally. Empirical functions proposed by
Leopold (1953) and Carlston (1969) have been used to determine flow velocity
(V?r?smarty et al, 1989). However, these equations were intended to define the
variation of flow velocity with river stage at a given location. Their utility in
characterizing flow velocity over a general domain is thus brought into question.
A method of deriving velocity from observed flow measurements is proposed
later on in this section.
As water flows through the landscape, it incurs losses in the form of
seepage into the ground and evaporation into the subsurface. As shown in section
4.1.3, a first order loss function can be included in the first passage times
distribution. The first order loss coefficient, K, can be expressed either in terms of
fraction lost per unit of travel time or fraction lost per unit of travel distance. The
former is sensitive to changes in flow velocity while the latter is not. It is
consequently more convenient to adopt the fraction lost per unit of travel distance
81
approach when performing calibrations since a mass balance performed at the
beginning of the calibration remains valid at the end of the run irrespective of
changes in velocity. The exponential function associated with first order losses
determines the fraction of input flow that makes it through the flow path to the
sink. The presence of side flow and spring flow in natural flow systems makes it
difficult to define a physical basis for the computation of the loss coefficient. It
must be determined from observed flow data.
After reviewing several calibration approaches, the method of moments
(Dooge, 1973, pp.254) was adopted. The method is based on the premise that if
the application of a given transformation to one time series yields another then the
difference between scalar moments of the two time series is equal to the moment
of the transformation. The first moment of a flow time series taken about the
origin is a scalar quantity, and it corresponds to the flow time. The second
moment taken about the mean is a measure of variance, and dispersion coefficient
can be calculated from the variance using Equation 4-12. Consequently, given
two flow time series from gauging stations located a given distance apart on the
same river reach, the flow velocity and dispersion coefficient of the first passage
times distribution linking them can be determined from Equations 4-11 and 4-12.
Whenever possible two gauge stations with no major tributaries between them are
selected. The method can be applied when tributaries are present though the
values of the parameters thus obtained are less reliable.
The loss coefficient is used to achieve a mass balance. When time series
from two gauging stations are used, the value of the loss coefficient is not a
82
meaningful quantity. The presence of side flow from adjacent sources and
recharge due to springflow in natural systems ensures that losses incurred along
the flow path cannot be attributed to the usual mechanisms of evaporation and
infiltration. When a soil water balance model is used in conjunction with a routing
model, the loss coefficient can be viewed as a calibration factor whose value can
be determined analytically from mass balance considerations.
The main advantage of using the method of moments is that it can only
yield one correct solution for a given pair of time series. Optimization routines
may yield alternative solutions for the same set of data. The value of the resulting
parameters may also depend on the type of objective function used to arrive at a
solution. In addition, the convolution integral often results in highly non-linear
relationships for which conventional non-linear optimization routines such as the
Generalized Reduced Gradient (GRG2) algorithm (Lasdon and Waren, 1978)
used in the Excel Solver are unsuitable.
The accuracy of parameters determined by this approach is heavily
dependent on the temporal resolution of the flow data. There is currently no
public domain database of daily flow time series available for global applications.
Monthly flows produced by V?r?smarty et al (1996) are the most complete global
time series of runoff data available for this study. The effect of inaccuracies
arising from using monthly time series may be minimized by using a long time
series and by obtaining parameters from as many different pairs of stations as are
available for each river basin.
83
The key conceptual limitation of this method lies in the assumption that
peak discharge at the downstream gauging station is associated with the arrival of
flow from the upstream gauging station. The assumption may not hold in basins
where there is very high rainfall in the drainage area between the two stations.
Peak discharge at the downstream discharge station could occur independent of
the arrival of flow from the upstream station. Wherever possible two stations with
similar quantities of flow are selected. This ensures that the contribution of any
additional flow entering the reach between the two stations can only augment or
shift slightly the time of occurrence of peak flow at the downstream station. This
method of calibrating flow parameters is especially useful because it lends itself
to large-scale implementation on a river network.
In conclusion, velocity dispersion coefficient and loss coefficients for the
source to sink model must be derived from flow data. Methods for deriving these
values have been suggested in this section. The watershed based approach of
HMS is similarly discussed in section 4.2 below.
4.2 THE HYDROLOGIC MODELING SYSTEM APPROACH
The modeling approach used by HMS differs from the source-to-sink
approach in that the river basin is viewed as a series of hydrologic elements linked
to each other in series. HMS performs routing by passing flow from one element
to the next until it gets to a terminal point. Any local flow generated within each
element is added to flow entering it from upstream elements, and the total flow is
passed on to the next downstream element. The assumptions underlying an HMS
84
model, its component modules and methods as well as a calibration procedure are
described in this section.
4.2.1 HMS Model Assumptions
The assumptions of linearity and time invariance discussed in section 4.1.1
are also made in an HMS model using the Soil Conservation Service method for
watershed transformation and the Muskingum method for in-stream routing.
These assumptions allow an HMS model to deal with temporal variations within
each hydrologic element before moving to the next downstream element. Flows
from all upstream elements can subsequently be superimposed in the downstream
element. For a non-linear system, the response in each element is a function of the
flow.
Another important assumption is made about the sequential nature of the
connection between hydrologic elements. Subbasins are assumed to be connected
to river reaches via junctions. (The term subbasin is used instead of watershed in
HMS and is maintained in discussions about HMS in this dissertation). The
sequential connectivity assumption implies that there is no side flow into a routing
reach. In reality, river reaches as a whole act as drainage outlets for the subbasin
associated with them. Water from the subbasin can enter the river network at any
point within the subbasin. The definition of a subbasin connected to the river
reach at the subbasin outlet is therefore a simplification of the real life system.
The river basin cannot be divided into separate overland and in-stream flow
components without allowing for side flow in the river reaches. However, it is a
85
simplification that allows the river basin to be modeled in a vector (as opposed to
a raster) environment. The definition of a subbasin in which flow from each
location enters the river network at a particular point in along the reach requires a
fully distributed channel routing model. Such a model is computationally
intensive and difficult to implement in a vector environment. The partially
distributed approach of HMS allows parameters within hydrologic elements to be
lumped while defining different responses and parameters for each subbasin and
river reach.
4.2.2 HMS Model Components
Figure 4-2: Diagram of HMS Model Components
As shown in Figure 4-2, a typical HMS model is made up of three
components; a basin model, a precipitation model and a control specification
model. The basin model contains information describing the hydrologic elements
86
in the basin being modeled. It contains a list of the routing parameters and
element connectivity as well as the methods for computing losses, transforming
runoff and routing flow in each element. The classes of hydrologic elements
include subbasins, river reaches, junctions, diversions, reservoirs, sources and
sinks. Figure 4-3 shows the symbols used to represent these hydrologic elements
in HMS.
Figure 4-3: Elements of an HMS Basin Model
CRWR-Prepro, a preprocessor developed at the Center for Research in
Water Resources of the University of Texas at Austin, allows the element
connectivity to be derived from Digital Elevation Models (DEMs). This GIS
based preprocessing results in the creation of a basin file which can be interpreted
by HMS. The preprocessing procedure is described in section 3.4 of this
dissertation.
A portion of the precipitation falling on a subbasin is lost through
infiltration and evaporation. Methods of computing these losses in HMS include
the Soil Conservation Service curve number method, an initial plus constant loss
approach, and the Green and Ampt infiltration method. After accounting for these
losses, the remaining excess precipitation is transformed to the subbasin outlet
where it enters the stream system. HMS allows the user to specify one of several
87
runoff transformations for each subbasin. These include the Clark unit
hydrograph, Snyder?s synthetic unit hydrograph and the Soil Conservation
Service transformation method. When the runoff enters a river reach, it is routed
to the next downstream element using a user defined routing method. The routing
options available in HMS include the Muskingum method, the Modified Puls, the
Kinematic Wave and the Muskingum-Cunge method. Junctions in HMS serve as
control points where flow from multiple reaches or subbasins are aggregated
before being passed on to the next downstream element. Sinks on the other hand
serve as terminal points where flow is aggregated without being released. The
presence of control structures in real river systems is allowed for by the inclusion
of diversion points and reservoirs in the HMS basin description. Storage
dependent diversions and releases can be specified for the respective diversion
and reservoir components. The information contained in the basin model does not
vary with the quantity of flow.
The precipitation model on the other hand contains all information
describing time varying inputs. In this model, the quantity of precipitation
occurring during each time interval is specified. There are several options for
providing precipitation data. These include precipitation gages with or without
gage weights, gridded precipitation or spatially averaged precipitation with a
single value provided for each subbasin. The model also contains standard
precipitation events typical of certain geographical areas. However, these data are
only available for the United States, east of the Rocky Mountains. Discharge time
series for flow gauging stations are also stored in the precipitation module.
88
The control specification model is used to specify simulation start and end
times, and routing intervals. The control model plays an important role in HMS
modeling because in addition to defining simulation parameters it serves to tie the
routing process to reality. Actual dates and times are used in reporting simulation
results allowing for easy interpretation of flow data. The reference framework for
these dates is set in the control model. Routing intervals that can be specified in
an HMS model include 1, 2, 3, 4, 5, 6, 10, 15, 20 and 30 minutes as well as 1, 2,
3, 4, 6, 8, 12 and 24 hours. No other routing interval can be specified in HMS
though flow routed at one of these intervals can subsequently be resampled to
another routing interval.
For each HMS simulation, the user is allowed to specify one of each of the
component basin, precipitation and control modules. The routing methods used
for the simulations in this study are described in section 4.2.3 below.
4.2.3 HMS Routing Methodology
The routing methods defined for HMS subbasins and river reaches in this
study are described in this section along with the justification for their selection.
The range of methods available for this study was limited not by the methods
available in HMS but rather by the methods available in the preprocessor used in
the study. The benefits of using a method not available in the preprocessor would
be easily outweighed by the tedious nature of the task of piecing together and
parameterizing a hydrologic network of a large basin by hand.
89
The HMS preprocessor, CRWR-Prepro, currently allows the user to select
one of two methods for determining subbasin losses, namely the Initial plus
Constant loss approach, and the Soil Conservation Service curve number
approach. The Initial plus Constant loss approach was selected for this study
because it relies on currently available data while the curve number approach does
not. There is currently no global database of curve numbers, and the development
of such a database would be beyond the scope of this research. The initial and
constant losses can also be set to zero when input runoff is supplied from a soil
water balance model, allowing runoff to be transported to the subbasin outlet
without losses.
Runoff generated in a subbasin is limped together and routed to the basin
outlet using the Soil Conservation Service unit hydrograph method. In this
method, the distribution of flow at the subbasin outlet is computed by stretching
the SCS curvilinear hydrograph to fit the dimensions of the input. The parameters
of the curvilinear hydrograph are stored as a table of values. However, the SCS
triangular hydrograph, shown in Figure 4-4 (reproduced from Chow et al, 1988,
pp.229), is often used to illustrate its shape. This is because the triangular
hydrograph results in the same distribution of area under the rising and falling
limbs as the curvilinear hydrograph. The lag time in each subbasin is computed as
the larger of 0.6 times the length of the longest flow path divided by the flow path
velocity and 3.5 times the routing interval. The lag time thus computed is equated
to the time to peak, t
p
, in the figure while the routing interval is equated to the
duration of excess rainfall, t
r
. The other parameters, time to peak, t
p
, peak
90
discharge, q
p
, and duration direct runoff duration, t
b
, can be determined from the
other two parameters.
Figure 4-4: Soil Conservation Service Triangular Hydrograph
The preprocessor allows the user to select either the Muskingum or the
Lag routing method for in-stream routing. In the Lag method, flow is transported
through the river reach by pure translation using a lag time computed by dividing
the length of the reach by flow velocity. Muskingum method was selected
because of its ability to include a dispersion process in its routing of flow. In
HMS, each river reach is subdivided into multiple sub-reaches, and the
Muskingum method is applied in sequentially to each sub-reach. The subdivision
of a routing reach into sub-reaches ensures that numerical stability considerations
are satisfied. If the routing reaches become too long, the magnitude of the error
due the approximate numerical integration used in the Muskingum method
91
becomes significant compared to the magnitude of the flows. On the other hand,
if the reaches are too short the routing process tends towards pure translation, and
the effect of flow redistribution due to the Muskingum method becomes
insignificant. The number of routing sub-reaches, n, defined for each river reach
in CRWR-Prepro is given by Equation 4-17 (Olivera and Maidment, 1999).
Equation 4-17 1
/
2int +
?
?
?
?
?
?
?
=
t
VsLs
xn
where x is the Muskingum x coefficient
Ls is the length of the reach
Vs is the flow velocity in the reach
t? is the routing interval
The discharge from each characteristic reach is computed using the
Muskingum storage function and the continuity equation. The method of
coefficients is employed by HMS for this routing as shown in Equation 4-18 to 4-
21 below (Maidment, 1993, pp.10.7).
Equation 4-18
j
Q
C
3
j
I
C
2
1j
I
C
1
1j
Q
++
+
=
+
Equation 4-19
txK
Kxt
c
?+?
??
=
)1(2
2
1
Equation 4-20
txK
Kxt
c
?+?
+?
=
)1(2
2
2
92
Equation 4-21
txK
txK
c
?+?
???
=
)1(2
)1(2
3
where I is the inflow into the reach
Q is the discharge from the reach
K is the Muskingum K often interpreted as lag time
x is the Muskingum x coefficient which controls dispersion
t? is the routing interval
Dooge (1973) points out that solving the Muskingum equation using the
method of coefficients is equivalent to applying the response function given in
Equation 4-22 below.
Equation 4-22 )(
1)1(
exp
)1(
1
)(
2
t
x
x
xK
t
xK
tu ?
?
?
?
?
?
?
?
?
?
?
?
=
where )(t? is the Kronecker delta function which assumes a value of 1 at t
equal to zero and a value of zero elsewhere
The parameters that must be supplied by the user include the routing
interval, river reach and watershed lengths, flow velocities, and Muskingum x for
river reaches. The choice of routing interval takes into account the temporal
resolution of the input runoff and the desired resolution of the output flows. The
implication of the choice of interval is discussed in section 5.3.4. The lengths of
the routing reaches and watershed are automatically computed from flow length
grids by the preprocessor, CRWR-Prepro. Muskingum K is equated to travel time
and is computed by dividing the river reach length by the flow velocity supplied
93
by the user. The derivation of the user defined flow velocity and Muskingum x
parameters is discussed in section 4.2.4.
In summary, the choice of routing methods in this study was dictated by
the availability of data and by the limitations of the HMS preprocessor used. The
methods selected required two DEM based parameters (river and watershed flow
length) and two additional parameters (flow velocity and Muskingum x). Flow
losses were not taken into account in this study.
4.2.4 Estimating HMS Routing Parameters
Dooge (1973, pp.254) proposes a procedure for determining the
Muskingum x parameter based on channel characteristics for uniform channels.
The method requires the specification of channel slope, Froude number and a
reference discharge which allows the hydraulic radius to be calculated. The
approach cannot be used in this study since models requiring an assumption about
the quantity of flow were similarly rejected for the STS method in section 4.1.4.
The method of moments, applied in the source to sink modeling application, could
not be adopted for HMS because of difficulties in determining a single response
function for the multiple subreach Muskingum routing method applied in HMS.
The method of moments proposed by Dooge (1973, pp.254) corresponds to the
single reach case and is not applicable in HMS.
In this study, the Muskingum equation was employed to route the flow
observed at an upstream station to a downstream gauging station while using
HMS?s optimization routine to calibrate the parameters in the reaches in-between.
94
The optimization routine adjusts the routing parameters until a match is found
between the routed and observed flows at the downstream station. A spreadsheet
solution similar to this optimization computation is difficult to implement because
of the multiple routing subreaches involved. The Muskingum equation impulse
response function Equation 4-19 above is for the single routing reach case. From
the nature of the equation, it is easy to see that the equation is not linearly
decomposable. The presence of the conditional second term ensures that
convolving the function with itself will result in at least three terms, and hence
multiple convolutions will result in an equation that is not similar in form to the
original equation. The impulse response function for the multiple reach routing
case is not available in the literature. Consequently, the multiple reach case would
have to be solved by manually reconvolving the output of preceding subreaches.
An optimization run with multiple reaches quickly becomes too complex to
implement in a spreadsheet.
HMS has an inbuilt calibration program in which the values of specified
parameters can be calibrated from observed rainfall and corresponding stream
flow measurements. In this program, the user can specify the initial values of the
parameters, whether a parameter is to be changed or held constant during a
particular run, constraints on the values of parameters, the nature of the objective
function and the search method to be employed during an optimization run. In
addition, the program has inbuilt constraints (listed in Appendix 4 of this
dissertation) which limit the range of values a parameter can assume. The user
specified parameter constraints are ignored if they fall outside those imposed by
95
HMS. The method of parameter calibration proposed in the HMS documentation
(Peters, 1998) combines the rainfall-runoff transformation process with the
routing process in determining the routing parameters. It is therefore not ideal for
independent calibration of the routing process.
For independent calibration of in-stream routing parameters, the
optimization routine in HMS could be used to calibrate flow between two stations
in a manner similar to that described for the Source to Sink model in section 4.1.4.
The procedure adopted in this study was to create a separate model using the river
reaches located between the two stations. A source with an associated flow gauge
was inserted at the upstream station while a sink with an associated flow gauge
was inserted at the downstream end. Observed flows at the respective upstream
and downstream stations are introduced into the model by editing the gauge data
from the edit menu in the project definition window of HMS. An optimization is
then run to calibrate first the Muskingum K parameter and then the Muskingum x.
96
CHAPTER 5
Model Applications
5.1 THE APPLICATION BASINS
As outlined in section 1.2 of this dissertation, the testing of the source to
sink and watershed based routing models under various modeling scenarios is an
important aspect of this study. Two drainage basins were selected for this testing
based on the three criteria described below. The first criterion was that the basin
had to be large enough to justify the application of a continental scale routing
model. Effects, such as the spatially variability of parameters, are relatively
insignificant in small basins, but take on much more significance in large basins
because of their capacity to influence flow. The selection of large basins ensures
that such effects can be taken into account. Secondly, the application basins had
to have a variety of topographic and climatic regions. Such variety is necessary if
the models are to be applied globally. The availability of runoff data for model
calibration was the third criterion for selection. The Congo and Nile basins were
selected because they fulfilled these requirements and because of their location in
Africa, a continent of special interest to the author. The basins are described in
sections 5.1.1 and 5.1.2. This is followed in sections 5.2 and 5.3 by a discussion
of the results of model runs for the respective STS and HMS models. In section
5.4, STS and HMS routed flows are compared with observed flows in the Nile
basin.
97
5.1.1 The Congo Basin
Along with the Nile, the Niger and the Zambesi, the Congo is one of the
four major basins of Africa as shown in Figure 5-1.
Figure 5-1: The Four Major Basins of Africa
The Congo basin is located on the Equator in Central Africa and has a drainage
basin which extends through 7 countries including the Democratic Republic of
Congo, the Republic of Congo, the Central African Republic, Cameroon,
Tanzania, Angola and Zambia. From its sources in the East African lakes of
Tanganyika, Mweru and Bangweulu, it flows predominantly north for the first
half of its course. Near the Congolese city of Kisangani, the river makes a 90-
98
degree turn and heads westward towards its mouth at the Atlantic Ocean. Figure
5-2 shows a map of the Congo with its main tributaries, Sangha, Kasai, Ubangi
and Lualaba highlighted.
Figure 5-2: Area Map of the Congo Basin
The mean annual rainfall in the basin is about 1600 mm, and the mean
annual runoff is about 370 mm. With a basin surface area of over 3.8 million
square kilometers, this results in a mean annual discharge of approximately 45000
m
3
/s at its mouth in the Atlantic Ocean. The Congo ranks second only to the
99
Amazon in the quantity of discharge at the mouth. The high rainfall in the region
allows flow in the basin to be simulated without having to deal with the problems
of runoff generation in arid basins which are best dealt with in soil water balance
models.
Its location in the equatorial belt ensures that different parts of the Congo
basin receive substantial rainfall throughout the year. The northern and central
portions of the basin have two major rainfall seasons which begin in March and
October each year (Kazadi, 1996). In the south, the two rainfall seasons gradually
merge into a single season beginning in December and lasting for six months each
year. These rainfall peaks are associated with the passage of the Inter-Tropical
Convergence Zone (ITCZ). This is a large zone of low pressure caused by
excessive heating from an overhead sun. During the northern summer, the midday
sun is directly overhead in the tropical regions of the Northern Hemisphere. This
results in higher temperatures and consequently, lower air pressures at the surface.
Moist air flows from the oceans towards these low pressure areas. The moisture is
released as rainfall on the land surface when the air is forced to rise on entering
the convergence zone or by orographic effects. The ITCZ causes heavy rainfall in
the areas it passes over as it moves north and south between the tropics during the
respect northern and southern summers. The Congo basin is thus representative of
a large river basin in which the spatial distribution of input varies significantly
with time.
100
5.1.2 The Nile Basin
Located in the northeastern part of Africa, the Nile basin is another of
Africa?s large river basins. Figure 5-3 shows the Nile basin with it major
tributaries and gage locations highlighted.
Figure 5-3: Major Tributaries and Flow Gages of the Nile Basin
101
Its two longest tributaries have their sources in Lake Albert and the
Luvironza River in the East African Highlands. The Albert Nile and Victoria Nile
as the two tributaries are called emerge together from Lake Albert to form the
White Nile which subsequently flows northward through northern Uganda and the
Sudan. Near the Sudanese capital of Khartoum, the White Nile is joined by
another major tributary, the Blue Nile, which rises out of Lake Tana in the
Ethiopian Highlands. The merged Nile then flows northward through Egypt to its
mouth at the Mediterranean Sea. From its most distant source to its mouth, the
Nile flows a distance of 6680 kilometers, making it the world?s longest river. Its
drainage basin extends through 8 countries and covers an area of 3.25 million km
(Revenga et al, 1998).
Rainfall in the basin varies from 1200 mm per year in the Highland Lake
region to almost zero in the desert region near its mouth (Hurst and Phillips,
1938). This rainfall distribution provides a unique setting in which to study the
behavior of water as it flows over long distances without significant input from
the watersheds closest to its mouth. In addition, there is a wide variety of
topography which influences the flow of water in the Nile basin. The Blue Nile
basin is characterized by steep slopes on the order of 1m per km between Lake
Tana and Khartoum. In the Sudd marshes in southern Sudan, the White Nile flows
through a very flat region with slopes on the order of 0.01 m per km or less (Hurst
and Phillips, 1938). Flow velocities in these marshes are an order of magnitude
lower than in the rest of the basin, and the river loses about half of its flow on its
way through the marshes. The East African Highlands are also characterized by
102
lakes and marshes though the influence of these on stream flow is not as
pronounced as those of the Sudd. The Nile is thus representative of a large river
basin with several zones of significantly different topography and hydrology.
5.2 SOURCE TO SINK MODEL RUNS
In section 3.3, the development of a global STS model was described.
Drainage basins were delineated from stretches of coastline in that description. In
the watershed based modeling approaches of HMS, drainage basins are delineated
from a single 1 km cell. The boundaries of drainage basins in the two approaches
are thus inconsistent and their results cannot be compared. To allow for
intercomparison of model results, the drainage basin boundaries used in the STS
and HMS models must coincide exactly. The subbasins from the HMS delineation
were merged together to form a single drainage basin. Sources for the STS model
were defined by intersecting the merged basin with a 0.5 by 0.5 degree mesh.
Average travel distance from each of source to the basin outlet was determined as
the mean of all flow length grid cells underlying each source. The resulting STS
model consists of 1378 sources with flow parameters determined from the same
DEM as the HMS model. Figure 5-4 below shows the sources defined in the
Congo Basin model used for the model intercomparison runs. The darker areas
represent areas of increasing flow length.
103
Figure 5-4: Source to Sink Model of the Congo Basin
5.2.1 Longitudinal Decomposability of the STS Model
In developing models of hydrologic systems, methods and conceptual
models which yield the same value irrespective of the size of the modeling unit
selected are more reliable. The control volume in a source to sink model is the
flow path. Hence it is important to affirm that a source will yield the same
104
response at its sink irrespective of the number of subreaches into which its path is
divided. Olivera (1996) identifies the property of flow paths or other control
volumes to yield the same response irrespective of the number of constituent
segments as longitudinal decomposability. Given the longitudinal
decomposability of the first passage times distribution (Olivera, 1996), a source to
sink model using this equation is expected to be longitudinally decomposable
along the flow path from the source to the sink. A set of runs conducted to verify
this property is described in this section.
Figure 5-5: Setup for Evaluating Longitudinal Decomposability
The simplest approach to testing the longitudinal decomposability of a
flow path is to divide it into a number of subreaches and determine if the resulting
flow path response to an impulse is the same as that of the single reach. In a river
basin, flows routed directly from the point of runoff generation to a sink could be
compared with flows routed sequentially through intermediate flow elements. A
source to sink model was setup with a set of four rectangular cells, each with a
surface area of 1000 km
2
. As shown in Figure 5-5, the distance between the
centroids of adjacent cells varied from 2000 km between cell 1 and 2, to 1000 km
105
and 1200 km between cells 2 and 3 and cells 3 and 4, respectively. From cell 4,
the most downstream, to the outlet, sink 5, a further distance of 800 kilometers
was allowed such that cells 1, 2 and 3 are 5000, 3000 and 2000 kilometers away
from sink 5, respectively. An instantaneous impulse of 1 millimeter of water was
applied simultaneously to all four cells and routed to the sink by two different
approaches. In the first approach, source to sink routing, flow is routed directly
from the respective cells to the sink. In the second approach, cell to cell routing,
the input is routed from each cell to the next downstream cell until it arrives at the
sink. The diffusion wave equation is used to route flow in both cases. A flow
velocity of 0.3 m/s and a dispersion coefficient of 1500 m
2
/s are assumed
throughout the setup. Figure 5-6 below illustrates the difference between the two
routing approaches.
Figure 5-6: Source to Sink versus Cell to Cell Routing
106
The discharge arriving at the sink from the approaches are shown in the
Figure 5-7 below. The discharge hydrographs indicate that for the setup used in
this run, the difference between the runs is negligible. The cell to cell model
results in a slightly longer tail but the magnitude of this difference is very small
compared to the long distance over which flow was routed.
Figure 5-7: Comparison of Diffusion Type STS and CTC Discharge
The first moment of the source to sink model response is approximately
0.36 hours higher than that of the cell to cell model. A comparison of the second
moments of the respective responses indicates that the cell to cell model does in
fact exhibit more dispersion with a 0.24% increase in second moment. The
107
magnitude of this difference is very small and could be ignored given
uncertainties in the routing parameters.
5.2.2 Effect of Spatial Resolution of the STS Model
Another set of runs were setup to study the effect of spatial resolution on
model outputs. In this set of runs, the effect of additional modeling detail
achieved through the definition of additional sources on the basin response is
examined. A suitable modeling spatial resolution can be defined, for a given
modeling time step, by determining the spatial resolution below which no
significant change occurs in the basin response. In other words, the task is to
determine the coarsest spatial resolution at which a basin can be represented with
a STS model without significantly distorting the basin response. If the modeling
resolution used is finer than required, the model will generate the same result but
at a higher cost of computer resources and processing time. A coarser resolution
model on the other hand would begin to exhibit some of the properties of scale
dependence present in other models.
The Congo basin was modeled with sources defined at different
resolutions. The resolutions used include half a degree, 10 arc-minutes and 5 arc-
minutes. At each of these resolutions, the flow length from each source to its
corresponding sink was calculated as the mean value of the cells in the underlying
30-arcsecond flow length grid. The total drainage area as well as velocity and
dispersion coefficient were kept constant for each of the runs. Figure 5-8 shows
the sources defined in the Congo basin at each of the three resolutions.
108
Figure 5-8: Congo Basin STS Models at 30, 10 and 5 arcminute Resolutions
109
In response to this question, the STS model was run with the terrain being
subdivided into different resolutions. The resolutions used include 30, 10 and 5
minutes. The results of this simulation (shown in Figures 5-9 (a) and (b))
indicated that much of the noise seen in the STS hydrograph is due to the
coarseness of the spatial discretization. A relatively smooth curve can be obtained
by using a finer discretization of the drainage area. The improved smoothness is
however to be differentiated from the change in dispersion characteristics
observed with cell to cell models. The definition of finer resolution sources results
in less deviation of a given flow value from preceding and subsequent values. It
does not however result in changes in the degree of spreading out experienced by
inflow as would be the case if the dispersion characteristics of the model were
altered. The overall effect of going from a 30 arc-minute to a 10 arc-minute
resolution (Figure 5-9 (a)) is much more pronounced than that resulting from the
change from a 10 arc-minute to a 5 arc-minute resolution (Figure 5-9 (b)). This
implies that as the resolution of sources becomes finer, the basin response
becomes better defined though its dispersion characteristics remain unchanged.
110
Figure 5-9(a): Comparison of Congo Basin Responses for 30 and 10 arcminute
Resolution STS Model with Pure Lag Routing
Figure 5-9(b): Comparison of Congo Basin Responses for 10 and 5 arcminute
Resolution STS Model with Pure Lag Routing
111
When the dispersion effect is included in the routing process, the effect of
resolution is much less pronounced. While a slight deviation may be observed at
the points of inflexion, the basin responses shown in Figure 5-10 indicated that
the overall basin is relatively insensitive to the resolution of the sources. The
spreading out of inflow along each flow path is sufficient to obsure the effect of
spatial discretization.
Figure 5-10: Effect of Spatial Resolution on Dispersed STS Hydrograph
In Table 5-1, a comparison of the lag times and variance of the pure lag
and dispersed basin response for the three resolutions is presented. For each set
of runs (pure lag and dispersed), the 5 arc-minute resolution is used as the
112
baseline for determining magnitude of change in lag time and standard deviation.
The results indicated that for both sets of runs, the magnitude of changes in
dispersion characteristics with changes in resolution were small as indicated by
the high correlation values. The higher correlation for the dispersed routing case
indicates that the inclusion of the dispersion effect makes the basin response much
less sensitive to the resolution of the sources. Comparing the first set of runs to
the second also reveals that the introduction of the dispersion effect results in a
reduction of the overall basin lag time by a full routing interval from 90.2 days to
89.2 days. There is a corresponding increase in standard deviation though of a
much smaller magnitude.
Table 5-1: Effect of Spatial Resolution on Lagged and Dispersed Basin Response
Spatial
Resolution
First
Moment
(days)
Second
Moment
(days
2
)
Correlation Change in
Lag Time
(hrs)
Change in Std.
Deviation (hrs)
Lagged Basin Response
30? 90.223 1521.525 0.9533 0.480 3.549
10? 90.244 1528.636 0.9975 0.024 1.364
5? 90.243 1533.083 BASELINE BASELINE BASELINE
Dispersed Basin Response
30? 89.135 1557.380 0.9997 1.272 4.379
10? 89.189 1567.439 0.9999 0.024 1.325
5? 89.188 1571.815 BASELINE BASELINE BASELINE
113
5.2.3 Effect of Spatial Distribution of STS Parameters
In the previous set of runs in section 5.2.2, the effects of varying the
spatial resolution of STS model sources were explored. This section describes a
series of runs in the Nile basin to examine the effects of spatially varying routing
parameters on the basin response. The goal in this section is to determine whether
or not the assumption of uniform routing parameters made in several source to
sink models (Naden, 1993, Lohmann et al, 1996) results in a significant alteration
of the impulse response of a large drainage basin. The routing parameters
examined in this study are the velocity, dispersion coefficient and loss coefficient.
Given a field of spatially distributed parameters, a process of determining
equivalent uniformly distributed parameters is required. The process adopted
depends on whether the spatially distributed parameter is scalar or vector in
nature. Scalar quantities can be added together and subtracted from each other
because they are linearly decomposable. If the quantities are vector in nature then
both the magnitude and direction must be taken into account. Vector quantities
are therefore not linearly decomposable. For example, the mean value of velocity
in a field of spatially varying velocities is a meaningless quantity since each
individual velocity value in the field has a magnitude and a direction. The
arithmetic mean value of a field of parameters is hence only meaningful if the
parameters in the field are scalar quantities. To determine a representative
uniform value for a field of vector quantities, two approaches may be adopted.
One approach would be to first resolve all the values in the field into a single
114
direction and determine the magnitude of the parameter in the resolved direction.
The second approach would be to convert the vector parameter into an associated
scalar parameter, determine a single value of the scalar quantity by summation or
averaging and use the single scalar value to determine a representative value of
the vector parameter for the field. For example, a velocity field could be
converted into a field of flow time, a scalar quantity. A mean or total value of
flow time could be determined from the resulting field and used to derive a single
value of velocity for the entire field. By so doing, the total lag time in the basin,
the first moment of the basin response, is preserved allowing the effect of spatial
variability to be analysed independently.
The spatial variability runs were conducted in the Nile basin because of
the availability of flow data throughout the basin (Hurst and Phillips, 1938).
Figure 5-11 shows the STS model of the Nile basin with key parameteric zones
highlighted. As described in section 5.1, the Blue Nile basin is characterized by
high flow velocities and low dispersion while the Sudd marshes are characterized
by low velocities and high dispersion. Part (b) of the figure shows the sources
defined in the basin with longer flow lengths represented by darker shades. The
model consists of 1808 sources attached to a single sink.
115
Figure 5-11: Nile Basin (a) Parametric Zones and (b) Sources
The STS model was run in several iterations. In the first set of runs, the
effect of hydrodynamic dispersion was held constant by applying a dispersion
coefficient of 1500 m
2
/s for all runs. The model was set up first with a uniform
flow velocity and then with spatially variable velocity. For the variable velocity
case, three velocity zones were defined. In the Blue Nile basin, water flows from
an elevation of 3000 m above sea level in the Ethiopian Highlands to an elevation
of 1400 m at its confluence with the White Nile at Khartoum, a distance of about
1700 km. Observed flows published by the Nile Control Department in Egypt
(Hurst and Phillips, 1938) indicate that peak discharges take aproximately 4 days
to travel the 418 kilometers from Wad el Aies to Khartoum. This indicates a flow
116
velocity of approximately 1.2 m/s. In the Sudd Marshes in the Sudan, the White
Nile takes about three months to travel the 711 kilometers between Mongalia and
Malakal resulting in a velocity of about 0.1 m/s. Similarly, velocity in the White
Nile outside of the marshes is estimated based on a flow length of 794 km and a
travel time of 20 days between Malakal and Mogren just south of Khartoum as
0.5m/s.A1km
2
grid of these three velocities was first developed, and a grid
containing its inverse was calculated. Flow time from each grid cell to the basin
outlet was computed by using Arc Info?s weighted flow length function with the
flow direction grid and the grid containing the inverse of velocity as inputs. The
resulting grid contains cumulative flow time along the flow path from each grid
cell to the basin outlet. The flow time for each source was then determined as the
mean of all 1 km
2
cells in the flow time grid. The flow time for each source can
be used for computing the discharge hydrograph directly from Equation 4-11 as
described in section 4.1.3 of this dissertation. It can also be used to compute a
new flow path velocity which preserves total flow time along each flow path from
source to sink. The latter option was adopted in this study because it allows the
diffusion wave equation to be used in its original form without having to
introduce additional parameters. The basin response for this run was compared
with the response for the equivalent uniform velocity case. The equivalent
uniform velocity was calculated by dividing the sum of flow time along all flow
paths in the basin by the sum of flow length for the whole basin. The equivalent
uniform velocity therefore preserves the average basin lag time, allowing the two
simulations to be compared. The basin response for the respective non-uniform
117
and uniform velocity cases are illustrated in Figure 5-12. The light grey line
represents the basin response for the non-uniform velocity case while the darker
line represents the uniform velocity case.
Figure 5-12: Effect of Spatial Distribution of Velocity with Uniform Dispersion
Coefficient
For the Nile basin, the difference between using a uniform velocity as
opposed to a non-uniform velocity is significant. The basin response is modal for
uniform velocity and bimodal for non-uniform velocity. A routing method which
could not account for the spatial variability of flow velocity could not possibly
represent the hydrology of a basin like the Nile. The run serves to highlight the
importance of accounting for spatially variable velocities in large basins.
The next set of runs was conducted to study the effect of spatially
distributed dispersion coefficient on the basin response. In this set of runs, the
118
flow velocity was kept constant while uniform and non-uniform dispersion
coefficients dispersion coefficients were adopted for respective runs. For the non-
uniform case, dispersion coefficients of 4500 and 300 m
2
/s were assigned for the
Sudd and Blue Nile basins respectively while the rest of the basin was assigned a
coefficient of 1500 m
2
/s. These dispersion coefficients were order of magnitude
estimates based on the measured values reported in Fischer (1979) and Seo and
Cheong (1998). They were not intended to duplicate the actual basin response but
rather as estimates to allow the effect of spatial variability to be studied. Table 5-2
summaries the flow routing parameters used in the Nile basin.
Table 5-2: Flow Routing Parameters used in the Nile Basin
Region Velocity (m/s) Dispersion Coefficient (m
2
/s)
Blue Nile Basin 1.2 300
Sudd Marshes 0.1 4500
Rest of the Nile Basin 0.5 1500
A1km
2
grid containing the dispersion coefficient values was created for
the Nile Basin. From Equation 4-10 in section 4.1.3, variance along each flow
path can be estimated given flow path velocity and length as well as the local
values of dispersion coefficient for all intervening segments. The flow path
variance for the flow path from each source to its sink was computed by
averaging the flow path variance for all grid cells within each source. A
representative dispersion coefficient was then computed for the flow path
119
associated with each source using Equation 4-12. This coefficient results in the
same amount of variance along the flow path from source to sink as would result
from routing flow sequentially through the respective dispersion zones. As with
the flow path velocity, the dispersion coefficient was inserted into the diffusion
wave equation and used in the routing. The results from this run were compared
with those from the uniform dispersion coefficient case. A uniform dispersion
coefficient which resulted in the same total variance when applied uniformly to all
the flow paths was then determined from the sum of flow path variance in the
non-uniform case. For the coefficients zones defined in the Nile basin, the
equivalent uniform dispersion coefficient was determined as 907.55 m
2
/s.
Figure 5-13: Effect of Spatial Distribution of Dispersion Coefficient with Uniform
Velocity
120
The results of the two simulations are shown in Figure 5-13. The lighter
gray line represents the basin response with both dispersion coefficient and
velocity uniform while the darker line represents the non-uniform dispersion
coefficient case. Peaks and troughs are slightly less defined in the non-uniform
dispersion case, particularly in the latter part of the hydrograph. This seems to
imply that there is more complete dispersion of flow in the non-uniform flow
case. However, the difference between the observed values in the two runs is
minimal compared to the obvious disparity observed in similar runs with the
velocity. This would seem to imply that velocity is a much more significant
parameter in determining the response of a basin and that the influence of
spatially distributed dispersion coefficients could be ignored without significantly
altering the basin response. However, before such an assertion could be made, the
combined effect of spatially variable velocity and dispersion coefficient had to be
examined. A simulation run was consequently performed with both of the
parameters non-uniformly distributed.
Figure 5-14 shows a comparison of the non-uniform velocity, uniform
dispersion coefficient case with the case in which both parameters are non-
uniform. Non-uniform dispersion coefficient did not have a marked influence on
the basin response when applied with uniform velocity. However, when velocity
is non-uniform, the effect of a non-uniform dispersion coefficient on the Nile
basin response was significant. The second of the two discharge peaks observed
in the Nile was substantially attenuated, resulting in a 30% reduction in peak
flow. The reason for this difference is that the topographic conditions that result in
121
low velocities are the same as those associated with high dispersion coefficients.
The peak attenuation associated with low velocities is therefore compounded by
the high dispersion coefficients in the same region. Areas with higher velocities
are likewise more likely to have low dispersion coefficients. Water travels quickly
through such areas and the effect of dispersion is not very pronounced.
Figure 5-14: Effect of Spatial Distribution of Dispersion Coefficient with Non-
Uniform Velocity
The key inference from the runs with spatially variable parameters is that
if variable velocity zones exist in a basin then it is important to considered the
effect of non-uniform dispersion coefficient as well. On the hand, a uniform
dispersion coefficient can be assumed if the assumption of a uniform velocity can
be shown to reasonably represent flow conditions in the basin.
122
5.2.4 Effect of Temporal Resolution of STS Model
Having examined the effect of spatial resolution and spatial variability on
basin response, another series of runs was conducted to examine the effect of
temporal resolution of the routing model on the basin response. The goal in this
set of runs was to determine how the choice of a routing time interval impacts the
flow time and dispersion characteristic of a drainage basin. This determination is
important because it provides useful information on whether or not the source to
sink model can be used to establish a global database of velocity and dispersion
coefficients. A model that is temporally scale independent will yield the same
results for a given set of parameters irrespective of the routing interval. It can
hence be used as a standard for the determination of routing parameters for use in
other routing models. STS models of the Congo Basin were used to route inflow
at routing intervals of 6, 12 and 24 hours respectively using a uniform flow
velocity of 0.3m/s and a dispersion coefficient of 1500 m
2
/s. Figure 5-15 shows
the resulting time series plotted with a single horizontal time scale.
123
Figure 5-15: Effect of STS Routing Interval on Congo Basin Response
The graphs of the respective routing intervals are indistinguishable from
one another. A comparison of their respective first and second moments in Table
5-3 reveals that there is a minimal shift whose magnitude is much smaller than the
routing interval.
Table 5-3: Effect of STS Routing Interval on Congo Basin Response
Routing
Interval
First Moment
(days)
Second Moment
(days
2
)
Change in Lag
Time (hrs)
Change in Std.
Deviation (hrs)
6 hours 89.169 1549.488 0.432 0.431
12 hours 89.174 1549.930 0.312 0.297
24 hours 89.187 1550.903 BASELINE BASELINE
124
The results of this simulation indicate that the source to sink model is
relatively insensitive to choice of routing time interval. Routing parameters
determined at one temporal resolution would therefore be applicable to other
temporal resolutions.
5.3 HYDROLOGIC MODELING SYSTEM MODEL RUNS
Having examined the effect of various routing and model parameters on
the response of a source to sink model, a similar set of simulations were
performed on the watershed based model for comparison. As discussed in section
3.4, validation of the source to sink routing process in this study is performed by
comparison with an HMS model of the same basin. A series of simulations
performed with the watershed based approach of HMS are described in this
section. The change in the response of a sample HMS reach due to longitudinal
decomposition is examined in section 5.3.1. The Congo basin response is
compared for different models resolutions and routing intervals in sections 5.3.2
and 5.3.4 respectively while the ability of an HMS model to represent spatially
distributed parameters in the Nile Basin is examined in section 5.3.3.
125
Figure 5-16: Sample HMS Schematic for Congo Basin
Figure 5-16 shows an HMS schematic for the Congo Basin delineated at a
10000-cell threshold. The process of deriving an HMS basin file containing a
description of the hydrologic elements and the connectivity among these elements
is described in section 3.4 of this dissertation. In that delineation, the threshold for
stream delineation was specified as 10000 km
2
. A more detailed model of the
Congo basin was created for the model intercomparisons. The parameters for the
second delineation of the Congo basin are summarized in Table 5-4. In the Nile
basin, spatially distributed velocities, Muskingum x and loss parameters were
used as described in section 5.3.3.
126
Table 5-4: Parameters for the Congo Basin Delineation
Process Inputs Parameter Value Outputs
Stream
Definition
Flow
Accumulation
Grid
Stream Definition
Threshold (km2)
1000 Stream
Grid
Stream
Linking
Stream Grid Number of Links 1791 Stream
Link Grid
Outlet
Definition
Stream Link
Grid
Number of Outlets 1791 Outlet
Grid
Watershed
Definition
Flow Direction
Grid, Stream
Link Grid
Number of
Watersheds
1791 Watershed
Grid
Stream
Vectorization
Stream Grid Minimum Stream
Length
710 Stream
Coverage
Watershed
Vectorization
Watershed Grid Number of
Watersheds
1791 Watershed
Coverage
SCS Unit
Hydrograph
watershed velocity 0.3 m/s
initial losses 0
Calculating
Watershed
Attributes
Watershed
Coverage, Flow
Direction Grid
constant losses 0
Watershed
Parameters
Muskingum Method:
in-stream velocity 0.3 m/s
Calculating
Stream
Attributes
Stream
Coverage, Flow
Direction Grid
Muskingum x 0.3
Stream
Routing
Parameters
Subbasins 1791
River reaches 943
Junctions 894
Basin File
Creation
River coverage,
watershed
coverage
Sinks 1
HMS
Basin File
127
5.3.1 Longitudinal Decomposability of the HMS Model
Longitudinal decomposability along each flow path linking a source to its
sink in the diffusion type STS model was verified in section 5.2.1 of this
dissertation. It is imperative that a similar verification be conducted for the HMS
model. A model was set up in HMS consisting of a 162 kilometer river reach
linked to a subbasin at its upstream end and a sink at its downstream end as shown
in Figure 5-17. The watershed served as a source of input for the stream, and
input runoff was translated instantaneously through it and without redistribution
or losses.
Figure 5-17: HMS Schematic for Longitudinal Decomposability Runs
The Muskingum routing method was used to route flow through the river
reach. The routing parameters used include a flow velocity of 0.3 m/s, and a
Muskingum x of 0.3. This resulted in a Muskingum K of 150 hrs and a
Muskingum n of 4 for the reach. HMS routes flow by applying the Muskingum
equation successively to all n sub-reaches until the flow arrives at the end of the
river reach. In built HMS parameter constraints limit the number of sub-reaches
into which a river reach can be divided. Hence the task was to determine the
extent to which longitudinal decomposability is maintained in a reach when the
HMS imposed limitations on the number of subreaches, n, are adhered to. An
128
impulse of 100m
3
/s was applied to the subbasin, and the discharge at the sink was
observed for various values of n between 4 and 8. HMS returned warning
messages for values of n outside this range. These messages indicated that the
Muskingum equation was unstable for the routing parameters defined.
Muskingum n values of 4 and 8 correspond to sub-reach lag times of 18.75 hrs
and 37.5 hours, respectively. This implies that HMS allows the user to subdivide
the 150-hour lag time into four 37.5 hour segments, eight 18.75 hour segments or
other whole number of uniform segments between 4 and 8. The problem could be
complicated even further by subdividing the single reach into smaller reaches by
introducing junctions and then calculating a Muskingum n for each reach.
However, HMS generates error messages when the lag time along the routing
reach is less than the routing interval, 24 hours in this case. The largest number of
Muskingum routing reaches that could be created from the 150 hour reach is six.
A routing reach with a lag time of just over 24 hour would require a Muskingum n
value of 1. Thus the introduction of additional junctions would not lead to a
Muskingum n value outside the 4 to 8 range, assuming a 24-hour routing interval.
Figure 5-18 below shows the river reach response for different values of n with
the Muskingum K and x values unaltered.
129
Figure 5-18: Changes in River Reach Response with Muskingum n
The figure indicates that for the same value of Muskingum x and total
reach length, routing flow sequentially over several small distances produces less
dispersion than would result from routing the same flow over larger stretches of
distance. The Muskingum routing method is consequently not linearly
decomposable. The significance of this result is that it provides an indication of
the degree of variability permissible in HMS. Table 5-5 provides an analysis the
moments of the time series generated for each of the 5 values of Muskingum n
used in the simulations. The equivalent values of dispersion coefficient are
determined by equating the second moments of the Muskingum and diffusion
wave equations for a routing subreach. The dispersion coefficient can thus be
determined from Equation 5-1 (Montes, 1998, pp. 574).
130
Equation 5-1 )5.0( xLvD ?=
where D is the equivalent dispersion coefficient in m
2
/s
L is the length of the reaching in m
x is the dimensionless Muskingum x
v is the velocity of translation through the reach in m/s
Table 5-5: Effect of Muskingum n on River Reach Response
Number of
Routing
Subreaches
First
Moment
(days)
Second
Moment
(days
2
)
Equivalent
Dispersion
Coefficient (m
2
/s)
Change in
Standard
Deviation (hrs)
4 6.250 3.905 2429 BASELINE
5 6.250 3.125 1944 5
6 6.250 2.603 1619 8.7
7 6.250 2.233 1389 11.6
8 6.250 1.953 1215 13.9
Even within the HMS imposed limitations, a 100% increase in dispersion
coefficient can be achieved by simply altering the number of subreaches. The
result indicates that the response of a routing reach is very sensitive to spatial
resolution. If the Muskingum x parameter is allowed to assume different values,
there is a limited range of dispersion coefficients that the reach can assume given
the routing interval and velocity.
131
Table 5-6: Range of Dispersion Coefficients of a Muskingum Routing Reach with
24 hour Routing Interval and Velocity of 0.3m/s
Muskingum
x
Maximum
Routing Reach
Lag Time (hrs)
Minimum
Routing Reach
Lag Time (hrs)
Maximum
Dispersion
Coeff. (m
2
/s)
Minimum
Dispersion
Coeff. (m
2
/s)
0 150.0 12.2 24300 1976
0.1 120.0 13.4 15552 1737
0.2 60.0 15.0 5832 1458
0.3 40.0 17.2 2592 1389
0.4 30.0 20.0 972 648
0.5 UNSTABLE UNSTABLE N/A N/A
Table 5-6 lists the range of Muskingum K and the corresponding
dispersion coefficients that are available to an HMS model with a routing interval
of 24 hours and a velocity of 0.3 m/s. For each Muskingum x value, the maximum
and minimum values of Muskingum K were determined by repeatedly altering the
K value for a single reach with n =1. The respective maximum and minimum
Muskingum K values were determined as the values greater and less than which
error messages regarding the numerical stability of the Muskingum solution are
triggered. Flow velocity for the reach was then determined by dividing the lag
time by the actual length of the reach. The significance of this result is that it
demonstrates that by choosing a particular value of the Muskingum x, the modeler
limits the range of dispersion coefficients. However, the lag time for a particular
132
reach also plays a role in determining the amount of dispersion the reach
experiences. The application of a uniform Muskingum x coefficient does not
therefore result in a model with uniform dispersion properties. To achieve such a
model, each river reach would have to be assigned a different value of
Muskingum x. In section 5.3.3, an HMS model of the Nile basin is set up using
equivalent Muskingum K and x parameters calculated from a spatially distributed
field of velocity and dispersion coefficients.
5.3.2 Effect of Spatial Resolution of the HMS Model
To examine the effect of spatial resolution on the basin response, two
HMS models of the Congo were created using stream delineation thresholds of
1000 and 10000 km
2
, respectively. These thresholds resulted in HMS models
containing 1791 and 214 subbasins, and 943 and 152 river reaches respectively.
River reaches with lag times greater than the maximum permissible lag time of
150 hours in HMS were subdivided into shorter reaches using junctions. The
basin responses from the respective delineations are shown in Figure 5-19.
133
Figure 5-19: Effect of Delineation Threshold on Congo Basin Response for HMS
Model with Muskingum Routing
The specification of a higher cell threshold implies that input runoff
spends more time on the land surface as opposed to within the river channel. The
effect of the approximate unit hydrograph process on the overall basin response
thus becomes more pronounced. Conversely, the reduction in time spent in the
river channel, results in an overall reduction in the amount of dispersion the flow
undergoes for the same set of routing parameters. Consequently, the larger
delineation threshold, the less dispersed the resulting hydrographs. Also note that
since the Muskingum x is scale dependent, the value of 0.3 used for the model at
the finest resolution is no longer applicable at the coarser scale. A procedure for
134
predicting the appropriate value of Muskingum x to use for a change in resolution
isdescribedinsection5.3.3.
The HMS basin responses at the respective 1000 and 10000 cell
resolutions were also compared with STS basin responses using various values of
dispersion coefficient. Figure 5-20 shows a comparison of the finer resolution
HMS and STS model response of the Congo basin. For both models, a flow
velocity of 0.3m/s was assumed. A dispersion coefficient of 2000 m
2
/s was
assumed for the STS model while a Muskingum x of 0.3 was assumed for the
HMS model. The basin responses from the two models would map almost
perfectly on top of each other if the STS model response were shifted 2 days
ahead. Possible explanations for this shift are discussed in section 5.2.4.
Figure 5-20: Comparison of 1000 Cell Threshold HMS Model Response with the
STS Model Responses in the Congo Basin
135
The Congo basin response for the HMS model delineated at the 10000 cell
threshold was similarly compared with an STS model response with a dispersion
coefficient of 1500 m
2
/s. This value of dispersion coefficient was selected based
on the range of possible values of the coefficient for a Muskingum routing reach
with a velocity of 0.3 m/s, a Muskingum x of 0.3 and a routing interval of 24
hours. The results of this comparison are shown in Figure 5-21.
Figure 5-21: Comparison of 1000 Cell Threshold HMS Model Response with the
STS Model Responses in the Congo Basin
The match between the two basin responses is not as close as that obtained
with the 1000 cell threshold model. The likely explanation for this discrepancy is
the increased role of the subbasin response function in the overall basin response.
136
The SCS synthetic unit hydrograph used for this analysis plays a greater role in
the overall basin response at this threshold than it did at the 1000 cell threshold.
In summary, the comparison of responses of two HMS models of the same
basin delineated at different thresholds revealed the basin response is very
sensitive to spatial scale. When the scale of the model is changed, the basin
response changes even if the routing velocity and Muskingum x parameters
remain unchanged. Each HMS model must be parameterized separately making
the method unsuitable for establishing scale independent routing parameters. The
next logical step is therefore to determine whether routing parameters established
using a more scale independent model like the diffusion type source to sink model
can be used to parameterize an HMS model. An analysis of this nature is
describedinsection5.3.3.
5.3.3 Effect of Spatial Distribution of HMS Parameters
The manner in which hydrologic elements are defined in an HMS model
allows it to incorporate spatially distributed routing parameters. Hence the task in
this section is not to define a procedure for representing spatially distributed
parameters as was done for the source to sink model in section 5.2.3. Instead, a
procedure is defined for computing equivalent Muskingum K and x parameters
from a field of spatially distributed dispersion coefficients, D, and flow velocities,
v. The same dispersion coefficients and velocity fields used in the Nile basin STS
model in section 5.2.3 were provided as inputs, thus allowing the basin responses
from the respective model runs to be compared.
137
The first step in determining equivalent Muskingum routing parameters
was to extract a flow direction grid containing only in-stream cells by multiplying
the stream grid with the flow direction grid. The resulting grid was then
multiplied by the outlet grid to insert sinks, in the flow direction grid, at the end of
each river reach. As in section 4.1.3, weight grids of (1/v) and (D/v
3
) were applied
for the computation of first and second moment using the weighted flow length
function. Equivalent values of Muskingum x were computed for the longest flow
length in each stream reach by equating the second moments of the Muskingum
equation with those of the diffusion wave equation using Equation 5.2 which is
Equation 5.1 rearranged to make x the subject. The procedure for performing this
computation was written in the Arc Macro Language (AML) of GIS software, Arc
Info, and is included in Appendix 4.
Equation 5-2
?
?
?
?
?
?
?=
vL
D
x 5.0
In this first estimate of Muskingum x, a single routing reach (Muskingum
n = 1) is assumed. The determination of the number of routing reaches however
depends on the value of x. Hence the two parameters must be estimated by
iteration. For the first and second iterations, the respective values of Muskingum x
are computed as
Equation 5-3
?
?
?
?
?
?
?
?
?=
1
1
5.0
vL
D
x
and
138
Equation 5-4
?
?
?
?
?
?
?
?
?=
2
5.0
2
vL
D
x
where
1
x and
2
x are the respective values of Muskingum x
1
L and
2
L are the respective values of routing reach length
However, the only parameters that changes between iterations is the
number of routing subreaches, the Muskingum n, and hence, nLL /
12
=
Substituting for
2
L in Equation 5-4 gives
Equation 5-5
?
?
?
?
?
?
?
?
?=
1
2
5.0
vL
D
nx
which is equivalent to
Equation 5-6 )5.0(5.0
12
xnx ??=
For the more general case in which n is not equal to one in the first
iteration, the value of the Muskingum x can similarly be shown to be
Equation 5-7 ()
i
i
i
i
x
n
n
x ?
?
?
?
?
?
?
?
?
?=
+
+
5.05.0
1
1
where i is the number of iterations
Figure 5-22 shows the respective basin responses for first and second
iterations of the Muskingum x and n parameter computations. HMS basin models
of the Nile were created for the Muskingum parameters resulting from the
respective first and second moments. Negative Muskingum x values which
139
occurred in 14 of the 184 routing reaches were set to zero since HMS does not
allow the parameter x to assume negative values.
Figure 5-22: Basin Response for Iterations of Muskingum x and n Parameters
The second iteration shows an increase in dispersion properties over the
first. Local peaks and troughs are consequently less pronounced in the basin
response after the second iteration of the Muskingum x and n computation. The
basin response also exhibits the bimodal nature of observed for the distributed
parameter case of the STS model in section 5.2.3. Figure 5-23 shows a
comparison of the basin response from the respective Muskingum and dispersion
wave routing approaches using the HMS and STS models.
140
Figure 5-23: Equivalent HMS and STS Basin Response for Distributed
Parameters
A comparison of HMS and STS basin responses in Figure 5-23 reveals
that while it is able to capture the essential characteristics of the response, the
method of moments does not result in an exact duplication of the basin response.
The equivalent Muskingum parameters result in a basin response that is more
dispersed for the shorter reaches (those associated with the first peak) and less
dispersed for the longer reaches. In the Nile basin, the reaches closest to the outlet
have lower dispersion coefficients compared to those transmitting inflow from the
areas upstream of the Sudd Marshes. Hence the equivalent Muskingum
parameters tend to underestimate high dispersion coefficients while
141
overestimating low coefficients. This result is also significant because it reveals
that despite differences in their spatial scale dependence, STS and HMS models
parameterized from a common spatially distributed field of velocity and
dispersion coefficient yield very similar basin responses. The problems of spatial
scale dependence of HMS models can consequently be avoided by determining
routing parameters for various hydrologic elements from spatially distributed
fields of velocity and dispersion coefficient. An approach for deriving these
parameters from observed flow data is demonstrated in section 5.4.1 of this
dissertation.
5.3.4 Effect of Temporal Resolution of HMS Model
The response of a basin in an HMS model is sensitive to the temporal
resolution specified for routing. One source of this sensitivity is the subbasin lag
time which is computed as a function of the routing interval. If Muskingum
method is specified as the routing option within the river reaches, then the number
of routing subreaches, Muskingum n, into which a river reach is divided also
varies with temporal resolution as indicated by Equation 4-14 in section 4.2.3.
While the lag time within each river reach is independent of temporal resolution,
changing the number of routing subreaches results in changes in the dispersion
properties of the reach. The magnitude of this change is discussed in section 5.3.3
of this dissertation. An analysis of the relative importance of subbasin lag time
and Muskingum n in influencing basin response is presented in this section.
142
Additional sources of temporal scale dependence are also studied by performing
runs in which the above mentioned factors are kept constant.
Subbasin lag times in HMS are computed as the larger of 0.6 times the
length of the longest subbasin flow path divided by the flow path velocity and 3.5
times the routing interval, as described in section 4.2.3. This method of
computation ensures that subbasin lag times do not fall below the specified
minimum of 3.5 times the routing interval. The smaller the subbasins defined in
the study area, the more frequently the minimum lag time rule is enforced.
Defining larger subbasins does not solve the problem either since HMS allows a
maximum subbasin lag time of 500 hours. Rather than trying to strike a balance
between the number of subbasins and the subbasin lag time during delineation, it
may be easier to select a routing interval which results in a lower minimum
subbasin lag time and hence a lower number of enforced subbasin lag times. The
routing interval should consequently be determined not by the desired temporal
resolution of the output but rather the number of enforced lag times (and
consequently the accuracy of routing).
To quantify the effect of routing interval on overall subbasin response, a
spreadsheet was created using the subbasin lag function at different routing
intervals. The HMS model of the Congo basin, described in section 5.2, was used
in this analysis. For each temporal resolution, a tally was kept of the number of
subbasins in which the minimum lag condition was enforced. The overall mean
lag time per subbasin was computed and compared with the mean lag time for the
case in which no minimum subbasin lag time is enforced. The difference between
143
the two cases is calculated as the effective enforcement. Table 5-7 shows the
effect of enforced lag times at various routing intervals in the Congo basin.
Table 5-7: Effect of HMS Enforced Minimum Subbasin Lag
Routing
Interval
Enforced
Minimum Lag (hrs)
Number of
Enforcements
Mean Lag
Time (hrs)
Effective
Enforcement (hrs)
24 hours 5040 1550 87.9 33.0
18 hours 3780 1265 71.2 16.3
12 hours 2520 614 59.8 4.9
8 hours 1680 263 56.6 1.7
6 hours 1260 179 55.7 0.8
3 hours 630 58 55.0 0.1
1 hours 210 15 54.9 0.0
0.1 hours 21 0 54.9 0.0
It is evident from the analysis above that it would take a very small routing
interval to completely eliminate subbasin lag time enforcements. Such a routing
interval would be computationally inefficient for performing long term routing.
The resulting flows would also not be reliable given that the temporal resolution
of the input runoff is much coarser than that of the output flows. There are at least
two approaches to getting around this limitation. The first would be to create a
basin file using the actual lag times for all subbasins instead of the HMS enforced
minimum lag times. The warning messages issued by HMS would then have to be
ignored when the basin file is subsequently used in simulations. The second
approach would be to accept a few of lag time enforcements in exchange for
improved computational efficiency. An objective approach to determining a
144
suitable routing interval would be to perform an analysis like the one shown in
Table 5-7. The routing interval could be chosen such that the difference between
the resulting subbasin lag time and that of the zero enforcements case is much less
than the temporal resolution of the input data. This difference is labeled as the
effective enforcement in Table 5-7 above. For the example of the Congo basin, an
HMS model with a routing interval of 12 hours would be suitable because it
results in an effective enforcement of 4 hours which is much less than the
resolution of the daily input available for running the model. No significant delay
is expected when a model of this resolution is used. To verify this hypothesis,
HMS basin files of the Congo were created, at different temporal resolutions, for
the respective enforced minimum and unbounded lag time cases. Routing
intervals 12 hourly and 24 hourly were specified for the models. For the first set
of runs, pure lag routing was specified for river reach routing while the
Muskingum method was specified for the second set of runs to allow the scale
dependence of the methods to be studied. Following is an analysis of the results of
the temporal resolution runs. The hydrographs resulting from the run are
described in terms of their respective first moment and second moments. Figures
5-24 (a) and (b) show the results of the pure lag routing runs at 12 and 24 hour
intervals respectively. The figures show that the minimum subbasin rule has the
effect of shifting the whole basin response forward by an almost uniform amount.
145
Figure 5-24(a): Effect of Minimum Subbasin Lag for 24hr-Interval Lag Routing
Figure 5-24(b): Effect of Min. Subbasin Lag for 12hr-Interval Lag Routing
146
From Table 5-8, the difference in first moments between the basin
response with and without minimum subbasin lag enforcement is 0.47 days at a
12-hour routing interval and 0.99 days at the 24-hour interval. A much larger
difference of about 23 days is observed between the basin response at the
respective 12 and 24 hour intervals. The results in Table 5-8 also show a slight
change in dispersion characteristics between runs with minimum subbasin lag
enforcements and runs without it. However, a much larger increase in second
moment, approximately 200%, occurs with the change in routing interval, from 12
to 24 hours.
Table 5-8: Effect of Min. Subbasin Lag Rule on Basin response, Lag Routing
Routing
Interval (hrs)
Min. Subbasin
Lag Enforced
First Moment
(days)
Second Moment
(days
2
)
12 Yes 68.71 494.17
12 No 68.24 505.57
24 Yes 91.25 1546.90
24 No 90.26 1538.11
In the second set of runs, the Pure Lag routing was replaced with the
Muskingum method in river reach elements in order to determine whether the
features observed in the first set of runs were specific to the routing method. As
with the Lag routing runs, routing was performed at 12 and 24 hour intervals with
and without enforcing the minimum subbasin lag rule. Figure 5-25(a) and (b)
show the effect of the minimum subbasin lag rule on basin response with
Muskingum routing in the river reaches.
147
Figure 5-25(a): Effect of Min. Subbasin Lag for 24hr-Interval Muskingum
Routing
Figure 5-25(b): Effect of Min. Subbasin Lag for 12hr-Interval Muskingum
Routing
148
The basin responses for the Muskingum runs, Figure 5-25(a) and (b),
display similar characteristics to those observed for the Lag routing runs. As the
routing interval is increased, there is increased dispersion and the effect of the
subbasin lag rule becomes more pronounced. From Table 5-9, the lag times in the
Muskingum runs are approximately equal to those in the Lag runs and the second
moments are generally higher in the former.
Table 5-9: Effect of Subbasin Lag Rule on Basin response, Muskingum Routing
Routing
Interval (hrs)
Minimum
Subbasin Lag
Enforced
First Moment
(days)
Second Moment
(days
2
)
12 Yes 68.45 493.68
12 No 68.43 494.31
24 Yes 91.25 1580.71
24 No 90.26 1572.00
However, the difference between lag first and second moments at the 12 and 24
hour routing intervals remains the most significant feature of the runs. Figures 5-
26(a) and (b) show the effects of routing interval on the basin response for the
respective Lag and Muskingum runs.
149
Figure 5-26(a): Effect of Routing Interval on Response with Lag Routing
Figure 5-26(b): Effect of Routing Interval on Response with Muskingum Routing
150
For both the Muskingum and Lag routing runs, an increase in routing
interval leads to increased dispersion. This increased dispersion can be attributed
to the stretching of flow at junctions where one flow element transfers flows to
the next downstream element. The discretization of time into routing intervals
ensures that flow is stretched at such intersections. For example, a flow time of
3.4 days could be interpreted in one of three ways when routing flow at a daily
time interval. The flow time could be converted to a whole number of routing
intervals by rounding off to the nearest whole number, 3 days in this example.
The flow time could also be rounded off to the next higher whole number, 4 days
in this example, with the implicit assumption that flows arriving after the third
day can only contribute to the fourth day. The third approach, which is adopted by
HMS, involves resolving the uncertainty in definition of arrival interval by
dividing the flow between the two adjacent intervals, the 3
rd
or the 4
th
days in this
example. For the 3.4-day flow time assumed in the example above, 60% of the
flow is assumed to arrive within the 3
rd
day while the remaining 40% arrives
within the 4
th
day. This subdivision results in the stretching out of flow, and the
magnitude of the stretching out is proportional to the magnitude of the routing
interval over which the flow is split. Larger routing intervals therefore lead to
more dispersed flows.
There is a secondary mechanism that acts to increase dispersion with
increasing routing interval when the Muskingum method is applied. The number
of routing subreaches, the Muskingum n, in each river reach decreases with
increasing routing interval. Consequently, keeping all other factors constant, the
151
amount of dispersion experienced by flow within each routing subreach also
increases. This increase must however be balanced with the decrease in dispersion
due to the reduction in the number of points at which flow is rerouted within each
river reach. (Dispersion at change over points occurs within a river reach every
time flow is passed from one routing subreach to the next). The change in
dispersion with change in routing interval is slightly higher for the Muskingum
runs than pure lag runs in Table 5-9. This indicates that the length of the routing
sub-reaches is a more important factor in determining routing reach dispersion
characteristics than the number of internal change points within the reach.
5.4 COMPARISON WITH OBSERVED FLOWS
The concept of a basin response, used in sections 5.2 and 5.3 of this
dissertation, allows the effect of modeling parameters to be studied. Changes in
the basin response provide an indication of how changes in a given model
parameter influence flow characteristics. However, a basin response is not a
realistic representation of the flow transformation process in a continental scale
basin because it incorporates an assumption of uniformly distributed input. In
large basins, the occurrence of rainfall is non-uniform both in space and in time.
The validation of the transformation process of a continental scale routing model
must consequently be done by comparison with observed flow data. Simulation
runs to compare routed flows with observed flows are described in sections 5.4.1
and 5.4.2 for the respective STS and HMS models.
152
5.4.1 Comparing STS Simulated Flows with Observed Flows
The distributed nature of sources in the STS model makes it more suited
for routing spatially distributed input. This would in turn require the use of a land-
atmosphere interaction or other water balance model for input runoff generation.
The use of such a runoff model would result introduce an additional source of
errors making it difficult to determine which set of errors to attribute solely to the
routing model. A single source, whose location was designated by a river gauging
station, was defined for the comparison in this study. Mongalla, a gauging station
located just upstream of the Sudd marshes in the Nile basin, was defined as the
source. Lake Nyong, located about 409 km downstream of Mongalla, was defined
as the sink location. In-stream flows recorded at Mongalla were routed using the
diffusion wave equation to Lake Nyong. Observed flows at Lake Nyong were
then compared with simulated flows at the same station. Because of the presence
of the Bahr el Zeraf and other braided channels between the two stations, the
magnitude of the routed and simulated flows at the station were not compared.
However, the translatory characteristics of flow including velocity and dispersion
coefficient can be deduced from the intercomparison.
Appropriate values of velocity and dispersion coefficient were determined
using the method of moments previously mentioned in section 4.1.4 of this
dissertation. 10 day observed flows for the period 1953 to 1957 (Hurst and
Phillips, 1931-) were used in this analysis. The flow values used in this analysis
are also given in Appendix 5 of this dissertation. As shown in Figure 5-27, both
the inflow hydrograph at Mongalla and the discharge at Lake Nyong were broken
153
into 3 periods each representing a full hydrograph cycle with a duration of one
year.
Figure 5-27: Inflow and Discharge Hydrographs in the Sudd Marshes
Within each of these periods, a constant discharge baseflow separation was
performed to allow for better isolation of the input flow signal. The lowest flow
recorded during the period was defined as the baseflow, and all flows above this
baseflow were defined at runoff for the source at Mongalla. Likewise, a baseflow
for each period was defined for the discharge at Lake Nyong, and flows in excess
of this were defined as sink discharges. Figures 5-28 (a), (b) and (c) show the
excess inflows and discharges for periods 1, 2 and 3, respectively.
154
Figure 5-28: Excess Input and Discharge Runoff after Baseflow Separation for (a)
Period 1, (b) Period 2 and (c) Period 3
155
Table 5-10 below shows the duration of the periods, the 1
st
and 2
nd
moments of the excess flows and the velocity and dispersion coefficients resulting
from the moments. Velocity is computed by dividing the difference in 1
st
moments of the input and discharge by the distance between the source and the
sink (409 km). Variance along the flow path is computed as the difference
between the 2
nd
moments of the input and discharge about their respective means.
Dispersion coefficient is then computed from the variance using Equation 4-12 in
section 4.1.4.
Table 5-10: Routing Parameters Computed by Method of Moments
Period Period 1 Period 2 Period 3
Start 440 800 1160
End 790 1150 1510
Input Baseflow 43.9 47.2 47
Discharge Baseflow 38.5 42.7 40.8
Input 1st Moment about origin 173.5046 191.4872 185.4032
Discharge 1st Moment about origin 217.6234 234.3801 227.1073
Input 2nd Moment about mean 5184.534 3208.575 5134.602
Discharge 2nd Moment about mean 5532.649 3803.898 5584.91
Velocity (m/s) 0.107297 0.110363 0.113509
Dispersion Coefficient (m
2
/s) 3924.246 7302.953 6010.043
Figures 5-29 (a), (b) and (c) show a comparison of observed excess
discharge at Lake Nyong (shown as black bars) with the simulated discharge
(shown as gray bars) using the velocity and dispersion parameters computed for
the respective periods in Table 5-10. The input at Mongalla for each period is
represented by a continuous gray line.
156
Figure 5-29: Routing Excess Runoff using Computed Velocities and Dispersion
Coefficients for (a) Period 1, (b) Period 2 and (c) Period 3
157
The mean values of the velocity and dispersion coefficient were computed from
the three period values as 0.11 m/s and 5746 m
2
/s, respectively. Figure 5-30
shows the effect of routing the entire input at Mongalla to Lake Nyong using the
diffusion wave equation with the computed mean velocity and dispersion
coefficient. In the figure, the light gray line represents the input at Mongalla while
the thicker gray line is the simulated discharge at Lake Nyong. The black line
represents the observed flows at Lake Nyong.
Figure 5-30: Effect of Routing Input Flow using Mean Velocity and Dispersion
Coefficients
158
From the diagram, it can be seen that despite differences in magnitude, the
transformation with mean parameters is able to capture key features of the
observed flow at Lake Nyong. The simulated and observed hydrographs peak at
about the same time indicating that the mean velocity is representative of the
translation velocity in the field. The width of the hydrograph cycles of the
simulated and observed cases are also approximately the same indicating that the
dispersion coefficient results in the right amount of spread. It can there be
concluded that the method of moments results in parameters that ensure a
representative transformation of inputs at a source to discharges at a sink. Hence a
STS model for which velocity and dispersion coefficient are estimated by the
method of moments results in a good representation of the flow process in a real
hydrologic system.
5.4.2 Comparing HMS Simulated Flows with Observed Flows
Having compared performance of the respective STS and HMS models
under various parametric conditions, a set of runs was required to compare
simulated flows from the model with observed data. For these comparisons, the
routing process had to be isolated from the soil water balance process which is
typically used to generate input runoff for routing models. Without such isolation
of the two processes, it would be impossible to determine from the resulting flows
which errors to attribute to the routing process and which to attribute to the water
balance process. Hence model verification in this study was done by routing
159
observed flows at an upstream gauging location to a downstream location using
an HMS model and comparing the results with observed flows at the downstream
location. The quantity of simulated flow in this test is not expected to equal the
observed flow exactly since flow entering the river system between the two
gauging locations is not taken into account. However, the flows simulated by the
model should duplicate the major features of the observed hydrograph at the
downstream gauge if upstream stations with comparatively high flows as chosen
as the input sources. Comparisons of this nature were performed for the HMS
model and the results are described in this section. The simulations were
conducted in the Nile basin because of the availability of flow data at various
locations throughout the basin. Upstream source stations were located at
Mongalla and Roseires in the White and Blue Nile respectively. Downstream
observation points were located at Malakal, Khartoum and Aswan. Figure 31 (a)
shows a map of the Nile Basin with the gauge and source locations indicated. Part
(b) of the figure shows the HMS representation of the basin.
160
Figure 5-31: Nile Basin (a) Gauge Locations and (b) HMS Model Components
Malakal is the first major gauging station that captures the flows coming
out of the Sudd Marshes in Southern Sudan. The high dispersion coefficients and
low velocities encountered in the region ensure that flow undergoes significant
redistribution between Mongalla and Malakal. In addition, flows from the small
Bahr el Ghazel and Sobat tributaries of the White Nile are neglect in this analysis.
161
Nevertheless, the simulated hydrograph at Malakal generally arrives at the same
time and has similar flow distribution as the observed flows at the same station.
Figure 5-32 shows the match between the simulated and observed flows at
Malakal.
Figure 5-32: Observed and HMS Simulated Flows at Malakal
Simulated flows match observed flows almost exactly at Khartoum
because of the relatively short time of translation for flow between Roseires and
Khartoum. Steep slopes in the region result in high flow velocities with low in-
stream dispersion rates. The routing model was consequently able to simulate
162
flow between the two stations almost exactly. Figure 5-33 shows the match
between observed and HMS simulated flows at Khartoum.
Figure 5-33: Observed and HMS Simulated Flows at Khartoum
Aswan is located downstream of the intersection between the White and
Blue Nile tributaries. Routed flows from the two tributaries are superimposed at
the station. Figure 5-34 shows a comparison between the simulated and observed
flows at Aswan. The simulated hydrographs capture the rise and fall of the
observed flow hydrographs well though the simulated flows are generally smaller
in magnitude.
163
Figure 5-34: Observed and HMS Simulated Flows at Aswan
The results of these simulations indicate that an HMS model can be used
to route flow in a continental scale basin like the Nile. It allows flows to be routed
over long distances through various types of topographic conditions, and routed
flows from various upstream sources can be superimposed at location further
downstream. The use of HMS as a baseline model for comparisons with the STS
model is consequently justifiable.
164
CHAPTER 6
Conclusions and Recommendations
6.1 CONCLUSIONS
The measure of the success in any study is the extent to which the research
objectives have been met. In this section, the objective of the study as laid out in
section 1.2 of this dissertation are revisited to determine the extent to which they
have been met. The lessons learnt in meeting each objective are subsequently
summarizedinsections6.1.1to6.1.3.
The first research objective was stated as follows:
?to develop a GIS database to support large-scale surface water routing
globally?
The global terrain analysis performed with 1 kilometer Digital Elevation
Models of the earth resulted in corrected terrain grids, flow direction grids, flow
accumulation grids, flow length grids and drainage basins of the whole world.
These data were subsequently used to parameterize two surface water models for
continental scale.
The second research objective was stated as follows:
?to implement a modeling framework which incorporates basin
boundaries in a grid based model while maintaining computational efficiency by
only performing routing at desired locations?
165
A routing model which subdivided the earth into a series of 0.5 by 0.5
degree boxes was implemented globally. The modeling framework allowed for
the incorporation of drainage basin boundaries and modeling unit boundaries
defined in water balance model by intersecting vector representations of the
respective boundaries to define new modeling units called sources. The modeling
framework was implemented in the routing of runoff generated by a Global
Circulation Model (GCM) to terminal points (also called sinks) at the continental
margin and in internally draining catchments. Computational efficiency was
ensured by routing flows from directly from the sources to their respective sinks.
The third objective was stated as follows:
?to examine of the robustness of the source to sink approach as compared
to the watershed based approach in continental scale applications?
The Source to Sink (STS) routing model was compared with the
watershed modeling approach of the Hydrologic Modeling System (HMS) in the
Congo and Nile basins, two of the largest river basins in Africa. The basin
responses computed by the two modeling approaches were compared for various
spatial and temporal resolutions, and for uniform and non-uniformly distributed
routing parameters.
The three objectives of the study have therefore been met. The key
inferences made during each phase of the research are described in sections 6.1.1
to 6.1.3 below.
166
6.1.1 Inferences from Global Terrain Analysis
Due to the high number and varied characteristics of pits in DEMs, the
differentiation of inland catchments from erroneously created pits is best done by
referring to other secondary sources of information such as cartographic maps and
databases of lakes. Inland catchments identified from such secondary sources can
subsequently be preserved in the DEM.
In creating DEMs which correctly represent flow direction, the insertion
of a single null value cell at the lowest elevation of an inland catchment is
sufficient to preserve the boundaries of the inland catchment as well as the flow
direction of all cells that drain towards it. Large inland water bodies such as the
Caspian Sea can also be preserved in the DEM using a single null value cell.
Flow direction grids can be checked for the presence of pits by computing
a histogram of the cell values. The presence of cell values other than powers of 2
from 0 to 8 implies that there are pits in the DEM.
The eight direction pour point model used in determine flow direction in
this study is biased in favor of the four cardinal directions. However, this bias
does not adversely impact the location of river networks and basin boundaries
delineated from the resulting flow direction grid.
The definition of the study area in a large scale terrain analysis must take
into account the location of major drainage basin boundaries as well as the
computer processing resources available for the study. In this study for example,
the inclusion of the drainage basin of the Caspian Sea in Asia would result in a
167
DEM with well over 50 million grid cells. A DEM of this size would be too large
to be processed using the computer resources available for the study.
After eliminating erroneously created pits, delineations performed using
the 30 arc-second DEM resulted in river networks and drainage basins that
corrected represented the actual locations of the respective features. While there
were some minor deviations between delineated and actual drainage features, no
blunders which would severely impact the results of hydrologic modeling
applications developed using these derived datasets were identified. The 30 arc-
second DEMs are therefore suitable for continental scale hydrologic delineation
and parameter derivation.
6.1.2 Inferences from Global STS Model Implementation
If the modeling units of a source to sink model are defined in a vector
environment (as opposed to a raster environment) then they can assume any
shape. Hence a source defined on a rectangular mesh (in geographic coordinates)
can be split up into smaller sources if the drainage basin boundary passes through
it. Thus a source to sink model constructed in a vector environment can retain the
rectangular land surface subdivision used in other grid-based models while
incorporating the exact location of basin boundaries delineated from the DEM.
After the processing of underlying DEMs, source to sink models of any
resolution can be developed with relative ease. The process involves dividing the
continental margin into desired segments, delineating drainage basins from these
segments and subdividing the basins into sources whose parameters are then
168
determined from the processed DEM using zonal averaging and summation
functions. This is a much less involved process than that required to change the
resolution of a cell to cell model or a watershed based model.
The use of the diffusion wave equation instead of the linear reservoir
equation in continental scale cell to cell routing models could significantly reduce
the scale dependence of these models.
Links between routing models and atmospheric or soil water balance
models which provide input runoff for the routing process can be established by
spatially intersecting the respective modeling units. Such an intersection results in
modeling units that are uniquely linked to both the input and routing models by
key fields.
Links between routing models and ocean models which receive routed
flows can be established by using ocean modeling units located at the continental
margin as outlets in the drainage basin delineation. All modeling units defined
within each drainage basin can subsequently be linked to the single receiving
ocean model unit used to delineate the basin.
The use of non-uniform velocity and dispersion parameters in continental
scale routing applications can significantly alter the response of a river basin. It
may expose important hydrograph features that cannot be depicted with a uniform
parameter model. The STS model implemented in this study is an improvement
over other STS models because it can account for spatially distributed velocity
and dispersion parameters.
169
6.1.3 Inferences from HMS Model Development
The preprocessor used in the delineation of HMS hydrologic elements,
CRWR Prepro, was able to accurately classify and derive the connectivity among
3150 hydrologic elements at the continental scale. The preprocessor was also able
determine hydrologic routing parameters for each of the elements based on the
DEM and additional parameters supplied.
While the preprocessor ensures that subbasin routing parameters are
within the ranges permitted by HMS, it does not perform a similar check for river
reach routing parameters. River reach parameters that are outside the constraints
set by HMS can result from the preprocessing and may require manual editing to
correct.
The process of deriving HMS basin files from DEMs at the continental
scale is memory intensive and time consuming. Consequently, care must be taken
in defining the modeling parameters such as delineation threshold, flow velocity
and other routing parameters to be used in creating the model.
In working with HMS at the continental scale, it was determined that some
of its imposed constraints on routing parameters, particularly those relating to the
maximum river reach lag time and the maximum number of routing subreaches
for a Muskingum reach, could be removed without affecting the accuracy of the
routed flows.
Given a field of spatially distributed flow velocities and dispersion
coefficients, the weighted flow length function can be used in conjunction with
170
and the method of moments to determine routing parameters for hydrologic
elements in an HMS model.
An HMS model is sufficient to represent the processes observed in flow at
the continental scale. It is able to perform flow routing over long distances,
through various topographic conditions and aggregate flow from various sources.
It is therefore a suitable model to use for model intercomparisons at the
continental scale.
6.1.4 Inferences from Comparing STS and HMS Models
The hydrologic processes exhibited by an HMS model were successfully
duplicated using a STS model. A better match was obtained when a finer
delineation threshold was specified for the HMS model. At the finer threshold
input runoff spends more time within the flow channels than on the land surface,
making it similar to the single flow path approach of the STS model.
Any model that involves subdividing the earth into a series of linked
segments as in the watershed based and cell to cell models, will result in
numerical dispersion. This is because at the end of each segment, an essentially
continuous response must be discretized into time steps. The discretization results
in the redistribution of flow between two or more adjacent time steps at the end of
each segment. The larger the number of segments, the larger the net dispersion
experienced by flow arriving at the discharge outlet.
The discretization of continuous flow at the end of a segment occurs only
once in a source to sink model irrespective of the resolution of the model. In a
171
watershed based model, discretization occurs every time flow is passed from one
flow element to another. The choice of spatial resolution of a watershed based
model greatly influences the dispersion properties of the model while the
dispersion properties of a source to sink model are relatively unaffected by the
choice of spatial scale.
Redistribution of flow due to discretization is also influenced by the
choice of routing interval. A large routing interval results in greater redistribution
every time flow is discretized. The dispersion characteristics of a watershed based
model are consequently impacted more by the choice of temporal resolution than
those of a source to sink model in which redistribution only occurs once.
Velocity and dispersion parameters in a river basin can be derived from
observed flow data by equating the difference between the moments of flow at
two gauging stations to the moments of the flow routing method used in the
model. The method of moments is consequently a good approach for deriving
spatially distributed velocity and dispersion parameters.
Grids of spatially distributed velocity and dispersion coefficients can be
used to parameterize both a watershed based model using the Muskingum
equation for in-stream routing and a dispersion type source to sink model,
irrespective of resolution. The resulting basin responses for the STS and HMS
models are very similar. Hence, parameterizing watershed based models from
spatially distributed grids of velocity and dispersion coefficient can reduce their
spatial scale dependence.
172
6.2 RECOMMENDATIONS
The development of a global grid of flow velocity and dispersion
coefficients would greatly advance the development of continental scale routing
models. The derivation of these parameters from observed flow data using the
method of moments offers a promising avenue to proceed with this development.
The approach lends itself to automated implementation using network algorithms
since distances and hierarchy along a river network can be determined
automatically. Further research could be conducted to investigate the possibility
of implementing this approach globally using a river network model with flow
gauges and associated flow time series along it.
In testing the longitudinal decomposability of the diffusion type source to
sink model in this study, a diffusion type cell to cell routing model was
developed. The model offers all the advantages of a conventional cell to cell
model while simultaneously affording many of the advantages of diffusion wave
routing such as longitudinal decomposability as well as spatial and temporal scale
independence. The implementation of the diffusion type cell to cell routing and its
comparison with the linear reservoir type cell to cell model would provide
valuable insight into its viability at the continental scale.
The assumption of time invariant flow velocities and dispersion
coefficients may be insufficient to represent the flow regime in some regions.
Further research is required to study the effect of temporally varying routing
parameters on the basin response. In regions where time invariant parameters are
shown to be insufficient, alternate flow paths between a source and a sink could
173
be defined to allow flow to be routed using different response functions at
different flow stages. Further research is required to develop ways of partitioning
flow into different flow paths to allow the continued application of the source to
sink method in such regions.
Flow regimes in nature are often interrupted by the presence of reservoirs
and other control structures with non-linear release schedules which cannot be
included in the linear source to sink model. A secondary source to sink model
could be used to route upstream flow into the reservoir, followed by the
application of the reservoir regulation rules to determine reservoir discharges.
Such a nested source to sink approach would allow the reservoir and the sources
located upstream of it can be treated a single source whose input could
subsequently be routed in the primary source to sink model as before. Research
towards the implementation of such an approach and its testing in a large basin
would further enhance the ability of the source to sink model to represent flow
regimes encountered at the continental scale.
The use of the diffusion wave equation for routing was found to result in
significantly less scale dependent models. Its inclusion as a routing option in a
watershed based model like HMS would result in a much more versatile model.
Further research should be conducted to determine the scale dependence of a
watershed based model employing the diffusion wave equation for in-stream
routing.
174
APPENDICES
Contents
Appendix 1: AML Code for Source to Sink Model Preprocessing
Appendix 2: Parameters of the T42 mesh
Appendix 3: FORTRAN Code for Source to Sink Routing Model
Appendix 4: AML Code for Computing Equivalent Muskingum Parameters from
a Field of Velocities and Dispersion Coefficients
Appendix 5: 10 day flows in the Nile basin, 1953-57
175
Appendix 1: AML Code for Source to Sink Model Preprocessing
/* -------------------------------------------------------------
/* STSPREPRO.AML
/* ARC MACRO LANGUAGE (AML) CODE FOR PREPROCESSING STS MODELING UNITS
/* -------------------------------------------------------------
/* BY FRANCISCO OLIVERA AND KWABENA ASANTE
/* CRWR, UNIV. OF TEXAS AT AUSTIN
/* WITH EXTRACTS FROM CODE BY FERDINAND HELLWEGER
/* -------------------------------------------------------------
/* THIS PROGRAM TAKES USGS FLOW DIRECTION GRIDS (HYDRO1K DATASET) AND
/* USES THEM TO DEFINE SINKS, DELINEATE RIVER BASINS, DEFINE AND
/* PARMETERIZE SOURCES FOR THE STS MODEL.
/* -------------------------------------------------------------
/* INPUTS INCLUDE:
/* USGS FLOW DIRECTION GRID -> AF_FD
/* OCEAN 5-DEGREE MESH COVERAGE -> SEALAMB
/* LAND SURFACE MESH (0.5-DEGREE) COVERAGE -> MESHAFRLAM
/* SOIL WATER BALANCE UNITS MESH COVERAGE -> T42AFR
/* -------------------------------------------------------------
/************************************************************
/* INITIAL SETTINGS
SETWINDOW AF_FD AF_FD
SETCELL AF_FD
/* ASSIGNING NODATA TO THE SINK CELLS (INTERNAL AND EXTERNAL)
FDR0=CON(AF_FD <> 55537, AF_FD)
SINKS = SINK(FDR0)
SINKSZERO = ISNULL(SINKS)
FDR = FDR0/SINKSZERO
/* CREATING THE CONTINENT BORDER (VECTOR AND RASTER)
BORDER = GRIDPOLY(CON(FDR,1))
QUIT
CLEAN BORDER BORDERLN ##LINE
GRID
/* INITIAL SETTINGS
SETWINDOW AF_FD AF_FD
SETCELL AF_FD
SEASHORE1=CON(LINEGRID(BORDERLN),1)
SEASHORE_SW = CON(SEASHORE1(0,-1) == 1 AND SEASHORE1(1,0) == 1,1)
SEASHORE_SE = CON(SEASHORE1(0,-1) == 1 AND SEASHORE1(-1,0) == 1,1)
SEASHORE_NE = CON(SEASHORE1(0,1) == 1 AND SEASHORE1(-1,0) == 1,1)
SEASHORE_NW = CON(SEASHORE1(0,1) == 1 AND SEASHORE1(1,0) == 1,1)
176
SEASHORE2=
MERGE(SEASHORE1,SEASHORE_SW,SEASHORE_SE,SEASHORE_NE,SEASHORE_NW)
SEAMESH = POLYGRID(SEALAMB, SEALAMB-ID)
SEASHORE = SEASHORE2*SEAMESH
SEA_WSH = WATERSHED(FDR,SEASHORE)
SEA_WSH_PC = GRIDPOLY(SEA_WSH)
QUIT
/* WHAT ABOUT BUILDING THE TOPOLOGY OF SEA_WSH
TABLES
SELECT SEALAM.PAT
ALTER AREA
AREACOARSE
~
~
~
~
SELECT MESHAFRLAM.PAT
ALTER AREA
AREAFINE
~
~
~
~
QUIT
INTERSECT MESHAFRLAM SEA_WSH_PC MESHFINEWSH
QUIT
TABLES
SELECT MESHFINEWSH.PAT
RESELECT GRID-CODE <0
PURGE
Y
QUIT
/* CREATE THE ZONE GRID FOR THE ZONALSTATS COMMAND
/* USE 500 INSTEAD OF 1000 FOR BETTER POLYGON DEFINITION
POLYGRID MESHFINEWSH FINEWSHGRID MESHFINEWSH-ID
500
Y
GRID
/* INITIAL SETTINGS
SETWINDOW AF_FD AF_FD
SETCELL AF_FD
/* COMPUTE MEAN FLOWLENGTH FROM FLDS USING FINEWSHGRID TO DEFINE ZONES
177
MEANLENGRID = ZONALMEAN ( FINEWSHGRID , FLDS , DATA )
/* CONVERT TO AN INTEGER GRID
MEANLEN = INT ( MEANLENGRID )
BUILDVAT MEANLEN
/* ASSIGN 999999999999 TO NODATA CELLS
MEANLEN2=CON (ISNULL(MEANLEN), 999999999999, MEANLEN)
BUILDVAT MEANLEN2
QUIT
/* MOVE POLYGON LABELS TO CENTROID BUT ENSURE ALL LABELS WITHIN POLYGON
CENTROIDLABELS MESHFINEWSH INSIDE
/* ADD X AND Y COORDINATES TO PAT AND CREAT FIELD SPOT
ADDXY MESHFINEWSH
DROPITEM MESHFINEWSH.PAT MESHFINEWSH.PAT SPOT
ADDITEM MESHFINEWSH.PAT MESHFINEWSH.PAT SPOT 16 16 B 3
/* GO INTO ARC PLOT TO READ VALUES FROM GRID TO POLYGON
/* DISPLAY 1040 DOES NOT REQUIRE A VISUAL DISPLAY, ONLY AN OUTPUT FILE
(DISTORE)
/* USEFUL FOR REMOTE EXECUTION
AP
DISPLAY 1040
DISTORE
MAPE MESHFINEWSH
/* RESEL MESHFINEWSH POLY OVERLAP BORDER POLYGON WITHIN
/* THIS LINE IS NOT REQUIRED WHEN USING MEANLEN2
/* DECLARE AND OPEN A CURSOR TO READ AND WRITE TO THE PAT
CURSOR PTCUR DECLARE MESHFINEWSH POLY RW
CURSOR PTCUR OPEN
CURSOR PTCUR FIRST
&S OLD$ECHO [SHOW &ECHO]
/*START A LOOP TO GO THROUGH THE PAT, FIND THE CELL VALUE AT EACH
/*POINT LOCATION, AND WRITE IT TO THE SPOT ITEM IN THE PAT.
&DO &WHILE %:PTCUR.AML$NEXT%
&S :PTCUR.SPOT := [SHOW CELLVALUE MEANLEN2[SHOW SELECT MESHFINEWSH POLY
1 XY]]
CURSOR PTCUR NEXT
&END
CURSOR PTCUR CLOSE
178
Q /*QUIT FROM ARCPLOT
&DATA ARC INFO; ARC
SEL MESHFINEWSH.PAT
ALTER SPOT
FLOWLENGTH,,,,,,,,
Q STOP
&END
&RETURN /*TO CALLING AML
179
Appendix 2: Parameters of the T42 mesh
(SUPPLIED BY WILLIAM G. LARGE,CLIMATE SYSTEM MODEL GROUP,
NATIONAL CENTER FOR ATMOSPHERIC RESEARCH)
THE FOLLOWING ARE THE PARAMETERS OF THE T42 MESH ON WHICH THE COMMUNITY
CLIMATE MODEL (CCM3) GENERATES INPUT RUNOFF.
I AND J ARE COORDINATE INDICES FROM INDICATING THE LOCATION OF THE RUNOFF
GENERATING UNIT
XEDGE AND YEDGE ARE THE LONGITUDE AND LATITUDE OF THE LOWER LEFT CORNER
OF EACH RUNOFF GENERATING UNIT (VALUES OF ?1 ARE DUMMY VALUES)
XMID AND YMID ARE THE LONGITUDE AND LATITUDE OF THE CENTROID OF EACH
RUNOFF GENERATING UNIT (VALUES OF ?1 ARE DUMMY VALUES)
I/J XEDGE XMID YEDGE YMID
1 -1.406 0 -90 -87.864
2 1.406 2.813 -86.48 -85.097
3 4.219 5.625 -83.705 -82.313
4 7.031 8.438 -80.919 -79.526
5 9.844 11.25 -78.131 -76.737
6 12.656 14.063 -75.342 -73.948
7 15.469 16.875 -72.553 -71.158
8 18.281 19.688 -69.763 -68.368
9 21.094 22.5 -66.973 -65.578
10 23.906 25.313 -64.182 -62.787
11 26.719 28.125 -61.392 -59.997
12 29.531 30.938 -58.602 -57.207
13 32.344 33.75 -55.811 -54.416
14 35.156 36.563 -53.021 -51.626
15 37.969 39.375 -50.23 -48.835
16 40.781 42.188 -47.44 -46.045
17 43.594 45 -44.649 -43.254
18 46.406 47.813 -41.859 -40.464
19 49.219 50.625 -39.068 -37.673
20 52.031 53.438 -36.278 -34.883
21 54.844 56.25 -33.487 -32.092
22 57.656 59.063 -30.697 -29.301
23 60.469 61.875 -27.906 -26.511
24 63.281 64.688 -25.115 -23.72
180
25 66.094 67.5 -22.325 -20.93
26 68.906 70.313 -19.534 -18.139
27 71.719 73.125 -16.744 -15.348
28 74.531 75.938 -13.953 -12.558
29 77.344 78.75 -11.162 -9.767
30 80.156 81.563 -8.372 -6.977
31 82.969 84.375 -5.581 -4.186
32 85.781 87.188 -2.791 -1.395
33 88.594 90 0 1.395
34 91.406 92.813 2.791 4.186
35 94.219 95.625 5.581 6.977
36 97.031 98.438 8.372 9.767
37 99.844 101.25 11.162 12.558
38 102.656 104.063 13.953 15.348
39 105.469 106.875 16.744 18.139
40 108.281 109.688 19.534 20.93
41 111.094 112.5 22.325 23.72
42 113.906 115.313 25.115 26.511
43 116.719 118.125 27.906 29.301
44 119.531 120.938 30.697 32.092
45 122.344 123.75 33.487 34.883
46 125.156 126.563 36.278 37.673
47 127.969 129.375 39.068 40.464
48 130.781 132.188 41.859 43.254
49 133.594 135 44.649 46.045
50 136.406 137.813 47.44 48.835
51 139.219 140.625 50.23 51.626
52 142.031 143.438 53.021 54.416
53 144.844 146.25 55.811 57.207
54 147.656 149.063 58.602 59.997
55 150.469 151.875 61.392 62.787
56 153.281 154.688 64.182 65.578
57 156.094 157.5 66.973 68.368
58 158.906 160.313 69.763 71.158
59 161.719 163.125 72.553 73.948
60 164.531 165.938 75.342 76.737
61 167.344 168.75 78.131 79.526
62 170.156 171.563 80.919 82.313
63 172.969 174.375 83.705 85.097
64 175.781 177.188 86.48 87.864
65 178.594 180 90 -1
66 181.406 182.813 -1 -1
181
67 184.219 185.625 -1 -1
68 187.031 188.438 -1 -1
69 189.844 191.25 -1 -1
70 192.656 194.063 -1 -1
71 195.469 196.875 -1 -1
72 198.281 199.688 -1 -1
73 201.094 202.5 -1 -1
74 203.906 205.313 -1 -1
75 206.719 208.125 -1 -1
76 209.531 210.938 -1 -1
77 212.344 213.75 -1 -1
78 215.156 216.563 -1 -1
79 217.969 219.375 -1 -1
80 220.781 222.188 -1 -1
81 223.594 225 -1 -1
82 226.406 227.813 -1 -1
83 229.219 230.625 -1 -1
84 232.031 233.438 -1 -1
85 234.844 236.25 -1 -1
86 237.656 239.063 -1 -1
87 240.469 241.875 -1 -1
88 243.281 244.688 -1 -1
89 246.094 247.5 -1 -1
90 248.906 250.313 -1 -1
91 251.719 253.125 -1 -1
92 254.531 255.938 -1 -1
93 257.344 258.75 -1 -1
94 260.156 261.563 -1 -1
95 262.969 264.375 -1 -1
96 265.781 267.188 -1 -1
97 268.594 270 -1 -1
98 271.406 272.813 -1 -1
99 274.219 275.625 -1 -1
100 277.031 278.438 -1 -1
101 279.844 281.25 -1 -1
102 282.656 284.063 -1 -1
103 285.469 286.875 -1 -1
104 288.281 289.688 -1 -1
105 291.094 292.5 -1 -1
106 293.906 295.313 -1 -1
107 296.719 298.125 -1 -1
108 299.531 300.938 -1 -1
182
109 302.344 303.75 -1 -1
110 305.156 306.563 -1 -1
111 307.969 309.375 -1 -1
112 310.781 312.188 -1 -1
113 313.594 315 -1 -1
114 316.406 317.813 -1 -1
115 319.219 320.625 -1 -1
116 322.031 323.438 -1 -1
117 324.844 326.25 -1 -1
118 327.656 329.063 -1 -1
119 330.469 331.875 -1 -1
120 333.281 334.688 -1 -1
121 336.094 337.5 -1 -1
122 338.906 340.313 -1 -1
123 341.719 343.125 -1 -1
124 344.531 345.938 -1 -1
125 347.344 348.75 -1 -1
126 350.156 351.563 -1 -1
127 352.969 354.375 -1 -1
128 355.781 357.188 -1 -1
129 358.594 -1 -1 -1
183
Appendix 3: FORTRAN Code for Source to Sink Routing Model
!ROUTDISP VERSION 1.0
! *********************************************
!THIS PROGRAM TAKES RUNOFF GENERATED BY A VERTICAL BALANCE AND ROUTES
!IT TO THE CONTINENTAL MARGIN.THIS VERSION OF THE PROGRAM ACCOUNTS FOR
!BOTH GEOMORPHOLOGICAL AND HYDRODYNAMIC DISPERSION.THIS VERSION ALLOWS
!FOR THE USE OF SPATIALLY VARIABLE FLOW VELOCITY, DISPERSION COEFFICIENT
!AND LOSS COEFFICIENT
!
!INPUTS :
! SINKCELL.TXT - LIST OF GRID CODES OF THE OUTLET CELLS
!
! PARAMETER.TXT - LIST OF PARAMETERS (COMMA DELIMITED, NO HEADER)
! PARAMETER(1) = NUMBER OF TIME STEPS TO BE SIMULATED
! PARAMETER(2) = INTERVAL FOR RUNOFF GENERATION IN SECONDS
! PARAMETER(3) = NUMBER OF OUTLETS OR SINK CELLS
! PARAMETER(4) = NUMBER OF DRAINAGE POLYGONS OR SOURCES
!
! LANDCELL.TXT - LIST OF ATTRIBUTES OF SOURCES (COMMA DELIMITED)
! WITH A SINGLE LINE ASSIGNED FOR EACH SOURCE.
! ATTRIBUTES INCLUDE
! SOURCE-ID DEFINED BY THE SOURCE COVERAGE POLYGON ID
! AREA OF THE SOURCE IN METERS SQUARED
! SINK-ID DEFINED BY GRID-CODE OF DRAINAGE BASINS
! FLOW LENGTH FROM SOURCE TO SINK IN METERS
! VELOCITY ALONG THE FLOW PATH IN METERS PER SECOND
! DISPERSION COEFFICIENT IN METERS SQUARED PER SECOND
! LOSS COEFFICIENT IN FRACTION OF FLOW LOST PER SECOND
! KEY FIELD LINKING SOURCES TO INPUT RUNOFF UNITS
!
!
! INPUT RUNOFF FILES - H0001.TXT TO H3700.TXT
! ONE FILE FOR EACH SIMULATION TIME STEP CONTAINING
! AN X-COORDINATE INDEX FROM 1 TO 128
! AY-COORDINATE INDEX FROM 1 TO 64
! INPUT RUNOFF IN MILLIMETERS PER SECOND
!
! INFILES.TXT - LISTS NAMES OF INPUT RUNOFF FILES
!
!OUTPUTS :
!
! INFLOW.TXT - RUNOFF VOLUME (M3) ENTERING EACH BASIN PER TIME STEP
!
! SINKFLOW.TXT - DISCHARGE (M3/S) TO EACH SINK PER TIME STEP
!
!
184
! AUTHORSHIP:KWABENA ASANTE, CRWR, UT AUSTIN,JULY, 1999
!
! @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!
!OPENFILESANDREADDATAINTOARRAYS
REAL,DIMENSION(:),ALLOCATABLE::SINKCELL,FLOWTIME
REAL, DIMENSION(:), ALLOCATABLE : : ADESUM,LOSSFACTOR
REAL, DIMENSION(:,:), ALLOCATABLE : : SINKFLOW,RUNOFF,LOSSES
REAL, DIMENSION(:,:), ALLOCATABLE : : LANDCELL,INFLOW,RUNOFFIN
REAL, DIMENSION(:,:), ALLOCATABLE : : COEFFICIENTS,ADECOEFF,LOSSES
REAL :: INTERVAL, FLOWTIMEB,MAXVELOCITY,MAXFLOWTIME,MAXDISPCOEFF
REAL :: MAXLENGTH,ADEPARTA,ADEPARTB, EXPPARTA
REAL, PARAMETER :: PI = 3.1415927
INTEGER, DIMENSION(:), ALLOCATABLE : : T42RUN
INTEGER :: XCOORD(8192) , YCOORD(8192) , T42NUM(8192)
INTEGER :: ALLOCATESTATUS,OPENSTATUS,SINKCELLS,OUTSTEPS,LTIME
INTEGER :: SOURCES,STEPS,INPUTSTATUS,MAXWIDTH,CTIME
CHARACTER(11), DIMENSION(:), ALLOCATABLE : : FILENAME
PRINT *, "OPENING INPUT FILE PARAMETER.TXT.........."
OPEN (UNIT = 10, FILE = "PARAMETERS.TXT", STATUS = "OLD", ACTION
.= "READ", POSITION = "REWIND", IOSTAT = OPENSTATUS)
IF (OPENSTATUS > 0) STOP "***CANNOT OPEN PARAMETERS.TXT FILE***"
READ (10, *, IOSTAT =INPUTSTATUS)STEPS,INTERVAL,SINKCELLS,SOURCES
IF (INPUTSTATUS > 0) STOP "****INPUT ERROR IN PARAMETERS****"
IF (INPUTSTATUS < 0) STOP "****NOT ENOUGH PARAMETERS****"
! PARAMETERS
! PARAMETER(1) = STEPS - NUMBER OF TIME STEPS TO BE SIMULATED
! PARAMETER(2) = INTERVAL - INTERVAL FOR RUNOFF GENERATION
! PARAMETER(3) = SINKCELLS - NUMBER OF SINK CELLS
! PARAMETER(4) = SOURCES - NUMBER OF DRAINAGE POLYGONS CELLS
! ALLOCATE MEMORY FOR ALLOCATEABLE VARIABLES
ALLOCATE (SINKCELL(SINKCELLS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: SINKCELL()***"
ALLOCATE (INFLOW(SINKCELLS,STEPS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: INFLOW()***"
ALLOCATE (RUNOFF(SOURCES,STEPS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: RUNOFF()***"
185
ALLOCATE (RUNOFFIN(8192, STEPS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: RUNOFFIN()***"
ALLOCATE (T42RUN(SOURCES), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: T42RUN() ***"
ALLOCATE (FILENAME(STEPS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: FILENAME()***"
ALLOCATE (ADESUM(SOURCES), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: ADESUM()***"
ALLOCATE (LOSSFACTOR(SOURCES), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "**NOT ENOUGH MEMORY: LOSSFACTOR()**"
ALLOCATE (LOSSES(SINKCELLS,STEPS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: LOSSES()***"
ALLOCATE (LANDCELL(SOURCES,8),STAT=ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: LANDCELL()***"
!LANDCELL(SOURCES,1)=POLYGON IDENTIFICATION NUMBER
!LANDCELL(SOURCES,2)=POLYGON AREA
!LANDCELL(SOURCES,3)=SINK WHERE FLOW IS DISCHARGED(GRID CODE)
!LANDCELL(SOURCES,4)=FLOW LENGTH FROM SOURCE TO SINK
!LANDCELL(SOURCES,5)=VELOCITY FROM SOURCE TO SINK
!LANDCELL(SOURCES,6)=DISPERSION COEFFICIENT ALONG THE FLOWPATH
!LANDCELL(SOURCES,7)=LOSS COEFFICIENT ALONG THE FLOWPATH
!LANDCELL(SOURCES,8)=SOIL WATER BALANCE UNIT FOR INPUT RUNOFF
ALLOCATE (FLOWTIME(SOURCES), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: FLOWTIME()***"
PRINT *, "OPENING INPUT FILE SINKCELL.TXT.........."
OPEN (UNIT = 11 , FILE = "SINKCELL.TXT" , STATUS ="OLD" , ACTION
.= "READ", POSITION = "REWIND" , IOSTAT = OPENSTATUS )
IF (OPENSTATUS > 0) STOP "***CANNOT OPEN SINKCELL.TXT FILE***"
READ (11, *, IOSTAT = INPUTSTATUS)SINKCELL
IF (INPUTSTATUS > 0) STOP "****INPUT ERROR IN SINKCELL****"
IF (INPUTSTATUS < 0) STOP "****NOT ENOUGH SINK CELLS****"
PRINT *, "OPENING INPUT FILE LANDCELL.TXT.........."
OPEN (UNIT = 12 , FILE = "LANDCELL.TXT" , STATUS ="OLD" , ACTION =
. "READ" , POSITION = "REWIND" , IOSTAT = OPENSTATUS )
IF (OPENSTATUS > 0) STOP "****CANNOT OPEN LANDCELL.TXT FILE****"
DO Z = 1, SOURCES,1
READ (12, *, IOSTAT = INPUTSTATUS)LANDCELL(Z, 1),
.LANDCELL(Z, 2), LANDCELL(Z, 3), LANDCELL(Z, 4),
.LANDCELL(Z, 5), LANDCELL(Z, 6), LANDCELL(Z, 7), LANDCELL(Z, 8)
186
END DO
IF (INPUTSTATUS > 0) STOP "****INPUT ERROR IN LANDCELL****"
IF (INPUTSTATUS < 0) STOP "****NOT ENOUGH LAND CELLS****"
!OPENFILESANDREADDATAINTOARRAYS
PRINT *, "OPENING INPUT FILE INFILES.TXT.........."
OPEN (UNIT = 41, FILE = "INFILES.TXT", STATUS = "OLD", ACTION =
."READ", POSITION = "REWIND", IOSTAT = OPENSTATUS)
IF (OPENSTATUS > 0) STOP "***CANNOT OPEN FILES4.TXT***"
PRINT *, "OPENING INPUT FILE AFRUNOFF.TXT.........."
OPEN (UNIT = 43, FILE = "AFRUNOFF.TXT", STATUS = "OLD", ACTION =
."READ", POSITION = "REWIND", IOSTAT = OPENSTATUS)
IF (OPENSTATUS > 0) STOP "***CANNOT OPEN AFRUNOFF.TXT***"
READ (43, *, IOSTAT = INPUTSTATUS)(T42RUN(X), X =1, SOURCES,1)
IF (INPUTSTATUS > 0) STOP "****INPUT ERROR IN T42NUMS****"
IF (INPUTSTATUS < 0) STOP "****NOT ENOUGH T42NUMS****"
!DETERMINE THE MAXIMUM FLOW LENGTH, TIME AND DISPERSED FLOW WIDTH
MAXLENGTH =0
MAXDISPCOEFF =0
DOX=1,SOURCES,1
FLOWTIME(X) = INT(LANDCELL(X,4) / (LANDCELL(X, 5)*INTERVAL))
LOSSFACTOR(X) = EXP(-LANDCELL(X,7)*FLOWTIME(X)*INTERVAL)
IF (FLOWTIME(X) .GT. MAXFLOWTIME)THEN
MAXLENGTH =LANDCELL(X, 4)
MAXFLOWTIME =MAXLENGTH /(LANDCELL(X, 5)*INTERVAL)
MAXVELOCITY =LANDCELL(X, 5)
MAXDISPCOEFF =LANDCELL(X, 6)
ENDIF
END DO
!THIS GIVES THE WIDTH OF THE DISPERSED FLOW IN METERS
!(FISCHER ET AL, 1979, MIXING IN INLAND LAKES AND RESERVOIRS)
MAXWIDTH = INT((SQRT(32*MAXDISPCOEFF*MAXFLOWTIME*INTERVAL))/
.(MAXVELOCITY*INTERVAL))
ALLOCATE (COEFFICIENTS(SOURCES,MAXWIDTH), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "***NOT ENOUGH MEMORY: COEFFS.()***"
ALLOCATE (ADECOEFF(SOURCES,MAXWIDTH), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "**NOT ENOUGH MEMORY: ADECOEFF(,)**"
187
!CREATE AN OUTPUT FILE FOR STORING ROUTED FLOWS TO SINK (OUTLET) CELLS
PRINT *, "CREATING OUTPUT FILE DISPCOEF.TXT.........."
OPEN (UNIT = 44 , FILE = "DISPCOEF.TXT" , STATUS ="UNKNOWN" ,
.ACTION = "READWRITE" , POSITION = "REWIND" , IOSTAT = OPENSTATUS)
IF (OPENSTATUS > 0) STOP "***CANNOT CREATE DISPCOEF.TXT FILE***"
PRINT *, "COMPUTING DISPERSION COEFFICIENTS.........."
DOX=1,SOURCES,1
ADESUM(X) = 0
IF (LANDCELL(X,4).GT.0)THEN
DOY=1,MAXWIDTH,1
FLOWTIMEB = (INT(FLOWTIME(X)) - INT(MAXWIDTH/2)+ Y)* INTERVAL
IF (FLOWTIMEB .GT. 0) THEN
TEMPCALC =4*PI *LANDCELL(X, 6) * FLOWTIMEB ** 3
ADEPARTA =LANDCELL(X, 4)/( TEMPCALC ** 0.5)
EXPPARTA=-((LANDCELL(X,5)*FLOWTIMEB)-LANDCELL(X, 4))**2
ADEPARTB =EXP(EXPPARTA/(4*LANDCELL(X,6)* FLOWTIMEB))
ADECOEFF(X,Y) = ADEPARTA *ADEPARTB
ELSE ! IF (FLOWTIMEB .LE. 0)
ADECOEFF(X,Y) = 0
ENDIF ! (FLOWTIMEB .GT. 0)
ADESUM(X) = ADESUM(X) + ADECOEFF(X, Y)
ENDDO ! Y = 1, MAXWIDTH,1
DO Z = 1, MAXWIDTH,1
COEFFICIENTS(X, Z) = ADECOEFF(X, Z) / ADESUM(X)
ENDDO
ELSE ! IF (LANDCELL(X,4).LE.0)
COEFFICIENTS(X,(INT(MAXWIDTH*0.5)+1)) = 1
!DOZ=2,MAXWIDTH,1
!COEFFICIENTS(X,Z) = 0
! ENDDO
ENDIF ! IF (LANDCELL(X,4).GT.0)
188
ENDDO ! DO X=1,SOURCES,1
!LANDCELL(X,4) IS THE FLOWLENGTH FROM THE UNIT TO THE OUTLET IN M
!COMPUTE FIRST PASSAGE TIMES COEFFICIENTS FOR EACH INTERVAL IN MAXWIDTH
!START AND END AT HALF A MAXWIDTH ON EITHER SIDE OF THE ELEMENT FLOWTIME
!MAKEALLOWANCEFORCASEWHERETHISRANGEGOESPASTTHETIMEOFINPUT
PRINT *, "WRITING OUTPUT FILE DISPCOEF.TXT.........."
DOY=1,SOURCES,1
WRITE (44, 200) (LANDCELL(Y, 1)), (COEFFICIENTS(Y, X), X = 1,
.MAXWIDTH)
END DO
DOX=1,STEPS,1
READ (41, *, IOSTAT = INPUTSTATUS)FILENAME(X)
IF (INPUTSTATUS > 0) STOP "**INPUT ERROR IN FILENAMES**"
IF (INPUTSTATUS < 0) STOP "****NOT ENOUGH FILENAMES****"
!PRINT *, FILENAME(X)
OPEN (UNIT = 42, FILE = FILENAME(X), STATUS = "OLD", ACTION =
. "READ", POSITION = "REWIND", IOSTAT = OPENSTATUS)
IF (OPENSTATUS > 0) STOP "***CANNOT OPEN INPUT FILES***"
DO Y = 1, 8192, 1
READ (42, *, IOSTAT = INPUTSTATUS)XCOORD(Y) , YCOORD(Y),
.RUNOFFIN(Y, X)
IF (INPUTSTATUS > 0) STOP "**INPUT ERROR IN PARAMETERS**"
IF (INPUTSTATUS < 0) STOP "****NOT ENOUGH PARAMETERS****"
T42NUM(Y) = ((YCOORD(Y)*1000) + (129-XCOORD(Y)))
END DO ! Y = 1, 8192, 1
END DO ! X = 1, STEPS,1
DOX=1,SOURCES,1
DO Y = 1, 8192, 1
IF (T42RUN(X).EQ.T42NUM(Y)) THEN
!PRINT*,"MATCH FOUND"
DO Z = 1, STEPS,1
RUNOFF(X,Z) = RUNOFFIN(Y,Z)*86.4
189
END DO ! Z = 1, STEPS,1
END IF
END DO ! Y = 1, 8192, 1
ENDDO !X=1,SOURCES,1
!CREATE AN OUTPUT FILE FOR STORING ROUTED FLOWS TO SINK (OUTLET) CELLS
PRINT *, "CREATING OUTPUT FILE SINKFLOW.TXT.........."
OPEN (UNIT = 15 , FILE = "SINKFLOW.TXT" , STATUS ="UNKNOWN" ,
.ACTION = "READWRITE" , POSITION = "REWIND" , IOSTAT = OPENSTATUS)
IF (OPENSTATUS > 0) STOP "***CANNOT CREATE SINKFLOW.TXT FILE***"
!CREATE OUTPUT FILE FOR SOIL WATER BALANCE RUNOFF PER BASIN PER INTERVAL
PRINT *, "CREATING OUTPUT FILE INFLOW.TXT.........."
OPEN (UNIT = 16 , FILE = "INFLOW.TXT" , STATUS ="UNKNOWN" , ACTION
.= "READWRITE" , POSITION = "REWIND" , IOSTAT = OPENSTATUS )
IF (OPENSTATUS > 0) STOP "***CANNOT CREATE INFLOW.TXT FILE***"
!CREATE AN OUTPUT FILE FOR STORING LOSSES FOR BASINS
PRINT *, "CREATING OUTPUT FILE LOSSES.TXT.........."
OPEN (UNIT = 17 , FILE = "LOSSES.TXT" , STATUS ="UNKNOWN" , ACTION
.= "READWRITE" , POSITION = "REWIND" , IOSTAT = OPENSTATUS )
IF (OPENSTATUS > 0) STOP "***CANNOT CREATE LOSSES.TXT FILE***"
PRINT *, "ROUTING FLOWS TO SINK.........."
DOX=1,SOURCES,1
DOY=1,STEPS,1
!LANDCELL(X,2) IS THE AREA OF THE POLYGON IN M2
!COMPUTE FLOW FROM RUNOFF, AND CONVERT FROM M/DAYTOM3/S
!COMPUTE LOSSES FROM FLOW TIME AND LOSS COEFFICIENT
! UPDATE FLOW AVAILABLE FOR ROUTING
RUNOFF(X,Y) = ((RUNOFF(X,Y) * LANDCELL(X,2))/86400)
DO M = 1, SINKCELLS,1
IF (SINKCELL(M) .EQ. LANDCELL(X,3)) THEN
INFLOW(M, Y) = RUNOFF(X,Y)*86400 + INFLOW(M, Y)
LOSSES(M,Y)= RUNOFF(X,Y)*86400*(1 - LOSSFACTOR(X))
.+OSSES(M,Y)
190
GOTO 100
END IF
END DO ! M = 1, SINKCELLS,1
100 IF (FLOWTIME(X).GT.0)THEN
RUNOFF(X,Y)=RUNOFF(X,Y)*LOSSFACTOR(X)
ELSE
RUNOFF(X,Y) = RUNOFF(X,Y)*1
END IF
ENDDO !Y=1,STEPS,1
!PRINT*,"WRITING FLOWS TO LANDFLOW.TXT.........."
! WRITE (14,70) LANDCELL(X,1), (LANDFLOW(X,Y), Y = 1, STEPS)
ENDDO !X=1,SOURCES,1
PRINT *, "WRITING TOTAL RUNOFF TO INFLOW.TXT.........."
WRITE (16, 90) (0.0), (SINKCELL(Z), Z = 1, SINKCELLS)
WRITE (17, 90) (0.0), (SINKCELL(Z), Z = 1, SINKCELLS)
DO M = 1, STEPS,1
WRITE (16, 80) (M-1), (INFLOW(Z, M), Z = 1, SINKCELLS)
WRITE (17, 110) (M-1), (LOSSES(Z, M), Z = 1, SINKCELLS)
END DO
OUTSTEPS =(STEPS +INT(MAXFLOWTIME)+INT(MAXWIDTH/2))
ALLOCATE (SINKFLOW(SINKCELLS,OUTSTEPS), STAT = ALLOCATESTATUS)
IF (ALLOCATESTATUS /= 0) STOP "****NOT ENOUGH MEMORY****"
!INITIALIZE SINKFLOW()
DO R = 1, OUTSTEPS,1
DO S = 1, SINKCELLS,1
SINKFLOW(S, R) = 0
END DO ! S = 1, SINKCELLS,1
191
END DO ! R = 1, OUTSTEPS,1
PRINT *, "AGGREGATING FLOWS AT OUTLET.........."
DOX=1,SOURCES,1
DOY=1,STEPS,1
DO Z = 1, SINKCELLS,1
IF (SINKCELL(Z) .EQ. LANDCELL(X,3)) THEN
LTIME=INT(FLOWTIME(X) )
DO M = 1, MAXWIDTH,1
CTIME =Y+LTIME +M-INT(MAXWIDTH/2) - 1
IF (CTIME .GT. 0) THEN
SINKFLOW(Z,CTIME)= SINKFLOW(Z,CTIME) +
.(RUNOFF(X,Y)*COEFFICIENTS(X, M))
ELSEIF (CTIME .LE. 0) THEN
CTIME =Y+LTIME
SINKFLOW(Z,CTIME)= SINKFLOW(Z,CTIME) +
.(RUNOFF(X,Y)*COEFFICIENTS(X, M))
END IF ! IF (COEFFICIENTS(X,M).GT.0)THEN
END DO ! M = 1, MAXWIDTH,1
END IF ! IF (SINKCELL(Z) .EQ. LANDCELL(X,3)) THEN
END DO ! Z = 1, SINKCELLS,1
ENDDO ! Y=1,SOURCES,1
ENDDO ! X=1,SOURCES,1
PRINT *, "WRITING OUTPUT FILE SINKFLOW.TXT.........."
WRITE (15, 90) 0.0, (SINKCELL(Z), Z = 1, SINKCELLS)
DOY=1,OUTSTEPS,1
WRITE (15, 90) (Y-1), (SINKFLOW(Z, Y), Z = 1, SINKCELLS)
END DO
192
80 FORMAT(I5, 1000F15.1)
90 FORMAT(F7.2, 1000F15.5)
110 FORMAT(I5, 1000F20.1)
200 FORMAT(F10.1, X, 1000F12.10)
PRINT *, "PROCESSING COMPLETE"
END
193
Appendix 4: AML Code for Computing Equivalent Muskingum
Parameters from a Field of Velocities and Dispersion Coefficients
/* AML FOR COMPUTING MUSKINGUM RIVER REACH PARAMETERS FROM FIELDS OF
/* FLOW VELOCITIES AND DISPERSION COEFFICIENTS
/* INPUTS INCLUDE GRIDS INCLUDE THE FOLLOWING
/* FLOWDIR: A FLOW DIRECTION GRID
/* STREAMGRID: A STREAM GRID
/* OUTLETGRID: AN OUTLET GRID
/* LINKGRID: A STREAM LINK GRID
/* VELOCITY: A GRID OF LOCAL VELOCITY
/* DISPCOEFF: A GRID OF LOCAL DISPERSION COEFFICIENT
GRID
SETWINDOW FLOWDIR FLOWDIR
SETCELL FLOWDIR
/* COMPUTE RIVER FLOWDIR FROM STREAM GRID AND BASIN FLOWDIR
RIVERFDR = FLOWDIR * STREAMGRID
/* INSERT NODATA CELLS AT STREAM JUNCTIONS BY MULTIPLYING BY OUTLET GRID
LINKFDR = CON(ISNULL(OUTLETGRID), RIVERFDR)
/* COMPUTE LINK FLOW LENGTH
LINKLENGTH = FLOWLENGTH ( LINKFDR )
/* COMPUTE MAXIMUM LINK FLOW LENGTH
LINKLENMAX = ZONALRANGE (LINKGRID, LINKLENGTH)
/* COMPUTE THE INVERSE OF VELOCITY
INVVEL =1/VELOCITY
/* COMPUTE WEIGHTED FLOW LENGTH GRID USING THE INVERSE OF VELOCITY AS
WEIGHT
LINKTIME = FLOWLENGTH ( LINKFDR, INVVEL )
/* COMPUTE THE FLOW TIME FOR EACH STREAM LINK
LINKTIMEMAX = ZONALRANGE (LINKGRID, LINKTIME)
/* COMPUTE THE FLOW VELOCITY FOR EACH STREAM LINK
LINKVELOCITY = LINKLENMAX / LINKTIMEMAX
/* COMPUTE THE DISPERSION WEIGHT, (D/V3)
DISPWEIGHT = DISPCOEF * POW (INVVEL,3)
/* COMPUTE WEIGHTED FLOW LENGTH GRID USING (D/V3) AS WEIGHT
LINKVARIANCE = FLOWLENGTH ( LINKFDR, DISPWEIGHT )
194
/* COMPUTE LINK VARIANCE
LINKVARMAX = ZONALRANGE (LINKGRID, LINKVARIANCE)
/* COMPUTE LINK MUSKINGUM X
LINKDISP = LINKVARMAX * (POW (LINKVELOCITY, 3)) / LINKLENMAX
/* COMPUTE LINK MUSKINGUM X
LINKMUSKX =0.5-(LINKDISP /(LINKVELOCITY * LINKLENMAX))
/* CONDITIONAL STATEMENT TO REMOVE NEGATIVE MUSKINGUM X VALUES
LINKMUSKX1=CON ( LINKMUSKX <0,0,LINKMUSKX )
/* COMPUTE LINK MUSKINGUM N, 86400 SECONDS ASSUMING DAILY ROUTING
INTERVAL
LINKMUSKN =(INT (2 * LINKMUSKX1*LINKTIMEMAX / 86400)+1)
/* RECOMPUTE MUSKINGUM X BASED ON MUSKINGUM N
LINKMUSKX2 = 0.5 - LINKMUSKN * ( 0.5 - LINKMUSKX1)
/* CONDITIONAL STATEMENT TO REMOVE NEGATIVE MUSKINGUM X VALUES
LINKMUSKX3=CON ( LINKMUSKX2<0,0,LINKMUSKX2)
/* COMBINE LINKGRID AND LINKLAG
JOINGRID = COMBINE ( LINKGRID, LINKTIMEMAX, LINKVELOCITY, LINKMUSKX3)
&RETURN
195
Appendix 5: 10 day flows in the Nile basin, 1953-57
OBSERVED FLOWS AT MONGALLA, 1953-57,
(REPORTED IN HURST AND PHILLIP, 1931-)
PERIOD FLOWS PERIOD FLOWS PERIOD FLOWS PERIOD FLOWS PERIOD FLOWS
0 65.7 360 49.9 720 53.6 1080 53.9 1440 57.9
10 64.3 370 48.9 730 53 1090 52.9 1450 57.4
20 62.7 380 47.4 740 52.3 1100 50.5 1460 56.8
30 61.8 390 46.1 750 52.3 1110 49.9 1470 56.7
40 59.4 400 45.7 760 51.1 1120 49.1 1480 55.9
50 56.7 410 45.5 770 50.4 1130 48.6 1490 55.4
60 55.9 420 45 780 49.4 1140 48.6 1500 55.7
70 54.6 430 44.5 790 48.9 1150 47.9 1510 56.6
80 53.2 440 43.9 800 48.1 1160 47.9 1520 59.5
90 52.6 450 51.2 810 47.2 1170 49.7 1530 63.6
100 52.4 460 54.4 820 52.7 1180 54.1 1540 62.6
110 56.2 470 46.8 830 49.9 1190 60.2 1550 67.6
120 62.6 480 64.9 840 58.3 1200 59.6 1560 68.9
130 59.8 490 58 850 58.9 1210 62.5 1570 67.4
140 62.8 500 61.6 860 51.8 1220 66.4 1580 92.1
150 63.4 510 62.9 870 50.2 1230 67.4 1590 98.5
160 61.9 520 60.5 880 49.9 1240 60.4 1600 83.3
170 69.6 530 56.8 890 53.1 1250 56.9 1610 87.4
180 63.8 540 63.4 900 56.8 1260 62.2 1620 74.8
190 78.7 550 66.4 910 52.7 1270 58.4 1630 68.6
200 64.5 560 66 920 54.4 1280 66.8 1640 66.5
210 65.6 570 75.9 930 59.3 1290 66.2 1650 82.1
220 84.4 580 82.9 940 79.4 1300 79.8 1660 87
230 76.4 590 101 950 74.2 1310 94.8 1670 72.1
240 60.3 600 107 960 87.1 1320 105 1680 72.9
250 64.1 610 112 970 101 1330 109 1690 71.2
260 63.2 620 82.1 980 98.4 1340 91.8 1700 64.8
270 55.3 630 78.2 990 116 1350 102 1710 65.9
280 58.5 640 70 1000 98 1360 98.7 1720 63.8
290 70.2 650 68.3 1010 83.5 1370 91.6 1730 71.3
300 67.5 660 68.3 1020 82 1380 78.3 1740 65.5
310 64.4 670 62 1030 86.6 1390 68.4 1750 67
320 58.6 680 56.9 1040 68.1 1400 63.1 1760 63.1
330 55.6 690 55.6 1050 59.9 1410 61.8 1770 62.6
340 53.6 700 54.5 1060 57 1420 60.2 1780 61.7
350 51.1 710 54.1 1070 55.7 1430 58 1790 60.7
196
OBSERVED FLOWS AT LAKE NYONG, 1953-57
(REPORTED IN HURST AND PHILLIP, 1931-)
PERIOD FLOWS PERIOD FLOWS PERIOD FLOWS PERIOD FLOWS PERIOD FLOWS
0 49.1 360 44.8 720 45.4 1080 48.6 1440 55.1
10 46.9 370 44.3 730 44.8 1090 48.2 1450 53.1
20 44.8 380 43.4 740 44.2 1100 47.9 1460 50.6
30 45.2 390 42.4 750 43.7 1110 47.3 1470 48.2
40 45.2 400 41.8 760 43.4 1120 46.5 1480 47
50 45.2 410 41.5 770 43.2 1130 45.4 1490 46
60 45.6 420 41.1 780 42.9 1140 44.3 1500 45.2
70 45.4 430 41 790 42.9 1150 43.2 1510 45.3
80 44.4 440 40.9 800 42.9 1160 42 1520 45.6
90 42.5 450 40.6 810 42.8 1170 40.8 1530 46
100 42.2 460 40 820 42.8 1180 41.2 1540 46
110 42.9 470 39.2 830 42.8 1190 42.5 1550 46
120 44.2 480 38.5 840 42.7 1200 43.7 1560 46
130 43.5 490 38.6 850 42.7 1210 43.8 1570 46.8
140 42.7 500 39.1 860 42.7 1220 43.6 1580 47.8
150 43.6 510 39.7 870 42.7 1230 43.5 1590 48.7
160 43.9 520 40.3 880 42.7 1240 43.7 1600 50.1
170 44.2 530 41.3 890 43 1250 43.9 1610 51.6
180 44.9 540 42.2 900 44 1260 44.1 1620 52.7
190 45.5 550 43.1 910 44.9 1270 44.5 1630 52
200 46.3 560 44.2 920 45.5 1280 45.1 1640 50.9
210 47.3 570 45.1 930 45.7 1290 45.8 1650 50.1
220 48.2 580 46.1 940 45.8 1300 47 1660 51.3
230 48.8 590 47 950 46.3 1310 48.7 1670 53
240 48 600 48 960 47.3 1320 50.3 1680 54.4
250 47.4 610 49.2 970 48.4 1330 51.1 1690 53.9
260 47.1 620 50.5 980 49.4 1340 51.6 1700 52.9
270 46.6 630 51.6 990 50.5 1350 52.2 1710 52
280 46.5 640 51.2 1000 50.6 1360 52.9 1720 51.5
290 46.5 650 50.5 1010 49.6 1370 53.7 1730 51
300 46.2 660 49.7 1020 48.7 1380 54.4 1740 50.5
310 46 670 48.9 1030 48.6 1390 54.5 1750 49.9
320 45.6 680 48.6 1040 48.8 1400 54.4 1760 49.2
330 45.5 690 48.3 1050 48.9 1410 54.4 1770 48.6
340 45.2 700 47.5 1060 48.9 1420 54.7 1780 48.8
350 45 710 46.5 1070 48.7 1430 55 1790 49.2
197
REFERENCES
Barkstrom, B., E. Harrison, R. Lee, (1990), Earth Radiation Budget Experiment,
Preliminary Seasonal Results, Transactions - American Geophysical
Union, vol. 71, pp. 279, 299, 304-305
Bartholomew (John) and Sons, ltd., (1975), The Times Altas of the World,
Quadrangle/The New York Times Book Company, NY, NY
Birkett,C.M. and I.M.Mason, (1995), A new Global Lakes Database for a semote
sensing programme studying climatically sensitive large lakes, Journal of
Great Lakes Research, Vol. 21, no. 3, pp. 307-318
Bonan, G.B., (1996), The NCAR land surface model (LSM version 1.0) coupled
to the NCAR Community Climate Model, NCAR Technical Note
NCAR/TN-429+STR, National Center for Atmospheric Research,
Boulder, Colorado
Branstetter, M. L., and J. S. Famiglietti, (1998), An Investigation of the Effects of
Continental Runoff on Climate Dynamics Using a Parallel Earth System
Model, Proceedings of the American Geophysical Union Spring Meeting,
Boston, MA
Carlston, C.W., (1969), Downstream variations in the hydraulic geometry of
streams: Special emphasis on mean velocity, American Journal of Science,
vol. 267, pp. 499-509
Chow, V.T., D.R. Maidment, L.W. Mays, (1988), Applied Hydrology, McGraw-
Hill, Inc., NY, NY
Clark, C.O., (1945), Storage and the Unit Hydrograph, Transactions, America
Society of Civil Engineers, vol. 110, pp. 1419-1488
Coe, M., (1994) Simulating Continental Surface Waters: An Application to
Holocene Northern Africa, Journal of Climate, v10, pp. 1680 -1689
Cunge, J.A., (1969), On the subject of a flood propagation computation method
(Muskingum Method), Journal of Hydraulic Research, 7, no. 2, pp. 205-
230
Dooge, J.C., (1973), Linear Theory of Hydrologic Systems, USDA Technical
Bulletin no. 1468, USDA, Washington, DC
198
Easterling, D. R., T. Peterson, and R. K. Thomas, (1996), On the development
and use of homogenized climate data sets, Journal of Climate, 9, 1429-
1434, ?http://www.ncdc.noaa.gov/wdcamet.html#GPCP?
Environmental Systems Research Institute Inc. (ESRI), (1992), ArcWorld 1:3M,
User?s Guide & Data Reference, ESRI, Redlands, CA
Environmental Systems Research Institute Inc. (ESRI), (1993), Data Dictionary:
The Digital Chart of the World, ESRI, Redlands, CA
Food and Agriculture Organization (FAO) of the United Nations, (1970-78), Soil
map of the world, scale 1:5,000,000, volumes I- X: United Nations
Educational, Scientific, and Cultural Organization, Paris,
?http://daac.gsfc.nasa.gov/CAMPAIGN_DOCS/FTP_SITE/INT_DIS/read
mes/soils.html?
Fischer, H.B., E.J. List, R.C.Y. Koh, J. Imberger and N.H. Brooks, (1979),
Mixing in Inland and Coastal Waters, Academic Press, New York
Fischer, H.B., (1967) The Mechanics of Dispersion in Natural Streams, Journal of
the Hydraulic Division, ASCE, v.93, 187 - 216
Gupta, V.K., I. Rodriguez-Iturbe, and E.F. Wood. (eds.), (1986), Scale Problem in
Hydrology : runoff generation and basin response, Dordrecht, Holland
Gupta, V.K., O.J. Mesa, D.R. Dawdy, (1994), Multiscaling theory of flood peaks:
Region quantile analysis, Water Resources research, vol. 30, no. 12, pp.
3405-3421
Hagemann, S. and L. Dumenil, (1998), A parameterization of the lateral
waterflow for the global scale, Climate Dynamics 14: 17 ? 31
Hellweger, F., and D.R.Maidment (1996), A GIS Preprocessor for Lumped
Parameter Hydrologic Modeling Programs, CRWR Online Report 97-8,
The University of Texas at Austin
Hortal, M., and A.J. Simmons, 1991: Use of reduced Gaussian grids in spectral
models, Monthly Weather Review, 119, pp. 1057-1074.
Huffman, G.A., R.F. Adler, P.A. Arkin, A. Chang, R. Ferraro, A. Gruber, J.
Janowiak, R.J. Joyce, A. McNab, B. Rudolf, U. Schneider, and P. Xie,
(1997), The Global Precipitation Climatology Project (GPCP) Combined
199
Precipitation Data Set, Bulletin of American Meteorological Society, 78,
pp. 5-20, ?ftp://ftp.ncdc.noaa.gov/pub/data/gpcp/version1/?
Hurst, H.E. and Phillips, (1931-), The Nile Basin, Volumes I ? V, Nile Control
Department
Kiehl, J.K., J.J. Hack, G.B. Bonan, B.A. Boville, B.P. Briegleb, D.L. Willianson,
P.J.Rasch, (1996), Description of the NCAR Community Climate Model
(CCM3), Tn-420+STR, NCAR, Boulder, Colorado
Kite, G.W., A. Dalton, K. Dion, (1994), Simulation of streamflow in a macroscale
watershed using general circulation model data, Water Resources
Research, vol. 30, No. 5, pp. 1547-1559
Korzoun, V.I., A.A. Sokolov, M.I. Budyko, K.P. Voskresensky, G.P. Kalinin,
E.S. Konoplyantsev and M.I. Lvovich (edited by), (1977), Atlas of World
Water Balance, The UNESCO Press, Paris
Kull, D.W., Feldmean, (1998), Evolution of Clark?s Unit Graph Method to
Spatially Distributed Runoff, Journal of Hydrologic Engineering, v3, no.1,
pp. 9-19
Lasdon, L. and A.D. Waren, (1978), Generalized Reduced Gradient Software for
Linearly and Nonlinearly Constrained Problems, in Design and
Implementation of Optimization Software, H. Greenberg, ed., Sijthoff and
Noordhoff, Pubs., 1978, pp. 363-397
Legates, D.R., and C.J. Willmott (1990), "Mean Seasonal and Spatial Variability
in Gauge-Corrected Global Precipitation", International Journal of
Climatology, v.10, pp. 111-127
?http://www.scd.ucar.edu/dss/datasets/ds236.0.html?
Leopold, L.B., (1953), Downstream Change of Velocity in rivers, American
Journal of Science, vol.251, pp. 606-624
Liston, G.E., Y.C. Sud, E.F. Wood (1994), Evaluating GCM land surface
hydrology parameterizations by computing river discharges using a runoff
routing model: Application to the Mississippi Basin, Journal of Applied
Meteorology, v33, n3, pp.394 - 406
Loaiciga, H.A., (1997), Runoff scaling in large rivers of the world, Professional
Geographer, vol.49, n3, pp.356-365
200
Lohmann, D., E. Raschke, B. Nijssen, D. P. Lettenmaier, (1998), Regional Scale
Hydrology: 1. Formulation of VIC-2L Model Coupled to a Routing
Model, Hydrological Sciences Journal, vol.43, n. 1, pp.131-141
Lohmann, D., R. Nolte-Holube, E. Raschke, (1996), A Large-Scale Horizontal
Routing Model to be Coupled to Land Surface Parameterization Schemes,
Tellus, vol. 48A, pp.708-721
Kazadi, S., F. Kaoru, (1996), Interannual and long-term climate variability over
the Zaire River Basin during the last 30 years, Journal of Geophysical
Research, v.101, no. D16, 21,351 - 360
Maidment,D.R.(ed.) (1993), Handbook of Hydrology, McGraw-Hill, Inc., NY,
NY
Maidment, D.R., F. Olivera, A.Calver, A. Eatherall, W.Fraczek, (1996), Unit
Hydrograph Derived from a Spatially Distributed Velocity Field,
Hydrological Processes, vol. 10, 831-844
Meyer, O.H., (1941), Simplified Flood Routing, Civil Engineering, v.11, no. 5,
pp. 306 - 307
McCuen, R.H. (1998), Hydrologic Analysis and Design, Prentice-Hall Inc., Upper
Saddle River, New Jersey
Miller, J.R., G.L.Russell, G.Caliri (1994) Continental-Scale River Flow in
Climate Models, Journal of Climate, v7, pp. 914 - 928
Montes, S., (1998), Hydraulics of Open Channel Flow, ASCE Press, Reston, VA
Naden, P.S., (1993), A routing model for continental-scale hydrology, Macroscale
modeling of the the Hydrosphere, Proceedings of Yokahama Symposium,
July 1993, International Association of Scientific Hydrology (IASH)
Publ., 214, 67-79
Naden, P.S., (1992), Spatial variability in flood estimation for large catchments:
the exploitation of channel network structure, Hydrological Science
Journal, v37, pp.53-71
Nash, J.E., (1959), A note on the Muskingum Flood-Routing Method, Journal of
Geophysical Research, vol. 64, pp. 1053 ? 1056
201
Olivera, F., (1996), Spatially Distributed Modeling of Storm Runoff and Non-
Point Source Pollution Using Geographic Information Systems, Doctoral
Dissertation, Submitted to the Graduate School, The University of Texas
at Austin
Olivera, F., and D.R.Maidment (1999), GIS Tools for HMS Modeling Support,
Proceedings of the 19
th
ESRI Users Conference-CDROM, San Diego,
California
Olivera, F., and D.R.Maidment (1999), Geographic Information Systems (GIS)
Based Spatially Distributed Model for Runoff Routing, Water Resources
Research, vol. 35, no. 4, pp. 1155-1164
Olivera, F., J.S.Famiglietti and K.O.Asante (2000), Global-Scale Flow Routing
Using a Source-To-Sink Algorithm, Submitted to Water Resources
Research
Perry, G.D., P.B. Duffy, and N.L. Miller, (1996), An extended data set of river
discharges for validation of general circulation models, Journal of
Geophysical Research, v.101, no. D16, pp. 21,339-349
?ftp://kosmos.agu.org/? in the directory ?/apend/96Jd00932?
Peters, J.C., (1998), HEC-HMS, Hydrologic Modeling System, US Army Corps
of Engineers Hydrologic Engineering Center, Davis, CA
Pilgrim, D.H., (1981), Effects of catchment size on runoff relationships, J.
Hydrology, 58: 205-221
Pilgrim, D.H., (1977), Isochrones of Travel Time and Distribution of Flood
Storage from a Tracer Study on a Small Watershed, Water Resources
Research, vol. 13, no. 3, pp. 587-595
Ponce, V.M., R.M.Li, D.B.Simons, (1978), Applicability of Kinematic and
Diffusion Models, J.Hydraulic Engineering, vol.104, no.hy3, pp. 353-
360,1978
Reed, S.M., J. Patoux, D.R. Maidment (1997), CRWR Online Report 97-1:
Spatial Water Balance of Texas, CRWR, The University of Texas at
Austin
Reed, S.M., D.R. Maidment (1997), CRWR Online Report 95-3: A GIS
Procedure for Merging NEXRAD Precipitation Data and Digital Elevation
202
Models to Determine Rainfall-Runoff Modeling Parameters, CRWR, The
University of Texas at Austin
Revenga, C., S. Murray, J. Abramovitz, A. Hammond, (1998), Watersheds of the
World, A joint publication of the Water Resources Institute and
WorldWatch Institute, Washington, D.C.
Rodriguez-Iturbe, I. and J.B. Valdes, (1979), The geomorphological structure of
hydrologic response, Water Resources Research, vol. 15, no. 6, pp. 1409-
1420
Sherman, L.K. (1932), Streamflow from rainfall by unit-graph method,
Engineering News record, Vol. 108, 501-505
Snyder, F.F. (1938), Synthetic Unit-Graphs, Trans. Am. Geophysical Union, Vol.
19, 447-454
Sausen, R., S.Schubert, L.Dumenil, (1994), A model of river runoff for use in
coupled atmosphere-ocean models, Journal of Hydrology, 155 (1994), pp.
337 - 352
Seo, I.W., T.S. Cheong, (1998), Predicting Longitudinal Dispersion Coefficient in
Natural Streams, Journal of Hydraulic engineering, vol. 124, no.1, pp. 25 -
32
Shahin, M., (1985), Hydrology of the Nile Basin, Developments in Water Science
No.21, Elsevier, Amsterdam, Netherlands
Simmons, A.J., and D.M. Burridge, (1981), An energy and angular-momentum
conserving vertical finite difference scheme and hybrid vertical
coordinates, Monthly Weather Review, 109, pp.758-766.
Sposito, G, (1998), Scale Dependence and Scale Invariance in Hydrology,
Cambridge University Press, Cambridge, UK
Storm, B. and J. Kjelds, (1999), Integrated River Basin Management Modeling
Water Use and Water Quality Simulation, Proceedings of the ASCE Water
Resource Conference, Seattle, WA
USGS, (1996), GTOPO30 Documentation, Sioux Falls, South Dakota
USGS, (1997), HYDRO1K Documentation, Sioux Falls, South Dakota
203
V?r?smarty, C.J., B.Fekete, and BA Tucker (1996), River Discharge Database,
Version1.0 (RIvDIS v1.0) Volumes 0-6, A contribution to IHP-V Theme
1, Technical Documents in Hydrology Series, UNESCO, Paris
?http://www.rivdis.sr.unh.edu/?
V?r?smarty, C.J., B.Moore III, A.L. Grace, M.P. Gildea, M.Melillo, B.J.
Peterson, E.B. Rastetter, P.A. Steudler, (1989), Continental Scale Models
of Water Balance and Fluvial Transport: An Application to South
America, Global Biogeochemical Cycles, v3, no.3, pp. 241 - 265
Vose,R.S.,R.L.Schmoyer,P.M.Steurer,T.C.Peterson,R.Heim,T.R.Karl,and
J. Eischeid, (1992), The Global Historical Climatology Network: long-
term monthly temperature, precipitation, sea level pressure, and station
pressure data. ORNL/CDIAC-53, NDP-041, Carbon Dioxide Information
Analysis Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee
Zoch, R,T., (1934), On the relation between rainfall and stream flow, Monthly
Weather Review, vol. 62
204
Vita
Kwabena Oduro Asante was born in Kintampo, Ghana, on November 23,
1971, the son of Beatrice Oduro Asante and Kofi Oduro Asante. After completing
his work at Starehe Boys Centre and School, Nairobi, Kenya, in 1989, he studied
computer programming at the Kestrel College, Nairobi. He entered the University
of Nairobi in August of 1990 and received the degree of Bachelor of Science in
Civil Engineering in December, 1994. He served as an Engineering Trainee at the
Water Resources Research Institute in Accra, Ghana, before entering the Graduate
School of The University of Texas in August, 1995. He was awarded the degree
of Master of Science in Civil Engineering in December, 1997. During his tenure
at the university, he served as a teaching assistant in the Department of Civil
Engineering and as a graduate research assistant at the Center for Research in
Water Resources.
Permanent address: P.O.Box 3188, Accra, Ghana
This dissertation was typed by the author.