DIFFUSION W AVE MODELING OF DISTRIBUTED CATCHMENT DYNAMICS

By Stefano Orlandini 1 and Renzo Rosso2

ABSTRACT: A diffusion wave model of distributed catchment dynamics is presented. The effects of catchment
topography and river network structure on storm-flow response are incorporated by routing surface runoff in
cascade throughout a digital elevation model (OEM) based conceptual transport network, where the Muskingum-
Cunge scheme with variable parameters is used to describe surface runoff dynamics. Dynamic scaling of hy-
draulic geometry is also incorporated in the model by using the "at-a-station" and "downstream" relationships
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

by Leopold and Maddock. Numerical experiments indicate that the model is more than 98% mass conservative
for possible slope and roughness configurations, which may occur for hillslopes in a natural catchment. Fluc-
tuations in the simulated discharge may occur in response to discontinuities in rainfall excess representation if
Courant number Cu during the simulation exceeds a threshold of about 3. Catchment scale simulations with
different temporal resolution show that the model response is independent of structural parameters (model
consistency). Also, the overall accuracy is preserved for computationally inexpensive space-time discretizations
(for which Cu > 3) because fluctuations that may occur at the local scale are dampened when propagating
downstream. Comparison of model results with observed outlet hydrographs of the Rio Missiaga experimental
catchment (Eastern Italian Alps) show this approach to be capable of describing both overland and channel
phases of surface runoff in mountainous catchments.

INTRODUCTION scale to basin scales, experience spreading mechanisms from

molecular to geomorphological diffusion in which the full de-
Recent advances in remote sensing, geographic information velopment of any transition overwhelms the previous ones.
systems, and computer technology make the use of distributed This indicates that models based on accurate description of the
hydrologic models an attractive approach to flow simulation geometry and topology of the river network and a simplified
and prediction. The linkage of a distributed hydrologic model dynamics can capture the features of travel time distributions
with the spatial data handling capabilities of digital elevation in a broad range of dispersivities.
models (OEMs) and digital terrain models (DTMs) offers ad- In this paper a simple and computationally inexpensive
vantages associated with the full utilization of spatially dis- model for distributed catchment dynamics is developed. The
tributed data to analyze hydrologic processes. The major areas control of downslope water movement by catchment topog-
of application of distributed models are in forecasting the ef- raphy and channel network topology is treated by processing
fects of land-use change, the effects of spatially variable inputs each of the DEM cells, from upland areas to the basin outlet,
and outputs, the movement of pollutants and sediments, and through a conceptual transport network extracted from DEM
the hydrologic response of ungauged catchments (Beven and data. Each elemental network channel is assumed to have bed
O'Connell 1982). slope and length that depend on location within the DEM-
Surface runoff constitutes an important component of the extracted transport network, and rectangular cross section
hydrologic response of a catchment to climate and land use. whose width varies dynamically according to the scaling prop-
A pioneering approach to the physical modeling of surface erties of stream geometry as described by the "at-a-station"
runoff was initiated by the simplified model of a catchment and "downstream" relationships first introduced by Leopold
idealized by an open-book geometry (Wooding 1966). Also, and Maddock (1953). Distinction between overland and chan-
numerical schemes to solve the governing equations of the nel flow is based on the "constant critical support area" con-
physical processes have been given increasing importance, cept [see Montgomery and Foufoula-Georgiou (1993)]. The
both for overland flows (Kibler and Woolhiser 1972; Ross et use of grid-based OEMs instead of vector elevation data [see
al. 1979; Brandford and Meadows 1990) and for channel flows Moore and Grayson (1991)] is motivated by a series of inter-
(Szymkiewicz 1991). Despite the abundance of modeling acting factors. On one hand, grid-based digital models are sim-
schemes, the computational resources needed to run detailed ple to use and elaborate and are still the more widespread
models may be excessive for large-scale catchment simula- support available for applications in operational hydrology. On
tions; hence there is a need for simpler conceptual techniques, the other hand, for large-scale catchment simulations, the dy-
which embody the essential concepts of component subpro- namic scaling of channel geometry (see section titled "Chan-
cesses while also allowing space-time variability to be consid- nel Network Geometry") on a quasi-two-dimensional network
ered. extracted from grid-based data allows a physically realistic de-
Using the geomorphological-dispersion approach Rinaldo et scription of surface flow phenomena without incurring more
al. (1991) found that in channel networks the dispersive char- complicated (and not always physically realistic) two-dimen-
acter of the transport processes is determined by geometry and sional sheet flow representations.
topology rather than the dispersive processes dominant at the A Muskingum-Cunge method with variable parameters is
mesoscale and microscale. Transport processes, from micro- applied to control the amount of numerical diffusion in such
a way that it matches the diffusion of the physical problem.
I Dr., OllAR, Politecnico di Milano, Piazza Leonardo da Vinci, 32,
Milano 20133, Italy. (e-mail: stefano@idral.iar.polimi.it) The reliability and robustness of the routing model in terms
'Prof., OllAR, Politecnico di Milano, Piazza Leonardo da Vinci, 32, of mass conservation, accuracy, and consistency are investi-
Milano 20133, Italy. (e-mail: rr@idral.iar.polimi.it) gated through numerical experiments for a wide range of pos-
Note. Discussion open until December I, 1996. To extend the closing sible configurations of basin hillslopes (see section titled "Nu-
date one month, a written request must be filed with the ASCE Manager merical Experiments at Hillslope Scale"). Model application
of Journals. The manuscript for this paper was submitted for review and
possible publication on November I, 1995. This paper is part of the
to the 4.35-km2 Rio Missiaga catchment (Eastern Italian Alps)
Jour""l of Hydrowgic Engineering, Vol. I, No.3, July, 1996. @ASCE, shows that, even for computationally inexpensive space-time
ISSN 1084-0699/96/0003-0103-0113/$4.00 + $.50 per page. Paper No. discretizations, mass conservation, accuracy, and consistency
11936. are significantly preserved (see section titled "Catchment
JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996/103

J. Hydrol. Eng., 1996, 1(3): 103-113

 Scale Modeling"). Finally, as an example of model application of that point and the exponent b' is a characteristic of the
for investigative purposes, Eulerian and Lagrangian analysis channel network as a whole. When considering discharge with
of a Rio Missiaga catchment flood event is presented. Future the same frequency, Q" at different points moving downstream
work will consider large-scale catchment simulations to further in the basin the corresponding channel width B r is given by
test the capability of the model to reproduce the distributed
flood dynamics and investigate the space-time variability of Br = a"Qt; (2)
flow characteristics throughout the basin. where a" and b" = characteristics of the river network. Al-
though this is not necessarily equivalent to considering dis-
MODEL DESCRIPTION charge resulting from the same event at different points, it can
Local contributions to infiltration excess runoff produced by be assumed as a reasonable indicator of what may occur in
that situation (Bras 1990).
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

the time compression approximation (TCA) model described

in Orlandini et al. (in press, 1996) are routed onto a DEM- At-a-station and downstream relationships, (1) and (2), re-
based conceptual transport network via a Muskingum-Cunge spectively, can be coupled into a single equation to represent
scheme with variable parameters. Surface runoff is not allowed the dynamics of flow width at any point of the channel net-
to infiltrate in downslope cells even when the soil profile of work in response to a rainfall excess forcing event. From (1)
these cells is not saturated. The runoff-runon problem is likely one obtains, for a given point in the channel network
to be important where broad sheet flow occurs but relatively
unimportant when the flow concentrates in defined rills, chan- (3)
nels, and streams, as assumed in this model. There are two
main aspects of the routing model presented in this work. The where Br and Qr = values of flow width and discharge at that
first aspect concerns the description of a watershed in terms point as given by (2). Eq. (2) yields
of hydrologic features; the second concerns the representation
of the flood wave phenomena that significantly characterize Br (Qr)b"
the catchment dynamics. These aspects are described in the B ro = Qro (4)
following two sections, respectively.
where B and Qro = values of the flow width and discharge at
IO

the basin outlet, for the considered frequency. By coupling (3)

Channel Network Geometry and (4), one obtains
The drainage network for a given catchment is automati-
cally extracted from the DEM. Each grid cell is characterized B =B
ro (Qr )b"(Q)b' (5)
by a maximum-slope pointer, and the network links are or- Qro Qr
ganized into a stream-ordering system extracted from the DEM As discharges Qr in the various reaches corresponding to an
data (Band 1986). Ordering is defined by assigning to each outflow Qro are not known, estimation from a geomorpholog-
cell two numbers: the cell order and a link number. The first ical relation is necessary. Leopold et al. (1964) report a rela-
is the sum of the orders of the neighboring upslope cells from tionship of wide applicability as
which the cell can receive water; source cells are assigned
order 1. The link number identifies the first cell of each link (6)
in the network (Orlandini 1995). Distinction between overland
where Qb = bank-full discharge; and A = drainage area. The
and channel flow is based on the "threshold drainage area
exponent in (6) may be closer to unity for small watersheds,
concept" as described in Montgomery and Foufoula-Georgiou
especially in humid climates. To economize on the number of
(1993). This is shown in Fig. 1 with reference to the Rio Mis-
model parameters, the exponent in (6) is assumed equal to
siaga catchment application reported in a later section. A single
unity and the value of Qr at a point of the river network is
formulation is used to describe the hydraulic geometry of the
expressed using the formula
channel network for both overland and channel flow phenom-
ena, hut different parameter values can be taken to characterize (7)
cells in which overland or channel flow is assumed to occur.
Runoff over hillslopes or agricultural watersheds initially where qr = corresponding reference value for rainfall excess;
starts as sheet flow, then it concentrates into a series of small and A = upstream drainage area. Future adjustments of this
channels. The flow concentrations are due to either topo- assumption may be easily incorporated in the context of the
graphic irregularities or differences in soil erodibility. As run- model structure. From (7) the following is obtained:
off continues, erosion progresses and the channels deepen and
Qr A
widen as a function of slope steepness, runoff characteristics, -=- (8)
and soil erodibility. Such erosion-formed microchannels are
called rills or rivulets (Emmett 1978; Li et al. 1980). To min- where A o = catchment area, which is substituted for QrlQro in
imize the computational effort and economize on the number (5) to obtain
of model parameters, in the model presented here the rill for-
mations are lumped at the DEM elemental scale into a single -b' (A )b"-b' b'
conceptual wide channel with bed slope and geometry de- B = BroQro A
o
Q (9)
pending on DEM cell location and flow characteristics.
The elemental channels of the network are taken as wide This can be written as
rectangular channels with width depending on upstream drain- (10)
age area and flow discharge. According to Leopold and Mad-
dock (1953), if one considers discharge of various frequencies where ~ = constant depending on location of the point in the
Q at a point of the river network, the flow width B scales with river network and is given by
Qas
(lA = BroQ;"b'(AIAol"-b' (11)
B =a'Qb' (1)
The hydraulic geometry of the channels in the river network
where the proportionality constant a' depends on the location is lumped into a single parameter, that is, BroQ;"b'. During a
104/ JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996

J. Hydrol. Eng., 1996, 1(3): 103-113

 Q1:: = c1Q1+ 1 + c2Q1 + C3 Q1+1 + C4qt~1 (13)

where Q1::
= discharge at network link point (i + 1)t:..s and
time (j + 1)t:..t; and q{~ll = lateral inflow rate at the (i + l)th
space interval and (j + l)th time interval (see Fig. 2). The
routing coefficients, C" C2 , C3 , and C4 , are given by
C _ ck(t:..tlt:..s) - 2X
(14)
1 - 2(1 - X) + ck(t:..tlt:..s)
ck(t:..tlt:..s) + 2X
(15)
C
2
= 2(1 - X) + ck(f:..tlt:..s)
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

2(1 - X) - ck(t:..tlt:..s)
(16)
C
3
= 2(1 - X) + ck(t:..tlt:..s)

2c t:..t
C4 = - - - _ : : .k . . - _ - - (17)
2(1 - X) + ck(t:..tlt:..s)

where Ck = kinematic wave celerity; and X = weighting factor

used for discretizing the temporal derivative of the kinematic
flow equation
aQ aQ
FIG. 1. 50 x 50 m 2 Resolution OEM of Rio Missiaga Experi- - + Ck - = CkqL (18)
mental Catchment Showing Cells in Which Overland Flow at as
(White Cells) and Channel Flow (Dark Cells) Occurs
X is used to match the numerical diffusion coefficient of the
scheme
flood the active channel width scales dynamically with dis-
charge with the same exponent of the at-a-station relationship, Dn = c,t:..s(1I2 - X) (19)
and with the upstream basin area with an exponent equal to
the difference of those of downstream and at-a-station rela- and the hydraulic diffusivity D h in the convection-diffusion
tionships. If one lets Qro be equal to 1, B ro is the only structural flow equation
parameter of the channel network. aQ aQ a2 Q
Eq. (10) can be applied to represent the network dynamics - + Ck - = Dh -2 + CkqL (20)
at as as
in response to a rainfall excess event. For b ' = 0 hydraulic
geometry varies in space but not in time, while for b ' > 0 it If one incorporates (IO) into the Manning-Gauckler-Strick-
varies both in space and time. This structure is used to describe ler friction equation, (see Appendix I), one obtains
both overland and channel flow. For all those catchment cells
Q= 0il-2/(3+2b')k~/(3+2b')S}I[2(3+2b')]05/(3+2b') (21)
with upstream drainage area A less than a given threshold A *,
it is assumed that only overland flow occurs, whereas channel
where k s = apparent Gauckler-Strickler roughness coefficient;
flow is assumed to occur in cells for which A exceeds A *. This
Sf = friction slope; and 0 = flow area (see Appendix I). For
simple criterion for catchment cell classification is represented
b ' = 0 one has the case of channel with fixed width 0il. The
in Fig. 1 for the Rio Missiaga catchment, assuming A * = 0.50 kinematic flood wave celerity Ck = dQ/dfl is expressed by
km 2 • Different values of b', b", and B ro (Qro = 1) may be as-
sumed to characterize overland and channel flow cells. Ck = 51(3 + 2b')Q/O (22)
Diffusion Wave Formulation that expands into the expression
The same routing scheme is applied for both overland and Ck = 5/(3 + 2b')0il-215k'J5S~1l0Q2I5(1-h') (23)
channel flow, but different distributions of the Gauckler-Strick-
ler roughness coefficient may be considered to account for the in which one has assumed Sf = So, where So = sin ~ is the
different processes that characterize overland and channel channel bed slope; and ~ = channel bed inclination angle (see
flow. The routing model developed in this study is based on Appendix I). The effect of the dynamic scaling of channel
the Muskingum-Cunge method (Cunge 1969). This method is width with flow on the flood wave celerity is represented in
simple to formulate and use, and it allows physical damping Fig. 3 for a set of parameters taken as representative of channel
of the flood wave to be artificially described through the nu- flow for the Rio Missiaga catchment, by varying b' in (23).
merical diffusion in the scheme (due to truncation errors). The From the friction equation (21), one obtains the following
model routes surface runoff downstream, link by link, from
the uppermost cell in the basin to the outlet. A given grid cell
will receive water from its upslope neighbors and discharge to Q!+1 ai+ 1
I 1+1
its downslope neighbor according to the network ordering sys- G+ 1)t.t
tem described earlier. At any catchment cell the lateral inflow 1+1
t.t qLl+1
rate qL (L 2 r- 1) is given by
j t.t j
qL = q!:ut:..ylt:..s (12) Q!I
t.s Q i+1
1
where q (Lr- ) = local contribution to infiltration excess run-
off; t:..x and t:..y = DEM-cell sizes in the direction of the x and
y horizontal coordinates; and t:..s = channel length within the o
cell. o j t.s (i+ 1 )t.s s
Inflow hydrographs and lateral inflows qL are routed onto FIG. 2. Space-Time Computational Grid of Musklngum-Cunge
each individual channel via the routing equation Method

JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996/105

J. Hydrol. Eng., 1996, 1(3): 103-113

 3-.--------------, Simulation of backwater effects in catchment dynamics can be
b·-O.OO- used to allow flow over "digital dams" (due to errors in DEM
b·-0.10-
data); for example, in the model presented by Julien et al.
....
i ..
b·-0.20 ---
b·=0.30····....
b'=0.40 _.• 00
(1995), surface water accumulates behind the barrier until the
.§.2 b·=0.50···· depth exceeds the barrier height and then spills over. In the
.. -..---- model presented here, "digital pits" are filled in a prepro-
.- ---- . cessing step, before extracting the channel network from the
catchment DEM, so that a certain degree of accuracy in the
;fii::j;;;'?;::-:.:::.:::::.:::' catchment topography representation is forfeited in exchange
for simplicity and robustness in the flow dynamics description.
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

Numerical Experiments at Hillslope Scale

Numerical experiments over hypothetical hillslopes are per-
o 2 4 6 8 10 fonned in order to investigate mass conservation, accuracy,
discharge. Q (m~ .-1) and consistency of the presented model. Analysis of the nu-
FIG. 3. Effect of Dynamic Scaling of Channel Width on Flood merical scheme is carried out for the case b' = 0: this simplifies
=
Wave Celerity Ck [E:Jj (23)], for 8 ro 7 m (Oro 1 m 3 S-I), A A o , = = notation without too much loss of generality. In fact, numerical
= =
b" 0.50, k s 5 m' s-" and So 0.10 = experience shows that the numerical behavior of the model is
not very sensitive to b ' and that hi = 0 represents the most
expression for the hydraulic diffusivity in the convection-dif- severe condition. In the context of this study, accuracy refers
fusion equation (see Appendix I) to the avoidance of physically unrealistic outflows and fluc-
tuations in discharge that may occur in response to disconti-
3QI-b' cos 13
(24) nuities in rainfall excess representations. Consistency refers to
D h
= 2(3 + 2b')'1R,So the ability of the routing method to produce the same results
regardless of grid size. This concept is to be used within the
and thus, by matching D n and D h given by (19) and (24),
framework of a simplified routing method in which the objec-
respectively, the weighting factor X can be expressed as a
tive is to optimize the numerical diffusion by matching the
function of channel and flow characteristics, that is
physical diffusion, and it should not be confused with the
x= 1/2(1 - D) (25) property of consistency of a conventional finite difference
scheme, in which the objective is to minimize the amount of
where numerical error in order to properly solve the differential equa-
tion under consideration (Ponce and Theurer 1982).
3QI-b' cos 13
D =---=-----'-- (26) The examined hillslope consists of 20 square cells with grid
(3 + 2b )'1R,Soc !ls
' k spacing !lx, arranged to fonn a single plane (Fig. 4). Hillslope
response to a synthetic step function of rainfall excess is in-
Wave celerity Ck and weighting factor X are varied at each vestigated over a range of rainfall excess rates q, grid spacings
computational cell according to (23) and (25), in which dis- !lx, time-step resolutions !It, mean slopes So, Gauckler-Strick-
charge Q is estimated via a three-point average discharge Q = ler roughness coefficients ks , and channel widths at the outlet,
(Q{ + Q{+I + Q{+1)/3 (Ponce and Yevjevich 1978). By var- Bro. The base case parameter set for the runs and the parameter
ying X with the flow the numerical diffusion coefficient D n is ranges investigated are reported in Table 1. Parameters q, t1x,
used to simulate the hydraulic diffusivity D h of the actual flood t1t, So, ks , and Bro , respectively, are varied, keeping all other
wave. For X = 1/2 there is no numerical diffusion. For X > parameters fixed.
1/2 the numerical diffusion coefficient is negative and the nu-
merical scheme is therefore unstable; however, such occur-
rence is avoided by matching numerical and physical diffusiv-
ities and this ensures the unconditional stability of the scheme.
In the Muskingum-Cunge method, X is interpreted as a dif·
fusion-matching factor and thus negative values of X are pos-
sible.
The Muskingum-Cunge method offers several advantages
over the standard kinematic wave methods. First, the solution
is obtained through a linear algebraic equation, (13), instead
of a finite difference or characteristic approximation of a par-
tial differential equation; this allows the entire hydrograph to FIG. 4. DEM of Hypothetical Hillslope Used to Test Mass Con-
be obtained at required cross section without solving over the servation and Accuracy of Routing Scheme (Darker Shades
entire length of the channel for each time step. Second, the Represent Higher Cell Elevations)
solution using (13) will pennit a more flexible choice of time
and space increments in order to obtain accurate and stable TABLE 1. Model Parameters for HiIIslope Scale Tests
computations as compared to the kinematic wave method. Fi-
nally, the diffusion method provides grid independence for a Parameter BCPS· Parameter ranges
wide range of resolution levels, while the kinematic method (1 ) (2) (3)
does not (Ponce 1986). q.mmh I 10 1, 2. 3, 5, 10, 20
Some of the disadvantages of the Muskingum-Cunge Ax,m 50 10, 20, 30, 50, 100, 200, 300, 500
method are that it cannot handle downstream disturbances that At, h 1/12 1/240, 1/120, 1/60. 1/36, 1/24, 1/12, 1/6
propagate upstream, and that it does not accurately predict the
So 0.084 0.017, 0.084, 0.336
k s , m l /3 S-I 10.0 0.7, 1.4,2.8,5.6, 10.0,20.0,40.0,80.0
discharge hydrograph at a downstream boundary when there B ro , m 2.0 0.5, 1.0, 2.0, 4.0, 8.0, 16.0
are large variations in the kinematic wave speed, such as those a Base case parameter set
that result from inundation of large floodplains (Natural 1975).
106/ JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996

J. Hydrol. Eng., 1996, 1(3): 103-113

 Mass conservation is monitored by integrating the outflow iments performed in this study, it is found that fluctuations
hydrographs for all runs reported in Table 1. The relative mass occur when, in the computational domain, Courant numbers
balance errors are expressed as exceed a threshold value Cu*; an estimate of Cu* is 3. The
hillslope response for the base case parameter set of Table 1
(27)
is shown in Fig. 5. Courant numbers during the simulation
where range from 0.59 to 4.73 and this produces fluctuations in dis-
charge in response to the discontinuities in rainfall excess rep-

Vq = Li q dA dt (28)
resentation. The effect of the Courant number on numerical
dispersion is empirically investigated by varying all the factors
of (32) as reported in Table 1. For example, the effect of in-
is the cumulative volume of rainfall excess produced over the
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

15 - , - - - - - - - - - - - - - r 0.20B
entire hillslope area A during the simulation period T, and
rainfall exce........ (
Idlacharae
Va = i Q dt (29)
-(

0.138 f.'
is the cumulative volume of outlet discharge Q during the 1a
simulation period T. Relative errors in mass balance are found
to be less than 2% for all simulations.
To obtain an accurate solution of the Muskingum-Cunge
method it is generally recommended to respect the empirical
I
!oS;
5 0.08ll f
'V
criterion E!

(30)
whereas negative values for C2 and C3 do not affect the overall o 5 8
accuracy of the method (Ponce and Theurer 1982). Condition
(30) is verified if, and only if,
Cu+D~l (31)

where Cu = ct t1tlt1s is the Courant number; and D = 1 - 2X. -,-------------r 0.208

For wide channels with rectangular cross sections, one has rainfall exc..a ·..···1
Idlacharve -I
Cu = 5/3B-21~k~/~s'lfl0(atlt:lS)o.2" (32)

0. " cos ~
3 0.138 f.'
D=-""::;'''--'':='::-::--
3
5/3k f~S ~3/1OB " as
(33) 1a
Thus, (30) can be expressed in terms of discharge as

Q ~ [5(13 cos- Cu) k~"S~3/10B3"as]513

~
(34)
0.08ll f
When overland flow is examined, (34) may not be attained,
especially in upland areas. In this case the value of Q in (33) -h..,...,r-r-r-r-M-...,.'r'TT'T'T"T........,...,..,...,r-r-1- 0.0
can be forced to meet (34), which is equivalent to the condi- 023 .158
time. t (hr)
tion FIG. 6. HlIIslope Response to Step Function of Rainfall Ex-
D= 1 - Cu (35) ce.. for Base Case Parameter Set of Table 1, Except with 4t =
1/24 h; Fluctuations In Discharge Evident In Fig. 5 Are Signifi-
for calculations of X = 1/2(1 - D). This causes the numerical cantly Reduced (During Simulation, Courant Number Cu
Ranges from 0.29 to 2.33)
diffusion coefficient of the computational scheme to exceed
the hydraulic diffusivity for small values of discharge. Com- 15 - , - - - - - - - - - - - - - r O . 8 3 3
putational results show that for the range of situations that may
occur in a natural catchment, (35) improves the accuracy of
the scheme without significantly changing the features of
catchment response. 0.555 f.'
Since no general criteria for accuracy (in the sense specified "e
previously) are available for the Muskingum-Cunge method
.....
a
with variable parameters, the results of numerical experiments
~
are reported as an alternative indicator for accuracy and con-
sistency. In the linear analysis of the scheme, the numerical
dispersion is minimized when the Courant number is kept
0.277 i
close to 1 (Cunge 1969). In the nonlinear case the Courant
number is allowed to vary at each computational cell of Fig.
2; therefore, it is not possible to establish a priori whether the o +'..,...,r-r-r-r-M-TTT;=.,.,."'I""T"T"'l""'r"T"'I"'T'"+- 0.0
accuracy is preserved. Nevertheless, accuracy is still affected o 2 J 4151
time, t (hr)
by the factors of (32), and the influence of these factors is
FIG. 7. HlIIslope Response to Step Function of Rainfall Ex-
related to the corresponding exponent. For instance, accuracy
is affected by space and time resolution (order 1) more than
ce.. for Base Case Parameter Set of Table 1, Except with 4% =
100 m and B", = 4 m (During Simulation, Courant Number Cu
by topographic slope (order 3/10). From the numerical exper- Ranges from 0.18 to 2.49)

JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996/107

J. Hydrol. Eng., 1996, 1(3): 103-113

 creasing the time resolution (smaller at) is represented in Fig. overland flow that occurs in rivulets, and this value is in agree-
6. For the same set of parameters used for the run plotted in ment with the values derived from the data of Emmett (1978).
Fig. 5, except with at = 1/24 h, the Courant number ranges The Gauckler-Strickler roughness coefficient used in the con-
from 0.29 to 2.33. This caused numerical dispersion to be text of this study is to be assumed as an overall parameter to
reduced and fluctuations in discharge to be significantly de- which many flow processes that occur in a natural hillslope
creased. The same effect can be obtained by varying each of may contribute. This may explain the small values of k s re-
the factors in (32) so as to reduce the Courant number during quired to describe flows in natural channels with respect to
the simulation. In particular, the effects of increasing the grid the case of artificial channels.
spacing ax encourage the model application for large-scale The model was calibrated using data for the flood event
catchment analysis (Fig. 7). shown in Fig. 8. Model validation was then carried out using
data for the events shown in Figs. 9 and 10, where calibration
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

CATCHMENT SCALE MODELING values of soil and channel network hydraulic properties are
maintained; initial conditions of catchment wetness are esti-
The network routing model described earlier was tested on
mated on the basis of precipitation prior to the simulation
the Rio Missiaga experimental catchment, located on the west-
events. A 1-h time step is used for infiltration excess calcu-
ern side of the Cordevole Valley (Belluno, Italy). The catch- lations and a 5-min time step (at = 1/12 h) is used to control
ment's main geological and geomorphological features are re- the accuracy of the routing scheme. The calibration event
ported in Friz et al. (1983). A 4.35-km2 extension of the Rio shown in Fig. 8 is reproduced by the model simulation in both
Missiaga catchment was horizontally discretized into 1,740 the infiltration excess and kinematic subsurface runoff com-
cells with a 50-m grid spacing (ax = 50 m, see Fig. 1). The ponents (see Appendix II). The rising limb and first peak of
Rio Missiaga catchment ranges in elevation from 1,100 m the hydrograph reflect surface runoff response, whereas suc-
above sea level at the outlet to more than 2,400 m above sea
cessive peaks and the tail reflect kinematic subsurface storm-
level. Identification of the channel network from the DEM of flow response. As shown in Figs. 9 and 10, the validation
the catchment was carried out through the critical support area events are reasonably well simulated by the model. The hy-
with constant threshold A * = 0.50 km2 • This threshold value drographs for the two validation events represent predomi-
was selected on the basis of visual similarity between the ex- nantly surface runoff response, as indicated by the single peaks
tracted network and the blue lines depicted on topographic of shorter duration and shorter tails.
maps. Structural parameters of the drainage network and As shown in Fig. 9, fluctuations in discharge that may occur
Gauckler-Strickler roughness are reported in Table 2. Expo- at the elemental cell-channel scale for at = 1/12 h (see Fig.
nents hi and h" introduced earlier are assumed to be hi = 0.26
5) are dampened when propagating downstream, and this ex-
and h" = 0.50, as found by Leopold and Maddock (1953) for plains why the outlet hydrographs do not significantly differ
the semiarid catchments of the United States. These exponent from the hydrographs obtained for extremely high temporal
values provide flow widths that are in reasonable agreement resolution, at = 1/60 h, for which fluctuations do not occur
with the hydraulic geometry of the catchment network, al- (see Fig. 6). The consistency of the method is also confirmed.
though more extensive field data must be collected in order to
fit and verify relationships (1) and (2), and to conduct a com-
TABLE 2. Model Parametera for Rio Mlnlaga Catchment Sim-
prehensive model calibration and parameter optimization. ulations
Surface cover and soil properties are assigned to each DEM
cell on the basis of the available information on catchment Values
attributes and a TCA water balance model is applied to cal- Parameter Overland flow Channel flow
culate rainfall excess at a 1-h time-step resolution. Complete (1 ) (2) (3)
details of the TCA water balance model formulation and par- b' 0.26 0.26
ameterization are given by Orlandini et al. (in press, 1996). b" 0.50 0.50
The nonuniform distribution of soil and vegetation parameters Bro , m 25 7
produces a distributed rainfall excess response of the catch- Q,o, m' s-' 1 1
ment that is qualitatively in agreement with field observations. ks : m l/3 S-l 0.7 5.0
Infiltration excess is the dominant surface runoff production "lin, where n is Manning's roughness.
mechanism for lowland areas, where fine-grained sedimenta-
tion results in relatively low surface hydraulic conductivities, 3.0,.--------------,
although strong subsurface kinematic storm-flow response
from the coarse-grained upland areas can also occur during
extreme storm events. The subsurface storm-flow response of
the catchment is described in this study through the diffusion
wave formulation outlined in Appendix II.
The apparent Gauckler-Strickler roughness coefficient ks is \2.0
the only (distributed) parameter of the routing model. For the E
......
Rio Missiaga catchment, uniform distributions for ks are as-
sumed and, as reported in Table 2, different values of ks are ~
c
.s:.
considered for cells in which overland or channel flow is as- .~ 1.0
sumed to occur (see Fig. 1). This allows one to take into ac- ...
count the different dissipative processes that characterize over-
land and channel flow. Small values of k s assumed to calibrate
the model are in agreement with observed data reported in the
literature. Kouwen and Li (1980) found, for flow through
grass-lined channels, values of Manning roughness n = 0.05- o 16 32 48 64
0.5 m- 1I3 s (k s = 2-20 m ll3 S-I). Newson and Harrison (1978) time (hr)
found that for sheet flow n = 25 m- 1I3 s (ks = 0.04 m ll3 S-I). FIG. 8. Comparison between Simulated and Observed Outlet
Bathurst (1986) used an intermediate value of Gauckler-Strick- Storm Flow for October 9-12, 1987 Rio Mlsslaga Flood Calibra-
ler roughness, ks = 0.75 m l/3 S-l (n = 1.33 m -1/3 s), to describe tion Event

108/ JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996

J. Hydrol. Eng., 1996, 1(3): 103-113

 2.0.,...-------------,
simulated (61-1/12 hr)-
simulated (At= 1/60 hr) _.•..
observed ·Go·

..... 1.5
I
..
1
0 1. 0
~
.
Cl
.c.
o
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

'ii 0.5

o 4 8 12 16
time, t (hr)
FIG. 9. Comparison between Simulated and Observed Outlet
Storm Flow for September 29, 1991 Rio Mlssiaga Flood Valida-
tion Event

6.0 ...,-----....,.~----:-)
- ~-------, FIG. 11. DEM of Rio Missiaga Catchment Showing Selected
simulated b'-0.26 -
simulated b'''O) •.•.. Locations (Black Cells, A-E) and Path (Dark Cells) for Eulerian
observed ...... and Lagrangian Analysis of September 14, 1993 Flood Event

..... TABLE 3. Characteristics of Path Shown in Fig. 11 and Peak

I Discharges at Selected Locations·
.. 4.0
E
...... Sb A '!Ac '!Ad Qpe.k
a Location (km) (km2 ) (m) (m) (m 3 s-')
oS (1 ) (2) (3) (4) (5) (6)
~
Cl
.c. A 1.450 1.3825 3.94 5.32 0.27
lil 2.0 B 1.842 3.1125 5.92 6.46 1.10
'ii C 2.192 3.4675 6.24 6.63 2.02
D 2.592 4.2000 6.87 6.94 3.77
E 3.062 4.3500 7.00 7.00 4.41
aRio Missiaga September 14, 1993 flood event (see Fig. 10).
bDownstream spatial coordinate.
o 4 8 12 16 c Static description (h' = 0).
time. t (hr) dDynamic description (h' = 0.26).
FIG. 10. Comparison between Simulated and Observed Outlet
Storm Flow for September 14, 1993 Rio Mlssiaga Flood Valida-
tion Event ..... 6 . 0 . . . , - - - - - - - - - - - - - - - - ,
'i

In Fig. 10 the effect of the dynamic scaling of flow width on

.!.
In location
location
location
A·····
B·······
C···
the catchment response is investigated. As expected from (23) a location 0-
location E-
and Fig. 3, the change in wave celerity in the case of dynam- g4.0-
ically scaled flow width (h' = 0.26) produces a slower catch- ~
o
ment response and a smaller hydrograph peak in comparison .2
to the case of static flow width (h' = 0). This effect is enhanced c:
CIl
.~
when the discharge increases, and thus further investigations ClI
for larger scale catchment simulations are suggested.

!I'\l
c 2.0 J',
As an example of model application for investigative pur- ~ I,
poses both Eulerian and Lagrangian analysis of the September
14, 1993 flood event are presented. The five locations A-E
~
~
(channel sections) and the flow path shown in Fig. 11 are
selected for the Eulerian and Lagrangian analysis, respectively. ...J.1\·.
. ., 0.0 +.,.-,.-r-:r--1-.,.-r-,..e;,h"::';:''';~., .:\~.~.~
Ir-=;~-I
Some hydrologically relevant features along the selected path o 4 8 12 16
are reported in Table 3. As shown in Fig. 12, the Eulerian time. t (hr)
analysis of the flood event is carried out by plotting the surface FIG. 12. Surface Runoff Hydrographs at Locations A-E
runoff hydrographs at five locations along the path. The La- Shown In Fig. 11, for September 14, 1993 Rio Missiaga Flood
grangian analysis of the flood event is carried out by plotting Event
surface flow characteristics along the selected path, for four
given instants in time. Flow discharges, depths, and velocities variability when using the geomorphologic instantaneous unit
are shown in Figs. 13, 14, and 15, respectively. For the con- hydrograph to model basin response [see Agnese et al. (1988)].
sidered flood event, average velocities increase in a downslope The effect of topography on catchment dynamics can be
direction along the path, and the rate of increase is qualita- expressed in terms of the Peelet number. The Peelet number
tively in agreement with published data [e.g., Carlston (1969),
Pilgrim (1977)]. Therefore, one must properly account for this Pe = 4UYIDh (36)

JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996/109

J. Hydrol. Eng., 1996, 1(3): 103-113

 ~6.0
t=9 hr -
"E
...... t=10 hr---
t=11 hr .......
o t=12 hr -_..

I: ,I
I
..' I

-ID ,}
,I .'
•..""./
"'0 2 . 0 .' ,/
"
~ ~; fi
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

~ ,',' ......)
:ii ,'
I'"
..... ,1
/,.l
i'."·...·•·
",,'

~ 0.0 +.,.--.,.--.,.--..,....:;.i"'."'i~~."".:~:;
. ::::.:•...,.__•... ,_.__.~-F::r-_r__r__r_-.-l
o 1 2 :5 4
space along channel, B (km)
FIG. 13. Surface Runoff Discharge along Path Shown In Fig.
11, for Different Instants in Time of September 14, 1993 Rio Mis-
siaga Flood Event

0.6.,.----------------,
FIG. 16. Map of Hydraulic Diffusivities Dh over Rio Missiaga
t=9 hr - Catchment (See Table 2), as Obtained from Eq. (38) in Which g,
E
...... t-l0 hr--- = 10 mm h-'; Darker Shades Represent Higher Diffusivities: Dh
t=11 hr .......
>- t=12 hr-'-"
Ranges from 0 to 1.93 m' s·' with Average Value of 4.86 x 10-'
m's·'
~ 0.4
I: values of Peelet number, (20) will be mostly of a convective
~ nature. From (36) and (24) the Peelet number can be expressed
'e;,
c as
....c
.r: 0.2
.... Fe = 8So/cos 13 (37)
0..
Gl
"0 Although the steeply sloping Rio Missiaga catchment is char-
~
o
;;:::
acterized predominantly by convection rather than diffusion,
the variability of hydraulic diffusivity of surface flows in re-
sponse to a hypothetical rainfall excess forcing of intensity q,
o 1 2 3 4 = 10 mm h -1 and duration greater than the time of concentra-
space along channel, s (km) tion of the basin
FIG. 14. Surface Runoff Depth along Path Shown In Fig. 11, I-b' b"-b'
for Different Instants in Time of September 14, 1993 Rio Mlssi- - 3q, Ao cos 13 AI-b"
aga Flood Event Dh = 2(3 + 2b I b'
)B,oQ,o So
(38)

,.... 2 . 0 . , . . . - - - - - - - - - - - - - - - , is plotted in Fig. 16 to show how the application of the Mus-

'..§... t=9 hr -
t=10 hr---
kingum-Cunge scheme to the considered channel network is
capable of distinguishing hillslope areas, where convection is
t=11 hr .......
t=12 hr _._.. the predominant physical mechanism, and lowland areas,
~ 1.5 where physical damping may be important. Comparison of

..
~
5 1.0 I-"
" I

/ ......
, ....... ,-
.
Figs. 1 and 16 emphasizes how the model can automatically
identify catchment areas characterized by different dynamics.
The kinematic wave description applies for steeply sloping
.~ I ,
o ~ ~t hillslopes, while the diffusion-convection wave description ap-
.. ' ,r•••••

~ t'" ".l'·· ,., ... plies in lowland areas, where the conditions of slope and flow
o ,,' / ._i depth make physical damping an important effect.
=E 0.5 '
,/(.•./'"
r""

~ SUMMARY AND CONCLUSIONS

~ /~. ·-1" The Muskingum-Cunge routing scheme with variable pa-
;;::: O.0 -+-...,.....,..-I,F~l:;>;=.,=;=::;::~:...._..,.._..,..__r__r__r__r_....__l rameters was applied to DEM-based conceptual transport net-
o 1 2 :5 4 works in order to describe distributed catchment dynamics in
space along channel, s (km) response to storm events. Each elemental channel of the net-
FIG. 15. Flow Velocity along Path Shown in Fig. 11, for Differ- work has bed slope that depends on the DEM cell location,
ent Instants in Time of September 14, 1993 Rio Missiaga Event and cross-section width that is allowed to vary dynamically
according to the Leopold and Maddock (1953) at-a-station and
with U = mean flow velocity and Y = flow depth, playing for downstream relationships. Overland and channel flow is as-
the diffusion-convection equation, (20), the same role as the sumed to occur in those cells for which the upstream drainage
Reynolds number (R = 4UY/v, where v is the kinematic vis- area A is less or greater than a given threshold value A *, re-
cosity) plays for the momentum equation, as a measure of the spectively. The appraent Gauckler-Strickler roughness coeffi-
ratio of the convective flux to the diffusive flux. For very high cient is assumed as an overall distributed parameter to describe
110/ JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996

J. Hydrol. Eng., 1996, 1(3): 103-113

 the different processes that characterize overland and channel width and depth; sand t = spatial and temporal coordinates;
flow. qL = lateral inflow; ~ = channel bed inclination angle; So = sin
Numerical experiments were performed using hypothetical ~ = channel bed slope; and Sf = friction slope (Hayami 1951).
hillslopes in order to test mass conservation and accuracy of Since B = B[Q, (Y)], (41) can be written in the form
the method for possible slope and roughness configurations
that may occur in a natural catchment. Errors in mass balance aQ + B* ay =q (43)
did not exceed 2%. Numerical experiments indicate that ac- as at L

curacy is preserved for a wide range of slope and roughness where it has been defined that
characteristics, space-time resolutions, and rainfall excess forc-
ing rates. Fluctuations in discharge may occur in response to aB
discontinuities in rainfall excess representation if the Courant B* =B + y- (44)
aY
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

number Cu exceeds a threshold of about 3. However, model

application to Rio Missiaga catchment shows these fluctua- Eqs. (43) and (42) can be coupled into a single equation by
tions to be dampened when propagating downstream, so that deriving the first with respect to s and the second with respect
accuracy is generally preserved for catchment scale simula- to t in order to eliminate the mixed derivative a2 YI(asat). On
tions, and also for computationally inexpensive space-time dis- applying the equations
cretizations (Ax = 50 m, At = 1/12 h), for which Cu > 3 may
occur. (45)
The combined use of the DEM-derived transport network
structure, of the scaling concept for hydraulic geometry of and, from (43)
channels, and of the diffusion wave equations for both over-
land and channel flow can help in investigating the nonlinear
dynamics of basin response to storm rainfall, where basin to-
aY = J....
at B*
(q _aasQ)
L
(46)
pography and stream channel geometry play an important role.
Computational experience indicates this approach as an effi- one obtains
cient alternative to the Saint-Venant equations in modeling sur-
face runoff in natural basins, where these equations could be (47)
hardly applied due to insufficient detail in the knowledge of
boundary conditions. Further theoretical work and field in- where
spections are suggested to investigate the representation of hy-
draulic geometry of reaches in natural drainage systems, es- _= _J.... (aSf _ cos ~ aB*)jaSf (48)
pecially at the mid- and large-size catchment scales. Ck B* ay B* as aQ

ACKNOWLEDGMENTS and

This research was jointly supported by the Ministero dell'UniversitA e

della Ricerca Scientifica e Tecnologica of Italy through the grant
D -
h -
lj(~ asf )
cos ~ aQ
(49)
MURST 40% Processi Idrologici Fondarnentali, and by the Gruppo Na-
zionale per la Difesa dalle Catastrofi Idrogeologiche, contribution The wave celerity Ck can be expressed as sum Ck = c; + cZ,
95.00273.PF42. Vigilio Villi (CNR, Padova, Italy) is gratefully acknowl- where
edged for providing data for the Rio Missiaga catchment. Discussions
with Paolo La Barbera (Universitll degli Studi di Genova, Italy), Marco C' = _J.... aSflasf
Mancini (Politecnico di Milano, Italy), and Claudio Paniconi (CRS4, (50)
k B* ay aQ
Cagliari, Italy) are acknowledged and appreciated. The writers thank the
anonymous reviewers for their helpful comments. and
APPENDIX I. CHARACTERISTICS OF FLOW IN WIDE
(51)
RECTANGULAR CHANNELS WITH DYNAMICALLY
VARIABLE WIDTH
One obtains
The Manning-Gauckler-Strickler equation for wide channels
with rectangular cross section (52)

(39) _" - 3b' cos ~ aQ

C (53)
k - 2(3 + 2b')SoOOQb' -as
where Q = discharge; B = flow width; ks = apparent Gauckler-
Strickler roughness coefficient; Sf = friction slope; and Y = and (24) for D h • From computational results, cZ « c; and
flow depth, can be generalized to the case in which the flow thus, neglecting cZ, Ck equals the kinematic celerity Ck given
width varies with discharge by incorporating (10) into (39). If by (23).
one considers the obtained equation
APPENDIX II. DIFFUSION WAVE FORMULATION FOR
(40)
SUBSURFACESTORMFLOW
and (to) to express the flow area 0 = BY, one can solve for To incorporate the rapid subsurface storm-flow response of
Q to obtain the generalized friction (21). The kinematic wave the catchment [see Beven (1981)], a diffusion wave model
celerity Ck = dQ/dO given in (23) is directly derived from (21). similar to the model described earlier is applied. Local con-
The hydraulic diffusivity D h in (20) is derived from the con- tributions to subsurface runoff that result from the soil water
tinuity equation and a simplified momentum equation
balance model described in Orlandini et al. (in press, 1996)
aQ + aO = q' a(Y cos ~) =S _ S are routed on a conceptual transport network extracted from
as at L, as 0 f
(41,42) DEM data. Wave celerity and diffusivity are estimated by in-
corporating Darcy's law into the model instead of the Man-
where Q = discharge; 0 = BY = flow area; B and Y = flow ning-Gauckler-Strickler equation. Subsurface wave celerity is
JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996/111

J. Hydrol. Eng., 1996, 1(3): 103-113

 described by applying the kinematic wave approximation of Li, R. M., Ponce, V. M., and Simons, D. B. (1980). "Modeling rill den-
saturated subsurface flow and assuming that the hydraulic gra- sity." J. lrrig. and Drain. Div., ASCE, 106(1),63-67.
dient equals the topographic slope So. From the flow equation Montgomery, D. R., and Foufoula-Georgiou, E. (1993). "Channel net-
work source representation using digital elevation models." Water Re-
Q+ = 0+ KhS O (54) sour. Res., 29(12), 3925-3934.
Moore, I. D., and Grayson, R. B. (1991). "Terrain-based catchment par-
where Q+ = subsurface storm-flow discharge; n+ = apparent titioning and runoff prediction using vector elevation data." Water Re-
flow area; and K h = hydraulic conductivity in a direction par- sour. Res., 27(6),1177-1191.
allel to the flow, one obtains Mosley, M. P. (1979). "Streamflow generation in a forested watershed,
New Zealand." Water Resour. Res., 15(4),795-806.
(55) Natural Environment Research Council (NERC). (1975). Flood Studies
Rep., Flood Routing Studies Ill, Inst. of Hydro., Wallingford, England.
and Newson, M. D., and Harrison, J. G. (1978). "Channel studies in the
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

Plynlimon experimental catchments." Rep. 47, Inst. of Hydro., Wal-

D: = Q+/(6,B+So) (56) lingford, England.
Orlandini, S. (1995). "Space-time dependence of catchment-scale hydro-
where c: and D; = kinematic wave celerity and diffusivity of logic processes: comparison between physically based distributed mod-
subsurface stormflow; 6, = soil porosity; and B+ =flow width. els at different levels of conceptualization," PhD thesis, DIIAR, Poli-
High hydraulic conductivities are indicated by the small re- tecnico di Milano, Milano, Italy (in Italian).
sponse times of many catchments to subsurface flow, partic- Pilgrim, D. H. (1977). "Isochrones of travel time and distribution of flood
ularly forested catchments [see Mosley (1979), Sloan and storage from a tracer study on a small watershed." Water Resour. Res.,
13(3),587-595.
Moore (1984)]. This may be due to water flowing through
Ponce, V. M. (1986). "Diffusion wave modeling of catchment dynam-
interconnected macropores created by roots, other organic mat- ics." J. Hydr. Engrg., ASCE, 112(8), 716-727.
ter, and cracks between aggregates and rocks rather than Ponce, V. M., and Theurer, F. D. (1982). "Accuracy criteria in diffusion
through the soil matrix. Hence the soil porosity 6, and hy- routing." J. Hydr. Div., ASCE, 108(6), 747 - 757.
draulic conductivity K h used to characterize the subsurface Ponce, V. M., and Yevjevich, V. (1978). "Muskingum-Cunge method
flow may be different from those measured from small soil with variable parameters." J. Hydr. Div., ASCE, 104(12), 1663-1667.
Rinaldo, A., Marani, A., and Rigon, R. (1991). "Geomorphological dis-
cores. For fields application of the model, 6, and K h should be persion." Water Resour. Res., 27(4), 513-525.
considered as effective parameters of the soil profile rather Ross, B. B., Contractor, D. N., and Shanholtz, V. (1979). "A finite ele-
than being derived from the measured values of soil porosity ment model of overland flow and channel flow for assessing the hy-
and hydraulic conductivity for the soil matrix [see Moore and drologic impact of land-use change. " J. Hydro., Amsterdam, The Neth-
Grayson (1991)]. erlands, 41, 11-30.
Sloan, P. G., and Moore, I. D. (1984). "Modelling subsurface storrnflow
APPENDIX III. REFERENCES on steeply sloping forested watersheds." Water Resour. Res., 20(12),
1815-1822.
Agnese, C., D' Asaro, F., and Giordano, G. (1988). "Estimation of the Szyrnkiewicz, R. (1991). "Finite element method for the solution of the
time scale of the geomorphologic instantaneous unit hydrograph from Saint Venant equations in an open channel network." J. Hydro., Am-
effective streamflow velocity." Water Resour. Res., 24(7), 969-978. sterdam, The Netherlands, Vol. 122, 275-287.
Band, L. E. (1986). "Topographic partition of watersheds with digital Wooding, R. A. (1966a). "A hydraulic model for the catchment-stream
elevation models." Water Resour. Res., 22(1), 15-24. problem. I: Kinematic wave theory." J. Hydro., Amsterdam, The Neth-
Bathurst, J. C. (1986). "Physically-based distributed modelling of an up- erlands, Vol. 3, 254-267.
land catchment using the Systeme Hydrologique Europeen." J. Hydro., Wooding, R. A. (1966b). "A hydraulic model for the catchment-stream
Amsterdam, The Netherlands, Vol. 87, 79-102. problem. II: Numerical solutions." J. Hydro., Amsterdam, The Neth-
Beven, K. (1981). "Kinematic subsurface storrnflow." Water Resour. erlands, Vol. 3,268-282.
Res., 17(5), 1419-1424. Wooding, R. A. (1966c). "A hydraulic model for the catchment-stream
Beven, K. J., and O'Connell, P. E. (1982). "On the role of distributed problem. Ill: Comparison with runoff observation." J. Hydro., Am-
models in hydrology." Rep. 81, Inst. of Hydro., Wallingford, England. sterdam, The Netherlands, Vol. 4, 21-37.
Brandford, G. E., and Meadows, M. E. (1990). "Finite element simula-
tion of nonlinear kinematic surface runoff. " J. Hydro., Amsterdam, The APPENDIX IV. NOTATION
Netherlands, 119, 335-356.
Bras, R. L. (1990). Hydrology: an introduction to hydrologic science. The following symbols are used in this paper:
Addison-Wesley, Reading, Mass.
Carlston, C. W. (1969). "Downstream variations in the hydraulic ge-
ometry of streams: special emphasis on mean velocity." Am. J. Sci., A = upstream drainage area;
Vol. 267,499-509. Ao = catchment drainage area;
Cunge, J. A. (1969). "On the subject of a flood propagation computation A* = threshold drainage area;
method (Muskingum method)." J. Hydr. Res., 7(2), 205-230. a' = coefficient of "at-a-station" width-discharge relationship,
Emmett, W. W. (1978). "Overland flow." Hillslope hydrology, M. J. Eq. (1);
Kirkby, ed., John Wiley & Sons, Inc., New York, N.Y., 145-176. a" = coefficient of "downstream" width-discharge relationship,
Friz, C., Gatto, G., Villi, v., and Caleffa, G. (1983). "Groundwater re- Eq. (2);
sources of a typical catchment in the Dolomites area: The Rio Missiaga
Catchment (Belluno, Italy)." Memorie Scienze Geologiche, Vol. 36,
B = channel width;
293 - 315 (in Italian).
Br = reference value of channel width at-a-station, for discharge
Hayami, S. (1951). "On the propagation of flood waves." Disaster Pre- Qr;
vention Res. Inst. Bull. 1, Kyoto, JlI-pan. Bro reference value of channel width at the catchment outlet,
Julien, P. Y., Saghafian, B., and Ogden, F. L. (1995). "Raster-based hy- for discharge Qro;
drologic modeling of spatially-varied surface runoff." Water Resour. b' = exponent of at-a-station width-discharge relationship, Eq.
Bull., 31(3), 523-536. (1);
Kibler, D. F., and Woolhiser, D. A. (1972). "Mathematical properties of b" = exponent of downstream width-discharge relationship, Eq.
the kinematic cascade." J. Hydro., Amsterdam, The Netherlands, Vol. (2);
15,131-147.
Kouwen, N., and Li, R. M. (1980). "Biomechanics of vegetative channel
Cu = Courant number, Cu = ck!i.t/!i.s;
Ck = kinematic wave celerity, Ck = dQ/dO;
linings." J. Hydr. Div., ASCE, Vol. 106, 1085-1103.
Leopold, L. B., and Maddock Jr., T. (1953). "The hydraulic geometry of c; = kinematic wave celerity for subsurface flow;
steam channels and some physiographic implications." Profl. Paper D = cell Reynolds number, D = 1 - 2X;
252, U.S. Geol. Surv., Washington, D.C. D h = hydraulic diffusivity, D h = Q/(2BSo);
Leopold, L. B., Wolman, M. G., and Miller, J. P. (1964). Fluvial processes D: hydraulic diffusivity for subsurface flow;
in geomorphology. W. H. Freeman, San Francisco, Calif. D. = numerical diffusion coefficient, D. = c !i.s(1I2
k - X);
112/ JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996

J. Hydrol. Eng., 1996, 1(3): 103-113

 = spatial discretization index; R = Reynolds number, R = 4UYlv;
j = temporal discretization index; T = simulation time period;
ks = Gauckler-Strickler roughness coefficient (ks = lin); 1 = time variable;
n = =
Manning roughness coefficient (n lIks); Va = cumulative volume of outlet discharge;
Pe = P6clet number, Pe = 4UYIDh ; Vq = cumulative volume of rainfall excess;
Q = surface flow rate, discharge; X = Muskingum-Cunge weighting factor;
Q+ = subsurface flow rate; x = horizontal space coordinate;
Qr = reference value for discharge at-a-station; Y = flow depth;
Qro = reference value for discharge at the catchment outlet; y = horizontal space coordinate;
Q = three-point average discharge in the computational cell; .1s = space interval, along channel;
q = rainfall excess rate, per unit of catchment area; .11 = time interval;
qr = reference value for rainfall excess rate, qr = QrIA; .1x = DEM cell size in the x-direction;
Downloaded from ascelibrary.org by UFRJ - Universidade Federal Do Rio De Janeiro on 05/25/23. Copyright ASCE. For personal use only; all rights reserved.

qL lateral inflow rate, per unit of channel length; .1y = DEM cell size in the y-direction;
Sf = friction slope; EM = relative error in mass balance;
So = channel bed slope; 9, = soil porosity; and
s = spatial coordinate, along channel; n = flow area.

JOURNAL OF HYDROLOGIC ENGINEERING / JULY 1996/113

J. Hydrol. Eng., 1996, 1(3): 103-113