EPA/600/B-18/256 | August 2018 |
United States
Environmental Protection
Agency
Flow Routing Techniques for
Environmental Modeling
-------
EPA/600/B-18/256
August 2018
Flow Routing Techniques for
Environmental Modeling
by
Jan Sitterson1, Chris Knightes2, Brian Avant3
1Oak Ridge Associated Universities
2 U.S. Environmental Protection Agency, Office of Research and Development
3 Oak Ridge Institute for Science and Education
Chris Knightes - Project Officer
Office of Research and Development
National Exposure Research Laboratory
Athens, GA, 30605
l
-------
Notice
The U.S. Environmental Protection Agency (EPA) through its Office of Research and
Development funded and managed the research described here. The research described herein
was conducted at the Computational Exposure Division of the U.S. Environmental Protection
Agency National Exposure Research Laboratory in Athens, GA. Any mention of trade names,
products, or services does not imply an endorsement by the U.S. Government or the U.S.
Environmental Protection Agency. The EPA does not endorse any commercial products,
services, or enterprises. This document has been reviewed by the U.S. Environmental Protection
Agency, Office of Research and Development, and approved for publication.
2
-------
Abstract
This report describes a few ways to simulate the movement of water through a network of
streams. Flow routing connects excess water from precipitation and runoff to the stream to other
surface water as part of the hydrological cycle. Simulating flow helps elucidate the transportation
of nutrients through a stream system, predict flood events, inform decision makers, and regulate
water quality and quantity issues. Three flow routing techniques are presented and discussed in
this paper. Constant volume and changing volume techniques use the continuity equation, while
the third technique, the kinematic wave approximation, uses the continuity and momentum
equations. Inputs and outputs differ for each flow routing technique, with kinematic wave being
the most complex model. Each option has a set of applications it is best suited for. Numerical
errors such as distortion and instability are potential errors that can affect the accuracy of model
outputs. Routing methods should be chosen based on input data availability and scope of the
problem being addressed. This report informs modelers of algorithms, input data needed, and the
applicability of a few flow routing techniques used in environmental modeling.
3
-------
Foreword
The U.S. Environmental Protection Agency (EPA) is charged by Congress with protecting the
Nation's land, air, and water resources. Under a mandate of national environmental laws, the
Agency strives to formulate and implement actions leading to a compatible balance between
human activities and the ability of natural systems to support and nurture life. To meet this
mandate, EPA's research program is providing data and technical support for solving
environmental problems today and building a science knowledge base necessary to manage our
ecological resources wisely, understand how pollutants affect our health, and prevent or reduce
environmental risks in the future.
The National Exposure Research Laboratory (NERL) Computational Exposure Division (CED)
develops and evaluates data, decision-support tools, and models to be applied to media-specific
or receptor-specific problem areas. CED uses modeling-based approaches to characterize
exposures, evaluate fate and transport, and support environmental diagnostics/forensics with
input from multiple data sources. It also develops media- and receptor-specific models, process
models, and decision support tools for use both within and outside of EPA.
The goal of the Hydrologic Micro Services (HMS) project is to develop an ecosystem of inter-
operable water quantity and quality modeling components. Components are light-weight and can
be integrated to rapidly compose work flows to address water quantity and quality related
questions. Each component may have multiple implementations ranging from macro (coarse) to
micro (detailed) levels of modeling the physical processes. The components leverage existing
internet-based data sources and sensors. They can be integrated into a work flow in two ways:
calling a web service or downloading component libraries. For light-weight components, it is
generally more efficient to call a web service, however, it is more efficient to have local copies
of components if the component requires large amounts of input/output data.
Tom Pierce, Acting Division Director for CED
4
-------
Table of Contents
Notice 2
Abstract 3
Foreword 4
List of Figures 5
List of Tables 6
List of Equations 6
Introduction 7
Stream Routing Options 9
Constant Volume Routing 11
Changing Volume Routing 12
Kinematic Wave Routing 13
Channel Geometry 17
Potential Errors 19
Discussion 20
References 20
Appendix A 22
List of Figures
Figure 1: A simplified diagram of the water cycle within a watershed 7
Figure 2: The momentum terms and equation used for each model 9
Figure 3: Graphs of flow moving through instantaneous and non- instantaneous models 11
Figure 4: Instantaneous flow schematic 12
Figure 5: The flow-depth regression graph 13
Figure 6: The computational grid for solving the finite difference numerical approximation 15
Figure 7: Stream location visualization 16
Figure 8: Channel Cross Sectional Geometry Descriptions 18
Figure 9: Outputs from models that show numerical errors 19
5
-------
List of Tables
Table 1: Input parameters for flow routing schemes 10
Table 2: Outputs for flow routing schemes 10
Table 3: Table of formulas for channel shapes 18
List of Equations
Equation (1) 8
Equation (2) 8
Equation (3) 11
Equation (4) 11
Equation (5) 12
Equation (6) 12
Equation (7) 13
Equation (8) 13
Equation (9) 13
Equation (10) 14
Equation (11) 14
Equation (12) 14
Equation (13) 14
Equation (14) 15
Equation (15) 15
Equation (16) 16
Equation (17) 16
Equation (18) 17
Equation (19) 17
Equation (20) 19
6
-------
Introduction
As part of the hydrological cycle, stream flow plays an important role in the movement of
water through surface water systems by connecting precipitation and runoff to downstream water
bodies, like lakes and oceans (Perlman, 2016). Figure 1 shows the water cycle within a
watershed and the major interconnected components such as precipitation, runoff, soil water,
evapotranspiration, groundwater, and streamflow. A watershed is a topographic area where all
surface waters and shallow groundwaters drain to a specific point (Edwards et al., 2015). Figure
1: A simplified diagram of the water cycle within a watershed. (Brewster, 2017; ESRI, 2017)As
precipitation falls onto the surface it can evaporate, be taken up by vegetation, infiltrate into the
ground, and create surface and subsurface runoff to surface waters such as streams. These stream
systems allow water to travel long distances and connect to larger watersheds or bodies of water.
Surface water systems have great importance because they provide drinking water, electricity,
and recreation for people. Being able to model this movement allows one to predict and trace the
flow of water through a hydrologic system (Chow et al., 1988; O'Sullivan et al., 2012).
Simulating flow throughout a stream network of a watershed, not only at a pour point of a
watershed, is essential for studying water quality and quantity issues throughout the watershed.
Determining the flow rate through a stream segment is a central task of surface water hydrology
studies looking at erosion, pollution concentrations, and flooding occurrences, among many
things (Chow et al., 1988).
Soil Moisture
Infiltration
The Water Cycle
Precipitation
Evapotranspiration
Figure 1: A simplified diagram of the water cycle within a watershed. (Brewster, 2017; ESRI,
2017)
1
-------
Flow in an open stream network is governed by the continuity equation, in which mass is
conserved, and the momentum equation, in which forces are at equilibrium.
Continuity Equation:
^ + ^£=n (1)
dx at
Where (?is discharge [m3/s], Ac is cross-sectional area [m2], fis time [s], and xis length [m]. The
continuity equation states that the change in discharge over a set length of stream must offset the
change in cross-sectional area over a set time to maintain the same amount of mass pushing
through the system. In the conservation of mass, or continuity equation, volumetric flow rates
between channel cross sections are equal since the volume of the channel is unchanged over time
and mass is neither created or destroyed (Chaudhry, 2008).
Momentum Equation:
dy v dv 1 dv (2)
^ ^0 dx g dx g dt
Where 5"ois bed slope, 5>is frictional slope, y\s surface elevation [m], xis length [m], v\s water
velocity [m/s], g\s gravitational acceleration [m/s2], fis time [s]. The names of grouped variables
are the pressure term rr, the convective acceleration term ancj the local acceleration term
ox g ox
The momentum equation describes the equilibrium of forces acting on a fluid (UCAR,
2006). The conservation of momentum equation is based on Newton's second law of motion
where the rate of change in momentum of liquid volume is equal to the result of external forces
acting on the liquid volume (Chaudhry, 2008).
Simple flow models use the conservation of mass for an inflow-outflow relationship.
More complex models use a physics-based method that incorporates the conservation of mass
and conservation of momentum equations (Chow et al., 1988; UCAR, 2006). The kinematic,
diffusive, steady dynamic, and full dynamic wave models incorporate both equations at different
complexities to solve multiple streamflow conditions. Figure 2 shows how the momentum
equation can be broken down into specific terms for several modeling purposes. The kinematic
wave approximation assumes the frictional and bed slope forces dominate, and the rest of the
terms are assumed to be negligible. In the diffusive wave approximation, the pressure term is
included but the convective and local acceleration is negligible. The approximation for the steady
dynamic wave model includes all the terms except the local acceleration term. The full dynamic
wave approximation uses all terms in the full momentum equation.
8
-------
Kinematic
friction slope=Channel slope
Diffusive
Includes press jre term
Steady Dynamic
r vp _ v vv
f "0 Sx g Sx
dy v 3v
Includes convective acceleration
Foil Dynamic
« « i/;. i.' yy jl vt'
f 0 Iix g 8x g dt
includes local acceleration
dy v Ov 1Ov
Figure 2: The momentum terms and equation usedfor each model.
Stream Routing Options
Three flow routing methods for steady uniform flows are described here in order of
increasing complexity: constant volume, changing volume, and the kinematic wave routing
method.
1. Constant volume routing is a simple method using only the continuity equation where
flow going out of the segment is equal to the flow going into the segment. Volume,
velocity, and depth remain constant with changing flows in this scheme.
2. Changing volume routing is another conservation of mass method in which, flow coming
in is equal to the flow going out of the segment but volume, velocity, and depth change
with flow based on channel hydrogeometry.
3. The kinematic wave model describes one-dimensional steady uniform flow waves
moving through a stream channel. It is a simplification of the continuity and momentum
The constant volume and changing volume methods are easily calculated by using only a
simplified continuity equation. They are suitable for preliminary estimates of the time and shape
of a flood wave (O'Sullivan et al., 2012). Methods for simulating flow with the continuity and
the momentum equations, like the kinematic wave option, describe the temporal and spatial
variations in flows throughout the stream network but at a computational cost. Table 1 shows the
different input parameters needed for each of the stream routing options, while Table 2 shows the
outputs of each model. Increasing model complexity requires increased numbers of parameters
and computation costs. Figure 2 shows the differences in the hydrographs produced by the
different methods described above. See Appendix A for the source code used in running these
flow routing options.
equations.
9
-------
Table 1: Input parameters for flow routing schemes.
Inputs
Boundary
Flow
[m3/s]
Length
[m]
Bottom
Width
[m]
Z-
slope
[2]
Depth
Multiplier
y
Depth
Exponent
c
Manning's
Roughness
n
Channel
Slope
So
Initial
Depth
[m]
Constant
Volume
V
Changing
Volume
V
V
V
V
Kinematic
Wave
V
V
V
V
V
S
*Z-slope needed for triangle and trapezoid channel geometries
Table 2: Outputs for flow routing schemes
Outputs
Outflow [m3/s]
Velocity [m/s]
Volume [m3]
Depth //;//
Constant
Volume
V
Changing
Volume
V
V
V
V
Kinematic
Wave
Both the constant volume and the changing volume routing methods are examples of
instantaneous flow models. This means the entire system responds immediately to a change in
upstream inflow. Figure 3 A shows the output of an instantaneous model with no time delay in
movement of the flow wave downstream. The kinematic wave routing option is a non-
instantaneous flow model. Figure 3B shows the model output where the flow wave is translated
through the stream system at different times. When there is an increase in flow at 20 minutes
upstream, 1500 meters downstream has not seen the implications of this increased flow yet. This
time delay in flow is more representative of real world processes.
10
-------
Instantaneous
Non-Instantaneous
20 ¦
x=500m
x= 1000m
x= 1500m
20 ¦
60 ¦
x=500m x= 1000m x=15Q0m
0
20 <10 60
lime (min)
BO
20 40 60 BO
Time (min)
Figure 3: Graphs showing how flow moves through instantaneous and non- instantaneous
models. A) The model with instantaneous wave propagation for all segment location and B) A
model showing a time delay in flow downstream.
Constant Volume Routing
Constant volume flow routing is the simplest of the three options presented. Equation 3
shows the application of the conservation of mass where outflow (Qout) °f the stream segment
of interest is equal to the inflows (Qin) into the stream segment of interest. Constant volume
routing is an instantaneous model so any changes in inflows anywhere translates throughout the
entire system without a time lag (see Figure 3 A). Volume, velocity, width, and depth remain
constant. To work within a stream network, outflow of a segment of interest is equal to
the inflow from the segment before (Qiit) and the sum of any boundary flows (Qs,t) int0 the
segment of interest (Equation (4)).
Qout ~ Qin
(3)
n
(4)
11
-------
Boundary 1
<3si I
Qj Qi = Qbi
Ql = Qb2
Q3 = Qi + Q2
Qa, = Q3= Qi + Qz
Figure 4: The schematic above shows how incoming flows are instantaneous for all stream
segments within a flow path.
Changing Volume Routing
For the changing volume routing option, flow is propagated instantaneously as with
constant volume routing using Equations (3) and (4). However, in this method volume, velocity,
area, and depth change as flow changes. To solve the continuity equation a relationship must be
made between the two unknowns; flow (Q) and cross-sectional area (Ac). Flow [m3/s] can be
described as the product of the cross-sectional area [m2] of the stream channel and the velocity
(v) [m/s] of water moving through the system.
Q = vAc (5)
As flow changes, the velocity and/or depth must change to maintain mass balance; the
velocity is therefore inversely proportional to the cross-sectional area. A flow-depth regression
model based on the log of empirical observations is used to find depth assuming the power
relationship below (Leopold; & Maddock, 1953).
d = yQc (6)
Where t/is depth [m], (?is discharge [m3/s], _Kand care constants calculated by the regression
model using observed data (Equation (5)). Taking the logarithm of Equation (6), gives:
Boundary 2
12
-------
log(,d) = c * log(Q) + log(y) (7)
Figure 5 shows an example of the regression model. Volume, velocity, depth, and area of the
segment are then calculated based on the geometry of the stream. Ac is calculated using equations
provided in Table 3. Velocity (v) is calculated as given by Equation (8), and volume (V)is
calculated as given in Equation (9).
Q
V = TC
V = L*AC
0.75
0.50
¦o
o> 0 25
o
0.00
-0.25
Figure 5: The flow-depth regression gi'ciph to solve for the slope, c, and the y-intercept, log(y).
Kinematic Wave Routing
The kinematic wave routing scheme does not assume that flow changes instantaneously
throughout the system. This routing option uses the Continuity and Momentum equations often
called the Saint Venant equations when used together (Miller, 1984). The kinematic wave
method adopts the idea that the natural movement of flood waves is governed mostly by the
friction and bed slopes. Kinematic wave models assume the pressure term, the convective
acceleration, and the local acceleration terms of the Momentum Equation to be insignificant, thus
only gravity and friction forces are included (Chapra, 1997; Miller, 1984). Some general
13
(8)
(9)
1.0
1.5
2.0
2.5
log Q
-------
assumptions for using the kinematic wave flow routing technique are that the flows are one-
dimensional with no backwater effects, and the velocity is constant during the model time step
(MacArthur & Devries, 1993). The changing flow wave propagates downstream, allowing an
increased flow at the boundary to travel down the stream network with a lag in time (M.J.
Lighthill & Whitham, 1955; Miller, 1984). The kinematic wave routing scheme translates a flood
wave through the system with no attenuation of the peak flow (Nwaogazie, 1986). The resulting
changes in flow, velocity, width, and depth associated with the flow wave are calculated in this
routing option (Ambrose & Wool, 2017). The derivation of the kinematic wave equation is
presented below.
Continuity Equation:
dQ_ £^c = Q (10)
dx dt
Where Q is discharge [m3/s], Ac is cross-sectional area [m2], fis time [s], and xis segment length
[m]. In this equation we have two unknowns, Qand Ac; therefore, a second equation is
needed. This second equation is obtained by combining, the Momentum Equation and
Manning's Equation (described later).
Momentum Equation:
dy v dv Idv (11)
^ ^0 dx g dx g dt
Qy fly
Where 5>is frictional slope, 5^is bed slope, — is the pressure term, is the convective
acceleration term, and is the local acceleration term. Since the pressure, convective, and
acceleration terms are assumed to be insignificant, the Momentum Equation simplifies to Sf =
So-
Manning's Equation:
q = ±aJJ^ ('2
Where n is the Manning's roughness number, Ac is cross-sectional area [m2], and 5>is the
frictional slope of the channel. Because the frictional slope is equal to bed slope in the
Momentum Equation, Manning's Equation can be used to solve for one of the unknowns from
the Continuity Equation. Rearranging Manning's Equation to solve for area, with the hydraulic
radius substituted as R = Ac/P and Sf = S0.
Ac =
2 \ 5
nPzQ
yfSo
(13)
14
-------
Where /"is the wetted perimeter [m] described in the Channel Geometry section. By comparing
Equation (13) to the power function:
Ac = ocQP
(14)
np2/3
jlj^ / q /r
we can set (3 = 3/5 and a = ( ,— ) ' . Equation (14) is differentiated with respect to time, and
V ^0
the result is inserted as the second term in Equation (10) to obtain an equation with Q as the only
unknown:
dQ R 1 dQ
-^ + aPQp~1^ = 0
dx dt
(15)
To solve this equation, we use the finite difference numerical approximation. There are
many methods for the approximation. An explicit method is described by Chapra (1997), a
partially implicit method is described by Chow et al. (1988), and a fully implicit version is used
in USEPA's Water Quality Analysis Simulation Program (WASP) (Ambrose & Wool, 2017).
We chose to implement a solution like the one presented in WASP. The finite difference method
solves the partial differential equation on a grid representing the space (i) time (t) plane
(Figure 6). An implicit scheme uses the points along the unknown coordinate for a solution.
Location /is the location of interest and t +1 is moving forward in time.
t+1 --
I
¦ —i
H
Qi-l, t+1
Ax
| d
. Q-i, t+i
w %
L
%
Qi,t
_ Qj,t Qi-i.t+i
dx Ax
dQ _ Qi,t+1 ~ Qi,t
dt
At
i-1
Distance
Figure 6: The computational grid for solving the finite difference numerical approximation. The
red dot represents the unknown value for space and time, the black dots represent known values.
15
-------
Using the equations for the discharge variation terms ^ and ^ shown in Figure 6 to solve the
kinematic equation implicitly leaving two unknown values:
{Qi,t ~ Qi-l,t+l\ onP-1 (Qi,t+1 ~~ Qi,t\ _ n (16)
(Vi.t Qi-l,t+l\ , r>r.0-1 /vi,t+1 Vi,t\
( £ J+ (—Zi—J=
The subscript i, t+1 is the location and time we are interested in solving for, thus Q,it+i is the
unknown flow. <2i_ljt+1is the flow from the upstream segment at the time we are interested in
and Qit is the flow from the segment of interest at the previous time step. Figure 7 is a
visualization of the stream segment location.
Ax
t+1
Qii, t+i
Figure 7: Stream location visualization. Qi^+i is the unknown flow and Qi_lt+1is the flow fi'om
the upstream segment.
Rearranging Equation 16 to find the unknown flow:
_ ~(Qi,t ~ Qi-i,t+1) A „ (17)
Ql't+1 ~ Axa/3(Qlt+1)P-i Ql*
This equation cannot be solved algebraically because the unknown, Qiit+1 is on both sides of the
equation in its simplest form. The equation must be solved using an iteration. For the initial
calculation we set one of the unknowns to equal the outflow from the current time and location
(Qi,t+1 = Qi,t) > then the unknown, Qiit+1 is solved in an iterative loop where is equal to
the new answer for Qi t+1arid solved again until it converges with an error of less than 10"3.
16
-------
error =
Qjt+1 ~ Qi,t+1
(18)
Once the solution converges, stream characteristics are updated. Area is calculated by Equation
(19) below. The depth variable is calculated using the equations in Table 3 for the respective
channel geometry. Alpha is then recalculated using the new depth. If the absolute difference
between the two calculated alphas is greater than a 10"3 error, then alpha is updated to the new
value and area is recalculated. Using the updated area; depth, velocity, and volume are
recalculated using equations in Table 3, and Equations 8, 9, respectively.
Natural channels can have a wide range of shapes but models use simplified channel
shapes (Moglen, 2015). Three shapes are presented to represent natural channel geometries;
rectangular, trapezoidal, and triangular (see Figure 8). The shape of the channel segment is used
to relate area, depth, volume, velocity, and flow. Alpha (a) will change with the shape and depth
of the channel due to the wetted perimeter calculation for each channel shape, as presented in
Ac = aQp
(19)
Channel Geometry
Table 3.
17
-------
Rectangular
~
d
w
Trapezoidal
07
<4 ~
W
Triangular
Figure 8: Channel Cross Sectional Geometry Descriptions where d=Depth (m), w = Width (m),
z=run of side slope (default = 2).
Table 3: Table of formulas for channel shapes
Channel
Shape
Cross Sectional
Area
Wetted
Perimeter
Alpha
Rectangular
w * d
2 d + w
( n 2\3/5
te(2d+w)I)
Trapezoidal
w * d + z * d2
w + 2 dV 1 + z2
/ 2\3/s
(wSw+2d^i+z2y)
Triangular
z * d2
2 d-J 1 + z2
18
-------
Potential Errors
Numerical approximations may introduce errors through their calculations. As natural
flood waves travel downstream, they may attenuate and disperse. However, the kinematic wave
solution is not structured to incorporate these processes, it only translates waves through the
system. Due to numerical approximations in the solution of the kinematic wave equation, these
processes may be artificially introduced and create distortion of the flood wave (Nwaogazie,
1986). Numerical distortion is a type of discretization error that occurs due to the numerical
approximations of an implicit solution which varies from the true solution (Fread, 1974).
Distortion can be seen in the model output hydrograph as dispersion of the flood wave or as
dampening (attenuation) of the wave peak. The attenuation of the peak flow is introduced by the
limitations of the numerical solution (Chaudhry, 2008). Another concern when modeling flood
waves is numerical instability. Instability in model output is defined as growth without bound of
numerical errors from truncation and roundoff (Fread, 1974; Strelkoff & Falvey, 1993). These
numerical errors can be amplified during the successive computations which can mask the true
solution or crash the model. Nwaogazie (1986) describes in detail how large time steps can cause
instability and how the chosen explicit or implicit solution scheme can cause attenuation of the
peak flow and dispersion in a hydrograph. Figure 9 shows the resulting hydrograph from
kinematic wave models containing numerical errors. In order for a model to maintain stability
the time step has to meet the Courant Condition (Chapra, 1997, Equation (20)). Many models,
including USEPA's WASP, incorporate code to make sure the Courant Condition is not violated
(Ambrose, 2017). The Courant Condition states that the time step (At) must be less than the time
for a wave to travel the segment distance (Ax).
A.^Ax „ (2°°)
AO -
30 -
1/1
m
E
J 20 -
U_
10 -
o-
Figure 9: Outputs from models that show numerical errors of instability (green), distortion
(orange), and the expectedflow wave (blue).
Ax dQ
At < —, ck = —-
Numerical Errors
ft Expected
JJ Distortion
/
Instability
f
//
/ /
I
1 /
I
// 1
//
jj
JJ J
0 5 10 15 20 25
Time (hrl
19
-------
Numerical errors can significantly affect the accuracy of the solution (Fread, 1974). The
magnitude of the two modeled peak flows differ from the expected flow wave as shown in the
figure. The rising and falling limbs of the modeled hydrographs also diverge from the expected
hydrograph. This can have substantial implications on flood prediction and regulatory decision
making.
Discussion
Modeling flow throughout a stream network provides useful information about the
volume, depth, discharge rate, and velocity in the network. Having detailed information within a
watershed is useful for water quality and quantity studies. The routing methods presented show
the effects of channel geometry and frictional forces on the discharge rate. Modeling flow is
important for planning, designing, regulating, and managing our water resources (Nwaogazie,
1986). Routing methods should be chosen based on input data availability and scope of the
objective. Streamflow routing methods simulate flow moving through a stream network, not just
discharge at a pour point. The volumetric methods presented here require fewer input parameters
than more complicated methods. Among the three methods presented in this report, the
kinematic wave flow routing option simulates the most detail within the stream segment but is
more computationally intensive.
Simple routing schemes like constant and changing volume should be used when the
specific data needed for complex routing schemes are not available (Nwaogazie, 1986). Constant
volume routing models can be applied anywhere to study the movement of water through a
system, because stream characteristics are not needed. Changing volume routing methods are
useful when the model time step is smaller than the time it takes for a flow wave to move
through the system. For example, if the model time step is hourly and the flow wave moves
through the system in one day, specific details from more complex models are unnecessary,
because the volume changes only slightly in the model time step. The complex kinematic wave
routing option is useful when detailed information about the stream system is available and when
travel time of flow is important. The prediction of the a resulting flood wave from a dam failure
is an example of a situation in which the more complex kinematic wave model is needed in
environmental modeling (UCAR, 2006).
References
Ambrose, R., & Wool, T. (2017). WASP8 Stream Transport-Model Theory and User's Guide.
Brewster, R. (2017). Paint.net 4.0.17 (Version 4.0.17). Retrieved from https://www.getpaint.net/
Chapra, S. C. (1997). Surface Water-Quality Modeling: The McGraw-Hill Companies, Inc. .
Chaudhry, M. H. (2008). Open-Channel Flow (Second Edition ed.): Springer.
Chow, V. T., Maidment, D. R., & Mays, L. W. (1988). Applied Hydrology : McGraw-Hill Book
Company.
20
-------
Edwards, P. J., Williard, K. W. J., & Schoonover, J. E. (2015). Fundamentals of Watershed
Hydrology. Journal of Contemporary Water Research & Education, J54( 1), 3-20. doi:
10.1 lll/j,1936-704X.2015.03185.x
ESRI. (2017). ArcGIS Pro Release 2.0. Redlands, CA: Environmental Systems Research
Institute.
Fread, D. L. (1974). Numerical properties of implicit four-point finite difference equations of
unsteady flow. (NWS HYDRO-18). National Oceanic and Atmospheric Administration
Leopold;, L. B., & Maddock, T. (1953). The hydraulic geometry of stream channels and some
Physiographic implications. U.S. Geologic Survey Prof. Pap, 252, 56.
M.J. Lighthill, & Whitham, G. B. (1955). On Kinematic Waves I. Flood Movement in Long
Rivers Kinematic Waves (pp. 281-316). London: Proc. R. Soc.
MacArthur;, R., & Devries, J. (1993). Introduction and Application of Kinematic Wave Routing
Techniques Using HEC-1.
Miller, J. (1984). Basic Concepts of Kinematic-wave Models
Moglen, G. E. (2015). Fundamentals of Open Channel Flow. Boca Raton, FL: CRC Press.
Nwaogazie, I. L. (1986). Comparative analysis of some explicit-implicit streamflow models.
Adv. Water Resources, 10, 69-77. doi: 030917088702006909
O'Sullivan, J. J., Ahilan, S., & Bruen, M. (2012). A modified Muskingum routing approach for
floodplain flows: Theory and practice. Journal of Hydrology, 470, 239-254. doi:
10.1016/j.jhydrol.2012.09.007
Perlman, H. (2016, December 15). The Water Cycle- USGS Water Science School. Retrieved
12/18, 2017, from h ttp s: //water. uses, gov/edu/watercy et e. h tm 1
Strelkoff, T. S., & Falvey, H. T. (1993). Numerical Methods Used to Model Unsteady Canal
Flow. Journal of Irrigation and drainage engineering, 119(4), 637-655.
UCAR. (2006). Streamflow Routing. Retrieved 1/20/2018, from
https://www.meted.ucar.edii/hvdro/basic/Routine/
21
-------
Appendix A
Source code written in Python v3.6 for the flow routing options. The GitHub repository of the
code can be found at
https://github.com/quanted/hms flask/blob/dev/modules/hms/hvdrodvnamics.pv
Constant Volume:
##loop through segments to calculate outflow and parameters for each time
step.
for t in nange((len(time_s))):
for i in nange(nsegments):
Q_out.iloc[tj i + 2] = Q_out.iloc[tj i + 1] # + tributary flow
HP_depth.iloc[tj i+2] = (model_segs['init_depth'][i]) ## This is the
difference between Flow Routing and Stream Routing
if model_segs['ChanGeom'][i] == 'tz':
HP_aneaC.iloc[tj i+2] = model_segs['bot_width'][i] *
HP_depth.iloc[tj i+2] + model_segs['z_slope'][i] * HP_depth.iloc[tj i+2] ** 2
elif model_segs['ChanGeom'][i] == 'r':
HP_aneaC.iloc[tj i+2] = model_segs['bot_width'][i] *
HP_depth.iloc[tj i+2]
elif model_segs['ChanGeom'][i] == 't':
HP_aneaC.iloc[tj i+2] = model_segs['z_slope'][i] *
HP_depth.iloc[tj i+2] ** 2
HP_vel.iloc[tj i+2] = Q_out.iloc[tj i + 2] / HP_aneaC.iloc[tj i+2]
HP_vol.iloc[tj i+2] = model_segs['length'][i] * HP_aneaC.iloc[tj i+2]
HP_nt.iloc[tj i+2] = HP_vol.iloc[tj i+2] / Q_out.iloc[t, i + 2] / 24
/ 60 / 60
Changing Volume:
##loop through segments to calculate outflow and parameters for each time
step.
for t in nange((len(time_s))):
for i in nange(nsegments):
Q_out.iloc[tj i + 2] = Q_out.iloc[tj i + 1] # + tributary flow
HP_depth.iloc[tj i+2] =
l0**(model_segs['depth_intencept'][i])*Q_out.iloc[tj i +
2]**model_segs['depth_exp'][i] ## This is the only difference between Flow
Routing and Stream Routing
22
-------
if model_segs['ChanGeom'][i] == 'tz':
HP_aneaC.iloc[tj i+2] = model_segs['bot_width'][i] *
HP_depth.iloc[tj i+2] + model_segs['z_slope'][i] * HP_depth.iloc[tj i+2] ** 2
elif model_segs['ChanGeom'][i] == 'r':
HP_aneaC.iloc[tj i+2] = model_segs['bot_width'][i] *
HP_depth.iloc[tj i+2]
elif model_segs['ChanGeom'][i] == 't':
HP_aneaC.iloc[tj i+2] = model_segs['z_slope'][i] *
HP_depth.iloc[tj i+2] ** 2
HP_vel.iloc[tj i+2] = Q_out.iloc[tj i + 2] / HP_aneaC.iloc[tj i+2]
HP_vol.iloc[tj i+2] = model_segs['length'][i] * HP_aneaC.iloc[tj i+2]
HP_nt.iloc[tj i+2] = HP_vol.iloc[tj i+2] / Q_out.iloc[t, i + 2] / 24
/ 60 / 60
Kinematic Wave:
def penimeten_calc(shapej 6, wide, z_slope):
if shape =='tz':#fon trapezoid
penim=wide +2 *d*(l + z_slope**2)**0.5
elif shape =='r':#for rectangle
penim=2*d + wide
elif shape =='t':#fon triangle
penim=2 *d*(l + z_slope**2)**0.5
return penim
def depth_calc (shape, area, wide, z):
if shape =='r':#for rectangle
depth=anea/wide
elif shape =='tz':#fon tnapezoidshape =='tz':#fon trapezoid
depth=(-wide+np.sqnt((wide**2)-4*z*anea))/2*wide
elif shape =='t': #for triangle
depth == np.sqnt(anea/z)
return depth
##Loop through segments to calculate outflow and parameters for each time
step.
for t in nange(len(time_s)):
23
-------
for i in nange(nsegments):
p_calc=perimeter_calc(model_segs['ChanGeom'][i],HP_depth.iloc[t,i+2],
model_segs['bot_width'][i], model_segs['z_slope'][i])
alpha_data.iloc[tJi]=((model_segs['mannings_n'][i]/model_segs['chan_slope'][i
]**0.5)*(p_calc**(2/3)))**beta
Q_hat = Qout.iloc[tj i+2]
Qout.iloc[t+lj i+2] = Qout.iloc[tj i+2] + ((Qout.iloc[t+lj i+1] -
Qout.iloc[tj i+2])*(delta_t))/(model_segs['length'][i]*alpha_data.iloc[t,
i]*beta*Q_hat**(beta-l))
While abs((Q_hat - Qout.iloc[t+lj i+2])/Q_hat ) > le-3: # implicit
soln for Qout.iloc[t+lj i+2]
Q_hat = Qout.iloc[t+lj i+2]
Qout.iloc[t+lj i+2] = Qout.iloc[tj i+2] + ((Qout.iloc[t+lj i+1] -
Qout.iloc[tj i+2])*(delta_t))/(model_segs['length'][i]*alpha_data.iloc[t,
i]*beta*Q_hat**(beta-l))
HP_anea.iloc[tji+2]=alpha_data.iloc[t,i]*(Qout.iloc[t+l , i+2])**beta
HP_depth.iloc[t+lji+2]=HP_anea.iloc[t,i+2]/model_segs['bot_width'][i]
#necalculating alpha using new depth
p_calc2=penimeten_calc(model_segs['ChanGeom'][i],HP_depth.iloc[t+lj i+2],
model_segs['bot_width'][i], model_segs['z_slope'][i])
alpha_stan=((model_segs['mannings_n'][i]/model_segs['chan_slope'][i]**0.5)*(p
_calc2**(2/3)))**beta
while abs((alpha_stan-alpha_data.iloc[tji])/alpha_data.iloc[tji]) >
le-3:
alpha_data.iloc[tj i]=alpha_stan
HP_anea.iloc[tji+2]=alpha_data.iloc[tji]*(Qout.iloc[t+l ,
i+2])**beta
HP_depth.iloc[t+lJi+2]=depth_calc(model_segs['ChanGeom'] [i],
HP_anea.iloc[tJi+2]Jmodel_segs['bot_width'][i], model_segs['z_slope'][i])
alpha_stan=((model_segs['mannings_n'][i]/model_segs['chan_slope']
[i]**0.5)*(p_calc2**(2/3)))**beta
HP_velocity.iloc[tji+2]=Qout.iloc[t+l , i+2]/HP_area.iloc[t,i+2]
HP_vol.iloc[tj i+2]=model_segs['length'][i]*HP_anea.iloc[t, i+2]
24
-------
oEPA
United States
Environmental Protection
Agency
Office of Research
and Development
(8101R)
Washington, DC
20460
Official Business
Penalty for Private Use
$300
PRESORTED
STANDARD POSTAGE
& FEES PAID EPA
PERMIT NO. G-35
25
------- |