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 |