vvEPA
EPA/600/R-09/100 | September 2009 | vwvw.epa.gov/athens
United States
Environmental Protection
Agency
WASP7 Stream Transport
Model Theory and User's
Guide
Supplement to Water Quality Analysis
Simulation Program (WASP) User
Documentation
Ecosystems Research Division, Athens, GA 30605
National Exposure Research Laboratory
Office of Research and Development
-------
EPA/600/R-09/100
September 2009
WASP7 Stream Transport
Model Theory and Userfs
Guide
Supplement to Water Quality Analysis
Simulation Program (WASP) User
Documentation
by
Robert B. Ambrose, Jr., P.E.
U.S. EPA, Office of Research and Development
National Exposure Research Laboratory
Ecosystems Research Division
Athens, Georgia
Tim A. Wool
U.S. EPA, Region 4
Water Management Division
Atlanta, Georgia
U.S. Environmental Protection Agency
Office of Research and Development
National Expsoure Research Laboratory
Ecosystems Research Division
Athens, GA 30605
-------
NOTICE
The U.S. Environmental Protection Agency (EPA) through its Office of Research and
Development (ORD) funded and managed the research described herein. It has been
subjected to the Agency's peer and administrative review and has been approved for
publication as an EPA document. Mention of trade names or commercial products does
not constitute endorsement or recommendation for use.
11
-------
Abstract
The standard WASP7 stream transport model calculates water flow through a branching
stream network that may include both free-flowing and ponded segments. This
supplemental user manual documents the hydraulic algorithms, including the transport
and hydrogeometry equations, the model input and output, and a series of model
verification tests.
For one-dimensional, branching streams or rivers, flow routing can be calculated for free-
flowing stream reaches, for ponded reaches, and for backwater or tidally-influenced
reaches. Kinematic wave flow routing is a simple but realistic option to drive advective
transport through free-flowing segments. The kinematic wave equation calculates flow
wave propagation and resulting variations in flows, volumes, depths, and velocities
resulting from variable upstream inflow. This well-known equation is a solution of the
one-dimensional continuity equation and a simplified form of the momentum equation
that considers the effects of gravity and friction. Advective transport through ponded
segments is controlled by a sharp-crested weir equation. This equation calculates outflow
based on water elevation above the weir. Ponded reaches can include a flowing surface
water segment (epilimnion) as well as stagnant underlying segments (hypolimnion).
For dynamic flow through backwater segments, the momentum equations from
DYNHYD provide a simple solution for calculating outflows and resultant changes in
velocity, surface elevation, depth and volume. Driven by variable upstream flows and
downstream heads, the dynamic flow routine solves one-dimensional equations
describing the propagation of a long wave through a shallow water system while
conserving both momentum and volume. This approach considers the effects of gravity,
friction, and convective inertia, assuming that flow is predominantly one-dimensional,
that accelerations normal to the direction of flow are negligible, that channels can be
adequately represented by a constant top width with a variable hydraulic depth (i.e.,
"rectangular"), and that bottom slopes are moderate. This option can be used to represent
simple two-dimensional (x-y) water bodies with a branching link-node network.
To run the WASP7 stream transport module, the user must supply segment information
and flow information. Required segment information includes lengths, widths, and
depths for average flow conditions, as well as bottom slopes and Manning friction
coefficients. The hydrogeometric depth exponents may also be specified to control the
channel shape. Minimum channel depths for zero-flow conditions may be specified, with
a default value of 0.001 m. Bottom slopes less than 10"6 signify ponded or backwater
segments. For ponded segments, the minimum channel depth is interpreted as the outlet
weir height. If weir height is not specified, then the segment is treated as a backwater
using the dynamic flow option. Required flow information includes the flow pathways
for the main channel and each of the tributaries, as well as the inflow time functions for
the simulation period.
in
-------
Transport variables are provided as output for each of the WASP7 kinetic modules.
Standard output includes segment outflows, depths, velocities, and widths.
IV
-------
Table of Contents
Abstract iii
List of Figures vii
List of Tables vii
Introduction 1
2Background 1
2.1 WASP Transport Fields 1
2.2 WASP Surface Water Flow 2
2.2.1 WASP Surface Water Descriptive Flow Options 2
2.2.2 WASP Dynamic Stream Flow Option 3
3Development of Equations 4
3.1 Hydrogeometry 4
3.2 Governing Flow Equations 7
3.2.1 Kinematic Wave Flow 8
3.2.2 Ponded Weir Flow 9
3.2.3 Dynamic Flow 9
3.3 Implementation of Equations 10
3.3.1 Kinematic Wave Flow 11
3.3.2 Ponded Weir Flow 12
3.3.3 Dynamic Flow 14
4Stream Transport Model Inputs 16
4.1 Data set screen 16
4.1.1 Restart Options 16
4.1.2 Date and Times 17
4.1.3 Non-Point Source File 18
4.1.4 Hydrodynamics 18
4.1.5 Solution Technique 20
4.1.6 Time Step Definition 20
4.2 Segments screen 21
4.2.1 Description 21
4.2.2 Volume 21
4.2.3 Velocity (multiplier and exponent) 22
4.2.4 Depth (multiplier and exponent) 22
4.2.5 Segment Type 23
4.2.6 Bottom Segment 23
4.2.7 Length 23
4.2.8 Width 23
4.2.9 Minimum Depth 23
4.2.10 Slope 23
4.2.11 Segment Roughness 24
4.3 Flows screen 24
4.3.1 Flow Field 25
4.3.2 Flow Function 26
v
-------
4.3.3 Flow Path Function 27
4.3.4 Flow Time Function 28
SStream Transport Model Outputs 28
6References 29
TAppendix 1: Derivation of Equations 31
7.1 Hydraulic Exponents for Kinematic Wave Flow 31
SAppendix 2: Model Verification Tests 31
8.1 Kinematic Wave Tests 32
8.1.1 Stream Transport Test 1 32
8.1.2 Stream Transport Test 2 32
8.1.3 Stream Transport Test 3 32
8.1.4 Stream Transport Test 4 32
8.1.5 Stream Transport Test 5 33
8.2 Weir Overflow Verification Tests 33
8.2.1 Weir Overflow Test 1 - Steady Flow 33
8.2.2 Weir Overflow Test 2 - Variable Flow 33
8.2.3 Weir Overflow Test 3 -Long Term Dynamics 33
8.3 Dynamic Flow Verification Tests 33
8.3.1 Dynamic Flow Test 1 - Steady Backwater 33
8.3.2 Dynamic Flow Test 2 - Steady Flow, Elevation 34
8.3.3 Dynamic Flow Test 3 - Variable Stream Flow 34
8.3.4 Dynamic Flow Test 4 - EFDC Stream 34
8.3.5 Dynamic Flow Test 5 - Tidal Stream 34
8.3.6 Dynamic Flow Test 6 - 2-D Tidal Stream 34
VI
-------
List of Figures
Figure 1 - Model network with advective transport pathways 1
Figure 2 - WASP Descriptive Flow Options 3
Figure 3 - Channel hydraulic cross-sections 6
Figure 4 - Definition sketch for ponded flow 9
Figure 5 - WASP Main Screen Toolbar, data input buttons 16
Figure 6 - Dataset Parameterization Screen (WASP7.3) 17
Figure 7 - WASP Segment Definition Screen 21
Figure 8 - Flows screen 24
Figure 9 - Example WASP flow input 25
List of Tables
Table 1 - Hydraulic Exponents for Figure 3 6
Table 2 - Comparison of Empirical Hydraulic Exponents 6
Table 3 - Transport output variables 29
vn
-------
1 Introduction
The Water Quality Analysis Simulation Program, WASP7 (Di Toro et al. 1983,
Ambrose et al. 1988, Wool et al. 2001, Wool et al. 2006) is a general dynamic mass
balance framework for modeling contaminant fate and transport in surface waters. Based
on the flexible compartment modeling approach, WASP can be applied in one, two, or
three dimensions with advective and dispersive transport between discrete physical
compartments, or "segments." WASP provides a selection of modules to allow the
simulation of conventional water quality variables as well as toxicants.
The WASP kinetic models are based on a set of transport and transformation
equations. Advective transport is driven by water flow through a specified computational
network (e.g. Figure 1). Inflows bring boundary concentrations into the network, and
internal flows advect most constituents along specified flow paths through the network
and out the downstream boundaries.
Model Network
Figure 1 - Model network with advective transport pathways
2 Background
2.1 WASP Transport Fields
Advective transport in WASP is divided into six distinct types, or "fields." The
first transport field is advective flow in the water column. Advective flow carries water
quality constituents "downstream" with the water and accounts for instream dilution. The
second transport field specifies the movement of pore water in the sediment bed.
Dissolved water quality constituents are carried through the bed by pore water flow. The
-------
third, fourth, and fifth transport fields specify the transport of particulate pollutants by the
settling, resuspension, and burial of solids. Water quality constituents sorbed onto solid
particles are transported between the water column and the sediment bed. The sixth
transport field represents evaporation or precipitation from or to surface water segments.
2.2 WASP Surface Water Flow
Advective water column flows directly control the transport of dissolved and
particulate pollutants in many water bodies. In addition, changes in velocity and depth
resulting from variable flows can affect such kinetic processes as reaeration,
volatilization, and photolysis. In WASP, water column flow is input via transport field 1.
Circulation patterns may be described (flow options 1 and 2) or simulated by a
hydrodynamic model, such as DYNHYD, EPDRIV1, or EFDC (flow option 3).
Beginning with version 6.2, WASP is capable of internally calculating flow through a
one-dimensional, branching network using the kinematic wave formulation (flow option
4). In version 7.2, flow option 4 added ponded segments with weir overflow. In version
7.3, flow option 4 added backwater segments with dynamic, bi-directional flows.
2.2.1 WASP Surface Water Descriptive Flow Options
Two descriptive flow options are available in WASP - Net Flows (flow option 1)
and Gross Flows (flow option 2). For these descriptive flow options, WASP tracks each
separate inflow specified by the user from its point of origin through the model network.
For each inflow, the user must supply both a continuity (i.e., unit flow response) function
and a time function. The time function describes the inflow as it varies in time. The
continuity function describes the unit flow response throughout the network. The actual
flow between segments that results from each inflow is the product of the time function
and the continuity function.
If several inflow functions are specified, then the total flow between segments is
the sum of the individual flow functions. If the total flow into a segment is not equal to
the total flow out of a segment, then that segment volume is adjusted by WASP to
maintain continuity. Changes in inflows are instantly propagated throughout the model
network. While the steady-state effect of several tributaries can be described with these
flow options, unsteady flow transients are not calculated.
These descriptive flow options are illustrated in Figure 2. In the Net Flow option,
WASP sums all the flows at a segment interface to determine the direction of net flow,
and then moves mass in that ONE direction. In the Gross Flow option, WASP moves
mass with each separate flow at a segment interface. If opposite flows are specified at an
interface, WASP will move mass in BOTH directions. This option allows the user to
describe large dispersive circulation patterns. Note that if all flows are directed
downstream, then the mass transport in the Gross Flow option is equivalent to the mass
transport in the Net Flow option.
-------
I.Net I
2. Gros
7 low Option
10 cms 5 cms
5 c
ms
>
s Flow Option
10 cms 5 cms
10 c
5 c
<
ms
ms
Figure 2 - WASP Descriptive Flow Options
2.2.2 WASP Dynamic Stream Flow Option
The dynamic stream flow option was implemented to provide a more realistic
simulation of flow dynamics in branching, one-dimensional networks. This option can
include segments governed by kinematic wave flow, ponded weir overflow, and
backwater flow. WASP simulates flow through the network in response to time-variable
inflows and withdrawals.
As in the descriptive flow options, the user must supply both a continuity function
and a time function for each inflow (or withdrawal). Flow paths may diverge (branch)
and then re-join. For surface water segments, the user must specify bottom slopes and
roughness factors, as well as widths and depths for average flow conditions and a depth
hydraulic exponent for non-rectangular channels. The model uses the inflow and flow
path functions, along with specified channel geometry and hydraulic coefficients to
calculate time-variable water movement (flows and velocities) and channel
hydrogeometry (top widths, cross-sectional average depths, and volumes).
Kinematic wave routing was the first dynamic stream flow option introduced to
WASP6. Kinematic flow is controlled by bottom slope and bottom roughness. The
kinematic wave formulation can be used for most stream and small river reaches. For
these segments, the user must specify bottom slope greater than 10"6 and nonzero
roughness coefficients.
Beginning with version 7.2, the stream network can include ponded flow
segments along with kinematic flow segments. Flow through ponded segments is
controlled by a downstream low-head dam, weir, or natural sill. For these segments, the
user must specify bottom slope of 0 (or less than 10"6) and the downstream weir height.
Ponded segments may include stagnant underlying water layers.
-------
Beginning with version 7.3, the stream network can include dynamic flow (or
"backwater") segments along with kinematic flow and ponded flow segments. Dynamic
flow is controlled by gradients in surface elevation and velocity, as well as bottom
roughness. For these segments, the user must set bottom slope to 0 and specify bottom
elevation in reference to a downstream control point.
3 Development of Equations
3.1 Hydrogeometry
A good description of segment hydrogeometry as a function of flow can be
important in properly using WASP to simulate streams and rivers. For the hydrodynamic
linkage flow option, velocities and depths computed by the hydrodynamic model are used
in WASP. For the internal flow options (Net Flow, Gross Flow, Dynamic Stream Flow),
a set of user-specified hydraulic discharge coefficients defines the relationship between
velocity, depth, and stream flow in surface water segments. This method follows the
implementation in QUAL2E (Brown and Barnwell, 1987). For the descriptive flow
options (Net Flow, Gross Flow), segment velocities and depths do not influence the
transport scheme; they are only used in calculations of reaeration and volatilization rates.
For the Dynamic Stream Flow option, segment velocities, widths, and depths are integral
to the transport calculations.
Discharge coefficients giving depth and velocity from stream flow are based on
empirical observations of the stream flow relationship with velocity and depth (Leopold
and Maddox, 1953). The equations relate velocity, channel width, and depth to stream
flow through power functions:
Equation 1
v = vmult Qvexp
Equation 2
R = dmult- Q**
Equation 3
B = bmult-Q^
where v is velocity [m/sec], R is hydraulic radius, or cross-sectional average depth [m], B
is top width [m], vmult, dmult, and bmult are empirical coefficients, and vexp, dxp, and
bexp are empirical exponents. Cross-sectional area, A is the product of top width and
average depth, and from continuity, flow is given by:
Equation 4
Q = v-A = v-R-B = (vmult Qvexp ) (dmult Q**)- (bmult Qbexp)
= (vmult dmult bmult) Q^+^+b^
From inspection, the following hydraulic relationships hold:
-------
Equation 5
vmult dmult bmult = 1
Equation 6
v exp + dxp + b exp = 1
The Net Flow and Gross Flow options in WASP require specification of the hydraulic
relationships for velocity and depth; the width coefficients are calculated internally from
Equation 5 and Equation 6. The Kinematic Wave Flow option requires specification of
the hydraulic depth exponent dxp, along with depth Dm and width Bm under average flow
conditions. Manning's equation (rearranged) is used to calculate velocity vm under
average flow conditions, and then average flow Qm from depth, width, and velocity:
Equation 7
v =
) . V1/2
}m *f
n
Equation 8
A consistent set of hydraulic exponents are set (see Appendix, Section 7.1):
Equation 9
v exp = (%J dxp
Equation 10
b exp = 1 - dxp - v exp
Finally, a consistent set of hydraulic multipliers are then derived from mean flow width,
the hydraulic geometry equations, and Manning's equation:
Equation 11
bmult = Bm.Qm-"°*
Equation 12
vmult = 11 a
Equation 13
dmult = a I bmult
Equation 14
a =
n bmult
-------
Channel Profile
-505
Width
Channel
1. Rectangular
2. U-Shape
3. V-Shape
4. Shallow
Velocity
0.40
0.32
0.26
0.20
Depth
0.60
0.48
0.39
0.30
Width
0.00
0.20
0.35
0.50
Figure 3 - Channel hydraulic cross-sections
Table 1 - Hydraulic Exponents for Figure 3
Channel cross-sections for representative hydraulic geometry coefficients in Table
lare illustrated in Figure 3. Under mean flow, these channels are 10m wide and 0.5 m
deep. Leopold et al. (1964) have noted that stream channels in humid regions tend
towards a rectangular cross-section because cohesive soils promote steep side slopes
whereas noncohesive soils encourage shallow sloped, almost undefined banks.
Table 2 - Comparison of Empirical Hydraulic Exponents
Channel
Rectangular
Average of 158 U.S. Gaging Stations
Average of 10 Gaging Stations on Rhine River
Ephemeral Streams in Semiarid U.S.
Velocity
0.40
0.43
0.43
0.34
Depth
0.60
0.45
0.41
0.36
Width
0.00
0.12
0.13
0.29
Table 2 compares hydraulic exponents for a rectangular channel with data
reported by Leopold et al. (1964). Note that the average velocity exponent is relatively
constant for all channel cross sections. The major variation occurs as a decrease in the
depth exponent and concomitant increase in the width exponent as channel cross-sections
change from the steep side slopes characteristic of cohesive soils to the shallow slopes of
arid regions with noncohesive soils.
For site-specific river or stream simulations, hydraulic coefficients and exponents
must be estimated. Brown and Barnwell (1987) recommended estimating the exponents
(b and d) and then calibrating the coefficients (a and c) to observed velocity and depth.
The exponents may be chosen based on observations of channel shape noted in a
reconnaissance survey. If cross sections are largely rectangular with vertical banks, the
first set of exponents shown should be useful. If channels have steep banks typical of
areas with cohesive soils, then the second set of exponents is appropriate. If the stream is
-------
in an arid region with typically noncohesive soils and shallow sloping banks, then the last
set of exponents is recommended.
The key property of the channel that should be noted in a reconnaissance survey is
the condition of the bank slopes or the extent to which width would increase with
increasing stream flow. Clearly the bank slopes and material in contact with the stream
flow at the flow rate(s) of interest are the main characteristics to note in a reconnaissance.
This gives general guidance but it should be noted that values are derived for bankful
flows. Even in streams with vertical banks, the low flows may be in contact with a sand
bed having shallow sloped, almost nonexistent banks that are more representative of
ephemeral streams in semi-arid areas.
3.2 Governing Flow Equations
The WASP stream flow model consists of a set of one-dimensional equations
solving water flow and water volume in a branching stream or shallow river network.
This network can include free-flowing stream reaches (kinematic wave flow), ponded
reaches (weir overflow), and backwater or tidally influenced reaches (dynamic flow). The
equation of motion, based on the conservation of momentum, predicts water velocities
and flows. The equation of continuity, based on the conservation of volume, predicts
water heights (heads) and volumes.
The one-dimensional continuity equation is given by:
Equation 15
dx dt
where Q is volumetric flow, [m3/sec] and^4 is cross-sectional area [m2]. For rectangular
channels, where width is constant, Equation 15 becomes:
Equation 16
dO dH
+ B - = 0
dx dt
where B is channel width [m] and H is water surface elevation [m]. As presently
implemented in WASP, kinematic flow reaches have shapes described by the
hydrogeometric relationships described in Section 3.1, while ponded reaches and
dynamic flow reaches have rectangular channel shapes.
The equations of motion implemented in the three reach types are described in the
following sections.
-------
3.2.1 Kinematic Wave Flow
For one-dimensional, free-flowing stream reaches, kinematic wave flow routing is
a simple but realistic option to drive advective transport. The kinematic wave equation
calculates flow wave propagation and resulting variations in flows, velocities, widths, and
depths throughout a stream network. This well-known equation is a solution of the one-
dimensional continuity equation and a simplified form of the momentum equation that
considers the effects of gravity and friction:
Equation 17
g(s0-sf) = o
where g is acceleration of gravity [m/sec2], S0 is the bottom slope, and Sf is the friction
slope. Manning's equation expresses the friction force as a function of water velocity
and hydraulic radius:
Equation 18
where n is the Manning friction factor, v is water velocity [m/sec], and R is hydraulic
radius [m], which is equivalent to the cross-sectional average depth, D. From the
simplified momentum equation, So can be equated to Sf. Hydraulic radius can be
expressed as cross-sectional area divided by width, B [m]. Substituting these into the
Manning's equation and rearranging terms gives flow as a function of bottom slope,
cross-sectional area, and width:
Equation 19
1 A5'3
= -s12
n B23
Substituting this expression into the continuity equation and differentiating^ with respect
to time gives the kinematic wave differential equation:
Equation 20
= o
dx dt
where, for rectangular channels:
Equation 21
0 = 3/5, a =
1 V
For channels where width varies with flow, a and P are functions of hydraulic
coefficients:
-------
Equation 22
P = 0.6 + OA-dxp,
a =
n bmult
where the hydraulic coefficients dxp and bmult are defined in Section 3.1.
3.2.2 Ponded Weir Flow
For flow through ponded segments controlled by a downstream low-head dam or natural
sill (Figure 4), the sharp-crested weir overflow equation is a simple solution for
calculating outflows and resultant changes in depth and volume. Weir height Hw [m] and
width Bw [m] are specified by the user, and hydraulic head Hh is the difference between
ponded depth H and Hw.
H,
H
H
w
\
Q,
Figure 4 - Definition sketch for ponded flow
For a sharp-crested weir where Hh/Hw < 0.4, velocity and flow are related to head by
(Finnemore and Franzini 2002):
Equation 23
= 1-83 H[h'
Qo = V,
= 0
A =
1.83 flw
when H
-------
The dynamic flow routine solves one-dimensional equations describing the
propagation of a long wave through a shallow water system while conserving both
momentum and volume. The equation of motion calculates water velocities and flows.
The equation of continuity calculates surface elevations, along with associated depths and
volumes. This approach assumes that flow is predominantly one-dimensional, that
accelerations normal to the direction of flow are negligible, that channels can be
adequately represented by a constant top width with a variable hydraulic depth (i.e.,
"rectangular"), that the wave length is significantly greater than the depth, and that
bottom slopes are moderate. Reaches with steep bottom slopes are best solved with the
kinematic wave equations.
Equation 16 gives the equation of continuity implemented for dynamic flow
reaches. The equation of motion calculates local acceleration, the velocity rate of change
with respect to time [m/sec2]:
Equation 24
dv dv
= -v + a + af
dt dx g J
where v is water velocity [m/sec], ag is gravitational acceleration along the axis of the
channel [m/sec2], and a/is frictional acceleration [m/sec2]. The first term, convective
inertia or Bernoulli acceleration, represents the rate of momentum change by mass
transfer, [m/sec2]. The second term, gravitational acceleration, is driven by the slope of
the water surface:
Equation 25
dY_
where 7is surface elevation [m], and g is the acceleration of gravity [9.81 m/sec2]. The
third term, frictional acceleration, can be expressed using the Manning equation for
steady uniform flow:
Equation 26
af
R'
where n is the Manning friction factor, v is water velocity [m/sec], and R is hydraulic
radius [m], which is equivalent to the cross-sectional average depth, D.
Local wind acceleration is not included in this implemented of the dynamic flow
equations.
3.3 Implementation of Equations
WASP7 solves the kinematic flow, ponded flow, and dynamic flow equations for
appropriate surface water segments in a stream network using finite-difference
formulations for flow and for continuity.
For each segment, a maximum numerical time step DTmax is calculated from the
segment length and characteristic velocity, as described in the sections below. The overall
10
-------
time step is the product of the minimum DTmax in the network and a user-specified
fraction, DTP (default = 0.9) that is set to ensure stability:
Equation 27
DT = DTF.mm(DTmJ
This time step DTis divided into two half time steps. For kinematic wave reaches
and ponded weir reaches, flows are calculated sequentially for each half time step and
then averaged for subsequent use by the water quality module. Final velocities, depths,
volumes, and surface elevations at the end of the full time step are passed along to the
water quality module.
For dynamic flow reaches, the two half time steps are used in a predictor-
corrector scheme as described in Section 3.3.3 below. Flows, velocities, volumes, depths,
and surface elevations are updated following each of the numerical steps. Final flows,
velocities, volumes, depths, and surface elevations at the end of the full time step are
passed along to the water quality module.
3.3.1 Kinematic Wave Flow
To solve for flow in kinematic wave reaches, Equation 20 is expressed in finite difference
form:
Equation 28
Qot~ QDO _ Qot
T-.rp r> T U DO
DT a/3 L
where Qu is the upstream inflow, QDO is the outflow from the preceding time step, Qot is
the outflow for this time step, and DTis the time step [days]. Equation 28 is solved using
a Newton-Raphson approach:
Equation 29
DT
f(QDt) = Qnt+c.Q1^ +C2 = 0, c,= -- (QV-QD,\ C2=-QDO
a /3 L
Equation 30
f'(QDt) = l + Cl(l-/?)<2Df
Given an initial estimate of Qot, an updated estimate, Qot2, is calculated by:
Equation 31
/ (QDt)
f'(QDt}'
Qna ~ Q
Dt
Or
11
-------
Equation 29, Equation 30, and Equation 31 are solved in an iterative loop where Qnt it set
equal to Qot2 until err is less than 10"5. Given the new set of flows for the water column
network Qm,ij, volumes for all water segments "/" are updated using the continuity
equation:
Equation 32
DV - V(9 DT V - V + DV
^^i ~ Z-i^Dt,ij ^1 ' y t,i ~ ^ Qj ~t- -LS r j
j
Segment widths are updated with the new flows using Equation 3. Associated cross-
sectional areas and depths are then calculated:
Equation 33
Equation 34
A
To prevent slow numerical drift in calculated volumes during lengthy simulations, small
adjustments are made to flows based on differences between hydraulic radius calculated
from Equation 2 and cross-sectional average depth calculated through continuity.
Applying Equation 19:
Equation 35
n, °''
Equation 36
Qotj = QotJ + Qerr,
For each kinematic flow segment, a maximum stable numerical time step DTmax [days] is
calculated from the segment length L [m] and celerity c [m/s]:
Equation 37
c = vl p
Equation 38
nr = °-5 (L/\
86400 ' vc*
where 0.5 is a safety factor.
3.3.2 Ponded Weir Flow
For ponded reaches "/'", weir overflow (Equation 23) must be solved along with
continuity such that:
12
-------
Equation 39
and
Equation 40
rr TT
Hht,, = Hho,,
L
where Qut.t and Qnt.i are the upstream inflow and outflow for this time step [mVsec], H/
is the head for this time step [m], H/^- is the head from the previous time step [m], and
DT is the time step [days]. These equations are solved using a Newton-Raphson
approach where:
Equation 41
_J = Got ~ 1.83 flfl^ = 0
Equation 42
f(QDt] = 1 - 1.835(|#£ DT^
,2 nv B,L
z~z-
1.5-1.83-ZXT
Given an initial estimate ofQDt,j, a consistent value for Hhtis calculated using Equation
40. An updated estimate Qot2, is then calculated by:
Equation 43
- n _
~ *
f (QDt)
f'(QDt)'
Qna ~ Q
Dt
Or
Equation 41, Equation 42, Equation 43, and Equation 40 are solved in an iterative loop
where Qot and Hht are set equal to Qot2 and Hht2 until err is less than 10"5.
WASP7 solves this weir overflow equation for each ponded segment in a stream
network. In each of the numerical steps, calculated values ofdQ/dt and Q0 are used to
update volumes, depths, and Hh, which are used in the next numerical step to calculate
For each weir overflow segment, a maximum stable numerical time step DTmax
[days] is calculated from the segment length and overflow velocity v0:
Equation 44
DT = °'5
86400
where 0.5 and 1.5 are safety factors.
13
-------
3.3.3 Dynamic Flow
Equation 16 and Equation 24 form the basis of the hydrodynamic model
DYNHYD5, which is implemented within this version of WASP. These equations are
integrated numerically on a flexible, computationally efficient "link-node" network
(Feigner and Harris, 1970), which solves the equations of motion and continuity at
alternating grid points. At each time step, the equation of motion is solved at the links (or
"channels"), giving velocities for mass transport calculations, and the equation of
continuity is solved at the nodes (or "junctions"), giving heads for pollutant concentration
calculations. Link-node networks can treat fairly complex branching flow patterns and
irregular shorelines with acceptable accuracy for many studies. They cannot handle
stratified water bodies, small streams, or rivers with a large bottom slope. Link-node
networks can be set up for wide, shallow water bodies if primary flow directions are well
defined.
In WASP, nodes correspond to segments, and links correspond to segment
interfaces. For every dynamic flow segment in a WASP network, a distinct channel
number "ich" is defined for each of its downstream segments. Channel ich is defined by
upstream segment y and downstream segment /'. Positive flow in channel ich is outflow
from segment y' and inflow to segment /'. Negative flow in channel ich is outflow from
segment /' to segment y'.
In finite difference form, Equation 24 is given by:
Equation 45
Vtich - V,.
ich
= v
~
V,ch
where vf;!C/, is the velocity for this time step [m/sec], vich is the velocity from the preceding
time step [m/sec], Ax!C/j is channel ich length [m], Av!C/, lkxich is the velocity gradient in
channel ich with respect to distance [sec"1], &HiCh I^Xtch is the water surface gradient in
channel ich with respect to distance [m/m], and DTis the time step [days]. All values on
the right hand side of equation 20 are referenced to the previous time step.
The water surface gradient can be computed from the junction heads at either end
of the channel:
Equation 46
Axlch (Lj+Lt)/2
where Hj and Lj are the water surface elevation and length of the upstream segment[m],
and Hi and Lt are the water surface elevation and length of the downstream segment[m].
The velocity gradient cannot be computed directly from upstream and downstream
channel velocities because of possible branching in the network. If branching does occur,
there would be several upstream and downstream channels, and any computed velocity
gradient would be ambiguous. An expression for the velocity gradient within a channel
can be derived by applying the continuity equation to the channel and substituting
forO:
14
-------
Equation 47
dA
~dt
dO dA
=- = -v
(J jC (J jC
dx
Rearranging terms gives the channel velocity gradient:
Equation 48
dv 1 dA v dA
dx A dt A dx
Writing this in finite difference form and substituting BxR for A and BxAHfor AA gives
the velocity gradient term:
Equation 49
Av,.
ich
1
ich
V,ck
ich
Ax,
ich
R
ich
Ax
ich
The term AHich/At is computed as the average water surface elevation change
between the segments at each end of channel ich during time step t. Substituting
Equation 49 into Equation 45 and rearranging gives the explicit finite difference equation
of motion applied to each channel:
Equation 50
Vt,,ch = V,ch
DT
v.
ich
ich
RKh ^
R
Vich A-T/id
"zc/j ) ^-Xich
iiiL -y
,4/3 ich
V,
ich\
Writing the equation of continuity in finite difference form and rearranging terms gives:
Equation 51
H.
BJ LJ
Equation 50 and Equation 51 are solved using a 2-step predictor-corrector routine. Based
on initial velocities, surface elevations, and depths from the previous time step, new
velocities and flows are solved for the half time step, along with new surface elevations
and depths. Using these predicted half-time step values, velocity and flow derivatives are
recalculated for the half time step. These corrected derivatives are then used with the
initial velocities, depths, and surface elevations to calculate velocities and flows for the
full time step. Finally, surface elevations, depths, and volumes are calculated for the full
time step.
For each dynamic flow segment, a maximum stable numerical time step DTmax [days] is
calculated from the segment length L [m] and celerity c [m/s]:
15
-------
86400 /c/
where 0.5 is a safety factor.
4 Stream Transport Model Inputs
To implement stream flow routing, the user must specify information in the Data
set screen, the Segments screen, and the Flows screen, accessed from the gears, the cube,
and the faucet on the main WASP toolbar (Figure 5). Each of these is briefly described in
the sections below.
= e ~'3,ect Pre-processor ?»'ode! Post-Processor He:p
Qi&KjsJJL
ess?
Figure 5 - WASP Main Screen Toolbar, data input buttons
4.1 Data set screen
Figure 6 shows the Data set screen for version 7.3. In this screen, the user must
select the Model Type. In the Time Range section, the user must specify the simulation
Start Date, Start Time, End Date and End Time. In the hydrodynamics section, the user
must select the flow option. Finally, the Time Step information should be specified.
Default values are supplied for Fraction of max time step (0.9), Max time step (1.0 day),
and Min time step (0.0001 day). These values should work well for most cases, but if
numerical instability is encountered, lowering the Fraction of max time step (to 0.5 or
even 0.1) could help. In some cases, the user may want to specify a maximum times step
of less than a day. If diurnal output is desired, then a maximum time step of 0.1 or 0.05
days should give the necessary precision.
When creating a new input dataset the input parameterization data entry form is
the first one that needs to be completed. This form provides basic information that is
needed by the program to parameterize the other data entry forms that follow. This
screen informs the program what type of WASP file you are going to be creating.
4.1.1 Restart Options
The methods used by WASP to read and create restart files have changed
substantially in this version. In previous versions the user would have selected Create
Restart File, for WASP to write the final conditions of the simulation to an output file.
16
-------
This is true for the current version as well. If the user wants to restart a simulation with
the final conditions of previous simulation this radio must be set. At the end of the
WASP simulation a restart file with the same name as the WIF except with the extension
*.RST will be saved. With the current release of WASP if the user wants to use a restart
file they simple click on the Load Restart File button, this will allow the user to browse to
whatever restart file they want to use. Once the file is selected and the user clicks on the
Okay button, the restart file is opened up and segment volumes and state variable initial
conditions are reset to the values in the user selected *.RST file.
ftarameterS'.
S'PWfK.7' ""'i i"~'.,i&,-">?- V~?
Description
J Floyd's Fork TMDL Model
Comments
1/31/2QQ4-- End time
I Time Range
Start Date
I 1/1/1999
I Start Time
I n-nn
End Date
j ci -ly-iqqq
End Time
n-nn
| Skip Ahead to Date
| 1/1/1999
Skip Ahead Time
0:00
Model Type
JEutrophication
'Jon Point Source File
P Use NFS file r^ :.
MDC Clip K]=mp *v~±"
Hydrodynamics
(~~ Net Flows
f Gross Flows
/"* 1 Pi M > Iv If' f'
f~* Hydrodynamic Linkage
Hydrodynamic Linkage File
&:....
Solution Technique
EULER jrj
f* Disable WASP to WASP linkage
r Enable WASP to WASP linkage
^^"^^^^^^^^^^^^^*- 1
--. Restart Option
-^1 ff No Restart File
i P Create Restart File
I (J/ Load restart file now
Bed Volumes
G Static
1 i Time Step
j | 0.90
j { 1.0000
I Min time step
{ | 0.0030
! Solution Options
[ f~ Negative Solution Allowed
| ! ^7 OK i| x Cancel
Figure 6 - Dataset Parameterization Screen (WASP7.3)
4.1.2 Date and Times
The previous versions of WASP did not require that the model time functions be
represented in Gregorian date format. WASP requires all time functions be represented
in Gregorian fashion (mm/dd/year hh:mm:ss).
17
-------
Start Date and Start Time - The start time dialog is used to define the date and time for
the start of the simulation or time period being considered in the model input files. This
date and time correspond to time zero within the model.
End Date and End Time - The end time dialog is used to define the date and time when
the simulation will end.
Skip Ahead to Date and Time - This new addition to the WASP interface allows the
user to skip to any portion of the simulation and/or the selected loaded hydrodynamic
linkage file. When the user selects a hydrodynamic linkage file the start time and end
time of the file is read and the interface automatically sets the beginning and end time to
these values. It is best that the user build all of the time series (environmental, boundary
and dispersion) to cover this full range of time. Once the WIF is built the user can set a
date and time to skip the simulation from the start time set in the hydrodynamic file. This
is handy for using the whole hydrodynamic linkage file for calibration and verification,
and then using a small portion of the hydrodynamic linkage file for scenario analysis. It
could be the critical time period that will be used for the waste load allocation or TMDL.
The start time of the simulation should still be set the beginning time in the
hydrodynamic linkage file. The user can change the end time of the simulation by
changing the last date/time pair in the Time Step screen.
4.1.3 Non-Point Source File
The non-point source file is an external file that contains a time-series of loads
(kg/day) for a given segment and system. This file is typically created either by the user
manually or using other software like the Stormwater Management Model (SWMM) in
conjunction with the Linked Watershed/Waterbody Model. This file can be used to
provide loading information to WASP on virtually any time scale, from time step to time
step, to year average loads.
4.1.4 Hydrodynamics
There are currently four surface flow options available for WASP:
1. Net Flows - This descriptive flow option is used for stream systems. WASP
uses the input main stem and tributary flow paths to calculate net transport
across all segment interfaces. If the user specifies opposing flow paths across
a segment interface, WASP will net the flows across the interface and move
water and constituent mass in the net direction of flow. This option is best
used under steady flow conditions, as changes in flow do not result in changes
in volume.
2. Gross Flows This descriptive flow option is used for simple lake or
estuarine systems. WASP uses the input flow paths to calculate transport
across all segment interfaces. Each flow separately participates in the mass
balance calculations. If the user specifies opposing flow paths across a
segment interface, WASP will move water and constituent mass in both
directions of flow. Again, changes in flow do not result in changes in volume.
18
-------
3. Dynamic Stream Flow For one-dimensional, branching streams or rivers,
dynamic flow is calculated internally using kinematic wave flow routing,
ponded weir overflow, or backwater flow equations. Networks can be
converging or diverging. These hydrodynamic equations calculate flow wave
propagation and resulting variations in flows, volumes, depths, width, and
velocities throughout a stream network. For WASP7.3 and earlier, this option
is titled "1-D Network Kinematic Wave."
4. Hydrodynamic Linkage Realistic simulations of unsteady transport in
rivers, reservoirs, and estuaries can be accomplished by linking WASP7 to a
compatible hydrodynamic simulation. This linkage is accomplished through
an external "hyd" file chosen by the user at simulation time. The
hydrodynamic file contains segment volumes at the beginning of each time
step, and average segment interfacial flows during each time step. WASP7
uses the interfacial flows to calculate mass transport, and the volumes to
calculate constituent concentrations. Segment depths and velocities may also
be contained in the hydrodynamic file for use in calculating reaeration and
volatilization rates. Before using hydrodynamic linkage files with WASP, a
compatible hydrodynamic model must be set up for the water body and run
successfully, creating a hydrodynamic linkage file with the extension of
*.hyd. This is an important step in the development of the WASP input file
because the hydrodynamic linkage file contains all necessary network and
flow information. When Hydrodynamic Linkage is selected in the Data Set
Parameters screen, the user cannot provide any additional surface flow
information. When you are ready to begin the development of a WASP input
deck, simply open the hydrodynamic linkage file from within the data
preprocessor. The hydrodynamic linkage dialog box allows the user to browse
and select a hydrodynamic linkage file. The data preprocessor will open the
hydrodynamic interface file and extract the number of segments, the starting
and ending time. The data processor will also determine the set of boundary
segments (segments that receive flow from outside the model network) and set
the boundary concentrations to 1.0 mg/L. Once a hydrodynamic linkage file
is selected in the data preprocessor, WASP has enough information to execute
a simple test run with no loads or kinetics enabled. This step is recommended
to test the network and transport integrity. If the simulation is run for a
sufficient duration, concentrations should approach 1.0 mg/L throughout the
network. If you are getting a number other than 1 mg/L, you may have to use
a different time step in the hydrodynamic model. This is especially true if the
concentrations are oscillating between large and small numbers, a clear
indication of numerical instability. WASP has the ability to get hydrodynamic
information from a host of hydrodynamic models. If a hydrodynamic model
does not support the WASP linkage it is relatively straightforward to create a
hydrodynamic linkage file. The hydrodynamic models that currently support
the WASP7.X file format are: EFDC (three dimensions), DYNHYD (one
dimension, branching), RIVMOD (one dimension, no branching, CE-QUAL-
RIV1 (one dimension, branching), SWMM/Transport (one dimension,
branching, SWMM/Extran (one dimension, branching)
19
-------
4.1.5 Solution Technique
The user now has the ability to select the model solution technique to be used by
the water quality module during the simulation. Currently there are 3 solution techniques
that can be selected: 1) Euler - which is the traditional solution technique that has been in
WASP since its inception, 2) COSMIC Flux Limiting - this solution technique is
typically used when WASP is linked to multi-dimensional hydrodynamic models like
EFDC, 3) Runge-Kutta 4 step solution technique used for diurnal simulations.
4.1.6 Time Step Definition
Starting with WASP Version 7.3 the user no longer has control over the
computational time step. Time step optimization routines have been refined to the point
where the model can determine what the most appropriate time step should be used next.
This assures the most efficient run time as well as minimizing numerical dispersion
caused by too small of a time step. While the user can not set the time step directly, they
do have some control over what would be an acceptable time step.
Fraction of Maximum Time Step - This dialog box specifies what fraction of the model
calculated time step will be used for the next time step. Its primary purpose is aid the
user in keeping the model stable. The default is 0.9 (or 90%) of the optimal time step.
Maximum Time Step - This specifies the maximum time step that will be used. If the
time step optimizer calculates a time step larger than this value, this value will be used.
This could be important in constraining the time step for diurnal or daily calculations.
Minimum Time Step - This specifies the minimum time step that will be used. The
default minimum time step is defined in the model as 0.0001 days. Use this dialog to
raise the minimum time step.
20
-------
4.2 Segments screen
In the Segments screen, the user must enter a row for each segment in the model
network (e.g., by pressing "Insert," by pressing the down arrow from the bottom row or
by copying from spreadsheet). For each segment, the user must enter a minimum amount
of information, which depends on the selected flow option. For the Net Flows and Gross
Flows options, the user must enter Depth (i.e. "Depth Multiplier"), Segment Type,
Bottom Segment, and either Volume or Length and Width. If an input segment volume is
0, then WASP will calculate an initial volume as the product of Length, Width, and
Depth. For the Dynamic Stream Flow option, the user must enter Depth, Segment Type,
Length, Width, Minimum Depth, Slope, and Bottom Roughness. Some default values are
provided, including Description ("Wasp Segment"), Segment Type ("Surface"), Bottom
Segment ("none"), Depth Exponent (0.6), Minimum Depth (0.001 m), and Bottom
Roughness (0.05).
4.2.1 Description
Descriptive segment names may be entered, but are not required. If the user wants
to paste segment descriptions/names that include spaces, the segment description must be
placed in quotes (i.e. "120 Bridge").
4.2.2 Volume
Segment volumes [m3] should be specified when using the Net Flow or Gross
Flow options. If a segment volume is not entered (or is 0), then WASP will calculate that
volume from specified segment length, width, and depth.
For the Dynamic Stream Flow option, input volumes are only used for benthic
segments. WASP calculates initial water column segment volumes from length, width,
and depth under initial flow conditions.
When using the hydrodynamic linkage flow option, initial water column volumes
are read from the external hydrodynamic file. Only benthic segment volumes must be
entered.
Segments Parameter initial Cun entration; ] F
| Segment Description Volume
1
' .. > asp Segment 0
'V Jf 'fyiat 0
, 4 «, asp Segment 0
5 Wajp Segment 0
6 Wasp Segment 0
7 Wasp Segment 0
8 Wasp Segment 0
9 Wasp Segment 0
10 Downstream 0
.
raction Dissolved [
Velocity
Multipliei
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
Velocity
Exponent
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
Depth
Multiplier
02890
0.2890
0.2890
0.2890
0.2890
0.2890
0.2890
0.2890
0.2890
0.2890
Volume Scale Factor
1.0000000
+ Insert [ Delete | | J OK ]
X Cancel |
Depth
Exponent
0.4500
0.4500
0.4500
0.4500
0.4500
0.4500
0.4500
0.4500
0.4500
0.4500
Segment
Type
Surface
Surface
Surface
Surface
Surface
Surface
Surface
Surface
Surface
Surface
Bottom
Segment
None
None
None
None
None
None
None
None
None
None
Length
1000.0000
1000.0000
1000.0000
1000.0000
1000.0000
1000.0000
1000.0000
1000.0000
1000.0000
1000.0000
Width
2.0000
20000
2.0000
2.0000
2.0000
2.0000
2.0000
2.0000
2.0000
2.0000
Minimum
Depth
0.5000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
Slope
0.0000
0.0010
0.0010
0.0010
0.0010
0.0010
0.0010
0.0010
0.0010
0.0010
^^^1
Bottom
Roughness
0.0400
0.0400
0.0400
0.0400
0.0400
0.0400
0.0400
0.0400
0.0400
0.0400
Volume Conversion Factor
| 1.0000000
Figure 7 - WASP Segment Definition Screen
21
-------
4.2.3 Velocity (multiplier and exponent)
The velocity hydraulic multipliers [(m/sec) / (mVsec)] and exponents should be
specified when using the Net Flow or Gross Flow options. For the Dynamic Stream Flow
option, velocity exponents should be left at 0. Initial velocities [m/sec] should be
specified in the velocity multiplier column for backwater segments only. For segments
governed by kinematic wave flow or ponded weir overflow, the velocity multipliers and
exponents are internally calculated from the input depth multipliers and exponents along
with the input width. Any input velocity multiplier or coefficient will be ignored when
using these options.
4.2.4 Depth (multiplier and exponent)
The depth hydraulic multipliers [m / (m3/sec)] and exponents should be specified
when using the Net Flow, Gross Flow, or Dynamic Stream Flow options. Depth
multipliers are required for all segments. For benthic segments, the depth multipliers are
interpreted as segment depths [m].
For the Net Flow and Gross Flow options, the depth multipliers and exponents are
used along with initial segment flows to calculate initial segment depths. If a depth
multiplier is left at 0, it is reset internally to 1.0 and a message is issued to the screen. If a
depth exponent is left at 0, then the depth multiplier is equal to the initial segment depth
[m]. During simulations using these descriptive flow options, changing flows do not
directly change segment depths, even if the hydraulic exponent is nonzero. Depths are
recalculated along with volumes based on flow continuity. If a segment outflow
continuity multiplier is equal to the inflow continuity multiplier, then changing flows will
not alter that segment's volume or depth.
For the kinematic wave segments in the Dynamic Stream Flow option, the depth
multiplier is taken to be the cross-sectional average segment depth under average flow
conditions [m]. The depth exponent is a value generally between 0.3 and 0.6. If a
segment depth exponent is left at 0, a rectangular cross-section is assumed and the
exponent is reset internally to 0.6. The average depths and depth exponents are used
along with segment widths, slopes, and roughness factors to calculate consistent
hydraulic coefficients, which are then used to calculate segment depths under initial flow
conditions. During simulations, changing flows directly change hydraulic depths based
on the hydraulic coefficients. Total segment depth is equal to the hydraulic depth plus
the user-input zero-flow minimum depth.
For the ponded segments in the Dynamic Stream Flow option, the depth
multiplier is taken to be the initial cross-sectional average depth. During simulations,
depths are recalculated along with volumes based on flow continuity. Total segment
depth is equal to the calculated depth plus the minimum depth, which is set internally to
0.001 m.
For backwater segments in the Dynamic Stream Flow option, the depth multiplier
is taken to be the initial cross-sectional average depth. The depth exponent should be left
at 0, indicating a rectangular cross-section. During simulations, depths are recalculated
along with volumes based on the momentum and flow continuity equations. Total
segment depth is equal to the calculated hydraulic depth plus the minimum depth.
22
-------
4.2.5 Segment Type
Segment type is entered using a pick list. Four segment types are available:
"Surface," "Subsurface," "Surface Benthic," and "Subsurface Benthic." The default
segment type is "Surface," which represents upper water column segments in contact
with the atmosphere. "Subsurface" represents underlying water column segments.
"Surface Benthic" represents the upper benthic sediment segments in contact with the
water column. "Subsurface Benthic" represents underlying benthic segments.
4.2.6 Bottom Segment
Bottom segment is the segment immediately underneath the current segment. The
bottom segment is entered by typing in the segment number or by using a pick list. If no
segments are underneath the current segment, then the bottom segment is designated
"none."
4.2.7 Length
Segment length [m] is the bottom length along the center of the flow line from the
upstream end to the downstream end of the segment.
4.2.8 Width
Segment width [m] is the top width averaged along the length of the segment. If
no input volume is specified, then width is used along with length and depth (multiplier)
to calculate an initial volume.
For the Dynamic Stream Flow option, top width should be specified for average
flow conditions. For kinematic wave segments, average top widths are used along with
average depths, depth exponents, slopes, and roughness coefficients to back-calculate a
consistent set of hydraulic coefficients. For ponded weir overflow segments, top width
values should be specified for the weir. For backwater flow segments, top widths are
constant in time, assuming rectangular channel cross-sections.
4.2.9 Minimum Depth
Minimum depth [m] is the average segment depth under zero-flow conditions,
used only in the Dynamic Stream Flow option. If this cell is left blank, a default
minimum depth of 0.001 m is assigned internally. Total depth is hydraulic depth plus
minimum depth. For ponded segments, the minimum depth field is used to enter the weir
height.
At present, the minimum depth field is also used to input bottom elevations for
the dynamic flow option. Bottom elevations are negative values representing the vertical
distance from the segment bottom (cross-sectional average) to the bottom of the
downstream control segment, which is either a weir or a boundary.
4.2.10 Slope
Segment slope [m/m] is the elevation drop divided by length averaged over the
segment length. This is usually calculated as the upstream elevation minus the
downstream elevation divided by segment length. Slope is used only in the Dynamic
Stream Flow option. Segments with slope greater than 10"6 are considered free-flowing,
23
-------
and flow is computed using the kinematic wave equation. Segments with slope less than
or equal to 10"6 are considered ponded, and flow is computed using the weir overflow
equation.
4.2.11 Segment Roughness
Segment roughness is the Manning's roughness coefficient n. Roughness is used
in the Dynamic Stream Flow option for kinematic wave flow and dynamic flow
segments. Roughness coefficients should usually be between 0.01 and 0.15. If a
coefficient of 0 is input for a free-flowing segment, WASP will reset the coefficient to
0.05 and issue a message to the screen.
4.3 Flows screen
The Flows screen is used to define advective transport, including surface water
and pore water flow, as well as solids settling and resuspension, precipitation and
evaporation. The Flows screen is also used to define downstream boundary elevations
and two-dimensional channel networks for the Dynamic Flow option.
The flow input screen is a complex screen that contains four tables (Figure 8).
The upper left quadrant is used to select the transport field, such as "Surface Water" flow.
For each transport field selected, the upper right quadrant is used to define a set of
transport functions, including upstream and tributary inflows. For each transport function,
the bottom two quadrants are used to define the flow path and the flow time function.
" *- 1 »
- 1 -
V--*. -t ,"! ~
l>l, ,,1
nTHa*
t i,rs
.,.
61 ""la
* * . (
r*,r t^
Figure 8 - Flows screen
24
-------
For surface water and pore water transport, the upstream inflow, each tributary inflow,
pore water inflow, and any flow withdrawals must be described by continuity path
functions and inflow time functions. An example is shown in Figure 9.
Figure 9 - Example WASP flow input
4.3.1 Flow Field
The transport field must be selected. Six transport fields are available:
1. Surface Water - This transport field is used to describe surface water flows. These
flows transport both the particulate and dissolved fractions of a constituent. If the
user has selected the hydrodynamic linkage option they will not be able to enter
information here.
2. Pore Water - This transport field is used to describe pore water flows. These flows
transport only the dissolved fraction of a constituent.
25
-------
3. Solids 1 - This transport field is used to describe solids type 1 settling and
resuspension. These flows transport only the particulate fraction of a constituent that
is mapped to solid type 1 in the Systems screen.
4. Solids 2 - This transport field is used to describe solids type 2 settling and
resuspension. These flows transport only the parti culate fraction of a constituent that
is mapped to solid type 2 in the Systems screen.
5. Solids 3 - This transport field is used to describe solids type 3 settling and
resuspension. These flows transport only the parti culate fraction of a constituent that
is mapped to solid type 3 in the Systems screen.
6. Evaporation/Precipitation - This transport field subtracts/adds water from the model
network. No constituent mass is added, removed, or transported.
Scale Factor - The scale factor for a transport field multiplies all flows associated with
that field by the input value. This is generally used to scale flows in sensitivity tests. The
default value is 1.0.
Conversion Factor - The conversion factor for a transport field multiplies all flows
associated with that field by the input value. This is generally used to adjust input flow
units to the internal units of mVsec. If flows are specified in ft3/sec, the conversion factor
should be 0.02832. The default value is 1.0.
4.3.2 Flow Function
The user can define several flow functions for the selected transport field. Each
flow function must have its own flow path function (lower left table) and flow time
function (lower right table). Normally, a Flow Function defines a discrete inflow, such as
upstream flow, tributary flow, or pore water flow. Special flow functions are also used in
conjunction with the Dynamic Flow option to define downstream boundary elevations or
two-dimensional (x-y) channel networks, as described in the sections below.
To insert a flow function, first highlight the Surface Water flow field in the upper
left table, then move the cursor to the upper right quadrant and click on the insert button.
The resulting flow function cell, labeled "Flow Function," can be edited to provide a
descriptive name.
To insert additional flow functions, either click on "insert" or highlight the last
flow function and press the down arrow. To delete a flow function, select the function by
highlighting the row and click on the delete button. Deleting a flow function will delete
the corresponding flow path function (lower left table) and flow time function (lower
right table).
Function Name - When a flow function is inserted, it is given the default name "Flow
Function." The Function cell can be edited to provide a descriptive name, such as
"upstream inflow," "tributary inflow," "downstream elevation," or "channel network."
Interpolation Option - The default interpolation option for the flow time function
associated with a flow function is "Linear." This can be changed to "Step" to provide for
26
-------
a step function. To change options, click in the Interpolation cell, press the down arrow,
and select the interpolation option for this flow function.
Once a flow function is selected and named, the user must define the associated
flow path function and flow time function. Be sure that the correct flow field and flow
function are highlighted before entering these next screens.
4.3.3 Flow Path Function
The flow path function traces this flow from its point of entry into the model
network to its point of exit either from the model network or to an alternate pathway
associated with another flow function. The flow path consists of a set of rows,
corresponding to segment interfaces. Each row will have a set of segment pairs and a
fraction of flow multiplier.
Segment Pairs - The segment pairs consist of a "From" segment and a "To" segment,
and define the direction of flow across this segment interface. Either the "From" or the
"To" segment can be defined as "Boundary." Normally the first row will define the
inflow from "Boundary" to the upstream segment and the last row will define the outflow
from the downstream segment to "Boundary." If this flow path is a tributary, then the last
row will define the outflow from the downstream tributary segment to a segment in
another tributary or the main stem of the river.
Positive values of flow transport water and constituent mass in the defined
direction from the first segment to the second segment. Negative flows transport water
and constituent mass from the second segment to the first segment. For example, if
"From" is segment 1 and "To" is segment 2, then negative values of flow in the time
function will cause transport from 2 to 1. Note: While the kinematic wave option checks
to make sure that all flow paths are ultimately connected to outflows, neither the
preprocessor nor the model can assure that the segments are connected properly.
Connectivity is the responsibility of the user.
Fraction of Flow - The fraction of flow column defines what fraction of the total flow in
this pathway moves between these segment pairs. For surface water flow, the fraction of
flow is normally 1.0. This allows the user to split flows from one segment into two or
more downstream directions. This can be used to define diverging and converging flows,
but must be used carefully. The sum of all fractions entering each segment must normally
equal the sum of all fractions leaving. If the sum is greater than 1.0, then that segment's
volume will continually increase. If the sum is less than 1.0, then that segment's volume
will continually decrease; if the volume reaches 0, the simulation will end badly.
Note for downstream boundary elevation - If a downstream boundary elevation
function is being defined for a dynamic flow network, then the flow path should consist
of a single segment pair from "Boundary" to the downstream segment number. The
fraction of flow should be set to 0.
Note for channel network - If a two-dimensional channel network is being defined for
the Dynamic Flow option, then each segment pair defines a unique flow channel. The
27
-------
fraction of flow multiplier is used to enter the initial channel cross-sectional area.
Channel hydraulic radius is calculated internally as the average of the upstream and the
downstream segment depths. Channel width is calculated internally as the cross-sectional
area divided by hydraulic radius. The channel network must not include any boundaries.
Channels connect two segments within the model network. WASP distinguishes channel
networks from traditional flow paths by the absence of boundaries and flow path
multipliers greater than 1.0. The flow time function associated with a channel network is
not used, and flow values can be left at 0.
4.3.4 Flow Time Function
The flow time function is a table consisting of dates, times, and inflow values
[mVsec]. Each row in the table represents a single point in time. During a simulation,
inflows are interpolated between these points based on the flow function interpolation
option selected (see Section 4.3.2). At least two rows must be entered in the flow time
function to allow for interpolation.
Date - As with other WASP time functions, the date must be entered as mm/dd/year
(e.g., 01/01/2004). The first date in the time function should correspond with the Start
Date specified in the Data Set Screen (Section 4.1.2). The last date in the time function
normally will correspond with the End Date.
Time - The time must be entered as hh:mm (e.g., 14:30).
Value - The inflow for this date and time is specified in units of [m3/sec]. Different units
can be used if a conversion factor is provided with the Flow Field (Section 4.3.1). Note
that if a downstream boundary elevation function is being defined, then the value entered
will be surface elevation [m].
The time function table allows the user to enter time variable flow information.
For constant flows, two rows should be specified with the simulation start and end dates,
and the constant flow value. The user can enter the information by hand, paste in from a
spreadsheet, or query in from database/spreadsheets.
5 Stream Transport Model Outputs
A set of output display variables is defined for each water quality module. Calculated
values of selected display variables are stored for every segment at every output time
interval through a simulation. Display variables are selected in the Output control screen.
Variables checked in the "Output" box will be available to the WASP graphical post-
processing software MOVEM. For each variable with a checked "CSV" box, WASP will
produce a separate comma-delimited file containing output for all segments and all output
times.
A subset of transport variables is included in the output for all of the water quality
modules. These are listed in Table 3. Users are encouraged to explore patterns and
28
-------
relationships among these variables to better understand the dynamics controlling
transport in their water body.
Table 3 - Transport output variables
Variable
Flow In
Flow Out
Ave Segment Flow
Water Velocity
Volume
Depth
Surface Elevation
Surface Width
Residence Time
Maximum Time Step
Actual Time Step
Units
m3/sec
m /sec
m /sec
m/sec
m3
m
m
m
days
days
days
Module
Heat
X
X
X
X
X
X
X
X
Eutr
X
X
X
X
X
X
X
X
X
Adv
Eutr
X
X
X
X
X
X
X
X
X
X
Toxi
X
X
X
X
X
X
X
Adv
Toxi
X
X
X
X
X
X
X
Merc
X
X
X
X
X
X
X
X
6 References
Ambrose, R. B., J. L. Martin, and T. A. Wool, 2006. Wasp? Benthic Algae - Model
Theory and User's Guide. U.S. Environmental Protection Agency, Washington, DC,
EPA/600/R-06/106 (NTIS PB2007-100139), 2006.
Wool T.A., R.B. Ambrose, J.L. Martin, and E.A. Comer, 2001. The Water Quality
Analysis Simulation Program, WASP6; Part A: Model Documentation. U.S.
Environmental Protection Agency, Center for Exposure Assessment Modeling, Athens,
GA.
Ambrose, R.B., Jr., T.A. Wool, and J.L. Martin, 1993. The Water Quality Analysis
Simulation Program, WASPS; Part A: Model Documentation. Internal Report
Distributed by USEPA Center for Exposure Assessment Modeling, U.S. Environmental
Protection Agency, Athens, GA.
Ambrose R.B., T.A. Wool, J.P. Connolly, and R.W. Schanz, 1988. WASP4, A
Hydrodynamic and Water Quality ModelModel Theory, User's Manual, and
Programmer's Guide. EPA/600/3-87-039, U.S. Environmental Protection Agency,
Athens, GA.
29
-------
Di Toro DM, JJ. Fitzpatrick, and R.V. Thomann, 1983. Water Quality Analysis
Simulation Program (WASP) and Model Verification Program (MVP) - Documentation.
Contract No. 68-01-3872, Hydroscience Inc., Westwood, NY, for U.S. EPA, Duluth,
MN.
Chapra, S.C. 1997. Surface Water-Quality Modeling, McGraw-Hill, New York, New
York, 844 pp.
Chapra, S.C. 2003. QUAL2K: A Modeling Framework for Simulating River and Stream
Water Quality (Beta Version): Documentation and Users Manual. Civil and
Environmental Engineering Dept, Tufts University, Medford, MA.
Chapra, S.C. and GJ. Pelletier. 2004. QUAL2K: A Modeling Framework for Simulating
River and Stream Water Quality, Version 1.3: Documentation and Users Manual. Civil
and Environmental Engineering Dept., Tufts University, Medford, MA.
Feigner, K.D. and H.S. Harris, 1970. Documentation Report - FWQA Dynamic Estuary
Model. U.S. Department of Interior, Federal Water Quality Administration.
Wool, T.A., R.B. Ambrose, and J. L. Martin. 2001. "The Water Analysis Simulation
Program, User Documentation for Version 6.0," Distributed by USEPA Watershed and
Water Quality Modeling Technical Support Center, Athens, GA.
30
-------
7 Appendix 1 : Derivation of Equations
7. 1 Hydraulic Exponents for Kinematic Wave Flow
Manning's formula (Equation 19) provides the basis for deriving relationships among the
hydraulic exponents. Rearranging terms gives cross-sectional area as a function of flow:
Equation 54
A ,/,3/5 0-3/10 n2/5 ^3/5
A = n o0 B ()
Cross-sectional average depth (hydraulic radius) is given by:
Equation 55
d = - = n3'5 S0mo B-3'5 Q3/5
B °
Velocity is given by:
Equation 56
v = Q- = n-3/5 S3/w B-2'5 Q2/5
Substituting Equation 3 for width as a function of flow gives:
Equation 57
d = n315 S03/w
Equation 58
3/w
v = n-35 S30
Equation 57 and Equation 58 give the depth and velocity hydraulic exponents as a
function of the width exponent:
Equation 59
dxp = 0.6 (l-6exp) vexp = 0.4 (l-6exp)
Comparing the flow exponents confirms that the velocity exponent is 2/3 of the depth
exponent, confirming Equation 9.
8 Appendix 2: Model Verification Tests
Model verification tests were designed to assure that the equations are
implemented correctly in the model code. Results are stored in separate folders at:
\WASP7\QA\Stream Transport 4-Stream Kinematic Wave Flows\
Verification tests are outlined below. Results are detailed in a companion document.
31
-------
8.1 Kinematic Wave Tests
8.1.1 Stream Transport Test 1
This series tests the kinematic wave flow routines in WASP using the Simple
Toxicant module with intermediate slope and steady inflow with an example problem
taken from Chapra, Example 14.6, pp. 253, 254. Results are compared with analytical
solutions implemented in a spreadsheet, and by simple hand calculations.
Test la - Upstream Inflow Only
Test Ib - Upstream and Pore Water Inflow
Test Ic - Upstream and Precipitation Inflow
Test Id - Upstream Inflow and Evaporation Outflow
8.1.2 Stream Transport Test 2
This series tests the kinematic wave flow routines in WASP using the Heat
module with shallow slope and step changes in inflow. Results for each flow step are
compared with analytical calculations calculated in a spreadsheet.
Test 2a - Rectangular Cross-Section
Test 2b - U Cross-Section
Test 2c - V Cross-Section
8.1.3 Stream Transport Test 3
This series tests the kinematic wave flow routines in WASP using the
Eutrophication module with shallow slope, U-shape cross-section, step changes in
upstream inflow, and inflow or withdrawal at Segment 3.
Test 3a - Constant Tributary Inflow
Test 3b - Constant Flow Withdrawal
8.1.4 Stream Transport Test 4
This series tests the kinematic wave flow routines in WASP using the Mercury
module with a branching stream system. A medium stream with moderate slope is
connected to a small tributary with shallow slope.
Test 4a - Step Inflows
Test 4b - Long-Term, Variable Inflows
Test 4a specifies step changes in upstream and tributary inflows. Results for each flow
step are compared graphically and to analytical solutions from a spreadsheet. Test 4b uses
variable upstream and tributary inflows repeating in a 2-year pattern over a long
simulation period. Results are examined graphically for stationary (repeating) output.
32
-------
8.1.5 Stream Transport Test 5
This tests the kinematic wave flow routines in WASP using the Heat module with
a diverging-converging stream system. A small river with steep slope diverges into two
branches receiving 40% and 60% of the upstream flow. These branches converge
downstream. This test uses step changes in the upstream inflow. Results for each flow
step are compared graphically and to analytical solutions from a spreadsheet.
8.2 Weir Overflow Verification Tests
Model verification tests were designed to assure that the equations are
implemented correctly in the model code. Results are stored in separate folders at:
\WASP7\QA\Stream Transport\4-Ponded Weir Flows\
Tests are outlined below.
8.2.1 Weir Overflow Test 1 - Steady Flow
This tests the ponded weir overflow routine in WASP using the Mercury module
with steady inflow and sequentially increasing weir heights. Results are compared with
analytical calculations.
8.2.2 Weir Overflow Test 2 - Variable Flow
This tests the ponded weir overflow routine in WASP using the Mercury module
with sequentially increasing inflow and sequentially increasing weir heights. Results are
compared with analytical calculations.
8.2.3 Weir Overflow Test 3 - Long Term Dynamics
This tests weir overflow routines in WASP for long-term performance with
variable inflows repeating in a 2-year pattern. Results are examined for long-term drift in
depths, volumes, and velocities.
8.3 Dynamic Flow Verification Tests
Model verification tests were designed to assure that the equations are
implemented correctly in the model code. Results are stored in separate folders at:
\WASP7\QA\Stream Transport\4-Dynamic Flows\
Tests are outlined below.
8.3.1 Dynamic Flow Test 1 - Steady Backwater
This tests the dynamic flow routine in WASP with mild slope, steady inflow, and
downstream pond with weir. Weir height is set to provide downstream ponded depth
33
-------
equal to the kinematic flow depth. Results are compared with output from an equivalent
kinematic wave simulation.
8.3.2 Dynamic Flow Test 2 - Steady Flow, Elevation
This tests the dynamic flow routine in WASP with mild slope, steady inflow, and
constant downstream boundary elevation. Results are compared with output from an
equivalent DYNHYD simulation linked with WASP.
8.3.3 Dynamic Flow Test 3 - Variable Stream Flow
This tests the dynamic flow routine in WASP using the Mercury module with
sequentially increasing flow. Results are compared with output from a DYNHYD
simulation linked with WASP.
8.3.4 Dynamic Flow Test 4 - EFDC Stream
This tests the dynamic flow routine in WASP with mild slope, steady inflow, and
constant downstream boundary elevation. Results are compared with output from an
equivalent EFDC simulation linked with WASP.
8.3.5 Dynamic Flow Test 5 - Tidal Stream
This tests the dynamic flow routine in WASP with flat slope, no inflow, and
sinusoidal tidal downstream boundary elevation. Results are compared with output from
equivalent DYNHYD and EFDC simulations linked with WASP.
8.3.6 Dynamic Flow Test 6 - 2-D Tidal Stream
This tests the dynamic flow routine in WASP on a 2-dimensional network with
flat slope, no inflow, and sinusoidal tidal downstream boundary elevation. Results are
compared with output from an equivalent EFDC simulation linked with WASP.
34
-------
vvEPA
United States
Environmental Protection
Agency
Office of Research
and Development (8101R)
Washington, DC 20460
Official Business
Penalty for Private Use
$300
EPA/600/R-09/100
September 2009
Please make all necessary changes on the below label,
detach or copy, and return to the address in the upper
left-hand corner.
If you do not wish to receive these reports CHECK HERE
D; detach, or copy this cover, and return to the address in
the upper left-hand corner.
PRESORTED STANDARD
POSTAGE & FEES PAID
EPA
PERMIT No. G-35
Recycled/Recyclable
Printed with vegetable-based ink on
paper that contains a minimum of
50% post-consumer fiber content
processed chlorine free
------- |