United States
Environmental Protection
Agency
Municipal Environmental Research EPA-600/2-78-148
Laboratory August 1978
Cincinnati OH 45268
Research and Development
Stream Models
for Calculating
Pollutional Effects
of Stormwater Runoff
-------
EPA-600/2-78-148
August 1978
STREAM MODELS FOR CALCULATING
POLLUTIONAL EFFECTS OF STORMWATER RUNOFF
by
Robert Smith and Richard 6. Eilers
Wastewater Research Division
Municipal Environmental Research Laboratory
Cincinnati, Ohio 45268
MUNICIPAL ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
CINCINNATI, OHIO 45268
-------
DISCLAIMER
This report has been reviewed by the Municipal Environmental Research
Laboratory, U.S. Environmental Protection Agency, and approved for publica-
tion. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
-------
FOREWORD
The Environmental Protection Agency was created because of increasing
public and government concern about the dangers of pollution to the health
and welfare of the American people. Noxious air, foul water, and spoiled
land are tragic testimonies to the deterioration of our natural environment.
The complexity of that environment and the interplay between its components
require a concentrated and integrated attack on the problem.
Research and development is that necessary first step in problem
solution, and it involves defining the problem, measuring its impact, and
searching for solutions. The Municipal Environmental Research Laboratory
develops new and improved technology and systems to prevent, treat, and
manage wastewater and solid and hazardous waste pollutant discharges from
municipal and community sources, and to minimize the adverse economic,
social, health, and aesthetic effects of pollution. This publication is
one of the products of that research—a most vital communications link
between the research and the user community.
The information presented here describes the use of computer models to
simulate the biological, physical, and chemical reactions that occur in a
flowing stream. Hydraulic aspects are also considered. These models make
it possible to estimate the pollutional effect of stormwater runoff on
receiving streams. As such, this study fulfills the need for continuing
technology assessment in emerging areas.
Francis T. Mayo
Director, Municipal Environmental
Research Laboratory
m
-------
ABSTRACT
Three related studies are described that provide the means to quantify
the pollutional and hydraulic effects on flowing streams caused by storm-
water runoff. Mathematical stream models were developed to simulate the
biological, physical, chemical, and hydraulic reactions that occur in a
stream. Relationships take the form of differential equations with the two
independent variables of time and distance. The differential equations can
be solved directly by means of calculus or by digital computer using
numerical methods. The solution would be the concentration of species of
pollutional interest, such as BOD and dissolved oxygen, within the stream
as a function of distance and time. The solution can be steady-state or
transient. There is sufficient information presently available for finding
steady-state solutions. However, when the pollution loads and/or the
initial conditions for the flowing stream vary with time, the problem
becomes much more difficult, and the technology for handling the transient
situation has not been adequately developed. The purpose of this report
is to show how the solution can be found for the case where the pollution
loading is a transient, especially as it applies to stormwater overflow.
The work described in this report was done inhouse by the Systems and
Economic Analysis Section of EPA and covers a period from April 1977 to
November 1977 and the work was completed as of December 1977.
-------
CONTENTS
Foreword i i i
Abstract iv
Figures vi
Tables viii
Metric Conversion Table ix
1. Introduction 1
2. Summary and Conclusions 3
3. Recommendations 5
4. Stormwater Overflow Pollution Stream Model (SWOPS) 6
Background 6
Closed-Form Solutions 8
Solution by Numerical Integration 11
Accuracy of Numerical Integration Solution 19
Summarizing the Impact of Square Pollution Pulses ... 22
References for Secti on 4 31
Appendix for Section 4 32
5. Stormwater Overflow Hydraulic Stream Model (SWOHS) 43
Background 43
Stream Conceptualization 43
Physical Laws Used 48
Systems of Uni ts Used 48
Conservation-of-Mass 49
Conservation-of-Momentum 49
Numerical Integration Scheme Used 50
Description of Input Required 51
Logical Structure of Program .... 53
Solution for Hypothetical Problems 56
Program Listing 60
References for Section 5 72
Appendix for Section 5 73
6. Effect of Stormwater on Stream Dissolved Oxygen 81
Background 81
Problem of Deoxygenation and Reaeration in Tubular
Streams 83
Stream Response with Dispersion 84
Non-Dimensional Groupings Used to Summarize Computed
Resul ts 86
Summary and Conclusions 92
References for Section 6 95
Appendix for Section 6 96
-------
FIGURES
Number Page
1 Diagrams for explanation of numerical integration used
in stream model for study of storm overflow events 12
2 Diagram illustrating how a storm overflow hydrograph
can be divided into discrete flow rates for use in the
stream model 13
3 Computed BOD profiles in the stream caused by 24-hour
storm overf 1 ow event 20
4 Dissolved oxygen deficit profiles in the stream caused
by 24-hour storm overflow event 21
5 Computed BOD profiles in the stream caused by 1-hour
storm overf 1 ow event 23
6 Dissolved oxygen deficit profiles in the stream caused by
1-hour storm overf1ow event 24
7 Computed standard deviations for dissolved oxygen deficit
profiles at time from entry of leading edge of pulse for
dispersion coefficient volumes of 0.1, 0.2, and 1.0 sq mi/day 25
8 Correction factor for peak dissolved oxygen deficit to be
applied to the value computed with the Streeter-Phelps
model 27
9 Time between entry of leading edge of square stormwater
pulse and time of peak DO deficit in terms of *
non-dimensional groupings 28
10 Diagram of solution reach used in SWOHS showing cross-
sections and basic equations used for depth, cross-sectional
area, wetted perimeter, hydraulic radius, and surface area 44
11 Computed wave profiles caused by a square hydrograph
entering a steady-state stream with 598.2 cfs volume flow 57
VI
-------
FIGURES (continued)
Number Page
12 Computed flow over terminal weir minus steady-state flow
resulting from a square hydrograph introduced 11.4 miles
upstream of the wei r 58
13 Diagram of second hypothetical problem solved with
the SWOHS program 61
14 Steady-state variation of water surface elevation in
the second hypothetical problem solved with the SWOHS
program 62
15 Typical stormwater overflow event divided into one-half
hour sub-loads 82
16 Comparison between computed dissolved oxygen deficit
peak values using the program SWOPS and the closed
form solution 89
17 Correction factor for peak dissolved oxygen deficit to
be applied to the value computed with the Streeter-Phelps
model 93
18 Relationships for finding the time to the peak dissolved
oxygen deficit from the non-dimensional grouping N 94
-------
TABLES
Number Pa9
1 Computed Results 29
2 Definitions of FORTRAN Symbols Used in SWOPS 32
3 FORTRAN Source Listing for the SWOPS Program 35
4 Sample Printout for SWOPS (DL = 0, PRNT =1) 41
5 Sample Printout from SWOPS (DL = 1, PRNT = 1) 42
6 Definitions for Variables Used in the SWOHS Program 45
7 Printout for First Hypothetical Problem Solved with SWOHS 63
8 Printout for Second Hypothetical Problem Solved with SWOHS .... 65
9 Fortran Listing for the SWOHS Program 66
10 Definitions of FORTRAN Symbols Used in HYSS 75
11 FORTRAN Source Listing for HYSS 77
12 Printout for the First and Second Hypothetical Problems Using
HYSS 80
13 Listing for the Program SWOCFS which Applies the Superposition
Principle to the Closed-Formed Solution for an Impulsive BOD
Load 87
14 Sample Printout from SWOCFS 90
15 Computed Values for Peak Dissolved Oxygen Deficit and Time of
Occurrence for Square Pollutographs 91
-------
METRIC CONVERSION TABLE
Customary Units Multiplier
cubic foot, cu ft .02832
cubic foot, cu ft 28.32
cubic feet per minute, cfm .04719
cubic feet per second, cfs .02832
foot, ft .3048
gallon, gal 3.785
mile, mi 1.609
miles per hour, mi/hr 1.609
million gallons, mil gal 3785.0
million gallons per day, mgd 43.81
million gallons per day, mgd .04381
pound, Ib .4536
square miles, sq mi 2.590
Metric Units Reciprocal
cubic meter, m 35.31
liter, 1 .03531
liters per second, I/sec 2.119
3
cubic meters per second, m /sec 35.31
meter, m 3.281
liter, 1 .2642
kilometer, km .6215
kilometers per hour, km/hr .6215
cubic meters, m .0002642
liters per second, I/sec .02282
cubic meters per second, m /sec 22.82
kilogram, kg 2.205
2
square kilometers, km .3861
-------
SECTION 1
INTRODUCTION
This report consists of separate descriptions of the following three
related studies that were performed inhouse by the U.S. Environmental
Protection Agency:
1. Stormwater Overflow Pollution Stream Model (SWOPS)
2. Stormwater Overflow Hydraulic Stream Model (SWOHS)
3. Effect of Stormwater on Stream Dissolved Oxygen
These studies were made in order to develop better methods for quantifying
the pollutional effects in streams caused by Stormwater overflows.
Mathematical stream models were developed to simulate the biological,
physical, and chemical reactions that occur in a flowing stream. Hydraulic
aspects are also taken into account. The relationships take the form of
differential or difference equations with two independent variables: time
and distance along the stream. The differential equations can be solved
with analog or digital computers or, in some cases, calculus can be used
to find analytical (closed-form) solutions. The solution of the differential
equations is the concentration of all species of pollutional interest, such
as BOD and dissolved oxygen, within the stream as a function of distance
and time resulting from the assumed pollution load on the stream.
The solution is divided into the transient regime and the steady-state
regime. The steady-state solution is simply the solution that will apply
when time becomes infinitely large, and therefore the steady-state solution
is a function of distance only. The transient solution is the concentration
of all pollution species as a function of distance and time beginning at
time zero and ending when the transient solution equals the steady-state
solution.
The pollution loads on the stream can also be specified as steady-
state or transient. A steady-state load on the stream is one where the
rate of pollution entering the stream is constant. A transient load on the
stream is one where the rate of pollution entering varies with time. It
should be pointed out that the solution is required only for a specified
segment of the stream, and the range of interest for time is from zero to
the time at which the solution reaches the steady-state regime.
-------
When the initial conditions for the river at distance zero are constant
for all time and the pollution loads on the stream are steady-state, the
solution of interest is the steady-state solution. This set of assumptions
is made in most stream modeling, and many analytical expressions and computer
programs are available for finding the steady-state solution. However, when
the pollution loads and/or the initial conditions for the river at
distance zero vary with time, the problem is much more difficult and the
technology for treatment of this transient situation is poorly developed.
The purpose of these reports is to show how the solution can be found for
the case where the pollutional load is a transient, especially as it
applies to storm overflow events.
-------
SECTION 2
SUMMARY AND CONCLUSIONS
Since stormwater overflows from urban areas served by combined sewers
contribute a mass of pollutants to the stream roughly equivalent to the
sanitary sewage treated to a secondary level, environmental planners must
estimate the effect of stormwater on the dissolved oxygen level in the
receiving stream. Because of the transient nature of the problem, computa-
tional methods available were complex and time consuming. Although the
Streeter-Phelps model neglects the very important dispersion term, it is
often used for estimating the impact of stormwater overflows because of the
simplicity of the solution.
In an effort to improve the technology for assessing the importance of
stormwater overflows a digital computer stream model (SWOPS) was developed
for computing the dissolved oxygen (DO) deficit in the stream as a function
of time and distance along the stream caused by any specified overflow. The
accuracy of the computation was judged by comparsion with the known closed-
form solution for the BOD concentration in the stream resulting from an
impulsive load. The superposition principle was used to find the BOD response
to any overflow from the closed-form solution for the impulsive load. It was
discovered that the computational accuracy was greatly improved by the use of
the Lagrange coordinate system instead of the commonly used Euler coordinate
system. A second digital computer program (SWOHS) was developed to study the
hydraulic effect on the stream of large volume overflows. Since the reaera-
tion coefficient is known to be related to the stream velocity and water depth
the hydraulic aspects of the problem can affect the reaeration coefficient and
thus the dissolved oxygen levels in the stream. Computations made with SWOHS
show that the hydraulic effect is likely to be minimal in all except the most
extreme cases.
It was found that a time interval of 15 minutes gave sufficient accuracy
for computing the DO deficit profiles with SWOPS. Dissolved oxygen deficit
profiles were computed with SWOPS using a wide range for all variables and
a method for presenting the results in a succinct form was sought. It was
found that the ratio of the peak DO deficit in the stream to the peak DO
deficit computed with the Streeter-Phelps model could be presented as a single-
valued function of a non-dimensional grouping. It was also found that if the
time-to-the-peak was expressed as a non-dimensional grouping this grouping
could also be presented as a single-valued function of the first non-dimen-
sional grouping. Thus, a method was found for summarizing the results of the
computations in an easily used manner.
-------
In performing the computations with SWOPS it was noted that as the width
of the square pollution pulse was decreased the form of the DO deficit pro-
files approached the shape of the normal probability density function. Using
this fact as a clue, a closed-form solution for the DO deficit response to an
impulsive BOD load was then discovered. This closed-form solution can now be
used with the superposition principle to find the DO response in the stream
to any stormwater pollution event. A short digital computer program (SWOCFS)
was then developed to apply the closed-form solution for the impulsive load
with the superposition principle although the computations can be made by
hand if necessary. Thus, the technology for estimating the impact of any
stormwater overflow on the dissolved oxygen resources of the stream has been
advanced to a point where making the computations is feasible even by hand
computation.
-------
SECTION 3
RECOMMENDATIONS
The effect of storm and combined sewer overflows on the quality of the
receiving stream is poorly understood. One reason for this is that the
effect of time varying or transient concentrations of dissolved oxygen
deficit or toxicants on aquatic life is difficult and expensive to investi-
gate. Another reason is the difficulty of calibrating and validating a time
dependent stream model. The prevailing opinion at this time is that in
flowing streams of reasonably large volume the impact of periodic dissolved
oxygen depressions caused by storm and combined sewer overflows is probably
not critical to aquatic life. It is also understood that in lakes or even
estuaries the impact of storm overflows can be critical. There is some
evidence that storm overflow solids can settle out on the stream bottom and
build up benthic deposits which can have a significant oxygen demand in dry
weather or low flow periods which may seriously deplete the oxygen resources
of the stream. These benthic deposits can also become resuspended in high
flow periods to cause significant BOD loads in the stream. This aspect of
the storm water pollution problem needs to be investigated first by locating
and investigating sites where this mechanism occurs and second by devising
a stream model capable of predicting the magnitude of the effect based on
the physical parameters of the stream and the pollution sources.
-------
SECTION 4
STORMWATER OVERFLOW POLLUTION STREAM MODEL (SWOPS)
BACKGROUND
The Stormwater Overflow Pollution Stream Model (SWOPS) is a mathematical
model for a natural flowing river or stream developed primarily to study the
effect of dispersion within the stream on transient changes in water quality
caused by storm and combined sewer overflow events. The stream is assumed to
be prismatic, that is, the cross-sectional area and velocity are fixed for
all distances along the stream.
The length of the river to be studied, the solution reach, is divided
into equal intervals of distance (AX) along the stream. Processes such as
advection, dispersion, depletion of oxygen reserves by BOD degradation, and
reaeration from atmospheric air are expressed as difference equations with
time (t) and distance (X) as the independent variables. These difference
equations can be expressed as partial differential equations when the incre-
ments of time and distance are allowed to approach zero.
The point along the stream where the storm overflow occurs is specified
and the characteristics of the overflow are expressed as a volume flow (cfs)*
and a concentration of pollutant (BOD) for each time increment beginning at
time equals zero. The solution to the differential equations over the domain
of time and distance along the stream is found by means of numerical integra-
tion using the Crank-Nicolson method. Lagrange coordinates, using a distance
reference point moving downstream at the stream velocity, are used instead of
the more conventional Euler coordinates, using a distance reference point
fixed on the stream bank, to improve the accuracy of the numerical solution.
The accuracy of the numerical integration scheme is checked against a closed
form solution for BOD concentration in the stream caused by a square p^lse
of BOD pollution entering the stream. Finally, the computed results are
presented as functions of non-dimensional groupings to allow the analyst to
estimate the effect of dispersion without use of the digital computer program.
The solution is roughly divided into two regimes; the transient regime
and the steady-state regime. The steady-state solution is simply the solution
*For simplicity, customary units of measurement are being used in this report
instead of the required metric units, and a metric conversion table is given
on page ix.
-------
which will apply when time becomes infinitely large and, therefore, the steady
state solution is a function of distance only. The transient solution is the
concentration of all pollutional species as a function of distance and time
beginning at time zero and ending when the transient solution equals the
steady-state solution.
The pollutional loads on the stream can be specified as steady-state or
transient. For example, a steady-state load on the stream is one where the
rate of pollution entering the stream (Ib/day) is constant. A transient
load on the stream is one where the rate of pollution entering varies with
time. An example of a transient load is a storm overflow in the form of a
square pulse of four hours duration with a fixed concentration of BOD.
It should also be pointed out that the solution is required only for a
specified segment of the river, the solution reach corresponding to the
distance from X-0 to X=XMAX. Similarly, the range of interest for time is
restricted from t=0 to the time at which the solution reaches the steady-state
regime.
When the initial conditions for the river upstream of the pollution
sources are constant for all time and the pollutional loads on the stream are
steady-state the solution of interest is the steady-state solution. This
set of assumptions is made in most river modeling and many analytical expres-
sions and digital computer programs are available for finding the steady-state
solution. However, when the pollutional loads and/or the initial conditions
for the river vary with time the problem is more difficult and the technology
for treatment of this transient problem is not well documented. The purpose
of this report is to show how the solution can be found for the case where
the pollutional load is a transient especially as it applies to storm over-
flow events.
The solution for storm overflow events will be found by means of
numerical integration of the basic differential equations. However, errors
of various kinds are indigenous to numerical integration and it is, there-
fore, necessary to compare the numerical integration solution to a solution
known to be mathematically accurate in order to develop a knowledge of the
limitations of the numerical solution. Analytical solutions (closed-form)
for a number of simple stream modeling problems are available. First of
these is the Streeter-Phelps (1) steady-state solution for BOD and dissolved
oxygen profiles in a prismatic stream assuming that no dispersion exists.
Analytical expressions are also available (2) for steady-state BOD and dis-
solved oxygen profiles in a prismatic stream including dispersion. However,
it was shown by Dobbins (2) that for the range of dispersion coefficients
normally experienced in natural streams the Streeter-Phelps solution is a good
approximation to the steady state solution including dispersion. A closed-
form expression is available (3) for calculating the concentration of BOD in
a stream resulting from an impulsive load of BOD entering the stream as a
function of time and distance from the dump point. By an impulsive load we
mean a load which enters the stream instantaneously. The Streeter-Phelps
solution and the solution for an impulsive load were used to test the accuracy
of the numerical integration scheme used in SWOPS. These two closed form
solutions will be derived in the following section.
-------
CLOSED-FORM SOLUTIONS
The pioneering work of H. W. Streeter and Earle B. Phelps of the U. S.
Public Health Service is presented in Public Health Bulletin No. 146 (1)
dated February, 1925. They visualized the Ohio River as many segments
divided by planes perpendicular to the flow of the river which move in the
direction of flow with no mixing across the dividing planes. Biological
activity, using the BOD pollutional load, depletes the dissolved oxygen
reserves of the river and additional oxygen enters each segment from the
atmosphere at a rate directly proportional to the dissolved oxygen deficit.
The Streeter-Phelps model is described by the following two differential
equations.
dL/dt = -KrL (1)
dD/dt = KdL - KaD (2)
where L = 5-day BOD concentration, mg/1
t = time of travel, days
K = rate constant for deoxygenation
-1
and sedimentation, day
D = dissolved oxygen deficit, mg/1
«d = deoxygenation rate constant, day
Ka = reaeration rate constant, day
LQ = BOD concentration at X=0 for
all time, mg/1
Equation (1) can be integrated simply as follows:
L = L^V (3)
If the volume flow of the pollutional load is small compared to the volume
flow of the river LQ is computed simply as follows:
Lo = QI Li/Q (4)
Q.J = volume flow of pollutional stream, cfs
Q = volume flow of river, cfs
L.J = 5-day BOD concentration in pollutional
stream, mg/1
-------
Solving equation (2) is slightly more involved than solving equation (1).
Equation (2) can be written in the standard form for a linear first order
differential equation as follows:
dD/dt + KD = K,L
a u
(5)
K t
If we multiply both sides of equation (5) by the integrating factor, e a
K t
it can be seen that the left-hand side is an exact derivative of De a
Therefore, equation (5) can be written as follows:
D
-,t A
= Kd LQ ye(Ka"Kr)
— '0 o
dt
(6)
This equation is readily integrated to give the following result:
D - K
-Kr)
H- D
(7)
In equation (7), D is the initial dissolved oxygen deficit in the
stream as it enters the solution reach. Equations (3) and (7) can be used to
check the numerical integration solution in certain cases. If it is assumed
that the cross-sectional area of the stream and the stream velocity are
constant with distance, time and distance are related by X = time x velocity
and equations (3) and (7) can be used to find the steady-state solution.
It will be shown later that as the length of the pollution pulse
approaches 24 hr in length the effect of dispersion is minimized and the
Streeter-Phelps solution gives a good approximation to the true transient
solution. For shorter pulses, the solution can only be found by means of
numerical integration. A partial check on the accuracy of the numerical
solution can be had by comparison with the closed-form solution for BOD
concentration in the stream resulting from an impulsive BOD load.
The following expression for the BOD concentration in the stream caused
by an impulsive or spike load of BOD entering the stream can be derived by
Laplace transforms as shown by Thomann (3).
L(x,t) ='C exp(-%((x-vt)/s)) exp(-Kft)
(8)
L = BOD concentration in stream, (cu mi)
x = distance from spike input, miles
t = time from spike input, days
v = stream velocity, miles/day
~
-------
K = rate constant for BOD removal, days"
C = I/(As /2II )
s = (2Et)h
E = dispersion coefficient, sq miles/day
A = cross-sectional area of stream, sq miles
Equation (8) expresses the concentration of BOD in the stream as Ib/cu
mi per pound of BOD introduced in the spike input.7 Pounds per cubic mile
can be converted to mg/1 by the factor 1.008 x 10 . The stream cross-
sectional area (A) equals Q/v where Q is the stream volume flow in cubic
miles/day and v is the stream velocity in miles/day. Cubic feet per second 7
(cfs) for stream flow can be converted to cu mi/day by the factor 5.87 x 10~ .
Thus, equation (8) can be expressed as mg/1 of BOD per pound of BOD introduced
in the spike as follows:
L(x,t) = (0.07395/(sQ/v)) exp(-%((x-vt)/s)2) exp(-Kft) (9)
To find the expression for BOD in the stream (L) as a function of time
and distance equation (9) must be multiplied by the mass of pollutant intro-
duced. If the combined sewer overflow is a square pulse with a volume flow
of QIN (cfs) and a BOD concentration of LOIN (mg/1) and a duration of TAU
(days) the amount of pollutant is found as follows:
Ib 5-day BOD - QIN*0.646*LOIN*TAU*8.33 (10)
Multiplying equation (9) by equation (10) finally gives an expression for BOD
in the stream as a function of time and distance.
L, mg/1 = (0.3979/(Qs/v))*QIN*LOIN*TAU exp(-M(x-vt)/s)2)exp(-Krt) (11)
It can be seen from equation (11) that the BOD profile in the stream at any
time (t) has the form of the normal probability density function with a mean
of stream velocity times time from the spike and a standard deviation of
(2Et)*s.
Equation (11) cannot be used directly to check the results of the
numerical integration because it applies only to a mass of pollutant intro-
duced instantaneously. However, since the differential equations for the
processes occurring within the stream are linear any square pulse can be
divided up into sub-pulses and the responses to the sub-pulses can be summed
to find the response to the square pulse. A simple computer program named
STREAM-4 (4) has been devised to perform this type of computation. The
results from STREAM-4 can then be used to check the accuracy of the numerical
integration scheme.
10
-------
SOLUTION BY NUMERICAL INTEGRATION
The processes of advection, dispersion, deoxygenation caused by the
presence of BOD, and reaeration at the water surface can all be expressed as
difference equations. The corresponding partial differential equations can
then be found by letting the increments of time and distance approach zero.
The difference equations are solved simultaneously by means of a
numerical integration scheme to find the solution which consists of computed
values for each of the two dependent variables, BOD (L) and dissolved oxygen
deficit (D) at each of the lattice points DX apart in the X direction and DT
apart in the time direction. This is shown diagrammatically in Figure Ic.
Values for the dependent variables must be supplied for all points along the
time=0 axis and for all points along the X=0 axis which represents the
upstream boundry of the solution reach. At the downstream boundry of the
solution reach, shown in Figure Ic, the values for L and D are assumed to be
zero. In addition, the characteristics of the pollutional pulse and the point
along the X axis at which it enters must be specified. The pollutional pulse
is assumed to begin at time=0.
The first step is to divide the solution reach into equal increments
(AX) by planes perpendicular to the stream velocity vector. The volumes
between the dividing planes are sometimes called control volumes. This is
shown in Figure la and three adjacent control volumes are shown in Figure Ib.
If the volume of each increment is V and the volumetric flow rate for the
stream is Q, advection of BOD (L) and dissolved oxygen deficit (D) into and
out of control volume (2) is expressed by the following two difference
equations:
V AL2/At = Q (L1 - L2) (12)
V AD2/At = Q (D1 - D2) (13)
These two expressions apply only when the Euler coordinate system is used.
The Euler coordinate system assumes that the incremental volumes are fixed
with respect to the shore and the stream flows through each volume.
When the Lagrange coordinate system is used the X=0 point is fixed with
respect to the contents of one of the control volumes at time=0. Therefore,
in the Lagrange system no advection will occur in the tubular stream, and the
relationships for advection shown in equations (12) and (13) can be dropped
from the set of equations to be solved.
When the Euler coordinate system is used the pollutional pulse will
enter a single control volume and will cause a positive rate of change of BOD
in that control volume. For example, assume that an outfall is located at
control volume (2). When the hydrograph for the outfall is known, such as
the unit hydrograph shown in Figure 2, the hydrograph is divided into incre-
ments of At and the flow assigned to each increment is the average of the
11
-------
to
to
•r—
X
CD
DX
3
-o-
Ib
la
) (.
Be
\
\ t
) (.
) f
•\ f
) (,
&f
) (,
\ f
> (.
} (•
~\ r
> (
\ t
) (.
) f
1
(.
\ f
i \
Df
\.
\ C
1 DX -
5f
\
V f
) I
\ /
} \
» £
•< —
I
1
h
c
1 1
1
t
\ f
) (.
i c
3
<
N
)
\ —C*i
) \)
\ C i C*-
Distance axis, miles
Figure 1. Diagrams for explanation of numerical integration used in
stream model for study of storm overflow events.
12
-------
T
\
-------
logicalractivity combined. The rate constant for loss of BOD by biological
flow at the boundries of At. Thus, the hydrograph is replaced with a bar-_
graph as shown in Figure 2. When a plot of BOD concentration versus time is
known, this is also divided into At increments in a similar way. The average
volume flow from the hydrograph for each At is called Qin- Similarly,
the average BOD concentration over each At is called L, . Using the Euler
coordinates, all input of pollution will occur at segment (2) and the rate
of change of BOD can be written as follows where Qin and Lin are functions
of At.
V AL2/At = Qin L-n (14)
When Lagrange coordinates are used, the ratio of AX to At must always
equal the stream velocity. When this is true, it can be seen that pollution
associated with the first At interval of the hydrograph will enter increment
(N), pollution associated with the second At will enter increment (N-l), etc.
where increment (N-l) is adjacent to and upstream of increment (N).
Equations for simulating the loss of BOD and the loss and replenishment
of dissolved oxygen in the stream are the same for both sets of coordinates.
For example, the loss of BOD in the stream in written simply as follows:
V AL2/At = -KrL2 V (15)
Here, K is the rate constant for the loss of BOD by sedimentation and bio-
logicalractivity comt
activity alone is K..
The use of dissolved oxygen by biological activity is expressed as
follows:
V AD2/At = Kd L2 V (16)
Reaeration of the stream through the water surface is expressed as
follows:
V AD2/At = -Ka D2 V (17)
Similar terms are included in the stream model for the oxygen consumption of
nitrogenous compounds and the benthic demand of sludge deposits on the stream
bottom but these are omitted here to simplify the discussion.
The final term to be considered is dispersion of BOD and dissolved
oxygen across the planes which separate the incremental stream volumes.
Although dispersion is known to consist of both eddy mixing and molecular
diffusion, the form of the equation used to simulate dispersion has the same
form as Fick's law for molecular diffusion. Pick's law states that the rate
of diffusion across the cross-sectional area (A) of the stream is equal to a
constant (E, sq miles/day) times the product of stream cross-sectional area
and the gradient of the concentration of the diffusion material along the
stream.
14
-------
Again, referring to the second segment, the rate of diffusion of BOD
upstream from the segment is written as:
rate of upstream diffusion, Ib/day = A E (L2 - L^/AX (18)
Similarly, the rate of diffusion downstream is:
rate of diffusion downstream, Ib/day = A E (L9 - U)/AX (19)
i— O
Since the Area A, (sq miles) can be expressed as V/AX the net rate of diffu-
sion into segment 2 is:
Net rate of diffusion 9 9
into segment 2 = EV(L] - L2)/AX + EV(I_3 - L2)/AX (20)
V AL/At = EV(L1 - 2L2 + L3)/AX2 (21)
The following analogous expression can be used to describe the diffusion
of dissolved oxygen deficit.
V AD/At = EV(D1 - 2D2 + D3)/AX2 (22)
If the diffusion coefficient (E) is different for BOD and dissolved
oxygen deficit we can use (E ) for BOD and (E .) for dissolved oxygen deficit.
By summing the right-hand sides of equations (12), (14), (15), and (21) the
complete difference equation for the rate of change of BOD in increment (2)
is written as follows:
AL2/At - Q(L1 - L2)/V + Qin Lin/V - KfL2 + E^L-, - 2L2 + L3)/AX2 (23)
Similarly, the right-hand sides of equations (13), (16), (17), and (22) can
be summed to find the rate of change of dissolved oxygen deficit in incre-
ment (2).
AD9/At = Q(D, - D9)/V + K,L9 K D9 + E .(D, - 2D9 + DJ/AX2 (24)
L I L.
-------
A crude solution of equations (23) and (24) can be found by taking the
values for L and D as those at the I'th time row to estimate values for L
and D at the (I+l)'th time row. Since the solution for L does not involve
D the solution for L should be computed first. If this approach is taken,
it is known from experience that the time and distance intervals must be
very small to achieve any accuracy and this results in extreme amounts of
computer time as well as magnification of computational error. A more
efficient approach is to use the sum of values at the I'th and the (I+l)'th
time points divided by two as the values to be used in equations (23) and
(24). To implement this type of solution, called the Crank-Nicolson method
by Smith (5), a linear equation must be written for each point along the X
axis within the solution reach and these equations must be solved to advance
the solution one time interval. Fortunately, the linear equations contain
only three unknowns except for the first and last which contain two unknowns.
Thus, a square tridiagonal matrix with a dimension equal to the number of
distance intervals in the solution reach must be solved for each time point.
It will be shown that the solution to the set of simultaneous linear equa-
tions is very simply accomplished with the digital computer.
To facilitate the following discussion, equations will be written in
the FORTRAN language. For example, the value for BOD concentration (mg/1)
at the I'th distance point from the upstream boundry of the solution reach
and J'th time point from time=0 will be written as L(I,J). Also, the
asterisk (*) will indicate multiplication and the slash (/) will indicate
division. Since in equations (23) and (24) the dispersion coefficient
always appears divided by delta X squared we will redefine E /AX2 and E./AX2
as EC and ED. c a
The time increment (days) will be called DT and DT/2 will be called
DT2. Similarly, the ratio of increments on the left-hand side of equation
(23) will be called DLDT(I.J) and the left-hand side of equation (24) will
be called DDDT(I,J). Since the solution is built up one row at a time pro-
gressing towards greater values of time, the following two equations will
apply:
.J) + DLDT(I, J+1))/2)*DT (25)
D(I,J+1) = D(I,J) + ((DDDT(I.J) + DDDT(I, J+1))/2)*DT (26)
Combining equations (23) and (25) it can be seen that a linear equation of
the following form will result for each point within the solution reach:
A*Ld-l, J+D + B*L(I, J+l) + C*L(I+1, J+l) = D (27)
The constants A, B, C, and D can all be found from known values of L
(I,J) along the J'th time line, from the characteristics of the incoming
pollutional pulse, or from the assumed coefficients. L(l, J) is the value
of BOD along the vertical left-hand boundry of the solution plane and is
supplied as input to the program.
16
-------
It can be seen from Figure Ic that as the
time row to the next the solution contains one
a new line is computed. Therefore, if the solution
tion reach containing NX distance intervals for all
increments the number of distance points must be MT
number of distance points must then be reduced by one for
computed to insure that only valid points on the (J-l)'th
compute points on the J'th line. The number of distance points which enter
into the solution is NINC+1 which has an initial value of MT+NX+1. The
values for L and D at the NINC+2 interval are assumed to be zero. This is
shown diagrammatically in Figure Ic.
solution progresses from one
less valid point each time
is required for a solu-
times out to MT time
+ NX + 1 initially. The
each time line
line are used to
The index (I) in
Figure 1. Values for
as follows:
equation (27) ranges from (2)
the constants A-D in equation
to NINC+1 as shown in
(27) are calculated
A = -Q*DT2/V - EC*DT2
B = Q*DT2/V + 1 + EC*DT
C = -EC*DT2
Q*DT2*(L(I-1,J)-L(I,J))/V
+ DT2*(QIN(I,J)*LIN(I,J)
+ EC*DT2(L(I-1,J) - 2*L(I,J)
- KR*DT2*L(I,J)
(28)
KR*DT2
When (I) has a value of (2) an additional term of -A*L(1,J+1) should
be added to the expression for D(I,J) if the initial concentration of BOD
in the stream upstream of the pollutional pulse is not zero. The equations
for A, B, C, and D as given apply to the Euler coordinate system. Cor-
responding equations for the Lagrange system are obtained by simply deleting
the first term in the expressions for A, B, and D.
Equations (24) and (26) can similarly be combined to form one linear
equation for each point along the X axis all with the same form as equation
(27) except that the D(I,J) is substituted for L(I,J). Relationships for
A, B, and D which apply to dissolved oxygen deficit (D) are given below:
A = -Q*DT2/V - ED*DT2
B = Q*DT2/V + 1 + 2*ED*DT2 + KA*DT2
C - -ED*DT2
17
-------
= Q*DT2*(D(I-1,J) - D(I,J)/V
+ ED*DT2*(D(I-1,J) - 2*D(I,J)
+ D(I,J)*(1 - KA*DT2)
+ KD*DT2*(L(I,J) + L(I
If the dissolved oxygen deficit (D) is not zero along the left-hand vertical
boundry of the X-T plane, an additional term should be added to the expres-
sion for D(I,J) when 1=2. This term is -A*D(1,J+1).
The coefficients as they appear above apply to the Euler coordinate
system. The first term of the equations for A, B, and D should be deleted
when the Lagrange coordinate system is used.
The number of simultaneous linear equations to be solved to advance
the solution one time increment is NINC. Each equation has three terms
centered along the diagonal of the matrix except the first and last which
have two. This tridiagonal of the matrix can be solved by Gauss Elimination
but a more direct and convenient method is given by Bunce and Hetling (6).
To apply this method, first, compute two vectors R(N) and S(N) as follows:
S(2) = B (30)
S(N) = B - C*A/S(N-1) (31)
In equation (31) the index ranges from 3 to NINC+1. The vector R(N) is com-
puted as follows:
R(2) = L(2,J) (32)
R(N) = L(2,N) - R(N-1)*A/S(N-1) (33)
In equation (33) the index (N) ranges from 3 to NINC+1. Notice that the
vector S(N) must be computed first in order that it can be used in generating
vector R(N). The values for L on the (J+l)'th row are then computed as
follows beginning with a value of NINC for N and working backward down the
numerical scale until N=2.
L(N,J+1) = (R(N) - C*L(N+1,J+1))/S(N) (34)
Using equation (34) values for BOD on the J+l row are calculated between an
index of 2 and an index of NINC. The BOD value in the last (NINC+1) distance
point is computed as follows:
L(NINC+1, J+l) = R(NINC+1)/S(NINC+1) (35)
18
-------
To solve for the values of dissolved oxygen deficit on the (J+l)'th
row the vectors S(N) are computed in a completely analogous way. The value
of D(NINC+1,J+1) is found as follows:
D(NINC+1,J+1) = R(NINC+1)/S(NINC+1) (36)
Letting the index (N) take on all values between NINC and 2 working backward
down the numerical scale, the remaining values for dissolved oxygen deficit
are computed as follows:
D(N,J+1) = (R(N)-C*D(N+1,J+1))/S(N) (37)
This method of solution is also discussed on page 136 of the book by
R. V. Thomann (3). The report (6) by R. E. Bunch and Leo J. Hetling entitled
"A Steady State Segmented Estuary Model" gives a particularly clear presenta-
tion of the method. A diagram of the stream divided into segments is shown
in Figure la. The solution reach is between segment 2 and segment NINC+1.
Initial conditions in segment (1) are assumed to be fixed for all time and
both the BOD and dissolved oxygen deficit are assumed to be zero in segment
NINC+2. The pollutional pulse is shown entering segment (3). This segment
number is called POINT in the computer program and it must be far enough
downstream so that the dispersion will not carry the pollutant upstream
beyond the second segment.
A number of computer programs have been developed based on the process
relationships and the computational scheme described here. Listings of
these programs and instructions for their use are given in references 7-11.
ACCURACY OF NUMERICAL INTEGRATION SOLUTION
The computational accuracy of these programs was compared to the
Streeter-Phelps solution and to the closed-form solution for an impulsive
load. It was found that the use of Lagrange coordinates allowed much greater
accuracy of computation than that possible with Euler coordinates. When
the Lagrange coordinates are used the distance interval must equal the
time interval times the stream velocity. Using the Lagrange coordinates,
it was found that a 15 minute time interval gave sufficient accuracy when
the duration of the square pollutional pulse was one hour. Similarly, it
was found that when the duration of the pollutional pulse was 24 hours or
greater the solution approximated the Streeter-Phelps solution. A one-hour
time interval gives sufficient accuracy when the duration of the pulse is
24 hours.
For example, Figure (3) shows computed points connected by a dashed
line for a 24 hour pulse entering a stream with a velocity of 13 miles
per day. The corresponding Streeter-Phelps solution is shown by the solid
line. It can be seen that the accuracy is adequate. The corresponding
values for dissolved oxygen deficit are shown in Figure (4) which also shows
that the solution for a 24 hour pulse agrees well with the Streeter-Phelps
solution.
19
-------
7 -
6
-T+ — ±± irrhtt-r!: ±t~r:o:~r~ "'-^fftt--'--^- - +
:4g^:: |E|£Ed ;|::|"::$:::Wi"::g::|
T1 — THT -1-F^-TTJ-TJ-'""'''^'^'"""^"^'""''''''^'
::|±#|||: :|:::|:::::::I:::|::::::|
•-"kb~+ ~ i __i:t-i Ti--i — 4
±:3:fa t4ffi Vill 1 1 1 HI 1 1 1 1-^rH^^Frt
::i:!e*:E=->E£--]=g-!i-L streeter-Phelps 5
:EE:,::: EF| ±| :::|___:.± -j-.---.-i
:::::|:|::|f::: :::|:||;||:E|EEEEEE|E|EEE
:::::±-S "i I i^ti--:::::::::]
FaBSBiBBft 1 1KB
L: '111 ! ill
- o H Jp± ic ::::::--:::::::::
;^::::;:E^p:|::-:::;::::^;i:-;::::::::
- Oj 'r T
- ^ ---+ T T —
j-j M i M Ml r> ' -'- N44--1 ' 4
- c --T (^ ::::--:
: g j|p:: ±j|,::::::g :::::::::::::::::::::::
;g:|:|:::i:p^:::::::::E:::::::::::::|::E
: ° ::::::::: ::::::::?•+::: --S"= Numerica
n o f ffu. L L-jjJ |J L yj-l ] \\\\+ • I
;^-EE~EpE|EEE:||-~EiE~EEEEE~EEE
:::::::p:::::::±::::::|::::::±::::i:::^::Sj--:
p .
i---P i 5-;
::Srl:::|:ffi|:::::::::::J::::::"::|:l!:2»
::::::::: :::±:r:: ::: + ::::: ::::::::::;:::;:::^:5:
|:::::|:|| : :I ::::::::::: ::::;|:::::::::|::
1.5 days —::-::!::
::::::::::::::::::::::::::::::: V = 13.1 ml/ds
:::::::::::::::|:::::::::::::: L = 6.25 mg/]
::::|:~|::i::::::::::::g: K = 1.0 day a
::|:|:::::g:::;|::::::i:: ^ = 0.23 day
::::::::::|:::l::|:::::::g K - 0.4 day -1
. a.
:::::i:E:::;::^::|::::::: E = 0.1 sq mi/d<
Solution :::::::::: + :::: :::::::::::::: :± ::::::
1 Integration Solution + ^
1 1 j j i M j I n 1 1 1 1 II i M 1 1 j j 1 1 1 1 1 1 1 1 1 M 1 1 1 :::::::::::::::::: :±
::::| :::::::::: :::|::::: :::;::::::::::::::::::::;::::::::
-__. ± ± r± ...±
na__
:EEEEEEEEEEEEEEE^:|j;EEEEEEEE|;EEEEEEEEiEE?!E||;EE
*ir- i L*Ua''j:e '"
-i 2.5 days --s :.:•=. s:;;:-:::::::::::::
j_._. 3.5 days
LY H4:
i [jj|j
iy::::
E|!=i =
0 4
12 16 20 24 28 32 36 40 44
Distance downstream from outfall, miles
48
Figure 3. Computed BOD profiles in the stream caused by 24 hour storm
overflow event
20
-------
1.0
.9
.7
.6
.5
.4
.3
.2
.1
.0
I" V =
J L =
o
J1 r K
• H "I
O t~~ V
LH tE d
:d a
13
6
1
0
0
,1 mi/day
.25 mg/1
0 day'1 Ht
23 day
-1
4 day
-1 :::
E = 0.1 sq mi/day
\
Streeter-Phelps Solution
Numerical Integration Solution!:
12 16 20 24 28 32 36 40 44
Distance downstream from outfall, miles
48 52
Figure 4. Dissolved oxygen deficit profiles in the stream caused by 24 hour storm
overflow event
-------
The computer program described in reference (4) was used to find BOD
profiles in the stream caused by a square pollutional (BOD) pulse of one
hour duration. The one hour duration was broken up into ten increments and
each increment was taken as an impulsive load for which a closed-form solu-
tion is available. The solution for the one hour pulse was then taken as
the sum of the solutions for the ten increments. BOD profiles in the stream
caused by the one hour BOD pulse are shown by the solid lines in Figure (5)
which have the shape of the normal probability density function. The
computed points are circled in Figure (5). The corresponding Streeter-Phelps
solutions are shown by the vertical bars with dashed lines. It can be seen
that the computed points fit the closed-form solution reasonably well. The
time interval used was 15 minutes. The corresponding profiles for dissolved
oxygen deficit are shown in Figure (6). The solid lines represent the
numerical integration solution and the bars enclosed by dashed lines repre-
sent the corresponding Streeter-Phelps solution. It can be seen that
dispersion can have a significant effect on the shape of the BOD and dissolved
oxygen deficit profiles caused by storm overflows.
SUMMARIZING THE IMPACT OF SQUARE POLLUTION PULSES
The digital computer program described in this report can be used to
precisely compute the BOD and dissolved oxygen deficit (DO) profiles in the
receiving stream caused by any defined storm overflow event.
Because of other uncertainties associated with the storm overflow
problem it is seldom of interest to compute the stream impact with great
precision. Therefore, it is clearly desirable to present the computed
results in terms of graphical relationships which can be easily used to
find a rough solution to any problem.
In general, when a square pollutant pulse is introduced into the stream
a short segment of the stream is polluted and this segment is carried down-
stream at the velocity of the tubular stream. The initial length of the
polluted segment is simply the duration of the square pollutant pulse times
the stream velocity. When dispersion is assumed to be zero the length of
the polluted segment is fixed and the Streeter-Phelps solution can be used.
If dispersion is introduced into the problem the length of the polluted
segment will increase with time and the shape of the BOD and DO profiles
(at any time) will have a shape which approaches the normal error function
as time becomes great. The shorter the width of the pollution pulse and the
greater the value of the dispersion coefficient the less time is required
for the shape of the DO profile to reach the normal curve shape. The
standard deviation of both the BOD and DO profiles will equal (2Et)% where
t is the time from entry of the initial portion of the pulse. The standard
deviations of the profiles from a one hour square pulse with values of
dispersion coefficient (E) of 0.1, 0.2, and 1.0 sq mi/day are shown in
Figure (7). It can be seen that the profiles are normal after about 150
15 minute time intervals or 1.56 days. Since the maximum peak occurred at
1-1.2 days it can be seen that the profiles did not become normal until after
the maximum peak. The effect of dispersion is minimized when the initial
22
-------
ro
GO
LI
v = 13.1 mi/day
-1.4 mg/1
LQ = 6.25 mg/1
K = 1.0 day
r -ir
= 0.23 day
K = 0.4 day"1
E = 0.1 sq ml/day-
2.5 days J
21 22 30 31
17 18
33
34 35 43 44 45
46
47 48
Distance downstream from outfall, miles
Figure 5. Computed BOD profiles in the stream caused by one hour storm overflow event. Solid line
= closed-form solution, dashed line = Streeter-Phelps solution, circled points = numerical
solution
-------
rv
txO
C
o
•H
4->
rt
f-i
4-1
c
O
CJ
C
O
O
•H
O
•H
C
CD
fc>0
X
X
o
I/)
•H
Q
0.8
0.7
0.6
0.5
0.4
0.5
0.2
0.1-
v = 13.1 mi/day:
;LQ = 6.25 mg/1 E:
;K = 1.0 day'1
:E = 0.1 sq mi/day:;
= 0.23 day
K = 0.4 day
a
-1
3.5 days:
1;
r
47 48
17 18 19 20 21 22 30 31 32 33 34 35 43 44 45 46
Distance downstream from outfall, miles
Figure 6. Dissolved oxygen deficit profiles in the stream caused by one hour storm overflow event.
Solid line = numerical solution dashed line - Streeter-Phelps solution
-------
to
1—
fd
s_
OJ
-M
-M
(T3
fl3
•o
c
GO
5 6 7 8 9 10
10
Figure 7.
100 1000
Number of 15 min time intervals after entry
Computed standard deviations for dissolved oxygen deficit
profiles at time from entry of leading edge of pulse for
dispersion coefficient values of 0.1, 0.2, and 1.0 sq mi/day,
-------
length of the polluted segment is large or the dispersion coefficient is
small. The BOD profile will be square initially and the DO profile will be
zero initially.
The peak DO values within the polluted segment will increase from zero
to a maximum (DOMAX) and then descrease to zero. It was found that the
ratio of the maximum peak to the Streeter-Phelps peak DO value could be
plotted as a single valued function of the following non-dimensional grouping:
E = dispersion coefficient, sq miles/day
W = pulse width in the stream (duration x
stream velocity), miles
KA = reaeration rate coefficient, days"
KD = deoxygenation rate coefficient, days"
The relationship between this non-dimensional grouping and the ratio of
DOMAX to the Streeter-Phelps DOMAX is shown in Figure (8).
The time from the introduction of the initial part of the square
pollutant pulse to the point where DOMAX occurs is also of interest. It
was found that if this time (t) in days is expressed in terms of the
following non-dimensional grouping it will plot as a single valued function
of the first non-dimensional grouping:
Et/W2
E = dispersion coefficient, sq miles/day
t = time from initial introduction of pollutant
to time at which DOMAX occurs, days
W = initial pulse width, miles
A plot of this second non-dimensional grouping as a function of the first
is shown in Figure (9). Computed results used to plot Figures 8*and 9 are
shown in Table 1 .
In summary, it has been shown that by the use of non-dimensional
relationships a reasonably complete picture of the impact on the dissolved
oxygen resources of the stream caused by a storm overflow can be estimated.
However, if a closed-form expression for the time-distance profiles of
dissolved oxygen deficit caused by an impulsive or spike load could be
found a much more straight-forward and satisfactory method for estimating
the impact of storm overflows would be possible.
26
-------
1.0
0.1
9
8
7
6
6
4
3
2
10
9
8
7
6
5
4
3
2
1
).
— P-
-I*
F=^
X
o
, Q
' Pi
i X!
' a
. 0>
. 0
" X
: Q
-
H-
— 1
-®£"'-2sI
— H — r
1
~i it
— H-: = ::qp*
I
+ IIP11
2
1
^+-TT
3r— T
S»
3 ^
r '^
_Li— i — L_ -_<4:
::;E|EE|E|
r
T
T
+fn
T
1
.1.
r
•A
tii ^
i
—
i.
E
E - dispersion coefficient, sq mi
W = pulse width in stream, miles
KA = reaeration coefficient, days
KD = deoxygenation coefficient, da
/da
1
ys
y
i
DOMAX = peak dissolved oxygen deficit, mg/1
~rri l i \ i i i i r rvi i \ 1 1 1 1 rr i i i i fi 1 1 M r n \ i > i i i i i \ i i i i 'i M i i i i M t
N
*
i
i
\
^
T
-
N
--
VW/(KA*KD)
™
?
~]~
i
T
T,
4
1 i
^
:i
ON
I
k
' L
3
r-U -
I
^
<;
==
--
r
— i
b
—
\
~~
—
~. i
~~:
|
r
fru
o
^j.
*
3 45678910 2 3 45678910 2
1.0 10.0
0.01
Figure 8. Correction factor for peak dissolved oxygen
deficit to be applied to the value computed
with the Streeter-Phelps model
27
-------
3 4 567891
0.01 ,
0.1
1.0
10
100
Figure 9 . Time (t, days) between entry of leading edge of square
stormwater pulse and time of peak DO deficit in terms of
non-dimensional groupings where: E=dispersion coef,
sq mi/day W=pulse width in stream, miles KA=reaeration
coef, I/day KD=deoxygenation coef, I/day.
-------
Table 1
Computed Results
Velocity, mi/day
Pulse width (W) , mi
I K
V V
Kr Kd
E
SP-t*
SP-DOMAX/LO
n K
V ~- V
Kr - ^d
E
SP-t*
SP-DOMAX/LO
in v
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
IV K
Ka
Kr = Kd
E
SP-t*
SP-DOMAX/LO
V
Ka
Kr = Kd
. E
SP-t*
SP-DOMAX/LO
vi K
Kr = Kd
E
SP-t*
SP-DOMAX/LO
VII K
Kr = Kd
E
SP-t*
SP-DOMAX/LO
VIII
i\a
K\r
y» l\^4
E
SP-t*
SP-DOMAX/LO
IX .,
Ka
K — V
r ~ M
E
SP-t*
SP-DOMAX/LO
4.68
0.7
= 0.1
= 0.4774
0.1071
4.68
0.7
0.2
0.4774
0.1071
4.68
0.7
0.3
0.4774
0.1071
-,4.68
0.7
1.0
0.4774
0.1071
2.0
0.7
0.1
0.8076
0.1989
2.0
0.7
1.0
- 0.8076
0.1989
0.4
0.2
1.0
3.466
0.25
0.4
0.4
1.0
2.5
0.368
0.4
0.7
1.0
- 1.865
0.474
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
2.4
0.1
= 0.229
= 0.01658
2.351
2.99
0.229
0.01178
3.324
5.98
0.229
0.00963
4.071
8.97
0.229
0.00528
7.433
29.9
0.406
0.02316
2.907
4.06
0.354
0.00722
9.193
35.40
0.844
0.00344
18.803
84,40
0.719
0.00637
15.811
71.90
0.604
0,01014
13.747
60.40
6.0
0.25
.250
.03965
.940
.400
.240
.02878
1.330
.768
.240
.02371
1.628
1.152
.229
.01316
2,973
3.664
.427
.05647
1.163
.683
.406
.01835
3.677
6.469
1.635
.01104
7.521
26.160
1.210
.01910
6.324
20.160
.937
,02860
5.499
14.992
12.0
0.5
.312
.06911
.470
.125
.271
.05339
.665
.217
.260
.04503
.814
.312
.240
.02588
1.487
.960
.490
.10405
.581
.196
.417
.03636
1.839
1.668
1.729
.02213
3.761
6.916
1.281
.03811
3.162
5.124
.948
,05698
2.749
3,792
24.0
1 .0
.417
.09602
.235
.042
.354
.08435
.332
.071
.323
.07569
.407
.097
.260
.04868
.743
.260
.635
.16146
.291
.635
.448
.07021
.919
.448
1.760
,04389
1.880
1.760
1.312
.07536
1.581
1.312
,979
,11222
1.375
.979
32.15
1 .34
.458
.10201
.175
.026
.406
.09452
.248
.045
.375
.08781
.304
.063
.292
.06146
.555
.163
.708
.17908
.217
.039
.469
.09081
.686
.261
1,792
.05827
1.403
.998
1.333
.09976
1,180
.742
1,010
,14796
1.026
.563
48.0
2.0
.490
.10569
.118
.012
.458
.10276
.166
.023
.437
.09940
.204
.033
.344
.07970
.372
.086
.781
.19247
.145
.020
.531
.12407
.460
.133
1.865
.08496
.940
.466
1.406
.14428
.791
.352
1,083
,21182
.687
.271
29
-------
The FORTRAN listing, input and output definitions, and sample
printouts for the SWOPS computer program are given in the Appendix
30
-------
REFERENCES FOR SECTION 4
1. Streeter, H. W. and Phelps, E. B. "A Study of the Pollution and
Natural Purification of the Ohio River III", Public Health Bulletin
No. 146, 1925.
2. Dobbins, William E. "BOD and Oxygen Relationships in Streams",
Jour. Sanitary Engineering Div.s ASCE SA3, June, 1964 pp 53-77.
3. Thomann, Robert V. "Systems Analysis and Water Quality Management",
McGraw-Hill Book Co., 1972.
4. Eilers, Richard G. "Stream Model (STREAM4): Closed-Form Solution
for Simulation of BOD in Streams", Internal Memo, May 17, 1977.
5. Smith, G. D. "Numerical Solution of Partial Differential Equations",
Oxford University Press, New York, 1965.
6. Bunce, Ronald, E. and Hetling, L. J. "A Steady State Segmented
Estuary Model", Tech. Paper No. 11, FWPCA, U.S. Dept. Interior.
7. Eilers, Richard G. "Stream Model (STREAM1): Dresnack-Dobbins
Numerical Analysis of BOD and DO profiles", Internal Memo, May 17,
1977.
8. Eilers, Richard G. "Stream Model (STREAM2): Closed-Form Solution
for Stream Simulation Using Euler Coordinates", Internal Memo, May 17,
1977.
9. Eilers, Richard G. "Stream Model (STREAMS): Crank-Nicolson Implicit
Method Solution Using Euler Coordinates", Internal Memo, May 17, 1977.
10. Eilers, Richard G. "Stream Model (STREAMS): Crank-Nicolson Implicit
Method Solution Using Lagrange Coordinates with a Single Influent
Point", Internal Memo, May 17, 1977.
11. Eilers, Richard G. "Stream Model (STREAM6): Crank-Nicolson Implicit
Method Solution Using Lagrange Coordinates with Multiple Influent
Points", Internal Memo, May 17, 1977.
31
-------
APPENDIX FOR SECTION 4
Table 2
Definitions of FORTRAN Symbols Used in SWOPS
Input Variables
NCASE number of data cases to be calculated by the program.
LIST alpha-numeric identification for each data case.
KR rate constant for removal of BOD by deoxygenation and
sedimentation, days" .
KN rate constant for removal of ammonia nitrogen by
deoxygenation, days' .
KA rate constant for reaeration from the atmosphere, days' .
KD rate constant for removal of BOD by deoxygenation, days' .
DX distance interval for the numerical integration
procedure, miles.
DT time interval for the numerical integration procedure, days
XN number of distance increments in the solution reach.
TM number of time increments to be calculated.
EC dispersion coefficient for BOD, sq mi/day.
EN dispersion coefficient for ammonia nitrogen, sq mi/day.
ED dispersion coefficient for dissolved oxygen, sq mi/day.
VEL velocity of the stream, miles/day.
BEN benthic oxygen demand, mg/l/day.
32
-------
Table 2 (continued)
POINT program control: distance point (number of increments from
the first upstream point of the solution reach) at which
the storm overflow enters the stream.
Q total stream flow, cfs.
ZM program control: number of time increments in the storm
hydrograph.
DAY program control: time point at which the program begins
printing output, days.
Y1-Y5 set concentrations for the BOD or DO increment calculations,
mg/1.
DL program control: 0 = BOD calculation, 1 = DO calculation.
PRNT program control: Q = increment calculation printout,
1 = selected concentrations printout.
XI distance increment at which program printout begins.
X2 distance increment at which program printout stops
(note that X2-X1 must equal 20).
QNT(M) volume of the storm overflow stream at time increment M, cfs
LOT(M) BOD concentration of the storm overflow stream, mg/1.
LNT(M) ammonia nitrogen concentration of the storm overflow
stream, mg/1.
Output Parameters
N distance increment along the stream.
DIST distance along the stream at increment N, miles.
TT time of travel downstream, days.
LO(N) carbonaceous oxygen demand (5-day BOD) at N, mg/1.
LN(N) nitrogen ultimate oxygen demand (exclusive of LO) at N,
mg/1.
DO(N) dissolved oxygen deficit at N, mg/1.
33
-------
Table 2 (continued)
M time increment of calculation.
IT time interval at which the maximum BOD or DO occurs.
TTRAV time of travel (I*DT), days.
DOMAX maximum DO, mg/1.
LOMAX maximum LO, mg/1.
SUM1-SUM5 sum of DO or BOD at time increment M based on Y1-Y5 inputs,
IM1-IM5 number of increments exceeding the Y1-Y5 inputs at time
increment M.
DOT(I) if DL=0 - DO deficit in the stream at distance interval I/
initial BOD concentration in the stream.
if DL=1 - BOD concentration in the stream at distance
interval I/initial BOD concentration in the stream.
34
-------
Table 3
FORTRAN Source Listing for the SWOPS Program
C (SWOPSl
C CRANK-MTCOLSON IMPLICIT METHOD SOLUTION (LAGRANGE COORDINATES)
C SINGLE INFLUENT POINT
C ROD OR DO SUMS APE CALCULATED
C R, G. FILERS SEPTEMBER 1P77
PEAL KR,KN,KA,KD,LOf700) , LN ( 700 ) , LOS ( 700 ) , LOT (30) ,LNT(30) ,
. LOTA,LOMAX
DIMENSION LISTf40),DO(700) , QNT ( 30 ) ,DE ( 700 ) ,R(700),S(700),
. DOT(700),NA(21)
OPEN( UN IT*lfNAMEB» SWOPS, DAT' , TYPE a 'OLD' , FOP MB » FORMATTED ' ,
. READONLY)
TNal
10*6
PEAD(IN,10) NCASE
10 FORMATd?)
DO 1000 TITsI, NCASE
DO 20 Nsl ,700
LOfN)*0,
LNfN)aO.
DO(N)aO.
LOSCN)aO.
DE(N)aO.
R(N)sO.
S(N)aO.
DOT(N)aO.
20 CONTINUE
DO 30 Mai ,30
QNT(M)aO.
LOT(M)aO.
LNTfM)aO.
30 CONTINUE
READ(IN,40) LIST
40 EORMAT(40A2)
READ(IN,50) KR,KN,KA,KD,DXeDT,XN,TM
RE AD (IN, 50) EC,EN,ED,VEL,BEN,POINT,Q,ZM
READ(IN,50) DAY,Y1,Y2,Y3SY4,Y5»DL,PRNT
PEADfTN,50) X1,X2
50 FOPMAT(8E10,0)
NX»XN
MT»TM
IP«POINT+1 .
N1aX1
N?«X?
NTsN1
35
-------
Table 3 (continued)
DO 55 Nalf21
NA(N)sNT+N
55 CONTINUE
PEAr>(IN,6Q) rQNTfM) ,M«1,30)
REftDf TNjfiQ) fLQTfM),Ms1 ,30)
PEADfIN,60) (LNTCM),Ms1 ,301
60 FORMATf 10F8.0)
WRTTE(IO,70) LIST
70 FORMATdH! ,/////, 20X,40A2,//)
WRITE (10, 80) KR,KN,KA,KDtOX»DT,XN,TM
80 FOR" AT CSX, ' KR • , 1 OX „ • KN ' , 1 OX , ' K A ' , 1 OX , ' KD ' , 1 OX , ' DX • , \ OX , ' DT ' ,
. 10X, 'XN' ,1 OX, *TM' ,/,8F12.4,/)
WRTTFC 10,90) EC,FN,ED,VEL,PEN,PQINT,O,ZM
90 FORMS TC8X, ' EC ' , 1 OX , ' EN ' , 1 OX , ' FD ' , 9X , ' VEL ' , 9X , 'REN' ,7X, 'POINT' ,
. 11X, '0' ,10X, 'ZW ,/,fiF12.4,/)
WRITEf 10,95) DAY,Y1 , Y2 , Y3 , Y4 , Y5 , DL ,PRNT
95 FORMATf7X, 'nAY',10X,'Yt',1pX,'Y2',10X,'Y3',lOX,'Y4',lOX,'Y5',
. 10X,'DL',fiX,'PRNT«,/,8F12'. 4,/)
WRITFdO,97) XI, X?
97 FORMA Tf flX, (Xl',10X,«X?',/,?F12,4,/)
WPITEf 10, 100) fQNT(M) ,M«1 ,30)
100 FORM AT (7X, 'ONT',/,t OFi9t4,/f10Fl2.4,/,lOFl2.4,/)
WRTTFdO, 110) ftOTfM),v»l f 10)
1 10 FOPMATC7X, 'LOT'f/,tOF12,4f/,10F12.4f/f10F12.4,/)
WRITEdO, 120) CLNTfM),Msl ,30)
120 FORMAT (7X, 'LNT'f/.10F12.4,/,10Fl?.4,/,10F12.4,////)
IF fPRNT) 128,121,128
121 IF (DL) 124,122,124
122 WRITE(IO, 123)
123 FORMAT f 8X, 'M',9X,'T',5X, 'TTRAV ,5X, 'LOW AX' ,6X, 'SUM ' ,6X, 'SUM2' ,
. 6X, 'SUM 3' ,6X, 'SUM4' ,6X, 'SUM?' ,9X,'IM1',3X,'I.M2»,3X,'TM3»,3X,
, 'IM4' ,3X, 'IM5« ,/)
GO TO 12«
124 WPITEf 10,125)
125 FORM AT (8X, 'M',9X,'I',5X, iTTPAV' ,5X,'DOMAX',6X, 'SUVl ' ,6X, 'SUM2> ,
. 6X, 'SUM?' ,6X, 'SUM4' ,6X, ' SIJM5 ' ,9X,'IMl',3X.'lM2',3X»'IM3',3X
GO TO 137
128 IF (DL) 134,130,134
130 WRTTEfIO,132) (N A ( II ) , I lal , 2 11
132 FOPMftTfl2X, 'BOD IN STREAM / INITIAL BOD IN STREAM',//,
. 4X, 'M' ,21(2X,I3,1X),/)
GO TO 137
134 WRTTFf 10,136) (M A ( T I ) , I I«l , 2 1 )
13ft FORMATM2X, 'DO DEFICIT IN STREAM « 100 / INITIAL BOD ',
. 'IN STREAM', //,4X. 'Mi, 21(2X,T3,lin,/)
137 DT*DT/24.
DXaVEL«DT
36
-------
Table 3 (continued)
TOAYaDAY/DT+l
FN«FJN/DX##2
FDaFD/DX**2
NINCBMT+NX+2
DT?sDT/2.
TlsO
I2aO
T3aO
T4sO
I5aO
SUMlaO.
SUM2»0.
SUM4rO.
Tf n«QNTt
DO 500 Mel,MT
D 0 M a X x 0 .
TTPAVsO.
TTaO
TRAV1 BO.
TT1*0
TPAV?sO.
1
NINCsNINC
As-EC«DT?
Ro1 ,+Er#DT2*2.+KP*DT2
LOf 1 )»0.
IPaIP-1
K«NINC+1
DO ISO N«2,K
LOS(N)aLOfN)
IF (N-IP> 140,138,140
139 IF (M-MZ1 139,119,140
139 ONTA«ONTfM)
LOTAainTfM)
GO TO 145
140 QNTAaO.
LOTA*0.
1 45 DEf Nl
. LOfN)*fl .-
150 CONTINUE
P(2')»*DEC2)
S(2)aB
no 160 Ne3,K
-2 . *LO f N 5 +LO f N-f 1 5 )*HT2*
37
-------
Table 3 (continued)
160 CONTINUE
DO 170 N*3,K
R(N)sDF, fM)-R(N-n*A/S(N-1
170 CONTINUE
KsNINC-1
DO IfiO Nsl ,K
I»NINC+1-N
LO(T)"fPf I)
IF (LOMAX-L
175 LOMAXaLOflf
TRSVlaI*DT
ITlsI
180 CONTINUE
ron )»o.
175,17SilSO
Pal ,+2.*En#nT24KA*DT2
DO 190 N«2,K
) )*Fn*DT24-DO C N ) * f 1 ,-KA*DT2)*
190 CONTINUE
R(2)sDFf2)
DO 200
?00 CONTINUE
DO 210 N»3,K
210 CONTINUE
KsNTNC-1
DO 220 Nsl.K
TsNTNC»N4t
IF
215,215,220
TPAV?sI#DT
IT2*I
220 CONTINUE
IF (DM 924,222,224
222 TTPAVsTRAVl
TTsTTt
38
-------
Table 3 (continued)
GO TO 276
724 TTPAVsTRAV?
TT»IT?
226 TMlsO
IM?SBO
IV5sO
no 370 Nsl,NlNC
IF (DM 230,228,250
728 VAR«l,0(N)
GO TO 21?
230 VARsDOfN)
73? TF fVAP-Yl) ?40,??5,?15
235 IM1«TW1+1
240 TF (VAR»Y2^ 260,250,250
250 TM?sTM24-1
260 IF (VAP-Y3) 780,270,270
270 TM3slM3f1
280 IF fVAP-Y4) 300,290,290
790 IM4=TM4+1
300 IF (VAR-Y5) 320,310,310
310 IM5sTM5+1
320 CONTINUE
iisn +TMI
I2BT24.IM?
T3sI3+IM^
T4sI4*TM4
SUM4sT4#nX*DT
SUM5aT5*nX*DT
IF (PPNTT 340,370,340
340 DO 350 N»N1 ,N2
IF fflM 345,342,345
34? HOTf N)=LOfN)/BnD
GO TO 350
345 DOTfN)«nO(N)/POD#lOO.
350 CONTINUE
360 FORMftTf I5(,21F6.31
GO TO 500
370 IF (DM 390,380,390
380 WRTTEf 10,400) M , IT , TTR AV , LOW AX , SUM 1 , SUM2 , 5UM3 , SUM4 , SUM5 ,
, IM1 ,TM2, IM?
GO TO 405
39
-------
Table 3 (continued)
390 WRTTF. (10 s 4001 "• , TT , TTRA V , DOM AX , SUM 1 , SUM 2 , SUM 3 , SUM 4 , SUMS ,
, TM! f\W1, TM3, TM4, TM5
400 FORMAT (5X, T5,5X,I5»7F10,5,5X»5I61
40? Ttr (M-IDAY) 500f410»"iOO
410 WRTTF(IO»4?0) M,DAY
420 FOPMATr///,5X, 'M a ' , T5 f 5X , ' DAY S',F7. ?,///)
WRTTFf TO, 4301
DO 450 Nsl ?NX
TTB(N*POTNT1#DT+DAY
WPTTFf TO, 4401 N , n I ST , TT „ LO f N 1 , L« f N ) , DO (N )
440 FOPMAT(6X,I6,5F12e5)
450 CHMTTNUE
500 rOMTTNUE
TF fPRNT) lOOOsBOS^ 1000
505 WRITF(Tnf 51 0) T 1 , T 2 , T .1 , I 4 , TS
510 FOPMAT(///
1000 fONTTNUF
FND
40
-------
Table 4
Sample Printout from SWOPS (DL = 0, PRNT = 1)
. OXYGFN DEMAND - DISTANCE INCREMENTS 40 TO 60
9-71-77
K"
0.2500
FC
o. 1000
niY
7S.oooo
xi
40.0000
OUT
1 294 .ooon
0.0000
n . ooon
LOT
Ho. ooon
0.0000
o.oooo
LNT
0.0000
0.0000
o.oooo
KN
0. 0000
FN
0. 1 000
Y1
6.0000
X7
60.0000
1704 .ooon
o . oono
o.ooon
1 1 n.nnno
o.oooo
o.oooo
o.oooo
o . oonn
o.oooo
K6
0.9800
Fn
O.inno
Y?
4 .0000
1294.0000
0.0000
o.nooo
1 1 o.ooon
0.0000
0.0000
0.0000
o. noon
o.oooo
n
7.4
7
1794
0
0
1 10
0
0
0
0
0
KD
.7500
VFL
.0000
Y3
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
DX
0.0000
RFN
0.0000
Y4
1.0000
0 . 0000
0.0000
o.oooo
o.oooo
0.0000
o.oooo
0.0000
0.0000
0.0000
DT
0.7500
POINT
50.0000
Y5
0.5000
o.nooo
0.0000
o.oooo
0.0000
0.0000
0.0000
o.oooo
0.0000
0.0000
100
950
f,
0
0
n
0
0
0
0
0
0
XN
.0000
0
.0000
DL
.onoo
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.oono
.oono
TM
50.0000
ZM
4.0000
PPNT
1 .0000
0.0000
0.0000
0.0000
0.0000
o.nooo
o.oooo
0,0000
0.0000
0.0000
o.oooo
o.oooo
0.0000
0.0000
o.oooo
0.0000
0.0000
0.0000
0,0000
o.oooo
0,0000
0.0000
o.oooo
0.0000
0,0000
0.0000
0.0000
0.0000
POD IN STPFUM / INITT6L POD IN STREAM
41
42
45
46
47
4B
50
51
57
53
54
55
56
57
58
59
60
1
7
3
4
5
6
7
R
9
1 0
1 1
17
1 3
1 4
15
o .onn n .000
o.ooo n . ooo
0,000 0,000
o '. n o o o.nnn
0 .000 0.000
o .oon n. noo
0. 000 0.000
0.000 0 . 000
O.OOO 0.000
n.ooo r, ,noo
o.ooo o.ooo
o . noo o.ooo
0.000 0.000
o.ooo o.onn
0. ooo o .000
0.
0.
0.
o.
o.
o.
0.
o .
o.
n.
0,
0.
o.
o.
0.
ooo
000
000
ooo
ooo
ooo
ooo
000
ooo
000
000
000
000
orm
ooo
0
n
0
0
0
0
0
0
n
0
0
0
0
0
0
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
. 000
.000
.000
.000
.000
0
0
0
0
0
n
0
0
0
0
0
0
0
0
0
.000
.000
.000
.000
.000
.000
. 000
.000
.000
. onn
.000
.000
.001
.001
.001
0
0
0
0
0
0
0
0
0
0
n
0
0
0
0
.000
.000
.000
.000
.000
.001
.002
.003
.004
.005
.006
.OOR
.010
.011
.01 3
0
0
0
0
0
0
0
0
0
0
0
0
0
n
0
.000
.000
.000
.008
.024
.040
.054
.068
.080
.093
. 1 04
.115
.125
.135
.144
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
.000
.000
.008
.007
.988
.970
.953
.936
.970
.905
.B91
.877
.R64
.851
.839
r,
0
1
0
0
0
0
0
0
0
0
0
0
n
0
. 000
.no8
.007
.996
.993
990
.987
.994
.980
.976
.977
.968
.964
.960
.955
0
1
0
0
0
n
0
0
0
n
0
0
0
n
0
.008
.006
.994
.991
.987
.983
.979
.975
.971
.966
.967
.957
.953
,948
.943
0
0
0
0
0
0
0
0
0
0
0
0
0
r.
0
.982
.956
.939
.973
.90P
. 893
.879
.866
.R53
.840
.879
.817
.R06
.796
.786
0
0
0
0
0
0
0
0
0
0
n
0
0
0
0
.008
.074
.038
.052
.065
.078
,090
.101
.117
.172
.131
.140
. 149
.'57
.1*5
0.
0.
0 ,
0 .
0.
o.
0.
o.
o.
0.
0.
0 .
o.
0.
0.
ooo
000
001
007
002
004
005
006
008
009
01 1
013
015
017
019
0
0
0
0
0
0
n
0
0
0
0
0
0
0
0
.000
.000
,000
.000
.000
. 000
,000
.000
.000
.000
.001
.001
.001
.001
.001
0.
0.
0.
o .
o.
o.
0.
0.
0.
0.
0 ,
o.
o.
o .
0 .
ooo
000
000
000
000
000
000
000
000
000
000
ooo
000
000
000
o.
0.
o.
o.
0.
0.
o.
o.
0.
0.
o.
n.
0.
0.
0.
ooo
000
000
000
000
ooo
000
ooo
ooo
000
000
000
ooo
000
000
0.
0.
o.
0.
0.
o.
0.
0,
o.
0,
0.
o.
o.
0.
0.
ooo
000
000
000
000
000
000
ooo
000
000
000
000
000
000
000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
n
.000
.000
.000
.000
.000
,000
.000
.000
.000
,000
,000
.000
.000
.000
.000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,000
.000
.000
.000
.000
.000
,000
.000
.000
.000
.000
,000
.000
.000
.000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
.000
,000
,000
.000
.000
,000
.000
.000
.000
.000
.000
.000
,000
.000
.000
0,000
o.ooo
o.ooo
0,000
o.ooo
0 ,000
0,000
0.000
o.ooo
0,000
0.000
o.ooo
0.000
o.ooo
o.ooo
-------
ro
Table 5
Sample Printout from SWOPS (DL = 1, PRNT = 1)
L r)XY(?FN r>E«AWl - DISTANCE INCBF.MFNTS 40 TCI 60
1-77
KB
0.7500
FC
0 . 1 000
DSY
78 . oonn
X1
4n.nonn
ONT
704. nnoo
o . nnon
o, nonn
LOT
1 i n.nnnn
o, noon
o.onon
LNT
n.nnoo
0 . oonn
n.nnnn
KM
o.oooo
FN
0. 1000
Y1
6 .onoo
XT.
60.0000
1794
0
0
1 1 0
0
n
0
0
0
. 0000
. 0000
.0000
.0000
. ooon
.onoo
.nnno
,0000
.0000
0.9800
FD
o . 1 nno
Y7
4.onno
1 294.0000
O.onon
o.nooo
1 1 o.nnoo
o.nnno
n, onno
o.onno
o.nnno
o.onoo
KD
0.7500
VFL
24.0000
Y3
7.0000
1794
0
0
1 10
0
0
0
0
0
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
0
0
1
0
0
0
0
0
0
0
n
0
DX
.0000
BF.N
.0000
Y4
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
.0000
t)T
0.7500
POINT
50.0000
Y5
0.5000
o.oooo
0,0000
0,0000
0,0000
o.oooo
0.0000
0.0000
o.oooo
O.nooo
too.
950.
1 .
o.
0.
0.
0.
0.
0.
0.
o.
o.
XN
nonn
0
0000
ni,
0000
0000
0000
0000
nnoo
0000
ooon
0000
0000
0000
TM
50.0000
ZM
4,0000
PPNT
i .0000
0.
o.
0,
o.
0.
0,
o.
0.
o.
0000
0000
0000
0000
0000
0000
0000
oooo
0000
0,0000
0,0000
o.oooo
o.oooo
o.oooo
0.0000
0,0000
0,0000
0.0000
o.oooo
0.0000
0,0000
0,0000
0.0000
o.oooo
o.oooo
0,0000
0,0000
DEFICIT TN 5TPFJM . 100 / INITIAL BOO IN STREAM
"1
1
7
3
4
5
(,
7
n
9
10
1 1
17
1 3
1 4
15
40
0 .000
o.ooo
n. non
o.oon
n.oon
n.nno
n . nno
0.000
o.non
0 .oon
n . oon
n.oon
n.nnn
n .nnn
n .nno
41
0.000
o.ono
n.nno
n .nnn
1.000
n.nno
n.oon
0.000
n.ooo
n.noo
0.000
o.oon
o.ooo
o.nnn
n.nnn
42
n.nno
o.onn
o.nnn
o.nno
n.oon
n.oon
n.ooo
n.ooo
o.nno
n.nnn
n.nno
n.oon
O.onn
n.oon
o.nno
43
0.000
o.ono
o.nno
o.noo
o.onn
o.nno
n. non
n. oon
o.oon
n.nnn
n, noo
o.nno
O.ooo
n.nnn
n.nno
44
0,000
n.nnn
n . oon
O.nno
O.n *n
n. nnn
n.nnn
n.onn
n .non
O.nnn
n.nnl
n.nnl
n.nnl
o.nnT.
n . on?
45
n.ono
o.non
o.ono
0.000
n.oon
0,001
n.oo?
0.003
0.005
o.oon
0,012
n.017
n.073
o.n3o
0.039
46
o.ooo
o.noo
0.000
0.007
0.011
0.027
0.050
0.079
0.115
0.155
n.700
0.249
0.307
n.359
0.41S
47
0.000
0,000
0.007
0.1 16
0.393
0.639
0.874
t .099
1.311
1 .57?
1 .721
1 .913
2.097
7.275
7.4*7
48
0, 000
o.oo?
0. 1 36
0.395
0.648
0.899
t .1*5
1 .388
1 .6?7
1 .867
2.093
7.320
7.543
2.762
7.976
0
0
0
0
0
1
1
1
1
7
7
7
7
2
3
49
.002
.»"
.393
.646
.894
.138
.378
.614
.845
.073
.296
.515
.730
.941
.147
50
1.175
1.368
T.599
5.820
.037
.235
.431
.670
.801
.976
?.H6
!.309
2.468
S.671
J.769
51
0.002
0.010
0.075
0.047
0.076
0.109
0.1 48
0.191
0.238
0.289
0.342
0.399
0.459
0.571
0.584
57
n.ono
o.ooo
0.001
0,001
0.003
0.005
o.ons
0.012
0.017
n.n?7
0.079
0,037
0.046
0,056
0.067
53
0.000
0,000
0,000
0.000
0.000
0.000
n.noo
n.ooi
r. , o o i
o.not
0.007
0.002
0.003
o.no4
o.oos
54
o.ooo
0,000
o.noo
0,000
0.000
o.ono
n.ooo
o.ofln
0 . 000
o.ooo
o.oon
o.ooo
o.ooo
o.ooo
0,000
55
o.nno
o.ooo
o.ooo
0,000
0,000
0,000
o.noo
o.ono
n.onn
o.non
o.ooo
o ,000
o.ooo
o.noo
n.noo
56
0,000
0.000
0.000
o.noo
0,000
o.non
n.oon
o.noo
o.ooo
o.oon
o.ooo
o.nno
o ,nno
o.onn
n .nno
57
0.000
0.000
O.ooo
o.ooo
0,000
0,000
o.oon
0.000
o.non
o.ooo
o.ooo
o.oon
n.oon
o.ooo
o.oon
58
n.ooo
o.ooo
o.ooo
o.ooo
0,000
o.ooo
n.ooo
0.000
o.ooo
0.000
n.ooo
o.ooo
0.000
n.ooo
o.ooo
59
0.000
0.000
o.ooo
0.000
0,000
0.000
0.000
0,000
0 .000
0.000
0.000
0,000
0.000
0.000
0.000
60
0,000
o.ooo
0.000
o.ooo
0.000
0,000
o.ooo
0.000
0,000
0,000
0,000
o.ooo
0,000
0,000
0,000
-------
SECTION 5
STORMWATER OVERFLOW HYDRAULIC STREAM MODEL (SWOHS)
BACKGROUND
A digital computer program (SWOHS) has been developed to study the
dynamic (time-dependent) hydraulic response of the receiving stream to
large stormwater overflow hydrographs entering the stream. This problem
can be important in stream pollution studies because the reaeration co-
efficient is known to depend on both stream velocity and stream depth. The
two physical laws used in the program are conservation of mass and conserva-
tion of momentum. When these two laws are expressed as partial differential
equations they are sometimes (1) referred to as Saint-Venant's equations
after an early (1848) hydrologist who expressed them for the first time. The
program numerically integrates the St. Venant equations.
STREAM CONCEPTUALIZATION
In order to perform the integration, the stream is first divided into
a number of increments by planes perpendicular to the direction of flow. A
stream diagram with trapezoidal cross-sectional areas is shown divided by the
solid lines in Figure 10. Terms used in Figure 10 are defined in Table 6.
It is convenient to visualize the stream as a series of nodes (or junctions)
connected by channels. The nodes are visualized as having the property of
mass storage and the channels are visualized as having the properties assoc-
iated with flow of water. For example, properties such as elevation (above
sea level) of the stream bottom, elevation of the water surface, and stream
surface area are associated with the nodes. Channels, on the other hand, are
visualized as having the properties of length, width, velocity, roughness,
and side-slope if the trapezoidal shape is assumed. A rectangular channel
can be assumed by setting the value for the side-slope equal to zero. In
Figure 10, stream increments bounded by the solid lines are represented by
nodes and the increments bounded by the dashed lines are represented by
channels. The node is visualized as being located at the center of the
increment bounded by the the solid lines and the length of the channel which
connects two adjacent nodes is, therefore, one-half the length of the two
adjacent increments bounded by solid lines. Elevation of the stream bottom
(Z) and elevation of the water surface (Y) are associated with the node. The
depth of the node is simply Y-Z. The surface area (A) is also associated
with the node. The depth of the channel connecting two nodes is taken as
the average depth of the two nodes as shown by the first equation in Figure
10. The inter-connection between nodes and channels is described by the
43
-------
QIN
AA
Y(3)
Z(3)
1. CD(I)
2. CA(I)
3. WP(I)
4. HR(I)
5.
6.
7. A(5)
CD (3)
BB
CA(3)
CW(3)
4-
\
CD(3)*SS(3)
CC
Z(4)
J
//////
CD(I)*CW(I) + CD(I)**2*SS(I)
CW(I) + 2*CD(I)*(1 + SS(I)**2)**0.5
CA(I)/WP(I)
(CW(I-l) + 2*CD(I-l)*SS(I-l))*CL(I-l)/2 + (CW(I) + 2*CD(I)*SS(I))*CL(I)/2
CW(1) +«2*CD(1)*SS(1))*CL(1)
(CW(5) + 2*CD(5)*SS(5))*CL(5)
Figure 10 . Diagram of solution reach used in SWOHS showing cross-sections and basic equations
used for depth, cross-sectional area, wetted perimeter, hydraulic radius, and surface
area.
-------
Table 6
DEFINITIONS FOR VARIABLES USED IN THE SWOHS PROGRAM
NCASE
LIST
XI = NI
TN = NT
DT
XNCH = NCH
PRIN = NPRIN
TIDE = NTIDE
FMC1
FMC2
Tl, T2, T3, T4,
T5, T6, T7
TP
STR = JSTRT
STO = JSTOP
XSTRM = NSTRM
QSTRM(I)
number of cases to be computed by the program
alpha-numeric title of the data case
number of nodes
number of time increments to be calculated
length of time increment, seconds
number of channels
time increment at which the program printout
is to begin
identifying number of the tidally forced node
a constant having a value of 7.285 when English
units (cfs and sec) are used and 4.9 when metric
units are used
a constant having a value of 16.1 when English
units (cfs and sec) are used and 4.9 when metric
units are used
coefficients in sinusoidal expression for elevation
of tidally forced node as a function of time
tidal period, hours
time increment at which the storm input begins
time increment at which the storm input ends
number of node where storm flow enters
volume of storm flow at time increment I during
the duration of the storm, cfs (cu m/sec)
45
-------
Table 6 (continued)
1(1) elevation (above some reference point) of the
stream bed at the I'th node, feet (m)
Y2(I) elevation of the water surface at the I'th
node at time point T+DT, feet (m)
QSS(I) volume flow of the steady stream into the
I'th node, cfs (cu m/sec)
CV2(I) velocity of the water in the I'th channel at
time point T+DT, feet/second (m/sec)
CL(I) length of the I'th channel, feet (m)
CW(I) width of the flat bottom of the stream bed
for the I'th channel, feet (m)
RN(I) Manning roughness coefficient for the
I'th channel
SS(I) side-slope of the I'th channel
NDAM(I) vector to indicate the location of the dam/weir
controlling the flow from the I'th node - if the
integer I appears in the I'th row of NDAM a
dam/weir controls the flow out of the I'th node
W1(I) first weir coefficient for the I'th node
W2(I) second weir coefficient for the I'th node, ft (m)
W3(I) third weir coefficient for the I'th node
NCN(I,1) number of the upstream node for the I'th channel
NCN(I,2) number of the downstream node for the I'th channel
NCN(I,3) set equal to one if the I'th channel enters
the downstream channel laterally
CV1(I) velocity of the water in the I'th channel at
time point T, feet/second
Y1(I) elevation of the water surface at the I'th node
at time point T, feet (m)
CD1(J) depth of J'th node at time T, feet (m)
46
-------
Table 6 (continued)
CD2(J) depth of J'th node at time T+DT, feet (m)
Al one-half surface area of the channel at time T,
sq ft (sq m)
A2 one-half surface area of channel at time T+DT,
sq ft (sq m)
CAT cross-sectional area of channel at time T,
sq ft (sq m)
CA2 cross-sectional area of channel at time T+DT,
sq ft (sq m)
Q01 flow leaving a node through all downstream channels
at time T, cfs (cu m/sec)
Q02 flow leaving a node through all downstream channels
at time T+DT, cfs (cu m/sec)
QI1 flow entering a node from all upstream channels at
time T, cfs (cu m/sec)
QI2 flow entering a node from all upstream channels at
time T+DT, cfs (cu m/sec)
QIN(I) steady flow entering the I'th node over and above the
flow entering through numbered channels and is equal
to QSS(I) + QSTRM(I), cfs (cu m/sec)
NO number of upstream node
NI number of downstream node
WP1 wetted perimeter of channel cross-sectional area
at time T, feet (m)
WP2 wetted perimeter of channel cross-sectional area
at time T+DT, feet (m)
HR1 hydraulic radius of channel equal to the cross-sectional
area/wetted perimeter at time T, feet (m)
HR2 hydraulic radius of channel equal to the cross-sectional
area/wetted perimeter at time T+DT, feet (m)
47
-------
matrix NCN(NCH,3) with one row for each channel. The first column of NCN
contains the number of the upstream node (IN) and the second column contains
the number of the downstream node (NO). The third column refers to the
direction at which side streams enter the main stream and this will be
discussed later in this report.
PHYSICAL LAWS USED
The two independent variables used in the problem are time and distance
along the stream. The Euler coordinate system which uses a point attached
to the river bank as the zero distance point is used for the computation
because stream properties such as roughness coefficients and shape of the
trapezoidal cross-section are associated with points along the stream bed.
As mentioned earlier, the purpose of the program is to apply the principle
of conservation-of-mass to each node and the principle of conservation-of
momentum to each channel. For example, each node can be thought of as a
control volume and if the net flow entering the node (flow entering - flow
leaving) is positive the water depth in the node will increase. Alterna-
tively, if the net flow into the node is negative the water level in the
node will decrease. The principle of conservation-of-momentum applied to
each channel requires that if the resultant of the force system acting on
the water within the channel acts upstream the flow in the channel will be
decelerated and if the net force acts downstream the water will be acceler-
ated. The forces acting on the mass within the channel consist of the
gravity force, the pressure acting across the vertical dividing planes, the
frictional force exerted by the stream bed and the force equivalent to the
net momentum flux entering and leaving the control volume. For clarifica-
tion, these concepts will be discussed with respect to the stream diagram
shown in Figure 10.
SYSTEMS OF UNITS USED
Definitions given in the following table are expressed in the British
gravitational system which used pounds as the unit of force, feet as the
unit of length, and seconds as the unit of time. When the British system
is used the value of 7.285 for FMC1 and a value of 16.1 for FMC2 are
supplied as input. When the metric system is used, the values for FMC1 and
FMC2 are both 4.9. The metric gravitational system uses kilograms as the
unit of force, meters as the unit of length and seconds as the unit of time.
Corresponding units for other variables used in the program are Iftted below:
British Metric
1. Force Ibf kgf
2. Length ft m
3. Time sec sec
4. Acceleration of ?
Gravity, g 32.2 ft/sec 9.8 m/sec
48
-------
5. Mass (force/g) slug metric slug
6. Volume cu ft cu m
7. Velocity ft/sec m/sec
8. Volume Flow cfs cu m/sec
9. Area sq ft sq m
CONSERVATION-OF-MASS
To apply the conservation-of-mass principle to each node the surface
area of the node must first be computed. The nodal surface area is taken
as one-half the surface area of the two adjacent channels except for terminal
nodes with only one adjacent channel. For terminal nodes the nodal area is
taken as equal to the area of the one adjacent channel. This is demonstrated
by equations 5-7 in Figure 10. The net flow into each node is then computed
and divided by the nodal area to find the rate of change of the elevation of
the water surface in the node. When this rate is multiplied by the time
increment (DT) an incremental change in water surface elevation is found.
When the flow out of a node is controlled by a dam or a weir the weir
equation must be used to compute the flow out of the node. No distinction
between dams and weirs is made in the program because the flow is computed
from the following relationship in both cases.
QOUT = W1*(Y(5)-W2)**W3 (2)
In this relationship QOUT is the flow over the weir (cfs), Y(5) is the
elevation of the water surface in the 5'th node and Wl, W2, and W3 are weir
constants supplied as input to the program.
CONSERVATIQN-OF-MOMENTUM
To apply the principle of conservation-of-momentum to each channel it
is first necessary to compute the average channel depth, the channel cross-
sectional area, the wetted perimeter, and the hydraulic radius as shown by
equation 1-4 in Figure 10. Having computed these variables the mass of water
in the channel is known and the principle of conservation-of-momentum can be
applied by dividing the net force acting on the water by the mass to find the
time rate of change of velocity within the channel. For example, the gravity
force on the I'th channel can be computed as follows:
FG = 62.4*CA(I)*(Z(NI)-Z(NO)) (3)
In this equation NI is the number of the upstream node and NO is the number
of the downstream node.
49
-------
The force caused by pressure acting across the upstream and downstream
boundry planes can be written as follows:
FP = 62.4*CA(I)*((Y(NI)-Z(NI))-(Y(NO)-Z(NO))) (4)
It can be seen that the gravity and pressure forces can be added to give the
following simplier result.
FPG = FG+FP = 62.4*CA(I)*(Y(NI)-Y(NO)) (5)
The frictional force which acts in a direction opposite to the direction of
flow can be written as follows:
FF = 62.4*CA(I)*CL(I)*CV(I)**2*RN(I)**2/2.21/HR(I)**(4/3) (6)
Finally, a momentum force must be computed if the momentum flux entering the
channel differs from the momentum flux leaving the channel. For the simple
case where momentum flux enters the I'th channel only from the upstream
(I-l)'th channel the momentum force (FM) can be computed as follows:
FM - 62.4/32.2*(CA(I-l)*CV(I-l)**2 - CA(I)*CV(I)**2) (7)
The incremental change in velocity in the I'th channnel (DV) can now be found
from the following relationship. The mass of water (M) in the channel is
computed simply as 62.4*CA(I)*CL(I)/32.2.
M*(DV/DT) = FPG + FM - FF (8)
A number of complications are involved and these will be discussed later
on, but the equations presented here show in a rudimentary way the
physical laws used in the program.
NUMERICAL INTEGRATION SCHEME USED
Because of the quadratic terms which appear in equation (6) for the
frictional force (FF) an implicit numerical integration scheme such as the
Crank-Nicolson method could not be used for this problem. Instead, an itera-
tive scheme was used which is based on making the derivatives used for
advancing the solution one time increment (DT) equal to the mean of their
values at T and T+DT, Also, the conservation-of-mass equations and the
conservation-of-momentum equations must be solved simultaneously. Therefore,
the value of each variable at T and T+DT must be separately identified. For
example, Y1(I) is the water elevation in the I'th node at time = T and Y2(I)
is the water elevation in the I'th node at time = T+DT. All necessary vari-
ables are similarly marked.
Values for all variables at time = 0 must, of course, be supplied as
input to the program. To advance the solution by one time increment (DT)
the value of all variables at time = T+D are set equal to their values at
time = T. This constitutes the first guess at the value of the variables
at time = T+DT.
50
-------
All of the channels are then considered one-at-a-time. The upstream
node of the I'th channel (NI) is identified and value for Y2(NI) at the
later time point (T+DT) is found from equations for conservation-of-mass
using the most up-to-date estimates for all variables at the T+TD time
points. The downstream node (NO) for the I'th channel is identified and a
new value for Y2(NO) is found in the same way. The difference between the
new computed values for Y2(NI) and Y2(NO) and the previous values is noted.
If they differ by more than a set tolerance a flag (NTEST) is set equal
to one.
Continuing our consideration of the I'th channel, all forces acting on
the water within the channel are evaluated using the mean of values at
time = T and the most up-to-date values at time = T+DT. The mean mass of
water in the channel between T and T+DT is computed and equation (8) is used
to compute a new value of CV2(I) for T+DT. The newly computed value is
compared to the previous value and if they differ by more than some set
tolerance the flag (NTEST) is set equal to one. This procedure is performed
for each channel in turn.
After all channels have been considered and new values for Y2(J) and
CV2(I) have been found the flag NTEST is examined and if it has been set
equal to one the whole procedure is repeated using values for Y2(J) and
CV2(I) computed during the previous iteration. This iteration is repeated
until NTEST no longer equals one indicating that the computation has
converged. The final values for Y2(J) and CV2(I) are then assumed to be
the correct values and are printed. The numerical integration scheme is
then advanced one time increment and the whole procedure described above
is repeated until the iteration converges. These values are printed and
the solution advances to the next time point until the specified number (NT)
of time increments have been computed.
DESCRIPTION OF INPUT REQUIRED
The time increment to be used (DT) is read in as seconds and the number
of time increments over which a solution is to be found is read in as NT.
The number of channels in the stream system is read as NCH. The placement
of nodes and channels is described by a matrix (NCN) with NCN rows and three
columns. The I'th row of NCN, corresponding to the I'th channel, contains
the identifying number of the upstream node (NI) in the first column, the
number of the downstream node (NO) in the second column, and a one in the
third column if the I'th channel enters the main stream laterally. If a zero
appears in the third column the direction of flow in the I'th channel is
assumed to be in the same direction as the flow in the downstream channel.
This distinction is made to properly handle the momentum flux term.
Two constants, FMC1 and FMC2, are read-in to provide the flexibility
of changing the system of units. Hhen the British system (Ibf, ft, sec) is
used the values for FMC1 is 7.285 (32.2/2/2.21 )and the value for FMC2 is
16.1 (32.2/2). When the metric system (kgf, meter, sec) is used the value
for both FMC1 and FMC2 is 4.9 (9.8 m/sec2 divided by 2).
51
-------
Before the descriptive parameters for the stream are read in, the
analyst should assure himself that the variables used to describe the stream
are feasible. That is, they must be capable of corresponding to some steady-
state solution. Otherwise, the numerical integration procedure can drive
Y(J) below Z(J) corresponding physically to the node running dry. A program
HYSS described in the Appendix has been developed to check the input to
guard against this problem. When a feasible steady-state solution is known
the elevation of the stream bed, Z(J), and the elevation of the water surface,
Y2(J) is read into the program for each node at time = 0. For each channel
the length, CL(I), the width of the trapezoidal stream bottom, CW(I), the
side-slope, SS(I), the Manning roughness coefficient, RN(I), and the stream
velocity, CV2(I) are read into the program.
A steady flow into each node is provided for by the input variables
QSS(J) where J is the node number. Normally, QSS(J) is provided only for the
terminal nodes although this restriction is unnecessary. The hydrograph of
the stormwater overflow is divided into time intervals (DT) of the same
length as the interval used in the integration scheme. The volume flow (cfs)
read from the hydrograph at the mid-point of each interval is read into
the program as the vector QSTRM which provides for a maximum of 50 values.
The number of the node to which QSTRM is applied is input as NSTRM. The
variable L indexes time in the program. When L equals JSTRT the first
value of QSTRM is applied to node NSTRM. Successive values of QSTRM are
applied for all L values greater than JSTRT until L exceeds JSTOP. Thus,
the first and last values of L over the storm duration are JSTRT and JSTOP.
The number of the node which is tidally controlled is input as NTIDE.
The elevation of the node numbered NTIDE is computed as follows:
Y(NTIDE = Tl + T2*SIN(WT) + T3*SIN(2*WT) + T4*SIN(3*WT) + T5*COS(WT)
+ T6*COS(2*WT) + T7*COS(3*WT) (9)
WT = 6.283185*L*DT/TP/3600 (10)
Thus, the constants T1-T7 must be supplied as input and the period in hours
(TP) must also be supplied as input. In this expression for the elevation
of the NTIDE node; L is the number of time increments since time = 0.
When flow into or out of a node is controlled by a weir or a 4am
(hydraulically they are both treated the same) this is specified by the
input vector NDAM with J rows; one for each node. When the flow out of the
J'th node is controlled by a weir/dam this is indicated by placing the number
J in the J'th row of NDAM. If flow into the J'th node is controlled by a
weir/dam the number of the node (KD) which controls the flow into the J'th
node is placed in the J'th row of NDAM. The three weir coefficients are
provided by the vector W1(J), W2(J), and W3(J) where J is the node number
of the node controlled by the weir/dam. The flow out of the J'th node is
computed as follows:
52
-------
QO = W1(J)*(Y(J) - W2(J))**W3(J) (11)
Flow into the J'th node controlled by the KD'th node is computed as follows:
QI = W1(KD)*(Y(KD) - W2(KD))**W3(KD) (12)
If a dam/weir is located in a non-terminal node no channel is assigned
downstream of the node. The logic of the program is arranged to assume that
the momentum flux of the channel upstream of the dam/weir is dissipated by
the dam/weir and no incoming momentum flux is assigned to the channel down-
stream of the node on the downstream side of the dam/weir.
LOGICAL STRUCTURE OF THE PROGRAM
When input to the program is completed, the first step in advancing the
solution from time = 0 to time = DT is to set the values of node water surface
elevation and channel velocity at time = DT equal to their values at time = 0.
This constitutes the first guess at the solution at time = DT. At any time
other than zero this procedure involves setting Y1(J) equal to Y2(J) and
CV1(I) equal to CV2(I). Thus, in order to use the same logic in the program,
initial values for Y and CV are input as Y2(J) and CV2(I).
The iterative numerical integration procedure for advancing the solution
by one time increment continues by considering each channel (I) one-at-a-time
and first identifying the node (K) downstream of the I'th channel. The first
column of NCN is scanned to identify channels which flow out of node K
and the second column of NCN is scanned to identify channels which flow into
node K. When a channel which flows into or out of node K is identified the
other terminal node of the channel is identified from NCN and called M.
The depth, cross-sectional area, and surface area of any channel flowing
into or out of node K is then computed from the properties of nodes K and
M. The variable NSEG is given the value of the total number of channels
which flow into or out of node K.
A test is made to determine if a dam/weir controls flow into or out of
node K. If NDAM(K) equals K the flow out of node K over a weir is computed
from the inputted weir coefficients and the elevation of water surface in
node K. If NDAM(K) is not zero, and not K, it must equal KD and the flow
into node K is controlled by the water level in node KD. The flow from node
KD into node K is computed from the inputted weir coefficients Wl(KD),
W2(KD), and W3(KD). If NSEG is greater than one, the surface area of node
K is taken as one-half the sum of the surface areas of all incoming or
outgoing channels. If NSEG is less than two, the surface area of the one
channel entering or leaving node K is taken as the surface area of node K.
The net flow into node K is then computed. Flow into node K includes
the sum of the flow in all channels entering node K, the steady flow QSS(K),
the flow entering over a weir from an upstream node, and stormwater overflow.
If K equals NSTRM and the time index (L) falls between JSTRT and JSTOP
inclusive, the volume flow of stormwater entering node K is given by the
53
-------
vector QSTRM. The flow out of node K is found by summing the flow in all
channels leaving node K with flow over a weir, if one exists at node K.
Subtracting the flow out of node K from the flow into node K gives the net
in-flow which can then be divided by the surface area of node K to find the
change in water surface elevation at node K. The new value for Y2(K) is
then compared to the previously computed value YTEMP(K) and if they differ
(in absolute value) by more than a specified tolerance (0.01 ft is used in
the program) the flag NTEST is set equal to one. In either case YTEMP(K)
is set equal to the newly computed value for Y2(K).
If NTIDE equals K the whole procedure outlined above is bypassed and the
level of water in node K is found from the tidal equation (9). The same
procedure is then repeated for the node upstream of the I'th channel unless
that node is tidally controlled.
Conservation of momentum in the I'th channel is considered next. To
consider conservation of momentum in channel I, the depth, cross-sectional
area, wetted perimeter, and hydraulic radius of channel I are computed from
the most up-to-date values of variables involved. The mass of water in the
I'th channel is computed as follows:
M = CA(I)*CL(I)*62.4/32.2 (13)
The derivative of velocity over the time interval DT is expressed as follows:
DVDT = (CV2(I) - CV1(I))/DT (14)
The principle of conservation of momentum as applied to the I'th channel is
expressed as follows:
M*MVDT = FPG + FM - FF (15)
Where FPG is the sum of the gravity and pressure forces acting on the mass
M, FF is the frictional force acting on M, and FM is the equivalent momentum
force caused by the net momentum flux entering and leaving the I'th channel.
To assure good accuracy in the numerical integration, all variables must
be averaged over the interval DT. That is, variables used in applying the
principle of conservation of momentum must be the mean of values existing at
T and T+DT. In terms of the average values of all variables, the three
forces acting on the I'th channel are expressed as follows:
FPG = 62.4*CA(I)*(Y(NI) - Y(NO)) (16)
FF = 62.4*CA(I)*CL(I)*ABS(CV(I))*CV(I)*RN(I)*^2/2.21/HR(I)**(4/3) (17)
FM = 62.4/32.2*(CA(J)*ABS(CV(J))*CV(J) - CA(I)*ABS(CV(I))*CV(I)) (18)
The absolute values are used in the expressions for FF amd FM to make sure
that the FF force acts in a direction opposite to the direction of flow and
54
-------
that the momentum force acts downstream when the momentum flux entering from
upstream exceeds the momentum flux leaving the I'th channel. Also, in the
expression for FM the first term in the brackets must be summed for all
channels (J) upstream of the I1th channel.
When equations 13, 14, 16, 17, and 18 are substituted into equation 15
it can be seen that the resulting equation is quadratic in terms of the
unknown variable CV2(I). This can be readily solved if care is taken to
select the proper root. If we let X be equal to CV2(I) the quadratic equa-
tion can be put into the following form:
AX2 + BX + C = 0 (19)
The values for the constants in this equation can be shown to be as follows:
A = 1/DT
B = 7.285*RN(I)**2/HR2(I)**(4/3) + 0.5/CL(I)
C = CP1 + CP2
CP1 = 7.285*ABS(CVl(I))*CVl(I)*RN**2/HRl**(4/3) (20)
- 16.1*(Y1(NI)+Y2(NI)-Y1(NO)-Y2(NO))/CL(I)
- CV1(I)/DT
CP2 = (ABS(CVHD) - FMI1 - FMI2)*CV1 (I)/CL(I)/2
FMI1 = CA1(J)*ABS(CV1(J))*CV1(J)/CA1(I)
FMI2 = CA2(J)*ABS(CV2(J))*CV2(J)/CA2(I)
The index (J) which appears in the expression for FMI1 and FMI2 must be summed
over all channels which flow into the I'th channel. If the number which
appears in the third column of NCN of the row corresponding to the J'th
channel is one, the direction of flow in the J'th channel is perpendicular
to the flow in the receiving channel and the momentum flux (FMI1 or FMI2) is
set equal to zero. If the value FMI1 is found to have a value of zero the
I'th channel is assumed to be a terminal channel and the values of FMI1 and
FMI2 are computed as follows: The upstream node of the I'th channel is K.
FMI1 = QSS(K)*CV1(I)/CA1(I) (21)
FMI2 = QSS(K)*CV2(I)/CA2(I) (22)
FMI1 will also be zero if the I'th channel is downstream of a dam/weir but
since QSS(K) will be zero the value for FMI1 and FMI2 will remain zero. Care
should be taken to make QSS(K) zero at any node downstream of a dam/weir.
This type of calculation is made for each of the NCH channels and if the
55
-------
newly computed value for CV2(I) differs from the previously computed value
VTEMP(I), by more than the tolerance (0.001 ft/sec) the flag NTEST is set
equal to one and a second iteration is necessary. If NTEST is found to be
zero, indicating that the iteration has converged, Y1(I) is set equal to Y2(I)
and CV1(I) is set equal to CV2(I) and the numerical integration is advanced
one time increment. This procedure is continued until the solution has been
found for NT time increments. The results of the computation are printed
after each time step.
SOLUTION FOR HYPOTHETICAL PROBLEMS
The first problem selected for testing SWOHS involved a prismatic
stream flowing at a constant velocity into which a square hydrograph of
15 minutes duration and a volume flow roughly equal to the stream flow was
introduced. The trapezoidal cross-section of the stream had a width of
30 ft and a side-slope of four. The solution reach was divided into 20
lengths of 4000 ft. each. The Manning roughness coefficient was set at
0.023 and the bed slope at 0.0005 ft/ft. Before entry of the square hydro-
graph, the steady stream flow was 598.2 cfs, the velocity of flow was 3.036
ft/sec and the depth in each node was 4.207 ft. The length of the time
interval used in the numerical integration was 300 seconds.
To simulate the effect of a large stormwater overflow a square hydro-
graph with a duration of 15 minutes and a volume flow of 600 cfs was
introduced into the 5'th node. The computed elevations of the water
surface in each of the 20 nodes above the steady-state elevation is shown
for one-hour intervals after entry in Figure 11. From Figure 11 it can be
seen that the maximum increase in water surface of 1.2 ft occurred at the
time corresponding to the end of the square hydrograph. One hour later
the pulse has spread significantly and the amplitude has reduced to about
0.35 ft. In general, the square hydrograph produced a wave which moved
downstream at a velocity of about 4 ft/sec and attenuated rapidly.
The computed flow over the weir in node 20 minus the steady-state
flow of 598.2 cfs is shown plotted in Figure 12. Thus, in about 11.4 miles
of stream the peak of 600 cfs introduced at node 5 has been reduced to
about 46 cfs. The volume of flow introduced (540,000 cu ft) equals the
area under the hydrograph at node 20 as shown by the cumulative curve in
Figure 12. The interval over which the 600 cfs overflow was introduced is
also shown in Figure 12 between the 10'th and 12'th time interval^. Thus,
about 5 hours elapsed between the time the overflow was introduced and the
time at which the flow over the weir began to increase. The fact that the
integrated flow over the weir equaled the integrated flow into the reach
shows that the conservation-of-mass principle is being satisfied. A
simple closed-form expression for the computed result is not available to
check the numerical integration scheme.
One of the reasons for simulating the hydraulic response of the stream
to large stormwater hydrographs was to consider the effect on the reaeration
coefficient. Many relationships are available for estimating the reaeration
56
-------
1234567
8 9 TO 11 12 13 14 15 16 17 18 19 20
Node Number
Figure 11 . Computed wave profiles caused by a square hydrograph
(600 cfs-15 min duration) entering a steady-state
stream with 598.2 cfs volume flow.
57
-------
50
40
30
20
10
i.
O)
o
t
O)
o
>
600 qfs
'400
«
(D
500
300
200
;|ioo
-;;8:::;;;<
0
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140
Time, number of 15 min intervals
Figure 12. Computed flow over terminal weir minus steady-state flow resulting from a
square hydrograph (600 cfs-15 min duration) introduced 11.4 miles upstream
of the weir.
-------
coefficient as a function of the depth and velocity of the stream. The
work of O'Connor and Dobbins (2) is widely used for this purpose. The two
relationships developed by O'Connor and Dobbins, one for deep and one for
shallow streams are shown by equations (23) and (24) below:
Ka = 1440 (DU^/H1-5 (H greater than 5 ft) (23)
K = 1016.4 D^/H1'25 (H less than 5 ft) (24)
Ka = reaeration coefficient, days
D = diffusivity, sq ft/hr (80 x 10"6 sq ft/hr typical value)
H - water depth, ft
S = slope of stream channel, ft/ft
U = stream velocity, ft/sec
If the typical value of 80 x 10" sq ft/hr for diffusivity is substituted
into equations (23) and (24) the following two equations result:
K = 12.9 l/VH1'5 (25)
a
K - 48.5 SWH1'25 (26)
a
In the first example the steady state parameters for the stream were
3.036 ft/sec for velocity, 4.207 ft for depth, and 0.0005 for slope.
Substituting these values into equations (25) and (26) yields an estimate
of 2.6 days" from equation (25) and 1.2 days' from equation (26). Since
the depth of the stream is less than 5 ft equation (26) should be used
according to O'Connor and Dobbins.
The greatest change in the stream characteristics occurred just after
the 600 cfs event in node 5 where the depth increased by 1.3 ft and the
velocity increased by 0.73 ft/sec. Substituting these into equation (25)
gives an estimate of K of 1.94 instead of 2.6 days" . From equation (26)
1 -1
the estimate for K is 0.86 days instead of 1.2 days . Thus, the
a
change is about 25% in both cases and this is the peak variation which
lasts a very short period of time. Since the correct value of the reaeration
coefficient is probably not known within 25% it can be concluded that for
the problem selected the effect of the hydraulic transient on the value
of K, is negligible.
a
The second hypothetical problem was taken from the report (3) on the
RECEIVE-II model written by the Raytheon Company. This problem which is
59
-------
described in Section V of the Raytheon report involves one tidally controlled
terminal node and one dam/weir located in the center of the solution reach.
A diagram for the problem is given in Figure 13. The solution reach is
divided into ten 2000 meter lengths. The dam/weir is located in the 5'th
node and the 10'th node is tidally controlled. The steady flow into the
upstream terminal node is 5.66 cu m/sec and a second steady flow of 0.566
cu m/sec enters the 3'rd node. The solution showed that the elevation of
the water surface in nodes 1-5 and the stream velocity in channels 1-4
approach a steady-state condition. For nodes 6-10 the water surface
elevation as a function of time beginning at 750 minutes and ending at
1500 minutes from time = 0 is shown in Figure 14. The water surface
elevation in node 6 is essentially constant and the water surface elevation
in node 7 varies only slightly. The water surface elevation in nodes 8
and 9 varies significantly in phase with the variation in the tidally
controlled node number 10. These results are in substantial agreement with
the results of the RECEIV-II model. One significant difference between the
two models is that the surface areas are supplied as input in the RECEIV-II
simulation. The input for the two problems and a sample output at one time
point are shown in Tables (7) and (8).
PROGRAM LISTING
The FORTRAN listing for the program SWOHS is shown in Table 9 and
the definitions for variables used in the program is given in Table 6.
Examples of input and output are shown in Tables 7 and 8.
60
-------
O)
.(->
QJ
OJ
-------
LoJ
800
900
1000
1100 1200
Time, minutes
1300
1400
1500
Figure 14. steady-state variation of water surface elevation of nodes 6-10 in the second hypothetical
problem solved with the SWOHS program.
-------
Table 7
Printout for First Hypothetical Problem Solved with SWOHS
STREAM MODFT, TfST CASE - DYN 6 M 1C /HYDR AUL 1C (1)
(SWflHS) 9-12=1977
XT
70.nnn
T1
o.ooo
STB
l n.non
OSTBM
(inn , nonn
o.nnnn
n .ooon
o.onon
n.ooon
T
1
7
1
4
5
6
7
p
9
1 n
1 1
17
1 3
t 4
15
t ft
17
1 0
1 0
7n
T
1
7
1
4
S
A
7
0
<)
1 n
1 1
1?
1 •>
1 4
1 S
1 6
\ 7
1 R
1 9
2n
NT
25.ono
T7
n.non
sin
i 2,000
6nn.nnoo
o.oooo
n .nnon
n, nnnn
o.ooon
ZCT1
Rnn.oooo
7
-------
Table 7 (continued)
10
JTs 50
JTi 55
60
I
t
2
3
4
5
*
7
4
9
10
1 1
(7
1 3
14
15
1 6
1 7
IB
19
20
I
1
2
3
4
5
6
7
8
9
10
1 1
1 1
1 3
14
15
16
1 7
1?
19
?0
T
1
T.
3
4
5
A
7
A
9
10
1 1
12
13
1 4
15
IS
17
in
19
30
CV2(I)
3.0365
3.0365
3.0185
2.8044
3.3330
3.1076
3.0497
3.0386
3.0369
3.0366
3.0365
3.0365
3.0365
3.0365
3.0365
3.0365
3.0365
3.0365
3.1223
0.0000
CV2CT)
3.0365
3.0366
3.0354
2.4361
1.6837
3.2849
3.0980
3.0499
3.0393
3.0371
3.0366
3.0365
3.0365
3.0365
3.0365
3.0365
3.0365
3.0365
3.1224
o.oooo
CV2CI)
3.0365
3.0369
3.0127
2.22J6
3.7682
3.4653
3. 1883
3.0797
3.0474
3.0391
3.0371
3.0366
3.0365
3.0365
3.0365
3.0365
3.0365
3.0365
3. 1224
O.ooon
Y2CT5
904.3070 ,
802.5070
800.2070
798.5009
796.8307
794.2712
792.2174
790.2097
788.2073
786.2070
784.2070
782.2070
780.2070
779.2070
776.2070
774.2070
772.2170
770.2070
768.2070
766.0357
Y2CI)
804.2070
802.2070
800.2069
798,2142
797.J55J
794.4501
792.2638
790.2192
798.2095
786.2075
784.2071
782.2070
780.2070
778,2070
776.2070
774.2070
772.2070
770.2070
768 .2070
766.0356
Y2CD
804.2070
802.2070
800.2057
798.2685
797.5016
794.681 5
792,3601
790.2491
788,2174
786.2094
784.2075
782.2071
780.2070
778,2070
776.2070
774.2070
772.2070
770.2070
768.2070
766.0356
QTNf I)
599.2000
0,0000
0.0000
0.0000
600.0000
0.0000
0.0000
o.oooo
o ,0000
0.0000
0.0000
o.oooo
0.0000
0.0000
o.oooo
o.oooo
0.0000
0.0000
o.oooo
0.0000
GTNCI)
598.2000
0.0000
o.oooo
0.0000
600.0000
o.oooo
0.0000
0.0000
o.oooo
o.oooo
o.oooo
o.oooo
o.oooo
0.0000
0,0000
0,0000
0,0000
0,0000
o.oooo
o.oooo
QINf 11
598.2000
o.oooo
0.0000
o.oooo
600.0000
o.oooo
0.0000
o.oooo
o.oooo
o , oooo
0 . 0000
0,0000
o.oooo
0.0000
o.oooo
0.0000
0.0000
0 .0000
0.0000
0.0000
CO
598,2054
598,2047
598.0182
608.6623
731 , 1630
619.6157
601 .7906
598.8104
598 .3096
590.2323
599.2027
599.2008
598 . 2006
598.2006
59S.2006
599.2006
599.2006
598,2006
598,1823
0.0000
CO
598, 1995
598.2045
598,6759
564.4816
883.2436
678.7817
617.1514
602.2786
599.0403
598.3716
598.2379
598.2034
598.2009
598.2007
598.2007
598.2007
598,2007
598.2007
599.1871
0.0000
cn
598.2005
598,1810
599,3221
537,8660
966.3295
753.2759
648.048S*
611.8713
601 .6057
599.0009
599.3951
598.2395
598,2034
598.2009
598.2006
590.2006
598.2006
598.2006
598, 1906
0.0000
CA2(I1
197.0035
197.0074
196.8132
21 7.0403
219.3700
199.3869
197.3922
197.0696
197.0152
197.0074
197.0035
197.0035
197.0035
197.0035
197.0035
197.0035
197.0035
197.0035
191.5818
0,0000
CA2(I)
197.0035
196.9996
197.2289
231.7133
239.7731
206.6395
199.2073
197.4738
197. 1006
197.0229
197.0074
197.0035
197.0035
197,0035
197.0035
197.0035
1 97,0035
197.0035
191 .5818
0.0000
CA2(I1
197.0035
196.9685
198.9303
242.0036
256.4422
217.3754
203.2580
198.6908
197.4155
1 97. 1006
197.0268
197.0074
197.0035
1 97.0035
197.0035
1 97,0035
197.0035
197.0035
191.5819
0.0000
64
-------
Table 8
STREAM MnnFL TFST CASK - DYNAMTC/HYDRAULIC C2t
tSWOHS)
9-12-1977
CT>
cn
XT
1 0.000
T1
0.001
STP
o. ooo
Ov5TPM
0,0000
0 .0000
0 . OOOO
0. OOOO
n.oooo
T
1
7
3
4
5
6
7
8
9
10
NT
75.000
T7
-0.116
STn
0.000
0.0000
0.0000
o.oooo
0.0000
0.0000
zm
1 7.0000
1 0.0000
p. oooo
6,0000
4,7000
7.4666
0.7131
-1 .oooo
-2.0000
-3.0000
nr
300.000
Tl
-O.B78
XSTRM
0.000
0,0000
0.0000
0. OOOO
0.0000
0, OOOO
Y2f 11
1 3.0000
1 1 .0000
9.0000
7 .0000
6.7000
4.0000
2.0000
1 , oooo
1 .0000
1 .0000
XNCH
ft. 000
T4
-0.047
0.0000
0.0000
0.0000
0.0000
0.0000
OfiS(T)
5.6600
0.0000
0.5660
0.0000
0,0000
o.oooo
o.oooo
0.0000
0.0000
o.oooo
PRIN
0.000
T5
• 0.090
0.0000
0.0000
o.oooo
0.0000
o.oooo
CV?(I)
o. i ooo
0. 1 000
0. 1 000
0. 1000
0. 1000
0. 1 OPO
0. 1 000
0. 1 000
0.0000
o.oooo
TIDE
10.000
T6
-0.475
0.0000
o.oooo
o.oooo
0,0000
o.oooo
CL(I1
7000.0000
2000.0000
2000,0000
2000,0000
2000,0000
2000,0000
2000,0000
2000,0000
0.0000
0,0000
FMC1
4.900
T7
0.101
0,0000
0.0000
0.0000
0,0000
0,0000
CWfl)
40.0000
40.0000
40.0000
40.0000
60.0000
RO.OOOO
100.0000
170.0000
o.oooo
o.oooo
FMC7
4.900
TP
25.000
o.oooo
0.0000
0.0000
0,0000
o.oooo
RNf 11
o. 1000
0.1000
0, 1 000
0.1000
o.iooo
0.1000
0. 1000
0.1000
0,0000
o.oooo
0,0000
o.oooo
0,0000
0,0000
o.oooo
SS(I)
0,0000
0.0000
0,0000
0,0000
o.oooo
0,0000
0.0000
0,0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
NDAW(I)
0
0
0
0
5
5
0
0
0
0
Tj =
T
1
7
1
4
5
6
7
R
9
10
HI
0
0
0
0
73
0
0
0
0
0
,1T =
(11
.0000
.0000
.0000
.0000
.5000
.oooo
.0000
.0000
.0000
.0000
5 MTN
W2f T 1
o.onoo
o.oooo
0.0000
o.oono
6.7000
o.oooo
0 . OOOO
0. OOOO
o.oooo
O.oooo
t
1
7
1
4
*,
6
7
8
9
10
cv7r T
0.
o.
o.
0,
0,
0.
0.
0,
o.
o.
W3U)
o.oooo
o.oooo
o.oooo
o.oooo
t .5000
0.0000
o.oooo
o.oooo
0.0000
o.oooo
NCN
-1 NCN-2
1
7
3
4
6
7
8
9
0
0
) Y7U1
394H
30RO
3979
3000
4R79
3605
1074
5697
oooo
oooo
1 7
10
9
6
6
3
1
1
0
-0
.9R44
.9905
.0072
.9919
.7441
.9410
,9fi'i7
.0093
.8H1 1
.3710
OINdl
5.6600
0.0000
0.5660
0.0000
0,0000
0.0000
0.0000
0.0000
0.0000
0 . OOOO
2
3
4
5
7
R
9
10
0
0
1
NCN
CO
15.6629
15.9351
15.8684
1 R. 71 36
39,5048
4R ,7174
25.0308
180.8472
0.0000
0,0000
-3
0
0
1
0
0
0
0
0
0
0
CA?m
39,6769
40.0348
39.8820
60.7195
81 ,8022
I 30.46R4
244.5321
333.4991
0,0000
0,0000
-------
Table 9
FORTRAN Listing for the SWOHS Program
(SWOHSI STREAM MODFL --- DYNAMIC /HYDRAULIC
P. G. FILERS SEPTFMPEP 1977
DIMENSION LISTC40),7OO),Y1f30),Y?(30),QlN(30),CDJf30),CD2(30),
. CV1 MO),CV2f 30),CLMO),CWf30),RNOO),SS(30) ,YTEMPf30) ,VTEMP(30),
. CM (30),CA2MO).QSvSf 30).W1 f 30) ,W2( 30) f W3 f 30) ,NCNf 30, 3) ,
. NDAMf 30) , QSTRM(50)
OPEN f UNITs 1 ,NAVF.B» SWOHS.DAT I ,TYPFs'OLD» , FOPM« I FORMATTED ' ,
. READONLY)
TNsl
1086
READriN,10) NCASF.
10 FOPMATfT?)
DO 1000 IIT»lfNrASF
DO 20 I«l,30
Zf I)«0.
Yl fl)a0.
Y2(T)sO.
QTNf T)«0.
CD1 (T)«0.
Cn2(T)«0.
CV1 f T)«0.
CL(T)*0.
CWf T)sO.
RNfT)*0.
SSf I)*0.
YTFMp(T)=0
VTFMPf I)sO
CA1 f I)sO.
f A2(I)«0.
OSSf I )»0.
W1 (I)«0.
W2ri)*0.
W3f I)*0.
NCNf Tf 1)*0
l, 2)«0
20 CONTINUE
DO 2S Ts1 ,«50
OSTPM(I)sO.
25 CONTINUE
PFADfTN,40) LIST
40 FORMAT(40A2)
PFADfTN,50) XI,TN, DT,XNCHf PRIM f TIDE, F«C1 ,FMC2
PFADfTNf«>0) Tl,T?,T3,T4,Tsi.T6,T7,TP
66
-------
Table 9 (continued)
READfIN,50)
50 FnRMAT(8F10,
NXIaXT
NT*TN
NPPINaPPIN
NCHsXNfH
NTir>E = TIDE
OSUVz25.
JSTRTaSTR
60
70
100
1 10
1 15
116
119
120
STR,STO,XSTPM
0)
NSTPMsXSTRM
PFADfIN,60)
RFAHfIN,60)
PEAHfIN,60)
RFADfIN,60)
READfIN,60)
BFAOfIN,60)
PEAPfIN,60)
READfIN,60)
READfIN,70)
PEAOfIN,60)
READfIN,60)
PEADfIN,60)
fOSTPMfT) , Is1,50)
f Z f T ) , T a 1 , N X T )
fY2fI),T*1,NXT)
fQSSfI),lsl,NXT)
fCV2fT),T=1,NXT)
fCLf T),Isl,NXn
fRNfT)
,NXI)
,NXI)
fSSfT),I=1,NXT)
fNDAM(T),Isl,NYI)
fMlfT),T*1,NXT)
,NXT)
T)
fNCNfT,1 ) ,Ts1 ,NXT)
I,3),T*1,NXT)
PEAD(TN,70)
PFADf IN, 70)
PEADf IN, 70)
FOPMAT(IOIfl)
WRITFf in, 100) LIST
FORMAT (1 HI , /////, 20X.40A?,//)
WPTTEfTn,l10) XT,TN,DT,XNCH,PRJNf TIDE,FMC1
FORM AT f8X, 'XI » , 10X, 'NT' , 10X, 'DT1 ,8X, 'XNCH' ,8X,
•PRTNi),8Xr'TIOE',PX,'FMCl »,PX,|FMC2',/,8F12,3,/)
WRITE fid, 11 5) T1,T?,T3,T4,T5,T6,T7,TP
FORM AT f 8X, «Tl',10X,'T2',lOY,'T3',10X,1T4',10X,'T5',10X,'T6'>
10X,'T7', 10X,»TP',/»8P'12.^»'/)
WRTTFf in, 11 6) 5TR,STO,XSTRM
FORM AT f 7X, 'STR' ,9X, 'STH' ,7V, ' XSTRM • ,/,3Fl2.3,//,5X» 'QSTRM» )
WRITFf 10,1 1 8) fQSTPMf I),Tz1 ,50)
FDPMATf 10F12.4)
WRITFCIO, 120)
FORM AT f //, 11X, 'I' ,6X. 'Z(T) ' »7X» (Y2(I) ' ,6
. 'Ct.(Tl ' ,7X, 'CW(T) ' ,7X, "RNf T) ' »7X, 'SS(I)
HO 140 1*1, NX!
WRITE CIO, 130) T,Z(T),Y2f T),QSSf I ) , CV2 C I )
. NDAWfT)
130 FORMftT(8X, 14, 8F1 2.4,112)
i) ' ,6X, «CV2fI) ' ,7X,
Xi 'NDAM(I) ')
( I ) ,CW 1 1 ) »RN f I) ,SS f I ) ,
67
-------
Table 9 (continued)
140 CONTINUE
WRTTFflO, 150}
150 FOPMATf //,1 1X,'T<,5X, I,W1 m,W2t n.Wlf I).NCN(I, I ),NCNf T,2),NCN(T,3)
160 FOPMATf8X,T4,3F12.4,3H2)
J70 CONTTMUE
LL*0
KOIJTasO
PC 600 L=l, NT
TF fL-JSTRT) ?00,lflO,180
180 IF fL-JSTOP) 1^0,190,200
190 LL«l,L+t
200 ICNT*0
,TTaL#DT/60.
DO 210 I si, NX I
vi m=Y2m
CVl (T^sCV2f I)
210 CONTINUE
TF (NTIDF) 220,220,215
21s; Y1 fNTIDE)»T1*T5-fT6 + T7
Y2(NTIDE)*Y1 (NTTDE)
220 NTFSTsO
ICNTsTCNTf 1
DO 470 1*1, NCH
NIsNCNfl,!)
NOsNCNf I,?)
KsNO
230 A1aO.
A2»0.
001*0.
002=0.
NSFG*0.
DO 360 J«t,NCH
IF (NDAM(K)-K) 239,250,239
239 IF (NCN(J,1)-K) 270,240,270
240 NSFGasNSFC,+ l
MsNCNf J,?)
CD1 (J)s(Yl (MUY1 (K)»Z(M)-7fK))/2
CA1 C.nsCDl fJl#CWfJ)4-CDl(J)##2«SS(J)
68
-------
Table 9 (continued)
Ln#*2*SS( J)
A1 (J1*CV1 CJ>
A2d,n*CV2r
-------
Table 9 (continued)
400 K s N T
GO Tfl 2.30
410 NTswrNf TP 1 )
CA1 msCD] f T5*CWm*rD1 (T)#*2*SS(I
HR1 sCAl (T)/WP1
?f D/WP2
P=1
fYtfNI)4Y2(NT).Yl
FMTIsO.
KsNCNf T, 1 )
DO 440 Js! ,NCH
TF fNCNfJ,?)-K) 440,4?0,440
420 TF fNCN£J,;n-n 410,440,4^0
430 FMT1sFMIH>CAl ( J ) #CV1 f J) *APg f <" V! ( J ) ) /C A 1 ( I )
440 CONTINUE
IF (FMin 444,442,444
44? FMii*Qssno*cvi f n/CAi m
FMT2BOSS(K)#CV2(T)/C»2(I)
444 rp2«fFMT1fFMT2-CV1fI)«ABSfrvlfI)))/CLfI)/2.
IF fr) 446, 44?, 448
446 A«ABSfA)
GO TH 44Q
44P
44Q
TF fABS(VTEMP(I).CV2ttn-.OOn 460,460,450
450 NTFSTsI
460 VTEMP(nsCV2f I)
470 rONTTNUE
TF fTCNT«25) 4SO,,500,500
430 TF fNTEST-1) 500,220,220
SOO IF (!>NPPIN) 600,510,510
510 WPTTF(TO,5?0> L,,TT
520 FOBMRTf/^X,1!. s'»T4,3X,,MTss«,l4,1X,'MIN>,6X»'Il,5X,«CV2(T)ll.
. 7X,»Y2n)',6X,sOTW{T5!,10X,'CQf»6X,'CA2(I)')
KKKaNT+1
DO 540 1*1, KKK
T)«CA2(T)
70
-------
Table 9 (continued)
WPITF. f TOf 530) I,CV2(I),y?fT)»OIN(T)fCO,CA2fT)
530
540
IF (ROUT) 1000,600.1000
600
1000
CLOSFflJNIT»l
CMP
71
-------
REFERENCES FOR SECTION 5
Chow, V. T., "Open-Channel Hydraulics, McGraw-Hill Book Company, 1959.
O'Connor, D. J. and Dobbins, W. E., "Mechanism of Reaeration in Natural
Streams" ASCE Transactions, Vol. 123, 1958, pp. 641-684.
Raytheon Company, Oceanographic and Environmental Services, Portsmouth,
Rhode Island, New England River Basins Modeling Project, Volume III,
Part 1, RECEIV-II Water Quantity and Quality Model, December 1974, EPA
Contract No. 68-01-1890.
72
-------
APPENDIX FOR SECTION 5
DESCRIPTION OF HYDRAULIC STEADY-STATE (HYSS) PROGRAM
To avoid computational problems in the use of SWOHS initial values for
water surface elevation, Y(J), at all nodes and the velocity, CV(I) in all
channels should correspond to a feasible steady-state solution. The program
HYSS has been developed to find an initial feasible steady-state solution.
The variable names used in HYSS correspond, where possible, to variable names
used in SWOHS. Variables read into HYSS as input include QOUT (cfs) which
is the volume flow out of the downstream terminal node. The downstream
terminal node is assumed to be controlled by a dam/weir and the three weir
coefficients, Wl, W2, and W3 must be read in as input. A maximum nodal depth,
CDMAX, in feet must be supplied as input. The number of channels in the
stream is supplied as NCH. Two constants, VAR1 = 32.2 and VAR2 = 14.57
corresponding to th English system of unit are read in as input. The steady
flow into each node, QIN(J), is read in as input. The following additional
variables whose names correspond to input for SWOHS must be supplied for
each channel; CL(I), CW(I), RN(I), and SS(I).
The first thing to be computed is the water surface elevation in the
downstream terminal node, Y(KI), from the value for QOUT and the weir/dam
coefficients. The flow in all channels is found next. This is done by
taking each channel (I) in turn, identifying the upstream node (NI) and
then setting the flow in the I'th channel equal to the sum of QIN(NI) plus
the flow in all channels whose downstream node is NI. Care should be taken
to make the number of any channel which flows into the I'th channel less
than the number of the I'th channel.
Initial values for Y(J) and CV(I) are set equal to zero. Each channel
is then considered one-at-a-time beginning with the channel upstream of the
terminal node (KI) and working backward through the list of channel numbers.
The upstream node (NI) and the downstream node (NO) are identified for each
channel. The elevation of the water surface in the downstream node is known.
The elevation of the water surface in the upstream node is known to be
between Z(NI) and Z(NI)+CDMAX. A continuously halving or Bolzano iteration
is then performed to find the elevation of the water surface Y(NI) in the
upstream node. This iteration is performed for each channel in turn finding
the channel velocity and the elevation of the water surface for each node.
As a check, DVDT, the rate of change of velocity with respect to time is
computed and this should always be zero or a very small number.
73
-------
If the variables input to the program are incompatible with a feasible
steady state solution, this will be detected by noting that the elevation
of the water surface equals Z(J) (within the tolerance) or equals Z(J) plus
CDMAX. If this happens new values for RN, CW, or SS must be input and the
program rerun until a feasible solution is found. The definitions for the
variables used in the HYSS program are given in Table 10. The listing for
HYSS is shown in Table 11 and the output for the first and second hypothetical
problem is shown in Table 12.
74
-------
Table 10
Definitions of FORTRAN Symbols Used in HYSS
NCASE number of cases to be executed by the program
LIST alpha-numeric title of the data case
QOUT sum total of all flows entering the system, cfs
Wl first weir coefficient
W2 second weir coefficient
W3 third weir coefficient
CDMAX maximum feasible depth of channel, ft (M)
XNCH = NCH number of channels
VAR1 a constant having a value of 32.2 when English
units (cfs and sec) are used (9.8 metric units)
VAR2 a constant having a value of 14.57 when English
units (cfs and sec) are used (9.8 metric units)
QIN(I) steady flow entering the I'th node, cfs (M/Sec)
CL(I) length of the I'th channel, ft (M)
CW(I) width of the flat bottom of the stream bed
for the I'th channel, ft (M)
RN(I) Manning roughness coefficient for the I'th channel
SS(I) side-slope of the I'th channel
Z(I) elevation (above some reference point) of the
stream bed at the I1th node, ft (M)
NCN(IJ) number of the upstream node for the I'th channel
75
-------
Table 10 (continued)
NCN(I,2) number of the downstream node for the I'th channel
NCN(I,3) set equal to one if the I'th channel enters
the downstream channel laterally
76
-------
Table 11
FORTRAN Source Listing for HYSS
fHYSS) PROGRAM FOR PREPARING INPUT TO SWOHS
P,G. EILEPvS SEPTEMBER 1976
DIMENSION LTST(40)fZ(30),QTNf30),Ct.f30),Cwr30),RN(3
-------
Table 11 (continued)
Y(Kn=W2+(QOUT/Wl
no 75 jsi,NC"
DO 70 K=1 ,NCH
IF (NCNfK,2)-Nn 70,65,70
65 Of.TlsOfJUQfK)
70 CONTINUE
75 CONTINUE
78 DO 80 Is1 ,NCH
VTFMPd)sCV(I)
80 CONTINUE
no ?oo TT=I ,NCH
I«KT-TT
NTaNCNCI, 1 )
NO»NCN(I, ?)
YMAXxZfNT J+CDM&X
100
WP«CWd)+?.«CDd5#n ,+SSdl **?)*#. 5
HR*CAd)/WP
CVf I)aO( T)/CAd)
FMTsO.
KsNCNd, 1 )
DO 170 J»1,NCH
IF f NCNf J,?5-K) 120,115,120
115 IF (NCN(J,.1)-n 117,120,117
117 FMT»FMI+CAf J)#CVf Jl#«2.
120 COMTTNME
IF fFMI) 126,125,126
125 FMIsQTN(K)#CVd)
126 DVnTf T)BCONST4 f FMI-FMO VC A d 1 /CL ( T 1
TF fnVHTfin 140,140,130
130 YMAXBYfNIl
GO TO 150
140 YMINsY(NT)
150 y(NI^sfYWAx + VMI^3^/?.
TF fARSfYTEMP. YCNin-. 00001 ) 200,200,100
200 COMTTNUF
MTFSTaO
DO 220 Isl,NCH
TF f ARS(VTFMPd)-CV(m-.OOn 220,220,210
210 NTFSTsl
78
-------
Table 11 (continued)
220 CONTINUE
IF fNTFST) 2<30,?Qn,7P
290 WRTTFdO, 3005 LIST
300 FOPMATMH1,/////,20X,40A2,//)
WRTTFfin,310)
310 FORMATflSX, 'GOUT',10X, 'W1 ',10X,'W2',10X,'W3',7X,
. •CDMAX',8X,»XNCH',flX,»VAP1»,flX,'VAR2')
WRITEfTO,320) QOHT,WJ ,W2fW3,fnMAX,XNCH,VARl,VAP2
120 FOPMATf9X.8F12.4,///I
WRTTFfin,310)
330 FOPMATC8X,(Tl,5X,tOTN',8X(,'CLl,8X,'CW«,8X,'RNSRX(,lSv^l,9
. 3X, «MfN-1 ' , IX, 'NCN-?' ,1X, 'MfN-31 ,7X, 'Y' »8X, TD«,flX» «CVI
. 6X, 'DVDT1 ,/)
DO 3*10 1 = 1,KT
WPTTFfln,140) I,OIN(n,CL(n,rW(I),RN(I),SS(I),ZfI),
. NCNf T,1),NCN(T,?),NCNfT.3l,YfT),CD(T1»CVfT),PVDTfI)
340 FnRMATf3X,T6,6FtO,3,lT6,4F1083)
350 fONTTNT'F
1000 CHNTTNIJE
ClOSFCUNTTsl1
FND
79
-------
00
O
Table 12
Printout for the First and Second Hypothetical Problems Using HYSS
FIRST HYPOTHETICAL PROBLEM
(HYSS)
9-12-1977
OOIIT
Ml
W7
59R.7196 20.0000 756.4ooo
I
l
7
1
4
•>
6
1
8
g
m
1 1
t 2
1 1
1 4
1 5
1ft
1 7
t a
1 1
70
QTN
598.240
n.onn
0.000
0.000
n.ono
n.noo
0 .000
0.000
0.000
0. 000
o.ooo
0 . (100
0.000
0.000
0.000
0.000
0.000
0.000
o.ooo
0.000
CL
4000. 000
4000.000
4000 .000
4000.000
4000.000
4000. 000
4000.000
4000.000
4QOO.OOO
4000.000
4000.000
4000.000
4000.000
4000.000
4000.000
4000.000
4000.000
4000.000
4000.000
o.ooo
cw
'0.000
10.000
10.000
10.000
30.000
10,000
10.000
10.000
10.000
50.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
10.000
50.000
0.000
PN
0.023
0.023
0,023
0.023
0.023
0.023
0.023
0.023
0.023
0.023
0.023
0.023
0.023
0.023
0,023
0.023
0.023
0.023
0,023
0.000
W3
1 ,5000
SS
4.000
4 .000
4,000
4.000
4,000
4.000
4.000
4 .000
4.000
4.000
4.000
4.000
4.000
4.000
4.000
4.000
4.000
4.000
4.000
0.000
CDMAX
1 0,0000
Z
BOO. 000
79R.OOO
796.000
794.000
79? ,000
790.000
7RR.OOO
7fl6.ooo
7R4.000
7P2.000
700,000
778 .000
776.000
774.000
772.000
770.000
76B.OOO
766.000
764.000
767.000
KNCH
19,0000
NCM-1
1
2
3
4
5
6
7
R
9
10
1 1
1 2
1 3
14
15
16
I 7
IB
19
0
NCN-2
2
3
4
5
6
7
S
9
10
1 1
1 2
1 3
14
15
16
17
1 fl
19
20
0
VAP1 VflR2
32.
NCN-3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
2
0
,2000 14.5700
Y
R04 .707
802.207
R00.707
798.207
796.207
794.207
792.207
790.207
788.207
786.207
7R4.207
782.207
7R0.207
778.707
776.207
774.207
772.207
770.207
768.207
766.036
CD
4.207
4.207
4,207
4.207
4.707
4.207
4.207
4.207
4.707
4.707
4.707
4.707
4.207
4.207
4.207
4,207
4.207
4,207
4.122
o.ooo
cv
3,037
3.037
3.037
3.037
3.037
3.037
1.037
3,037
3,037
3.037
037
037
037
037
037
037
037
036
122
o.ooo
DVPT
o.ooo
0.000
0,000
0,000
o.ooo
o.ooo
o.ooo
0.000
0.000
o.ooo
0.000
0.000
0.000
0,000
o.ooo
o.ooo
o.ooo
0,000
0.000
o.ooo
HYPOTHETICAL PROBLEM
(HYSS)
q.12-1977
omit
21 .0000
i
t
7
3
4
5
6
7
a
f)TN
16
0
0
0
7
0
0
0
.000
,000
,000
.000
.000
.000
.000
.000
Wl
70.0000
W7
488.9000
CL 4 C»l
1 2700.
17500.
10200.
1*100.
37500.
3750.
3600.
0.
000
000
000
000
ooo
000
ooo
000
73
77
43
41
11
10
10
0
.000
.000
.000
.000
.000
,000
.000
.000
0
0
0
0
0
0
0
0
W3
1.5000
CDM6X
7.0000
XNCH
7 .0000
VAR1
17.7000
NCM-I NCN-2 NCN-3
VAR7
14.5700
cn
CV
DVDT
73.000
77.000
43.000
41 .000
11 .000
10.000
10,000
0.000
O.OSO
0.080
0.060
0.080
0.080
0.050
0.050
0.000
6.000
2.000
3.000
1 .000
1 .000
4 , 000
4.000
o.ooo
575.000
560.000
537.000
570.000
673.000
495.000
494 .000
488.000
1
2
3
5
4
6
7
0
2
3
4
4
6
7
8
0
0
0
0
0
0
0
2
0
575. 7Rt
56] .715
537.429
570.666
673.000
496.869
494.000
489.998
0.998
0.822
0.548
0.333
1 .768
0.935
0.999
0.000
0.553
0.679
0.655
0.508
0.562
0.729
0,677
o.ooo
0.000
o.ooo
o.ooo
0.241
0.000
0.000
0.016
o.ooo
-------
SECTION 6
EFFECT OF STORMWATER ON STREAM DISSOLVED OXYGEN
BACKGROUND
From an analysis of stormwater overflow measurements made in seven
cities, Heaney et al (5) concluded that the BOD load originating from
separate storm sewers in an urban area depends primarily on the population
density and ranges from 25-38 Ib BOD/acre-yr (28-43 kg/ha-yr). The cor-
responding load from combined sewers was found to be 3-4 times this amount.
Using the load estimating relationships developed by Heaney et al, the
pollution caused by combined sewer overflows in a community of 100,000
can be estimated at about 100 Ib BOD/acre-yr (112 kg/ha-yr). If the
community treats the municipal sewage to the secondary level (30 mg/1 BOD)
the load from the steady dry weather flow can be estimated at about 92 Ib
BOD/acre-yr (103 kg/ha-yr) assuming a population density of 10 persons per
acre. Thus, the amount of pollution from combined sewer overflows is
roughly equivalent to the load from the dry weather secondary treatment
plant. This type of reasoning has generated an interest in estimating
the impact of stormwater overflow events on the dissolved oxygen resources
of the receiving stream.
If dispersion in the receiving stream is neglected, the Streeter-
Phelps model (7) can be used to find the peak dissolved oxygen deficit
caused by any overflow event with minimal computational effort. However,
if the analyst wishes to include the effect of dispersion, no model of
equivalent simplicity has been available.
If the stormwater overflow event is divided into increments as shown
in Figure 15 and the stream is prismatic, each incremental load can be
approximated as an impulsive load for which a closed form solution is
available (1). Since the partial differential equations governing the
deoxygenation and reaeration process are linear the superposition principle
applies and the response of the incremental loads can be summed to find the
response of the entire stormwater overflow event. An alternative approach,
which must be used when the stream cannot be assumed to be prismatic is to
numerically integrate the governing partial differential equations. This
approach is tedious and time-consuming but when the pollutional pulse is
square and the stream is prismatic numerical integration results can be
conveniently summarized by means of non-dimensional groupings providing
preliminary estimates of the peak dissolved oxygen deficit occurring in the
stream as well as the time after entry of the pollutional pulse at which
the peak occurs.
81
-------
0
0 1
23456789
Time, hours
Figure 15. Typical stormwater overflow event divided into half hour sub-loads
-------
PROBLEM OF DEOXYGENATION AND REAERATION IN TUBULAR STREAMS
When the Euler coordinate system, which measures distance along the
stream (X) from a point on the stream bank is used, the following two
partial differential equations describe the concentration of BOD (L) and
the concentration of dissolved oxygen deficit (D) as functions of distance
and time.
3L/3t = -v(3L/3X) + E(32L/3X2) - Ky,,L (1)
3D/3t = -v(3D/3X) + E(32D/3X2) + K ,L - K D (2)
Q a
These equations can be simplified by converting to the Lagrange coordinate
system which measures distance (x) from a point which moves downstream with
the velocity of the stream. When the Lagrange coordinate system is used,
the first term on the right-hand side of both equations becomes zero and
the Lagrangian equations take the following simplier form:
3L/3t - E(32L/3x2) - KpL (3)
3D/3t = E(32D/3x2) + K,L - K D (4)
a a
If dispersion (E) in the stream is assumed to be zero, distance is no
longer an independent variable and equations (3) and (4) reduce to the
following two equation which govern the Streeter-Phelps model:
dL/dt = -KrL (5)
dD/dt = K.L K D (6)
u a
Equation (5) can be integrated simply as follows:
L/LQ = e-V (7)
If the hydrograph of the storm overflow entering the stream is square with
a flow rate of Q. and a BOD concentration of L..., the initial BOD concentra-
tion (L ) in the stream is calculated by the following relationship:
"0'
L0 = Q.L./Q (8)
83
-------
The Streeter-Phelps model was developed for steady state analysis of
deoxygenation and reaeration in streams but can be applied equally well to
a plane in the stream perpendicular to the direction of flow in which the
initial BOD concentration is given by equation (8). Equation (7) shows
that the BOD concentration in that plane falls off exponentially as the
plane moves downstream at the velocity of the stream. The DO deficit in
the plane will increase from the initial value (DQ) to some maximum value
and then approach zero as the time form entry becomes great. The relations
describing the behavior of DO deficit can be found by substituting equation
(7) into equation (6) and integrating to find the following result:
D/LQ = K
When disperison is assumed to exist in the stream values for L and D can
no longer be associated with planes and their values are functions of
both time and distance along the stream.
STREAM RESPONSE WITH DISPERSION
When BOD is introduced into the stream by a storm overflow event, the
BOD and DO deficit concentrations produced in the stream are functions of
time and distance from the point in the stream at which the load is intro-
duced and are referred to as the BOD and DO deficit responses in the stream
to the storm overflow load. The form of the BOD and DO deficit responses
must satisfy equations (3) and (4) and will also depend on the nature of
the load introduced.
Since equations (3) and (4) are linear partial differential equations,
the principle of superposition applies. The storm overflow load can be
divided into a number of sub-loads, the response of each sub-load solved
for, and the response of the total load found as the sum of the responses
of the sub-loads. Therefore, if the solution for a sub-load is known, the
solution for any overflow event can be found by superposition. The most
convenient kind of sub-load is the impulsive load. An impulsive load can
be visualized physically as a load which enters the stream instantaneously.
A following closed-form expression for the BOD response to an impulsive
load is given by Thomann (8).
L' = "/A , e- e- (10)
(2EtT "
This equation, which is a solution to equation (3) for an impulsive load,
shows that the shape of the BOD concentration profile at any time has the
shape of the normal probability density function with a standard deviation
\,
equal to (2Et)2. The peak value of the profile is infinite initially and
falls off quickly with time. Since equation (3) is analogous to the
84
-------
classical one-dimensional heat conduction equation, methods for deriving
the solution given by equation (10) can be found in textbooks such as
Operational Mathematics by Churchill (3).
To use equation (10) the storm overflow hydrograph and the corresponding
plot of BOD concentration versus time are both divided into a number of
equal time increments. The value of M for each time increment is then found
by multiplying the average flow (0^) over the interval by the average BOD
concentration (L.) and multiplying the product by AT. This procedure is
shown diagrammatically in Figure 15.
The cross-sectional area of the stream (A) is found by dividing the
stream flow (Q) after introduction of the stormwater by the stream velocity
after introduction of the stormwater. Adjusting the units gives the
following modification of equation (10) applied to a single impulsive load.
2nJ
4 (v AT/s) e'
(11)
In equation (11) s is defined as (2Et)
(8).
2 and L is computed using equation
A stormwater overflow pollution stream model (SWOPS) was developed
as described in ref. 6 using conventional numerical integration methods
as described by Bunce and Hetling (2). The DO deficit profiles in the
stream resulting from square stormwater overflow events were computed
using a wide range for the various variables involved. It was noted from
the computed results that, as the width of the square pollution pulse was
decreased, the shape of the DO deficit profiles at any time approached
the shape of the normal probability density function. Assuming that the
shape of the DO deficit response at any time has the shape of the normal
probability density function it follows that the second partial derivative
of D with respect to distance could be expressed as follows:
82D/3x2
= (D/2Et) ((x/2Et) -
(12)
Substituting equations (11) and (12) into the governing partial differential
equation (4) then produced the following expression for the DO deficit
response in the stream to an impulsive load of BOD when DQ is zero.
D/L,
J_
2n
(vAT/s) (Kd/(Ka-Kf))
(13)
The second term on the right side of equation (9) should be added to
equation (13) if DQ is positive.
85
-------
A rigorous derivation of equation (13) is contained in the referenced
(1) paper by Bennett. A recurrence formula derived by Di Toro (4) which
relates the concentrations of sequentially reacting substances to the
concentration of a single first order reacting substance can be used to
derive equation (13) from equation (11).
To compare the accuracy of the two methods of solution the DO deficit
response to a square pollution pulse was computed with the digital computer
stream model SWOPS (6) and with equation (13) using the superposition
principle. A listing of the program SWOCFS developed to apply equation (13)
using the superposition principle is shown in Table 13. A sample printout
is shown in Table 14.
A comparison of values computed by the two methods is shown in
Figure 16. The values plotted are the peak values for D/LQ at 15-minute
intervals after the leading edge of the square wave enters the stream.
Values computed with SWOPS are shown by the solid line and values computed
using equation (13) are shown by the circled points. In this analysis,
o
the hypothetical stream had a volume flow of 950 cfs (26.9 m /sec) and a
velocity of one mile/hr (1.61 km/hr). The values used for rate constants
were 0.25 days"1 for both Kf and Kd and 0.98 days"1 for K . The width of
the square pulse was taken as one hour and the time increment used was 15
minutes. The volume flow of the pulse was 1294 cfs (36.6 m /sec) and the
BOD concentration was 110 mg/1. The difference between the two solutions
can be attributed to computational error.
NON-DIMENSIONAL GROUPINGS USED TO SUMMARIZE COMPUTED RESULTS
Although equation (13) can be used with the superposition principle
to find the complete DO deficit response to any known stormwater overflow
with a minimum of computation effort, many problems require only an
estimate of the peak DO deficit concentration in the stream resulting
from a square pollution pulse. However, the peak DO deficit concentration
might be required for hundreds or even thousands of individual square
pollution pulses. Therefore, a method which would allow estimation of
the peak DO deficit concentration simply by reading a graph is needed for
many problems. To fill this need the program SWOPS was used to compute
the peak DO deficit concentration from square pollution pulses over a
wide range of all variables involved. These computed results are shown
in Table 15. It can be seen from equations (9) and (13) that the term
Kd/(Ka-Kr) appears in both expressions for D and will, therefore, not
appear in the ratio. Thus, only Kr and Ka need be specified in computing
the ratio of the two dissolved oxygen deficits.
By analyzing the computed results it was found that, if the computed
peak DO deficit concentration was divided by the corresponding DO deficit
concentration computed from the Streeter-Phelps model, this ratio will plot
as a single-valued function of the following non-dimensional grouping.
86
-------
Table 13
Listing for the Program SWOCFS which Applies the
Superposition Principle to the Closed-Form Solution
for an Impulsive BOD Load
c r swncrs ) CLOSED FOPM SOLUTION FOR STREAM SIMULATION
C P. G. FILERS SEPTEMBER 1077
PEAL KP,KA,KD.LO(?00),DOf200),LOIN(30),MDT,LOMAX
DIMENSION LIST f40),QIN(30)
nPENniNIT*J,NAMEB«swnCFS.DATi , TYPE* 'OLD ' , FORMn • FORMATTED ' ,
. READONLY)
TNs1
10 = 6
READ(TW,10) NCASE
10 EORMftTfT2)
DO 1000 111 = 1 , NCASE
DO 1 5 Tsl ,200
Lom*o.
DOC i)»o.
15 CONTINUE
DO 16 Js1 , 30
OTNf I)=0.
LQTNfn=0.
1ft CONTINUE
PEADfTN,20) LIST
20
PEADfIN,10) DT,DN,TM, FORK
30 EORMATf 9F10.0)
READfTN,35) f QTN ( I ) , 1 = 1 , 10 1
PEADfIN,35) fLOTNm, IB! ,30)
35 FORMATflOFB.O)
NDeDN
MTsTM
WRTTFfTO,40) LIST
40 EOPMAT( 1H1 ,/////, 20X,40A2, //)
wpiTEf 10,50) KR,EC.VEL,QSS.KD,KA
50 FORM ATf flX, ' KP • , 1 OX , ' EC ' , 9y , * VEL ' ,
. 9X, 'OSS' ,1 OX, «Kn'f10X,'KA'./.6F12.4,/)
WPITFfIO,60t DT,DN,TM,EOPK
60 FOpMftTfQXj'DT'.IOX.'DNi.lOXf'TMSSX
. 4F12.4,//)
WPTTEfTO,65) fQTN(T),I»1 ,301
6$ FOPM ATf 7X, 'OIN',/,tOEl2,4,/,10Fl2.4,/,lOFl2.4e/)
WRTTFfIO,66) (LOINf T) , T«l , 30)
f>6 FORMATr6X,'LOIN»f/,10Fl2.4,/,10E12.4./,10Et2.4,/)
DTaDT/24.
87
-------
Table 13 (continued)
DXsVFL*DT'
TsO.
DO ?00 Mal.MT
TsT+PT
NXsM
JF fM.ND) 90,80,70
70 MXeNP.
80 DO P5 1*1 ,21
Lom*o.
DOf USED,
8^ CONTINUE
MPTsT
DO 1!SO KstjNX
DO 11 0 Tal ,21
TF f VftP-1 8,1 90,110,110
90 IF fFOPK1) Q6f9^,
-------
o
£
0)
O)
CD
X
O
•o
O)
o
CO
to
to
0)
D_
100 200 300 400
Number of 15 min intervals since entry
500
Figure 16. Comparison between computed dissolved oxygen deficit peak values
using the program SWOPS shown by the solid line and using the
closed form solution shown by the circled points. Time is measured
as time intervals (15 min) from the entry of the first increment
of the one-hour square pollution pulse.
89
-------
Table 14
Sample Printout from SWOCFS
S»MPLF TF.ST r»SF FOP THE CLOSFD FORM SOLUTION
1-12-1977
KP Ff VFt, OSS
0.7511 0.1001 74,0000 510.0000
PT DM TM FOPK
1.7500 16.0000 600.0100 1.0000
QTN
3.0100 1.8000 1.B001 1.8000
4,7000 4.7000 9.4000 9.4000
0.0000 0.0001 0.0010 0.0000
LOTN
7.7010 8.1000 8.1000 8.3000
5.5000 5.5000 5.5100 5.5000
1
2
3
4
5
6
7
8
9
1 0
1 1
12
1 3
14
15
16
17
IS
1 0
70
21
72
73
74
75
76
77
78
79
10
0.0000
0 ,OO
1,00
0 .00
0.00
0,00
0,00
1.10
0.00
0.00
0.00
0^13
0.04
1.04
0.05
1.06
0.07
1,09
1. 00
0.10
1.11
0.17
1.13
0.14
0.15
0.16
1.1'
0.18
0.10
0.70
0.71
1.00
0.00
1.00
1. 00
0 .00
0.00
0.10
0.00
0.00
0.13
0.04
0.04
0.05
0.06
0.07
O.OR
1,09
0.10
0.1 1
0.17
0.13
1,14
0.15
0.16
0.17
0.18
1,19
1.70
1.71
0,77
0. 0000
0, 10
0,00
0.00
0.10
0.00
0,00
0,00
0,10
0.03
0,04
0.04
0.15
0.16
0.07
O.OP
0.09
0.11
0.1?
0.11
0.14
0.15
0,16
0.17
0,1 P
1.19
0.70
0.71
0.77
1.23
o . 74
0.00
0. 00
0.00
0.00
(••.00
0.01
0.10
0.03
1.05
1,06
0,17
O.OP
1.09
0.10
fy. 1 1
0.1 3
0.14
0.15
0.16
o. 1 a
0.19
1.70
0.71
1.27
0.73
1,74
1,25
0,76
0.77
0.78
0.0100
0.10
0,00
0.10
0.00
0.00
0.00
0.13
0,15
0,o6
O,o7
0.08
0.09
0.10
0.17
0.13
0.14
0.16
0.17
0.18
0.19
«.?1
0.77
0.71
0.74
0.76
0.77
0.78
0.79
0. 31
O.-M
0.00
0.00
0.00
f'.OO
0.00
0.03
0.05
0.06
0.07
0.08
0.09
0. 0
0. 2
0. 3
0. 4
0. 6
0. 7
0, 8
0.70
0.71
0.72
0.73
0,74
0,76
0.77
0,78
0,79
0,10
0.11
0.33
o.oooo
o.oo
0,00
0.00
o.oo
0,03
0.03
0.06
0.07
0.08
0.09
0.10
0.17
0.13
0.14
0.16
0.17
0.18
0.20
0.21
0.27
0.24
0,25
0.26
0.27
0.29
0.30
0.31
0.37
0.33
0,34
0,00
0,00
0.00
0.03
0.05
0,06
0.07
0.08
0,10
0.11
0.12
0.14
0.15
0,16
0.1 P
0.19
0.70
0.27
0.73
0.24
r..25
0,77
0.28
0,29
0. 30
0.31
0,37
0.34
0.35
0.36
0.
5.
o.
0.
5.
5.
o.
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
n
0
0
0
0
KD
7500
ROOO
4000
0000
2000
5000
0000
.00
.01
.03
.05
.06
.07
.08
.10
.' 1
.'2
.'*
.15
.16
.'B
.'9
.20
.22
.23
,24
.25
.77
.28
.29
.30
.31
.17
.33
.14
.I*
.37
0,00
0.03
0,05
0.06
0,07
0,08
0,09
0.11
0.17
1.13
0.14
0.15
0,17
0.18
0.19
0.70
0.21
0,22
0.23
0.74
0.75
0.26
0.77
0.7B
0,?9
0.30
0.31
0.12
0.33
0.13
K»
0.9800
5.8000
9.4000
0.0000
5.2000
5.5000
0,0000
0,03
0.14
0.04
0,05
0,06
0,06
0,07
O.OP
0,09
0,09
0.10
0.11
0,12
0,12
0.13
O.J4
0.15
0,15
0,16
0.17
0.17
0.18
0,19
0.19
0.20
0.20
0.21
0.27
0.77
0.23
0.00
0.10
0,00
0,00
0,00
0,10
0.01
0.01
0.01
0,02
0,02
0,02
0.03
0,03
1,13
0,04
0,04
0.05
0.05
0.05
0.06
o.os
0,06
0.07
0,07
0.08
0.08
0,08
0,09
0,09
5.8000
0.0000
0,0000
5.2000
0.0000
0.0000
0.00
0.00
0.00
o.oo
0.00
o.oo
0,00
0.00
0.00
0,00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.10
0.10
0.1)
1.01
0.01
0.01
0.01
0.11
0.01
0,01
1,01
0.07
0.07
5.8000
0.0000
0.0000
5.7000
o.oooo
0,10
0.00
0 .00
0.00
0.00
0.00
0.00
0.00
0.00
o.oo
0.00
0.00
0,10
0,10
0.00
0.00
o.oo
0.00
0,00
0.00
0.00
0.00
0.00
n.oo
0.00
0.00
r. .00
o.oo
0.10
0.00
0,0000
0,00
0,00
0,00
0.00
0,00
o.oo
0.10
o.oo
0.00
o.oo
0.00
0,00
0.00
0.00
o.oo
o.oo
0.00
0.00
o.oo
0.00
0.00
0.10
o.oo
0.00
0.00
0,00
o.oo
0,00
0.00
0,00
4.2000
o.ooofl
0.0000
5,5000
0.0000
0.00
0,00
0.00
0.00
0.00
0.00
o.oo
0.00
o.oo
o.oo
o.oo
0.00
0.10
0.00
0.10
o.oo
0,10
o.oo
0.00
0.00
o.oo
0,00
0,00
o.oo
0,00
0,00
0.00
0,00
0,00
0.00
o.oooo
0.00
0.00
o.oo
o.oo
0,00
0.00
0.00
0.00
o.oo
0.00
0.00
0.00
0.00
0.00
0.00
0.00
o.oo
0.00
o.oo
o.oo
0.00
o.oo
o.oo
0.00
0.00
0.00
0.00
o.oo
0,00
1.00
4.2000
o.oooo
0.0000
5.5000
0.0000
0.0000
o.oo
0.00
0.00
o.oo
0,00
0.00
O.Ofl
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0,00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0,00
0.00
0.00
o.oo
0,00
0.00
0,00
0.00
0,00
0.00
0,00
0,00
0.00
0.00
0,00
0.00
0,00
0.00
0.00
0.00
0.00
0.00
0.00
0,00
0.00
0.00
0.00
0.00
0.00
0.00
o.oo
0.00
o.oo
0.00
o.oo
o.oo
0.00
0.00
0.00
0,00
0,00
0.00
o.oo
0,00
o.oo
o.oo
0.00
0.00
0.00
o.oo
0,00
o.oo
0.00
0.00
0.00
0,00
0.00
0.01
0,00
0.00
0,00
0.00
0.00
0.00
0.00
0.00
0,00
0.00
0,00
0,00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
o.oo
0.00
0.00
0,00
0,00
0.00
-------
Table 15
Computed Values for Peak Dissolved Oxygen Deficit arid
Time of Occurrence for Square Pollutographs
Velocity, mi/day
Pulse width (W) , mi
I K
Kr Kd
1 c
SP-t*
SP-DOMAX/LO
n
Ka
Kr - Kd
E
SP-t*
SP-DOMAX/LO
in
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
IV K
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
v
Ka
Kr Kd
. E
SP-t*
SP-DOMAX/LO
VI
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
VII
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
VIII
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
IX .,
Ka
Kr Kd
E
SP-t*
SP-DOMAX/LO
4.68
0.7
0.1
= 0.4774
0.1071
4.68
0.7
0.2
0.4774
0.1071
4.68
0.7
0.3
0.4774
= 0.1071
-.4.68
= 0.7
= 1.0
= 0.4774
- 0.1071
= 2.0
= 0.7
= 0.1
0.8076
0.1989
2.0
0.7
1.0
0.8076
0.1989
0.4
0.2
1.0
3.466
0.25
0.4
0.4
1.0
2.5
0.368
0.4
0.7
1.0
- 1.865
0.474
t*
DOMAX/LQ
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LQ
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/LO
N
T*
t*
DOMAX/Ln
N
T*
2.4
0.1
= 0.229
- 0.01658
= 2.351
= 2.99
= 0,229
= 0.01178
= 3.324
5.98
= 0.229
= 0.00963
4.071
= 8.97
0.229
0.00528
7.433
= 29.9
- 0.406
0,02316
2.907
4.06
0.354
0.00722
9.193
35.40
0.844
0.00344
18.803
84.40
0.719
0.00637
15.811
71.90
0,604
0.01014
13.747
60.40
6.0
0.25
.250
.03965
.940
.400
.240
.02878
1.330
.768
.240
.02371
1.628
1.152
.229
.01316
2.973
3.664
.427
.05647
1.163
.683
.406
.01835
3.677
6.469
1.635
. Oil 04
7.521
26.160
1.210
.01910
6.324
20.160
.937
,02860
5.499
14.992
12.0
0.5
.312
.06911
470
.125
.271
.05339
.665
.217
.260
.04503
.814
.312
.240
.02588
1.487
.960
.490
.10405
.581
.196
.417
.03636
1.839
1.668
1.729
.02213
3.761
6.916
1.281
.03811
3.162
5.124
.948
,05698
2.749
3,792
24.0
1.0
417
.09602
.235
.042
,354
.08435
.332
.071
.323
.07569
.407
.097
.260
.04868
.743
.260
.635
.16146
.291
.635
.448
.07021
.919
.448
1.760
.04389
1.880
1.760
1.312
.07536
1.581
1.312
,979
,11222
1,375
.979
32.15
1.14
.458
.10201
.175
.026
.406
.09452
.248
.045
.375
.08781
.304
.063
.292
.06146
.555
.163
.708
.17908
.217
.039
.469
.09081
.686
.261
1.792
.05827
1.403
.998
1.333
.09976
1,180
.742
1,010
,14796
1,026
.563
48.0
2.0
.490
.10569
.118
.012
.458
.10276
.166
.023
.437
.09940
.204
.033
.344
.07970
.372
.086
.781
.19247
.145
.020
.531
.12407
.460
.133
1.865
.08496
.940
466
1.406
.14428
.791
.352
1.083
,21182
.687
.271
91
-------
(14)
In this relationship 'W is the pulse width in the stream (miles) defined
as the duration of the square pulse (days) times the stream velocity (miles
per day). The single-valued relationship between the ratio and N is shown
in Figure 17.
The time between entry of the leading edge of the square pollution
pulse and the time at which the peak DO deficit concentration occurs (t*)
is also of interest in some problems. By analysis of the same computed
data, shown in Table 14, it was found that the non-dimensional time (T*)
between entry and occurrence of the DO deficit peak when expressed as
follows also plotted as a single-valued function of the non-dimensional
group N.
T* = Et*/W2 (15)
The single valued relationship between T* and N is shown in Figure 18.
It can be seen from Figure 17 that the effect of dispersion in the
stream is to reduce the peak DO deficit concentration significantly. The
time at which the peak DO deficit occurs is shown by Figure 18 to be much
shorter than the corresponding time to the peak DO deficit as computed with
the Streeter-Phelps model.
SUMMARY AND CONCLUSIONS
A closed-form expression for the dissolved oxygen deficit response in
the stream to an impulsive BOD load can be used with the superposition
principle to estimate the response to any known stormwater overflow event.
The closed-form solution agrees well with corresponding numerical integra-
tion solution. Non-dimensional groupings have been developed to summarize
computed values for peak dissolved oxygen deficit and time of occurrence
for any square stormwater pollution pulse.
92
-------
1 0 9
8
6
5
4
3
2
01 in
. 1 10
6
5
4
3
2
0.01 '
0
r
1
—
>
e:
>
(-
f-
(-
r-
r
C.
^
4
C
c
r
4-
U
"•^
>
c-
5
C
—
j
"I T """ "
'
345
::iK
T" ' ~
|
6
S
1
7
i
i
r^
/
i
s <
?r
T
y
"
j i
.
s
r'
0
0
^
-
va - -- -
*, "
2
ffl|||[ttt|lj
* (•
C >^
__ S_ 4
* '
' .
S
3 4
j
— -e
5 6
7
1
8
S
<
)
r
\
10
).
-4
0
V
V
T
^r
...
>
2
Figure 17. Correction factor for peak dissolved oxygen deficit to be applied
to the value computed with the Streeter-Phelps model.
-------
0.01
0.1
Figure 18. Relationship for finding the time (t*) to the peak dissolved
oxygen deficit from the non-dimensional groupina N.
94,
-------
REFERENCES FOR SECTION 6
1. Bennett, James P., "Convolution Approach to the Solution for Dissolved
Oxygen Balance Equation in a Stream", Water Resources Research, Vol.7,
No. 3, (1971) pp 580-590.
2. Bunce, Ronald E. and Hetling, Leo J., "A Steady State Segmented Estuary
Model" Tech. Paper No. 11, FWPCA, U.S. Dept. Interior, NTIS
PB-217-469 (1969)
3. Churchill, Ruel V., "Operational Mathematics" McGraw Hill Book Co.
(1958).
4. Di Toro, Dominic M., "Recurrence Relations for First Order Sequential
Reactions in Natural Waters", Water Resources Research, Vol. 8, No. 1,
(1972) pp 50-57.
5. Heaney, James P. et al, "Nationwide Evaluation of Combined Sewer
Overflows and Urban Stormwater Discharges" EPA-600/2-77-064 (1977).
6. Smith, R. and Eilers, R. G., "Documentation for SWOPS: A computerized
Stream Model to Study the Effect of Dispersion on Storm Overflow
Events", EPA Internal Memorandum (June, 1977).
7. Streeter, H. W. and Phelps, E. B., "A Study of the pollution and
Natural Purification of the Ohio River III - Factors Concerned in
the Phenomena of Oxidation and Reaeration", Public Health Bulletin
No. 146 (1925).
8. Thomann, Robert V., "Systems Analysis and Water Quality Management",
McGraw-Hill Book Co. 1972 pp 140-144.
95
-------
APPENDIX FOR SECTION 6
The following symbols are used in this section:
t = time, days
AT = incremental time for segmented storm overflow load, days
v = stream velocity, miles/day (km/day)
2
A = stream cross-sectional area, sq miles (km )
X = distance measured downstream from a point on the stream
bank, miles (km)
x - distance measured downstream from a point in the stream, miles (km)
M = mass of impulsive BOD load Ibs (kg)
L = concentration of BOD in the stream, mg/1
L = initial BOD concentration in the stream, mg/1
o
L' = concentration of BOD in the stream, Ib/cu mile (kg/km )
L. = concentration of BOD in the I'th segment of the storm overflow, mg/1
Q. = volume flow of I'th segment of the storm overflow, cfs (m /sec)
Q = volume flow of receiving stream downstream of the storm
overflow, cfs (m /sec)
D = concentration of dissolved oxygen deficit in the stream, mg/1
DQ = initial concentration of dissolved oxygen deficit in the stream, mg/1
o
E - dispersion coefficient, sq miles/day (km /day)
Kr = rate of BOD removal by oxidation and sedimentation, days'
K, = rate of deoxygenation in the stream, days"
K - rate of stream reaeration, days'
H = the constant 3.1416
96
-------
APPENDIX (continued)
The following symbols are used in the program SWOCFS. Symbols
corresponding to those used in the paper are KR = K , KA = K , EC = E,
r a
VEL = v, Q = Q and DT = AT in hours. Symbols not used in the paper are
defined as follows:
DN = number of time increments in storm overflow
TM = number of time increments over which solution is required
QIN(I) = volume flow of storm overflow at center of I'th time increment,
cfs (m /sec)
LOIN(I) = BOD concentration in overflow at center of I'th time increment,
mg/1
DO(J) = computed dissolved oxygen deficit at J'th stream segment, mg/1
LO(J) = computed BOD concentration at J'th stream segment, mg/1
97
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
1. REPORT NO.
EPA-600/2-78-148
3. RECIPIENT'S ACCESSION>NO.
4. TITLE AND SUBTITLE
STREAM MODELS FOR CALCULATING POLLUTIONAL EFFECTS
OF STORMWATER RUNOFF
5. REPORT DATE
August 1978 (Issuing Date)
6. PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
8. PERFORMING ORGANIZATION REPORT NO
Robert Smith and Richard 6. Eilers
'9. PERFORMING ORGANIZATION NAME AND ADDRESS
Municipal Environmental Research Laboratory--Cin., OH
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, Ohio 45268
10. PROGRAM ELEMENT NO.
1BC611
11. CONTRACT/GRANT NO.
12. SPONSORING AGENCY NAME AND ADDRESS
Same as above
13. TYPE OF RE PORT AND PERIOD COVERED
Inhouse (Jan-Dec 1977)
14. SPONSORING AGENCY CODE
EPA/600/14
15. SUPPLEMENTARY NOTES
Project Officer: Richard G. Eilers 513/684-7618
16. ABSTRACT
Three related studies are described which provide the means to quantify
the pollutional and hydraulic effects on flowing streams caused by stormwater
runoff. Mathematical stream models were developed to simulate the biological,
physical, chemical, and hydraulic reactions which occur in a stream. Relation-
ships take the form of differential equations with the two independent variables
of time and distance. The differential equations can be solved directly by
means of calculus or by digital computer using numerical methods. The solution
would be the concentration of species of pollutional interest, such as BOD and
dissolved oxygen, within the stream as a function of distance and time. The
solution can be steady-state or transient. There is a sufficient amount of
information presently available for finding steady-state solutions. However,
when the pollutional loads and/or the initial conditions for the flowing stream
vary with time, the problem becomes much more difficult and the technology for
handling the transient situation has not been adequately developed. The purpose
of this report is to show how the solution can be found for the case where the
pollution loading is a transient, especially as it applies to the stormwater
overflow.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS
c. COSATI Field/Group
*Surface water runoff
*Mathematical models
*Stream pollution
*Stream flow
*Water pollution
Stormwater overflow
Streeter-Phelps
Dissolved oxygen
deficit
Stormwater hydrograph
13 B
3I5TRIBUTION STATEMENT
Release to public
EPA Form 2220-1 (9-73)
19. SECURITY CLASS (This Report)
Unclassified
21. NO. OF PAGES
108
20. SECURITY CLASS (This page)
Unclassified
22. PRICE
98
6U.S GOVERNMENT PRINTING OFFICE. 1978— 757-140/1424
------- |