JDY1J.
AMICS AND KINEMATICS OF OVERLAND FLOWS
by
ShiuMing 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
ShiuMing Huang
has been approved for the
Department of
Civil Engineering
by
Ill
Huang, ShiuMing (M.S., Civil Engineering)
Dynamics and Kinematics of Overland Flows
Thesis directed by Assistant Professor James C. Y. Guo
The dynamicwave model for an overland flow was
assessed using the implicit central finitedifference
method. Flow equations were derived from the continuity
and momentum principles. This study attempts to
delineate the limitations of kinematicwave approxi
mation by performing a series of numerical studies using
the flow equations. The nondimensional dynamicwave
model was combined with Horton's infiltration formula
and numerically solved to establish a model representing
twodimensional overland flow over a pervious surface.
The comparisons between kinematicwave solutions and
dynamicwave solutions provide general guidelines and
further establish new criteria for determining when
kinematicwave 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
dynamicwave 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, ChiI, 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 DynamicWave Model and
Kinematic Wave Model ... 63
4.3 Evaluation of KinematicWave
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 FiniteDifference 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 RainfallKinematicWave Number . 71
4.6 Effects of RainfallKinematicWave
Number on the Outlet Flow Depth.............74
4.7 Effects of DynamicWave 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 waterresources 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 highspeed digital computers,
numerical schemes have been developed and used to obtain
discrete solutions for complicated openchannel flow
situations. Prior to the use of highspeed computers,
various simplified forms of the St. Venant equations
were used to approximate different hydrological
problems.
The kinematicwave approximation is the most
common simplification of the St. Venant equation. It
assumes that the friction slope, Sf,is balanced by the
channelbottom slope, SQ, but it ignores the pressure
and inertial terms. Lighhill and Whitham (1955) have
found that the floodwave propagation in overland flow
is mainly kinematic.' Wooding (1965) applied the
kinematicwave approximation to solve the partial
differential governing equations describing the flow for
a vshaped chatchment flow. His solutions represented
complete runoff hydrographs and the time of concentra
tion. The general solution for the kinematicwave
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 (11)
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 kinematicwave
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 kinematicwave theory is
the simpliest form of the complete dyanamicwave model.
With the simplification, kinematicwave 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 kinematicwave theory has its
limitation. Woolhiser and Liggett (1967) determined the
limitations of kinematicwave 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 kinematicwave
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 kinematicwave 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 kinematicwave approximation with Horton's infiltra
tion formula to analytically solve the propagation of a
twodimensional 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 kinematicwave approximation limitations under a
nonuniform rainfall excess due to infiltration decay
characteristics.
6
This study attempts to determine the limitations
of kinematicwave approximation by performing a
numerical study which combines the dynamicwave model
with Horton's infiltration formula to numerically solve
the propagation of twodimensional overland flow on a
pervious surface. The results from comparisons between
kinematicwave solutions and complete dyanamicwave
solutions can provide general guidelines, and can
further establish new criteria for determining when
kinematicwave 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 dynamicwave
equations representing overland flow resulting from a
uniform rainfall and a variable infiltration rate. The
discrete computational solutions were then compared with
the kinematicwave solutions developed by Guo (1984A and
1984B) in order to delineate the limitations of the
applicability of kinematicwave solutions.
The kinematicwave method is an approach which
has been proven to be a very useful tool in stormwater
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 kinematicwave 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 kinematicwave 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 finitedifference
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 overlandflow 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 overlandflow
hydrograph were numerically tested. They proposed a
semitheoretical 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 kinematicwave solution gives very accurate
results.
Morgali (1970) did a followup study on laminar
and turbulent overland hydrographs. In setting up the
solution of finitedifference 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 overlandflow 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 overlandflow plane. Giburek and
Overton (1973) studied subcritical kinematic flow in a
stable stream. An evaluation of the shallowwater
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 kinematicwave theory with
10
Horton's infiltration model. The accuracy of solution
is confirmed by numerical computations based on an
implicit finitediffernece 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 decayinfiltration model to the unsteady
dynamic wave routing and to further investigate the
limitations of the kinematicwave approximation.
11
CHAPTER II
ANALYTICAL MODEL
2.1 Introduction
The governing equations of a dynamicwave model
for an overland flow are derived from the continuity and
momentum principles. The equations were first applied
to the overlandflow problem by Keulegan (1944), and
later Morgali and Linsley (1965) did a study on computer
analysis of overland flow. The schematic illustration
of twodimensional flow is shown in Figure 2.1. The
accumulation of runoff on a unitwidth 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,
openchannel 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 = watersurface 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:
(qjAtqQAt) + (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
gTSr
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
twodimensional 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 twodimensional 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 axtaxnal 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 unitwidth
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 8SQ; (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 uniformflow 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 crosssectional
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/ax.
21
Substituting tha force term/ F, in Equation (2.13) into
Equation (2.5), the complete momentum equation for
unsteady, twodimensional 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 SaintVenant 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 watersurface 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 kinematicwave
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 kinematicwave 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 kinematicwave
approximation cannot model the floodwave attenuation,
23
and it is not adequate when the backwater effects are
significant. The kinematicwave model is accurate for a
slowly rising flood wave on a steep slope such as
overland flow (Ponce et al., 1978).
The diffusivewave approximation assumes that
DU/Dt = 0 and is often called the equation of gradually
varied flow. Gradually varied flow is a steadyflow
situation where the depth gradually varies along the
length of the channel. The improvement of
diffusivewave 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 quasisteady dynamicwave model neglects
the local unsteadyvelocity term, ^U/^t = 0. This model
was found to be valid to describe a long floodwave
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 dynamicwave 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 dynamicwave 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 uniformflow 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 wellgraded fines,
the infiltration is low. If the soil has several layers
of horizons, the leastpermeable 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 semiarid 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 soilair 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 (ifc),
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 overlandflow 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 (ifc)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* dx* 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)
(ifc)/(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* 9x*
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 (12) 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
straightline water surface in this uppermost control
volume. The storage between stations 1 and 2,
therefore, is given by
Storage
12 (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
shallowwater 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 highspeed digital
computers, however, a number of very efficient numerical
methods have been developed in the application of
solutionsolving 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 finitedifference
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 twodimensional
secondorder 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 twodimensional secondorder equation
can be categorized as twodimensional elliptic equa
tions, parabolic and hyperbolic equations. The best
known among the twodimensional 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 steadystate 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 FiniteDifference Method
Finitedifference 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
singlevalued; 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 centraldifference 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 "backwarddifference method",
respectively. The formulae are given by
/ (x) =(l/k)[/ (x + k) / (x)] (3.8)
(3.9)
Figure 3.1 Illustration of FiniteDifference
Method
}
45
3.3 CentralDifference Method
From Equation (3.8), it is clear that the
forwarddifference 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 centraldifference 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 finitedifference 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 finitedifference 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. Trialanderror
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 finitedifference 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 centraldifference method
provides stable solutions for an unsteady overland
flow. The implicit centralfinite 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 fourpoint implicit scheme (Preissman, 1960) is
adopted.
In this study, the finitedifference 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 centraldifference method, the general
finitedifference 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 xaxis,
1 ^ i^ IE
3 = station number along yaxis,
1 < j < JE
IE = last station on xaxis, which is
the spatial axis
JE = last station on yaxis, which is
the time axis
Rearranging Equation (3.14), a more general form is
obtained:
52
/l (i/j) = /i (jfj1)
+ (At/2Ax)[/2
 /2 (i + 1, j) ]
+ (At/2)[/3 (i, j )
+ /3 (i,j1)]
(3.15)
Substituting the dynamicwave 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,jl) + (At/2Ax)[q(il,j) q(i+l,j)]
+(At/2)[ie(i,j) + ie(i,jl)
(3.16)
q(i/j) = q(i,jl) +(At/2Ax)[(q(i 1,j)/y(il,j)
+ gy(il,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,jl)]
(SoSf)
(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 (j1)
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*(il,j) 
q* (i+l,j)] + [ie*(i,j) +
ie* (i,j1)] (3.18)
and
q*(i,j) = q*(i,jl) + ^L. [P(il,j) 
2Ax*
P (i+l,j)] + Â£!.[R(i,j) +R(i,j1)]
(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 shallowwater flow on a subcritical slope.
Using Equations (2.25) and (2.27), the downstream
boundary conditions are thus formulated by a backward
finitedifference method as follows:
55
y*(IE,j) = y* (IE, j1) + JÂ£L [q* (IE1, j) 
AX*
q*(IE,j)] + *l[ie*(IE,j) +
ie* (IE,j1)]
(3.22)
and
q* (JE, j ) = q* (IE, j1) + [P(IE1, j ) 
Ax*
P(IE,j)] +.^L [R (IE, j) +
R(IE/jl)]
(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 finitedifference
method is used because it gives the best approximate
results among the forward, backward, and
centraldifference 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
watersurface profile at current time step. The new
watersurface 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 (jl)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 Tjr*.
Â£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* >3Tp 1 Ye s
* I S & 4
Ae iff 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 channelflow
conditions. As the kinematic wave formulation neglects
all terms in the momentum equation except the friction
slope and the channelbottom 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 rainfallkinematicwave number, K^,
and the dynamicwave 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)
(ifc)/(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 kinematicwave and dynamicwave
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 DynamicWave and
KinematicWave Models
Figures 4.1, 4.2, 4.3, and 4.4 show the effects
of the rainfallkinematicwave number, K^, on the peak
ratio of flow depth, Y/YE (as defined in Section
2.5). The kinematicwave solutions of each case are
plotted as the solid line, and the dynamicwave
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 dynamicwave solution is
approximately 70% of the kinematicwave solution when
the rainfallkinematicwave 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 kinematicwave solution departs from the
dynamicwave 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 (ifc)
when time reaches infinity. It is noted that the
dynamicwave model has shorter time to peak than that of
the kinematicwave 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, kinematicwave
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 highspeed
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 dynamicwave solution is about 80% of the
kinematicwave solution when the rainfallkinematic 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 rainfallkinematicwave 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 kinematicflow model rather than that of
the dynamicflow model.
Figure 4.3 (withoÂ£= 1.31, and TE = 1198 sec.)
shows that the dynamicwave solution is approximately
88% of the kinematicwave solution when the
rainfallkinematicwave number, K^, is equal to 65.
Figure 4.4 (withd = 1.4 and TE = 1,158 sec)
shows that the dynamicwave solutionis approximately 90%
of the kinematicwave solution when the rainfall
kinematicwave 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 rainfallkinematicwave number, K^, and the peak
ratio of water surface depth, Y/YE, is demonstrated.
As the term, Y/YE, approaches 1.0, the kinematicwave
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
Kinematicwave solution
Figure 4.5 Variation of the Peak Depth Ratio vs. RainfallKinematic
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 rainfallkinemaitcwave number increases,
the ratio of Y/Ys becomes higher and approaches 1=0
when the rainfallkinematic 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 rainfallkinematicwave number increases,
the point of separation where the dynamicflow 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 rainfallkinematicwave 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 rainfallkinematicwave
number, Kj_, to the peak ratio are summarized and shown
in Figure 4.6. The figure shows that when the
rainfallkinematicwave number is greater than 180, the
kinematic wave reaches almost 100% agreement with the
dynamicwave model.
It is noted that when the value of the rainfall
kinematic parameter, Kj_, is smaller than 50, the
73
kinematicwave 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 dynamicwave model merges closely with the
kinematic model and will give results within 0.5% of the
dynamicwave model (point A).
Another paramter, the dynamicwave 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
kinematicwave number) x (1/SQ) x YE/2L. The study
showed that when the dynamicwave number is 38, the
dynamicwave solution is approximately 70% of the
kinematicwave solution. When the dynamicwave numbers
are 33, 27, and less than 6, the dynamicwave solutions
are 80%, 90%, and 100%, respectively.
The results are depicted in Figure 4.7. It
shows that when the dynamicwave number is less than 6,
the kinematicwave approximation will give results
within 5% of the dynamicwave solutions. The parameter,
D, can also be used as a complementary criterion to
check when the dynamicwave 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 kinematicwave routing. The results show that as
the dynamicwave 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 rainfallkinematicwave 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 dynamicwave and kinematicwave
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
kinematicwave 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 kinematicwave number, Kq, in
Equation (1.2), is equal to the value of equilibrium
water depth (YE) used in calculating the rainfall
kinematicwav 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 watersurface depth), qE (equilibrium
flow rate) and Kq (kinematicwave number), and K^
(rainfallkinematicwave 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 kinematicwave 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
KinematicWave Number, Kg 25012 15257 1832 2042
RainfallKinematic 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
waverouting 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 rainfallkinematicwave 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
kinematicwave assumptions are adequate and acceptable.
When the rainfallkinematicwave number is greater than
180, the kinematicwave approximation will give results
within 5% of the dynamicwave approximation.
2. When the rainfallkinematicwave 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 rainfallkinematic
wave number, K^, increases. This indicates that the
flow regime behaves closer to the kinematic model when
increases.
3. Another parameter, the dynamicwave 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 kinematicwave
assumptions are adequate and acceptable. It was
concluded that when the dynamicwave number is less than
6, the kinematicwave approximation will give results
within 5% of the dynamicwave approximation.
4. The testing of the kinemaitcwave model by
laboratory data shows good agreement between this study
and that of Woolhiser and Liggett (1967). It is
concluded that the kinematicwave 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, (rainfallkinematicwave
number), can readily be calculated for a given
watershed. The parameter used by Woolhiser and Liggett
(1967) Kg (kinematicwaver number) cannot be
calculated easily for a given watershed. This study
offers a much better and more complete dynamicwave
modeling.
Further study could add several features to the
routing model. The comparison of kinematicwave routing
to full dynamicwave 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. 24812500 (1970).
Basco, D. R.
"Introduction to Numerical Method"
Verification of Mathematical and Physical Models in
Hydraulic Engineering, ASCE, p. 280302 (1978).
Beutner, E. L., Gaebe, R. R., and Horton, R. E.
"Sprinkledplat Runoff and Infiltrationexperiments on
Arizona Desert Soils"
Trans. Amer. Geophys. Union, p. 550558 (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"
McGrawHill Book Company, New York (1960).
Chen, C. W.
"Computer Simulation of Urban Storm Water Runoff"
J. Hyd. Div., ASCE Hy 2, p. 289301 (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"
McGrawHill 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. 344346,
McGrawHill 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. 918926 (1973).
Gburck, W. J., and Overton, D. E.
"Subcritical Kinematic Flow in a Stable Stream"
J. Hyd. Div. ASCE, HY9, p. 14331446 (1973).
Grace, R. A., and Eagleson, P. S.
"The Modeling of Overland Flow"
Water Res. Research, p^ 393403 (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. 124 (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.
485492 (1984A).
GUO, C. Y.
"Fluid travel time"
Water for Resources Development, Proceedings of ASCE
Hydraulic Conference in Coeur d'Alene, Idaho, p.
380384, 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. 15311540 (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. 340349 (1938).
Horton, R. E.
"Approach toward a Physical Intterpretation of
Infiltration Capacity"
Proc. Soil Sci. Soc. Amer., v. 5, p. 399441 (1939).
Isaacson, E., et al.
"Numerical Solution of Flood Prediction and River
Regulation Problems: Report III"
Report No. IMNNYU235, 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. 129146 (1946).
Kersten,l R. D.
"Engineering Differential Systems"
McGrawHill Book Company, New York (1969).
Kuelegan, G. H.
"Spatially variable discharge over a sloping plane"
Trans. Amer. Geophysical Union, Part VI, p. 956958
(1944) .
Labadie, J. W., Morrow, D. M., and Lazaro, R. C.
"Urban stormwater control package for automated
realtime 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 ShallowWater Equation"
Proc. ASCE J. Mech. Div. 9B (EM2), p. 3971 (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. 281316
(1955).
Lindsley, R. A.
"Hydrology for Engineering"
McGrawHill Book Company, New York (1958).
Meadows, M. E.
"Finite element simulation of overland flow"
Advances in irregation and drainage, surviving external
pressure, AASCE, p. 466475 (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. 441460 (1970).
Morgali, J. R., and Linsley, R.
"Computer Analysis of Overland Flow"
J. Hyd. Div. ASCE, HY3, p. 81100 (1965).
Musgrave, G. W.
"How Much of the Rain Enters the Soil"
Yearbook of Agriculture, p. 151159 (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. 215296 (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.
353363 (1978).
90
Preissman A.
"Propagation des intumescences dans les cannaux et
rivieres"
1st Congress de 1'Assoc. Francaise de Calcul, Grenoble,
p. 443442 (1960).
Price, R. K.
"Comparison of four numerical metnods for flood routing"
Journal of the Hydraulic Division, ASCE, v. 100, no,.
HY7, p. 879899 (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. 429446
(1972) .
Scheaffer, J. R., et al.
"Urban Storm Drainage Management"
Marcel Dekker Inc., New York (1976).
Simons, D. B., Chen, Y. H., and SaezBenito, J. M.
"A mathematical model of Pools 5 through 8 in the Upper
Mississippi River System"
Colorado State University Report CER 7989, DBSYHXJMS
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. 3156 (1962)
Wooding, R. A.
"Hydraulic Model for the Catchment"
J. Hydrology, 3, p. 254267 (1965).
Woolhiser, D. A., and Liggett, J. A.
"Unsteady, OneDimensional Flow Over a Plane The
Rising Hydrograph"
Water Res. Research, 3(3), p. 753771 (1967).
Yen, B. C.
"Methodologies for Flow Prediction in Urban Storm
Drainage System"
Research Report No. 72, University of Illinois,
ChampaignUrbana, Water Res. Center, January (1973).
Yen, B. C.
"Storm Sewer System Design"
University of Illinois, ChampaignUrbana, Illinois
(1978) .
Yevjevich, V. and Barnes, A. H.
"Flood Routing through Storm Drains"
Colorado State University, Hydrology Papers, p. 4346
(1970).

PAGE 1
AND KINEMATICS OF OVERLAND FLOWS by ShiuMing 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 ShiuMing Huang has been approved for the Department of Civil Engineering by 7
PAGE 3
Huang, ShiuMing (M.S., Civil Engineering) Dynamics and Kinematics of Overland Flows iii Thesis directed by Assistant Professor James c. Y. Guo The dynamicwave model for an overland flow was assessed using the implicit central finitedifference method. Flow equations were derived from the continuity and momentum principles. This study attempts to delineate the limitations of kinematicwave approximation by_performing a series of numerical studies using the flow equations. The nondimensional dynamicwave model was combined with Horton's infiltration formula and numerically solved to establish a model representing twodimensional overland flow over a pervious surface. The comparisons between kinematicwave solutions and dynamicwave solutions provide general guidelines anq further establish new criteria for determining when kinematicwave 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 dynamicwave 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, ChiI, 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 DynamicWave Model and viiPAGE 60 61 KinematicWave Nodel o 63 v. 4o3 Evaluation of KinematicWave 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 FiniteDifference 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 RainfallKinematicwave Number 71 4.6 Effects of RainfallKinematicWave Number on the outlet Flow Depth . . 74 4.7 Effects of DynamicWave 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 waterresources 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 highspeed digital computers, numerical schemes have been developed and used to obtain discrete solutions for complicated openchannel flow situations. Prior to the use of highspeed computers, various simplified forms of the st. Venant equations were used to approximate different hydrological problems. The kinematicwave approximation is the most common simplification of the st. Venant equation. It assumes that the friction slope, Sf1is balanced by the channelbottom slope, s0 but it ignores the pressure and inertial terms. Lighhill and Whitham (1955) have found that the floodwave propagation in overland flow is mainly kinematic; Wooding (1965) applied the kinematicwave approximation to solve the partialdifferential governing equations describing the flow for a vshaped chatchment flow. His solutions represented complete runoff hydrographs and the time of concentration. The general solution for the kinematicwave 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 kinematicwave 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 kinematicwave theory is the simpliest form of the complete dyanamicwave model. With the simplification, kinematicwave 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 kinematicwave theory has its limitation. Woolhiser and Liggett (1967) determined the limitations of kinematicwave theory on the basis of a
PAGE 13
series of computer simulation studies of the dynamicwave model. They concluded that for a given uniform rainfall excess, ie, the applicability of kinematicwave theory may be best judged by the kinematicwave 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 development of kinematicwave 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 kinematicwave approximation with Horton's infiltration formula to analytically solve the propagation of a twodimensional 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 kinematicwave approximation limitations under a nonuniform rainfall excess due to infiltration decay characteristics.
PAGE 15
6 This study attempts to determine the limitations of kinematicwave approximation by performing a numerical study which combines the dynamicwave model with Horton's infiltration formula to numerically solve the propagation of twodimensional overland flow on a pervious surface. The results from comparisons between kinematicwave solutions and complete dyanamicwave solutions can provide general guidelines, and can further establish new criteria for determining when kinematicwave 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 dynamicwave equations representing overland flow resulting from a uniform rainfall and a variable infiltration rate. The discrete computational solutions were then compared with the kinematicwave solutions developed by Guo (1984A and 1984B) in order to delineate the limitations of the applicability of kinematicwave solutions. The kinematicwave method is an approach which has been proven to be a very useful tool in stormwater 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 kinematicwave 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 kinematicwave 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 finitedifference 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 overlandflow 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 overlandflow hydrograph were numerically tested. They proposed a semitheoretical 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 kinematicwave solution gives very accurate results. 9 Morgali (1970). did_ a followup study on laminar and turbulent overland hydrographs. In setting up the solution of finitedifference 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 overlandflow 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 overlandflow plane. Giburek and Overton (1973) studied subcritical kinematic flow in a stable stream. An evaluation of the shallowwater 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 kinematicwave theory with
PAGE 19
Horton's infiltration model. The accuracy of solution is confirmed by numerical computations based on an implicit finitediffernece 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 decayinfiltration model to the unsteady dynamic wave routing and to further investigate the limitations of the kinematicwave approximation.
PAGE 20
11 CHAPTER II ANALYTICAL MODEL 2.1 Introduction The governing equations of a dynamicwave model for an overland flow are derived from the continuity and momentum principles. The equations were first applied to the overlandflow problem by Keulegan (1944), and later Morgali and Linsley (1965) did a study on computer analysis of overland flow. The schematic illustration of twodimensional flow is shown in Figure 2.1. The accumulation of runoff on a unitwidth 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, openchannel 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 twodimensional 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 twodimensional 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 state3 that F Til a = DU !\1Dt F = sum of all externa.l :forces in a given direction (in cu.!" ca:;e, the flow direction) (lb) m = mass ( slu'.j} c. acceleration in t:he flcw dir'":lic!:ion U : crosssectional 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 extarnal 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 unitwidth 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 uniformflow 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 crosssectional 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 SaintVenant 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 analytical 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 watersurface 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 kinematicwave 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 kinematicwave 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 kinematicwave approximation cannot model the floodwave attenuation,
PAGE 32
23 and it is not adequate when the backwater effects are significant. The kinematicwave model is accurate for a slowly rising flood wave on a steep slope such as overland flow (Ponce et al., 1978). The diffusivewave approximation assumes that DU/Dt = o and is often called the equation of gradually varied flow. Gradually varied flow is a steadyflow situation where the depth gradually varies along the length of the channel. The improvement of diffusivewave 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 quasisteady dynamicwave model neglects the local unsteadyvelocity term, = o. This model was found to be valid to describe a long floodwave 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 dynamicwave 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 dynamicwave 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 uniformflow 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 wellgraded fines, the infiltration is low. If the soil has several layers of horizons, the leastpermeable 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 semiarid 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 soilair 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 Ei ...: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 (ifc), 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 overlandflow 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, ifc), 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 (if ) L c ; ie = e (if ) c n2 (if )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 (ifc)/(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 (12) 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 straightline 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 downstream boundary condition was not necessary. Noncritical conditions gave good results for the shallowwater 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 highspeed digital computers, however, a number of very efficient numerical methods have been developed in the application of solutionsolving 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 finitedifference 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 twodimensional secondorder 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 twodimensional secondorder equation can be categorized as twodimensional elliptic equations, parabolic and hyperbolic equations. The best known among the twodimensional 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 steadystate 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 FiniteDifference Method Finitedifference 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 centraldifference approximation. 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 "backwarddifference 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 i1 Figure 3.1 AX/2 /2 i i+l X Illustration of FiniteDifference 44
PAGE 54
45 3.3 CentralDifference Method From Equation (3.8), it is clear that the forwarddifference 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): .6ft(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 centraldifference 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 finitedifference 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 finitedifference 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' i1 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. Trialanderror 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 finitedifference 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 centraldifference method provides stable solutions for an unsteady overland flow. The implicit centralfinite 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 fourpoint implicit scheme (Preissman, 1960) is adopted. 50 In this study, the finitedifference 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 centraldifference method, the general finitedifference 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 xaxis, IE j station number along = yaxJ.s, j JE IE = last station on xaxis, which the spatial axis JE = last station on yaxis, which the time axis Rearranging Equation (3.14), a more general form is obtained: 51 is is
PAGE 61
l1 = l1 (j,j1) + (.O.t/26x) [j2 ( i1, j) /2 ( i + 1, j) ] + (.6.t/ 2) [fo3 ( i 1 j ) + /3 (i,j1)] ( 3 .15) Substituting the dynamicwave 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,j1) + [q(il,j) q(i+l,j)] [ie(i,j) + ie(i,j1) ( 3. 16) q(i,j) = q(i,j1) +(.c:.t/2Ax) [ (q(i 1,j)jy(i1,j) + gy(i1,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,j1)] (SoSf) (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 (j1) 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,j1) + 26X* [ q*(il,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 shallowwater flow on a subcritical slope. Using Equations (2.25) and (2.27), the downstream boundary conditions are thus formulated by a backward finitedifference method as follows:
PAGE 64
and Y*(IE,j) q*(JE,j) = Y* (IE,j1) [q*(IEl,j) AX* q*(IE,j)] + + ie* (IE,j1)] = qk(IE,j1) Axx P(IE,j)] +.P.t* [R (IE,j) + 2 R(IE,j1)] (3.22) (3.23) The stability and accuracy, plus the con55 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 finitedifference meth6d is used because it gives the best approximate results among the forward, backward, and centraldifference 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 watersurface profile at current time step. The new watersurface 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 (jl)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 if ie = t solut1ons of prev1ous 1me step as initial ess Sweep the dis char within 103 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 channelflow conditions. As the kinematic wave formulation neglects all terms in the momentum equation except the friction slope and the channelbottom 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 rainfallkinematicwave number, Ki, and the dynamicwave 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 (ifc)/(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 kinematicwave and dynamicwave 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 DynamicWave and KinematicWave Models Figures 4.1, 4.2, 4.3, and 4.4 show the effects of the rainfallkinematicwave number, Ki, on the peak ratio of flow depth, Y/YE (as defined in Section 2.5). The kinematicwave solutions of each case are plotted as the solid line, and the dynamicwave 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 dynamicwave solution is approximately 70% of the kinematicwave solution when the rainfallkinematicwave 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 kinematicwave 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 . Kinematicwave solution Dynamicwave solution 6k__ ..,.__ . 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 (ifc) when time reaches infinity. It is noted that the dynamicwave model has shorter time to peak than that of the kinematicwave model. Inthe 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, kinematicwave 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 highspeed 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 dynamicwave solution is about 80% of the kinematicwave solution when the rainfallkinematic 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
Kinematicwave solution 1.0 Y/YE Dynamicwave 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 kinematicflow model rather than that of the dynamicflow model. Figure 4.3 (withcl= 1.31, and TE = 1198 sec.) shows that the dynamicwave solution is approximately 88% of the kinematicwave solution when the rainfallkinematicwave number, Ki, is equal to 65. Figure 4. 4 (with ol. = 1. 4 and TE = 1, 158 sec) shows that the dynamicwave solutionis approximately 90% of the kinematicwave solution when the rainfallkinematicwave 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 rainfallkinematicwave number, Ki, and the peak ratio of water surface depth, Y/YE, is demonstrated. As the term, Y/YE, approaches 1.0, the kinematicwave. 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 Kinematicwave solution Y/Y E _k_ Dynamicwave 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 ..... Dynamicwave 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 Kinematicwave solution Y/YE .. tJ
PAGE 81
72 (1) When the rainfallkinemaitcwave number increases, the ratio of Y/YE becomes higher and approaches 1.0 when the rainfallkinematic 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 dynamicflow 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 rainfallkinematicwave number, Ki, to the peak ratio are summarized and shown in Figure 4. 6. The figure sho"t.tTS that when the rainfallkinematicwave number is greater than 180, .the kinematic wave reaches almost 100% agreement with the dynamicwave model. It is noted that when the value of the rainfallkinematic parameter, Ki, is smaller than 50, the
PAGE 82
73 kinematicwave 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 dynamicwave model merges closely with the kinematic model and will give results within 0.5% of the dynamicwave model (point A). Another paramter, the dynamicwave 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 (rainfallkinematicwave number) x (1/50 ) x YE/2L. The study showed that when the dynamicwave number is 38, the dynamicwave solution is approximately 70% of the kinematicwave solution. When the dynamicwave numbers are 33, 27, and less than 6, the dynamicwave solutions are 80%, 90%, and 100%, respectively. The results are depicted in Figure 4.7. It shows that when the dynamicwave number is less than the kinematicwave approximation will give results within 5% of the dynamicwave solutions. The parameter, D, can also be used as a complementary criterion to check when the dynamicwave 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 RainfallKinematic 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 kinematicwave routing. The results show that as the dynamicwave 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 rainfallkinematicwave 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 dynamicwave and kinematicwave 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 kinematicwave 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 kinematicwave number, Kq, in Equation (1.2), is equal to the value of equilibrium water depth (YE) used in calculating the rainfallkinematicwavenumber, 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 watersurface depth), qE (equilibrium flow rate) and Kq (kinematicwave number), and Ki (rainfallkinematicwave 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 kinematicwave 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 1time 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) KinematicWave Number, Kq RainfallKinematicwave 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 waverouting 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 rainfallkinematicwave 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 kinematicwave assumptions are adequate and acceptable. Hhen the rainfallkinematicwave number is greater than 180, the kinematicwave approximation will give results within 5% of the dynamicwave approximation. 2. When the rainfallIdnematicwave 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 rainfallkinematicwave number, Ki, increases. This indicates that the flow regime behaves closer to the kinematic model when Ki increases . 3. Another parameter, the dynamicwave 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 dynamicwave number less than 6, the approximation give results within 5% of the dynamicwave approximation. 4. The testing of the kinemaitcwave model by laboratory data shows good agreement .between this study and that of Woolhiser and Liggett (1967). It is concluded that the kinematicwave 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 (rainfallkinematicwave number), can readily be calculated for a given watershed. The parameter used by Woolhiser and Liggett (1967) Kq (kinematicwaver number) cannot be calculated easily for a given watershed. This study offers a much better and more complete dynamicwave modeling. Further study could add several features to the routing model. The comparison of kinematicwave routing to full dynamicwave 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. 24812500 (1970). Basco, D. R. "Introduction to Numerical Method" Verification of Mathematical and Physical Models in Hydraulic Engineering, ASCE, p. 280302 (1978). Beutner, E. L., Gaebe, R. R., and Horton, R. E. "Sprinkledplat Runoff and Infiltrationexperiments on Arizona Desert Soils" Trans. Amer. Geophys. Union, p. 550558 (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" McGrawHill Book Company, New York (1960). Chen, c. w. "Computer Simulation of Urban Storm Water Runoff" J. Hyd. Div., ASCE Hy 2, p. 289301 (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" McGrawHill 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 344346, McGrawHill 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. 918926 (1973). Gburck, w. J., and overton, D. E. 11Subcritical Kinematic Flow in a Stable Stream11 J. Hyd. Div. ASCE, HY9, p. 14331446 (1973). Grace, R. A., and Eagleson, P. s. 11The Modeling of overland Flow11 Water Res. Research, p: 393403 (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. 124 (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. 485492 (1984A). Guo, c. Y. "Fluid travel time" Water for Resources Development, Proceedings of ASCE Hydraulic Conference in Coeur d'Alene, Idaho, p. 380384, 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. 15311540 (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. 340349 (1938). Horton, R. E. "Approach toward a Physical Intterpretation of Infiltration Capacity" Proc. Soil Sci. Soc. Amer., v. 5, p. 399441 (1939). Isaacson, E., et al. "Numerical Solution of Flood Prediction and River Regulation Problems: Report III" 88 Report No. IMNNYU235, 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. 129146 (1946). Kersten,! R. D. "Engineering Differential Systems" McGrawHill Book Company, New York (1969). Kuelegan, G. H. "Spatially variable discharge over a sloping plane" Trans. Amer. Geophysical Union, Part VI, p. 956958 (1944). Labadie, J. w., Morrow, D. M., and Lazaro, R. c. "Urban stormwater control package for automated realtime 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 ShallowWater Equation" Proc. ASCE J. Mech. Div. 9B (EM2), p. 3971 (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. 281316 (1955). Lindsley, R. A. "Hydrology for Engineering" McGrawHill Book Company, New York (1958). Meadows, M. E. "Finite element simulation of overland flow" Advances in irregation and drainage, surviving external pressure, AASCE, p. 466475 (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. 441460 (1970). Morgali, J. R., and Linsley, R. "Computer Analysis of Overland Flow" J. Hyd. Div. ASCE, HY3, p. 81100 (1965). Musgrave, G. W. "How Much of the Rain Enters the Soil" Yearbook of Agriculture, p. 151159 (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. 215296 (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. 353363 (1978).
PAGE 99
Preissman A. "Propagation des intumescences dans les cannaux et rivieres" 90 1st Congress de l'Assoc. Francaise de Calcul, Grenoble, p. 443442 (1960). Price, R. K. "Comparison of four numerical metnods for flood routing" Journal of the Hydraulic Division, ASCE, v. 100, no,. HY7, p. 879899 (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. 429446 (1972). Scheaffer, J. R., et al. "Urban Storm Drainage Management" Marcel Dekker Inc., New York (1976). Simons, D. B., Chen, Y. H., and SaezBenito, J. M. "A mathematical model of Pools 5 through 8 in the Upper Mississippi River System" Colorado State University Report CER 7989, DBSYHXJMS 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. 3156 (1962) Wooding, R. A. "Hydraulic Model for the Catchment" J. Hydrology, 3, p. 254267 (1965). Woolhiser, D. A., and Liggett, J. A. "Unsteady, OneDimensional Flow Over a Plane The Rising Hydrograph" Water Res. Research, 3(3), p. 753771 (1967). Yen, B.. c. "Methodologies for Flow Prediction in Urban Storm Drainage System" Research Report No. 72, University of Illinois, ChampaignUrbana, Water Res. Center, January (1973). Yen, B. c. "Storm Sewer System Design" University of Illinois, ChampaignUrbana, Illinois (1978). Yevjevich, v. and A. H. "Flood Routing through Storm Drains" Colorado State University, Hydrology Papers, p. 4346 (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
