Citation
Dynamics and kinematics of overland flows

Material Information

Title:
Dynamics and kinematics of overland flows
Creator:
Huang, Shiu-Ming
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
ix, 93 leaves : illustrations ; 29 cm

Subjects

Subjects / Keywords:
Hydrology -- Mathematical models ( lcsh )
Watersheds -- Mathematical models ( lcsh )
Runoff -- Mathematical models ( lcsh )
Hydrology -- Mathematical models ( fast )
Runoff -- Mathematical models ( fast )
Watersheds -- Mathematical models ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 86-91).
Thesis:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Civil Engineering
Statement of Responsibility:
by Shiu-Ming Huang.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
16855333 ( OCLC )
ocm16855333
Classification:
LD1190.E53 1985m .H82 ( lcc )

Downloads

This item has the following downloads:


Full Text
JDY1J.
AMICS AND KINEMATICS OF OVERLAND FLOWS
by
Shiu-Ming Huang
B.S., National Taiwan University, 1971
M.S., Pennsylvania State University, 1973
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Civil Engineering
1985


This thesis for the Master of Science degree by
Shiu-Ming Huang
has been approved for the
Department of
Civil Engineering
by


Ill
Huang, Shiu-Ming (M.S., Civil Engineering)
Dynamics and Kinematics of Overland Flows
Thesis directed by Assistant Professor James C. Y. Guo
The dynamic-wave model for an overland flow was
assessed using the implicit central finite-difference
method. Flow equations were derived from the continuity
and momentum principles. This study attempts to
delineate the limitations of kinematic-wave approxi-
mation by performing a series of numerical studies using
the flow equations. The nondimensional dynamic-wave
model was combined with Horton's infiltration formula
and numerically solved to establish a model representing
two-dimensional overland flow over a pervious surface.
The comparisons between kinematic-wave solutions and
dynamic-wave solutions provide general guidelines and
further establish new criteria for determining when
kinematic-wave assumptions are valid.
Results from this study indicate that the
overland runoff is mainly kinematic when the rainfall-
kinematic number is greater than 180 or when the
dynamic-wave number is less than 6.


This abstract is approved as to form and
content.
)
Faculty member in charge of thesis


V
ACKNOWLEDGMENTS
I would like to express my deep gratitude to Dr.
James C. Y. Guo for his advice, dedication, and con-
tinuous guidance throughout the research and the
duration of the thesis work.
Thanks are also due to other faculty members of
the University of Colorado, including Dr. David W.
Hubly, Dr. William C. Hughes, and Dr. John Liu, who
introduced me to the aspects of water resources
engineering.
Special thanks go to my husband, Chi-I, and my
sons, David and Samuel, for their encouragement,
support, and understanding.
This research was supported in part by the
Education Program of the Martin Marietta Denver
Aerospace.


vi
TABLE OF CONTENTS
CHAPTER PAGE
I. INTRODUCTION . . ............. .......... . 1
1.1 Description of Problem .............1
1.2 Goals and Objective ................... 6
1.3 Literature Review .......... 7
II. ANALYTICAL MODEL ........... ........ 11
2.1 Introduction ........................11
2.2 Continuity Equation ................. 12
2.3 Momentum Equation ................... 16
2.4 Horton's Infiltration Model ..... 25
2.5 Nondimensionalization
of Governing Equations ............... 33
2.6 Initial and Boundary Conditions ... 36
III. NUMERICAL STUDY .............. 38
3.1 Introduction..........................38
3.2 Finite Difference Method ...... 42
3.3 Central Difference Method ...... 45
3.4 Implicit and Explicit Methods .... 46
3.5 Numerical Formulation
of the Flow Equations ........ 50
3.6 Numerical Computation Procedures . 57


vii
CHAPTER PAGE
IV. CASE STUDY...................................60
4.1 Description of Study Aspects .... 61
4.2 Comparisons between Predictions
of Dynamic-Wave Model and
Kinematic Wave Model ... 63
4.3 Evaluation of Kinematic-Wave
Model Applicability with
Available Laboratory Data .... 76
V. SUMMARY AND CONCLUSIONS .......... 83
5.1 Conclusions .....................33
5.2 Recommendations for Further Study . 35
REFERENCES .................... ............. 86
APPENDIX I NOTATIONS ........ ............. 92


LIST OF TABLES
viii
TABLE PAGE
2.1 Recommended Horton's Equation Parameters . 31
4.1 Illustration of Laboratory
Watershed Analysis ................ 82


LIST OF FIGURES
ix
FIGURE PAGE
2.1 Illustration of Overland Flow Plane .... 13
2.2 Control Volume to Derive
Continuity Equation ....................... 15
2.3 Control Volume to Derive
Momentum Equation 18
2.4 Schematic Representation
of Horton's Equation ...................... 30
3.1 Illustration of Finite-Difference Method . 44
3.2 Network of Grid Points in
Computational Domain ...................... 48
3.3 Flow Chart of Computational Procedures ... 58
4.1 Illustration of Outlet Hydrograph A .... 64
4.2 Illustration of Outlet Hydrograph B .... 67
4.3 Illustration of Outlet Hydrograph C .... 69
4.4 Illustration of Outlet Hydrograph D .... 70
4.5 Variation of the Peak Depth Ratio
versus Rainfall-Kinematic-Wave Number . 71
4.6 Effects of Rainfall-Kinematic-Wave
Number on the Outlet Flow Depth.............74
4.7 Effects of Dynamic-Wave Number
on the Outlet Flow Depth ...................75
4.8A Hydrograph of Laboratory Data, Set 1 . . . . 78
4.8B Hydrograph of Laboratory Data, Set 1 . . . . 79
4.9A Hydrograph of Laboratory Data, Set 2 ... . 80
Hydrograph of Laboratory Data, Set 2 . .
4.9B
81


CHAPTER I
INTRODUCTION
1.1 Description of Problem
The surface runoff produced by rainfall is
hydraulically classified as unsteady, nonuniform flow,
meaning that the velocity and depth of flow vary with
respect to both time and space. In a catchment, the
lowest point, or the outlet of the basin,accumulates the
runoff from the entire drainage area. This results in
the nonuniformity of the flow movements. In addition,
the time variations in rainfall excess also cause
unsteady and nonuniform flow mechanism.
With such a complexity in nature, it is not
possible to use the assumptions of steady uniform
channel flow, inherent in Manning's equation or Chezy's
equation, to properly compute the runoff. In hydraulic
engineering practice, the prediction of runoff is very
fundamental for water-resources planning.


2
The theoretical hydrodynamic equations governing
overland flow are generally attributed to Barre de St.
Venant and were formulated in the late 19th century.
Due to their mathematical difficulty, however, the St.
Venant equations have not yet been solved analytically.
With the advent of high-speed digital computers,
numerical schemes have been developed and used to obtain
discrete solutions for complicated open-channel flow
situations. Prior to the use of high-speed computers,
various simplified forms of the St. Venant equations
were used to approximate different hydrological
problems.
The kinematic-wave approximation is the most
common simplification of the St. Venant equation. It
assumes that the friction slope, Sf,is balanced by the
channel-bottom slope, SQ, but it ignores the pressure
and inertial terms. Lighhill and Whitham (1955) have
found that the flood-wave propagation in overland flow
is mainly kinematic.' Wooding (1965) applied the
kinematic-wave approximation to solve the partial-
differential governing equations describing the flow for
a v-shaped chatchment flow. His solutions represented
complete runoff hydrographs and the time of concentra-
tion. The general solution for the kinematic-wave
approximation of overland flow is similar to Manning's


3
equation, in which, q, the flow rate per unit width, is
a power function of flow depth, y:
q = o(ym (1-1)
in which o( = watershed characteristic constant, for a
uniform flow, &< = (1/n) 1.49 (Sq)1/2
m varies betweem 5/3 and 3, depending upon where
the flow regime is turbulent or laminar.
Wooding (1965) also concluded that the kinematic-wave
method was applicable to overland flow when the Froude
number at outlet is less than 2. In his approach,
however, the infiltartion rate was assumed to be
negligible or constant throughout the runoff period.
As mentioned above, the kinematic-wave theory is
the simpliest form of the complete dyanamic-wave model.
With the simplification, kinematic-wave theory requires
the least computational efforts in predicting an
overland flow runoff hydrograph. Neglecting the dynamic
flow factors, including inertial force and backwater
effects, however, the kinematic-wave theory has its
limitation. Woolhiser and Liggett (1967) determined the
limitations of kinematic-wave theory on the basis of a


4
series of computer simulation studies of the dynamic-
wave model. They concluded that for a given uniform
rainfall excess, ie, the applicability of kinematic-
wave theory may be best judged by the kinematic-wave
number, Kq, which is defined as:
S L
kq=_o (1.2)
H F.
o o
in which SQ = slope of plane
Lq = length of plane (ft)
Hq = normal depth at equilibrium runoff flow
rate, q
Fq = Froude nubmer at equilibrium runoff flow
rate, q and depth, HQ, that is, FQ = Q
:f9H0)1/2
(where VQ = q/HQ)
The study shows that the kinematic wave approximation is
within 10% of the kinematic approximation if Kq is
greater than 10.
These two significant advances in the develop-
ment of kinematic-wave approximation were based upon the
assumption of constant infiltration loss and uniform
rainfall intensity. In reality, the infiltration rate
of a watershed has a major influence on the volumetric
rate of runoff generated from a watershed under a given


5
rainfall event. A constant infiltration rate occurs
only when the catchment is saturated before the rainfall
begins. Therefore, the assumption of the constant
infiltration rate established the limitation of
Wooding's conclusion. Even with a uniform design
rainfall, the rainfall excess is no longer uniform when
the time variation of infiltration rate is considered.
In general, the mathematical model to describe
the decay of infiltration rate with respect to time is
empirically derived from the observation of infil-
trometer tests. There are several infiltration models
available, such as those of Horton (1939), Holtan
(1968), Philip (1969), Green (1911), and Soil
Conservation Service (SCS). Guo (1984A, 1984B) combined
the kinematic-wave approximation with Horton's infiltra-
tion formula to analytically solve the propagation of a
two-dimensional overland flow on a pervious surface. He
reformulated the equation of time of concentration and
concluded that Wooding's assumption of constant infil-
tration rate is a special case of the general situation.
Nothing, however, was advanced regarding the establish-
ment of kinematic-wave approximation limitations under a
nonuniform rainfall excess due to infiltration decay
characteristics.


6
This study attempts to determine the limitations
of kinematic-wave approximation by performing a
numerical study which combines the dynamic-wave model
with Horton's infiltration formula to numerically solve
the propagation of two-dimensional overland flow on a
pervious surface. The results from comparisons between
kinematic-wave solutions and complete dyanamic-wave
solutions can provide general guidelines, and can
further establish new criteria for determining when
kinematic-wave assumptions are adequate and acceptable.
1.2 Goals and Objective
The purpose of this study was to develop a
computer simulation model to solve the dynamic-wave
equations representing overland flow resulting from a
uniform rainfall and a variable infiltration rate. The
discrete computational solutions were then compared with
the kinematic-wave solutions developed by Guo (1984A and
1984B) in order to delineate the limitations of the
applicability of kinematic-wave solutions.
The kinematic-wave method is an approach which
has been proven to be a very useful tool in storm-water
modeling. Consideration is given to flow unsteadiness
without solving the complete St. Venant equation.
Ignoring backwater effects, however, imposes limitations


7
upon its applicability. Kinematic waves occur when the
dynamic terms in the momentum equation are negligible.
Under these circumstances, flow rate can be expressed as
a function of flow depth at all distance (x) and time
(t).
The basic objectives for this study are the
following.
(1) To develop a numerical routing model which
made use of the complete St. Venant equation which
includes Horton's infiltration formula. An implicit
formulation is used to avoid the instability associated
with explicit schemes.
(2) To develop criteria for determining whether
the kinematic-wave method is applicable. As kinematic-
wave method is less complex and therefore^ easier to
model, it is desireable to use this method when it gives
acceptable results. Test runs are made for various
conditions on a hypothetical system, and results of the
two models are compared to determine under what
conditions kinematic-wave assumptions are acceptable.
1.3 Literature Review
Extensive literature reviews regarding overland
flow have been discussed by Linsley et al. (1982) and
Chow (1959). Only a few references pertaining to this


8
study were found, for instance, Richey (1954), Behlke
(1957), Liggett (1959), Izzard (1946), Woo and Brater
(1962), Isaacson et al. (1956), Morgali and Linsley
(1965), Woolhiser and Liggett (1967), Morgali (1970),
Schreiber and Bender (1972), Gburek and Overton (1973),
and Guo (1984A and 1984B).
Isaacson et al (1956) used the finite-difference
mathematical scheme on a digital computer for synthesiz-
ing flows. Morgali and Linsley (1965) then performed a
numerical analysis, using the computer for overland
flow. They concluded that an overland-flow hydrograph
can be synthesized by using an appropriate mathematical
model and numerical methods. The influences of various
physical parameters (such as slopes, roughness, plane
length, and rainfall intensities) on the overland-flow
hydrograph were numerically tested. They proposed a
semi-theoretical formula for estimating the time of
concentration under a uniform rainfall excess.
Woolhiser and Liggett (1967) used the finite-
difference intergration to solve the rising hydrograph
of the unsteady, overland flow. A dry channel was used
as an initial condition; the upstream and downstream
boundary conditions were, respectively, zero inflow and
critical flow depth (or no condition for supercritical
flow). The study concluded that when compared with


9
previous numerical, analytical, and experimental
results, the results show that in general there is no
unique dimensionless rising hydrograph for overland
flow, but that for most hydrologically significant
cases, the kinematic-wave solution gives very accurate
results.
Morgali (1970) did a follow-up study on laminar
and turbulent overland hydrographs. In setting up the
solution of finite-difference schemes, the approximation
of the resistance to flow, that is, the friction slope,
was found to be an important factor. Schriber and
Bender (1972) developed a mathematical model that can be
used to simulate an overland-flow hydrograph for a given
rain strom and watershed. A dimensional analysis is
used as an aid to determine predictive relationships for
the resistance of an overland-flow plane. Giburek and
Overton (1973) studied subcritical kinematic flow in a
stable stream. An evaluation of the shallow-water
equations in the lower 6.9 miles of Mahantago Creek,
Pennsylvania, showed that subcritical kinematic
conditions exist.
Guo's studies (1984A and 1984B) on effects of
infiltration on a runoff hydrograph established the
analytical solutions for a kinematic overland flow on a
pervious surface by combining kinematic-wave theory with


10
Horton's infiltration model. The accuracy of solution
is confirmed by numerical computations based on an
implicit finite-differnece scheme and experimental
data. Results indicate that the flood travel time for
overland flow on a pervious surface is not a constant
time interval and that the peak flood flow rate occurs
when rain ceases.
Recent trends in model developments of the
rainfall/runoff relationship recommended by the Urban
Drainage and Flood Control District (1982) indicate a
tremendous need for an efficient, but more complete,
runoff model. Such models can be valuable tools for the
planning, design, and management of practical problems
due to rainfall. It is the objective of this study to
incorporate a decay-infiltration model to the unsteady
dynamic wave routing and to further investigate the
limitations of the kinematic-wave approximation.


11
CHAPTER II
ANALYTICAL MODEL
2.1 Introduction
The governing equations of a dynamic-wave model
for an overland flow are derived from the continuity and
momentum principles. The equations were first applied
to the overland-flow problem by Keulegan (1944), and
later Morgali and Linsley (1965) did a study on computer
analysis of overland flow. The schematic illustration
of two-dimensional flow is shown in Figure 2.1. The
accumulation of runoff on a unit-width plane with a
slope, SQ and length, L, results from rainfall of
constant intensity, i.
The basic assumptions made in the flow equations
are as follows:
1. The plane has a unit width.
2. The hydraulic radius, Rp, is equal to the
depth of flow, Y.


12
3. The momentum coefficient is equal to 1, and
the energy correction coefficient associated
with velocity head is equal to unity.
4. The pressure distribution is hydrostatic at
any vertical plane in the flow.
5. "The slope of the plane is so small that
cos SQ = 1.
6. The momentum of incoming rainfall in the
direction of flow is negligible.
7. The friction loss in unsteady nonuniform,
open-channel flow is equal to friction loss
with steady uniform flow.
8. Surface slope and roughness are uniform
throughout the catchment.
9. Rainfall intensity is uniform in space.
10. Soil coverage is homogeneous in the entire
catchment so that infiltration rate varies
only with respect to time.
2.2 Continuity Equation
The governing equation of continuity for
unsteady flow may be derived by the conservation of mass
in a control volume, that is, a small space between the
channel sections. Considering the concept of control


\y\\
13
i, rainfall intensity*
urn 1

V
x = 0
upstream
boundary
urtrrr
-X = L
downstream
boundary*
---v x
i = rainfall intensity
y = water-surface depth.
SQ= slope of plane
V = velocity
ft= infiltration loss
Figure 2.1 Illustration of Overland Flow


14
volume between stations 1 and 2 as shown in Figure 2.2,
the principle of conservation of mass states that:
inflow outflow = rate of change in volume
stored
Therefore, the volume increased from time t to time (t
+ &t), can be represented by the area, ABCD, which is
equivalent to:
(q-j-At-qQAt) + (ieAt) ax = Ax Ay (2.1)
in which q^ = average inflow between time t and time
(t + At) (cfs/ft)
qQ = average outflow between time t and time
(t+At) (cfs/ft)
At = time increment (sec)
Ax = distance increment (ft)
Ay = average water surface depth increment (ft)
ie = average rainfall intensity excess
(ft/sec), which is the difference between
the rainfall intensity and infiltartion
loss


15
1
SQ = slope of plane
i = rainfall intensity
(T) = station number
= average inflow, between time t and timeCt+AtL
= average outflow between time t and time(t+At)
OBC = water surface at time t
OAD = water surface at time (..t+At)
ft = infiltration loss
Figure 2.2
Control Volume to Derive Continuity
Equation


16
Dividing Equation (2.1) by At ax yields
g-T-Sr
X
T Xe
(2.2)
Considering infinitesimal spacing in distance and time,
equation (2.2) becomes
M + M. = i-
3X 3t e
(2.3)
Equation (2.3) is the continuity equation for a
two-dimensional overland flow. When rainfall intensity
is neglected, equation (2.3) reverts to conventional
form (Chow, 1959).
+...££ =0 (2.4)
3X 9t
2.3 Momentum Equation
The momentum equation for a two-dimensional flow
may be derived from Newton's second law of motion.
According to that law, the change of momentum in the
waterbody is equal to the summation of all of the exter-'
nal forces that are acting on the body. The external
forces acting on a control volume are shown on the flow


17
profile in Figure 2.3. These external forces include
pressure gradient, body force, and the resistive force.
The most simple fora of Newtons law of motion
states that
.DU
F m a = si
Dt
(2.5)
where
sura of all external forces in a given
direction (in our case, the flow direction)
9 *1 t \
( j.O }
m = mass (slu*
~ accel a rat:
7T u - cross -seci
.1. u A H II (sec;
,uic J median ics,
m = J> AX Y
/here f density of fluid (slug/ft *J)
Ax * distance increment (ft)
y"~ average depth of control volume (ft)
Assuming that the momentum brought to the overland flow
by the rain is negligible, the axtax-nal forces acting or
the control volume between sections 1 and.2 in Figure
2.3 are composed mainly of four components: pressure,
gravity, friction force, and the momentum change due to
rainfall excess.


18
Fg = gravity force
Fp = pressure
i = rainfall intensity
T = specific weight of fluid (water)
^t = infiltration loss
Figure 2.3 Control Volume to Derive
Momentum Equation


19
The component of pressure, Fp, for a unit-width
overland flow can be represented by;
Fp = -TY Ay- T YAY (2.7)
where f = specific weight of the fluid (lb/ft3)
y" = average depth of control volume (ft)
Ay = increment of water depth (ft)
Equation (2.7) is legitimate when the channel bottom
slope is mild. For a very steep slope, the hydrostatic
force may be corrected by multiplying this term by
cos Sc.
The gravity force, Fg, is the component of the
weight of water in the flow direction. The value is
Fg = T sin 9 (y Ax) (2.8)
For a small slope,
sin 9 tan 8-SQ; (2.9)
therefore,
Fg T S0 (yAx) (2.10)
in which SQ = channel bottom slope
y == average depth of control volume (ft)
Ax = distance increment (ft)
The friction force, Ff, is the shear force on
the wetted perimeter.which may be expressed as
Ff = T Sf (yAx)
(2.11)


20
in which Sf = friction slope
Y = average depth of control volume (ft)
Ax = distance increment (ft)
It is assumed that the friction slope, Sf, computed by
the uniform-flow concept can be used for the nonuniform,
unsteady, and spatially varied flow. This assumption
will be discussed in more detail below.
The momentum change due to rainfall excess is
equal to ( -U ie) ax. Where U is the cross-sectional
average velocity, ie, is the average rainfall
intensity excess, and ax is the distance increment.
Summing these four external forces, the change
of momentum in the flow direction in Figure 3.2 can be
expressed as :
F = ma = (y*Axy) a
= AP + Ff + Fg (Uie)Ax (2.12)
Dividing Equation (2.12) by (y x) yields
F = 'JY + rso rSf - (2.13)
' y
The acceleration term, a in Equation (2.5), can be
further decomposed to local acceleration term, 9U/at,
and the convective acceleration term, UdU/a-x.


21
Substituting tha force term/ F, in Equation (2.13) into
Equation (2.5), the complete momentum equation for
unsteady, two-dimensional overland flow is
+ g 22. + = g (Srt S?) (2.14)
at 3x ax y v v '
where g = gravitational constant (ft/sec2)
SQ = channel bottom slope
Sf = friction slope.
When rainfall excess is negligible, Equation-
(2.14) is reduced to the Saint-Venant Equation, that is,
Equation (2.15).
!l + U^i + g!I + g (Sn S) = 0 (2.15)
9t 0X 3x 1
(1) (2) (3) (4) (5)
Equation (2,15) is a nonlinear, hyperbolic
partial differential equation. Due to tha complexity in
mathematics, this equation is usually solved numerically
with specified boundary conditions. No general analyti-
cal solution exists. When the vertical acceleration is
negligible, the hydrostatic pressure and frictional
resistance for an unsteady flow are the same as those in
a steady flow. The terms in Equation (2.15) are,
respectively:


22
(1) the acceleration due to time variation.
(2) the acceleration due to space variation.
(3) the force due to water-surface variation.
(4) the force due to the gravity, that is body
force
(5) the force due to the shear stress on wetted
perimeter.
Various approximations of Equation (2.15) were
used, depending on the terms of the momentum equation
which are considered (Yen, 1978). The kinematic-wave
approximation is the most simple form where the friction
slope is set to be equal to the channel bottom slope,
and the pressure and inertial terms are negligible. The
general format of the solution for a kinematic-wave flow
is similar to the Manning's equation.
q =o< ym (2.16)
in which c* = (1.49 / n ) (Sq)1/2 (Manning's
equation), n = plane roughness
or oi = C (Sq)1/2 (Chezy's equation)
The value of m varied from 5/3 to 3, depending on the
flow regime, turbulent or laminar. It has been found
that making the value of m = 2 provides a good agreement
with the observed data. The kinematic-wave
approximation cannot model the flood-wave attenuation,


23
and it is not adequate when the backwater effects are
significant. The kinematic-wave model is accurate for a
slowly rising flood wave on a steep slope such as
overland flow (Ponce et al., 1978).
The diffusive-wave approximation assumes that
DU/Dt = 0 and is often called the equation of gradually
varied flow. Gradually varied flow is a steady-flow
situation where the depth gradually varies along the
length of the channel. The improvement of
diffusive-wave approximation is that it includes the
consideration of a pressure term. In flood
computations, peak attenuation and backwater effects can
be simulated while two boundary conditions are
specified. The quasi-steady dynamic-wave model neglects
the local unsteady-velocity term, ^U/^t = 0. This model
was found to be valid to describe a long flood-wave
propagation in a channel when lateral inflow is
negligible.
The complete dynamic wave model includes all
terms shown in the momentum equation (equation 2.14).
Backwater effects, downstream effects, and flow reversal
can be simulated. To achieve numerical solutions for
the complete momentum equation is time consuming and
difficult in programming. Frequently, numerical
instability causes divergency in computation. The


24
complete dynamic-wave model has been used for both river
routing (Chen, 1973; Simons and Chen, 1979; Smith, 1980;
and Tucci, 1978) and sewer routing (Labadie et al.,
1978; Yen, 1973; and Yevjevich and Barnes, 1970). In
this study, the full dynamic-wave model is used for
numerical analysis to completely model the dynamic
behavior of overland flow.
For mathematical simplicity, the momentum
Equation (2.14) can further be transformed to a
differential equation in terms of flow rate, q, and flow
depth, y, only. Summing Equation (2.3) multiplied by U
and Equation (2.14) multiplied by y yields:
2.2
£2 + -L (3_+ 22_) = gy (SQ Sf) (2.17)
3t 9x y 2 o i
An examination of Equations (2.3) and (2.17)
shows that there are three unknowns: flow depth, y,
flow rate, q, and friction slope, Sf. In order to
obtain numerical solutions for the two unknowns, flow
depth, y, and flow rate, q, in the continuity equation
(2.3) and the momentum equation (2.17), an empirical
equation for estimating the friction slope, Sf, is


25
required. In this study, Manning's equation for a
steady nonuniform flow is used to evaluate the friction
slope Sf.
2.22y 3
in which n = Manning coefficient
q = flow rate (cfs/ft)
y = water surface depth (ft)
(2.18)
Equation (2.18) is use by assuming that the head loss at
a section for the nonuniform, unsteady flow is the same
as that for a steady nonuniform flow having the same
velocity and hydraulic radius of the section. This
assumption has never been precisely confirmed by either
experiment or theory, but errors arising from it are
believed to be small compared with those ordinarily
incurred in the use of the uniform-flow formula. Over
years of use, this assumption has been proven to be a
reliable basis for engineering designs.
2.4 Horton's Infiltration Model
For a given return period of rainfall, rainfall
intensity, i (in/hour), varies inversely with the
duration of the storm (Overton and Meadows, 1976):


a
26
1 =
a+t
(2.19)
in which Tr = return period (year)
= rainfall duration (minute)
a, b = parameters varying with return period and
geography.
The penetration of water through the soil surface is
known as infiltration. The infiltration rate or rate at
which water enters the soil at the surface is controlled
by surface conditions. Therefore, soil type is an
important factor in determining the infiltration rate.
When the soil has a large portion of well-graded fines,
the infiltration is low. If the soil has several layers
of horizons, the least-permeable layer will sometimes
control the steady infiltration rate.
The soil coverage also places an important role
in determining the infiltration rate. Vegetation, lawn
grass in particular, tends to increase infiltration by
loosening the soil near the surface. The antecedent
moisture condition is also important in determining the
initial infiltration rate. When rainfall occurs on an
areas that has little, antecedent moisture, that is, when
the ground is dry, the infiltration rate is higher than
that for a wet antecedent moisture such as from a
\
previous storm or from irrigation. Although antecedent
moisture condition is very important when calculating


27
runoff from less intense storms in nonurbanized areas,
the data collected in the urbanized basins by the Urban
Drainage and Flood Control District, Colorado, indicate
that the antecedent mositure condition has a limited
effect on runoff (Drainage Criteria Manual, 1982) .
Other factors affecting the infiltration rates include
slope of land, temperature, quality of water, age of
lawn, and soil compaction (Drainage Criteria Manual,
1982) .
The infiltration research has been conducted for
many years, but there still is no generally recognized
adequate quantitative model of natural infiltration.
The reason is that infiltration is affected by the com-
plex combination of soil characteristics, soil moisture
conditions, and other factors occurring in nature
(Overton and Meadow, 1976; Sheaffer et al., 1976.)
However, infiltration is extremely important, because it
can influence not only the volume and intensity of
rainfall excess rate but also the timing of the runoff
hydrograph. Some alternate represtations of infiltra-
tion are as follows;
Horton's Model
Holtan's Model
Soil Conservation Service (SCS) Method
Algorithms based on Darcy's law and/or


antecedent conditions
Simple Abstraction Method
These models range widely in complexity. Each of them
permits initial infiltration rates to depend on some
measure of antecedent soil moisutre and to vary with
time during the storm.
In the semi-arid region, such as the state of
Colorado, the infiltration model proposed by Horton
(1935) was found to provide a good balance between
simplicity and reasonable physical description of the
infiltration process by the Urban Drainage and Flood
Control District in Denver (Drainage Criteria Manual,
1982). Horton's infiltration model states:
ft fc + where ft = infiltartion rate at time t (in/hr)
fQ = initial infiltration rate (in/hr)
fc = final infiltration rate (in/hr)
e = natural logarinthm base
k = infiltration decay constant (1/sec)
t = time (sec)


29
Horton's equation is illustrated graphically in
Figure 2.4. In developing this equation, Horton
observed that infiltration is high at the beginning of
the storm and decays exponentially to a steady state
constant value, fc, as the pores in the soil become
saturated. The coefficients and initial and final
infiltration values are site specific and depend on the
soil and vegetative cover complex.
The infiltration rates are "capacity" rates,
indicating that the soil-air interface must be saturated
at all times. In practical terms, this means that it is
assumed that rainfall rate is always greater than
infiltration capacity rates, and hence some ponding will
always result. This is a major disadvantage in the use
of Horton's model, because rainfall rates are highly
variable and therefore will possibly fall below the
values of capacity infiltration rates (Linsley et al.
1982; Overton and Meadows, 1976).
The Urban Drainage and Flood Control District
has analyzed a considerable amount of rain/runoff data.
On the basis of this analysis, the values in Table 2.1
are recommended for use within the District in
conjunction with the 1982 version of the Colorado Urban
Hydrograph Procedure (CUHP). Most frequently found


INFILTRATION, f IN INCHES PER HOUR
30
Figure 2.4 Schematic Representation of
Horton*s Equation


31
TABLE 2.1 Recommended Horton's Equation Parameters
Decay
SCS Hydrologic Infiltration (in/hr) Coefficient
Soil Grouo* Initial fQ Final fc k (1/sec)
A 5.0 1.0 0.0007
B 4.5 0.6 0.0018
c 3.0 0.5 0.0018
D 3.0 0.5 0.0018
Remark: *SCS = Soil Conservation Service


32
within the District are SCS Hydrologic Soil Groups C and
D; however, areas of Group A and B soils may also be
found within the District (Drainage Criteria Manual,
1982). Hortons Model is therefore selected in this
study by using the available data source given by the
Urban Drainage and Flood Control District as shown in
Table 2.1.
Assuming a spatially and temporally uniform
rainfall intensity, i, over the entire watershed, and
assuming that the slope, surface roughness, and soil
texture are all uniform without any depression and
detention losses over the entire studying area, the
effective rainfall intensity, ie, for Horton's
infiltration model becomes
^e = ^ ^c
(fQ fc ) e "kt (2.21)
in which i = uniform rainfall intensity (in/hr)
fQ= initial infiltration rate (in/hr)
fc= final infiltration rate (in/hr)
It can thus be seen that the effective rainfall
intensity will approach a constant falue of (i-fc),
when time reaches infinity.


33
2.5 Nondimensionalization of Governing Equations
Through a normalization process, the governing
equations may be nondimensionalized to a more desirable
form. The nondimensionalized equations have the
advantages of reducing variables, simplifying the
computational algorithms, and comparing results of
different investigators. It is a better way to
represent the values in the numerical analysis
(Woolhiser and Liggett, 1967).
In formulating an overland-flow problem, the
dimensional units involved are time, length, and force.
Because accumulation of flow at a basin outlet eventu-
ally reaches its equilibrium state, the equilibrium
quantities were used to normalize governing equations to
produce more signficant results in terms of ratios to
the equilibrium state.
In this study, the equilibrium outlet depth,
Ye, is employed as the characteristic depth scale; the
watershed length, L, is used as the characteristic
length scale. The equilibrium time of concentration,
Te, determined by the equilibrium effective rainfall
(that is, i fc), is used as the characteristic time
scale. According to Guo (1984A, 1984B), when m = 2
(Equation 2.16), YE and TE can be computed by


[LAC ( i fc)]1/2
(2.22)
T
E
- (i fc)
(2.23)
where L = length of watercourse (ft)
oi = watercourse characteristic factor
The watercourse characteristic factor,oi is equal to
(1.49 (S0)ly/2)/n, which is constant for a given
homogeneous watershed.
Utilizing the above two normalization
quantities, YE and tE' the following dimensionless
parameters can be defined:
X*
x
L
t*
q* =
q
n2 (i-fc)L (q*)2
2.22 YJ>/'2 (Y*)10//3
Hi
(2.24)
When the above dimensionless variables in Equation
(2.24) are incorporated into the continuity equation
(2.3) and the momentum equation (2.17), the following
equations result:


35
ay* + sq* -
at* d-x* -Le
aq* + at* = So 9 rq*2 + ax* L y* ^y* (1 . V 2L2 -M) 1 Fi2 (y*)2]
Fi2 L So
where Fj_, the rainfall Froude
(2.25)
(2.26)
(i-fc)/(g Yg)1/2, which is constant for a given
rainfall intensity and watershed.
Equation (2.26) can further be simplified to the
following form:
aq* + 3P_ = R
at* 9-x*
where P = [ SH+ D(y*)2]
y*
R = [K y* (1 -§£)]
so
Yp2 i
D = E 1
2L2 Fi2
Kj
(2.27)
The parameter, D, reflects the effects of the watershed
length, L, the equilibrium depth, YE' and the constant


36
F^. The parameter, K^, reflects the effects of the
slope, SQ, the watershed length, L, equilibrium depth,
YE, and the constant F^.
2.6 Initial and Boundary Conditions
In the treatment of partial differential
equations, the determination of a solution is dependent
upon some given functional relationships which are
broadly termed as initial/boundary conditions. Solu-
tions must satisfy the governing equations on the
interior and also meet the requirements on the
boundaries of the computational domain simultaneously.
A necessary initial condition for the solution
of overland flow is that all depths and discharges per
unit width along the plane, 0 ^ x ^ L, must be known at
the initial given time, tQ. A dry bed initial
condition for which all depths and discharges are zero
at the time rainfall commences (t = 0) has been studied
and recommended by several investigators (Morgali and
Linsley, 1965; Woolhiser and Liggett, 1967). A dry bed
initial condition was adopted in this study.
The necessary boundary conditions are overland-
flow parameters at locations of x = 0, and x = L as
shown in Figure 2.1. Because there is no flow entering
the flow plane through the upstream boundry (x = 0), it
is assumed that the flow through the upstream boundary


37
of the plane is zero. This means that all flow will
enter the flow plane as rainfall. All flow entering the
upstream station (1-2) is accounted for by monitoring
the outflow for this section. The depth at the
uppermost boundary was allowed to rise as storage
increase in the section by the storage routing (Morgali
and Linsley, 1965). It is assumed that there is a
straight-line water surface in this uppermost control
volume. The storage between stations 1 and 2,
therefore, is given by
Storage
1-2 (Yl + y2) (2.28)
In this study, the storage volume method (Morgali and
Linsley, 1965; Morgali, 1970) was used to calculate the
depth at the most upstream boundary.
Regarding the downstream boundary condition, the
study performed by Morgali and Linsley (1965) concluded
that the consideration of critical depth at the down-
stream boundary condition was not necessary.
Noncritical conditions gave good results for the
shallow-water flow on a subcritical slope. The backward
finite difference method (discussed in Section 3.3) gave
good results, and it was therefore used in this study.


38
CHAPTER III
NUMERICAL STUDY
3.1 Introduction
Many of the differential equations for
engineering problems cannot be easily solved by
analytical methods. Consequently, many of these
equations have been of little value in application in
the past. With the advent of the high-speed digital
computers, however, a number of very efficient numerical
methods have been developed in the application of
solution-solving to sets of partial differential
equations. These numerical solutions are always
associated with special boundary conditions which
describe the physical constraints for the whole system.
Therefore, a good numerical modeling is not only a
method of computation but also a proper mathematical
description of some physical constraints.


39
According to Roache (1976), computational fluid
dynamics is certainly not a pure theoretical analysis.
If the initial and the boundary constraints of the whole
system are properly defined, computational fluid
dynamics is closer to natural system representations
than are physical experimentations or theoretical
analyses.
There are several basic approaches to transform
a partial differential equation into finite-difference
form, such as Taylor series expansion, polynomial
fitting, integral methods, and control volume method.
In general, it is assumed that a linear relationship
between two adjacent points is applicable to estimate
the time rate of change or the spatial rate of change in
the differential equation. The detailed derivation and
discussion of finite difference formulation for partial
differential equations can be found in many references
such as Rosenberg (1969), Smith (1971), and Roache
(1976).
The mathematical formulation of most problems in
science and engineering involving rates of change with
respect to two or more independent variables, usually
representing time, length, or an angle, leads either to
a partial differential equation or to a set of such


40
equations. Special cases of the two-dimensional
second-order equation.
a + b + C + d ^ + e^
3x2 axay ay2 ^y
+ f$ + g = 0 (3.1)
(where a, b, c, d, e, f, and g may be functions of the
independent variables x and y, and of the dependent
variable pf) occur more frequently than any other because
they are often the mathematical form of one of the
conservation principles of physics.
The basic two-dimensional second-order equation
can be categorized as two-dimensional elliptic equa-
tions, parabolic and hyperbolic equations. The best
known among the two-dimensional elliptic equations are
the Poisson's equation (3.2) and Laplace's equation
(3.3).
2 2
+ lA + G = 0 (3.2)
^ 2 2
9x Qy
=0 (3.3)
. 2 2 v '
in which 0 is a physical variable and
G is the source or sink term.


41
These equations are generally associated with
equilibrium or steady-state problems (Bascoi, 1970;
Kersten, 1969). For example, the velocity potential for
the steady flow of incompressible nonviscous fluid
satisfies Laplace's equation. It is the mathematical
way of expressing the idea that the rate at which fluid
enters any given region is equal to the rate at which it
leaves it.
Problems involving time, t, as one independent
variable usually lead to parabolic or hyperbolic
equations. The simplest parabolic equation is derived
from the theory of heat conduction. It depicts the
temperature, u, at a distance, x units of length, from
one end of a thermally insulated bar after t seconds of
heat conduction. In such a problem, the temperatures at
the ends of a bar of length, 1, are often known for all
times. In other words, the boundary conditions are
known.
The simplest hyperbolic equation is the one-
dimensional wave equation. An example for this type of
equation is used to describe the transverse displace-
ment, u, at a distance, x, from one end of a vibrating
string of length 1 after a time, t.


42
3.2 Finite-Difference Method
Finite-difference methods are approximate in the
sense that derivatives at a point are estimated by-
difference quotients over a small interval. The
solutions are not approximated in the sense of being
crude estimates, however, the method generally gives
solutions that are either as accurate as the data
warrant or as accurate as. necessary for the purposes in
which the solutions are required.
From Taylors theory (Kersten, 1969; Basco,
1970), for a function, ft, and its derivatives are
single-valued; that is, finite and continuous function
of x exists, then
X (x + k) = X (x) + k j/' (x) +
1/2 k2 / (x) +
1/6 k3 X"' (x) + ... (3.4)
/ (x k) =/ (x) k/' (x) +
1/2 k2 X (x) -
1/6 k3 X" (x) + (3.5)
in which k = infinitesimal fluctuation of variable x;
X' = first derivative of function etc.


By neglecting all terms of higher power of fluctuation
of k, equations (3.4) and (3.5) give
43
/" (x) = d2J2f/dx2
= (1/k2) [/(x + k) 2/ (x) + / (x k))
(3.6)
(x) = d//dx
= (l/2k) [/ (x + k) / (x k) ] (3.7)
According to Figure 3.1, Equation (3.7) approximtes the
slope of the tangent at point P by the slope of the
chord, LR. This is the central-difference approxi-
mation. If we approximate the slope of the tangent at
point P by the slope of the chord, PR, or by the slope
of the chord, LP, the methods are called the "forward-
difference method" and the "backward-difference method",
respectively. The formulae are given by
/ (x) =(l/k)[/ (x + k) -/ (x)] (3.8)
(3.9)


Figure 3.1 Illustration of Finite-Difference
Method
}


45
3.3 Central-Difference Method
From Equation (3.8), it is clear that the
forward-difference methods are defined by points to the
right of the point P. On the other hand, the backward-
difference methods are defined by the points to the left
of the point P. Equation (3.8) can alternatively be
expressed by Equation (3.10):
A/0'(x)/Ax = J?f (x + Ax) ft (x) (3.10)
Equation (3.9) can alternatively be expressed by
equation (3.11):
A^(x)/Ax = ft (x) / (x Ax) (3.11)
Similarly, the central-difference methods are,
therefore, defined by points located symmetrically with
respect to point P. Equation (3.7) can be expressed by
Equation (3.12):
£$(x) /Ax = $ (x + Ax/2)
- / (X AX/2) (3.12)


46
Using the finite-difference approximations in
solving partial differential equations, some error can
arise in the computation. The discretization error is
one of the errors which has been studied (Forsythe and
Wasow, 1960). According to the study, the central-
difference method gives an averaged estimation from
Equation (3.12); that is, it will give the approximation
closest to the true value. The backward and the forward
differences will give relatively poorer estimates.
3.4 Implicit and Explicit Methods
The complete St. Venant equations cannot be
solved analytically, and hence numerical schemes may be
used to obtain discrete solutions. As an unsteady flow
in nature, the flow domain in terms of computational
grids is composed of time and space steps, as shown in
Figure 3.2. The governing equations are approximated by
the algebraic finite-difference equations, which are
solved numerically. These equations can be solved by
one of the following two methods. (1) Implicit
schemes: the unknown parameters are expressed
implicitly in simultaneous equations, which are solved
by an appropriate solution technique. (2) Explicit
scheme: the unknown parameters are expressed explicitly
as functions of known parameters and solved directly.


47
Miller (1971) has reported a complete treatise
on the subject of the explicit and implicit schemes.
The work of Amein and Fang (1969) has also provided
notable computational examples. A complete description
of the numerical solutions of the governing St. Venant
equations of unsteady flows is given by Isaacson et al.
(1956).
An implicit and explicit grid network is shown
in Figure 3.2 to contrast the two approaches. A formula
which expresses one unknown pivotal value (point P,
Figure 3.2) directly in terms of known pivotal values
(points L, M, and R, Figure 3.2) is an explicit scheme.
The method applied to the governing equations usually
result in linear algebraic equations from which the
unknown can be evaluated directly without iterative
computations.
A method in which the calculation of an unknown
pivotal value (point P, Figure 3.2) necessitates the
solution of a set of simultaneous equations (points L',
M, and R'f Figure 3.2) is an implicit method. The
method involves nonlinear algebraic finite difference
equations whereby the solution is obtained by iterative
computation.
There are advantages and disadvantages to both
the explicit and the implicit methods. Explicit methods


48
Figure 3.2 Network of Grid Points in
Computational Domain


49
are only conditionally stable, meaning that errors will
grow as the solution progresses and are a function of
the step sizes of time and distance. Trial-and-error
solutions using computers are necessary to establish
stability criteria. According to Ponce (1978), the
explicit methods often result in a smaller computational
time step (At) in order to obtain the desired accuracy.
This limit on the length of At usually makes explicit
methods less efficient than implicit methods. The
explicit methods are easy to program, however, and are
generally easy to handle.
Implicit methods were claimed to be
unconditionally stable and more computationally
efficient than are explicit methods. When applied to
flood routing in open channels, the implicit methods
have been found to be more stable, faster, and more
accurate than other finite-difference methods (Amien and
Fang, 1969; Price, 1974). The implicit methods are
considerably more difficult to program, however.
Brakensik et al. (1966) and Amien and Fang
(1969) concluded that the central-difference method
provides stable solutions for an unsteady overland
flow. The implicit central-finite difference with
noncritical downstream boundary condition is, therefore,
selected to perform numerical computations in this


50
study. More details can be found in the study of
Morgali and Linsley (1965).
3.5 Numerical Formulation of the Flow Equations
The continuity and momentum equations are dis-
cretized in the continuous x t plane, the dependent
variables are depth and discharge, and the independent
variables are time and distance. The numerical scheme
of the four-point implicit scheme (Preissman, 1960) is
adopted.
In this study, the finite-difference equations
for computer programming are derived from Equations
(2.24), (2.25), and (2.26). Using the general format of
o/i/^t + 3/2/9X = /3
(3.13)
(in which jzfj,
flow variables)
and the central-difference method, the general
finite-difference form can be shown as follows:


51
Xidj) J2^,3)
At " 2Ax
=/3 (i,j) (3.14)
in which (i,j) = value of ^ at grid point (i,j),
etc.
At = time increment (sec)
Ax = distance increment (ft)
i = station number along x-axis,
1 ^ i^ IE
3 = station number along y-axis,
1 < j < JE
IE = last station on x-axis, which is
the spatial axis
JE = last station on y-axis, which is
the time axis
Rearranging Equation (3.14), a more general form is
obtained:


52
/l (i/j) = /i (jfj-1)
+ (At/2Ax)[/2
- /2 (i + 1, j) ]
+ (At/2)[/3 (i, j )
+ /3 (i,j-1)]
(3.15)
Substituting the dynamic-wave routing parameters into
equation (3.15), Equation (2.3) for the depth and
Equation (2.17) for the flow rate per unit width at
point (i,j) can then be solved by the following
equations:
Y(i/j) = y(i,j-l) + (At/2Ax)[q(i-l,j) -q(i+l,j)]
+(At/2)[ie(i,j) + ie(i,j-l)
(3.16)
q(i/j) = q(i,j-l) +(At/2Ax)[(q(i -1,j)/y(i-l,j)
+ gy(i-l,j)2/2) (q(i + l, j)/
y (i + l,j) + gy (i + l, j)2/2)]
+(At/2)g [y(i,j) + y (i,j-l)]
(So-Sf)
(3.17)


53
in which y(i,j) = the value of depth at ith point in
space and jth point in time, etc.
Sj = Average friction slope for time j and
time (j-1)
Using the nondimensionalized formula, as described above
in Section 2.5, Equations (3.16) and (3.17) can further
be transformed into:
y*(i,j) y*(i,j~l) + C q*(i-l,j) -
q* (i+l,j)] + [ie*(i,j) +
ie* (i,j-1)] (3.18)
and
q*(i,j) = q*(i,j-l) + ^L. [P(i-l,j) -
2Ax*
P (i+l,j)] + £|!.[R(i,j) +R(i,j-1)]
(3.19)
in which y* (i,j) = dimensionless depth at grid (i,j),
etc.
The numerical formulations of boundary
conditions are different from Equations (3.18) and
(3,19). Assuming that the water course is dry shortly
before rain starts, as discussed in section 2.7, the
initial conditions are:


54
Y (i,l) =0 where 1 >$ i ><: IE, IE is the outlet
station
and
g (i, 1) = 0 where 1 >< i >< IE.
Because there is no flow entering the upstream
extremity, the discharges at the most upstream boundary
of the plane is always zero (q = 0). Using the storage
volume method (Morgali and Linsley, 1965), the depth at
the upstream boundary, Equation (3.12), can be
transformed to:
y* (1, j ) = 2 [ie*(l, j) (Ax*) q* (2, j) ]
Ax* e
-y* (2,j) (3.21)
As to the downstream boundary conditions, Morgali and
Linsley (1965) have tried both critical conditions, that
is, free fall outlet and noncritical conditions. From
their study, noncritical conditions gave good results
for the shallow-water flow on a subcritical slope.
Using Equations (2.25) and (2.27), the downstream
boundary conditions are thus formulated by a backward
finite-difference method as follows:


55
y*(IE,j) = y* (IE, j-1) + J£L [q* (IE-1, j) -
AX*
q*(IE,j)] + *|l[ie*(IE,j) +
ie* (IE,j-1)]
(3.22)
and
q* (JE, j ) = q* (IE, j-1) + [P(IE-1, j ) -
Ax*
P(IE,j)] +.^L [R (IE, j) +
R(IE/j-l)]
(3.23)
The stability and accuracy, plus the con-
sistency and convergence of the results of the numerical
analysis have been studied (Basco, 1974; Meadows,
1984) In summary, the key points to reach the best
approximation of a differential equation are to use
(1) the "best" finite difference scheme,
(2) the optimal At/Ax ratio to give optimal
accuracy, and
(3) the numerical scheme of stable conditions.
In this study, the central finite-difference
method is used because it gives the best approximate
results among the forward-, backward-, and
central-difference methods (Kersten, 1969; Rosenberg,
1969). The optimal At/Ax ratio was determined by the
sentivity test in this study for different cases.


56
In computational procedure, the previous water
profile is used as the initial guess for the
water-surface profile at current time step. The new
water-surface profile can be obtained by successive
iteration. During updating, the flow variables (depth
and discharge) for the profile at jth time step and the
flow variables at (j-l)th time step are kept unchanged.
In this study, criteria similar to those used by
Meadows et al. (1984) are used for judging the
convergence.
* (n+1) . .
y________(1*3
i) -
*'(n) .. i -i \
y (1*3+1)
y ^ (i,j + l)
^10
-3
(3.24)
and
* (n+1) . .
q________(1*3
i) q
(n)
(i, j + 1)
<10"3 (3.25)
q*(n) in which y*(n+1)(i,j+l) = the value of y*(i,j+l) after
(n+l)th iteration, etc.


57
3.6 Numerical Computation Procedures
The method of successive iteration was used in
the computation. After each sweep of the computational
domain, the flow variables (depths and discharges) are
updated by new values until Equations (3.24) and (3.25)
are satisfied for the entire domain. Figure 3.3 shows
the detailed logic computational flow chart which
consists of the following steps:
1. Calculate the values of normalizing variables, YE
and T-jr*.
£j
2. Initialize the values of flow variables, such as
depth and discharge.
3. Advance time increment by t and compute the
corresponding rainfall excess, j
4. Adopt the solutions of previous time step as the
initial guess.
5. Sweep the computations of depth and discharge for
the entire plane, that is, i = l to i = IE,
6. Update the new values for depth and discharge after
each sweep.
7. Check convergence criteria. If Equations (3.24) and
(3.25) are satisfied, solutions are achieved for the
current instant, then computation goes back to step
3. Otherwise, it goes back to step 5.


|Start I
58
\ Calculate tt normalizinc le values of: r variables
1 ,
Initialize flow variable depth and discharge
' f
I Compute rainfall excess.
' '
' No 1- Check if t* >3T-p 1 Ye s
* I S & 4
Ae i-ff 1 ^ = -ft
1
(Adopt solutions of previous time step as initial cruess
Sweep the computat discharge for th r -i ion of depth and e entire plane t

Update the new values of depth and discharge after each sweep
^
Check convergence: j No
within 10J I
No
end
Figure 3.3 Flow Chart of Computational Procedures


59
In this study, the rainfall duration, t^, is
assumed to be as long as three times the equilibrium
travel time^ TE, of the basin. It is believed that
solutions developed under this rainfall condition should
provide enough information for general cases. Recession
of .hydrograph starts when rain ceases. In the computer
program, the rainfall intensity was set to be.zero after
rain ceases, that is, i = 0 for t > t^. The
computation stops when the eniter catchment goes back to
dry bed condition again, that is, depth and discharge
vanishes everywhere.


60
CHAPTER IV
CASE STUDY
As mentioned in Section 2.3, the kinematic wave
approximation is a simplifiecation.of the complete St.
Venant equations which have been recognized to compose
the most adequate model for describing open channel-flow
conditions. As the kinematic wave formulation neglects
all terms in the momentum equation except the friction
slope and the channel-bottom slope, the formulation is
less complicated than the complete dynamic equations,
and therefore it requires less computational effort.
Consequently, it is desirable to use the kinematic model
whenever possible. The main purpose of this study was
to use the numerical analysis to determine under what
conditions kinematic wave solutions would give results
close enough to those obtained using the complete
dynamic equations. The analysis of dimensionless
governing equations (Section 3.5) indicates that two
parameters, the rainfall-kinematic-wave number, K^,
and the dynamic-wave number, D, should be the two most


61
important parameters in the numerical study of overland
flow. These two parameters are defined as:
S Y_
Vi o E
Ki -
F. L
l
(4.1)
in which SQ =
slope of plane
equilibrium water surface depth (ft)
(i-fc)/(g Yg)1/2, where g is the
gravitational constant (ft/sec^)
length of plane (ft)
and D
? 2
2LF.
l
(4.2)
4.1 Description of Study Aspect
A hypothetical system was devised to evaluate
the performance of kinematic-wave and dynamic-wave
models in the watershed system. The purpose was to
determine the criteria that would require that the
complete St. Venant equations be used. This
hypothetical system can be exemplified by a shallow
sheet flow of excess rain on a homogeneous catchment
such as a parking lot.
The watershed area used in this study is
described as follows:


62
(1) The plane has unit width.
(2) The length of the plane, L, was set to be equal to
500 ft.
(3) The surface roughness, n, is uniform throughout the
catchment. The Manning's value was set to be equal
to 0.03.
(4) The surface slope, S0, is uniform throughout the
catchment. The value was ranged from 0.0005 to
0.1.
(5) The soil at the plane surface is homogeneous in
nature, and infiltration parameters are:
fQ = initial infiltration rate = 4.5 in/hr
fc = final infiltration rate = 0.6 in/hr
k = infiltration decay constant = 0.0018/sec
The values are taken from the recommendations of
the Urban Drainage and Flood Control District
(1982) for soil type B as shown in Table 2.1.
(6) The rainfall intensity, i, is uniform in space and
also in time. The rainfall intensity value was set
to be equal to 12 in/hr with rainfall duration, td,
of 3Te. The time, TE, is the equilibrium
travel time of the basin as defined in Section 2.5.


63
4.2 Comparisons between Predictions of Dynamic-Wave and
Kinematic-Wave Models
Figures 4.1, 4.2, 4.3, and 4.4 show the effects
of the rainfall-kinematic-wave number, K^, on the peak
ratio of flow depth, Y/YE (as defined in Section
2.5). The kinematic-wave solutions of each case are
plotted as the solid line, and the dynamic-wave
solutions obtained from this study are plotted as the
dotted line.
Several points of interest can be observed from
these figures. First, Figure 4.1 (with(.= 1.11, and
Te = 1301 sec) shows that the dynamic-wave solution is
approximately 70% of the kinematic-wave solution when
the rainfall-kinematic-wave number, K^, is equal to 55
at the rainfall duration of 3TE. It also shows that
when Kj_ = 55, there exists a "point of separation"
where the kinematic-wave solution departs from the
dynamic-wave solution. In Figure 4.1, the point of
separation is located at approximately 0.5TE (point
A). After the point of separation, the dynamic effects
of the flow become important. The reason for this
behavior is that at the beginning of the rainfall, the
flow is shallow, and thus there is no significant
backwater effects on the plane (that is, the flow shows
kinematic behavior). Only when the rainfall duration is




65
long enough (t^ > 0.5TE, in this case), will the
backwater effects be pronounced.
The third point of interest is the different
values of "time to peak" (defined as time of
l
concentration in this study) between the kinematic and
the dynamic solutions. The time to peak value for
kinematic routing is approximately 1.13TE (point B,
Figure 4.1), and that for dynamic routing is
approximately 1.09TE (point C, Figure 4.2). By using
the Horton's infiltration model, the effective rainfall
intensity, ie, levels off to a constant value (i-fc)
when time reaches infinity. It is noted that the
dynamic-wave model has shorter time to peak than that of
the kinematic-wave model. In the practice of runoff
prediction, time to peak is important. Generally the
longer the time to peak is, the lower the peak flow
rate, Qp, will be. Therefore, kinematic-wave
approximation would possibly underestimate Qp when
backwater effects are significant.
The fluctuation of the dynamic solutions as
shown in Figure 4.1 may be attributed to (a) the limited
accuracy of numerical computation due to high-speed
digital computer and (b) the dynamic effects of flood
routing. The peak discharge might fluctuate over the
routing period. When the downstream backwater effect


66
increases, the upstream water surface depth increases.
As a result, the velocity and the inertial force of the
flow increase, which causes the flow to change from
dynamic to kinematic behavior. These natural resonant
effects, changing from dynamic to kinematic and back to
dynamic behaviors could cause the numerical solutions to
fluctuate.
Figure 4.2 (with o<.= 1.22, and TE = 1243 sec.)
shows that the dynamic-wave solution is about 80% of the
kinematic-wave solution when the rainfall-kinematic wave
number, K^, is equal to 60 at a rainfall duration of
3Te. In comparison with Figure 4.1, Figure 4.2 shows
that when the rainfall-kinematic-wave number, K^,
increases (physically it means that the slope of the
plane becomes steeper, the equilibrium water depth
becomes greater and the constant, Fj_, becomes
smaller), the point of separation starts at 0.45 TE
(point A), which is earlier than that in Figure 4.1.
This indicates that the duration of kinematic behavior
is shorter because the higher Kj_ value causes an
earlier backwater effect in the flow regime. In
addition, the difference of the time to peak between the
dynamic (point C, Figure 4.2) and the kinematic (point
B, Figure 4.2) models is now reduced as compared with
Figure 4.1. The reason is that when value


Figure 4.2 Illustration of Outlet Hydrograph B. L = 500 ft,
d. = 1.22, T = 1243 sec, K_^ = 60, i = 12 in/hr,
f = 3.5 in/hr, f = 0.6 in/hr, k = 0.0018
o c


68
increases, the whole flow regime behaves more similarly
to that of the kinematic-flow model rather than that of
the dynamic-flow model.
Figure 4.3 (witho£= 1.31, and TE = 1198 sec.)
shows that the dynamic-wave solution is approximately
88% of the kinematic-wave solution when the
rainfall-kinematic-wave number, K^, is equal to 65.
Figure 4.4 (withd = 1.4 and TE = 1,158 sec)
shows that the dynamic-wave solutionis approximately 90%
of the kinematic-wave solution when the rainfall-
kinematic-wave number, K^, is equal to 70. Figures
4.3 and 4.4 also show the same effects of kinematic
versus dynamic behavior as those shown in Figure 4.2
when comparing the point of separeation and the values
of time to peak of the two solutions as discussed above.
Figure 4.5 integrates the data shown in Figures
4.1 through Figure 4.4. The effects of various values
of rainfall-kinematic-wave number, K^, and the peak
ratio of water surface depth, Y/YE, is demonstrated.
As the term, Y/YE, approaches 1.0, the kinematic-wave
model is in good agreement with the complete dynamic-
flow model. Based on Figure 4.5, the following
conclusions can be drawn:


Y/Y
Figure 4.3 Illustration of Outlet Hydrograph C. L = 500 ft,
t^= 1.31, T = 1198 sec, K. = 65, i = 12 in/hr
E 1
E
fQ = 3.5 in/hr, f = 0.6 in/hr, k = 0.0018


Figure 4.4 Illustration of Outlet Hydrograph D. L = 500 ft,
<- 1.40, T = 1158 sec, K. = 70, i = 12 in/hr,
jb> l
fQ = 4.5 in/hr, f = 0.6 in/hr, k = 0.0018
-u
o


1.0
B
Kinematic-wave solution
Figure 4.5 Variation of the Peak Depth Ratio vs. Rainfall-Kinematic
Wave Number. L = 500 ft, t = 3T i = 12 in/hr,
hj
= 4.5 in/hr, f =0.6 in/hr
o c


72
(1) When the rainfall-kinemaitc-wave number increases,
the ratio of Y/Ys becomes higher and approaches 1=0
when the rainfall-kinematic wave" number is greater than
130. This indicates that the flow behaves in a manner
more similar to that of the kinematic model when the
value of Kj_ increases.
(2) When the rainfall-kinematic-wave number increases,
the point of separation where the dynamic-flow behavior
becomes significant starts at an earlier time. This
shows that the duration of kinematic behavior is
becoming less at the beginning of the rainfall when
increases,
(3) When the rainfall-kinematic-wave number increases,
the values of time to peak of the kinematic and the
dynamic routing are closer to each other. This
indicates that the flow regime behaves more similar to
that of the kinematic model when increases.
The relationship of the rainfall-kinematic-wave
number, Kj_, to the peak ratio are summarized and shown
in Figure 4.6. The figure shows that when the
rainfall-kinematic-wave number is greater than 180, the
kinematic wave reaches almost 100% agreement with the
dynamic-wave model.
It is noted that when the value of the rainfall-
kinematic parameter, Kj_, is smaller than 50, the


73
kinematic-wave solution is a poor approximation (Figure
4.1), and the complete dynamic equation should be used.
The maximum error, however, in the outflow hydrograph
with = 50 is in the order of 30%, and it decreases
rapidly as increases.
Figure 4.6 also shows that when Kj_ increases
to 150, the dynamic-wave model merges closely with the
kinematic model and will give results within 0.5% of the
dynamic-wave model (point A).
Another paramter, the dynamic-wave number, D,
which is defined by YE2/ 2L2 F^2, is also
investigated. The parameter is equivalent to x
(1/S0) x Ye/2L, which is equal to the (rainfall-
kinematic-wave number) x (1/SQ) x YE/2L. The study
showed that when the dynamic-wave number is 38, the
dynamic-wave solution is approximately 70% of the
kinematic-wave solution. When the dynamic-wave numbers
are 33, 27, and less than 6, the dynamic-wave solutions
are 80%, 90%, and 100%, respectively.
The results are depicted in Figure 4.7. It
shows that when the dynamic-wave number is less than 6,
the kinematic-wave approximation will give results
within 5% of the dynamic-wave solutions. The parameter,
D, can also be used as a complementary criterion to
check when the dynamic-wave routing can be approximated


Ym|
ye
1.0


Figure 4.7 Effects of Dynamic Waye Number on the Outlet Flow Depth
Ln


76
by the kinematic-wave routing. The results show that as
the dynamic-wave number, D, becomes smaller (that is,
the slope of the watershed increases), the wave routing
behaves closer to the kinematic model. This result
agrees with that of the rainfall-kinematic-wave number,
4.3 Evaluation of the applicability of the Kinematic-
Wave Model by using Available Laboratory Data
To provide more comparisons between the
predictions of the dynamic-wave and kinematic-wave
modeling and verfication of the criteria established in
this study, the laboratory data obtained by Bentner et
al. (1940) were used. The experiments of Bentner et al.
were performed on a watershed located in Tucson,
Arizona. It is approximately 24 ft long and 6 ft wide.
For the most part, the area is sparsely vegetated or has
bare soils, and the soil is rather tight. Runoff
hydrographs in Figures 4.8A, 4.8B, 4.9A, and 4.9B were
collected by Bentner et al. (1940).
Woolhiser and Liggett (1967) concluded that the
kinematic-wave equation is a good approximation for most
overland flow situations of hydrologic interest, because
Kq (defined in section 1.1) rarely falls below 10. In
general, Kq values smaller than 10 can be expected only
for short, smooth planes with mild slopes and high rates


77
of lateral inflow. These conditions are not often found
in rural areas but are possible in urban areas.
At equilibrium state, the normal depth (H0)
used in calculating the kinematic-wave number, Kq, in
Equation (1.2), is equal to the value of equilibrium
water depth (YE) used in calculating the rainfall-
kinematic-wav enumber, K^, in Equation (4.1). Table
4.1 illustrates the parameters used, SQ (slope of the
plane), n (roughness of the plane), and ot(watershed
characteristic constant) and the calculated values YE
(equilibrium water-surface depth), qE (equilibrium
flow rate) and Kq (kinematic-wave number), and K^
(rainfall-kinematic-wave number) of the four cases
represented in Figures 4.8A, 4.8B, 4.9A, and 4.9B. It
is noted that for all laboratory cases the values of Kq,
which are the criteria used by Woolhiser and Liggett
(1967) are greater than 10 (within 10% of the kinematic
approximation). The values of K^, which are the
criteria established in this study, are greater than 70
(within 10% of the kinematic approximation). In
addition, for a highly kinematic wave flow, the two
values converge. All laboratory data show kinematic-
flow behavior. Therefore, this study also concluded
that the kinematic-wave equation is a good approximation
for most, but not all, overland flows.


Figure 4.8A Hydrograph of Laboratory Data, Set 1
1 = 24 ft, S = .089, n = 0.171
o


Figure 4.8B Hydrograph of Laboratory Data, Set 1,
L = 24 ft, SQ = .089, n = 0.136


Fiaure 4.9A Hydrograph of Laboratory Data, Set 2
L = 24 ft, SQ = 0.02, n = 0,039


81
Z
Figure 4.9B Hydrograph of Laboratory Data, Set 2,
L = 24 ft, Sq = 0.02, n = 0.045


TABLE 4.1
Illustration of Laboratory Watershed Analysis
Based on Four Cases in Figure 4.8A, 4.8B, 4.9A, and 4.9B
Parameters Case 1 Case 2 Case 3 Case 4
Slope, Sq 0.089 0.089 0.02 0.02
Roughness, n 0.171 0.136 0.039 0.045
Watershed Characteristic Constant, ot. 2.6 3.7 5.4 4.7
Equilibrium Depth, Ye (ft) 0.0202 0.0205 0.0170 0.0185
Equilibrium flow rate* qE (cfs/ft) 0.00108 0.00138 0.00156 0.00161
Equilibrium flow rate**, qE (cfs/ft) 0.00106 0.00138 0.00156 0.00161
Kinematic-Wave Number, Kg 25012 15257 1832 2042
Rainfall-Kinematic- Wave Number, Kj_ 25050 15268 1840 2057
Remark Kinematic Kinematic Kinematic Kinematic
laboratory Data
**Numerical Solution


83
CHAPTER V
SUMMARY AND CONCLUSIONS
5.1 Conclusions
The use of the complete St. Venant equation with
Hortons infiltration formula for unsteady overland flow
wave routing is a sophisticated and complicated
wave-routing model. The use of such models has not been
common among practicing engineers due to their
complexity and the requirements for large amounts of
computer time.
1. The rainfall-kinematic-wave number, K^,
which consists of the slope, SQ, the equilibrium flow
depth at outlet, YE, and the rainfall Froude number,.
Fj_, can serve as a criterion for determining when
kinematic-wave assumptions are adequate and acceptable.
When the rainfall-kinematic-wave number is greater than
180, the kinematic-wave approximation will give results
within 5% of the dynamic-wave approximation.
2. When the rainfall-kinematic-wave number
increases, the point of separation starts at an earlier


84
time. This shows that the duration of kinematic
behavior becomes shorter at the beginning of the
rainfall when Kj_ is larger. Furthermore, the values
of time to peak of the kinematic and the dynamic routing
are closer to each other, when the rainfall-kinematic-
wave number, K^, increases. This indicates that the
flow regime behaves closer to the kinematic model when
increases.
3. Another parameter, the dynamic-wave number,
D, which consists of the equilibrium flow depth at
outlet, Ye, watershed length, L, and the rainfall
Froude number, F^, serves as another checking
parameter for determining when kinematic-wave
assumptions are adequate and acceptable. It was
concluded that when the dynamic-wave number is less than
6, the kinematic-wave approximation will give results
within 5% of the dynamic-wave approximation.
4. The testing of the kinemaitc-wave model by
laboratory data shows good agreement between this study
and that of Woolhiser and Liggett (1967). It is
concluded that the kinematic-wave equation is a good
approximation for most, but not all, of the real world.
overland flows.


85
5.2 Recommendations for Further Study
This study incorporates the Horton infiltration
model into the complete dynamic St. Venant equations;
that model was not included in the study of Woolhiser
and Liggett (1967). Therefore, the effect of rainfall,
rather than the effect of runoff, is modeled. The
guideline parameter, (rainfall-kinematic-wave
number), can readily be calculated for a given
watershed. The parameter used by Woolhiser and Liggett
(1967) Kg (kinematic-waver number) cannot be
calculated easily for a given watershed. This study
offers a much better and more complete dynamic-wave
modeling.
Further study could add several features to the
routing model. The comparison of kinematic-wave routing
to full dynamic-wave routing should be expanded. For
example, tests should be conducted to consider the plane
of various lengths and various roughness, etc., to
determine if changes in these parameters would alter the
criteria developed.
Another recommendation is for additional
applications of the model. As different situations are
encountered, obvious needs for improvements will become
apparent. More study can be done to refine the
numerical modeling to serve as a reliable engineering
tool.


86
REFERENCES
Amein, M., and Fang, C. S.
"Implicit Flood Routing in Natural Channels"
J. Hyd. Div., ASCE, HY 12, p. 2481-2500 (1970).
Basco, D. R.
"Introduction to Numerical Method"
Verification of Mathematical and Physical Models in
Hydraulic Engineering, ASCE, p. 280-302 (1978).
Beutner, E. L., Gaebe, R. R., and Horton, R. E.
"Sprinkled-plat Runoff and Infiltration-experiments on
Arizona Desert Soils"
Trans. Amer. Geophys. Union, p. 550-558 (1940).
Behlke, C. E.
"The Mechanics of Overland Flow"
Ph. D. Dissertation, Stanford University, Stanford,
California (1957).
Brakensik, D. L., et al.
"Numerical Technique for Small Watershed Flood Routing"
U.S. Agric. Res. Serv. Publ., ARS 41:113 (1966).
Brater, E. F., and King, H. W.
"Handbook of Hydraulics"
McGraw-Hill Book Company, New York (1960).
Chen, C. W.
"Computer Simulation of Urban Storm Water Runoff"
J. Hyd. Div., ASCE Hy 2, p. 289-301 (1971).
Chen, Y. H.
"Mathematical Modeling of Water and Sediment Routing in
Natural Channel," Ph. D. Diss., Department of Civil
Engineering, Colorado State University, Fort Collins,
Colorado, June (1979).
Chow, V. T.
"Open Channel Hydraulics"
McGraw-Hill Book Company, New York (1959).
Drainage Criteria Manual
Published by Urban Drainage and Flood Control District,
Denver, Colorado (1982) .


87
Eagleson, P. S.
"Dynamic Hydrology," p. 344-346,
McGraw-Hill Book Company, New York (1970).
Forsthye, G. E., and Wasaw, W. R.
"Finite Difference Methods for Partial Differential
Equations"
John Wiley and Sons, Inc., New York (1960).
Fread, D.L.
"Technique for Implicit Dynamic Routing in Rivers with
Tributaries"
Water Res. Research, 9 (4), p. 918-926 (1973).
Gburck, W. J., and Overton, D. E.
"Subcritical Kinematic Flow in a Stable Stream"
J. Hyd. Div. ASCE, HY9, p. 1433-1446 (1973).
Grace, R. A., and Eagleson, P. S.
"The Modeling of Overland Flow"
Water Res. Research, p^ 393-403 (1966) .
Greebm, W. H., and Ampt, G. A.
"Studies on Soil Physics: I. Flow of Air and Water
through Soils"
Journal of Agr. Sci., v. 4, p. 1-24 (1911).
Guo, C. Y.
"Effect of Infiltration on Runoff Hydrograph"
Water today and tomorrow, preoceedings, ASCE, Irregation
and Drainage Conference in Flagstaff, Arizona, July, p.
485-492 (1984A).
GUO, C. Y.
"Fluid travel time"
Water for Resources Development, Proceedings of ASCE
Hydraulic Conference in Coeur d'Alene, Idaho, p.
380-384, August (1984B) .
Henderson, F.M.
"Open Channel Flow"
Macmillan Publishing Co., Inc., New York (1966).
Henderson, F. M., and Wooding, R. A. :
"Overland Flow and Groundwater Flow from a Steady
Rainfall of Finite Duration"
Geophysical Research, v. 69, No. 8, p. 1531-1540 (1964).


88
Holton, H. N., and Minshall, N. E.
"Plot Samples of Watershed Hydrology"
U.S. Agri. Res. Serv. Publ., ARS 41:33 (1968).
Horton, R. E.
"The interpretation and application of runoff plot
experiments with reference to soil erosion problems"
Soil Sci. Soc. Amer. v. 3, p. 340-349 (1938).
Horton, R. E.
"Approach toward a Physical Intterpretation of
Infiltration Capacity"
Proc. Soil Sci. Soc. Amer., v. 5, p. 399-441 (1939).
Isaacson, E., et al.
"Numerical Solution of Flood Prediction and River
Regulation Problems: Report III"
Report No. IMN-NYU235, Inst, of Mate, and Sciences, New
York University, New York, New York (1956).
Izzard, C. F.
"Hydraulics of Runoff from Developed Surfaces"
Proceeding, 26th Annual Meeting of the Highway Research
Board, V. 26, p. 129-146 (1946).
Kersten,l R. D.
"Engineering Differential Systems"
McGraw-Hill Book Company, New York (1969).
Kuelegan, G. H.
"Spatially variable discharge over a sloping plane"
Trans. Amer. Geophysical Union, Part VI, p. 956-958
(1944) .
Labadie, J. W., Morrow, D. M., and Lazaro, R. C.
"Urban stormwater control package for automated
real-time systems"
Colorado State University Report C 6174, prepared for
U.S. Depart, of Interior, Office of Water Research and
Technology, Decembrer (1978).
Lazaro, T. R.
"Urban Hydrology"
Ann Arbor Science Publishing Inc., Ann Arbor, Michigan
(1979).
Liggett, J. A.
"Unsteady Open Channel Flow with Lateral Inflow"
Ph. D. Dissertation, Stanford University, Stanford,
California (1959).


89
Liggett, J. A., and Woolhiser, D. A.
"Difference Solution of the Shallow-Water Equation"
Proc. ASCE J. Mech. Div. 9B (EM2), p. 39-71 (April,
1967) .
Lighthill, M. H., and Whithan, G. B.
"On Kinematic Waves, I, Flood Movement in Long rivers"
Proc. Royal Soc. of London, Ser. A., v. 229, p. 281-316
(1955).
Lindsley, R. A.
"Hydrology for Engineering"
McGraw-Hill Book Company, New York (1958).
Meadows, M. E.
"Finite element simulation of overland flow"
Advances in irregation and drainage, surviving external
pressure, AASCE, p. 466-475 (1983) .
Miller, W. A., Jr.
"Numerical solution of the equation for unsteady open
channel flows"
School of Civil Engineering, Georgia Institute of
Technology, Atlanta, GA (1971).
Morgali, J. R.
"Laminar and Turbulent Overland Flow Hydrographs"
J. Hyd. Div. ASCE, HY2, p. 441-460 (1970).
Morgali, J. R., and Linsley, R.
"Computer Analysis of Overland Flow"
J. Hyd. Div. ASCE, HY3, p. 81-100 (1965).
Musgrave, G. W.
"How Much of the Rain Enters the Soil"
Yearbook of Agriculture, p. 151-159 (1955).
Overton, D. E., and Meadow, M. E.
"Stormwater Modeling"
Academic Press, New York. (1976).
Philip, J. R.
"Theory of Infiltration"
Advances in Hydroscience, Academic Press, New York, v.
5, p. 215-296 (1969).
Ponce, V. M., Li, R. M., and Simons, D. B.
"Applicability of kinematic and diffusion models"
ASCE Journal of Hydraulics Division, v. 104, no. HY3, p.
353-363 (1978).


90
Preissman A.
"Propagation des intumescences dans les cannaux et
rivieres"
1st Congress de 1'Assoc. Francaise de Calcul, Grenoble,
p. 443-442 (1960).
Price, R. K.
"Comparison of four numerical metnods for flood routing"
Journal of the Hydraulic Division, ASCE, v. 100, no,.
HY7, p. 879-899 (1974).
Richey, R. P.
"The Fundamental Hydraulics of Overland Flow"
Ph D. Dissertation, Stanford University, Stanford,
California (1954).
Roache, P. J.
"Computer Fluid Dynamics"
Hermosa Publisher Inc., New Mexico (1976).
Schreiber, D. L., et al.
"Subcritical Kinematic Flow in a Stable Stream"
Journal of Hydraulic Division, ASCE, HY3, p. 429-446
(1972) .
Scheaffer, J. R., et al.
"Urban Storm Drainage Management"
Marcel Dekker Inc., New York (1976).
Simons, D. B., Chen, Y. H., and Saez-Benito, J. M.
"A mathematical model of Pools 5 through 8 in the Upper
Mississippi River System"
Colorado State University Report CER 79-89, DBS-YHX-JMS
21, prepared for U.S. Army Corps of Engineers, St. Paul
District, December (1979).
Smith, G. D.
"Numerical Solution of Partial Differential Equations"
Brunei College of Advanced Technology, London (1980).
Smith, R. E.
"Documentation for program KINEROS"
Unpublished report, Colorado State University, Fort
Collins, Colorado (1980).
Tucci, C. E. M.
"Hydraulic and water quality model for a river network"
Ph. D. Dissertation, Department of Civil Engineering,
Colorado State University, Fort Collins, Colorado
(1978) .


von Rosenberg, D. U.
"Methods for the Numerical Solution of Partial
Differential Equations"
Elsevier Publishing Company, Inc., New York (1969).
Woo, D. C., and Brater, E. F.
"Spatially Varied Flow from Controlled Rainfall"
Proc. ASCE, V. 88, p. 31-56 (1962)
Wooding, R. A.
"Hydraulic Model for the Catchment"
J. Hydrology, 3, p. 254-267 (1965).
Woolhiser, D. A., and Liggett, J. A.
"Unsteady, One-Dimensional Flow Over a Plane The
Rising Hydrograph"
Water Res. Research, 3(3), p. 753-771 (1967).
Yen, B. C.
"Methodologies for Flow Prediction in Urban Storm
Drainage System"
Research Report No. 72, University of Illinois,
Champaign-Urbana, Water Res. Center, January (1973).
Yen, B. C.
"Storm Sewer System Design"
University of Illinois, Champaign-Urbana, Illinois
(1978) .
Yevjevich, V. and Barnes, A. H.
"Flood Routing through Storm Drains"
Colorado State University, Hydrology Papers, p. 43-46
(1970).


Full Text

PAGE 1

AND KINEMATICS OF OVERLAND FLOWS by Shiu-Ming Huang ?' B.S., National Taiwan University, 1971 M.S., Pennsylvania State University, 1973 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of civil Engineering 1985

PAGE 2

This thesis for the Master of Science degree by Shiu-Ming Huang has been approved for the Department of Civil Engineering by 7

PAGE 3

Huang, Shiu-Ming (M.S., Civil Engineering) Dynamics and Kinematics of Overland Flows iii Thesis directed by Assistant Professor James c. Y. Guo The dynamic-wave model for an overland flow was assessed using the implicit central finite-difference method. Flow equations were derived from the continuity and momentum principles. This study attempts to delineate the limitations of kinematic-wave approximation by_performing a series of numerical studies using the flow equations. The nondimensional dynamic-wave model was combined with Horton's infiltration formula and numerically solved to establish a model representing two-dimensional overland flow over a pervious surface. The comparisons between kinematic-wave solutions and dynamic-wave solutions provide general guidelines anq further establish new criteria for determining when kinematic-wave assumptions are valid. Results from this stUdy indicate that the overland runoff is mainly kinematic when the rainfallkinematic number is greater than 180 or when the dynamic-wave number is less than 6.

PAGE 4

iv This abstract is approved as to form and content. Faculty member in charge of thesis

PAGE 5

v ACKNOWLEDGMENTS I would like to express my deep gratitude to Dr. James c. Y. Guo for his advice, dedication, and continuous guidance throughout the research and the duration of the thesis work. Thanks are also due to other faculty members of the University of Colorado, including Dr. David w. Hubly, Dr. William c. Hughes, and Dr. John Liu, who introduced me to the aspects of water resources engineering. Special thanks go to my husband, Chi-I, and my sons, David and Samuel, for their encouragement, support, and understanding. This research was supported in part by the Education Program of the Martin Marietta Denver Aerospace.

PAGE 6

vi TABLE OF CONTENTS CHAPTER PAGE I. II. III. INTRODUCTION . 0 1 1.1 Description of Problem 1 1.2 Goals and Objective Cl 0 6 1.3 Literature Review 7 ANALYTICAL MODEL 11 2.1 Introduction 11 2.2 Continuity Equation 12 2.3 Momentum Equation 16 2.4 Horton's Infiltration Model 25 2. 5 Nondimensionalj_zation of Governing Equations 33 2.6 Initial and Boundary Conditions 36 NUMERICAL STUDY . 3.1 Introduction 3.2 Finite Difference Method 3.3 Central Difference Method 3.4 3.5 Implicit and Explicit Methods Numerical Formulation of the Flow Equations 38 '3a 42 45 46 50. 3.6 Numerical Computation Procedures 57

PAGE 7

CHAPTER IV. CASE STUDY o 4o1 Description of Study Aspects 0 4o2 Comparisons between Predictions of Dynamic-Wave Model and vii-PAGE 60 61 KinematicWave Nodel o 63 v. 4o3 Evaluation of Kinematic-Wave Model Applicability with Available Laboratory Data SUMMARY AND CONCLUSIONS 5.1 Concltisions . . . 76 83 33 5.2 Recommendations for Further Study o o 35 REFERENCES . . 86 APPENDIX I NOTATIONS 92

PAGE 8

viii LIST OF TABLES TABLE PAGE 2.1 Recommended Horton's Equation Parameters 31 4.1 Illustration of Laboratory Watershed Analysis 82

PAGE 9

ix LIST OF FIGURES FIGURE PAGE 2.1 Illustration of Overland Flow Plane 13 2.2 Control Volume to Derive Continuity Equation . . . . . 15 2.3 Control Volume to Derive Momentum Equation . . . . . 18 2. 4 Schematic Representation of Horton's Equation . . . . 30 3.1 Illustration of Finite-Difference Method 44 3.2 Network of Grid Points in Computational Domain . . . . 48 3.3 Flow Chart of Computational 58 4.1 Illustration of Outlet Hydrograph A . . 64 4.2 Illustration of outlet Hydrograph B 67 4.3 Illustration of Outlet Hydrograph c 69 4.4 Illustration of Outlet Hydrograph D . . 70 4.5 Variation of the Peak Depth Ratio versus Rainfall-Kinematic-wave Number 71 4.6 Effects of Rainfall-Kinematic-Wave Number on the outlet Flow Depth . . 74 4.7 Effects of Dynamic-Wave Number on the Outlet Flow Depth . . . 75 4.8A Hydrograph of Laboratory Data, Set 1 . . 78 4.8B Hydrograph of Laboratory Data, Set 1 79 4.9A Hydrograph of Laboratory Data, Set 2 80 4.9B Hydrograph of Laboratory Data, Set 2 . . 81

PAGE 10

CHAPTER I INTRODUCTION 1.1 Description of Problem The surface runoff produced by rainfall is hydraulically classified as unsteady, nonuniform flow, meaning that the velocity and depth of flow vary with respect to both time and space. In a catchment, the lowest point, or the outlet of the basin,accumulates the runoff from the entire drainage area. This results in the nonuniformity of the flow movements. In addition, the time variations in rainfall excess also cause unsteady and nonuniform flow mechanism. With such a complexity in nature, it is not possible to use the assumptions of steady uniform channel flow, inherent in Manning's equation or Chezy's equation, to properly compute the runoff. In hydraulic engineering practice, the prediction of runoff is very fundamental for water-resources planning.

PAGE 11

2 The theoretical hydrodynamic equations governing overland flow are generally attributed to Barre de st. Venant and were formulated in the late 19th century. Due to their mathematical difficulty, however, the St. Venant equations have not yet been solved analytically. With the advent of high-speed digital computers, numerical schemes have been developed and used to obtain discrete solutions for complicated open-channel flow situations. Prior to the use of high-speed computers, various simplified forms of the st. Venant equations were used to approximate different hydrological problems. The kinematic-wave approximation is the most common simplification of the st. Venant equation. It assumes that the friction slope, Sf1is balanced by the channel-bottom slope, s0 but it ignores the pressure and inertial terms. Lighhill and Whitham (1955) have found that the flood-wave propagation in overland flow is mainly kinematic; Wooding (1965) applied the kinematic-wave approximation to solve the partialdifferential governing equations describing the flow for a v-shaped chatchment flow. His solutions represented complete runoff hydrographs and the time of concentration. The general solution for the kinematic-wave approximation of overland flow is similar to Manning's

PAGE 12

equation, in which, q, the flow rate per unit width, is a power function of flow depth, y: (1.1) in watershed characteristic constant, for a uniform flow, o1... = (1/n) 1. 49 (S0 ) 112 3 m varies betweem 5/3 and 3, depending upon where the flow regime is turbulent or laminar. Wooding (1965) also concluded that the kinematic-wave method was applicable to overland flow when the Froude number at outlet is less than 2. In his approach, however, the infiltartion rate was assumed to be negligible or constant throughout the runoff period. As mentioned above, the kinematic-wave theory is the simpliest form of the complete dyanamic-wave model. With the simplification, kinematic-wave theory requires the least computational efforts in predicting an overland flow runoff hydrograph. Neglecting the dynamic flow factors, including inertial force and backwater effects, however, the kinematic-wave theory has its limitation. Woolhiser and Liggett (1967) determined the limitations of kinematic-wave theory on the basis of a

PAGE 13

series of computer simulation studies of the dynamic-wave model. They concluded that for a given uniform rainfall excess, ie, the applicability of kinematicwave theory may be best judged by the kinematic-wave number, Kq, which is defined as: in which kq -so = Lo = Ho = S L 0 0 HF_2 0 d slope length normal of plane of plane depth at rate, q ( 1. 2) (ft) equilibrium runoff flow F0 = Froude nubmer at equilibrium runoff flow -vo rate, q and depth, H0 that F0 1/2 JgHo) (where V0 = q/H0 ) 4 The study shows that the kinematic wave approximation is within 10% of the kinematic approximation if Kq is greater than 10. These two significant advances in the develop-ment of kinematic-wave approximation were based upon the assumption of constant infiltration loss and uniform rainfall intensity. In reality, the infiltration rate of a watershed has a major influence on the volumetric rate of runoff generated from a watershed under a given

PAGE 14

5 rainfall event. A constant infiltration rate occurs only when the catchment is saturated before the rainfall begins. Therefore, the assumption of the constant infiltration rate established the limitation of Wooding's conclusion. Even with a uniform design rainfall, the rainfall excess is no longer uniform when the time variation of infiltration rate is considered. In general, the mathematical model to describe the decay of infiltration rate with respect to time is empirically derived from the observation of infiltrometer tests. There are several infiltration models available, such as those of Horton (1939), Holtan (1968), Philip (1969), Green (1911), and Soil Conservation Service (SCS). Guo (l984A, 1984B) combined the kinematic-wave approximation with Horton's infiltration formula to analytically solve the propagation of a two-dimensional overland flow on a pervious surface. He reformulated the equation of time of concentration and concluded that Wooding's assumption of constant infiltration rate is a special case of the general situation. Nothing, however, was advanced regarding the establishment of kinematic-wave approximation limitations under a nonuniform rainfall excess due to infiltration decay characteristics.

PAGE 15

6 This study attempts to determine the limitations of kinematic-wave approximation by performing a numerical study which combines the dynamic-wave model with Horton's infiltration formula to numerically solve the propagation of two-dimensional overland flow on a pervious surface. The results from comparisons between kinematic-wave solutions and complete dyanamic-wave solutions can provide general guidelines, and can further establish new criteria for determining when kinematic-wave assumptions are adequate and acceptable. 1.2 Goals and Objective The purpose of this study was to develop a computer simulation model to solve the dynamic-wave equations representing overland flow resulting from a uniform rainfall and a variable infiltration rate. The discrete computational solutions were then compared with the kinematic-wave solutions developed by Guo (1984A and 1984B) in order to delineate the limitations of the applicability of kinematic-wave solutions. The kinematic-wave method is an approach which has been proven to be a very useful tool in storm-water modeling. Consideration is given to flow unsteadiness without solving the complete st. Venant equation. Ignoring backwater effects, however, imposes limitations

PAGE 16

7 upon its applicability. Kinematic waves occur when the dynamic terms in the momentum equation are negligible. Under these circumstances, flow rate can be expressed as a function of flow depth at all distance (x) and time (t) The basic objectives for this study are the following. (1) To develop a numerical routing model which made use of the complete st. Venant equation which includes Horton's infiltration formula. An implicit formulation is used to avoid the instability associated with explicit schemes. (2) To develop criteria for determining whether the kinematic-wave method is applicable. As kinematicwave method is less complex and therefore;:. easier to model, it is desireable to use this method when it gives acceptable results. Test runs are made for various conditions on a hypothetical system, and results of the two models are compared to determine under what conditions kinematic-wave assumptions are acceptable. 1.3 Literature Review Extensive literature reviews regarding overland flow have been discussed by Linsley et al. (1982) and Chow (1959). Only a few references pertaining to this

PAGE 17

study were found, for instance, Richey (1954), Behlke (1957), Liggett (1959), Izzard (1946), Woo and Brater (1962), Isaacson et al. (1956), Morgali and Linsley (1965), Woolhiser and Liggett (1967), Morgali (1970), Schreiber and Bender (1972), Gburek and overton (1973), and Guo (1984A and 1984B). 8 Isaacson et al (1956) used the finite-difference mathematical scheme on a digital computer for synthesizing flows. Morgali and Linsley (1965) then performed a numerical analysis, using the computer for overland flow. They concluded that an overland-flow hydrograph can be synthesized by using an appropriate mathematical model and numerical methods. The influences of various physical parameters (such as slopes, roughness, plane length, and rainfall intensities) on the overland-flow hydrograph were numerically tested. They proposed a semi-theoretical formula for estimating the time of concentration under a uniform rainfall excess. Woolhiser and Liggett (1967) used the finitedifference intergration to solve the rising hydrograph of the unsteady, overland flow. A dry channel was used as an initial condition; the upstream and downstream boundary conditions were, respectively, zero inflow and critical flow depth (or no condition for supercritical flow). The study concluded that when compared with

PAGE 18

previous numerical, analytical, and experimental results, the results show that in general there is no unique dimensionless rising hydrograph for overland flow, but that for most hydrologically significant cases, the kinematic-wave solution gives very accurate results. 9 Morgali (1970). did_ a follow-up study on laminar and turbulent overland hydrographs. In setting up the solution of finite-difference schemes, the approximation of the resistance to flow, that is, the friction slope, was found to be an important factor. Schriber and Bender (1972) developed a mathematical model that can be used to simulate an overland-flow hydrograph for a given rain strom and watershed. A dimensional analysis is used as an aid to determine predictive relationships for the resistance of an overland-flow plane. Giburek and Overton (1973) studied subcritical kinematic flow in a stable stream. An evaluation of the shallow-water equations in the lower 6.9 miles of Mahantago Creek, Pennsylvania, showed that subcritical kinematic conditions exist. Guo's studies (1984A and 1984B) on effects of infiltration on a runoff hydrograph established the analytical solutions for a kinematic overland flow on a pervious surface by combining kinematic-wave theory with

PAGE 19

Horton's infiltration model. The accuracy of solution is confirmed by numerical computations based on an implicit finite-differnece scheme and experimental data. Results indicate that the flood travel time for overland flow on a pervious surface is not a constant time interval and that the peak flood,flow rate occurs when rain ceases. 10 Recent trends in model developments of the rainfall/runoff relationship recommended by the Urban Drainage and Flood Control District (1982) indicate a tremendous need for an efficient, but more complete, runoff model. Such models can be valuable tools for the planning, design, and management of practical problems due to rainfall. It is the objective of this study to incorporate a decay-infiltration model to the unsteady dynamic wave routing and to further investigate the limitations of the kinematic-wave approximation.

PAGE 20

11 CHAPTER II ANALYTICAL MODEL 2.1 Introduction The governing equations of a dynamic-wave model for an overland flow are derived from the continuity and momentum principles. The equations were first applied to the overland-flow problem by Keulegan (1944), and later Morgali and Linsley (1965) did a study on computer analysis of overland flow. The schematic illustration of two-dimensional flow is shown in Figure 2.1. The accumulation of runoff on a unit-width plane with a slope, S0 and length, L, results from rainfall of constant intensity, i. The basic assumptions made in the flow equations are as follows: 1. The plane has a unit width. 2. The hydraulic radius, Rp, is equal to the depth of flow, Y.

PAGE 21

12 3. The momentum coefficient is equal to 1, and the energy correction coefficient associated with velocity head is equal to unity. 4. The distribution is hydrostatic at any vertical plane in the flow. 5. -....The slope of the plane is s.o small that cos s0 = 1. 6. The momentum of incoming rainfall in the direction of flow is negligible. 7. The friction loss in unsteady nonuniform, open-channel flow is equal to friction loss th steady uniform flow. a. Surfaca slope and roughness are uniform throughout the catchment. 9. Rainfall intensity is uniform in space. 10. Soil coverage is homogeneous in the entire catchment so that infiltration rate varies only with respect to time. 2.2 Continuity Equation The governing of continuity for unsteady flow may be derived by the conservation of mass in a control volume, that is, a small space between the channel sections. Considering the concept of control

PAGE 22

i; rainfall intensity x=O upstream boundary ---;)1 X 11 ! i = rainfall intensity y = depth slope of plane s = 0 v = f = t velocity infiltration loss X= downst+eam boundary Figure 2.1 Illustration of overland Flow 13

PAGE 23

14 volume between stations 1 and 2 as shown in Figure 2.2, the principle of conservation of mass states that: inflow -outflow = rate of change in volume stored Therefore, the volume increased from time t to time (t can be represented by the area, ABCD, which is equivalent to: ( 2 .1) in which qX = average inflow between time t and time (t + At) (cfs/ft) q0 = average outflow between time t and time (t+At) (cfsjft) At = time increment (sec) Ax = distance increment (ft) AY = average water surface depth increment (ft) ie = average rainfall intensity excess (ftjsec), which is the difference between the rainfall intensity and infiltartion loss '

PAGE 24

5 = 0 i = Q) = qi = qo = OBC = OAD = ft = 1 slope of plane rainfall intensi.ty station number average inflow between i L L L @ timet and average outflow between time t and water surface at time t water surface at time (.t+At) infiltration loss lS t+At t t i.me Ct +At l_ time (t+At) Figure 2.2 Control Volume to Derive Continuity Equation

PAGE 25

16 Dividing Equation (2.1) by At AX yields ( 2. 2) Considering infinitesimal spacing in distance and time, equation (2.2) becomes ( 2. 3) Equation (2.3) is the continuity equation for a two-dimensional overland flow. When rainfall intensity is neglected, equation (2.3) reverts to conventional form (Chow, 1959). Jl9. + ... il = 0 at 2.3 Momentum Equation (2.4) The momentum equation for a two-dimensional flow may be derived from Newton's second law of motion. According to that law, the change of momentum in the waterbody is equal to the summation of all of the exter-' nal forces that are acting on the body. The external forces acting on a control volume are shown on the flow

PAGE 26

17 'l. .,... 2 prorL.e 1n .3. Thesa external forces include pressure gradient, body force, and the rasistive force. The :most simple fcr:m of N::rwton s la;,.; of mo1cion sta-te3 that F Til a = DU !\1-Dt F = sum of all externa.l :forces in a given direction (in cu.!" ca:;e, the flow direction) (lb) m = mass ( slu'.j} c. -accelera-tion in t:he flcw dir'":lic!:ion U : cross-sectional average velocity (ftjsec) = t:bn.H (sec) In fJ.,uid mechanics, ths c01:trol volurne mass is equal to: "t/here J 3 ""' of fluid (slug/f't.: ) .tsx = dist.anc.e incrl3ment ( ft j V== avsrage depth of vohrn-.e (ft:) 'l-... ,_,7 .. ; .;q1' 1 "'"'r .... ., h._.\..: .. :.._ o ..... e, t.ha ex-tarnal forces acting on '\:he control volume between sect i.ons l and. 2 in 2 3 are composed mainly of fou::: component.s: pressu::e, gravity, friction force, and the momentum change due to :=ainfall excess.

PAGE 27

(D,@= station number ft Ff = friction force F = gravity force g Fp = pressure i = rainfall intensity 0 = specific weight of fluid (\vater) ft = infiltration loss Figure 2.3 Control Volume to Derive Momemtum Equation 18 l"(y +-Ay) 2

PAGE 28

19 The component of pressure, Fp, for a unit-width overland flovr can be represented by: where Fp = -T Y 6. Y 1 Y 6 Y (2.7) y = specific weight of the fluid (lb/ft3 ) y = average depth of control volume (ft) 6 y = increment of water depth (ft) Equation (2.7) is legitimate when the channel bottom slope is mild. For a very steep slope, the hydrostatic force may be corrected by multiplying this term by The gravity force, Fg, is the component of the weight of water in the flow direction. The value is Fg = T sin e (y 6. x) For a small slope, sin e tan e s ; 0 therefore, Fg = TS0 (yAx) in s0 = channel bottom slope y = average depth of control 6.x = distance increment (ft) The friction force, Ff, is the ( 2. 8) ( 2. 9) ( 2 10) volume (ft) shear force the wetted perimeter. 'lrThich may be expressed as Ff = 1" sf (y !::. x) (2.11) on

PAGE 29

in which sf = friction slope y = average depth of control volume (ft) 6x = distance increment (ft) 20 It is assumed that the friction slope, sf, computed by the uniform-flow concept can be used for the nonuniform, and spatially varied flow. This assumption will be discussed in more detail below. The momentum change due to rainfall excess is equal to ( -u ie) 6X. Where U is the cross-sectional average velocity, ie, is the average rainfall intensity excess, and .6.X is the distance increment. Summing these four external forces, the change of momentum in the flow direction in Figure 3.2 can be expressed as : F = ma = Cf 6xy) a =6l? + Ff + Fg(Uie)6x (2.12) Dividing Equation (2.12) by (y x) yields F = 'dY + Ys ..,.Sf -Uie Y -0 I y ( 2 .13) The acceleration term, a, in Equation (2.5), can be further decomposed to local acceleration term, and the convective acceleration term,

PAGE 30

21 Substituting tha force term, F, in Equation (2.13) into Equation (2.5), the complete momentum equation for unsteady; two...:(Hmensional overland flo\-T is Ui oU + u -au -:-g + = g (So Sf) at dX ax y where g = gravitational constant (ftjsec2 ) S0 = channel bottom slope sf = friction slope. (2.14) When rainfall excess is negLigible, Equation-(2.14) is rsduced to the Saint-Venant that is, Equation (2.15). au + u cu + g + g cs0 -sf) = o ex ax ( 2 15) (1) (2) (3) (4) (5) Equation (2.15) is a nonlinear, hyperbolic partial differential equation. Due to the complexity in mathematics, this equation is usually solved numerically with specified boundary conditions. No general analyti-cal solution exists. When the vertical acceleration is negligible, the hydrostatic pressure and frictional resistance for an unsteady are the same as those in a steady flow. The terms in Equation (2.15) are, respectively:

PAGE 31

22 (1) the acceleration due to time variation. (2) the acceleration due to space variation. (3) the force due to water-surface variation. (4) the force due to the gravity, that is body force (5) the force due to the shear stress on wetted perimeter. Various approximations of Equation (2.15) were used, depending on the terms of the momentum equation which are considered (Yen, 1978). The kinematic-wave approximation is the most simple form where the friction slope is set to be equal to the channel bottom slope, and the pressure and inertial terms are negligible. The general format of the solution for a kinematic-wave flow is similar to the Manning's equation. q = ol ym ( 2 .16) in (1.49 1 n) (S0 )112 (Manning's equation), n =plane roughness or o{ = c (S0 ) 1 / 2 (Chezy' s equation) The value of m varied from 5/3 to 3, depending on the flow regime, turbulent or laminar. It has been found that making the value of m = 2 provides a good agreement with the observed data. The kinematic-wave approximation cannot model the flood-wave attenuation,

PAGE 32

23 and it is not adequate when the backwater effects are significant. The kinematic-wave model is accurate for a slowly rising flood wave on a steep slope such as overland flow (Ponce et al., 1978). The diffusive-wave approximation assumes that DU/Dt = o and is often called the equation of gradually varied flow. Gradually varied flow is a steady-flow situation where the depth gradually varies along the length of the channel. The improvement of diffusive-wave approximation is that it includes the consideration of a pressure term. In flood computations, peak attenuation and backwater effects can be simulated while two boundary conditions are specified. The quasi-steady dynamic-wave model neglects the local unsteady-velocity term, = o. This model was found to be valid to describe a long flood-wave propagation in a channel when lateral inflow is negligible. The complete dynamic wave model includes all terms shown in the momentum equation (equation 2.14). Backwater effects, downstream effects, and flow reversal can be simulated. To achieve numerical solutions for the complete momentum equation is time consuming and difficult in programming. Frequently, numerical instability causes divergency in computation. The

PAGE 33

24 complete dynamic-wave model has been used for both river routing (Chen, 1973; Simons and Chen, 1979; Smith, 1980; and Tucci, 1978) and sewer routing (Labadie et al., 1978; Yen, 1973; and Yevjevich and Barnes, 1970). In this study, the full dynamic-wave model is used for numerical analysis to completely model the dynamic behavior of overland flow. For mathematical simplicity, the momentum Equation (2.14) can further be transformed to a differential equation in terms of flow rate, q, and flow depth, y, only. Summing Equation (2.3) multiplied by U and Equation (2.14) multiplied by y yields: (2.17) An examination of Equations (2.3) and (2.17) shows that there are three unknowns: flow depth, y, flow rate, q, and friction slope, sf. In order to obtain numerical solutions for the two unknowns, flow depth, y, and flow rate, q, in the continuity equation (2.3) and the momentum equation (2.17), an empirical equation for estimating the friction slope, sf, is

PAGE 34

25 required. In this study, Manning's equation for a steady nonuniform flow is used to evaluate the friction slope sf. 2 2 sf = q n ( 2 .18) 10 2.22y 3 in which n = Manning coefficient q = flow rate (cfsjft) y = water surface depth (ft) Equation (2.18) is use by assuming that the head loss at a section for the nonuniform, unsteady flow is the same as that for a steady nonuniform flow having the same velocity and hydraulic radius of the section. This assumption has never been precisely confirmed by either experiment or theory, but errors arising from it are believed to be small compared with those ordinarily incurred in the use of the uniform-flow formula. over years of use, this assumption has been proven to be a reliable basis for engineering designs. 2.4 Horton's Infiltration Model For a given return period of rainfall, rainfall intensity, i (in/hour), varies inversely with the duration of the storm (Overton and Meadows, 1976):

PAGE 35

T a i = r b a+td in which Tr = return period (year) td = rainfall duration (minute) 26 (2.19) a, b = parameters varying with return period and geography. The penetration of water through the soil surface is kno1.-rn as infiltration. The infiltration rate or rate at which water enters the soil at the surface is controlled by surface conditions. Therefore, soil type is an important factor in determining the infiltration rate. When the soil has a large portion of well-graded fines, the infiltration is low. If the soil has several layers of horizons, the least-permeable layer \
PAGE 36

27 runoff from less intense storms in nonurbani z:ed areas, the data collected in the urbanized basins by the Urban Drainage and District, Colorado, indicate that the antecedent mositure condition has a limited effect on runoff (Drainage criteria Manual, 198.2). Other factors affecting the infiltration rates include slope of land, temperature, quality of water, age of lawn, and soil compaction (Drainage Criteria Hanual, 1982). The infiltration research has been conducted for many years, but there still is no generally recognized adequate quantitative model of natural infiltration. The reason is that infiltration is affected by the complex combination of soil characteristics, soil moisture conditions1 and other factors occurring in nature (Overton and Meadow, 1976; Sheaffer et al., 1976.) Hovrever, infiltration is extrgmely important, because it can influence not only the volurne and intensity of rainfall excess rate but also the timing of the runoff hydrograph. Some alternate represtations of infiltration are as follows: Horton's Model Holtan's Model Soil Conservation Service (SCS) Method Algorithms based on Darcys law andjor

PAGE 37

antecedent conditions Simple Abstraction Method These models range widely in complexity. Each of them permits initial infiltration rates to depend on some measure of antecedent soil moisutre and to vary with time during the storm. In the semi-arid region, such as the state of Colorado, the infiltration model proposed by Horton (1935) was found to provide a good balance between simplicity and reasonable physical description of the infiltration process by the Urban Drainage and Flood Control District in Denver (Drainage Criteria Manual, 1982). Horton's infiltration model states: (2.20) where ft = infiltartion rate at time t (injhr) f0 = initial infiltration rate (in/hr) fc = final infiltration rate (injhr) e = natural logarinthm base k = infiltration decay constant (1/sec) t = time (sec) 28

PAGE 38

29 Horton's equation is illustrated graphically in Figure 2.4. In developing this equation, Horton observed that infiltration is high at the beginning of the storm and decays exponentially to a steady state constant value, fc, as the pores in the soil become saturated. The coefficients and initial and final infiltration values are site specific and depend on the soil and vegetative cover complex. The infiltration rates are "capacity" rates, indicating that the soil-air interface must be saturated at all times. In practical terms, this means that it is assumed that rainfall rate is always greater than infiltration capacity rates, and hence some ponding will always result. This is a major disadvantage in the use of Horton's model, because rainfall rates are highly variable and therefore will_possibly fall below the values of capacity infiltration rates (Linsley et al. 1982; Overton and Meadows, 1976). The Urban Drainage and Flood Control District has analyzed a considerable amount of rain/runoff data. On the basis of this analysis, the values in Table 2.1 are recommended for use within the District in conjunction with the 1982 version of the Colorado Urban Hydrograph Procedure (CUHP). Most frequently found

PAGE 39

z H ... z 0 H E-i ...:I H l".:t! z H 30 INITIAL INFILTRATION f0 ... ...... .... ... .. .... -------FINAL INFILTRATION '. ... f c .. INFILTRATED BETWEEN tl & t 2 TIME --t IN SECONDS Figure 2.4 Schematic Representation of Horton ''s Equation

PAGE 40

TABLE 2.1 Recommended Horton's Equation Parameters SCS Hydrologic Soil Group* A B c D Decay Infiltration Cin/hr) Coefficient Initial f0 Final -fc k (1/sec) 5.0 4.5 3.0 3.0 1.0 0.6 0.5 0.5 0.0007 0.0018 0.0018 0.0018 Remark: *SCS = Soil Conservation Service 31

PAGE 41

32 within the District are scs Hydrologic Soil Groups C and D; however, areas of Group A and B soils may also be found within the District (Drainage Criteria Manual, 1982). Horton's Model is therefore selected in this study by using the available data source given by the Urban Drainage and Flood Control District as shown in Table 2.1. Assuming a spatially and temporally uniform rainfall intensity, i, over the entire watershed, and assuming that the slope, surface roughness, and soil texture are all uniform without any depression and detention losses over the entire studying area, the effective rainfall intensity, ie, for Horton's infiltration model becomes (2.21) in which i = uniform rainfall intensity (injhr) f0= initial infiltration rate (in/hr) fc= final infiltration rate (injhr) It can thus be seen that the effective rainfall intensity will approach a constant falue of (i-fc), when time reaches infinity.

PAGE 42

33 2.5 Nondimensionalization of Governing Equations Through a normalization process, the governing equations may be nondimensionalized to a more desirable form. The nondimensionalized equations the advantages of reducing variables, simplifying the computational algorithms, and comparing results of different investigators. It is a better way to represent the values in the numerical analysis (Woolhiser and Liggett, 1967). In formulating an overland-flow problem, the dimensional units involved are time, length, and force. Because accumulation of flow at a basin outlet eventually reaches its equilibrium state, the equilibrium quantities were used to normalize governing equations to produce more signficant results in terms of ratios to the equilibrium state. In this study, the equilibrium outlet depth, YE, is employed as the characteristic depth scale; the watershed length, L, is used as the characteristic length scale. The equilibrium time of concentration, TE, determined by the equilibrium effective rainfall (that is, i-fc), is used as the characteristic time scale. According to Guo (1984A, l984B), when m = 2 (Equation 2.16), YE and TE can be computed by

PAGE 43

(2.22) (2.23) where L = length of watercourse (ft) = watercourse characteristic factor The watercourse characteristic is equal to (1.49 (S0 ) 112}'n, which is constant for a given -"t."homogeneous watershed. Utilizing the above two normalization quantities, YE and TE, the following dimensionless parameters can be defined: X* = X Y* = y L YE t* = t q* = q YE (i-f ) L c ; ie = -e (i-f ) c n2 (i-f )L (q*)2 sf = c 2.22 YES/3 (Y*) 10/3. (2.24) When the above dimensionless variables in Equation (2.24) are incorporated into the continuity equation (2.3) and the momentum equation (2.17), the following equations result: 34

PAGE 44

(2.25) (2.26) where Fi, the rainfall Froude number, is equal to (i-fc)/(g YE)112 which is constant for a given rainfall intensity and watershed. 35 Equation (2.26) can further be simplified to the following form: aq* --+ = R at* where P = [ q*+ D(y*)2] y* (2.27) The parameter, D, reflects the effects of the watershed length, L, the equilibrium depth, YE, and the constant

PAGE 45

Fi. The parameter, Ki, reflects the effects of the slope, S0 the watershed length, L, equilibrium depth, YE, and the constant Fi. 2.6 Initial and Boundary Conditions 36 In the treatment of partial differential equations, the determination of a s9lution is dependent upon some given functional relationships which are broadly termed as initial/boundary conditions. Solu tions must satisfy the governing equations on the interior and also meet the requirements on the boundaries of the computational domain simultaneously. A necessary initial condition for the solution of overland flow is that all depths and discharges per unit width along the plane, 0 x L, must known at the initial given time, t0 A dry bed initial condition for which all depths and discharges are zero at the time rainfall commences (t = 0) has been studied and recommended by several investigators (Morgali and Linsley, 1965; Woolhiser and Liggett, 1967). A dry bed initial condition was adopted in this study. The necessary boundary conditions are overlandflow parameters at locations of x = o, and x = L as shown in Figure 2.1. Because there is no flow entering the flow plane through the upstream boundry (x = O), it is assumed that the flow through the upstream boundary

PAGE 46

37 of the plane is zero. This means that all flow will enter the flow plane as rainfall. All flow entering the upstream station (1-2) is accounted for by monitoring the outflow for this section. The depth at the uppermost boundary was allowed to rise as storage increase in the section by the storage routing (Morgali and Linsley, 1965). It is assumed that there is a straight-line water surface in this uppermost control volume. The storage between stations 1 and 2, therefore, is given by Storage 1 -2 (2.28) In this study, the storage volume method (Morgali and Linsley, 1965; Morgali, 1970) was used to calculate the depth at the most upstream boundary. Regarding the downstream boundary condition, the study performed by Morgali and Linsley (1965) concluded that the consideration of critical depth at the down-stream boundary condition was not necessary. Noncritical conditions gave good results for the shallow-water flow on a subcritical slope. The backward finite diffe,rence method (discussed in Section 3. 3) gave good results, and it was therefore used in this study.

PAGE 47

38 CHAPTER III NUMERICAL STUDY 3.1 Introduction Many of the differential equations for engineering problems cannot be easily solved by analytical methods. Consequently, many of these equations have been of little value in application in the past. With the advent of the high-speed digital computers, however, a number of very efficient numerical methods have been developed in the application of solution-solving to sets of partial differential equations. These numerical solutions are always associated with special boundary conditions which describe the physical constraints for the whole system. Therefore, a good numerical modeling is not only a method of computation but also a proper mathematical description of some physical constraints.

PAGE 48

39 According to Roache (1976), computational fluid dynamics is certainly not a pure theoretical analysis. If theinitial and the boundary constraints of the whole system are properly defined, computational fluid dynamics is closer to natural system representations than are physical experimentations or theoretical analyses. There are several basic approaches to transform a partial differential equation into finite-difference form, such as Taylor series expansion, polynomial fitting, integral methods, and control volume method. In general, it is assumed that a linear relationship between two adjacent points is applicable to estimate the time rate of change or the spatial rate of change in the differential equation. The detailed derivation and discussion of finite difference formulation for partial differential equations can be found in many references such as Rosenberg (1969), Smith (1971), and Roache (1976). The mathematical formulation of most problems in science and engineering involving rates of change with respect to two or more independent variables, usually representing time, length, or an angle, leads either to a partial differential equation or to a set of such

PAGE 49

40 equations. Special cases of the two-dimensional second-order equation. 2 2 aa+bil_+ ;}x2 + fJI + g = 0 ( 3 .1) (where a, b, c, d, e, f, and g may be functions of the independent variables x and y, and of the dependent variable f> occur more frequently than any other because they are often the mathematical form of one of the conservation principles of physics. The basic two-dimensional second-order equation can be categorized as two-dimensional elliptic equations, parabolic and hyperbolic equations. The best known among the two-dimensional elliptic equations are the Poisson's equation (3.2) and Laplace's equation (3.3). i 2fti + iP_ + G = 0 2 2 ax 'JY in which j is a physical variable and G is the source or sink term. (3.2) (3.3)

PAGE 50

41 These equations are generally associated with equilibrium or steady-state problems (Bascoi, 1970; Kersten, 1969). For example, the velocity potential for the steady flow of incompressible nonviscous fluid satisfies Laplace's equation. It is the mathematical way of expressing the idea that the rate at which fluid enters any given region is equal to the rate at which it leaves it. Problems involving time, t, as one independent variable usually lead to parabolic or hyperbolic equations. The simplest parabolic equation is derived from the theory of heat conduction. It depicts the temperature, u, at a distance, x units of length, from one end of a thermally insulated bar after t seconds of heat conduction. In such a problem, the temperatures at the ends of a bar of length, 1, are often known for all times. In other words, the boundary conditions are known. The simplest hyperbolic equation is the onedimensional wave equation. An example for this type of equation is used to describe the transverse displacement, u, at a distance, x, from one end of a vibrating string of length 1 after a time, t.

PAGE 51

42 3.2 Finite-Difference Method Finite-difference methods are approximate in the sense that derivatives at a point are estimated by difference quotients over a small interval. The solutions are not approximated in the sense of being crude estimates, however, the method generally gives solutions that are either as accurate as the data warrant or as accurate as. necessary for the purposes in which the solutions are required. From Taylor's theory (Kersten, 1969; Basco, 1970), for a function, Jf, and its derivatives are that is, finite and continuous function of x exists, then ;f (x + k) = f (x) + k ;/ (x) + 1/2 k 2 ;i (x) + 1/6 k 3 P"' (x) + I (x k) = / (x) -k ;t' (x) + 1/2 k 2 I (X) -1/6 k 3 Jl'' (x) + (3.4) ( 3. 5) in which k = infinitesimal fluctuation of variable x; ;( = first derivative of function p, etc.

PAGE 52

43 By neglecting all terms of higher power of fluctuation of k, equations (3.4) and (3.5) give ;/'' (x) = d 2 ,eifjdx2 = (l/k2 ) [p'(x + k) _2/ (x) + ;I (x -k)] ( 3. 6) f (x) = dp'jdx = (l/2k) [p (x + k) j (x k) ] (3.7) According to Figure 3.1, Equation (3.7) approximtes the slope of the tangent at point P by the slope of the chord, LR. This is the central-difference approxi-mation. If we approximate the slope of the tangent at point P by the slope of the chord, PR, or by the slope of the chord, LP, the methods are called the "forwarddifference method" and the "backward-difference method", respectively. The formulae are given by )I (x) = [fl (x + k) j (x) ] (3. 8) p (x) = (1/k) [p (x) f (x k) ] (3.9)

PAGE 53

i+l 0 i-1 Figure 3.1 AX/2 /2 i i+l X Illustration of Finite-Difference 44

PAGE 54

45 3.3 Central-Difference Method From Equation (3.8), it is clear that the forward-difference methods are defined by points to the right of the point P. On the other hand, the backwarddifference methods are defined by the points to the left of the point P. Equation (3.8) can alternatively be expressed by Equation (3.10): .6-ft(x) /boX =;I (x +box) ;I (x) (3 .10) Equation (3.9) can alternatively be expressed by equation (3.11): bop(x) /box = Ji (x) ;J (x box) (3 .11) Similarly, the central-difference methods are, therefore, defined by points located symmetrically with respect to point P. Equation (3.7) can be expressed by Equation (3.12): boft(X)/boX = jd (X +boX/2) -I (X 6X/2) ( 3 .12)

PAGE 55

46 Using the finite-difference approximations in solving partial differential equations, some error can arise in the computation. The discretization error is one of the errors which has been studied (Forsythe and Wasow, 1960). According to the study, the centraldifference method gives an averaged estimation from Equation (3.12); that is, it will give the approximation closest to the true value. The backward and the forward differences will give relatively poorer estimates. 3.4 Implicit and Explicit Methods The complete st. Venant equations cannot be solved analytically, and hence numerical schemes may be used to obtain discrete solutions. As an unsteady flow in nature, the flow domain in terms of computational grids is composed of time and space steps, as shown in Figure 3.2. The governing equations are approximated by the algebraic finite-difference equations, which are solved numerically. These equations can be solved by one of the following two methods. (1) Implicit schemes: the unknown parameters are 'expressed implicitly in simultaneous equations, which are solved by an appropriate solution technique. (2) Explicit scheme: the unknown parameters are expressed explicitly as functions of known parameters and solved directly.

PAGE 56

47 Miller (1971) has reported a complete treatise on the subject of the explicit and implicit schemes. The work of Amein and Fang (1969) has also provided notable computational examples. A complete description of the numerical solutions of the governing st. Venant equations of unsteady flows is given by Isaacson et al. (1956). An implicit and explicit grid network is shown in Figure 3.2 to contrast the two approaches. A formula which expresses one unknown pivotal value (point P, Figure 3.2) directly in terms of known pivotal values (points L, M, and R, Figure 3.2) is an explicit scheme. The method applied to the governing equations usually result in linear algebraic equations from which the unknown can be evaluated directly without iterative computations. A method in which the calculation of an unknown pivotal value (point P, Figure 3.2) necessitates the solution of a set of simultaneous equations (points L', M, and R', Figure 3.2) is an implicit method. The method involves nonlinear algebraic finite difference equations whereby the solution is obtained by iterative computation. There are advantages and disadvantages to both the explicit and the implicit methods. Explicit methods

PAGE 57

t'me,t j Figure 3.2 L' R' i-1 i+l distance, x Network of Grid Points in pomain 48

PAGE 58

49 are only conditionally stable, meaning that errors will grow as the solution progresses and are a function of the step sizes of time and distance. Trial-and-error solutions using computers are necessary to establish stability criteria. According to Ponce (1978), the explicit methods often result in a smaller computational time step (At) in order to obtain the desired accuracy. This limit on the length of usually makes explicit methods less efficient than implicit methods. The explicit methods are easy to program, however, and are generally easy to handle. Implicit methods were claimed to be unconditionally stable and more computationally efficient than are explicit methods. When applied to flood routing in open channels, the implicit methods have been found to be more stable, faster, and more accurate than other finite-difference methods (Amien and Fang, 1969; 1974). The implicit methods are considerably more difficult to program, however. Brakensik et al. (1966) and Amien and Fang (1969) concluded that the central-difference method provides stable solutions for an unsteady overland flow. The implicit central-finite difference with noncritical downstream boundary condition is, therefore, selected to perform numerical computations in this

PAGE 59

study. More details can be found in the study of Morgali and Linsley (1965). 3.5 Numerical Formulation of the Flow Equations The continuity and momentum equations are discretized in the continuous x -t plane, the dependent variables are depth and discharge, and the independent variables are time and distance. The numerical scheme of the four-point implicit scheme (Preissman, 1960) is adopted. 50 In this study, the finite-difference equations for computer programming are derived from Equations (2.24), (2.25), and (2.26). Using the general format of (3.13) (in which f 1 and / 2 and f 3 are flow variables) and the central-difference method, the general finite-difference form can be shown as -follows:

PAGE 60

'1 (i,j) -11 ,j2(i+l,j) = ( i 1 j ) ( 3 14) in which / 1 (i,j) =value of ;(1 at grid point (i,j), etc. = time increment (sec) 6X = distance increment (ft) i = station number along x-axis, IE j station number along = y-axJ.s, j JE IE = last station on x-axis, which the spatial axis JE = last station on y-axis, which the time axis Rearranging Equation (3.14), a more general form is obtained: 51 is is

PAGE 61

l1 = l1 (j,j-1) + (.O.t/26x) [j2 ( i-1, j) /2 ( i + 1, j) ] + (.6.t/ 2) [fo3 ( i 1 j ) + /3 (i,j-1)] ( 3 .15) Substituting the dynamic-wave routing parameters into equation (3.15), Equation (2.3) for the depth and Equation (2.17) for the flow rate per unit width at point (i,j) can then be solved by the following equations: 52 y(i,j) = y(i,j-1) + [q(i-l,j) -q(i+l,j)] [ie(i,j) + ie(i,j-1) ( 3. 16) q(i,j) = q(i,j-1) +(.c:.t/2Ax) [ (q(i -1,j)jy(i-1,j) + gy(i-1,j)2/2) -(q(i + 1, j)/ y (i + 1,j) + gy (i + 1, j)2/2)] +(.C:.t/2)g [y(i,j) + y (i,j-1)] (So-Sf) (3.17)

PAGE 62

53 in which y(i,j) = the value of depth at ith point in space and jth point in time, etc. sf = Average friction slope for time j and time (j-1) the nondimensionalized formula, as described above in Section 2.5, Equations (3.16) and (3.17) can further be transformed into: and .6tk Y*(ifj) = Y*(i,j-1) + 26X* [ q*(i-l,j) q* (i+l,j)J + [ie*
PAGE 63

54 Y (i, 1) = 0 where 1 i IE, IE is the outlet station and q ( i, 1) = 0 where 1 i IE Because there is no flow entering the upstream extremity, the discharges at the most upstream boundary of the plane is always zero (q = 0). Using the storage volume method (Morgali and Linsley, 1965), the depth at the upstream boundary, Equation (3.12), can be transformed to: Y* ( 1, j) = 2 [ ie ( 1, j) -q* ( 2, j) ] -y* ( 2, j) ( 3. 21) As to the downstream boundary conditions, Morgali and Linsley (1965) have tried both critical conditions, that is, free fall outlet and noncritical conditions. From their study, noncritical conditions gave good results for the shallow-water flow on a subcritical slope. Using Equations (2.25) and (2.27), the downstream boundary conditions are thus formulated by a backward finite-difference method as follows:

PAGE 64

and Y*(IE,j) q*(JE,j) = Y* (IE,j-1) [q*(IE-l,j) AX* q*(IE,j)] + + ie* (IE,j-1)] = q-k(IE,j-1) Axx P(IE,j)] +.P.t* [R (IE,j) + 2 R(IE,j-1)] (3.22) (3.23) The stability and accuracy, plus the con-55 sistency and convergence of the results of the numerical analysis have been studied (Basco, 1974; Meadows, 1984). In summary, the key points to reach the best approximation of a differential equation are to use (1) the "best" finite difference scheme, (2) the optimal At/AX ratio to give optimal accuracy, and (3) the numerical scheme of stable conditions. In this study, the central finite-difference meth6d is used because it gives the best approximate results among the forward-, backward-, and central-difference methods (Kersten, 1969; Rosenberg, 1969). The optimal At/AX ratio was determined by the sentivity test in this study for different cases.

PAGE 65

In computational procedure, the previous water profile is used as the initial guess for the water-surface profile at current time step. The new water-surface profile can be obtained by successive 56 iteration. During updating, the flow variables (depth and discharge) for the profile at jth time step and the flow variables at (j-l)th time step are kept unchanged. In this study, criteria similar to those used by Meadows et al. (1984) are used for judging the convergence. (n+l) (. 1 ) (rt) (. .+ 1 ) y 1,] y 1,] y*tnL(.i_, j+1 ) (3.24) and l *(n+1) ( .. l) *(n) ( .. + 1 ) q 1,] q 1,] -* (n) (. 1) q 1, J_ (3.25) in which y*(n+l} (i,j+l} = the value of Y*(i,j+l} after (n+l)th iteration, etc.

PAGE 66

57 3.6 Numerical computation Procedures The method of successive iteration was used in the computation. After each sweep of the computational domain, the flow variables (depths and discharges) are updated by new values until Equations (3.24) and (3.25) are satisfied for the entire domain. Figure 3. 3 the detailed logic computational flow chart which consists of the following steps: 1. Calculate the values of normalizing variables, YE and TE. 2. Initialize the values of flow variables, such as depth and discharge. 3. Advance time increment by t and compute the corresponding rainfall excess. 4. Adopt the solutions of previous time step as the initial guess. 5. sweep the computations of depth and discharge for the entire plane, that is, i = 1 to i = IE. 6. Update the values for depth and discharge after each sweep. 7. Chec]t convergence criteria. If Equations ( 3. 2 4) and (3.25) are satisfied, solutions are achieved for the current instant, then computation goes back to step 3. it goes back to step s.

PAGE 67

No I ie = !start I Check >3TE ] I i-f ie = t solut1ons of prev1ous 1me step as initial ess Sweep the dis char within 10-3 depth and sweep 58 Yef -f t No check: discharge = 0 No and depth = o Figure 3.3 Flow Chart of Computational Procedures

PAGE 68

59 In this study, the rainfall duration, td, is assumed to be as long as three times the equilibrium travel TE, of the basin. It is believed that solutions developed under this rainfall condition should provide enough information for general cases. Recession of _hydrograph starts when rain ceases. In the computer program, the rainfall intensity was set to be. zero after rain ceases, that is, i = 0 for t > td. The computation stops when the eniter catchment goes back to dry bed condition again, that is,. depth and discharge vanishes

PAGE 69

60 CHAPTER IV CASE STUDY As mentioned in Section 2.3, the kinematic wave approximation is a simplifiecation of the complete st. Venant equations which have been recognized to compose the most adequate model for describing open channel-flow conditions. As the kinematic wave formulation neglects all terms in the momentum equation except the friction slope and the channel-bottom slope, the formulation is less complicated than the complete dynamic equations, and therefore it requires less computational effort. Consequently, it is desirable to use the kinematic model whenever possible. The main purpose of this study was to use the numerical analysis to determine under what conditions kinematic wave solutions would give results close enough to those obtained using the complete dynamic equations. The analysis of dimensionless governing equations {Section 3.5) indicates that two parameters, the rainfall-kinematic-wave number, Ki, and the dynamic-wave number, D, should be the two most

PAGE 70

61 important parameters in the numerical study of overland flow. These two parameters defined as: in 'llhich so YE F l. J..r and D I< l. = = = = = s y = o E F. 2 L l. slope of plane equilibrium water surface (i-fc)/(g YE)l/2, g gravitational length of y 2 E plane constant (ft) 4.1 Description of Study Aspect ( 4 .1) depth (ft) is the (ftjsec2 ) (4.2) A hypothetical system was devised to evaluate the performance of kinematic-wave and dynamic-wave models in the watershed system. The purpose was to determine the criteria that would require that the complete st. Venant equations be used. This hypothetical system can be exemplified by a shallow sheet flow of excess rain ona homogeneous catchment such as a parking lot. The watershed area used in this study is described as follows:

