Group Title: methodology for selecting among water quality alternatives
Title: A Methodology for selecting among water quality alternatives
Full Citation
Permanent Link:
 Material Information
Title: A Methodology for selecting among water quality alternatives
Physical Description: Book
Language: English
Creator: Nicolson, Gilbert S., 1941-
Copyright Date: 1969
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Record Information
Bibliographic ID: UF00097770
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001027754
oclc - 09367602
notis - AFA9773


This item has the following downloads:

PDF ( 7 MBs ) ( PDF )

Full Text









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.




ACKNOWLEDGMENTS .............................................. ii

LIST OF TABLES ................................................ vi

LIST OF FIGURES ............................................... vii

ABSTRACT ...................................................... viii


I. INTRODUCTION AND OBJECTIVES ......................... 1

Introduction ........................................ 1
Historical Background ............................... 1
The Parameter of Major Interest ..................... 3
Water Quality Alternatives ............. ..... ... .. 5
Objectives .......................................... 7


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


Stream Temperature ................................... 51
Dissolved Oxygen .................................. 56
Determination of Allowable Waste Loading ............ 69

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


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


A. PROBABILISTIC DEFICIT FUNCTION ........................ 145

B. PROGRAM LISTING ........................................ 150

C. GLOSSARY OF PROGRAMMING TERMS ......................... 176

REFERENCES ................................ ..................... 180

BIOGRAPHICAL SKETCH ............................................. 185


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


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 Coefficie-nt 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. River-Reservoir System .................................... 80

13. Dynamic Progra inming Flo.-' Chart ............................. 98

14. Savannah River Basin ...................................... 116

15. P.cer-voir 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.



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 low-flow 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 low-flow periods. Alternatives at the

reservoir site include multi-level 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 least-cost alternative.


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 two-dimensional 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 obtaining-reliable cost comparison.



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 and-water 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 water-borne

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 benefit-cost

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 least-cost 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 least-cost 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


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 in-stream 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

low-flow periods has become a method of water quality control.

Bramhall and Mills7 criticized the use of low-flow 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 low-flow augmentation,

which would introduce a bias in favor of waste treatment. Young

pointed out the compatibility of flood control storage and low-flow

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 steady-state and stochastic units. The

effluent from the steady-state 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.


In planning to meet water quality standards, all alternative

schemes and combinations should be considered in attempting to find

the least-cost 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.


To develop alternative plans for meeting water quality standards

in a reservoir-river 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;

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 least-squares analysis

of the scream temperature data. Ward calculated these coefficients

for several rivers in Arkansas. Figure 1 illustrates the effect of



z W w
\w 4

+ -


z I '-

0 0
Co \ o s

Jo '3Nn1V8J]dLN31
! / l t
N / S =


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 forty-day 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-

and-after 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 short-wave solar radiation, H ; incoming long-wave radia-

tion H ; reflected short-wave and long-wave radiation, Hsr and Har; long-
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

relative magnitudes. (For a more detailed discussion, the reader is

referred to Edinger and Geyerl4 or to the Report of the 1960-61

Advanced Seminar of the Johns Hopkins University, Department of

Sanitary Engineering and Water Resources.15)

Absorbed Radiation

Short-wave 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 short-wave radia-

tion, but it is measured more easily than it is computed. The U.S.

Weather Bureau records short-wave 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

Long-wave atmospheric radiation is a function of cloud cover,

air temperature and air vapor pressure. Raphael presents the follow-

ing formula for long-wave 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 Stefan-Boltzman constant (4.15 x 10-8

Btu per square foot per degrees Rankine per day); and Ta is the air

temperature in degrees Fahrenheit. Long-wave radiation ranges from

2400 3200 Btu per square foot per day.

Reflected long and short-a'.'ave radiation are calculated as fractions

of incident long and short-wa: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


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 Stefan-Boltzman 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

vapor pressure of water in millimeters of mercury, mm Hg, determined

from the water surface temperature, T ; and e is the air-vapor 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)

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 at-which 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 wind-induced 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 three-dimensional 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
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




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 x-y plane in the positive z-direction


Ocp T OPp vx T 2 pc vx Sx

z = z+pz z = z

where v is velocity in the x-disection in feet per day.

The heat diffused across the x-y 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 x-direction 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 z-directions, 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

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








o I-




o w
0 Cw

0 0

M w

w cr
F- J


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, short-wave 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
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 be gr-atest

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 "turn-over" takes place.













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

"turn-over" 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

short-wave radiation is absorbed in the epilimnion and transfer of

heat through long-wave 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

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







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 low-level 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 this-low 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 multi-level 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



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 cross-sectional 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 ten-foot 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 twenty-five feet in July

and August. This layer would correspond to the epilimnion and would

have the same temperature throughout.

Second, the temperature of the ten-foot 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


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

interval. Data averaged over a seven-day 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


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




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)(eE-ea) + .26(a+bW)(TE-Ta) ........ (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)(e-eE)

+ (a+bW) 0.26 (Ts-TE)} ............................... (20)

Edinger and Geyer use a linear relationship to approximate

saturation vapor pressure and water temperature For ten-degree Fahren-

heit temperature ranges. The constants for the temperature increments

are shown in Table 1.



Temperature Range, OF B, mm-Hg/OF Intercept C (8)mm-Hg

40-50 .2910 -5.47

50-60 .4050 -11.22

60-70 .5553 -20.15

70-80 .6667 -27.80

80-90 .9900 -53.33

90-100 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:


H = [15.71 (T-TE) + .051(Ts2-TE2) + (a+bW) (Ts-TE)

+ 0.26 (a+bW)(Ts-TE) ... ........ ..................... (23)

Combining like terms yields:

AH = [15.7 + (0.26+8)(a+h'] (Ts-TE) +.051(T-T2) ..... (24)

From Equation (24), Edinger and Geyer define the exchange

coefficient, K in Btu per square feet per day per degree Fahrenheit,


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


eE = TE + C ( ) ......................................... (27)

Edinger and Geyer reduced Equation (19) to:

TE + .051TE2 = HR-1801 + K-15.7 [ea-C(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


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
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 wind-induced turbulence. The rate of oxygen absorption

is proportional to the difference between the saturation level and

the dissolved oxygen present:

dt K2 (Cs C) ................................ (29)

where C is the dissolved oxygen in milligrams per liter (mg/l); K is

the reaeration coefficient in days-l; 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 day-1 (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

saturation-dissolved 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 lower-layer 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


I '
10 -


30 -

40 JUNE 30,1967-
\\ \ AUGUST 16,1967---i
50 DECEMBER 7, 1967--.-c


S70 -

8 90_

100 -

S110 4


140 / 0




JUNE 3017 6--

AUGUST 16,107
DECEMBER 7,1 -----G




with TIME and DEPTH




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


Improving Release Water Ouality

Kittrell discussed several alternative methods for improving water

quality in released water. The most promising methods were multi-level

penstock intakes, turbine reaeration, and reservoir mixing.

Multi-level intakes permit the selection of that quality of water

which best meets downstream needs. The installation of multi-level

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 multi-level 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 steam-electric 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 low-flow 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


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.


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.

t=l a11 12 13

2 a21 22 23

3 a31 32 33


a a a
n anl n2 n3


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


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



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 reservoir-river 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.


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

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 x-direction is along the length of the stream, the z-direction

is vertical, and the y-direction 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


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 steady-state,

-T = A condition of this nature will negate diurnal temperature
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
f S dz is the net rate of heat transfer, AH. The term vx dz can
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)
where d will be taken as the mean depth of the stream. If in Equation
(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 cross-sectional 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) e-r2 + 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 steady-state 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 steady-state

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


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:

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 days-1

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 Streetcr-Phelps 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 days-1. 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
K9 = 3.3 1.33 ..................................... (46)

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:

dt K1 L K2 D .................................. (47)

Integrating Equation (48) yields the classical Streeter-Phelps


K1 La -K1t -K2t -K2t
D K La (e- t e ) + Da eK2 ............ (48)
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

0 U

I- IaJ
S I-

1/5W u! N39),XO 03A10SSIG

river to river and for different wastes, and estimates have ranged
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 day-1 and a standard deviation of 0.066 day-1

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)

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

load--for 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 Streeter-Phelps 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 Streeter-Phelps

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; (P-R), 0.55 to 2.9 mg/l/day; and K3, 0.07 to

0.11 day-1.

In many instances, one cannot neglect the terms represented in

the modification to the Streeter-Phelps, 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 Streeter-Phelps

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

dt = -K L + r ................................. (52)

dt = KIL-K2 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)

B = Csf Csi ..................................... (55)

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:

dt = K1 L + r ..................................... (56)

dt K1 I K2 D B + s ............................ (57)

Integration of Equations (56) and (57) for steady, uniform flow


L = (La- r) e- t + _r ............................... (58)
K1 K1

D K1La r -t _-+(Da BE) -Kt r+s ) t e-K 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

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+

(aD)2 o 2 + (;D)2 a 2
r r s-- s





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 days-1. 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) i-1 i- -
C r r- r-~ 00 r, 0h ON m 0%.

0 U

U <

] p u L


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)


1\0 C]
I' -,"

'-. C0

-- I- 0 r- Co c') co r---
N N c 1 r-- N H 1--

iL N- 0 O'. -. r-1 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') 4--4 cN -j -r -. C^

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~


0-. -. rN-

r'. C., CO

-0 %D0 0

Lr) 0". 'C'

'- N"- 0--.

L)n L'i N

ILn %0 C


.4 C')

Cf m mC C) CI C')
N- 'o. r- N- r^- el


L1 Li)
C'I C`]

-. ui


0 0

0 c

a E



o o



r? n)




pa H
< i-


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


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.


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

two-fold. 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
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
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


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


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

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.

Low-flow 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


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 low-flow 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.






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, in-stream 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.


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 river-reservoir system is shovn in Figure 12. The

system consists of a reservoir upstream from a lar-e 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.




z a

!- LU I-

k I-

1- LLI
g - - K ?s /(.9
z &- -

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 well-defined

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 Streeter-Phelps 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 Streeter-Phelps 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.
Liebman formulated the model as minimizing ci( i).. subject
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 Streeter-Phelps 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


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

Streeter-Phelps 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 non-zero 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 low-flow period. The low flow usually selected for sewage

treatment plant design is the lowest consecutive seven-day or ten-day

flow which has a frequency of once in ten-years. To construct a

plant based on this low-flow 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 steady-state treatment plant and an

auxiliary unit which is used during low-flow periods for providing

sufficient BOD removal to meet stream standards. He showed that

the combined cost of the steady-state and auxiliary units was less

than the cost of a single unit designed to satisfy the lou-flo\


Davis investigated the use of advanced vastc treatment

during critical flow periods as a substitute for lou-flow 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


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 steam-electric 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.
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


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
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:

Z qi Q .......................... (73)

The sequence of functions fn(q) is introduced as follows:

fN(q) = min C(ql,q2, ... .qt) ........ (74)

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 Q-q!l is used to obtain the minimum cost in the remaining !I-i

periods. By definition, the minimum cost of allocating Q-qll resource

to 1-1 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 (Q-qI } ............. (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.



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

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs