ABSTRACT
A diffusion wave method for catchment dynamics is presented. The
method has better convergence properties than kinematic wave models which
use offcentered derivatives in their numerical formulations. Unlike the offcentered
schemes, the diffusion wave scheme is formulated by matching physical
and numerical diffusivities. This results in an effective control of numeral diffusion
and leads to simulations which are essentially independent of grid size.
The diffusion method is further extended to the realm of dynamic waves by
including the Froude number dependence of the physical diffusivity. The resulting
formulation is believed to represent as complete a description of the
wave dynamics as is possible within the framework of diffusion wave theory.
Numerical experiments show that the diffusion wave method has better convergence
properties than existing kinematic wave methods featuring uncontrolled numerical diffusion.

1. INTRODUCTION
The mathematical modeling of catchment dynamics is well established
in hydrologic research and practice (1, 5, 6, 8, 20). Among the various
aspects of catchment modeling, the conversion of excess rainfall into
runoff is one that perhaps has received the most attention. Procedures to
accomplish this invariably resort to kinematic wave theory (11).
The kinematic wave is a valuable tool in the simulation of catchment
dynamics. With the steep slopes usually encountered in upland
watersheds, flow conditions are such that kinematic waves are reasonably
good approximations of the unsteady flow phenomena. Exceptions are
cases involving mild slopes coupled with fastrising hydrographs, for
which diffusion and dynamic waves are better representations than
kinematic waves.
Procedures for kinematic wave computations are either analytical or
numerical. Analytical solutions provide answers for a simplified class of
problems, while problems of a more general type are handled with
numerical solutions through attendant discretization of the solution
domain. The discretization, however, is not without its pitfalls. Numerical
solutions of kinematic waves are known to introduce variable amounts
of numerical diffusion, with the solution resembling a diffusion wave
rather than a kinematic wave (3). In schemes of first order, this
numerical diffusion is "uncontrolled," resulting in catchment response being
dependent on grid size.
In schemes of second order the numerical diffusion disappears, but a
certain amount of numerical dispersion (of the third order) still remains
(16). The search for a scheme exhibiting zero numerical diffusion and
dispersion may be illadvised because the solution itself stands to benefit
from these contributions. In effect, it has been shown that the diffusion
wave applies to a wider class of problems than the kinematic wave (13).
Therefore, a small amount of numerical diffusion in the kinematic wave
solution is certainly welcome.
The aim should be to "control" the amount of numerical diffusion in
such a way that it matches the diffusion of the physical problem. As will
be shown here, such a methodology has the decided advantage that it
results in catchment response being essentially independent of grid size.
This latter feature substantially enhances the credibility of the results,
leading to better modeling practice.
2. BACKGROUND
The mathematical modeling of distributed catchment dynamics has been
attempted in many ways (1, 5, 6, 8, 19, 20). A widely accepted formulation
is the selection of the kinematic wave equations together with a
simplified spatial representation of the catchment in the form of an "open book."
Such a configuration consists of two planes adjacent to one channel (Fig.
1). Rainfall is the inflow to the planes, and runoff from the planes is by
overland flow in the direction perpendicular to the channel alignment.
Inflow to the channel is by the lateral contribution from the planes, and
the outflow from the channel is the catchment response.
Fig. 1. Openbook catchment schematization (20).


A source of uncertainty in this type of modeling is the adequacy of
the frictional representation. Flow over the planes is usually of such small
depths that laminar flow may prevail. In contrast, flow in the channel
is likely to be of a turbulent nature. Yet in certain cases, especially near
the downstream end of the planes, the flow may be in the transitional,
or mixed laminarturbulent,
regime. For routine applications, friction coefficients are generally estimated
in such a way that they account for laminar flow in the overland
flow planes (6).
The kinematic wave equations are a common choice in the simulation
of flow over catchments. The equations are
well established and generally
lead to physically realistic solutions without the computational and
other complexities that usually plague full dynamic wave formulations.
To derive the kinematic wave equation, the usual statement of conservation
of mass in a control volume is coupled with a simplified form of
conservation of momentum which accounts only for frictional and body
forces. The statement of conservation of mass is the following:
∂Q ∂A
^{_____} + ^{______} = q_{L}
∂x ∂t
 (1) 
in which Q = flow rate, in cfs; A = flow area, in sq ft; q_{L} = lateral inflow, in cfs/ft.
The simplified momentum statement is S_{o} = S_{f} , in which S_{o} = bed slope and S_{f} = friction slope. This in fact a statement of uniform slow, which leads to an alternate way of expressing the simplied momentum statement:
in which α and β are coefficient and exponent, respectively. The values
of α and β contain information on frictional and crosssectional characteristics
[for instance, β = 3 for a wide channel with laminar friction; β
= 5/3 for a wide channel governed by Manning friction; β = 3/2 for a
wide channel with Chézy friction; and β = 4/3 for a triangular channel
with Manning friction (11)].
The speed of propagation of kinematic waves is obtained from Eq. 2:
∂Q Q
c = ^{_____} = β ^{____} = βv
∂A A
 (3) 
in which c = celerity of kinematic waves; and v = mean flow velocity.
The kinematic wave equation is obtained by multiplying Eqs. 1 and 3
and making use of the chain rule to obtain:
∂Q ∂Q
^{_____} + c ^{______} = cq_{L}
∂t ∂x
 (4) 
Equation 4 contains convection terms on the lefthand side and a source
term on the righthand side. Significantly, Eq. 4 does not contain a diffusion
term. This equation may be solved by analytical or numerical
means. A typical class of numerical solutions is that in which the values
of ∂Q/∂t, ∂Q/∂x, c
and q_{L} are expressed in terms of some or all of four
discrete adjacent values of Q, c and q_{L} in space and time (Fig. 2).
Alternatively, the value of c could be kept constant, leading to the linear
kinematic wave solution.
Fig. 2. Spatial and temporal discretization of kinematic wave equation.


A numerical solution of Eq. 4 following the aforementioned procedure
is not without its pitfalls. Unless the spatial and temporal derivatives are
perfectly centered in the rectangular grid, a certain amount of
diffusion will be generated in the numerical solution (Fig. 2). This numerical
diffusion can be eliminated by centering the spatial and temporal
derivatives, but a small amount of numerical dispersion will still remain (16).
The latter is most often simply a nuisance, but in certain cases it may
be of such magnitude as to invalidate the whole solution. Furthermore,
given the fact that there usually is a small amount of diffusion in the
physical problem, the complete elimination of numerical diffusion is
generally illadvised.
Fully offcentered schemes remain popular among practitioners (6, 7, 10).
Ponce et. al. (16) have shown that of the three possible fully offcentered
schemes (see Fig. 3), schemes I and II are conditionally stable and convergent
to the analytical solution as the Courant number (defined as the
ratio of physical celerity c to grid celerity Δs/Δt) approaches 1. Scheme
I, stable for Courant numbers less than or equal to 1, amounts to "upwind differencing."
Scheme II, stable for Courant numbers greater than
or equal to 1, is exactly the opposite image of scheme I. Scheme III is
unconditionally stable, but it is nonconvergent for any value of Courant
number (16).
Fig. 3. Three fully offcentered numerical schemes for solving that kinematic wave equation (16).


Schemes I and II are complementary, with versions of them in wide
use, for instance in the kinematic wave routing option of HEC1 and HECHMS (6).
Scheme III (and related versions), although nonconvergent, has also been
favored because of its features of unconditional stability (7, 10). In any
given practical application, local Courant numbers are likely to vary widely
and be different from 1. If such is the case, none of these three schemes
would solve the kinematic wave equation, certainly not scheme III given
its demonstrated nonconvergence. In theory, for Courant numbers close
to 1, as the discretization is made exceedingly small, schemes I and II
may asymptotically approach the analytical solution of the kinematic wave
equation. With today's powerful computers, this is certainly a possibility.
However, the question still remains whether the complete elimination
of numerical diffusion is the desirable strategy for this type of
problem.
In summary, the hydrologic modeler is faced with the following dilemma:
If he or she seeks to eliminate numerical diffusion, at whatever
cost, the solution may be purely kinematic and therefore lack the natural
diffusion property inherent in physically realistic flows. If, on the other
hand, the modeler uses a fully offcentered scheme such as I, II or III,
the end result will be an uncontrolled amount of numerical diffusion
only vaguely resembling the actual diffusion process.
The way out of this difficulty is to judiciously "match" physical
diffusion with numerical diffusion, a procedure which has been successfully used
in stream channel routing and other similar applications (9, 14, 18). A clear
advantage of this approach is that it renders the numerical solution
essentially independent of grid size. In other words, convergence is
satisfied at a much coarser grid resolution than with either schemes I or II.
Given the natural variation in Courant numbers likely to be present in
practical applications, a method that converges for a wide range of grid
sizes is certainly one to be reckoned with.
3. DIFFUSION WAVE FORMULATION
The soundness of the method of "matched diffusivity" is illustrated
here by comparing the results of two solutions for catchment dynamics.
The first is an application of scheme III to the kinematic wave equation,
Eq. 4, with the associated "uncontrolled" numerical diffusion effect. The
second is an application of the matched diffusivity concept. For
simplicity, both solutions are formulated in a linear mode, with celerity and
diffusivity parameters kept constant. The routing of linear waves visàvis nonlinear
waves has indicated that for small catchments the
nonlinear features are generally small and can be disregarded when
comparing the numerical accuracy of alternative schemes (14).
The application of scheme III to Eq. 4 in a linear mode leads to the
following routing equation:
_{n + 1} _{n + 1} _{n}
Q _{j + 1} = C_{0}Q_{j} +
C_{2}Q_{j + 1} + C_{0}Q_{L }
 (5) 
in which C_{0} = C/(1 + C); C_{2} = 1/(1 + C); and Q_{L} is the average lateral
inflow rate. C is the Courant number of the computational cell, C = c (Δt /Δs), where Δs is either Δx or Δy.
For overland flow computations, Q_{L} is the product of rainfall access
times the area of the overland flow cell. For channel routing, Q_{L}, is the
timeaveraged lateral contribution integrated along the length of a cell.
With an estimate of average wave celerity for plane and channel, together
with the proper specification of initial and boundary conditions,
Eq. 5 enables the calculation of catchment response for a given rainfall
excess excitation.
The derivation of the matched diffusivity method parallels that of
the MuskingumCunge method with the addition of lateral inflow
(3, 14, 15). The routing equations are the following:
_{n + 1} _{n + 1} _{n} _{n}
Q _{j + 1} = C_{0}Q_{j} + C_{1}Q_{j} + C_{2}Q_{j + 1 } + C_{3}Q_{L }
 (6) 
in which C_{0} = (1 + C + D)/(1 + C + D); C_{1} = (1 + C  D)/(1 + C + D); C_{2} = (1  C + D)/(1 + C + D); C_{3} = (2C)/(1 + C + D); and D is the cell Reynolds number defined as the ratio of physical diffusivity q/(2S_{o})
to numerical diffusivity (cΔs)/2, as follows (14):
q
D = ^{__________}
S_{o} c Δs
 (7) 
in which q = unitwidth discharge (q = vd; d = flow depth);
S_{o} =
bottom slope; c = kinematic wave celerity, and Δs = spatial grid
size.
Note that there are major differences between Eqs. 5 and 6 with regard
to their physical data requirements. Equation 5 requires only values of c for
both plane and channel (or alternatively: mean velocity v and a suitable
β value, see Eq. 3). On the other hand, Eq. 6 requires values of c, q and
S_{o} for both plane and channel (or alternatively: β, v, d, and S_{o} ).
Since the physical diffusivity is v = q /(2S_{o} ), it follows that unlike Eq.
6, Eq. 5 does not contain any physical information on diffusion. Not
surprisingly, however, the application of Eq. 5 does show a diffusionlike behavior.
This diffusion, unrelated to the physical problem, is directly traceable to the grid size and is a function of it.
4. DIFFUSIONCUMDYNAMIC WAVE FORMULATION
The preceding analysis has shown the advantages of using diffusion
waves in lieu of kinematic waves in the modeling of catchment response.
Aside from the assumption of linearity, the diffusion wave can
be further improved to include the dynamics of the wave phenomena.
Following Dooge (4), the diffusivity of the diffusioncumdynamic wave
formulation for the case of turbulent Chézy friction in sufficiently wide
channels is:
q F^{2}
v = ^{_____} ( 1  ^{_____})
2S_{o} 4
 (8) 
in which F is the Froude number, defined as F = v/(gd)^{1/2}, and g is the
gravitational acceleration. Equation 8 has the important property that the
physical diffusivity v is zero for F = 2, a threshold of neutral stability at
which both kinematic and dynamic waves propagate at the same speed
(for turbulent Chézy friction in sufficiently wide channels) (11, 12).
Equation 8 can be generalized for any frictional formulation and crosssectional
shape by recalling that the neutral stability Froude number is equal
to 1/(β  1). Therefore, the physical diffusivity expressed in terms of β is:
q
v = ^{_____} [1  (β  1)^{2}F^{2}]
2S_{o}
 (9) 
Given Eq. 9, a diffusioncumdynamic wave model can be formulated
based on Eq. 6, with Eq. 7 modified to account for the Froudenumber
and β dependence of the physical diffusivity, Eq. 9. In effect, for the
diffusioncumdynamic wave, Eq. 7 is modified to:
q
D = ^{___________} [1  (β  1)^{2}F^{2}]
S_{o} c Δs
 (10) 
With the cell Reynolds number defined as in Eq. 10, Eq. 6 ought to
represent as complete a description of the wave dynamics as is possible
within the context of a linear formulation. However, the usefulness of
Eq. 10 as compared to Eq. 7 will depend on whether the wave phenomena
being simulated can be properly construed as a dynamic wave (13).
Experience indicates that for a large class of problems of practical interest
in hydrology, the wave scale may well be in the diffusion (or kinematic)
range.
5. ILLUSTRATIVE EXAMPLE
The methodologies described in the previous section are compared here
by applying them to a typical problem of catchment dynamics. Following
widespread usage, the scheme of Eq. 5 is referred to as "kinematic
wave," although its appreciable and uncontrolled numerical diffusion is
recognized. The "matched diffusivity" method, Eqs. 6 and 7, is herein
referred to as "diffusion wave." The diffusioncumdynamic wave, Eqs.
6 and 10, although properly within the family of diffusion models, is
herein referred to as "dynamic wave" (see Tables 2 and 3).
The example consists of a catchment conceptualized as an open book
with two planes adjacent to one channel (Fig. 1). The planes are assumed to be
impermeable, with rainfall excess (net rainfall) equal to total
rainfall. The inflow to the planes is the rainfall excess, and runoff is by
overland laminar flow in a direction perpendicular to the channel alignment.
Inflow to the channel is by lateral contribution from the planes,
and the outflow from the channel is the catchment response.
The dimensions are: plane length (in the direction of plane flow) =
120 ft (36.58 m); plane width = 240 ft (73.15 m); channel length = 240
ft (73.15 m). The rainfall is 3 in/hr (76.2 mm/h), which, when multiplied
times the catchment area, leads to a flow concentration of 4 cfs
(0.113 m^{3}/s) at the outlet (assuming translation only and
neglecting diffusion). The bottom slopes are the following: slope of planes S_{op
} (in the
xdirection) = 0.01; slope of channel S_{oc} = 0.01.
Laminar flow in the planes and turbulent flow in the channel are calculated
following established hydraulic practice (2). Average flow conditions
in the planes are: discharge per unit width q_{p}= 0.004 cfs/ft (0.00037
m^{3}/s/m); mean velocity v_{p} = 0.5 fps (0.15 m/s);
mean flow depth d_{p} = 0.008 ft (0.0024 m); rating curve
exponent (laminar flow in wide channels) β_{p} = 3; and
mean wave celerity c_{p} = 1.5 fps (0.46 m/s). Average
flow conditions in the channel are discharge per unit width q_{c} = 1 cfs/ft (0.092 m^{3}/s/m);
mean velocity v_{c} = 3 fps (0.91 m/s); mean flow depth
d_{c} = 0.333 ft (0.1 m); rating curve exponent (turbulent flow in triangular
channel) β_{c} = 4/3; and mean wave celerity c_{c} = 4 fps (1.22 m/s).
These calculations lead to a timeofconcentration t_{c} = (120 ft/1.5 fps)
+ (240 ft/4 fps) = 140 s. Therefore, a rainfall duration of 3 min is chosen
to ensure the attainment of equilibrium at the catchment outlet. The total
rainfall volume is 0.15 in. (3.81 mm), or 720 cu ft (20.39 m^{3}).
5. RESULTS
The three methods were tested with regard to their behavior under
various levels of grid resolution. Five resolution levels were chosen,
varying from very coarse to very fine, as shown in Table 1. The results
of the simulations are shown in Tables 2 and 3 expressed in terms of:
(1) peak flow; and (2) timetopeak.
Figure 4 shows the calculated hydrographs for: (1) kinematic method;
and (2) diffusion method. As expected, the results of the kinematic method
varied substantially with the grid resolution. With a very coarse resolution
(A), the peak flow was 2.3507 cfs (0.0666 m^{3}/s), while for a very
fine resolution (E), the peak flow was 3.9641 of (0.1122 m^{3}/s). The
results of the diffusion method showed much better convergence properties
than the kinematic method. With a very coarse resolution (A), the peak
flow was 3.9336 cfs (0.1113 m^{3}/s), while for a very fine resolution
(E), it was 3.9971 cfs (0.1132 m^{3}/s). It is concluded that the diffusion
method provides grid independence for a wide range of resolution levels,
while the kinematic method does not. The demonstrated grid independence
of the diffusion method is clearly a significant asset of this
formulation.
TABLE 1. Tested resolution levels. 
Level (1) 
Δx (ft) (2) 
Δy (ft) (3) 
Δt (sec) (4) 
A 
120 
240 
60 
B 
60 
120 
30 
C 
30 
60 
15 
D 
15 
30 
7.5 
E 
7.5 
15 
3.75 
TABLE 2. Calculated peak flow by three methods (in cubic feet per second). 
Resolution level (1) 
Kinematic (2) 
Diffusion (3) 
Dynamic (4) 
A 
2.3507 
3.9336 
3.9390 
B 
3.0721 
3.9797 
3.9834 
C 
3.5921 
3.9932 
3.9954 
D 
3.8612 
3.9964 
3.9979 
E 
3.9641 
3.9971 
3.9985 
Note: See text for method description. 1 cfs = 0.0283 m^{3}/s. 
TABLE 3. Calculated time to peak by three methods (in seconds). 
Resolution level (1) 
Kinematic (2) 
Diffusion (3) 
Dynamic (4) 
A 
240 
180 
180 
B 
210 
180 
180 
C 
195 
180 
180 
D 
187 
180 
180 
E 
184 
180 
180 
Note: See text for method description. 
Fig. 4. Calculated outflow hydrographs as function of grid resolution (see Table 1 for explanation of curve parameter): (a) kinematic method; (b) diffusion method.

The results of the dynamic method were not plotted in Fig. 4 because
of their close resemblance to those of the diffusion method (see Tables
2 and 3). This resemblance underscores the lack of a significant dynamic
contribution in the tested example. The dimensionless criterion of Ponce
et. al. (13) may be used to prove whether the tested wave is indeed a
diffusion wave. According to this criterion, the wave is a diffusion wave
if the following inequality is satisfied:
g _{1/2}
TS_{o} (^{_____}) ≥ 30
d
 (11) 
in which T is the representative time scale of the wave phenomena. In
this case a time scale can be assumed to be equal to twice the rainfall
duration: T = 360 s. For the channel flow, with S_{o} = 0.01 and d = 0.333
ft (0.1 m), Eq. 11 leads to a lefthand side having a value of 35.4, satisfying
the inequality. This explains why the dynamic wave solution is
not very different from the diffusion wave solution for this particular
problem.
Mass conservation was monitored by integrating the outflow hydrographs
for all runs reported in Tables 2 and 3. All methods and runs
were 100% mass conservative. In addition, the diffusion method was
used with spatial resolution levels satisfying attendant accuracy criteria
(17).
In closing, it should be pointed out that had the discretization been
made exceedingly small, the peak flows calculated by the diffusion and
dynamic methods (see Tables 2 and 3) would have converged to the
peak flow value obtained by concentrating the flow at the outlet 4 cfs
(0.1132 m^{3}/s). This is attributed to the fact that translation (i.e., convection)
is the predominant physical mechanism for the tested example,
with diffusion playing a relatively minor role. Therefore, a diffusion
method (and diffusioncumdynamic extension) in which the diffusion
is small and controlled with physical parameters is certainly to be preferred
over a kinematic method in which the diffusion could be excessive
and bear no relation to the physics of the problem.
6. SUMMARY AND CONCLUSIONS
A diffusion wave method for catchment dynamics is formulated. The
method has better convergence properties than kinematic wave models
which use offcentered derivatives in their numerical formulations. Offcentering
the derivatives introduces numerical diffusion in uncontrolled
amounts, with the results of the simulation being highly dependent on
the chosen grid resolution.
Unlike the offcentered schemes, the diffusion wave scheme is formulated
by matching physical and numerical diffusivity. This results in
an effective control of numerical diffusion and leads to simulations which
are essentially independent of grid size.
The diffusion wave method is further extended to the realm of dynamic
waves by including the demonstrated Froudenumber dependence
of the physical diffusivity. The resulting formulation is believed
to represent as complete a description of the wave dynamics as is
possible within the framework of diffusion wave theory.
The three methods: (1) kinematic, with uncontrolled numerical diffusion;
(2) diffusion, with matched physical and numerical diffusivities;
and (3) dynamic, with matched and Froude numberdependent diffusivity,
are tested for convergence and overall accuracy. The diffusion
method is shown to have much better convergence properties than the
kinematic method. In particular, the demonstrated grid independence
of the diffusion method is an asset to be reckoned with. At least for the
example problem tested here, the solution of the dynamic method is
shown not to constitute a substantial improvement over that of the diffusion method.
APPENDIX I. REFERENCES
Agiralioglu, N., and Singh, V., "A Mathematical Investigation of Diverging Flow, 2. Numerical Solutions and Applications," Report No. MSSUEIRSCE804, Department of Civil Engineering, Mississippi State Univ., Mississippi State, Miss.
Chow, V. T. 1959. Open Channel Hydraulics, McGrawHill Book Co., New York, NY.
Cunge, J. A. 1969. "On the Subject of a Flood Propagation Computation Method
(Muskingum Method)," Journal of Hydraulic Research, Vol. 7, No. 2, 205230.
Dooge, J. C. I. 1973. "Linear Theory of Hydrologic Systems," Technical Bulletin No. 1468, Agricultural Research Service,
U.S. Department of Agriculture, Oct.
Harley, B. M., F. E. Perkins, and P. S.
Eagleson. 1970."A Modular Distributed Model of Catchment Dynamics," Report No. 133, Ralph M. Parsons Laboratory for Water Resources and Hydrodynamics, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Mass., Dec., 1970.
"HEC1, Flood Hydrograph Package, Users Manual. 1981. The Hydrologic Engineering Center, U.S. Army Corps of Engineers, Davis, Calif., Sept.,
revised Jan., 1985.
Huang, Y. H. 1978."Channel Routing by Finite Difference Method," Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY10, Oct., 13791393.
Kibler, D. F., and D. A.
Woolhiser. 1970. "The Kinematic Cascade as a Hydrologic
Model," Hydrology Paper No. 39, Colorado State Univ., Fort Collins, Colo., Mar.
Koussis, A. D. 1983.
"Unified Theory of Flood and Pollution Routing,"
Journal of Hydraulic Engineering, ASCE, Vol. 109, No. 12, Dec., 16521664.
Li, R. M., D. B.
Simons, D. B., and M. A. Stevens. 1975. "Nonlinear Kinematic Wave Approximation for Water Routing," Water Resources Research, Vol. 11, No. 2, Apr., 245252.
Lighthill, M. J., and G. B.
Whitham, G. B. 1955. "On Kinematic Waves I. Flood Movement in Long Rivers," Proceedings of the Royal Society of London, Vol. A229, May, 281316.
Ponce, V. M., and D. B.
Simons, D. B. 1977. "Shallow Wave Propagation in Open
Channel Flow," Journal of the Hydraulics Division, ASCE, Vol. 103, No. HY12, Dec., 14611476.
Ponce, V. M., R. M.
Li, R. M., and D. B. Simons, 1978. "Applicability of Kinematic and
Diffusion Models," Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY3, Mar., 353360.
Ponce, V. M., and V.
Yevjevich, V. 1978. "MuskingumCunge Method with Variable
Parameters," Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY12,
Dec., 16631667.
Ponce, V. M. 1979. "Simplified Muskingum Routing Equation," Journal of the Hydraulics Division, ASCE, Vol. 105, No. HY1, Jan., 8591.
Ponce, V. M., Y. H. Chen, and D. B.
Simons. 1979. "Unconditional Stability in
Convection Computations," Journal of the Hydraulics Division, ASCE, Vol. 105, No. HY9, Sept., 10791086.
Ponce, V. M., and F. D.
Theurer. 1982. "Accuracy Criteria in Diffusion Routing,"
Journal of The Hydraulics Division, ASCE, Vol. 108, No. HY6, June,
747757.
Price, R. K. 1978. "A River Catchment Flood Model," Proceedings, Institution of Civil Engineers, Part 2, Vol. 65, London, England, Sept., 655668.
Singh, V. P., and N. Agiralloglu. "A Mathematical Investigation of Diverging Flow, 1. Analytical Solutions," Report No. MSSUEIRSCE803, Dept. of Civil Engineering, Mississippi State Univ., Mississippi State, Miss.
Wooding, R. A. 1965. "A Hydraulic Model for the CatchmentStream Problem,"
Journal of Hydrology, Vol. 3, 254267.
APPENDIX II. NOTATION
The following symbols are used in this paper:
A = flow area;
C = Courant number, C = c Δt /Δs;
c = wave celerity;
D = cell Reynolds number, Eq. 7;
d = flow depth;
F = Froude number;
g = gravitational acceleration;
j = spatial discretization Index;
n = temporal discretization index;
Q = flow rate, discharge;
Q_{L} = average lateral inflow rate;
q = unitwidth discharge;
q_{L} = lateral inflow rate, per unit of channel length;
S_{f} = friction slope;
S_{o} = bed (or bottom) slope;
s = spatial variable, either x or y;
Δs = space interval, either Δx or Δy;
T = representative time scale of wave phenomena;
t = time variable;
t_{c} = time of concentration;
v = mean flow velocity;
x = spatial variable in the direction parallel to flow over planes;
Δx = space interval in xdirection;
y = spatial variable in the direction parallel to flow in channel;
Δy = space interval in ydirection;
α = coefficient of rating curve, Eq. 2; and
β = exponent of rating curve, Eq. 2.
Subscripts
c = channel; and
p = plane