PAGE 71

62 {1) The plane has unit width. (2) The length of the plane, L, was set to be equal to 500 ft. {3) The surface roughness, n, is uniform throughout the catchment. The Manning's value was set to be equal to 0.03. (4) The surface slope, s0 is uniform throughout the catchment. The value was ranged from 0.0005 to 0.1. (5) The soil at the plane surface is homogeneous in nature, and infiltration parameters are: f0 = initial infiltration rate = 4.5 in/hr fc = final infiltration rate = 0.6 injhr k = infiltration decay constant = 0.0018/sec The values are taken from the recommendations of the Urban Drainage and Flood Control District {1982) for soil type B as shown in Table 2.1. (6). The rainfall intensity, i, is uniform in space and also in time. The rainfall intensity value was set to be equal to 12 injhr with rainfall duration, td, of 3TE. The time, TE, is the equilibrium travel time of the basin as defined in Section 2.5.

PAGE 72

63 4.2 comparisons between Predictions of Dynamic-Wave and Kinematic-Wave Models Figures 4.1, 4.2, 4.3, and 4.4 show the effects of the rainfall-kinematic-wave number, Ki, on the peak ratio of flow depth, Y/YE (as defined in Section 2.5). The kinematic-wave solutions of each case are plotted as the solid line, and the dynamic-wave solutions obtained from this study are plotted as the dotted line. Several points of interest can be observed from these figures. First, Figure 4.1 (withQ(.= 1.11, and TE = 1301 sec) shows that the dynamic-wave solution is approximately 70% of the kinematic-wave solution when the rainfall-kinematic-wave number, Ki, is equal to 55 at the rainfall duration of 3TE. It also shows that when Ki = 55, there exists a "point of separation" where the kinematic-wave solution departs from the solution. In Figure.1, the point of separation is located at approximately 0.5TE (point A). After the point of separation, the dynamic effects of the flow become important. The reason for this behavior is that at the beginning of the rainfall, the flow is shallow, and thus there is no significant backwater effects on the plane (that is, the flow shows kinematic behavior). Only when the rainfall duration is

PAGE 73

Y/Y o.s._ 0 0 Figure 4.1 . Kinematic-wave solution Dynamic-wave solution ---6--------k-__ ..,.__ . 1 2 3 Illustration of Outlet Hydrograph A. L = 5000 ft, =. = 1301.3 sec, Ki = 55, i = 12 in/hr, f = in/hr, f = 0.6 in/hr, k = 0.0018 0 c .. t/TE 0'1

PAGE 74

long enough (td > O.STE, in this case), will the backwater effects be pronounced. The third point of interest is the different values of "time to peak" (defined as time of r concentration in this study) between the kinematic and 65 the dynamic solutions. The time to peak value for kinematic routing is approximately 1.13TE (point B, Figure 4.1), and that for dynamic routing is approximately 1.09TE (point c, Figure 4.2). By using the Horton's infiltration model, the effective rainfall intensity, ie, levels off to a constant value (i-fc) when time reaches infinity. It is noted that the dynamic-wave model has shorter time to peak than that of the kinematic-wave model. In-the practice of runoff prediction, time to peak is important. Generally the longer the time to peak is, the lower the peak flow rate, QP' will be. Therefore, kinematic-wave approximation would possibly underestimate Qp when backwater effects are significant. The fluctuation of the dynamic solutions as shown in Figure 4.1 may be attributed to (a) the limited accuracy of numerical computation due to high-speed digital computer and (b) the dynamic effects of flood routing. The peak discharge might fluctuate over the routing period. When the downstream backwater effect

PAGE 75

66 increases, the upstream water surface depth increases. As a result, the velocity and the inertial force of the flow increase, which causes the flow to change from dynamic to kinematic behavior. These natural resonant effects, changing from dynamic to kinematic and back to dynamic behaviors could cause the numerical solutions to fluctuate. Figure 4.2 (";litho<.= 1.22, and TE = 1243 sec.) shows that the dynamic-wave solution is about 80% of the kinematic-wave solution when the rainfall-kinematic wave number, Ki, is equal to 60 a.t a rainfall duration of JTE. In comparison with Figure 4.1, Figure 4.2 shows that when the number, Ki, increases (physically it means that the slope of the plane steeper, the equilibrium water depth becomes greater and the constant, F'i, becomes smaller), the point of separation starts at 0.45 TE (point A), which is earlier than that in Figure 4.1. This indicates that the duration of kinematic behavior is shorter because the higher Ki value causes an earlier backwater effect in the flow regime. In addition, the difference of the time to peak between the dynamic (point c, Figure 4.2) and the kinematic (point B, Figure 4. 2) models is reduced as compared with Figure 4.1. The reason is that when Ki value

PAGE 76

Kinematic-wave solution 1.0 Y/YE Dynamic-wave solution c ...... ----It ,1M 0.5 0 0 1 2 Figure 4.2 Illustration of OUtlet Hydrograph B. L = 500 ft, d= 1.22, T = 1243 sec, K. = 60, i = 12 in/hr, E 1 f = 3.5 in/hr, f = 0.6 in/hr, k = 0.0018 0 c .. 3 t/T E 0"1 ...,J

PAGE 77

68 increases, the whole flow regime behaves more similarly to that of the kinematic-flow model rather than that of the dynamic-flow model. Figure 4.3 (withcl= 1.31, and TE = 1198 sec.) shows that the dynamic-wave solution is approximately 88% of the kinematic-wave solution when the rainfall-kinematic-wave number, Ki, is equal to 65. Figure 4. 4 (with ol. = 1. 4 and TE = 1, 158 sec) shows that the dynamic-wave solutionis approximately 90% of the kinematic-wave solution when the rainfallkinematic-wave number, Ki, is equal to 70. Figures 4.3 and 4.4 also show the same effects of kinematic versus dynamic behavior as those shown in Figure 4.2 when comparing the point of separeation and the values of time to peak of the two solutions as discussed above. Figure 4.5 integrates the data shown in Figures 4.1 through Figure 4.4. The effects of various values of rainfall-kinematic-wave number, Ki, and the peak ratio of water surface depth, Y/YE, is demonstrated. As the term, Y/YE, approaches 1.0, the kinematic-wave. model is in good agreement with the complete dynamicflow model. Based on Figure 4.5, the following conclusions can be drawn:

PAGE 78

1.0 Kinematic-wave solution Y/Y E _k_ Dynamic-wave solution --6----_.a. -----'I 0.5 0 0 1 2 Figure 4.3 Illustration of Outlet Hydrograph c. L = 500 ft, 1.31, T = 1198 sec, K. = 65, i = 12 in/hr E 1 f = 3.5 in/hr, f = 0.6 in/hr, k = 0.0018 0 c 3 .. t/T E \0

PAGE 79

1.0 Y/YE 0.5 0 0 Figure 4.4 c ..... Dynamic-wave solution ... --k----____., -----_...,__ "' 'I / ., 1 2 Illustration of Outlet Hydrograph D. L = 500 ft, 1.40, T = 1158 sec, K. = 70, i = 12 in/hr, E 1 f0 = 4.5 in/hr, fc = 0.6 in/hr, k = 0.0018 < 3 t/TE -.J 0

PAGE 80

1.0 Kinematic-wave solution Y/YE .. tJ
PAGE 81

72 (1) When the rainfall-kinemaitc-wave number increases, the ratio of Y/YE becomes higher and approaches 1.0 when the rainfall-kinematic wavec number is greater than 180. This indicates that the flow behaves in a manner more similar to that of the kinematic model when the value of Ki increases. (2) When the number increases, the point of separation where the dynamic-flow behavior becomes significant starts at an earlier time. This shows that the duration of kinematic behavior is becoming less at the beginning of the rainfall when Ki increases. (3) When the number increases, the values of time to peak of the kinematic and the dynamic routing are closer to each other. This indicates that the flow regime behaves more similar to that of the kinematic model when Ki increases. The relationship of the rainfall-kinematic-wave number, Ki, to the peak ratio are summarized and shown in Figure 4. 6. The figure sho"t.tTS that when the rainfall-kinematic-wave number is greater than 180, .the kinematic wave reaches almost 100% agreement with the dynamic-wave model. It is noted that when the value of the rainfallkinematic parameter, Ki, is smaller than 50, the

PAGE 82

73 kinematic-wave solution is a poor approximation (Figure 4.1), and the complete dynamic equation should be used. The maximum error, however, in the outflow hydrograph with Ki = 50 is in the order of 30%, and it decreases rapidly as Ki increases. Figure 4.6 also shows that when Ki increases to 150, the dynamic-wave model merges closely with the kinematic model and will give results within 0.5% of the dynamic-wave model (point A). Another paramter, the dynamic-wave number, D, which is defined by yE2; 2L2 Fi2 is also investigated. The parameter is equivalent to Ki x (1/50 ) x YE/2L, which is equal to the (rainfallkinematic-wave number) x (1/50 ) x YE/2L. The study showed that when the dynamic-wave number is 38, the dynamic-wave solution is approximately 70% of the kinematic-wave solution. When the dynamic-wave numbers are 33, 27, and less than 6, the dynamic-wave solutions are 80%, 90%, and 100%, respectively. The results are depicted in Figure 4.7. It shows that when the dynamic-wave number is less than the kinematic-wave approximation will give results within 5% of the dynamic-wave solutions. The parameter, D, can also be used as a complementary criterion to check when the dynamic-wave routing can be approximated

PAGE 83

A Ymt 1.0j(----------Ye -. -_ .-.-.9 .8 .7 I I I .6 200 so 100 150 So Yr;: Figure 4.6 Effects of Rainfall-Kinematic Wave Number on the Outlet Flow Depth. L = 500 ft, i = 12 in/hr, f0 = 4.5 in/hr, = 0.6 in/hr, k = 0.0018 F2 L I --....1

PAGE 84

yml YE l.Oiat .9 .8 .7 .6 I 0 5 10 15 20 25 30 Figure 4. 7 Effects Dynamic Wa.ve Number on the Outlet Flow Depth 35 40 t YE .. 2L2 l. --.1 Ul

PAGE 85

, 76 by the kinematic-wave routing. The results show that as the dynamic-wave number, D, becomes smaller (that is, the slope of _the watershed increases), the wave routing behaves closer to the kinematic model. This result agrees with that of the rainfall-kinematic-wave number, Ki. 4.3 Evaluation of the applicability of the KinematicWave Model by using Available Laboratory Data To provide more comparisons between the predictions of the dynamic-wave and kinematic-wave modeling and verfication of the criteria established in this study, the laboratory data obtained by Bentner et al. (1940) were used. The experiments of Bentner et al. were performed on a watershed located in Tucson, Arizona. It is approximately 24 ft long and 6 ft wide. For the most part, the area is sparsely vegetated or has bare soils, and the soil is rather tight. Runoff hydrographs in Figures 4.8A, 4.8B, 4.9A, and 4.9B were collected by Bentner et al. (1940). Woolhiser and Liggett (1967) concluded that the kinematic-wave equation is a good approximation for most overland flow situations of hydrologic interest, because Kq (defined in section 1.1) rarely falls below 10. In general, Kq values smaller than 10 can be expected only for short, smooth planes with mild slopes and high rates

PAGE 86

77 of lateral inflow. These conditions are not often found in rural areas but are possible in urban areas. At equilibrium state, the normal depth (H0 ) used in calculating the kinematic-wave number, Kq, in Equation (1.2), is equal to the value of equilibrium water depth (YE) used in calculating the rainfallkinematic-wavenumber, Ki, in Equation (4.1). Table 4.1 the parameters used, s0 (slope of the plane), n (roughness of the plane), and tL(watershed characteristic constant) and the calculated values YE (equilibrium water-surface depth), qE (equilibrium flow rate) and Kq (kinematic-wave number), and Ki (rainfall-kinematic-wave number) of the four cases represented in Figures 4.8A, 4.8B, 4.9A, and 4.9B. It is noted that for all laboratory cases the values of Kq, which are the criteria used by Woolhiser and Liggett (1967) are greater than 10 (within 10% of the kinematic approximation). The values of Ki, which are the criteria established in this study, are greater than 70 (within 10% of the kinematic approximation). In addition, for a highly kinematic wave flow, the two values converge. All laboratory data show kinematicflow behavior. Therefore, this study also concluded that the kinematic-wave equation is a good approximation for most, but not all, overland flows.

PAGE 87

f =1 16 + ( 4 o 1 -1 16) eu eOt DRY RUN NO. 25 ::::s 0 .c. I I i:: .;.c
PAGE 88

-36.7ot r = 0.62 + .( 5. 36 -0.62)e .,.. __ ... .4,. ,t;Er 11UN NO. 26 I :J (: -,.10 _'t 0 :3. .Q) c. .A .5 -,-/.E 1 ca -a: !--' ... J.. I -. time in minutes \ 10 -20 30 Figure 4. 8B Hydrograph of Laboratory Data, Set 1, L = 24 ft, S0 = .089, n = 0.136 79

PAGE 89

I l I f = 0.50 +'(7.20 0.50)e3 L 6ot -i DRY RUN NO .t 8 g41 c I !31 i l ..L ...... & a .... .E .... Q) -ca 0: 2I 1-time in minutes 30 \ 2() Figure 4. 9 A Hydrograph of Laboratory Data, Set 2 L = 24 ft, 80 = 0.02, n = 0 . 039 80

PAGE 90

.. I I r (8 .14 0. 44) e -+'t = 1" H ...... "" 4 .. I '-WET RIJN NO. 49 ::J 0 i:..J .JJ .c '-Q) a. !.:. .5 3. --. fl .s r C1) ca a: .2. .... 1 30\ 1 time in _finutes 10 0 Figure 4. 9B Hydrograph of Laboratory Data, Set 2, L = 24 ft, s = 0.02, n = 0.045 o 81

PAGE 91

-..:. TABLE 4.1 Illustration of Laboratory watershed Analysis Based on Four Cases in Figure 4.8A, 4.8B, 4.9A, and 4.9B Parameters Slope, s0 Roughnes.s, n Watershed Characteristic Constant, cJEquilibrium Depth, 'LE (ft) Equilibrium flow rate* qE (cfs/ft) Equilibrium flow rate**, qE (cfsjft) Kinematic-Wave Number, Kq Rainfall-Kinematicwave Number, Ki Case 1 0.089 0.171 2.6 0.0202 0.00108 0.00106 25012 25050 Case 2 Case 3 case 4 0.089 0.02 0.02 0.136 0.039 0.045 3.7 5.4 4.7 0.0205 0.0170 0.0185 0.00138 0.00156 0.00161 0. 00138 0.00156 0.00161 15257 1832. 2042 15268 1840 2057 ' Remark Kinematic Kinematic Kinematic Kinematic *Laboratory Data **Numerical Solution

PAGE 92

83 CHAPTER V AND CONCLUSIONS 5.1 Conclusions The use of the complete st. Venant equation with Horton's infiltration formula for unsteady overland flow wave routing is a sophisticated and complicated wave-routing model. The use of such models has not been common among practicing engineers due to their complexity and the requirements for large amounts of computer time. 1. The rainfall-kinematic-wave number, Ki, which consists of the slope, s0 the equilibrium flow depth at outlet, YE, and the rainfall Froude number,. Fi, can serve as a criterion for when kinematic-wave assumptions are adequate and acceptable. Hhen the rainfall-kinematic-wave number is greater than 180, the kinematic-wave approximation will give results within 5% of the dynamic-wave approximation. 2. When the rainfall-Idnematic-wave number increases, the point of separation starts at an earlier

PAGE 93

time. This sho"tots that the duration of kinematic behavior becomes shorter at the beginning of the rainfall when.Ki is larger. Furthermore, the values 84 of time to peak of the kinematic and the dynamic routing are closer to each other, when the rainfall-kinematicwave number, Ki, increases. This indicates that the flow regime behaves closer to the kinematic model when Ki increases . 3. Another parameter, the dynamic-wave number, D, which consists of the equilibrium flow depth at outlet, YE, watershed length, L, and the rainfall Froude number, F i, serves as another checJdng parameter for determining when kinematic-'tvave assumptions are adequate and acceptable. It was concluded that when the dynamic-wave number less than 6, the approximation give results within 5% of the dynamic-wave approximation. 4. The testing of the kinemaitc-wave model by laboratory data shows good agreement .between this study and that of Woolhiser and Liggett (1967). It is concluded that the kinematic-wave equation is a good approximation for most, but not all, of the real world overland flows.

PAGE 94

85 5.2 Recommendations for Further Study This study incorporates the Horton infiltration model into the complete dynamic st. Venant equations; that model was not included in the study of Woolhiser and Liggett (1967). Therefore, the effect of rainfall, rather than the effect of runoff, .is modeled. The guideline parameter, Ki (rainfall-kinematic-wave number), can readily be calculated for a given watershed. The parameter used by Woolhiser and Liggett (1967) Kq (kinematic-waver number) cannot be calculated easily for a given watershed. This study offers a much better and more complete dynamic-wave modeling. Further study could add several features to the routing model. The comparison of kinematic-wave routing to full dynamic-wave routing should be expanded. For example, tests should be conducted to consider the plane of various lengths and various roughness, etc., to determine if changes in these parameters would alter the criteria developed. Another recommendation is for additional applications of the model. As different situations are encountered, obvious needs for improvements will become apparent. More study can be done to refine the numerical modeling to serve as a reliable engineering tool.

PAGE 95

REFERENCES Amein, M., and Fang, C. s. "Implicit Flood Routing in Natural Channels" J. Hyd. Div., AStE, HY 12, p. 2481-2500 (1970). Basco, D. R. "Introduction to Numerical Method" Verification of Mathematical and Physical Models in Hydraulic Engineering, ASCE, p. 280-302 (1978). Beutner, E. L., Gaebe, R. R., and Horton, R. E. "Sprinkled-plat Runoff and Infiltration-experiments on Arizona Desert Soils" Trans. Amer. Geophys. Union, p. 550-558 (1940). Behlke, c. E. "The Mechanics of overland Flow" Ph. D. Dissertation, Stanford University, Stanford, California (1957). Brakensik, D. L., et al. 86 "Numerical Technique for Small Watershed Flood Routing" u.s. Agric. Res. Serv. Publ., ARS 41:113 (1966). Brater, E. F., and King, H. w. "Handbook of Hydraulics" McGraw-Hill Book Company, New York (1960). Chen, c. w. "Computer Simulation of Urban Storm Water Runoff" J. Hyd. Div., ASCE Hy 2, p. 289-301 (1971). Chen, Y. H. "Mathematical Modeling of Water and Sediment Routing in Natural Channel," Ph. D. Diss., Department of Civil Engineering, Colorado State University, Fort Collins, Colorado, June (1979). Chow, v. T. "Open Channel Hydraulics" McGraw-Hill Book Company, New York (1959). Drainage Criteria Manual Published by Urban Drainage and Flood Control District, Denver, Colorado (1982).

PAGE 96

Eagleson, P. s. 11Dynamic Hydrology,11 P 344-346, McGraw-Hill Book Company, New York (1970). Forsthye, G. E., and Wasaw, w. R. 11Finite Difference Methods for Partial Differential Equations11 John Wiley and Sons, Inc., New York (1960). Fread, D.L. 11Technique for Implicit Dynamic Routing in Rivers with Tributaries11 Water Res. Research, 9 (4), p. 918-926 (1973). Gburck, w. J., and overton, D. E. 11Subcritical Kinematic Flow in a Stable Stream11 J. Hyd. Div. ASCE, HY9, p. 1433-1446 (1973). Grace, R. A., and Eagleson, P. s. 11The Modeling of overland Flow11 Water Res. Research, p: 393-403 (1966). Greebm, W. H., and Ampt, G. A. 11Studies on Soil Physics: I. Flow of Air and Water through Soils11 Journal of Agr. Sci., v. 4, p. 1-24 (1911). Guo, c. Y. "Effect of Infiltration on Runoff Hydrograph" 87 Water today and tomorrow, preoceedings, ASCE, Irregation and Drainage Conference in Flagstaff, Arizona, July, p. 485-492 (1984A). Guo, c. Y. "Fluid travel time" Water for Resources Development, Proceedings of ASCE Hydraulic Conference in Coeur d'Alene, Idaho, p. 380-384, August (1984B). Henderson, F.M. nopen Channel Flow11 Macmillan Publishing Co., Inc., New York (1966). Henderson, F. M. and Wooding, R. A. : "Overland Flow and Groundwater Flow from a Steady Rainfall of Finite Duration" Geophysical Research, v. 69, No. 8, p. 1531-1540 (1964).

PAGE 97

Holton, H. N., and Minshall, N. E. "Plot Samples of Watershed Hydrology" u.s. Agri. Res. Serv. Publ., ARS 41:33 (1968). Horton, R. E. "The interpretation and application of runoff plot experiments with reference to soil erosion problems" Soil Sci. Soc. Amer. v. 3, p. 340-349 (1938). Horton, R. E. "Approach toward a Physical Intterpretation of Infiltration Capacity" Proc. Soil Sci. Soc. Amer., v. 5, p. 399-441 (1939). Isaacson, E., et al. "Numerical Solution of Flood Prediction and River Regulation Problems: Report III" 88 Report No. IMN-NYU235, Inst. of Mate. and Sciences, New York University, New York, New York (1956). Izzard, c. F. "Hydraulics of Runoff from Developed Surfaces" Proceeding, 26th Annual Meeting of the Highway Research Board, v. 26, p. 129-146 (1946). Kersten,! R. D. "Engineering Differential Systems" McGraw-Hill Book Company, New York (1969). Kuelegan, G. H. "Spatially variable discharge over a sloping plane" Trans. Amer. Geophysical Union, Part VI, p. 956-958 (1944). Labadie, J. w., Morrow, D. M., and Lazaro, R. c. "Urban stormwater control package for automated real-time systems" Colorado State University Report c 6174, prepared for u.s. Depart. of Interior, Office of Water Research and Technology, Decembrer (1978). Lazaro, T. R. "Urban Hydrology" Ann Arbor Science Publishing Inc., Ann Arbor, Michigan (1979). Liggett, J. A. "Unsteady Open Channel Flow with Lateral Inflow" Ph. D. Dissertation, Stanford University, Stanford, California (1959).

PAGE 98

Liggett, J. A., and Woolhiser, D. A. "Difference Solution of the Shallow-Water Equation" Proc. ASCE J. Mech. Div. 9B (EM2), p. 39-71 (April, 1967) Lighthill, M. H., and Whithan, G. B. 89 "On Kinematic Waves, I, Flood Movement in Long rivers" Proc. Royal Soc. of London, Ser. A., v. 229, p. 281-316 (1955). Lindsley, R. A. "Hydrology for Engineering" McGraw-Hill Book Company, New York (1958). Meadows, M. E. "Finite element simulation of overland flow" Advances in irregation and drainage, surviving external pressure, AASCE, p. 466-475 (1983). Miller, w. A., Jr. "Numerical solution of the equation for unsteady open channel flows" School of Civil Engineering, Georgia Institute of Technology, Atlanta, GA (1971). Morgali, J. R. "Laminar and Turbulent Overland Flow Hydrographs" J. Hyd. Div. ASCE, HY2, p. 441-460 (1970). Morgali, J. R., and Linsley, R. "Computer Analysis of Overland Flow" J. Hyd. Div. ASCE, HY3, p. 81-100 (1965). Musgrave, G. W. "How Much of the Rain Enters the Soil" Yearbook of Agriculture, p. 151-159 (1955). Overton, D. E., and Meadow, M. E. "Stormwater Modeling" Academic Press, New York. (1976). Philip, J. R. "Theory of Infiltration" Advances in Hydroscience, Academic Press, New York, v. 5, p. 215-296 (1969). Ponce, V. M., Li, R. M., and Simons, D. B. "Applicability of kinematic and diffusion models" ASCE Journal of Hydraulics Division, v. 104, no. HY3, p. 353-363 (1978).

PAGE 99

Preissman A. "Propagation des intumescences dans les cannaux et rivieres" 90 1st Congress de l'Assoc. Francaise de Calcul, Grenoble, p. 443-442 (1960). Price, R. K. "Comparison of four numerical metnods for flood routing" Journal of the Hydraulic Division, ASCE, v. 100, no,. HY7, p. 879-899 (1974). Richey, R. P. "The Fundamental Hydraulics of Overland Flow" Ph D. Dissertation, Stanford University, Stanford, California (1954). Roache, P. J. "Computer Fluid Dynamics" Hermosa Publisher Inc., New Mexico (1976). Schreiber, D. L., et al. "Subcritical Kinematic Flow in a Stable Stream" Journal of Hydraulic Division, ASCE, HY3, p. 429-446 (1972). Scheaffer, J. R., et al. "Urban Storm Drainage Management" Marcel Dekker Inc., New York (1976). Simons, D. B., Chen, Y. H., and Saez-Benito, J. M. "A mathematical model of Pools 5 through 8 in the Upper Mississippi River System" Colorado State University Report CER 79-89, DBS-YHX-JMS 21, prepared for u.s. Army Corps of Engineers, St. Paul District, December (1979). Smith, G. D. "Numerical Solution of Partial Differential Equations" Brunel College of Advanced Technology, London (1980). Smith, R. E. "Documentation for program KINEROS" Unpublished report, Colorado State University, Fort Collins, Colorado (1980). Tucci, c. E. M. "Hydraulic and water quality model for a river network" Ph. D. Dissertation, Department of Civil Engineering, Colorado State University, Fort Collins, Colorado (1978).

PAGE 100

von Rosenberg, D. u. "Methods for the Numerical Solution of Partial Differential Equations" Elsevier Publishing Company, Inc., New York (1969). Woo, D. c., and Brater, E. F. "Spatially Varied Flow from Controlled Rainfall" Proc. ASCE, v. 88, p. 31-56 (1962) Wooding, R. A. "Hydraulic Model for the Catchment" J. Hydrology, 3, p. 254-267 (1965). Woolhiser, D. A., and Liggett, J. A. "Unsteady, One-Dimensional Flow Over a Plane The Rising Hydrograph" Water Res. Research, 3(3), p. 753-771 (1967). Yen, B.. c. "Methodologies for Flow Prediction in Urban Storm Drainage System" Research Report No. 72, University of Illinois, Champaign-Urbana, Water Res. Center, January (1973). Yen, B. c. "Storm Sewer System Design" University of Illinois, Champaign-Urbana, Illinois (1978). Yevjevich, v. and A. H. "Flood Routing through Storm Drains" Colorado State University, Hydrology Papers, p. 43-46 (1970). 91

PAGE 101

APPENDIX I NOTATIONS The following symbols are used in this thesis: a = acceleration F = force ft = infiltration rate f0 = initial infiltration rate fc = final constant rate i = rainfall intensity ie = rainfall excess K = exponent in Horton's infiltration model L = length of watercourse m = mass q = discharge S0 = bottom slope sf = friction slope td = rainfall duration = time spacing in the numerical computation t = time after rain starts TE = equilibrium travel time X = distance from upstream extremity = distance spacing y = depth variable Y = outlet depth YE = equilibrium flow depth 92

PAGE 102

APPENDIX I (continued) = watershed characteristic constant Subscripts: E = equilibrium travel time 1 = outlet parameter m = maximum depth 93