A METHODOLOGY FOR SELECTING AMONG
WATER QUALITY ALTERNATIVES
GILBERT S.
By
NICOLSON. JR.
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1969
ACKNOWLEDCMENTS
The author wishes to acknowledge the support furnished by the
U. S. Army Corps of Engineers. This dissertation was prepared under the
Corps' Graduate Fellowship Program, and any reproduction, in whole or in
part, for the purpose of the United States Government is authorized.
The author is very grateful to his chairman, Dr. E. E. Pyatt, who
was a source of encouragement and inspiration throughout this study. The
help and guidance of Dr. R. L. Patterson, Dr. W. C. Huber, and Dr. 1. G.
Keig during the course of this work were of inestimable value. The
suggestions of Dr. J. P. Heaney for improving the manuscript were appre
ciated. The author also acknowledges the help of Dr. J. C. Schaake and
Dr. D. H. Moreau, whose ideas and insights provided the foundation for
this study.
The author is grateful to the personnel of the Corps of Engineers:
Mr. Herb Rogers of the South Atlantic Division Office, Mr. Angelo Tabita
of the Jacksonville District Office, and Messrs. Frank Posey and Joe DeWitt
of the Savannah District Office. The assistance afforded the author during
the Fellowship Program by the personnel of the Project Planning Branch of
the Jacksonville District Office especially is appreciated.
This dissertation never could have been completed without the encour
agement and understanding of the author's wife. Her faith has been an
inspiration throughout the author's graduate career and for this he is
eternally grateful.
The author acknowledges the aid of Mrs. Linda Smith in the typing and
preparation of the manuscript. The typing of Mrs. Jean Weidner also is
greatly appreciated. To Mr. Bob Jandrucko special thanks are given for
his preparation of the drawings.
The help of the author's colleagues was gi itly appreciated. Their
many suggestions and encouragement added to the successfull completion of
this study.
The author acknowledges the financial support of a Public Health
Service Traineeship during his early graduate studies. Computer time was
furnished under FUPCA Grant No. WP01050 and by the University of Florida
Computing Center.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS .............................................. ii
LIST OF TABLES ................................................ vi
LIST OF FIGURES ............................................... vii
ABSTRACT ...................................................... viii
CHAPTER
I. INTRODUCTION AND OBJECTIVES ......................... 1
Introduction ........................................ 1
Historical Background ............................... 1
The Parameter of Major Interest ..................... 3
Water Quality Alternatives ............. ..... ... .. 5
Objectives .......................................... 7
II. RESERVOIR EFFECTS ON STREAM QUALITY ................. 8
Effects on Temperature .............................. 8
Heat Budget ........................................ 10
Internal Transfer Mechanisms .......... ...... ........ 14
The Heat Transfer Equation ........ .... .............. 16
Thermal Stratification of Lakes ..................... 18
Density Currents ................................... .. 24
Reservoir Prediction Models .......... .............. 28
Reservoir Temperature Model Used in This Study ...... 34
Dissolved Oxygen ...................................... .. 40
Improving Release Water Quality ..................... 46
Method Used in This Study .................. .......... 48
III. STREAM TDEPERATURE MAD DISSOLVED OXYGEN MODELS ...... 51
Stream Temperature ................................... 51
Dissolved Oxygen .................................. 56
Determination of Allowable Waste Loading ............ 69
IV. THE TRADEOFF BETWEEN STREAM TEMPERATURE AND
STREAM FLOW ........................................ 72
The Effects of Temperature on Waste
Assimilative Capacity ............................. 72
The Effects of Stream Flow on Waste Assimilative
Capacity ........................................... 75
The Stream Flow Stream Temperature Tradeoff ......... 76
Page
V. DEVELOPMENT OF METHODOLOGY ............................... 79
Description of the Problem ............................. 79
Water Quality Optimization Models ...................... 82
Dynamic Programming .................................. .. 88
VI. DATA REQUIREMENTS ..................................... 110
VII. TEST CASE .............................................. 114
Description of River Basin ................................ 114
Model Verification ..................................... 121
Waste Dischrage Model .................................. 124
Discharge Rates ........................................ 131
Alternatives Examined and Results ..................... 132
Summary of Test Case ................................... 139
VIII. SUMMARY ................................................. 140
APPENDICES
A. PROBABILISTIC DEFICIT FUNCTION ........................ 145
B. PROGRAM LISTING ........................................ 150
C. GLOSSARY OF PROGRAMMING TERMS ......................... 176
REFERENCES ................................ ..................... 180
BIOGRAPHICAL SKETCH ............................................. 185
LIST OF TABLES
Table Page
1. Proportionality Coefficients and Intercepts for the
Temperature, Vapor Pressure Approximation ................ 38
2. Sensivity Analysis ........................................... 67
3. Cost Matrix for Dynamic Programming Examole ................ 91
4. Evaluation of Equation (78) for N Equal T:.o in Example ..... 93
5. Evaluation of EQuation (78) for I Equal Three in Example ... 94
6. Evaluation of Equation (78) for II Equal Four in Example .... 95
7. Reservoir Dissolved Oxygen Constants and Epilimnion
Depth Variation ............................................. 117
8. River Hydraulic Data ........................................ 119
9. Distribution of Ungaged Runoff ............................. 120
10. Reservoir Release Temperatures ............................. 123
11. Base Flo .................................................... 126
12. Linear Regression Results .............. .................... 130
13. Data Listing for Trial Runs ................................ 133
14. Minimum Total Costs for Various Alternatives ............... 134
15. Optimum Decisions for Run A ................................ 135
16. Optimum Decisions for Run C ................................ 136
17. Optimum Decisions for Runs F and C ......................... 137
18. Results of Runs for Different 1l1 Values .................... 138
LIST OF FIGURES
Figure Page
1. Stream Temperature vs. Time .................................. 9
2. Coordinate S stem ............................................ 17
3. Annual Temperature Cycle for Lakes in Temperate Zones ...... 19
4. Eddy Diffusivity Coefficient and Temperature vs. Depth ..... 22
5. Temperature Profiles as Predicted by Dake and 11arleman ..... 25
6. Inflowu and Outflow. Velocity Profiles After Huber and
Harl man ................................................... 30
7. Assumed Velocity Profile .................................... 36
8a. Dissolved Oxy:gen Variation wich Time and Depth ............. 43
8b. Temperature Variation 'i th Time and Depth .................. 44
9. Macrix for Representing Dissolved Oxygen Variation
Within the Reser'.'oir ...................................... .49
10. Dissolved Oxygen Sag Curve ................................. 59
11. Hypothetical Curves for Waste Assimilative Capacity
vs. Release Temperature and Flo ............................... 77
12. RiverReservoir System .................................... 80
13. Dynamic Progra inming Flo.' Chart ............................. 98
14. Savannah River Basin ...................................... 116
15. P.cervoir Surface Temperatures ............................. 122
16. Downstream Tnimperature Profiles for August 1957 ............ 125
17. Wlastc Assimilative Capacity vs. TemperaLure and Flo': ....... 127
18. Waste /.ssimilative Capacity vs. Equilibrium Temperature
and Dissolved Oxygen ...... ............................... 128
Abstract of Dissertation Presented to the Graduate Council in
Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy.
A METHODOLOGY FOR SELECTING ANONG
WATER QUALITY ALTERNATIVES
By
Gilbert S. Nicolson, Jr.
June, 1969
Chairman: E. E. Pyatt
Major Department: Environmental Engineering
The Water Quality Act of 1965 requires each state to establish
stream quality standards. The major cause of violation of the standards
is municipal and industrial waste discharges. To maintain the standards
most state enforcement agencies require a high degree of waste treat
ment. Such requirements entail a large expenditure of private and
public funds for maintaining the standards during lowflow periods.
These funds could be used more.efficiently if provisions were made for
the examination of alternative schemes.
In the case of a municipality downstream from a reservoir, several
alternatives exist both at the reservoir and the waste treatment plant.
One of the alternatives at the treatment plant is an auxiliary unit
which could be operated during lowflow periods. Alternatives at the
reservoir site include multilevel outlet structures, reaeration of
releases, and reservoir mixing. This dissertation presents a methodology
for evaluation of these alternatives. For each alternative considered,
optimal decisions are made which minimize the total waste treatment costs
of the auxiliary unit. Comparison of the costs of providing the
alternative enable the selection of the leastcost alternative.
viii
The concept of a water quality tradeoff is introduced. The
tradeoff exists in reservoirs which thermally stratify during the summer
months. Stratification permits selective withdrawal from the reservoir
and a substitution of water at one temperature and flow rate for water
at another temperature and flow rate while still maintaining water
quality standards.
The methodology utilizes a twodimensional dynamic programming
algorithm for the optimization technique. Several water quality models
are needed to determine the reservoir release temperatures, the down
stream temperature profiles, and the dissolved oxygen response to
waste loads. Detailed analyses of the various models are presented.
A stochastic dissolved oxygen formulation is presented which permits
the determination of waste loads that satisfy the stream standard with
a given probability.
A test case is presented which illustrates the use of the
methodology. Included in the test case is the verification of the
temperature models used in the study. The existence of a water quality
tradeoff is demonstrated and the costs for several alternatives are
determined. The methodology is shown to be sensitive to the deoxygenation
coefficient, and it is concluded that a good estimate of the coefficient
is needed for obtainingreliable cost comparison.
I. INTRODUCTION A[TD OBJECTIVES
Introduction
During the rapid technological expansion of the twentieth century,
man's capacity to produce waste has caused deterioration of his environ
ment. The vast quantities of waste and the sometimes acute health
hazards and unpleasant esthetic conditions have produced a public
awareness of the need to protect the environment, not only for the
present generation but also for posterity.
One of the most important natural resources that must be preserved
for man's existence is water, which is used for recreation, navigation,
water supply, power production, maintenance of fish and wildlife,
irrigation, and assimilation of domestic and industrial waste. It is
this last use that is threatening the majority of the other uses and
that recent federal legislation is designed to control.
Historical Background
Federal interest in water was initially for the purpose of
navigation and flood control. Public Law 845 enacted by the Eightieth
Congress was the first federal legislation intended for pollution
abatement.I This law provided for monies to be spent for aiding in
design and construction of pollution abatement facilities. Since 1948,
several laws have been enacted, each with the general purpose of
providing more funds for abatement facilities. In 1963, Congress
directed that all new multipurpose reservoir projects should consider
additional storage for releases during periods of low flow.2 This
was the first step toward water quality control in reservoir projects.
However, Congress specified that such releases should not be made in
lieu of waste treatment.
With the Water Quality Act of 1965, Congress transferred the
Division of Water Supply and Pollution Control within the U. S. Public
Health Service to create the Federal Water Pollution Control Administra
tion (FW.PCA) in the Department of the Interior. The Act required that
each state set water quality standards subject to approval by the F.WPCA.3
The water quality standards should be designed to "enhance the quality
and value of our national resources... ."3 In the "Guidelines for
Establishing Water Quality Standards for Interstate Waters"4 the FUWPCA
stated that "no standard will be approved which allows any wastes
amenable to treatment or control to be discharged into any interstate
water without treatment or control regardless of the water quality
criteria andwater uses adopted."
The Water Quality Act of 1965, therefore, provides for the
establishment of quality standards for the preservation and maintenance
of one sector of man's environment. The maintenance of water quality
standards implies a vast outlay of public funds since all waterborne
waste must be subject to treatment and control.5 The Water Resource
Engineer must undertake the responsibility of providing the best
engineering plans for the expenditure of these funds.
Since water quality standards are established for an entire basin,
plans for meeting the standards should apply to the basin or to
independent parts of the basin. Water resource projects usually have been
designed according to the theory of efficient allocation of resources.
In the past, the efficiency criterion provided for the maximization of
net benefits. The engineer would list the purposes for which the project
was to be built. For each purpose he would perform a benefitcost
analysis. The purpose would be included in the project if the benefits
accruing to this purpose exceeded the cost of providing it. The scale
of each purpose would be increased until the benefits equaled the costs;
therefore, net benefits would be maximized.
With the establishment of water quality standards, the theory of
efficient allocation still can be satisfied. However, the constraint
is modified from requiring that benefits exceed costs to satisfying the
water quality standards. By changing the constraint, it is possible to
alter the efficiency criterion to selecting the leastcost method which
achieves the water quality objectives. This criterion does not require
that the standards be established under the principle of maximizing net
benefits. One can view water quality control under the stream standards
as a problem of maximizing net benefits. If water quality control is
assumed a benefit, the leastcost alternative which maintains the
standards will by necessity maximize net benefits. In either case, the
major difficulty is determining the costs of the alternatives.
The Parameter of Major Interest
As mentioned earlier, it was the deterioration of one of man's most
vital resources coupled with lethargic state programs which precipitated
4
federal legislation. The stream standards are established for the
purpose of preventing further pollution and improving water quality in
badly polluted streams. Therefore, the parameters of primary concern
should be the water quality parameter which is sensitive to pollution
and for which suitable mathematical expressions exist. These expressions
should trace the response of the parameter to waste loads and provide
a prediction of the response as a function of the waste loading. The
water quality parameter frequently used for pollution investigations
is dissolved oxygen. There are many reasons for the selection of
dissolved oxygen: first, it can be accurately measured; second, the
aquatic biota and bacteria which assimilate the organic content of
the waste require dissolved oxygen to maintain their metabolic processes;
third, the level of dissolved oxygen can be a determining factor as to
the species of bioca that can survive in the stream; fourth, the
absence of dissolved oxygen is readily associated with severe pollution.
For these reasons, dissolved oxygen will be the primary water quality
parameter considered in this study. Since dissolved oxygen is an
important resource with respect to the stream's capacity to assimilate
organic or inorganic waste, any action taken to increase the concentra
tion of dissolved oxygen or the total oxygen resources would also increase
the stream's waste assimilative.capacity.
A stream's waste assimilative capacity is affected by the temperature
of the stream. The stream's dissolved oxygen saturation level is
markedly affected by temperature,and the rate at which organisms decompose
the waste also is influenced by temperature. Placing a reservoir on a
stream alters the stream temperature regime below the reservoir. Since
the reservoir is a point of control in the basin, it is possible to
change the stream's assimilative capacity simply by controlling the
temperature of the release water. (The effects of reservoirs on stream
temperature and methods for controlling the releases will be detailed in
a later chapter.) The major concern of the present study deals with the
importance of temperature and methods of improving waste assimilative
capacity through temperature control. Other methods of water quality
control are reservoir mixing, turbine reaeration to increase the
dissolved oxygen in reservoir releases, and instream reaeration.
Water Quality Alternatives
There are several methods employed either to increase the allowable
waste loading or to alleviate the detrimental effects of pollution. The
method quite often employed is reducing the waste discharge into the
stream, which can be accomplished by installing waste treatment facilities
or by altering industrial processes to reduce waste. In the case of waste
treatment facilities, considerable capital investment by private and
public institutions is necessary to achieve high levels of waste treatment.
These monies should not be spent without some consideration to alternative
means of meeting water quality standards.
A second method of alleviating pollution problems is to provide
adequate dilution of wastes. TI most cases, water quality deteriora
tion is most severe during the summer and fall months when the natural
stream flow is lowest. Storing water in reservoirs for release during
lowflow periods has become a method of water quality control.
Bramhall and Mills7 criticized the use of lowflow augmentation as an
alternative method of meeting stream standards when compared to high
levels of waste treatment. However, in their analysis the authors
considered reservoir storage only for use as lowflow augmentation,
which would introduce a bias in favor of waste treatment. Young
pointed out the compatibility of flood control storage and lowflow
augmentation. Flow augmentation seems to be technically feasible as
a means of maintaining water quality objectives and should be con
sidered as an alternative to waste treatment.
A third method of water quality control is waste storage during
times when stream conditions are such that the waste discharge would
cause violation of stream standards. Loucks9 illustrated the use
of waste storage in the establishment of probabilistic stream standards.
Loucks did not consider any waste treatment in the storage facility,
but Moreau10 presented an analysis for the pulp and paper industry of
a treatment system composed of steadystate and stochastic units. The
effluent from the steadystate unit always received a constant degree
of treatment. When the receiving stream was unable to accept the con
stant effluent, a portion of the waste went to the stochastic unit
where it received sufficient treatment so that water quality standards
were not violated. Therefore, an alternative to a high degree of
treatment the year around would be the storing or treating part of the
waste in a variable waste treatment unit during low flow periods.
Objectives
In planning to meet water quality standards, all alternative
schemes and combinations should be considered in attempting to find
the leastcost alternative. A methodology or framework is needed
for comparing water quality alternatives in order to make the best
decision. The objective of this disseration is to develop a
methodology for selecting from among water quality control alternatives
for a reservoir river system. The methodology will employ mathematical
expressions for tracing the temperature and dissolved oxygen profiles,
an inventory system for reservoir operation and a dynamic programming
technique for selecting optimal reservoir releases.
In Chapter II the effects of a reservoir on temperature and
dissolved oxygen along with the models used to predict them will be
discussed. Chapter III details the stream models needed to calculate
downstream temperature and dissolved oxygen. A discussion of the
tradeoff between reservoir release temperatures and flows will be
presented in Chapter IV, and the development of the methodology to
make use of this tradeoff will be formulated in Chapter V. Data
requirements for the proposed methodology will be discussed in Chapter
VI. Chapter VII will consist of the test case illustrating the use of
the methodology and Chapter VIII will contain the summary and con
clusions. Flow charts of programs and a glossary of terms will be
placed in the Appendices.
II. RESERVOIR EFFECTS ON STREAM QUALITY
To develop alternative plans for meeting water quality standards
in a reservoirriver system, one must first understand the changes
imposed on the river with the introduction of a reservoir. This
chapter presents a discussion of the two quality parameters of partic
ular interest in this study: temperature and dissolved oxygen. The
discussion will include reservoir effects, methods used to predict or
describe the effects, and the model selected for this study.
Effects on Temperature
Ward11 has studied stream temperatures in regulated and unregulated
rivers and has proposed an expression for representing the temperature
as a function of time:
T = a [sin (bt + c)] + T ......................... (1)
where T is the stream temperature in degrees Fahrenheit; a is the
amplitude in degrees Fahrenheit;
360
b = 0= .987 degrees per day;
365 days
t is the number of days since October 1: c is the phase coefficient in
degrees; and T is the average annual temperature in degrees Fahrenheit.
The coefficients a, b and c are determined from a leastsquares analysis
of the scream temperature data. Ward calculated these coefficients
for several rivers in Arkansas. Figure 1 illustrates the effect of
9
I
<,
z W w
o0
\w 4
+ 
I C
r
z I '
0 0
Co \ o s
Jo '3Nn1V8J]dLN31
! / l t
N / S =

10
placing the Bull Shoals Reservoir on the White River: the average
temperature T was lower; the annual variation 2a decreased; and the
phase coefficient Icl increased, causing approximately a fortyday lag
in the sine curve. This represents the change in the time variation
of stream temperature for one specific reservoir, but, depending on
reservoir operation and outlet elevations, this is a typical before
andafter example. Churchill2, 13 has studied several cases of the
temperature change caused by reservoirs built by the Tennessee Valley
Authority. In each case, the effects were similar to those in Figure 1.
Before the models used to describe the changes in reservoir
temperatures can be discussed, it is necessary to understand the
relationship between the heat budget of a body of water and the
mechanisms of heat transfer in the water and the role each plays
during the annual cycle of reservoir temperatures. Most methods for
predicting water temperature depend upon the determination of the
water's heat budget, which is a summation of the available sources
and sinks of heat acting on the water. Calculations for the heat
budget are made to obtain the net rate of heat transfer at the
surface. The development of the equations for the heat budget will
follow the work of Edinger and Geyer.14
Heat Budget
There are seven major mechanisms of heat transfer at the water surface.
These are incoming shortwave solar radiation, H ; incoming longwave radia
tion H ; reflected shortwave and longwave radiation, Hsr and Har; long
a
wave back radiation, Hbr; heat conduction, Hc; and heat loss by evaporation,
H Following is a brief discussion of each of these terms and their
e
relative magnitudes. (For a more detailed discussion, the reader is
referred to Edinger and Geyerl4 or to the Report of the 196061
Advanced Seminar of the Johns Hopkins University, Department of
Sanitary Engineering and Water Resources.15)
Absorbed Radiation
Shortwave solar radiation reaching any point on the earth depends
upon the latitude, cloud cover, sun's altitude, and time of year.
Empirical formulas are sometimes used for estimating shortwave radia
tion, but it is measured more easily than it is computed. The U.S.
Weather Bureau records shortwave radiation and reports the measurements
in Climatological Data  National Summary. Monthly values for solar
radiation in the northern latitudes range from 400 2800 Btu per square
foot per day. The solar radiation incident on a body of water is
absorbed in the first ten to twenty feet, depending upon the optical
properties of the water.1
Longwave atmospheric radiation is a function of cloud cover,
air temperature and air vapor pressure. Raphael presents the follow
ing formula for longwave radiation:17
Ha = o o(Ta + 460)4 ........................... (2)
where Bo is an atmospheric radiation factor depending on cloud cover
and air vapor pressure; o is the StefanBoltzman constant (4.15 x 108
Btu per square foot per degrees Rankine per day); and Ta is the air
temperature in degrees Fahrenheit. Longwave radiation ranges from
2400 3200 Btu per square foot per day.
Reflected long and shorta'.'ave radiation are calculated as fractions
of incident long and shortwa:e radiation. Reflected radiations can
range from 100 300 Btu per square foot per day.
H = .03 H .................
ar a
H = .05 H .................
.................... (3)
.................... (4)
sr s
Edinger and Geyer combine H H H and H into one term called the
s a ar sr
absorbed radiation, HR.
H = H + H H H ....................... (5)
R s a sr ar
The absorbed radiation is independent of the water temperature and is,
therefore, conveniently calculated from meteorological data. Absorbed
radiation ranges from 2500 5700 Btu per square foot per day.
Heat Losses
Longwave
depends on the
radiation is:
back radiation is emitted from the body of water and
temperature at the water surface. The equation for back
H = y (T + 460)4 ............................
Br w s
where o is the StefanBoltzman constant; Y, is the emissivity of water
(0.97); and Ts is the surface water temperature in degrees Fahrenheit.
Back radiation is a substantial heat sink and ranges from 2400 3600
Btu per square foot per day.
Another large heat sink is evaporation, which is a function of
wind speed and the vapor pressure gradient between the air and the
water surface. Several empirical formulas are used for calculation of
heat loss due co the latent heat of vaporization. The general form
of the equation is:
H = (a + bW) (e e ) ......................
where a and b are constants depending upon the evaporation formula
employed; W is wind speed in miles per hour, mph; e is the saturation
5
vapor pressure of water in millimeters of mercury, mm Hg, determined
from the water surface temperature, T ; and e is the airvapor pressure
s a
in mm Hg. For the Lake Hefner Formula, a and b are 0 and 11.4, respec
tively; while in the Meyer Equation, for monthly evaporation, they are
73 and 7.3,respectively.14 Evaporative heat loss ranges from 2000 
8000 Btu per square foot per day.
Heat conduction is positive or negative depending upon whether
the difference between the air temperature and the water temperature
is positive or negative. The Bowen ratio is used to relate heat
conduction to evaporative heat loss.17
B= H
B = c ....................................... (8)
He
B is the Bowen ratio and is determined from the equation
C (T T) P
B= s a
B = C s a ................ ..... (9)
(e ea) 760
where T T e and e are as defined in preceding equations; P is
s a s a
the atmospheric pressure in mm Hg; and C is a coefficient determined
from experiments to be approximately 0.26. Substituting Equations (9)
and (7) into Equation (8) and assuming P is 760 mm Hg yields:
H = 0.26 (a + bW) (T T ) .................. (10)
c s a
Heat conduction can range from a minus 320 to a positive 400 Btu per
square foot per day.
Net Heat Transfer
The net rate of heat transfer, _H, is the algebraic sum of the
absorbed radiation, back radiation, evaporative heat loss, and
conductive heat losses. Assuming conductive heac transfer as positive
when the air temperature is greater than the water temperature, the
equation for net heat transfer is:
AH = H (Hb + H + H ) ....................... (11)
R br e c
The following discussion of Equation (11) is abstracted from Edinger and
Geyer. When the absorbed radiation, HR, is greater than the rate at which
heat is lost by a body of water through evaporation, conduction, and
back radiation, the net rate atwhich heat enters through the water surface
is positive (AH >0) and heat is added to the water. Heat leaves the water
(AH > 0) when absorbed radiation is less than the rate at which heat is
lost. When the net rate of heat transfer is zero (AH = 0), a special con
dition arises, which leads to the definition of the equilibrium temperature.
(The role of the equilibrium temperature in reservoir temperature predic
tion will be discussed later in this chapter.) Substituting Equations (6),
(7) and (10) into Equation (11) yields:
AH = H y (T + 460)4 (a + bW) (e e ) .26 (a + bW) (T T )
R s s a s a
............ (12)
Equation (12) is referred to as the "basic" heat transfer equation
because it can be derived without any assumptions as to the type of water
being studied or specification of an evaporation formula.14 Equation (12)
cannot be used to predict'water surface temperatures unless some other
method is used for obtaining AH. However, Equation (12) will be used to
derive the predictive equation used in this study.
Internal Transfer Mechanisms
The discussion of the heat budget was limited to the heat transfer
mechanisms at the water surface. The heat transfer mechanisms which
distribute the heat downward into the reservoir are the core of the
predictive reservoir temperature models. There are two modes of heat
transfer, other than absorption, in a body of water which distribute
the heat vertically.
One mode of heat transfer in a body of water is advection.
Advective heat transfer is the result of heat transferred in a flowing
fluid by macroscopic particles. The source of advected heat in
reservoirs is inflow, either surface or underground, outflows and
rainfall. In most reservoir studies, advected heat from inflow and
outflow is an important heat transfer mechanism.
The second mode of heat transfer is diffusion. Diffusion in water
can be either molecular diffusion or eddy diffusion. The heat flux
per unit area due to molecular diffusion is given by the product of
the temperature gradient and the molecular diffusivity coefficient a.
The diffusivity coefficient in water is very small in comparison to
eddy diffusivity (a = .00144 cm2/sec) and is considered independent of
water temperature.
Heat flux by eddy diffusion is the result of eddy diffusivity
caused by windinduced currents, convectional currents resulting
from surface waters cooling and sinking to lower depths, inflow
currents from tributaries to the reservoir, or outflow currents.
Wunderlich and Elder19 present a detailed discussion of what they term
turbulent diffusion mechanisms and their importance in reservoirs.
They indicate that convective currents are the most important source of
eddy diffusivity in the surface layer of reservoirs and that inflow
and withdrawal currents are important in the lower layers. The heat
flux per unit area by eddy diffusion is the product of the thermal
gradient and the eddy diffusivity coefficient (E). The eddy diffusivity
coefficient is not constant and varies from point to point in the
reservoir and from time to time. The eddy diffusivity coefficient is
large in comparison to the molecular diffusivity coefficient (E ranges
from .01 to 10 cm2/sec).
The Heat Transfer Equation
Basic mathematical expressions representing the modes of heat
transfer can be combined into a differential equation for the rate of
heat exchange within the water. The equation for heat transfer in a
fluid will be derived making use of the principle of conservation of
heat. This principle states that the rate of change of heat within
the body (stored heat) is equal to the difference between the rate of
heat inflow and rate of heat outflows by the mass transport mechanisms
plus the rate of external heat added or lost. In equation form this is:
[rate of change of stored heat] = [rate of heat inflow by mass
transport] [rate of heat outflow by mass transport] + [rate
of external heat added] .................................. (13)
The threedimensional coordinate system and unit cube shown in
Figure 2are used in the derivation. The heat content of any body is
given as
oc T
where 0 is the density of the substance in pounds per cubic foot: c
P
is the specific heat of the substance in Btu per pound per degree
Fahrenheit; and T is the temperature of the substance in degrees
Fahrenheit. For this derivation the substance will be water, in
/x
Y Z
FIGURE 2. COORDINATE SYSTEM
which the density and specific heac can be considered constants. The
rate of change of the heat content of the cube is
Pcp T
where t is time.
The heat advected across the xy plane in the positive zdirection
is:
Ocp T OPp vx T 2 pc vx Sx
z = z+pz z = z
where v is velocity in the xdisection in feet per day.
The heat diffused across the xy plane in the positive z direction
by molecular and eddy diffusion is:
P dx P x ax xP X x' x
z = z+tz z = z
where E is the eddy diffusivity coefficient in the xdirection in
square feet per day.
Let the external heat term be represented by S in Btu per
cubic feet per day. Taking the sum of the advected and diffused terms
in the x, y and zdirections, substituting the appropriate terms into
Equation (13), and dividing by c the final expression becomes:
aT ;3 T ; aT 3 aT
TE x[ (Ex + a) x] + [E + c) y + a) 
[vx aT + T + v T] + S ........................ (14)
ax y ay z az pc
p
Equation (14) represents the equation for temperature distribution
within a body of water. It entails the basic principles of heat transfer
used by other investigations, however, it is virtually impossible to
solve Equation (14) without making several simplifying assumptions.
Thermal Stratification of Lakes
The preceding discussions on the heat budget and heat transfer
mechanisms are germane to all studies of water temperature and models
developed for predicting water temperature. The annual variation of
lake temperature with discussion of the important heat transfer
mechanisms will now be presented.
Initial discussion on thermal stratification will be limited to
lakes since it is the most elementary body of water found in nature.
Figure 3 shows an annual cycle of temperature profiles for a lake.
Most lakes in the temperate zone approach isothermal conditions
3800100
A'ln.P
AV I'
La
A UVNr
B3QW3AOIN
Hld3G
w
z
0
N
w
o I
CC
c.
o
w
L.
O
0
o w
0
0 Cw
o
0 0
0
M w
w cr
F J
0
0 <
sometime around March 21. At this time, the slightest wind will
provide enough energy to mix the lake so that it is of constant
quality throughout. The eddy diffusivitycoefficient is constant from
top to bottom. As spring begins, shortwave solar radiation is
absorbed in the surface layers, warming the upper strata. At first,
wind on the lake will mix warmer surface water with the lower water;
but the rate of solar absorption soon exceeds the capacity of the
wind for mixing, and the surface water becomes farmer than the lower
water. During this early warming period, heat is transferred
vertically by the eddy diffusion resulting from wind currents. Since
warm water is less dense than colder water (above 40C), the warm
water will remain at the surface and form a layer of constant temper
ature called the epilimnion. The formation of a temperature gradient
and, hence, a density gradient creates a stable condition within the
lake. It would take considerable energy to mix the lake. Beneath
the epilimnion, the temperature profile decreases until it reaches
another layer of constant temperature. This layer is called the
hypolimnion and it contains the coldest water in the lake. Since
there is no temperature gradient in the hypolimnion, the water in it
is easily mixed and, therefore, the quality of the water is essentially
constant. The layer between the epilimnion and hypolimnion is called
the metalimnion and contains the thermocline (defined as the maximum
rate of decrease in temperature, 0T = 0). Once the lake stratifies,
3 z2
the effect of eddy diffusion is limited because of the stable condition
20
of the reservoir. Orlob discusses the change in the eddy diffusivity
with depth. Since the rate of heat e::change is greatest at the sur
face, the edd.' diffusivit.' coefficient would be expected.to be gratest
at this point. The coefficient would decrease somewhat in the
epilimnion, but would still be relatively large since the epilimnion
is completely mixed. The lake is most stable at the thermocline,
and the eddy coefficient would be expected to be lowest at this point.
It would then increase slightly to a second maximum in the hypolimnion.
Figure 4 shows the change of the eddy diffusivity coefficient with depth.
As summer approaches, the epilimnion absorbs more solar radiation
and the temperature of the surface layer rises, increasing the stabil
ity. The heat loss by evaporation and by back radiation prevents the
temperature of the epilimnion from exceeding the mean ambient air
temperature. As water evaporates at the surface, it loses its heat of
vaporization and thus a small layer of cooler water is formed at the
surface. During the summer evenings, when the air temperature falls
below the epilimnion temperature, more heat is transferred from the
water surface by evaporation, conduction, and back radiation than is
absorbed in the epilimnion. In both cases, an unstable condition
exists since the surface water is more dense than that below it. The
slightest wind will upset this condition and set up convective currents
which mix the cooler water into the epilimnion.
During the fall, the epilimnion begins to cool. The rapid
decrease in the air temperature causes large transfers of heat from
the epilimnion and thus sets up large convective currents. These
currents go deeper and deeper into the metalimnion as the temperature
gradient decreases and cause a deep layer of constant quality. The
lake continually cools until sometime in November when the lake becomes
isothermal and the fall "turnover" takes place.
V WATER SURFACE
EPIL
METAL
EDDY
DIFFUSIVITY
COEFFICIENT
/
LIMNION
IMNION THERMOCLIN
.TEMPERATURE
HYPOL IMNION
4. EDDY DIFFUSIVITY COEFFICIENT
TEMPERATURE vs. DEPTH
FIGURE
and
Depending upon the climate, the lake will remain essentially
isothermal until March or will stratify inversely if the temperature
of the water goes below 4C. Since water is most dense at about 4"C,
further cooling of the surface will cause the warmest water to be at
the bottom. When winter ends and solar radiation warms the colder,
upper layer to 4C, the lake will again be isothermal and the spring
"turnover" occurs. The annual cycle is completed at that time.
(Hutchinson presents a classification of lakes according to their
annual temperature properties and discusses the formation of additional
layers within lakes.21)
During the annual temperature cycle of lakes, the largest portion
of heat transferred is confined to the upper layer. Almost all the
shortwave radiation is absorbed in the epilimnion and transfer of
heat through longwave radiation, evaporation and conduction takes
place at the surface. The major modes of heat transfer are eddy
diffusion from wind currents and convective currents. Dake and
Harleman22 derived an equation similar to'Equation (15) for the
prediction of temperature profiles in lakes. Their equation is
(using the notation of this study):
T a2T
a S .............................. (15)
at a z2
p Cp
Since there are no velocity components in lakes, the advection
term is zero. Because of the difficulty of selecting the diffusivity
coefficient, E, a priori, Dake and Harleman assumed it to be zero in
their formulation. The source term, S, was a function of incident
solar radiation, (Ha Har) in (Etu per square foot per day) the
1
absorption coefficient of the fluid, n in centimeters and the
fraction of radiation absorbed at the water surface, 8. The source
term in equation form is:
S = n(l 61) nx
S(Ha Har) e ............... (16)
The linearity of Equation (16) permitted the superposition of
solutions for temperature distributions caused by surface solar
absorption, internal radiation absorption, and surface evaporative
loss. During the heating period, the evaporative loss at the surface
caused a negative thermal gradient in the upper layer. Ihen this
occurred, Dake and Harleman mixed the epilimnion to a uniform temperature.
This effect is shown in Figure 5. Even though eddy diffusivity was assumed
zero, the negative temperature gradient caused by evaporation and the
resultant mixing procedure could be interpreted as eddy diffusion
resulting from convective cooling.
The preceding discussion of the annual temperature cycle was limited
to lakes. Th'e major difference between lakes and reservoirs is that
advective heat transfer is not negligible in reservoirs whose inflows
and withdrawals can have considerable effect on the temperature distribu
tion. These effects are dependent upon several factors  the quantity
and time distribution of the inflow, the withdrawal rate, and the level
of the outlet structure in the reservoir.
Density Currents
A density current is a gravity flow of a fluid of a specific
density through fluid of another density. tlhen a reservoir is
,MIXED EPILIMNION
TEMPERATURE
URE PROFILE
5. TEMPERATURE PROFILES as PREDICTED
by DAKE and HARLEMAN
I
a.
w
0
FIGURE
thermally stratified, layers of different density will be present, and,
therefore, any flow of fluid will constitute a density current. Church
hill2 reported density currents in TVA reservoirs and their effects
on reservoir and stream quality. Cold water releases from deep reser
voirs were observed to underflow downstream reservoirsand to be dis
charged from lowlevel outlets with only a small detention time in the
reservoir. This was made possible because the incoming water was
colder than that in the reservoir and was, therefore, more dense.
Gravity pulled the heavier water into the lowest stratum in the
reservoir creating the underflow. Since outlets from the dam were
also at the level of thislow stratum, the underflow moved directly
through the reservoir.
Koberg and Ford2t, cited density currents at different depths
within the reservoir. These currents appeared to be caused by inflow
to the reservoir sinking below the surface and following the slope of
the reservoir until the layer of corresponding density was reached.
At this point the water spreads out horizontally causing a velocity
profile within that layer in the reservoir. This phenomenon was also
observed in a dye study on the Fontana Reservoir.1 The dye was mixed
into the inflow and monitored. The inflow was observed to seek a layer
of approximately the same temperature. The layer containing the dye
was found to be three degrees Fahrenheit higher than it was upon
entering the reservoir. This was caused by mixine with warm surface
water when the water entered the reservoir.
In view of these observations, most investigators assume that
water withdrawn from a reservoir under stratified conditions is
selectively taken from the layer at the elevation of the outlet. The
larger the density gradient near the outlet, the more confined the
withdrawal layer becomes. Inflow is assumed to be confined to the
layer in the reservoir of corresponding temperature.
With the background of the annual temperature cycle of lakes and
the concept of density currents, the change in the downstream temper
ature variation, as shown in Figure 1, can now be explained. If the
outlet structure is located in the hypolimnion, the colder water in
the reservoir will be released first. Incoming water will seek that
level in the reservoir corresponding to its temperature and will
remain in that layer until the water beneath is released. Thus,
incoming water is not released until it settles to the level of the
outlet. The overall effect is a prolonged release of cold water
from the reservoir. This has the combined effect of smoothing the
annual temperature variation of released water, decreasing the average
temperature of released water, and causing a shift in time of the
point of maximum temperature.1 The extent of these effects is dependent
upon the level of the outlet structure, the volume of vater above this
level, and the rate of inflow and outflow.
If the outlet level is deep in the reservoir, more hypolimnion
water must be released before the temperature will rise. The greater
the volume of water above the outlet, the longer it will take to release
the hypolimnion water. If the operating rule of the reservoir requires
a high discharge rate in the spring and early summer, hypolimnion water
will be discharged earlier and the phase shift of the regulated sine
curve in Figure 1 will be less. EB constructing multilevel outlets,
the reservoir operator could release water from any level in the reser
voir during the stratified period and thus select that temperature which
best meets downstream requirements. The methodology developed in this
dissertation is based on the assumption that water can be released from
any desired level in the reservoir.
Reservoir Prediction Models
Early models used to predict reservoir temperatures were quite
simple. Raphaell7 proposed a model for temperature prediction in
reservoirs that were well mixed during the summer months. He assumed
the reservoir to be completely mixed and added the heat budget to the
advected heat, in the form of inflow and outflow, to give the total
hear input. This total was then distributed throughout the reservoir
to produce an average reservoir temperature. Since most reservoirs
undergo thermal stratification during the heating period, Raphael's
model had only limited application.
Delay and Seaders2developed a procedure for predicting temperature
profiles in proposed reservoirs. Their method was based on the use
of temperature profiles from existing reservoirs of comparable size
and nearby locale. Again, a procedure of this nature would have
limited application.
Orlob and Selna26 developed a mathematical model which was used
to simulate reservoir thermal stratification. In their formulation
they combined molecular diffusivity and eddy diffusivity coefficients
into the "effective" diffusivicy coefficient. "Effective" diffusivity
can only be determined from thermal profile measurements made at the
reservoir. Since diffusion is a function of time and space, the
measurements must be made periodically during the heating cycle.
Therefore, considerable data must be obtained to determine the varia
tion of the "effective" diffusivity coefficient for a single year.
Since reservoir conditions may change from year to year because of
different sequences of reservoir inputs, the "effective" diffusivity
coefficient should be determined for several years in order to have
a satisfactory knowledge of its behavior for a particular reservoir.27
Orlob and Selna's model required considerable data input for applica
tion to a single reservoir. Their model could not be applied to
proposed reservoirs without considerable knowledge of the "effective"
diffusivity coefficient in several geographical regions.
Huber and Harleman28 have developed a mathematical model for
predicting thermal stratification in reservoirs. Their model requires
knowledge of the reservoir geometry, inflow volumes and temperatures,
meteorological conditions at the reservoir and outflow volumes. The
model is, to some degree, an extension of Dake and Harleman's work
on lakes as the eddy diffusivity is assumed to be zero. Included
in the model are consideration for mining of the inflow with reser
voir surface waters before diving into the reservoir as a density
current and a technique for handling the inflow and withdrawal
velocity profiles as Gaussian curves. Figure 6 illustrates these
velocity profiles.
Koh29 and Kao '0 have demonstrated in the laboratory that
withdrawals from a fluid under a density gradient do not come only
from the layer at the elevation of the outlet, but are a composite
PROF
FIGURE 6. INFLOW and OUTFLOW VELOCITY PROFILE
after HUBER and HARLEMAN
of the layers in the vicinity of the outlet. The velocity profile in
Koh's work uas very similar to that shown in Figure 6. The degree to
which the velocity profile extends into other layers is dependent upon
the density gradient at the outlet. If the gradient is large, the
discharge will be composed of a rather narrow swath of water from the
reservoir. However, when the gradient is small, the width of the layer
affected by the discharge increases.
Huber and Harleman combine the heat transfer equation with the
continuity equation and the equation for continuity of mass to derive:
n x
BT aT a 3T qi (Tin T) n(l Eo)o e
T +v (A )+ +
dt dz A 15z z A pc ..... (17)
where A is the crosssectional area of the reservoir and is a function
of depth; qi is the inflow rate per unit of vertical distance; Tin is
the temperature of the inflow; and all other terms are as defined earlier.
Equation (17) vas solved by finite difference techniques on a digital
computer. The equation yielded excellent laboratory agreement between
observed and predicted temperatures. The'model was tested on data
collected by the TVA at the Fontana Reservoir. The observed temper
atures at Fontana ranged from 7 18C. The model yielded very good
results until October when the predicted temperatures were about 2C
less than those observed. The model developed by Huber and Harleman
could be used to obtain predicted temperatures for proposed reservoirs.
Huber and Harleman's model was not chosen for this study due to the
amount of computer time necessary for an annual cycle.
A suitable mode'. should provide a satisfactory representation of
thermal stratification and the volume and temperature of water in the
three layers of the reservoir and should be easily adapted to com
putation on the digital computer.
Wunderlich and Elder19 have developed a graphical procedure for
determining reservoir outflow temperatures. Their technique is based
on several simplifying assumptions. First, the heat exchange between
the reservoir and the atmosphere occurs only in a tenfoot layer at
the surface of the reservoir. 'In the light of preceding discussion on
convective cooling, this assumption might be improved if the depth
were allowed to increase, as the heating season progresses, from, per
haps five feet in early April to twenty or twentyfive feet in July
and August. This layer would correspond to the epilimnion and would
have the same temperature throughout.
Second, the temperature of the tenfoot surface layer corre
sponds to the equilibrium temperature of the exposed body of water.
This is a convenient assumption as the equilibrium temperature can be
easily calculated from the heat budget, as will be explained in a later
section of this chapter.
Third, each water parcel entering the reservoir spreads into
a horizontal layer corresponding to its temperature. Inflows of
approximately equilibrium temperature enter just below the surface
layer. Field investigations show this to be approximately correct
and to simplify calculations, most investigators make this assumption.
The same assumption was made concerning outflows, i.e., the water
withdrawn from the reservoir comes from a layer between the top and
bottom of the intake.
Fourth, heat transfer by molecular and eddy diffusion is
neglected below the surface layer. Since molecular diffusion is
very slow and eddy diffusion is limited by the thermocline, this
assumption appears to be valid.
Wunderlich and Elder applied the graphical procedure to the Fon
tana Reservoir in Tennessee. The observed temperatures agreed very
well with the predicted temperatures except in the month of September
when the predicted temperatures were about 50F less than the observed
temperature.
Ross and MacDonald31 have applied a very similar model to
reservoir projects in Montana. The major differences were the
method of predicting surface water temperature and the assumed
thickness of the surface layer. Ross and MacDonald used an equation
similar to Equation (12) for the net rate of heat transfer at the
surface. Assuming that the epilimnion was completely mixed, they
expressed the time rate of change of epilimnion temperature as:
dTw (c pm/t)(Ti Tw)
dt = H + 1 i ) ... ............... (18)
c p(mw + Mi)
where T, is the temperature of the epilimnion; AH is the net rate of
heat transfer at the surface; A is the reservoir surface area; Cp is
the specific heat of water; mi is the change in mass of the epilimnion
during the time interval dt; mw is the mass of the epilimnion; and
T. is the temperature of the epilimnion at the beginning of the time
1
interval. Data averaged over a sevenday period were used to calculate
AH. The epilimnion temperature computed for the end of one time period
was used as the water temperature for calculations in the following
period. This technique for calculating the surface temperature
yielded results which agreed very . ell with observed temperatures.
Ross and MacDonald selected a surface layer depth of thirty
feet for Hungry Horse Reservoir. This value was selected because of
previous thermal observations at the site. They also assumed that
infloing water sought that level in the reservoir corresponding to its
temperature and that withdrawals were made from the layer at the outlet.
Their procedure yielded good agreement between observed and predicted
temperature profiles.
These last two temperature predictive techniques lend themsel.'es
;' to digital computation with a minimum of machine time. The model used
for predicting reservoir temperatures in this study, will now be
. presented employing features from each of the t.o techniques.
Reservoir Temperature Model Used in This Study
The purpose of the reservoir model is to provide an inventory of
the volume and temperature of the reservoir resulting from various
operating schemes. To fulfill this purpose, it is necessary to devise
a method for keeping account of the inflo'.. and outflo'. from the various
layers and the temperature within these layers.
Since reservoirs stratify into three layers, an inventory method
was selected for adding inflows and subtracting outflows from three
independent layers. Each layer was considered to have a constant
temperature and to be completely mi:.ed. Withdrawals from a particular
layer during any time interval were assumed to have the quality
parameter predicted for the layer at that time. Since the epilimnion
and hypolimnion have approximately constant temperature profiles, and
thus a small density gradient, releases made from these layers will
create constant velocity profiles as shown in Figure 7.
It seems logical to consider the hypolimnion and epilimnion as
being completely mixed. However, the metalimnion has a varying
temperature profile and thus a steep density gradient. Releases made
under a density gradient have a velocity profile as described by Koh.29
Withdrawals from the metalimnion will be a composite of the quality
parameters in that layer and will be approximately the average of
those parameters. Therefore, the release parameters from the metalim
Snion can be considered to be from a completely mixed layer, though
actually the layer has a steep temperature gradient.
Inflows to the reservoir ucre considered to enter the epilimnion,
metalimnion, or the hypolimnion depending upon their temperature. In
flows were added to the respective layers at the end of the time
interval.
Since hypolimnion temperatures change very little during the
heating season, they were assumed to remain constant throughout the
period of analysis. Any underflow currents in the hypolimnion were
mixed with the entire hypolimnion water. Inflows to the metalimnion
were mixed with that volume present at the end of a time interval,
and a new temperature determined.
Over the heating period, provision was made for increasing the
epilimnion depth, which could be taken as constant or changing
depending upon experience 'ith the reservoir being considered. As the
epilimnion volume increased, the "makeup" water was taken from the
metalimnion and mixed into the epilimnion. The temperature of the
U)
F
J I
epilimnion was assumed to equal the equilibrium temperature. When the
net rate of heat transfer at the water surface is zero (3H = 0), the
surface water temperature is considered to be at its equilibrium
temperature. When the equilibrium temperature is greater than the
surface water temperature, the net rate of heat transfer is positive
and the surface water temperature rises. When the equilibrium tem
perature is less than the surface temperature, the net rate of heat
transfer is negative and the water temperature decreases. Since
reservoir surface temperature increases in the spring, one would expect
that the equilibrium temperature is greater than the water temperature.
During the summer, the surface water temperature is approximately
constant, indicating that the equilibrium temperature is approximately
equal to the surface temperature. In the fall, the surface water
temperature decreases; therefore, the equilibrium temperature must
also be decreasing. Since the surface temperature is always striving
toward the equilibrium temperature, the equilibrium temperature should
be a relatively good predictor of the surface water temperature. The
following development of the equilibrium temperature is according to
Edinger and Geyer.111
If in Equation (12) the surface water temperature, Ts, is equal
to the equilibrium temperature, TE, then the net rate of heat transfer,
AH, is zero and the equation becomes:
HR = o w(TE+460)4 + (a+bU)(eEea) + .26(a+bW)(TETa) ........ (19)
where TE and eE have been substituted for T and e respectively.
Subtracting Equation (20) from Equation (12) eliminates HR:
AH = {y, o[(Ts+460)4 (TE+460)4] + (a+bW)(eeE)
+ (a+bW) 0.26 (TsTE)} ............................... (20)
Edinger and Geyer use a linear relationship to approximate
saturation vapor pressure and water temperature For tendegree Fahren
heit temperature ranges. The constants for the temperature increments
are shown in Table 1.
TABLE 1
PROPORTIONALITY COEFFICIENTS AND INTERCEPTS FOR
THE TEMPEPEATURE, VAPOR PRESSURE APPROXIMATION 14
Temperature Range, OF B, mmHg/OF Intercept C (8)mmHg
4050 .2910 5.47
5060 .4050 11.22
6070 .5553 20.15
7080 .6667 27.80
8090 .9900 53.33
90100 1.289 89.28
The term E in Table 1 is the proportionality coefficient in the
equation below:
(es e ) = (T TE) ............................. (21)
The expressions to the fourth power in Equation (20) were
approximated by the binomial expansion, keeping the first and second
order terms. The approximation is:
OYw(T s +460 Y460 [1 + 4Ts + 6( s)2] ............... (22)
T 460 460
Substituting for a and Y, and including the approximations
Equation (20) becomes:
39
H = [15.71 (TTE) + .051(Ts2TE2) + (a+bW) (TsTE)
+ 0.26 (a+bW)(TsTE) ... ........ ..................... (23)
Combining like terms yields:
AH = [15.7 + (0.26+8)(a+h'] (TsTE) +.051(TT2) ..... (24)
From Equation (24), Edinger and Geyer define the exchange
coefficient, K in Btu per square feet per day per degree Fahrenheit,
as:
K = 15.7 + (.026 + B)(a+bU) ................................. (25)
Using the binomial approximation for the fourth power term in
Equation (19), and the expressions:
K 15.7
a+b K 15.7.. ....... ............. ... .... ....... (26)
0.26 + S
and,
eE = TE + C ( ) ......................................... (27)
Edinger and Geyer reduced Equation (19) to:
TE + .051TE2 = HR1801 + K15.7 [eaC(e) + 0.26Ta ] ......... (28)
K K K 0.26+8 0.26+B
By estimating an initial temperature range for the constants B
and C (8), Equation (28) can be solved for the equilibrium temperature.
Ta, ea, HR and K are determined from meteorological data. If the
equilibrium temperature falls within the assumed range used for and
C (2), the calculation is complete. Otherwise, the calculation must
be repeated until TE falls within the estimated range.
Computer programs have been coded for calculating the epilimnion
temperature using Equation (28) and for keeping the inventory of water
volume and temperature in the three reservoir layers. The program list
ings are in the Appendix with appropriate flow charts and Glossary of terms.
The next section will discuss the effects of reservoirs on dissolved
oxygen and the technique employed in this study to account for these
effects.
Dissolved Oxygen
Dissolved oxygen levels in natural streams, with little or no
pollution, are generally near saturation, which depends upon the temp
erature and chloride content of the stream. However, with the intro
duction of a reservoir on a stream, dissolved oxygen concentration of
32
releases may range from 0 to 100% of saturation. This alteration
in dissolved oxygen levels is the result of many physical, chemical,
and biological processes taking place within the reservoir. The more
important processes will be discussed in this section along with recent
attempts at improving reservoir water quality.
Sources of Dissolved Oxygen
One of the principal sources of oxygen in natural waters is
atmospheric oxygen that is absorbed into the water across the air
water interface. Oxygen absorbed at the water surface is diffused
downward by windinduced turbulence. The rate of oxygen absorption
is proportional to the difference between the saturation level and
the dissolved oxygen present:
dC
dt K2 (Cs C) ................................ (29)
where C is the dissolved oxygen in milligrams per liter (mg/l); K is
the reaeration coefficient in daysl; and C, is the saturation level
of dissolved oxygen in milligrams per liter. The saturation values for
pure water can be expressed as:33
pure water can be expressed as:
Cs = 14.652 0.41022T + 0.0079910T 2 0.000077774T 3....... (30)
where T is temperature in degrees Centigrade'.
The reaeration coefficient is dependent upon the turbulence and
temperature of the water. Expressions for the reaeration coefficient
and its relation with temperature will be given in Chapter III. For
lakes, the reaeration coefficient has been reported as ranging from
0.05 to 0.10 day1 (see Babbitt and Baumann3).
Since dissolved oxygen is diffused into a body of water by wind
action, thermal stratification in reservoirs prevents atmospheric
oxygen from being transferred to the metalimnion and hypolimnion.
The absence of stratification in the epilimnion virtually insures
saturationdissolved oxygen levels within the epilimnion.
Photosynthesis by phytoplankton and aquatic plants is another
source of oxygen in water. Since sunlight is the catalyst for photo
synthesis, oxygen produced by this process is usually confined to the
layer absorbing the most sunlight. In reservoirs, this layer is the
epilimnion. Depending upon the concentration of phytoplankton, these
organisms can be an oxygen credit or debit. In the production of
oxygen, phytoplankton also respire, thus consuming oxygen. Large
concentrations of phytoplankton, "blooms," have been observed to almost
deplete oxygen from the water during night hours. Hutchinson21 points
out a case in which the oxygen varied from 18.7 mg/1 in the afternoon
to 2.2 mg/l in the early morning.
An important source of oxygen to the lowerlayer of a stratified
reservoir is inflowing water. If the oxygen content of the inflow is
larger than that at the level it enters the reservoir, the total oxygen
level will increase in that layer.
Dissolved Oxygen Sinks
The primary sinks of dissolved oxygen are chemical oxidation of
inorganic compounds and biological respiration. The presence of
oxidizable inorganic compounds may be from surface inflows, under
ground inflows, sediments and industrial pollution. In the presence
of oxygen, compounds such as sulfide and nitrite are oxidized to
sulfate and nitrate. Reservoirs which have large concentrations of
these compounds may have very little oxygen in the lower strata.
Biological respiration in reservoir waters may be from animals,
phytoplankton, zooplankton, or bacteria. Since fish require oxygen
to survive, they will certainly exert a demand on the oxygen resources
of a reservoir; however, the relative amount of oxygen consumed by fish
is probably very small.
Dissolved organic compounds are biochemically oxidized by
bacteria in the reservoir. The bacteria consume oxygen as they degrade
the organic material. The rate at which oxygen is consumed is directly
proportional to the amount of organic material present. In some lakes,
this rate has been observed to range from 0.07 to 0.28 mg/l/day.
The zooplankton and phytoplankton are sometimes interdependent.
Zooplankton consume oxygen as they feed on phytoplankton and bacteria
in the water. Naturally, large algal growths will spawn large pop
ulations of zooplankton and hence create an even larger oxygen sink
during night hours.
Dissolved Oxygen Profiles
The sources and sinks of oxygen combine to cause a variation in
dissolved oxygen with depth and time. This effect is shown in Figure 8
DISSOLVED OXYGEN in mg/I
I '
10 
20
30 
40 JUNE 30,1967
\\ \ AUGUST 16,1967i
50 DECEMBER 7, 1967.c
.60
S70 
8 90_
I
100 
S110 4
120
140 / 0
FIGURE 80. DISSOLVED OXYGEN VARIATION with
TIME and DEPTH
TEMPERATURE in C
CLARK HILL RESERVOIR
SAVANNAH RIVER
CORPS OF ENGINEERS
JUNE 3017 6
AUGUST 16,107
DECEMBER 7,1 G
I
d
FIGURE 8b. TEMPERATURE VARIATION
with TIME and DEPTH
/
/
I
9
I
/
with corresponding temperature profiles for data collected on the
Clark Hill Reservoir in Georgia. The dissolved oxygen is seen to
stratify as the reservoir undergoes thermal stratification. The
oxygen concentration is constant in the epilimnion, decreases in the
metalimnion, and increases in the hypolimnion before reaching a
minimum at the bottom. The high concentration in the epilimnion is the
result of surface aeration and phytoplankton. As mentioned earlier,
the oxygen sources are confined to the epilimnion with little, if any,
oxygen being diffused to the lower layer.
The point of lowest dissolved oxygen occurs in the metalimnion,
which is not uncommon and is termed the metalimnetic minimum.21 The
most widely known theory of the minimum is that it represents a region
into which dead phytoplankton, organic material, and zooplankton have
settled from the epilimnion. It is believed that the increasing
viscosity of the water in the metalimnion prevents this suspended
material from settling into the hypolimnion. The result is a large
growth of bacteria and zooplankton consuming the oxygen as they feed
on the dead phytoplankton. Churchill and Nicholas 35 report data
collected on the Boone Reservoir in Tennessee which seem to corrob
orate this theory. Profiles of zooplankton counts in the reservoir
show a maximum of zooplankton at the metalimnetic minimum.
The dissolved oxygen in the hypolimnion water of Figure 8a
decreases from June 30, 1967,to August 16, 1967. Evidently there is
oxidizable organic material being degraded by bacteria. At a depth of
one hundred ten feet, the rate of oxygen depletion is 0.077 mg/l/day,
which is in agreement with Hutchinson.
By December 7, 1967, the dissolved oxygen concentration is
uniform in the reservoir. This is to be expected after the fall turn
over. If another profile were shown for March or April, it probably
would show very little thermal gradient and uniform dissolved oxygen
concentration.
Improving Release Water Ouality
Kittrell discussed several alternative methods for improving water
quality in released water. The most promising methods were multilevel
penstock intakes, turbine reaeration, and reservoir mixing.
Multilevel intakes permit the selection of that quality of water
which best meets downstream needs. The installation of multilevel
penstock intakes adds considerable expense to the reservoir capital
costs. However, there is very little operating expense and benefits
from increased flexibility in reservoir management to meet downstream
quality requirements may offset the capital cost. It should be pointed
out that no methodology exists at this time for estimating the benefits
accruing from multilevel intakes, but the methodology developed in
this study will provide a basis for comparison with single outlets..
Turbine reaeration is easily employed at power installations where
the turbines are higher than the tailwater level, creating a substantial
vacuum at low flows. Vents opened to the atmosphere introduce large
quantities of air in the discharge water. The loss of power when the
vent is opened may range from 2.5 to 10.6 % of the total poL;er
being generated. The loss is inversely proportional to the head on
the turbines. In a test case of several hundred feet of head, the
power loss was negligible.37 Oxygen increases may range from 5000
to 10,000 lbs/day/1000 cfs.38
Since the work of Symons et al.39 and Irwin et al.,40 much
attention has been given to reservoir mixing as a means for controlling
water quality. It was shown that pumping water from the deepest portion
of small lakes would effectively modify the thermal gradient and in
crease the dissolved oxygen at lower depths. In short, the quality of
the lakes was improved. The Corps of Engineers in the South Atlantic
Division has completed one year of work on the first attempt to de
stratify a reservoir.32 Air diffusers were used at one location in
Alatoona Reservoir in Georgia to disrupt thermal stratification and
improve dissolved oxygen concentrations. The Corps is now in the
process of evaluating the data collected in the summer of 1968. No one
has yet determined the benefits of reservoir mixing, but a reasonable
comparison can be made from the methodology developed in this study.
Reservoir mixing causes warming of hypolimnion water in the
reservoir. If the reservoir is completely mixed for the entire heating
season, then the cool hypolimnion water will not be available for down
stream fisheries or for cooling water in steamelectric power plants.
This could be a severe detriment to overall basin economy and should be
included in analysis of reservoir mixing.
In the preceding discussion, the effects of reservoirs on
dissolved oxygen have been outlined and some methods of improving
water quality have been mentioned. The next section will outline
the method employed in this study to predict dissolved oxygen in
reservoir releases.
Method Used in This Study
In several water quality studies (e.g., Bramhall and Mills,
Davis,5 Worley et al.,41 Feigner42) little, if any, consideration
is given to temperature and dissolved oxygen variation in reservoir
releases. Bramhall and Mills, as well as Davis, make no mention at
all of reservoir water quality in their evaluation of lowflow aug
mentation. Worley et al. assume the dissolved oxygen of reservoir
releases to be at the saturation level. Feigner does evaluate the low
flow augmentation requirements due to arbitrary changes in temperature,
but does not attempt to control temperature in the releases at the
reservoir.
Probably the reason few attempts have been made to account for
the dissolved oxygen variation is that little work has been done to
develop a predictive model for dissolved oxygen in reservoirs. Due
to the complex biological processes taking place in reservoirs and the
variation in the processes from time to time and from reservoir to
reservoir, it is very difficult to formulate a dissolved oxygen model
which can be applied to reservoirs. In view of the previous dis
cussion on the effects of dissolved oxygen, assuming dissolved oxygen
saturation in the release water would be unrealistic and would
result in a bias in allowable :aste loads. Saturation levels in
releases probably would be maintained in the stream until waste loads
are discharged. With low dissolved oxygen releases, the stream will
receive oxygen through reaeration; but, depending upon the degree of
turbulence in the stream and time of travel to the point of waste
discharge, the dissolved oxygen level may be much less than saturation.
49
Therefore, the waste loadings in the saturated release may be much
larger because of total oxygen resources at the point of waste discharge.
In an attempt to better represent the dissolved oxygen content
of released water, a method for permitting dissolved oxygen to vary in
time and among layers is employed. This method consists of using a
matrix of coefficients which vary both vertically and horizontally.
The variation of the coefficients in the vertical represents the change
in the percentage of saturation for each time period. The horizontal
variation represents the change among layers. Figure 9 illustrates
the formation of this matrix.
E M H
t=l a11 12 13
2 a21 22 23
3 a31 32 33
o
a a a
n anl n2 n3
FIGURE 9. MATRIX FOR REPRESENTING DISSOLVED OXYGEN
VARIATION WITHIN THE RESERVOIR
An example of the use of this matrix might be that in the first
time interval the reservoir is completely mixed with dissolved oxygen
at saturation. For this interval al = 1.00, a2 1.00, and a 3
1.00. In the second time interval the dissolved oxygen might vary as
follows: 95% saturation in the epilimnion, 80,% saturation in the
metalimnion, 85' saturation in the hypollmnion. Then the second row
of coefficients is a21 = .95, a22 = .80, and a23 = .85. The matrix
is filled in this manner from the first time interval until the nth
interval, which in this study represents the last time interval.
The use of this matrix: permits the calculation of dissolved
oxygen within the reservoir given the temperature of the water in
each layer. A method of this type can be adapted very easily to data
taken on existing reservoirs and can be employed in the analysis of
proposed reservoirs to see what levels of dissolved oxygen in the
reservoir cause downstream water quality problems. As simple as it
may be, this method is an improvement over previous work in which no
variation in dissolved oxygen is assumed.
In this chapter, the effects of impoundments on temperature and
dissolved oxygen have been discussed. The marked variation of temper
ature and dissolved oxygen in time and space were illustrated by
Figures 3 and 8. Discussion was presented on the causes of these
variations, models used to describe them, and the methods employed
in this study to predict temperature and dissolved oxygen changes.
In addition, alternatives for improving water quality in reservoir
releases were mentioned, and it was stated that the methodology
developed in this study vill provide a means for comparing these
alternatives.
Upon release from the reservoir, temperature and dissolved oxygen
do not remain constant, but vary as the stream flows through the basin.
In Chapter 111, the models used in this study to describe the variation
in temperature and dissolved oxygen along the stream profile will be
presented.
III. STREAM TEMPERATURE AND DISSOLVED OXYGEN MODELS
In Chapter II, the models used in this study to predict the
effects of reservoirs on temperature and dissolved oxygen were
presented. Since this study deals with a reservoirriver system, it
is necessary to develop mathematical models to predict the downstream
variation of the water released from the reservoir. The changes in
stream quality will be discussed in this chapter, and mathematical
models used to describe these changes will be derived.
Stream Temperature
Since the water in the deeper layer of the reservoir is usually
at a low temperature, upon release it will approach the equilibrium
temperature determined by meteorological conditions. The rate at which
the water approaches equilibrium will depend upon the difference
between the water temperature and the equilibrium temperature, the
wind over the river, and the area of the stream water surface.
In the previous chapter, the concept of equilibrium temperature
was discussed and it was stated that water temperature is continually
approaching equilibrium. The difference between the equilibrium
temperature and the water temperature is the driving force for the
rate of change of water temperature. .When the difference between the
equilibrium and water temperature increases, the rate at which the
water temperature changes '.ill also increase.
51
The rate of temperature change also depends upon the wind. If
the water temperature is increasing, meteorological conditions are
such that the equilibrium temperature of the water is greater than
the in situ water temperature, causing the water to absorb heat.
Since the heat is provided from the air over the water, the more air
brought into contact with the water, the faster the water temperature
will increase. The source of air movement is wind; therefore, the
greater the wind velocity, the faster the rate of change of water
temperature. A large percentage of the heat exchange between the
stream and the overlying atmosphere takes place at the water surface.
If the water surface is small in comparison to the total volume of the
stream, the stream will exchange heat with the air at a slower rate
than if the surface were larger.
Heat exchange in streams is affected by thermal stratification
in the stream. Thermal stratification is caused by a rapid warming
of the surface and by high temperature discharges into the stream.
The thermal gradient may remain in sluggish streams, but it is
usually destroyed by turbulent mixing in swiftly flowing streams.
To make calculations easier, complete mixing usually is assumed.
Model Used in this Study
The derivation of the stream temperature model will be after
Edinger and Geyer.14 It is convenient to begin with Equation (14):
BT 3 3T a 3T 3 aT
Tt (E + a) 3 ]+ [97 (E + a) y + [ (Ez + a) ]
T + v, 3T + vT + S (31)
ax a' z z ] Cpc
p
where T is temperature in degrees Fahrenheit; t is time in days; Ex,
Ey,and E. are coefficients of eddy diffusivity in the x, y,and z
directions in square feet per day; a is the coefficient of molecular
diffusivity in square feet per day; v., vy,and vl are velocity components
in the x, y,and z directions in feet per day; and S is the external heat
source in Btu per cubic feet per day.
The xdirection is along the length of the stream, the zdirection
is vertical, and the ydirection is laterally across the scream. In
most stream temperature models, the first assumption is that the stream
can be considered to be completely mixed (e.g. Raphael,18 Duttweiler, 3
and Velz and Cannon ).
If a stream is completely mixed, there will be no lateral or
vertical gradients of temperature; i.e., i. OT = 0 Therefore,
two diffusion terms and two advection terms become negligible. This
leaves:
3T 3 T (T S
L=[S x (E + a) x] Vx x +p ...... (32)
For the purposes of this study, two additional conditions will be
required in the stream. First, the temperature at any point along the
stream will remain constant during the time interval for which the
temperature is predicted. This is the condition of steadystate,
T = A condition of this nature will negate diurnal temperature
at
variations. Secondly, the primary mode of heat transfer will be
advection. Under this condition, the heat transferred by diffusion
will be considered negligible when compared with advected heat. Thus,
the heat source term is balanced by the heat advection term:
x 3T S ..................................... (33)
ax pc
Taking the intergral of Equation (33) over the depth, z, yields:
d d
vx T dz = S dz.................... (34)
ax pc
0o p
For the purposes of this study, all heat exchanged into or out of the
body of water is assumed to go through the surface; hence, the term
ddI
f S dz is the net rate of heat transfer, AH. The term vx dz can
o
be thought of as the stream flow per unit width of stream. Then,
Equation (34) becomes:
Sr = r;H (35)
PCp dm vx 3 H ......................................... (35)
ax
where d will be taken as the mean depth of the stream. If in Equation
m
(24) of Chapter II the second order terms are neglected, the net rate
of heat exchanged becomes equal to the product of the exchange coeffi
cient and the difference between the water temperature and the equilib
rium temperature:
AH = K (Ts TE) .......................................... (36)
where K is the exchange coefficient defined in Equation (25) in Btu per
square feet per degree Fahrenheit per day:
K = 15.7 + (0.26 + 8) (a+bW) .............................. ... (37)
Substituting Equation (36) into Equation (35) and solving for the
rate of change in temperature in the downstream direction yields:
aT K (T T)
ax c38)
3x p d m v ....................................... (38)
pCp dm vx
Comparing Equation (38) with previous discussion on the change in
temperature reveals that the change in the longitudinal gradient is
a function of the driving force, (Ts TE), the exchange coefficient,
K, the mean depth and stream velocity. Further scrutiny shows that the
wind velocity is present in the exchange coefficient,and the surface
area of the stream is found in dm v by substituting in the continuity
equation for velocity:
dm v = d = d x ... ...... ............ (39)
A dm Ex As
where x is distance dounstream'in feet; As is stream surface area in
square feet; A is average crosssectional area in square feet; and Q is
the stream flow in cubic feet per second; and B is the average width
in feet. Thus it can be seen that all the factors which were said to
influence the rate of change of the release water are present in
Equation (38). Integrating Equation (38) between the limits x = 0 and
x = Xi yields the predictive equation used in this study.
Tf = (Ti TE) er2 + TE ............ ................ (40)
where Tf is the temperature at xl in degrees Fahrenheit; Ti is the
temperature at x = 0 in degrees Fahrenheit; and
K xl
r2 pcpd vx ...................................... (41)
If the river downstream from the reservoir is divided into several
reaches, Equation (40) can be used successively in each reach to
calculate Tf at the end of the reach,with Ti being the temperature at
the end of the previous reach. The reaches should be divided in such
a manner that the flow and hydraulic properties of the stream can be
considered as constants within the reach. To use Equation (40) one
must have available meteorological data to determine the equilibrium
temperature and the exchange coefficient. The temperature predicted
by Equation (40) is the average temperature during the time interval
for which the stream flow and meteorological data are considered
constant. The shorter the time intervals between data collection, the
closer the estimate becomes to observed stream temperature behavior.
However, temperature predicted from monthly data should be approximately
equal to the mean monthly temperature. It should be noted that pre
dictions for reaches whose travel times are longer than the intervals
between data collection can not be made because steadystate conditions
would not exist. The data averaging interval, the period for which the
data can be assumed approximately constant, can not be less than the
time for which the equation is assumed valid or the steadystate
condition is violated.
Dissolved Oxveen
The sources and sinks of dissolved oxygen in reservoirs were
discussed in Chapter II. It was pointed out that almost no methodology
exists for predicting dissolved oxygen in reservoirs. However, in the
case of streams, several attempts have been made to develop predictive
equations for dissolved oxygen. The equations estimate the dissolved
oxygen as a function of various parameters, which are determined from
data collected at the river or are assumed when data are unavailable.
The parameters generally are taken as constant throughout the reach of
the river to which the equation is applied. The.estimated values often
disagree with observed values with the result that little reliability
can be placed on the estimated oxygen sag curve. One of the major
reasons for this discrepancy is that biological regimes are assumed to
be constant, when actually they vary among rivers and even within the
same river. Hence, there is an increasing need for a predictive model
which reflects the uncertainties of the dissolved oxygen balance in
streams.
Streeter and Phelps 5 derived the most frequently used model for
predicting the stream dissolved oxygen profile. In their model, the
only oxygen sink :ias the utilization of oxygen by bacteria during the
decomposition of dissolved and suspended organic material. The rate
at which the organic material was oxidized was assumed proportional to
the concentration of the organic material. The equation used to
represent this process was the first order kinetic equation:
dL
d t 1 K L ........................................ (42)
where L is the biochemical oxygen demand (EOD) in milligrams per
liter at time t and K'1 is the deoxygenation coefficient in days1
Integration of Equation (42) yields:
L = La Kt (t3)
L = La e K It ......................................... (43)
here La is the ultimate BOD in milligrams per liter.
The only source of oxygen in the StreetcrPhelps equation was
surface reacration. The equation for the rate of dissolved oxygen
transfer was also expressed in a first order equation as:
d r
d = 2 D (
where D is the dissolved oxygen deficit in milligrams per liter and K2
is the reaeration coefficient in days1. The dissolved oxygen deficit
is expressed as the difference between the dissolved oxygen saturation
value, Cs in milligrams per liter, and the actual dissolved ox..ygen, C
in milligrams per liter:
D = Cs C .................... ......... ...... (45)
The reaeration coefficient was determined to be proportional to
the stream velocity and inversely proportional to the stream depth.
Their work on the reaeration coefficient has been confirmed by several
authors with only minor changesin exponents and constants. A recent
formula for the reaeration rate was presented by Langbein and Durum:46
V
K9 = 3.3 1.33 ..................................... (46)
dm
where d is the average depth in feet and v is the stream velocity in
feet per second.
Streeter and Phelps combined the oxygen source and sink into the
equation for the rate of change of dissolved oxygen deficit:
dD
dt K1 L K2 D .................................. (47)
Integrating Equation (48) yields the classical StreeterPhelps
Equation:
K1 La K1t K2t K2t
D K La (e t e ) + Da eK2 ............ (48)
K2KI
where Da is the initial dissolved oxygen deficit in milligrams per
liter. The dissolved oxygen sag curve, aS predicted by Equation (48),
is shown in Figure 10 The curve shown is typical of the oxygen
profile below waste discharges. The accuracy with which the Streeter
Phelps Equation predicts the sag curve depends upon many factors; among
them are: the selected. deoxygenation coefficient, KI; the variation
of the reoxygenation coefficient, K9; other sources of deoxygenation
such as phytoplankton respiration and benthol demand; and oxygen
production from photosynthesis.
Variation of K1 and K9
Streeter and Phelps determined that the deo::ygenation coefficient
was 0.23 days1 in their study. However, the coefficient varies from
/L
0 U
C:)
0
I IaJ
S I
0(fl
O
1/5W u! N39),XO 03A10SSIG
river to river and for different wastes, and estimates have ranged
l
from 0.1 to 0.7 days Because of the variation in K1 and K2,
Churchill and Buckingham47 used linear regression to develop a pre
dictive equation for dissolved oxygen. They felt that the results
of using K1 and K2 as constants for predicting dissolved oxygen under
varying conditions were highly questionable. Kothandaraman48 studied
K1 values calculated from data on the Ohio River and found them to be
indistinguishable from a random sample. He also found the K1 values
to vary randomly according to Gaussian probability with a mean of
0.173 day1 and a standard deviation of 0.066 day1
Variation in the reaeration constant has been found in several
investigations. 33, 49, 50 In each study, empirical expressions were
derived to predict K2 as a function of stream velocity and depth.
Langbein and Durum46 developed Equation (46) to predict K, over a
wide range of data collected in the field and laboratory. They
presented a graphical representation of the data points which showed
considerable scatter of the points about Equation (46). Kothandaraman8
examined data reported by the TVA33 and performed a regression analysis
to obtain an empirical equation for K2 of the form:
K2 = 5.827 v'924
2 = 5.327 .* 924 .............................. (49)
dm1.705
He found that percentage variation about Equation (49) was normally
distributed with a mean of zero and a standard deviation of 0.368.
In view of Kothandaraman's findings about the random behavior of
K1 and K2 it is easier to understand the discrepancy between predicted
and observed deficits. If K1 and K2 can be represented by a random
function, the deficit in Equation (49) is a function of random
variables and has a probability distribution of its own. In this
study, the dissolved oxygen deficit will be assumed to be a function
of random components and thus will have a mean and standard deviation.
This will permit a probabilistic expression for the effect of the waste
loadfor example, 100,000 lb BOD/day may cause a violation 5% of the
time. This means that if the same stream flow and waste loads were
combined on one hundred occasions, one would expect a violation of the
stream standard on five occasions.
Additional Sources and Sinks of Dissolved Oxygen
Other investigators have attempted to add additional sources and
sinks of dissolved oxygen to the StreeterPhelps Equation. Modifications
are warranted in situations where these additional sources and sinks play
significant roles in the stream's dissolved oxygen balance. Thomas5
presented a technique for stream analysis in which he added a rate
constant K3 to Equation (42) to represent the removal of organic material
by sedimentation. This was the first modification to the StreeterPhelps
Equation. Other investigators have added parameters such as Kb to account
for organic addition from runoff and scour, Db for benthol demand, R for
oxygen consumption by respiration of aquatic biota, and P for oxygen
produced by photosynthesis. When these parameters are included in
Equations (43) and (48), the result is a cumbersome equation for
dissolved oxygen deficit and organic removal:
d = (K1 + K L + K .......................... (50)
dt = iL K2D P + Db + ....................... (51)
An extensive data survey is necessary to estimate the values of
the parameters. It is often hard to obtain data from which reliable
estimates of the parameter can be made. For example, Camp52 found it
necessary to assume K3 = 0.1 in order to estimate K1 from the collected
data in a stream study. In his analysis, Camp found Db co range from
0.5 to 1.0 mg/l/day; (PR), 0.55 to 2.9 mg/l/day; and K3, 0.07 to
0.11 day1.
In many instances, one cannot neglect the terms represented in
the modification to the StreeterPhelps, but, on the other hand, to
try to determine the relative contribution from each term requires
extensive data collection. Also, the variance in K1 and K2 must be
taken into account when determining the other parameters. Because
of the difficulty in estimating the parameters, Moreau and Pyatt53
suggested lumping the added parameters to the StreeterPhelps
Equation into a composite term. These terms could be determined
from field data by subtracting the contributions of Kl and K2 in the
oxygen balance of the stream. Moreau and Pyatt represent the term in
the organic load by r in Equation (42) and in the deficit equation by
s:
dL
dt = K L + r ................................. (52)
dD
dt = KILK2 D s ................................ (53)
The signs of r and s terms are important. In Equation (52), r is
positive, which represents added organic loading, and in Equation (53),
s is positive, representing a dissolved oxygen sink. The combination
of these positive terms, therefore, causes the predicted dissolved
oxygen sag curve to be a conservative estimate. Equation (52) and (53)
will be used in the derivation of the dissolved oxygen equation used
in this study.
Dissolved Oxygen Equation Used in This Study
Since the stream temperature changes after discharge from the
reservoir, this factor must be included in the determination of the
dissolved oxygen saturation level. Liebman5t' proposes a method for
handling the changing saturation level by assuming a linear relation
ship within reaches:
Cs(t) = Csi + Bt ..................................... (54)
where
B = Csf Csi ..................................... (55)
Tt
Cs(t) is the saturation level at time t, Csi is the saturation level
at the top of the reach, Csf is the saturation level at the bottom of
the reach, Tt is the travel time in the reach, and B is the slope of
the saturation curve within the reach. This is a good approximation
since the saturation level is almost linear with temperature. The
slope B is positive uhen the water tenperatire drops within a reach
and thus represents a dissolved oxygen sink. The final representation
in differential form is:
dL
dt = K1 L + r ..................................... (56)
and
dD
dt K1 I K2 D B + s ............................ (57)
Integration of Equations (56) and (57) for steady, uniform flow
fields:
L = (La r) e t + _r ............................... (58)
K1 K1
and,
D K1La r t _+(Da BE) Kt r+s ) t eK t
D e e )(Da e  (I
2 '1 KI ":2
........... (59)
As pointed out by Kothandaraman,48 the "constants" K1 and K2 can
be represented by probability distribution functions each having its
own mean and variance. Due to the uncertainty of r and s, Moreau
and Pyatt53 suggest that each be assumed to be random variables which
may be represented by probability distribution functions. Therefore,
in this study K, K,, r and s are assumed to be random variables
which may be represented by a probability distribution function having
m 2 2 2
means P' rs and variances K 2 2 r, 2 One would
1 2 1 2
expect that K would be independent of K since K represents a
1 2 1
biological process and K2 is a physical process. For the same reason,
K1 is independent of r. However, K, and s represent biological
processes which, for purposes of this analysis, will be assumed to
be independent. It is felt that any dependence would contribute
only a small degree of variance to the dissolved oxygen deficit. In
a like manner, the remaining parameters are assumed independent.
Therefore, the dissolved oxygen deficit becomes a function of indepen
dent random variables and has a probability distribution of its own.
Meyer55 shows that a function Y of n independent random variables
l, x2 ..., x (for instance, Y = f (xl x ..., x )) has the follow
( o n I x n
ing approximations for its expected value and variance:
E(Y) = f ( ul' u2' ."' L ) + 1 2 ... (60)
2 i=l x2 i
i
Var(Y) 3 2 o 2 (61)
Var(Y) = 2 i 2 i2 ............ .......... (61)
i=l x
where i. and oi. are the mean and variance of the random variable x..
1 1 x1
whee wi ad a1 2arethemea ad varac fterndmvral
Taking the indicated partial derivatives in Equations (60) and (61),
the mean and variance of L and D are:
E(L) = L +1 {2L 5 2L
E2 1= 2 + 2 + 2L K 2) .........................
2 3Ka12 K12 + aK2 oK2
Var(L) = CL )2 o 2 + (aL2 a 2
3K K r
1 1 1 r
E(D) = D + 1 { 2D o 2 + ;2D a 2 + a2D 2 2D 2 }
2 K 2 K 1 K 2 K aYr r 3s2 s
Var(D) = (= D )2 a 2 + D
aK K ;K
1 1 2
o 2+
K2
(aD)2 o 2 + (;D)2 a 2
r r s s
(62)
(63)
(64)
(65)
The partial derivatives of Equations (62), (63), (64) and (65)
are given in Appendix A along with the expressions for E(L), Var(L),
E(D) and Var(D). The results of Kothandaraman48 indicate that the
distribution function of the dissolved oxygen sag is approximately
Gaussian. If such an assumption can be made, confidence limits of a
predetermined probability can be established for the deficit when the
expected value and variance are known. For instance, 90% of all
values of a Gaussian distribution lie within the interval p + 1.645 a,
where i is the expected value and a is the standard deviation
(2 = variance). In the case of dissolved oxygen, if the standard
required that dissolved oxygen should remain above 4.0 mg/l 95% of
the time, the standard would be met if the maximum deficit, calculated
as E(D) + 1.645 v'Var(D), did not cause the dissolved oxygen to fall
below 4.0 mg/l. When the minimum dissolved oxygen is larger than
4.0 mg/l, a maximum allowable waste load can be determined which
satisfies the standard.
Sensitivity Analysis
The expressions for the variance of L and D afford a convenient
method for determining relative sensitivities to the parameters K1,
K2, r,and s. The terms on the right hand side of Equation (65) are the
contribution of the respective parameters to the total variance. By
substituting nominal values of the parameters into Equation (65), one
can evaluate the partial derivatives and determine which parameters
cause the largest variance in D. Such an analysis w'as made for several
values of I:1, 12, r, and s. Kothandaraman found the deoxygenation
coefficient for the Ohio River to have a coefficient of variation, C ,
of 0,38 (Cv = i/o). If this value is assumed constant, for a K1 of
0.23 days', the standard deviation would be 0.0874 days1. This
relationship v\as assumed to hold for K1 in order to determine its
standard deviation. For the reaeration coefficient, the standard
deviation \'as assumed to be 0.368': of the mean, as found by Kothandaram.
The remaining nominal values are given in Table 2. A time of one day
was chosen to make the analysis since most river basins are subdivided
into reaches with travel times of about one day.
Examination of the results in Table 2 reveals that the variance
of D is caused entirely by K,' K2, and s however, the expected value
of D is dependent on all four parameters. The percentage of the variance
from K2 increases with increasing K2 indicating that in swift, shallow
streams variation in the deficit is due mostly to variation in K2
As :2 decreases the relative importance of KI in the variation of D
increases. Therefore, for deep, sluggish streams a reliable estimate
N rN r r C) tr o cY < mo
fl CM Ci N Ci .i C4
cc cn m m r mr) i1 i 
C r r r~ 00 r, 0h ON m 0%.
0 U
U <
u<
] p u L
0
U U
cc* ;S
va 4
r r r l
co Co CO 0o
c' CO co r.
iD j 3
a:" \Q % < co
I o LO L CO)
r'l
1\0 C]
I' ,"
'. C0
 I 0 r Co c') co r
N N c 1 r N H 1
iL N 0 O'. . r1 0 C)
C'" fC' c' j 1 c" ",r C
C 01 .C CO 0co 1 . e
 e n CI'n n L) L) .1 Ln)
0 o ) 0 0 0 L
m rIn 0rn o) Cp r~ Ln
0 rC' oCl co T' C'I 4 rN
t C' 0 o;. 1 .1 t 0 C
C') C') 44 cN j r . C^
:4
0 0 0 0 0 0 0 LC r
Co ". C) C CD C) L' )
C') C') C') C') ell C') e'l
Cen C C c'i (" C' fr"
rNl r C r] T r NC
N r
(, iO,
F~ r~
II
0. . rN
r'. C., CO
0 %D0 0
Lr) 0". 'C'
' N" 0.
Lr)
L)n L'i N
ILn %0 C
CO LO N
0
.4 C')
Cf m mC C) CI C')
N 'o. r N r^ el
r4
L1 Li)
C'I C`]
II II
. ui
O O
LC) LC)
0 0
II II
0 c
>,
,4
a E
OI
O
o o
II II
0
(N
r? n)
Q
CMI
:4
V)
4
<
Flj
pa H
< i
7
H
V,'
of K1 would be more important than streams with high :K2 values. An
increase in KI from 0.23 to 0.40 caused a large increase in the
variance of D. Such increases in K1 might occur during %,arm summer
months, the period of most pollution problems.
Large changes in r and s cause little change in their contributions
to the variation of P. Vlhen s is negative, corresponding to a photo
synthetic oxygen source, the expected value of D decreases and the
relative contribution of K2 in the variance decreases. This indicates
a decreasing importance in K2 for streams with a photosynthetic oxygen
source.
Also listed in Table 2 is the effect of K1 and r on L, One can
see that 1K has considerably more contribution in the variance of L
than r. Therefore, the estimation or determination of Ki is more
important than r. In summer the sensitivity analysis reveals that
most of the variation in D is caused by variation in KE and K2) while
in some cases s plays a more important role. It would seem, therefore,
that in field studies more time and data collection should be devoted
to determination of i1 and Y2 than to r and s.
Application of Eqilations (63) and (65)
The application of Equations (63) and (B5) requires special
consideration if the equations are used in piecemeal fashion from reach
to reach. In such cases the dissolved oxygen deficit and remaining
BOD at the end of a reach are used as the initial deficit. Da, and
ultimate DOD. La, in the next reach. If Da and La are deterministic
values at the beginning of the first reach and take on the values of
P and L at the beginning of the second reach, then Da and La are no
longer deterministic, but are stochastic. Since the deficit then becomes
a function of two more stochastic variables, the partial derivatives of
D with respect to these variables must be included in the determination
of the variance of D. Therefore, ( L)2 2 should be added to Equa
@La La
tion (63) and (. D)2, 2 + (2 D)2, 2 should be added to Equation (65).
aDg Da aLg a
The values for a 2 and o 2 are simply the variance of L and D at
La Da
the end of the previous reach. The expected values of L and D at the
end of a reach are used for La and Da in the next reach.
In this section, a stochastic model for dissolved oxygen has been
formulated as a function of various parameters which, because of the
uncertainty in assuming them constant, are taken to be independent
normally distributed random variables. The use of the stream temper
ature and dissolved oxygen models for determining the allowable waste
loading will be discussed in the following section.
Determination of Allowable Waste Loading
The thrust of many pollution studies is toward finding the maximum
allowable waste load the stream can assimilate under various conditions
ofstream flow and stream temperature. To determine waste loadings there
are several parameters in the stream and dissolved oxygen models which
must be known. In this study, stream velocity and depth are parameters
in the stream temperature model and are needed for determining the
reaeration rate in the dissolved oxygen equation; the equilibrium
temperature for the stream is fundamental to the stream temperature
model; the temperature at the upper and lower ends of each reach
determine the dissolved oxygen saturation slope and the values of the
deoxygenation and reaeration constants within the reach; and the maximum
deficit determines the allowable waste discharge.
Although the stream's velocity and depth are important in the
overall methodology, the model used to predict the relationship is
of little significance as long as it predicts reasonable values for
the particular stream being studied. In this study, velocity and depth
will be in terms of stream flow:
V = a 0Bv .................................... (66)
dm = h QBh . ................................... (67)
where a B,, ah, and Bh are constants in each reach, which are
determined from observed data. The first step in determining the allow
able waste load under given conditions is to calculate the mean velocity
and mean depth for each reach. The flow in any reach is assumed steady
and uniform. Tributary inflows to a reach occur at the upper end of
the reach.
The second step is the evaluation of the meteorological data to
determine the equilibrium temperature. (The procedure for calculating
the equilibrium temperature was discussed in Chapter II.) The
equilibrium temperature is assumed constant for all the downstream
reaches; however, when geographical situations indicate this assumption
to be false, another equilibrium temperature can be calculated for the
remaining stream sections.
After determining the equilibrium temperature, one calculates the
stream temperature profile by the stream temperature model. The
dcoxygenation and reaeration coefficients are calculated using the
average stream temperature within the reach.
The final step is to determine the minimum dissolved oxygen for
a given waste discharge. The sag curve is calculated at the end of
each reach, with the remaining organic load and dissolved oxygen deficit
becoming the initial conditions in the next reach. The minimum dissolved
oxygen is found when the slope of the sag curve becomes positive. Tracing
back to the beginning of the previous reach and calculating the dissolved
oxygen at the third points of each of the last two reaches will approxi
mately locate the minimum point if the change in slope is not too large.
The minimum dissolved oxygen is tested against the allowable minimum and
adjustments are made in the loading until the minimum is within 0.02 mg/l
of the standard. The waste loading which satisfies the dissolved oxygen
standard within 0.02 mg/l is the maximum allowable waste load for the
given stream flow and reservoir release temperature. The program listing
for the determination of allowable waste loading is included in
Appendix B.
In the next chapter a discussion of the factors which influence
a stream's waste assimilative capacity will be presented. The concept
of a tradeoff between release flows and temperatures will be introduced,
with hypothetical curves depicting the relationship between waste
assimilative capacity and the flow and temperature parameters.
IV. THE TRADEOFF BETWEEN STREAM TEMPERATURE AND STREAM FLOW
As mentioned in Chapter II, few investigators consider
discharge temperature variation from reservoirs in water quality
models. However, the variation in temperature may be quite significant
and provides the basis for considering a tradeoff between the temper
ature of the release and the quantity of discharge. The term "trade
off" in this study will mean that water at one temperature and
discharge rate may be exchanged for water at a different temperature
and discharge rate while still maintaining specified water quality
standards downstream. The exchange is possible because of the
variation in temperature within the reservoir.
The Effects of Temperature on Waste Assimilative Capacity
The effect of temperature on dissolved oxygen in streams is
twofold. First, temperature has a pronounced effect on the
biological activity in the stream, and secondly, the physical parameters
of the stream vary with stream temperature. In general, biological
activity increases with increasing temperature. This phenomenon is
caused by the dependence of enzyme reaction rates on temperature.
Increasing temperatures cause higher reaction rates, which in turn
stimulate metabolic activity and cause the oxygen requirements of the
organism to increase. Thus, increasing the temperature in a stream
causes faster oxidation of soluble organic material and a sharper race
of decrease in the oxygen sag curve. The variation of the deoxygenation
coefficient with temperature is often expressed as: 57
K T = K 046(T 20..) (68)
1 T 2 0 ..
where K, 20 is the value of the coefficient at 20C; T is the temper
ature in degrees centigrade; and K1, T is the value of the coefficient
at temperature T.
The physical properties affected by temperature relate primarily
to dissolved oxygen resources in the stream. Increasing temperature
causes the solubility of oxygen in water to decrease. This lowers the
saturation level of dissolved oxygen (Equation (30)). Also, increas
ing temperature decreases the viscosity of water, the rate of molecular
vibration, and the rate of oxygen dispersion. These factors cause the
reaeration coefficient to increase with increasing temperature. The
reaeration rate as a function of temperature is:33
.016(T20)
K K e 16 T 20) .......................... (69)
_, T 2, 20
where K2 20 is the reaeration coefficient at 200C; T is the temper
ature in degrees centigrade; and K, Tis the reaeration coefficient
T
at temperature T.
The high summer temperatures in streams cause lower levels of
dissolved oxygen saturation and higher coefficients of reaeration and
deoxygenation. However, comparing Equations (6S) and (69) reveals
that changing temperatures cause larger changes in K1 than in K2. A
change of 5"C from 200C causes a 20% change in K1 and a 9% change in
K2.
When the biological and physical effects are combined at high
temperatures, one observes that tlhe demand on the stream's dissolved
oxygen resources is highest at the time when the total available oxygen
is lowest. The net result is a decrease in the stream's waste assim
ilative capacity during the summer months. If cool water were avail
able in a reservoir and if it were of suitable quality, then releasing
the cool water would increase the assimilative capacity. However, the
degree of increase in assimilative capacity will depend upon the
distance of the waste outfall from the reservoir. Since the cool
water will warm as it passes through the basin, waste discharges far
do\.nstream will not have the benefit of the cool water.
Examples of the effects of temperature on the waste assimilative
capacity were reported by lWorley and Krenkel 5:' et al. Worley was
investigating low flo.. augmentation requirements in the Uillamette
River Basin. He found that increasing the temperature from 20'C to
25"C caused a 14" decrease in the total allowable waste load in the
basin for a streamflow., of 5500 cfs. Krenkel et al. studied the effects
of impoundment on the Coosa Piver in Georgia. Under free flo, condi
tions of 940 cfs, they reported large variations of allowable waste
loading with temperature. Decreasing the temperature from 20CC to 15cC
increased the allowable loading 40": from 45,000 Ibs BOD/day to 63,000
Ibs EOD/day, while a 5"C increase to 25"C decreased the allowable
loading 20'. to 32,000 lbs EOD/day. Such evidence supports the tradeoff
among temperatures for the purpose of improving waste assimilative
capacity.
The Effects of StreamFlow on Waste Assimilative Capacity
Streamflow influences the extent of dilution of the organic
waste and the reaeration coefficient of the stream. By increasing
flows, the concentration of the waste is decreased, and one would
expect a corresponding increase in the assimilative capacity of the
stream. Discussion will be presented in this section which questions
the validity of this assumption.
The equation for determining the reaeration coefficient was
given in Equation (46) as:
K = 3.3 v ........................... ......... (70)
2 d 1.3 3
m
Increasing flow increases the velocity of the stream as well as its
depth, but the reaeration coefficient is influenced to a greater extent
by the depth of the stream than by the velocity. Langbein and Durum46
point out that the reaeration coefficient at any location decreases
at about the 0.13 power of the discharge. Therefore, increasing flows
causes a decrease in the capacity of the stream to reaerate.
Lowflow augmentation as a means of water quality control is
usually employed in the summer months. If earlier releases at the
reservoir have emptied the colder water, the releases will be at
warmer temperatures, and by the time the flow reaches the waste
discharge point the temperature will probably have reached equilibrium.
The warmer temperature will cause high K1 values in the stream. The
sensitivity analysis in Chapter III indicated considerable variation
in the dissolved oxygen deficit for high K1 values and marked increases
in the expected value of the deficit. The sensitivity analysis also
76
revealed larger variations in the deficit with decreasing K2 values,
as well as an increase in the expected value of the deficit. In view
of these results, augmenting low flows with warm water releases would
cause an increase in K1 and a decrease in K2 with correspondingly larger
variations in the dissolved oxygen deficit. Therefore, considerable
case study should be made before lowflow augmentation is used as a
method for water quality control. Arbitrarily releasing water without
regard to temperature and influences on K2 may create water quality
problems instead of solving them.
The Stream Flow Stream Temperature Tradeoff
In the previous sections, the effects of temperature and flow
on waste assimilative capacity were discussed. With this background,
it is possible to infer the existence of a tradeoff between water
temperature and stream flow.
The effects of increasing temperature were said to decrease the
stream's assimilative capacity. With this in mind, hypothetical
curves for assimilative capacity versus temperature are shown in
Figure l1a. The curves are shown for various flow rates. For flow
rate Q1, the assimilative capacity decreases with increasing release
temperature. For the same temperature, increasing the flow increases
the assimilative capacity, but the effect of this increase is shown
to decrease as the curves slope to the right. The tradeoff can be
seen by noticing that the assimilative capacity for the coldest discharge
at flow rate Q1 exceeds the assimilative capacity at flow rate Q3. If
the waste discharge were at a distance from the reservoir such that
temperature had no effect on KI, then che curves '..ould become horizontal
lines as shown in Figure lla.
RELEASE TEM
FIGURE a. WASTE ASSIM.ILATIVE CAPACITY vo. TEMPERATURE
FLOW
FIGURE b. WASTE ASSIMrILATIVE CAPACITY vs. FLOW
FIGURE II. HYPOTHETICAL CURVES
for WASTE ASSIMILATIVE CAPACITY
vs. RELEASE TEMPERATURE and FLOW
The tradeoff is illustrated, also, in Figure lib. For the
warmest temperature, T3, the flow requirements become increasingly
large for small increases in waste assimilative capacity. This is
based on the assumption of diminishing returns for increasing flows.
If no tradeoff exists, the curves of Figure llb wouldd collapse to
the bottom curve.
For situations in vhich tie tradeoff in temperature does not
exist, the afterr quality variations in the reservoir have very little
effect on the stream's assimilative capacity. The only alternatives
to meeting the stream standard become waste treatment, with methods
such as variable waste treatment, instream aeration and increasing
flow's. Since the methodology developed in this dissertation is based
on the consideration of reservoir water quality variations, the
existence of a tradeoff is required for application of the methodology.
In Chapter V, earlier vater quality models %,ill be discussed and the
strategy employed in this study v'ill be d.'eeloped.
V. DEVELOPMENT OF ItETHODOLOCY
The construction of a reservoir in a river basin causes
water quality variation within the reservoir, which in turn modifies
the annual variation in the downstream reaches. Thermal stratifica
tion in the reservoir presents the opportunity for selectively with
drawing the quality of vater which best meets downstream require
ments. In this chapter, a methodology will be described which enables
the selection of the "best" water at successive time intervals. The
methodology is based upon the concept of a tradeoff among the layers
within the reservoir.
Descriotion of the Problem
The riverreservoir system is shovn in Figure 12. The
system consists of a reservoir upstream from a lare municipality,
which discharges its wastes into the stream. Water quality' standards
are established for the reaches downstream from the reservoir.
There are various alternative methods at the two control points,
reservoir and waste treatment plant, that can be employed to meet
the water quality standards. The objective is to determine which
alternative or combination of alternatives is the most economical.
This can be achieved by comparing the costs of various schemes.
80
C,
0/
J
/0
z a
! LU I
k I
G)
1 LLI
7>
g   K ?s /(.9
z & 
LLI^/
C,)/
\I7
Given any combination of alternatives (e.g., reservoir
mixing and waste treatment) and various constraints at the reser
voir (such as minimum releases for navigation and maximum flood
control storage), the problem is to minimize annual downstream waste
treatment costs while maintaining water quality standards and satisfy
ing the constraints. It is not sufficient to minimize downstream
treatment cost during each time interval because this policy would
release the best quality each month with no consideration to
the following months. Conceivably, this might lead to the point
that all water releases violate the stream standards. Nor is it
sufficient to release the minimum quality of water which, when com
bined with waste treatment, meets the stream standards. A policy' of
this nature conserves the higher quality water in the reservoir
instead of allocating the water among the various time intervals.
This policy also would encourage high levels of BOD removal in the
waste discharge and implies that even during severe water quality
periods the high quality water might remain in the reservoir. The
best policy is to allocate the water in the reservoir in such a
manner that annual downstream waste treatment is minimized. Since
the reservoir is thermally stratified during the heating cycle, and
most water quality problems occur during the summer and early fall,
the term "annual" will apply to the period from April to October.
Most reservoirs are only slightly stratified between Nov.ember and
March, with the result that there is no temperature tradeoff within
the reservoir.
To minimize annual waste treatment costs, an optimization
technique is needed which will enable the consideration of the tradeoff
between release temperature and flow. In each time period, decisions
must be made on the temperature of the release water and the amount
of water to release. Since the decisions made at each time interval
affect the state of the reservoir in the next time interval, the
technique must be able to incorporate the dynamic characteristics
of the reservoir. Before formulating the optimization technique,
discussion will be presented on the different techniques used in
water quality studies.
Water Quality Optimization Models
There are several features which all water quality optimi
zation models have in common. First, the models must be formulated
so that decisions can be translated into a ranking function or objec
tive function. The ranking function serves as a means for determining
the effectiveness of decisions, and should provide a welldefined
ordering of alternative decisions. The optimization technique is
applied to this function.
Second, the decisions are required to meet certain con
straints, such as maintaining water quality standards or releasing of
a minimum quantity of water from the reservoir at all times. Some
times the constraints are not imposed directly on the decision, but
rather on the consequences of the decision. For example, if a
decision were made on the quantity of waste discharge at a treatment
plant and the constraint uas a minimum dissolved oxygen level in the
stream, the decision would be disallowed if the response of the stream
to the waste violated the standard.
A third feature is the result of the constraints. Most water
quality models are formulated to investigate the effects of pollution
on the dissolved oxygen standard. Therefore, a predictive model for
the dissolved oxygen response to various waste loadings is necessary.
The model selected is usually the StreeterPhelps Equation.
Fourth an optimization technique is chosen which best
handles the formulated model. The technique should be amenable to
programming on a digital computer. In water quality studies, the
techniques fall into three categories: nonlinear programming, linear
programming, and dynamic programming. An example of the application
of each of the techniques will be presented.
Kerri 58 developed a model which might be used by regulatory
agencies to determine waste treatment requirements for several muni
cipalities and industries within a river basin. His objective func
tion was to minimize total waste treatment costs in the basin for
various dissolved oxygen standards. The degree of treatment required
at each treatment plant was determined at the optimal solution. He
used nonlinear programming to solve for the degree of treatment and
the StreeterPhelps Equation for predicting the dissolved oxygen
profile in the river. The model was tested on the Willamette river
basin for a critical low flow period to investigate several treatment
alternatives. The model provides a means for determining the economic
impact of setting uniform waste treatment standards. Kerri found that
in the Willamette basin, the minimum annual cost solution was
$2,999,000 for a dissolved oxygen standard of 5.0 mg/1. At the opti
mal solution, the muniicipalities were not required to provide waste
treatment, but four sulfite pulp mills were required to maintain 90'
EOD removal. By forcing all municipalities and industries to pro
vide 74% BOD removal, the cost increased to $4,733,000.
Liebman5 and P.eVelle et al.59 applied different optimizing
techniques to the same problem described by Kerri. The objective
function was that of minimizing total basin costs of waste treatment.
N
Liebman formulated the model as minimizing ci( i).. subject
i=l
m
to o s. (i=l,2,3,...,N): where pi is the percentage BOD removal
at the ith discharge point; om is the minimum dissolved oxygen
concentration in reach i resulting from the ih discharge of waste
and all preceding waste discharges; si is the minimum allowable
dissolved oxygen concentration in reach i. The function ci(pi) is
the total annual cost of providing pi percentage BOD removal at the
ith treatment plant. The function can be expressed in mathematical
form or in tabular form. Liebman used the StreeterPhelps Equation
to predict the dissolved oxygen response in the stream. He used
dynamic programming to optimally allocate the dissolved oxygen resources
in the stream. The model was tested on the Willamette River Basin and
yielded a minimum cost similar to that of Kerri; however, in Liebman's
study all municipalities were required to furnish primary treatment
of 35% BOD removal. In Liebman's model, the effects of changing the
dissolved oxygen standard in the downstream direction can be studied
permitting monetary evaluation of the standards in each reach.
ReVelle et al. assumed a linear relationship between organic
loading and dissolved oxygen response. They required the cost func
tion to be linear and minimized the total waste treatment cost in a
85
basin subject to the inequality constraints of the water quality
standards (e.g., dissolved oxygen required to be greater than or
equal to the minimum allowable dissolved oxygen) in each reach. The
StreeterPhelps Equation was linearized and formulated as equality
constraints in each reach. They used linear programming to solve
for the optimal solution. When solving the linear programming algo
rithm, the reach in which the water quality constraint is restrictive
can be found in the solution to the dual. For each water quality
constraint, a dual variable exists which will be nonzero if the con
straint is binding at the optimum. By relaxing that constraint, the
minimum cost can be decreased. ReVelle et al. applied the linear
programming technique to the Willamette River Basin and achieved
approximately the same results as Liebman; however, they were able
to identify the reach in which water quality standards could be
lowered to reduce the overall treatment costs.
The water quality models described above may aid the regula
tory agency in setting stream standards. In each model, the objective
is to lower waste treatment costs for the entire basin. The optimal
efficiency is determined for each plant to satisfy the standards
during the lowflow period. The low flow usually selected for sewage
treatment plant design is the lowest consecutive sevenday or tenday
flow which has a frequency of once in tenyears. To construct a
plant based on this lowflow criterion means that the plant will
actually "over treat" for the stream a large part of the time.
Moreau 10 suggests the use of a steadystate treatment plant and an
auxiliary unit which is used during lowflow periods for providing
sufficient BOD removal to meet stream standards. He showed that
the combined cost of the steadystate and auxiliary units was less
than the cost of a single unit designed to satisfy the louflo\
criterion.
Davis investigated the use of advanced vastc treatment
during critical flow periods as a substitute for louflow augmenta
tion. He found that in the Potomac Estuarv advanced waste treatment
was approximately the same cost as lov flow augmentation. Advanced
waste treatment provided more fle::ibilit, since the processes could
be operated only during the critical period. He also showed that
operating costs can be evaluated vis a vis capital costs in the
selection of an advanced 7aste treatment process. In his example,
step aeration plus micros training had a capital cost of $20,000,000
and an annual operating cost of $56,000 while chemical precipitation
had a capital cost of $200,000 and an annual operating cost of
$426,000. If the tvo processes are compared at an annual operation of
3.5 months, a 50 ,year project life and a i4' interest rate, the present
value of treatment costs are $33,700,000 for the step aeration plus
microstraining and $32,300,000 for chemical precipitation. At higher
interest rates chemical precipitation becomes more favorable.
For many municipalities, the establishment of stream stand
ards will require the construction of additional treatment plants.
In view of the work of iloreau and Davis, the addition of a lo', capital
cost, high operating cost auxiliary unit may be the most economical
alternative. If auxiliary. units w.ere installed, the objective could
be changed to minimizing the operating cost of the au:x:iliary treatment
unit.
60
Little formulated a model for optimizing reservoir
releases for hydroelectric power generation, which is similar to
minimizing auxiliary treatment costs. Since future river flows are
unknown and the release decision is made under uncertainty, he chose
to maximize the expected return, not the return itself. Little
considered the case of a hydroelectric dam and a steam plant operating
conjunctively to produce electric power for a region. The steam
plant was used to supply the demands which water could not satisfy.
His objective function was minimizing expected supplemental costs from
the steamelectric plant. He used stochastic dynamic programming to
determine the optimum releases for minimizing the objective function.
He assumed that the inflows to the reservoir were random variables
which could be represented by a Markov model of lag one. The optimi
zation technique required that the problem have a fixed end condition
from which to start. The expected value of the cost function was then
computed using the probability distribution of inflows for each time
interval beginning at the known point and working back to the present.
This technique is called backward looking stochastic dynamic programming.
61
Young commented on the various dynamic programming tech
niques and developed what is termed Monte Carlo Dynamic Programming.
Young used a forward looking dynamic programming model to determine
reservoir operating rules in flood control reservoirs. He assumed the
annual inflow to be known and then optimized the release quantity.
By generating many years of synthetic flows and optimizing the draft
for every additional flow, Young could find the optimal release rule.
Young also developed another method of determining the release rule.
ly assuming the optimal draft to be a linear function of storage,
inflow, and an inflow forecast, the optimal draft can be determined
for each time interval by making a regression analysis on the selected
variables.
In reviewing the previous models, one observes that dynamic
programming :as used for optimizing releases from reservoirs. This
technique is particularly useful for nonlinear cost functions, for
problems involving two or more decision variables, and for a system
which changes with time. The problem of interest in this study is
that of selecting the release water temperature and flow which mini
mizes downstream treatment costs during the summer months, Since
this involves two decision variables and a reservoir which changes
in accordance with the decisions, dynamic programming will be selected
as the optimizing technique. The formulation of the algorithm will
be presented in the ne.:t section.
Dynamic Programming
There are several general features of dynamic programming
which male the procedure applicable to water resource studies. The
technique is especially useful in determining the allocation of a
limited resource among various uses. In the problem at hand, it is
necessary to allocate water from a reservoir for succeeding time
intervals. Dynamic programming is very useful for multiple decision
problems which change over time in accordance with decisions made in
earlier time intervals. The dynamic nature of a reservoir is apparent
when one considers the marked effect releasing water from the various
layers has on the state of the reservoir. Another important feature
of dynamic programming is that nonlinear cost functions can be used.
This feature is particularly useful if the cost function itself
changes with time. If several constraints are required in the formu
lation of a dynamic programming algorithm, the constraints tend to
simplify the problem instead of making it more complex. This point
will be brought out in the development of the algorithm.
To facilitate the derivation of the algorithm, a one 
dimensional decision process will be considered. It is desired to
allocate a known quantity, Q, of water over several time periods, N.
The water is to be allocated in such a manner that the total cost, C,
of auxiliary waste treatment is minimized over the intervals. The
cost, ci, in the ith period is some function of the water, qi, re
leased in that period:
ci = gi(qi) ........................ (71)
The objective function is
N N
minimize C(ql,q2, ...qN) = Z ci = gi(qi) ..... (72)
i=l i=l
The minimization is subject to the constraint that the sum of the
qi's must be less than or equal to the total resource Q. This is
expressed as:
N
Z qi Q .......................... (73)
i=l
The sequence of functions fn(q) is introduced as follows:
fN(q) = min C(ql,q2, ... .qt) ........ (74)
{qil
The function fII(q) is the optimal solution for an allocation of the
quantity 0 of resource to N time intervals. It is seen that when N
is 1, fn(q) is simply the minimum of Equation (71) for i=l; i.e.,
fl(q) = g1(q) (75)
To develop the dynamic programming algorithm, it is necessary to
find a recurrence relationship connecting f1i(q) with fj_l(q) for an
arbitrary N and resource Q. If q, is the amount of resource allocated
to period II, then, regardless of the cost of qIj, the remaining re
source Qq!l is used to obtain the minimum cost in the remaining !Ii
periods. By definition, the minimum cost of allocating Qqll resource
to 11 periods is the function f,,_! (q,) An initial allocation of
qN. to period II has the cost:
N = gN (q 1 ) + f l (Q q **...........* (76)
The optimal allocation of ql, resource to the 1th period is that
allocation whichh minimizes Equation (76). The basic functional equa
tion can be expressed as:
f (Q) = min {g (qi)+rf_l( ) ... (77)
In the derivation of Equaition (77), a general principle
vas used: The Principle of Optiinality. This principle expressed
by Eellman 6 is: "An optimal policy has the property that whatever
the initial state and initial decision are, the remaining decisions
must constitute an optimal policy with regard to the state resulting
from the first decision."
To make the allocation of the resource Q amenable to digital
programmiing, the resource is divided into discrete units. Also, to
male the problem realistic, constraints are placed on the possible
states of Q. Such constraints might be a minimum release, qnmi' whichh
must be maintained at all times for navigation, and a maxitmum release,
qmax which might be the physical limit of the reservoir release
structures. These constraints change Equation (77) into the form:
f.m(Q) = min {gi (q )+f _l (QqI } ............. (78)
q min qqmax
Example Problem
To illustrate the use of the algorithm, an example alloca
tion problem is hypothesized. It is desired to allocate ten units
of water over four time periods in such a manner as to minimize the
total cost of treating waste to maintain water quality standards. A
minimum discharge of one unit must be maintained during each interval,
and the maximum discharge is four units. The cost of waste treatment
as a function of units of discharge is given in Table 3.
TABLE 3
COST LIATRIX FOR DYNiIIC PROGPA.1 II1IG EXAMPLE
Costs for Each Period
Units Released 1 2 3 4
1 3 4 5 7
2 2 2 3 4
3 0 1 1 1
4 0 0 0 0
Costs for waste treatment during the first time interval are
given by Equation (75). These values are simply the costs as given
for period one in Table 3. For period two, the function f2(Q) is
developed by finding the minimum cost of allocating from two to eight
units to the first two periods. A minimum of two units in the result
of the constraint that one unit must be released at all times and
