WATER POLLUTION CONTROL RESEARCH SERIES • 16130EXT12/69
MATHEMATICAL MODELS
FOR THE PREDICTION OF
THERMAL ENERGY CHANGES
IN IMPOUNDMENTS
ENVIRONMENTAL PROTECTION AGENCY • WATER QUALITY OFFICE
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes the
results and progress in the control and abatement of pollu-
tion of our Nation's waters. They provide a central source
of information on the research, development, and demon-
stration activities of the Water Quality Office, Environ-
mental Protection Agency, through inhouse research and grants
and contracts with Federal, State, and local agencies, re-
search institutions, and industrial organizations.
Inquiries pertaining to the Water Pollution Control Research
Reports should be directed to the Head, Project Reports
System, Office of Research and Development, Water Quality
Office, Environmental Protection Agency, Washington, B.C. 20242.
-------
MATHEMATICAL MODELS FOR THE PREDICTION OF THERMAL
ENERGY CHANGES IN IMPOUNDMENTS
by
Water Resources Engineers, Inc.
Walnut Creek, California
for the
WATER QUALITY OFFICE
ENVIRONMENTAL PROTECTION AGENCY
Project #16130 EXT
Contract # 14-12-422
December 1969
For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C., 20402 - Price $1.50
-------
EPA Review Notice
This report has been reviewed by the Water Quality Office,
EPA, and approved for publication. Approval does not signi-
fy that the contents necessarily reflect the views and poli-
cies of the Environmental Protection Agency, nor does mention
of trade names or commercial products constitute endorsement
or recommendation for use.
-------
TABLE OF CONTENTS
Page
LIST OF FIGURES v
LIST OF TABLES viii
I INTRODUCTION 1
BACKGROUND AND NEED 1
OBJECTIVES 6
ORGANIZATION OF THE RESEARCH EFFORT 8
GENERAL APPROACH 8
REPORT FORMAT 9
AUTHORITY AND RESPONSIBILITY 11
II CLASSIFICATION OF RESERVOIRS 13
PART I -- BASE SIMULATION STUDIES
III INITIAL MODEL APPLICATIONS 19
HUNGRY HORSE RESERVOIR 21
OBJECTIVES 22
MINOR MODEL REFINEMENTS 22
IV INTERNAL MIXING STUDIES 29
-------
MIXING PROPERTIES OF RESERVOIRS 30
EFFECT OF VARYING THE FUNCTIONAL REPRESEN-
TATION OF THE INTERNAL MIXING PROCESS 39
EXPONENTIAL REPRESENTATION 40
STEP FUNCTION 46
REPRESENTATION OF A(z,t) BY A
CONTINUOUS FUNCTION 50
PROPOSED FUNCTIONAL FORM FOR THE
EDDY CONDUCTIVITY COEFFICIENT 60
MODEL SENSITIVITY TO THE DIFFUSION
FUNCTION 63
VERIFICATION OF THE DEEP RESERVOIR MODEL 67
V SELECTIVE WITHDRAWAL STUDIES 75
THEORY 75
SIMULATION OF HUNGRY HORSE RESERVOIR
USING SELECTIVE WITHDRAWAL THEORY 85
VI MULTIPLE OUTLET STUDIES 93
TEST CASE I 94
TEST CASE II 98
VII SUMMARY AND CONCLUSIONS -- PART I 103
GENERAL 103
LIMITATIONS 104
RECOMMENDATIONS 106
-------
PART II -- SIMULATION OF WEAKLY STRATIFIED RESERVOIRS
VIII INTRODUCTION 109
IX SEGMENTING THE RESERVOIR 113
NOMENCLATURE 113
PROCEDURE 114
GEOMETRICAL CONSIDERATIONS 114
RESIDENCE Tim 115
TRIBUTARY INFLOWS 116
LONGITUDINAL RESOLUTION OF ISOTHERM 116
X INTERFACING ADJACENT RESERVOIR SEGMENTS 119
MOMENTUM AND BUOYANCY CONSIDERATIONS 120
INTERFLOW THICKNESS 121
MASS CONTINUITY 122
HEAT CONTINUITY 123
XI APPLICATION TO LAKE ROOSEVELT 127
GENERAL INFORMATION 127
PHYSICAL AND THERMAL DATA 128
METEOROLOGIC DATA 130
SEGMENTING THE RESERVOIR 133
SIMULATION RESULTS 138
GENERAL OBSERVATIONS 138
SPECIFIC OBSERVATIONS 144
REMARKS 149
111
-------
XII SUMMARY AND CONCLUSIONS -- PART II 151
GENERAL 151
RECOMMENDATIONS 152
CAVEAT 154
XIII LIST OF REFERENCES 155
IV
-------
LIST OF FIGURES
Page
FIGURE I The Columbia River Basin 4
FIGURE 2 Hungry Horse Reservoir 20
FIGURE 3 Observed Temperature Profiles,
Hungry Horse Reservoir, 1965 23, 24
FIGURE 4 Effective Diffusion and Temperature
Profiles for Selected Impoundments 32
FIGURE 5 Energy Distribution 34
FIGURE 6 Contribution to Vertical Heat Flux
by Individual Heat Transfer Mechanisms,
Hungry Horse Reservoir, 1965 41
FIGURE 7 Typical Computed Temperature Profile
with A(z,t) Described by Equation 10
and A(zT,t) - 5.0 x 10"5 m2 sec"1 44
FIGURE 8 Typical Computer Temperature Profile
with A(z,t) Described by Equation 10
and A(zT,t) - 1.0 x 10"4 m2 sec"1 44
FIGURE 9 Representation of A(z,t) with a Step
Function 47
FIGURE 10 Typical Computed Temperature Profile
with A(z,t) Described by Equation 11;
Simulation Day is 1 September 1965 50
FIGURE 11 Average Eddy Conductivity Profile for
Hungry Horse Reservoir, 1965 52
FIGURE 12 Log of Eddy Conductivity versus log
Stability -- Hungry Horse Data 54, 55
FIGURE 13 Typical Temperature and Eddy Conductivity
Profile Obtained by Use of Continuous
Function for A(z,t) 56
-------
FIGURE 14
FIGURE 15
FIGURE 16
FIGURE 17
FIGURE 18
FIGURE 19
FIGURE 20
FIGURE 21
FIGURE 22
FIGURE 23
FIGURE 24
FIGURE 25
FIGURE 26
Relationship between Hypolimnetic Value
of Eddy Conductivity Coefficient and
Deep Reservoir Hydraulics
Effect of Eddy Conductivity Parameter
Variation on Simulated Reservoir
Temperature Profiles; Hungry Horse
Reservoir, 1965
Reservoir Simulation with A(z,t)
Defined by Equation 12; Hungry
Horse Reservoir, 1965
Simulated and Observed Thermal Regime
for Hungry Horse Reservoir, 1965
Breakdown of Outlet Discharges and
Computed Outlet Temperatures for
Hungry Horse Reservoir, 1965
Simulated and Observed Temperature
Profiles for Hungry Horse Reservoir,
1965
Schematic Sketch of Flow Separation
in Stratified Flow
Schematic Diagram of the Application
of Selective Withdrawal Theory to a
Stratified Impoundment
Withdrawal Region for Case of Dividing
Streamline Falling Below Reservoir
Bottom
Withdrawal Region for Case of Computed
Dividing Streamline Crossing the
Thermocli ne
Effect of Convective Mixing on
Thermal Profile
the
Assumed Reservoir Density Structure
for Application of Craya's Criteria
Compa,rision of Outflow Temperatures
Computed from Debler's Criteria with
Measured Cooling Water Temperatures
(Fontana Reservoir)
59
64
66
68
69
70, 71
76
77
80
81
82
83
86
-------
FIGURE 27 Simulated Thermal Regime of Hungry
Horse Reservoir, 1965, Using Selective
Withdrawal Theory
FIGURE 28 Simulated Discharge Temperatures for
Individual Outlets as Computed with
Selective Withdrawal Theory; Hungry
Horse Reservoir, 1965
FIGURE 29 Simulated Thermal Regime of Hungry
Horse Reservoir, 1965, with Modified
Reservoir Release Scheme
FIGURE 30 Discharges and Discharge Temperatures
for Individual Outlets with Modified
Reservoir Release Scheme, Hungry Horse
Reservoir, 1965
FIGURE 31 Simulated Thermal Regime of Hungry
Horse Reservoir, 1965, with Addition
of Additional Turbine Intake; Hungry
Horse Reservoir, 1965
FIGURE 32 Discharges and Discharge Temperatures
for Individual Outlets with Addition
of Additional Turbine Intake; Hungry
Horse Reservoir, 1965
FIGURE 33 Schematic Illustration of a Reservoir
Divided into Segments and Elements
FIGURE 34 Relationship of the Interflow Properties
of the Model to Those of the Prototype
Reservoir
FIGURE 35 Plan View of Lake Roosevelt
FIGURE 36 Reservoir Segments in Plan View and
Typical Cross Sections, Lake Roosevelt
FIGURE 37 Reservoir Segments in Profile View,
Lake Roosevelt
FIGURE 38 Typical Results from Simulation of
the Thermal Regime of Lake Roosevelt,1967
FIGURE 39 Typical Shape of Computed Temperature
Profile during Cooling Period when
Convective Mixing is Suppressed
FIGURE 40 Effect of Linearization on the Temperature
Distribution of the Interflow
89
90
95
97
100
101
113
126
129
134
135
139, 140, 141, 142
146
147
VI1
-------
LIST OF TABLES
Page
TABLE 1 Physical and Meteorologic
Information, Hungry Horse
Reservoir 21
TABLE 2 Hungry Horse Reservoir --
Turbine Intake Temperatures;
and Detroit Reservoir --
Tail race Temperatures 87
TABLE 3 Physical and Thermal Data,
Lake Roosevelt 131 , 132
TABLE 4 Residence Times for Individual
Segments in Lake Roosevelt 136
-------
I. INTRODUCTION
BACKGROUND AND NEED
Water temperature has assumed an important role in water
resources planning and management because of its wide-ranging ef-
fects on aquatic flora and fauna and its influence on the physical
and chemical characteristics of the water itself. Some of the
better known effects of water temperature include: effects on fish,
particularly on various parts of their life cycle; effects on the
temporal and spacial variation in the dissolved oxygen content of
the water; effects on the toxicity level of various pollutants; and
effects on the taste and odor of drinking water. Details of thermal
effects can be found elsewhere (1, 2, 3, 4, 5, 6, 7).
In addition to the above mentioned effects, the water
temperature, through its influence on density, affects the hydro-
dynamic characteristics of the watercourse, particularly in impound-
ments which are subject to thermal stratification. Through its
influence on the hydrodynamic characteristics, the temperature can
exert a profound influence on the temporal and spacial variation of
water quality characteristics.
-------
Since the influence of temperature on the quality character-
istics of a watercourse is so broad, effective water quality planning
requires that there exist the capability to predict the temporal and
spacial variation of the water temperature under alternative develop-
ment plans in order to protect the needs of the water users within
the watercourse system. However, this is not a simple problem because
of the complexity of the mechanisms which govern heat transfer within
a waterbody and heat transfer between the waterbody and its surround-
ing environment. Before the time of World War I, it had been established
through the experiences of the industrial complexes in the eastern
and midwestern states that cooling water discharges could significantly
alter the thermal regime of a river (8); but it has also been recog-
nized, more recently, that a complex of reservoirs and the way in
which the water is discharged from them can have significant effects
on the thermal regime of the watercourse (9, 10).
In several areas of the country, the build-up of a complex
of dams within a watershed coupled with cooling water discharges from
industrial and thermal power plant cooling water discharges have
changed the thermal regime of the river to such an extent that the
native fish populations have either been eradicated, or their loss
is impending. The Columbia River Basin in the Pacific Northwest is
a typical case in which the existing complex of dams coupled with
the discharge of cooling water from the Hanford nuclear power plant
has altered the thermal regime of the Lower Columbia River to such
an extent that uncontrolled future development may very possibly
annihilate the existing anadromous fish run in the River.
- 2 -
-------
The Columbia River Basin, shown in Figure 1, encompasses a
drainage area of some 259,000 square miles, 15 percent of which is
in Canada. Thus the water resources management problems of the Basin
are international in scope. In contrast to the Columbia River in
the United States, which is highly developed, the river in its upper
reaches in Canada is almost as wild and untamed as it was when the
first white explorers saw it. Within the United States the Colubmia
River below Grand Coulee, and the Lower Snake River, are used pri-
marily for hydropower and navigation. In addition, this portion
of the basin supports a sizable commercial and sports fishery. Other
dams in the system, including Grand Coulee Dam, are used primarily
for hydropower, irrigation, and flood control and storage.
The present power system in the Pacific Northwest is al-
most all hydro. Consequently, the complex of dams in the Columbia
Basin is operated in a manner such that the base power demand plus
the peaking power needs are met. This type of operation, however,
is nearing an end. Future load growth in the area is expected to
triple in the next 20 years and, while hydro peaking capacity can
be increased at existing and planned projects, there are only very
limited opportunities to develop new projects that would provide
significant amounts of additional hydro-energy (11). Thus, future
plans call for the construction of thermal nuclear plants to meet
the base load demands with the hydro projects being used more and
more for peaking power only. According to a recent report (12),
practically all of the Columbia River mainstem above the Dalles Dam
is considered as preferred nuclear siting areas.
- 3 -
-------
LEGEND
MICA- 1,820 (storage)
DOWN1E CREEK- \,000-r
REVELSTROKE CANYON-630^
CANADA
WASHINGTON
U.S.
ARROW LAKES
(storage)—tt
.MURPHY CREEK 30_0
BOUNDARY-
BOX CANYON-60
NON
FED FED
5 I EXISTING OR UNDER CONSTRUCTION
^ § AUTHORIZED OR LICENSED
5 Q POTENTIAL SITE
NAMEPLATE RATING IN MEGAWATTS
UNCAN^
LAKE
(storage)
.CORRA LINN-4o{(storaqe)
.UPPER BONNINGTON -64
LOWER BONNINGTON -47
S. SLOGAN-47 t.
KOOTENAY CANAL,r270
BRILLIANT-82
CHIEF
CHELA*N-48
ROCKY REACH-7l2«,i
ROCK ISLAND-212
.•
COULEE
1944
MONTANA
HUNGRY HORSE-285
NOXON
RAPIDS
283
FALLS ^THOMPSON
II ^--
(storage) 30
^
KERR-168
BUFFALO RAPIDS No.2-120 I
BUFFALO RAPIDS No.4-120)
KNOWLES-256
MONUMENTAL. LITTLE GOOSE >405
405. _^_ ^,1 n\«/PD GRANITE-405
_W£NAPUM-83I
PRIEST RAPIDS-788
BEN FANKLIN-352-^
ASOTIN^2TO
CHINA GARDENS-180 r~
HIGH MOUNTAIN SHEEP--
McNARY-980
Columbia River
PNPC-1,225 USBR-I.5OO-,
-•SQO \
HELLS CANYON
OXBOW-190
BROWNLEE-360
SCALE IN
FIGURE 1. The Columbia River Basin
- 4 -
-------
As mentioned above, the present complex of dams in the
Columbia Basin, with its present scheme of operation plus the dis-
charge of cooling water from the Hanford nuclear power plant, has
already altered the temperature regime of the Lower Columbia suffi-
ciently to endanger the existing fisheries. By the construction
of additional dams in Canada and the United States under the Columbia
River Treaty (13), the possible addition of nuclear power plants
on the mainstem, and the alteration of the operation of the Columbia
system, there is a real possibility that the existing fisheries in
the Columbia will be lost unless temperature control is integrated
into the future planning and development of the Basin.
Recognizing the thermal problems attendant with the present
and future development of the Columbia River System, the Federal
Water Pollution Control Administration (FWPCA) has undertaken the
Columbia River Thermal Effects Project in order to determine the
"impact of hydroelectric projects and thermal loads on water tempera-
tures in the system from the Canadian border to the river mouth at
Astoria. Where appropriate, the Study will utilize established
methodology or proven tools to develop the facility for prediction
of temperature changes under various climatological, meteorological,
or operational conditions. In a number of instances, it will be
necessary to direct concurrent research efforts to provide the needed
capability.
It has been determined that a mathematical model or "simu-
lation package" for prediction of temperature changes in the river
system is an essential tool for the Study. This model or set of
models must have the capability to represent with reasonable
accuracy the thermal responses of each of the major types of hydraulic
- 5 -
-------
regimes in the system from comparatively simple uniform flow in
river reaches to the more complex stratified flows in reservoirs.
Capability is also needed to cope with highly unsteady flows asso-
ciated with power load variations or with tidal influences in the
river's estuary.
As an initial step in its program,the Study Team has
undertaken the development of a "run-of-the-river" temperature model.
This model is expected to be suitable for prediction of temperature
changes in reaches of the river system unaffected hydraulically by
the backwater of deeper impoundments and in certain of the more
shallow reservoirs which may be considered to be well-mixed verti-
cally. However, in order to provide more flexibility for prediction,
including those cases where flows are partially or strongly strati-
fied, the Study Team wishes to incorporate into its model package
certain components of a deep reservoir model developed by Water
Resources Engineers, Inc. (WRE). Further, it wishes that the model
be refined and restructured as necessary, so that it can cope with
the full range of practical problems to be studied in the many compli
cated reservoirs identified with the Columbia River System. To
this end, FWPCA has contracted with WRE to perform a research and
development study whose scope is sufficient to meet the required
objectives as stated below.
OBJECTIVES
A general objective of the research reported herein was to
devise a generalized mathematical model which would represent, with-
in the practical and reasonable limits of accuracy, the thermal
- 6 -
-------
changes which may be expected under alternative hydrologic, hydraulic
and climatologic conditions for operating reservoirs. It is in-
tended that this work support and complement that already completed
or underway by the Study Team in its investigation of the Columbia
River System.
Specific objectives of the work were to provide for the
Columbia River Thermal Effects Project the means for simulation of:
I. thermal energy changes which occur as a result
of internal mixing processes in deep, stratified
impoundments;
2. temperature changes which may result from the in-
stallation and operation of multiple outlets or other
physical devices for control of the temperature
regime; and
3. thermal energy distributions which result in weakly-
stratified reservoirs where so-called "thermal
underflows" may occur, i.e., where longitudinal
temperature gradients are significant.
These required capabilities were to be supplied to FWPCA in the form
of a specific model, or models, compatible with the FWPCA river-
run model and to be available for use at the end of the contract
period. Further details can be found in the contract proposal (14).
- 7 -
-------
ORGANIZATION OF THE RESEARCH EFFORT
GENERAL APPROACH
At the time that research work was initiated on this con-
tract, WRE had in hand a deep reservoir model which had been recently
developed under contract with the Department of Fish and Game of the
State of California (15). The model had been successfully verified
for deep water bodies having low and moderate discharge to volume
ratios by applications to Lake Washington in the State of Washington
(16) and to Fontana Reservoir in the TVA system, respectively. The
reader who is not familiar with this basic model is encouraged to
review reference 15 before proceeding, since the work reported herein
is specifically oriented toward the further development of this model.
A review of the initial applications of the model on Lake
Washington and Fontana Reservoir established that further refinement
of two features of the deep reservoir model would greatly enhance its
general applicability. These features were, specifically:
a. the description of the internal miixing process
within a reservoir; and
b. the methodology for introducing inflows and ex-
tracting outflows from the reservoir model.
The accomplishment of objectives 1 anvd 2 required that
the above refinements of the deep reservoir model be made in order
to provide FWPCA with a model displaying the capabilities required for
its general application. Since the development, of these capabilities
- 8 -
-------
can best be accomplished by working with real data, it was appro-
priate that a prototype test reservoir be selected which displayed
strong stratification characteristics and for which there existed
a data base adequate for the analyses contemplated. Hungry Horse
Reservoir in Montana was chosen as the prototype best suited to
the research needs. A set of test simulation runs was conducted on
this reservoir for the purpose of refining the description of the
internal mixing process and incorporating into the model the theory
of selective withdrawal from stratified impoundments. The base
simulations also provided opportunity for the identification and
correction of any weaknesses in the model. In addition, the reser-
voir itself served as a test case for verification of the model on
a deep reservoir typical of those found in the Pacific Northwest.
The accomplishment of objective 3 required that the gener-
alized model developed in meeting objectives 1 and 2 be restructured
to conform to the constraints imposed by a reservoir displaying the
characteristics of a weakly-stratified impoundment. This restructured
model was then tested on Lake Roosevelt which is a prototype of
the weakly-stratified reservoir class.
REPORT FORMAT
Since the successful accomplishment of objective 3 required
the meeting of objectives 1 and 2, it was expedient to break the work
under this contract into two phases or parts: Part I dealing with
the research and development pertaining to the base simulations of
the model on a deep reservoir; and Part II addressing itself to the
- 9 -
-------
restructuring of the generalized model developed in Part I of the
research, and the results of applying the new model to a weakly-
stratified reservoir. Consequently, the report contained herein is
also broken into two parts which are titled: Part I — BASE SIMULATION
STUDIES; and Part II--SIMULATION OF WEAKLY STRATIFIED RESERVOIRS.
Chapters III through VII constitute Part I of the report.
The first four chapters describe, sequentially, the steps taken in
meeting objectives 1 and 2. These steps were as follows:
I. An initial set of base simulations were conducted
on Hungry Horse Reservoir in order to: a) identify
(and correct) any specific weaknesses which might
exist in the model; b) establish a better criteria
for introducing tributary inflows into the model;
and c) explore the effect of varying the functional
representation of the internal mixing process.
2. On the basis of the knowledge gained from the base
simulations, efforts were concentrated on refining
the description of the internal mixing process.
3. The theory of selective withdrawal from stratified
impoundments was reviewed and incorporated into the
model.
4. Two test cases (one practical, the other hypothetical)
were simulated for downstream temperature control.
The summary and recommendations for Part I are contained
in Chapter VII.
- 10 -
-------
Part II of the report, which constitutes Chapters VIII
through XII, is more application oriented than Part I since the
basic unit of the weakly-stratified reservoir model is the deep
reservoir model. Part II contains, sequentially:
I. A conceptual description of the weakly stratified
reservei r model.
2. The process of adapting the model to the geometric
and hydraulic characteristics of the reservoir.
3. The details of restructuring the deep reservoir
model.
4. Application of the weakly-stratified reservoir
model to Lake Roosevelt.
The summary and recommendations for Part II are contained in Chapter
XII.
In addition to the material contained herein, a Supplement
to this report has been prepared under separate cover. The Supple-
ment is devoted to a detailed explanation of the use, utility, and
details of the computer program employed for simulating the internal
temperature structure of thermally stratified impoundments.
AUTHORITY AND RESPONSIBILITY
The work reported herein was performed under Contract
No. 14-12-422, dated June 28, 1968, between the Federal Water Pol-
lution Control Administration and Water Resources Engineers, Inc.
All aspects of the technical work were conducted by the staff of
WRE under the direct supervision of Dr. Gerald T. Orlob, President.
- 11 -
-------
Dr. Larry A. Roesner coordinated the project and was responsible for
the theoretical development; Mr. W. R. Norton developed the computer
simulation package and contributed substantially to the practical im-
plementation of the theory.
Dr. R. W. Zeller and Mr. John Yearsley, both of whom were
members of the FWPCA Study Team, were closely associated with the
project throughout its execution and were responsible for supplying
all the physical, hydrologic, and meteorologic data used in the develop-
ment and testing of the reservoir temperature models.
ACKNOWLEDGMENTS
WRE has enjoyed the close liaison and cooperative attitude
that existed between itself and the FWPCA during the performance of
this work. The establishment of this rapport is directly attributable
to Dr. Bruce A. Tichenor, project officer for the FWPCA in this study,
Dr. Zeller, and Mr. Yearsley. The attitude and enthusiasm of these
gentlemen materially aided the smooth and successful completion of
this study.
- 12 -
-------
II. CLASSIFICATION OF RESERVOIRS
In Chapter I, three different types of reservoirs were
mentioned and it was inferred that a different approach was required
to mathematically model the thermal regime of each type. It is
the purpose of the following discussion to define these reservoir
types or classes somewhat more precisely and to lay down some
general guidelines for reservoir classification.
The hydrodynamic and thermal behavior of a fluid system
is completely described by the simultaneous solution of three
differential equations, namely:
I. the Navier-Stokes equations
(.equations of motion);
2. the conservation of mass equation; and
3. the conservation of energy equation.
Unfortunately, however, these equations cannot be solved except for
a very limited number of steady-state cases. Consequently, real
world descriptions of the behavior of prototype systems require that
types or classes of solutions be formed which describe, in an approxi-
mate way, the actual hydrodynamic and thermal processes occurring
within the system. As applied to the specific objectives of the work
reported herein, these approximate solutions must be formulated in
- 13 -
-------
a manner such that they account for the principal hydrodynamic and
thermodynamic processes affecting the internal temperature structure
of the reservoir.
Experience with prototype reservoirs indicates that there
are three distinct classes of reservoirs, each of which requires a
different type of solution for the determination of its internal energy
(temperature) distribution. These classes, which were alluded to in
the introduction, are:
I. The deep reservoir which is characterized by horizon-
tal isotherms;
2. the weakly-stratified reservoir, characterized by iso-
therms which are tilted along the longitudinal axis
of the reservoir; and
3. the completely mixed reservoir whose isotherms are
vert i ca I ..
By their very definition, it is apparent that the deep reser-
voir solution and the completely mixed reservoir solution are one-
dimensional, because their thermal properties vary in one direction
only. The difference, of course, is that the solution is taken along
the vertical axis for the deep reservoir while the completely mixed
reservoir solution is taken along the longitudinal axis. In contrast,
the solution of the weakly-stratified reservoir requires that both
dimensions be considered.
- 14 -
-------
The single most important parameter determining into which
class a reservoir will fall, and hence the type of solution required,
is the densimetric Froude number which can be written for the reservoir
as
IF = ft 1^ (1)
where
L = reservoir length;
Q = volumetric discharge through the reservoir;
D = mean reservoir depth;
V = reservoir volume;
p = reference density;
3 = average density gradient in the reservoir; and
g = gravitational constant.
This number is the ratio of the inertial force of the horizontal flow
to the gravitational forces within the stratified impoundment; con-
sequently, it is a measure of the success with which the horizontal
flow can alter the internal density (thermal) structure of the reser-
voir from that of its gravitational static-equilibrium state.
In deep reservoirs, the fact that the isotherms are horizontal
indicates that the inertia of the longitudinal flow is insufficient
to disturb the overall gravitational static equilibrium state of the
reservoir except possibly for local disturbances in the vicinity of
the reservoir outlets and at points of tributary inflow. Thus, F
would be expected to be small for such reservoirs. In completely mixed
reservoirs, on the other hand, the inertia of the flow and its attendant
turbulence is sufficient to completely upset the gravitational structure
- 15 -
-------
and destratify the reservoir. For reservoirs of this class, F would
be expected to be large. In between these two extreme classes lies
the weakly-stratified reservoir in which the longitudinal flow possesses
enough inertia to disrupt the reservoir isotherms from their gravita-
tional static-equilibrium state configuration, but not enough to com-
pletely mix the reservoir.
For the purpose of classifying reservoirs by their Froude
number, Band p in equation 1 may be approximated as 10 kg m and
o 0
1000 kg m , respectively. Substituting these values and g into
equation 1 leads to an expression for F as
F - 320 fc $ (2)
where L and D have units of meters, Q is in cms, and V has units of
From this equation it is observed that the principal reservoir
m3.
parameters that determine a reservoir's classification are its length,
depth, and discharge to volume ratio ( H- ).
In developing some familiarity with the magnitude of F for
each of the three reservoir classes, it is helpful to note that
theoretical and experimental work in stratified flow indicates that
flow separation occurs in a stratified fluid when the Froude number
is less than about - , i.e. for F < - , part of the fluid will be
IT TT p
in motion longitudinally while the remainder is essentially at rest.
Furthermore, as F becomes smaller and smaller, the flowing layer be-
comes more and more concentrated in the vertical direction. Thus,
in the deep reservoir it is to be expected that the longitudinal flow
is highly concentrated at values of F « -while in the completely
IT n
mixed case F must be at least greater than — since the entire
77 1
reservoir is in motion and it may be expected in general that F »— .
- 16 -
-------
Values of F for the weakly-stratified case would fall between these
two limits and might be expected to be on the order of — . As an
illustration, five reservoirs are listed below with their Froude
numbers. It is known that Hungry Horse Reservoir and Detroit Reser-
voir are of the deep reservoir class and can be effectively described
with a one-dimensional model along the vertical axis of the reservoir.
Lake Roosevelt, which has been observed to fall into the weakly
stratified class is seen to have a Froude number on the order of —
which is considerably larger than F for either Hungry Horse or De-
troit Reservoirs. Finally, Priest Rapids and Wells Dams, which are
essentially completely mixed along their vertical axes, show Froude
numbers much larger than — , as expected.
RESERVOIR FROUDE NUMBERS
AVERAGE
RESERVOIR LENGTH DEPTH
me te rs me te rs
Hungry Horse 4.7x10 70
Detroit 1.5xl04 56
Lake Roosevelt 2.0xl05 70
A
Priest Rapids* 2.9x10^ 18
Wells* 4.6xl04 26
DISCHARGE TO
VOLUME RATIO F CLASS
sec'1
1.2 x 10"8 0.0026 Deep
3.5 x 10"8 0.0030 Deep
5.0 x 10"7 0.46 Weakly-
Stratified
-fi
4.6 x 10 2.4 Completely
Mi xe d
6.7 x 10"6 3.8 Completely
Mixed
*River run dams on the Columbia River below Grand Coulee Dam
- 17 -
-------
The above illustrations tend to confirm the validity of
using the densimetric Froude number to classify reservoirs; however,
just what values should be used to distinguish the transition points,
between the individual reservoir classes, is not clear at this time.
The determination of these points requires further investigation.
- 18 -
-------
PART I -- BASE SIMULATION STUDIES
-------
III. INITIAL MODEL APPLICATIONS
HUNGRY HORSE RESERVOIR
All base simulation studies in Part I were performed
using Hungry Horse Reservoir in Montana as the prototype. This
reservoir is located on the South Fork of the Flathead River in
west central Montana and is a major storage reservoir in the Colum-
bia River power system. It has a rated, useable storage capacity
of 2.98 million acre-feet and was constructed by the U.S. Bureau
of Reclamation during the period 1948-1953. A plan view of the
reservoir is illustrated in Figure 2 and information pertinent to
energy budget calculations is give in Table 1.
Time and money constraints required that the period of
prototype simulation be limited to a single year. Calendar year
1965 was recommended by the FWPCA Study Team as the period during
which the data base was most accurate. Although there was only suf-
ficient data to simulate the period from 11 May to 20 October, this
time period was of sufficient length to allow simulation of the
reservoir from the onset of stratification in the spring, through the
summer heating period, and well into the fall cooling portion of the
annual cycle. Since these periods are the most critical that the
deep reservoir model must be able to simulate, the data were regarded
as adequate for the required research and development. Average
center!ine temperature profiles for the main body of the reservoir,
- 19 -
-------
South Fork
Flat head River
Aurora Cr.
LEGEND
GAGING STATIONS,REC DING
GAGING STATIONS.
W/TEMP. REC., RECORDING
MISCELLANEOUS
MEASUREMENT SITES *
THERMAL SURVEY POINT
RAFT STATION, RECORDING
PRECIPITATION STATION,
TEMPORARY
PRECIR, TEMP., WIND, HUM.,
EVAPORIZATION STATION
RADIATION STATION,
RECORDING
GROUND WATER OBSERVA-
TION WELLS, RECORDING
RESERVOIR GAGE IN DAM
* PERODIC ESTIMATES OR
MEASUREMENTS MADE
ON 40 UNREFERENCED
TRIBUTARY SITES
Flossy Cr.
Clayton Cr.
Go/die Cr.
Murray Cr.
Me Inernie Cr.
Deep Cr.
Wheeler
Cr.
TANA
0123456
SCALE
IN MILES
.,£• SPOTTED BEAR
|U RANGER STATION
PERMANENT STATION
FIGURE 2. Hungry Horse Reservoir
- 20 -
-------
TABLE 1
PHYSICAL AND METEOROLOGICAL INFORMATION
HUNGRY HORSE RESERVOIR
PHYSICAL DATA
Location:
Elevation:
Intake Elevations:
Spillway:
Mean Gaged Inflow:
Mean Ungaged Inflow:
Average Reservoir Discharge
to Volume Ratio:
Latitude 48°20', Longitude 114°OT
1085neters (normal high water surface
elevation)
Lower Discharge Tubes - 975 meters
Penstock Intakes - 1013 meters
Type - Gloryhole, elevation 1082 meters
73.2 CMS from 9 tributaries
39.3 CMS
0.82 yr
-1
SURFACE AREA,m2xlO~8
10
1100
1.0
2.0
CD
- 1050
o
< 1000
LJ
_J
^ 950
AREA-
•^-VOLUME
0 1.0 2.0 3.0 4.0 5.0
RESERVOIR VOLUME, m3 x IO"9
METEOROLOGICAL DATA
Solar Radiation:
Atmospheric Radiation:
Cloud Cover:
Barometric Pressure:
Wind Speed:
Air Temperature
(wet bulb and dry bulb):
Solar Radiation Extinction
Coefficient in Water:
Evaporation Coefficient:
computed
computed
data from Kalispell
data from Kalispel1
data from Hungry Horse Reservoir
data from Hungry Horse Reservoir
-1
supplied by FWPCA, 0.345 m"
supplied by FWPCA, 2.1 x 10~y mps/mb/mps
-9
- 21 -
-------
observed at various times during the simulation period, are illus-
trated in Figure 3. The ability of the model to reproduce these
profiles was taken as the criteria for evaluating its performance.
OBJECTIVES
In the initial model applications a series of base simu-
lations was conducted. The purpose of these simulations was two-fold:
first, to identify any specific weaknesses in the model; and second,
to explore the effect of varying the functional representation of the
internal mixing process. Since the results of investigations on the
internal mixing process fall more logically into Chapter IV, this
chapter is devoted primarily to the identification and correction of
specific weaknesses found in the model.
MINOR MODEL REFINEMENTS
The initial applications of the deep reservoir model re-
vealed that the computer code contained deficiencies in the following
areas:
I. Computational efficiency;
2. Adequate detail in the description of the move-
ment of the water surface; and
3. Flexibility in the specification of the simulation
time step.
- 22 -
-------
MOO-
1050 —
1000-
(O
k_
o>
950
11 MAY
T
-26 MAY r-9 JUNE
V
I
- 23 JUNE
V
J
357 3579 357911 13 15 357911 13 15 17
TEMPERATURE,°C
h-
> MOO-
_j
LJ
1050
1000-J
950
7 JULY
\
-21 JULY
11 i
I
3 5 7 9 11 13 15 17 19 3579
TEMPERATURE,°C
13 15 17 19 21
FIGURE 3. Observed Temperature Profiles, Hungry Horse Reservoir, 1965
(Continued next page)
- 23 -
-------
\\QO-r4 AUGUST
1050-
1000-
CO
950
r 18 AUGUST
3 5 7 9 11 13 15 17 19 21 23 3 5 7 9 11 13 15 17 19 21 23
TEMPERATURE,°C
< 1100 -T- / SEPTEMBER
>
UJ
V
Ijj
1050
1000-
950
-16 SEPT.
V
I
r5 OCT.
\
- 2O OCT.
35791! 13 15 17 19 357911 13 15 35791! 357911
TEMPERATURE,°C
FIGURE 3 (continued).
Observed Temperature Profiles, Hungry Horse
Reservoir, 1965
- 24 -
-------
COMPUTATIONAL EFFICIENCY
At the time this work was initiated, WRE had only recently
completed the basic development and verification of the computer
code for deep reservoir simulation. As is often typical of develop-
mental work, computational efficiency was not strongly emphasized.
Consequently, a critical review of the computer code revealed that
many internal computational refinements could be made to shorten the
solution time. It is estimated that the reprogramming effort improved
computational efficiency approximately fifty percent which resulted
in an overall decrease in computer execution time of about 25 percent.
MOVEMENT OF THE WATER SURFACE
In the original version of the model, each reservoir slice,
or elementjWas of equal height. The reservoir surface was taken as
the top of the uppermost reservoir element which was completely filled
with water; partially filled elements were ignored. As a result of
this condition, changes in the water surface elevation were handled in
a stepwise fashion with complete elements being added or subtracted
from the reservoir as their upper surfaces became submerged or exposed,
respectively. With the height of an element typically taken as one
meter, changes in water surface elevation were accounted for in steps
of approximately three feet.
The step-wise approximation of the water surface elevation
has been eliminated by allowing the height of the surface element to
be variable (between one and two times the standard element thickness),
the upper boundary of the element coincides with the water surface
elevation. This refinement is believed to produce a better description
of the temperature structure near the reservoir's surface.
- 25 -
-------
VARIATION OF THE SIMULATION TIME STEP
Another refinement was made which allows the use of a simu-
lation time step that is different from the interval between successive
observations of the meteorologic and hydrologic input data. In the
original version of the model, the simulation interval was restricted
to the same duration as the interval between observations of meteorologic
input data. In the updated version of the computer code the simulation
time step may be either longer or shorter then the interval between
data observations; all necessary bookkeeping and time averaging
necessary for this feature are handled internally by code. This par-
ticular refinement was incorporated in order to facilitate reservoir
simulation in instances when advective heat transfer violates the basic
assumptions of the model. This problem occurs when, for a given time
step, the mass volume of flow through a reservoir element is greater
than the volume of the element.
Experience indicates that the most critical period with
regard to advection violations is in the spring of the year when the
tributary inflows to the reservoir are at their maximum; this situation
may also be expected to develop during periods when reservoir release
rates are high. One of the easiest ways to rectify this violation is
to shorten the simulation time step, a manipulation that is now trivial
since the required data conversions are handled internally.
TRIBUTARY INFLOWS
As stated in the introduction, the method of introducing
inflows to the reservoir was known to need refinement. The original
method was to place the entire flow of a tributary into that reservoir
- 26 -
-------
element whose density was nearest to that of the tributary flow. The
main objection to this criteria was that the inflow was artificially
constrained between the upper and lower planes bounding the element,
while in reality the vertical distribution of the flow varies with
both the degree of stratification and the quantity of flow.
In order to account for the effects of stratification and
the magnitude of inflow, the following scheme has been adopted.
First, the elevation, z , is located where the density of the reser-
voir is equal to that of the inflow. Next, the vertical thickness
of the inflow, D, is computed from equation 27-B in Chapter X. The
inflow is then distributed between the elevations z + 4- D and
1 ° 2
z - p" D, assuming a uniform velocity distribution.
The base simulations did not reveal a need for additional
refinements. It was concluded from these studies that the successful
refinement of the description for the internal mixing process plus
the incorporation of selective withdrawal theory into the model would
produce a general deep reservoir model with the capabilities required
to satisfy objectives 1 and 2.
- 27 -
-------
IV. INTERNAL MIXING STUDIES
Of the total effort expended in Part I of the research,
the major portion of it was devoted to the study and refinement
of the description of the internal mixing process in deep reservoirs.
The expenditure of this amount of effort was anticipated at the out-
set of the work because internal mixing is a basic feature of the
deep reservoir model. Thus, its proper description is imperative to
the proper representation of the vertical distribution of heat within
a stratified impoundment.
In this chapter, a qualitative overview of the experience
and theory relating to internal mixing in reservoirs is presented
in the first section. The second section reports the results of
varying the functional representation of the mixing process, and
describes the form that was chosen as best suited to the description
of the process. The last section describes the verification of the
deep reservoir model on Hungry Horse Reservoir as simulated with
the newly developed functional representation of the internal mixing
process.
- 29 -
-------
MIXING PROPERTIES OF RESERVOIRS
The two major factors which govern mixing in reservoirs
are gravity and turbulence. Gravity may cause or enhance mixing or
it may suppress mixing depending upon whether the density strati-
fication in the reservoir is unstable or stable, respectively.
Turbulence, on the other hand, can only promote mixing. When gravity
alone is responsible for mixing, the mixing process is called
"natural convection," while mixing in the presence of turbulence,
with or without gravitational effects, is referred to as "forced
convection."
Ordinarily, forced convection is the phenomenon of primary
concern in the evaluation of the mixing properties of reservoirs,
although at certain times, when the surface of the reservoir is
cooling, natural convection plays a primary role in epilimnetic
mixing. Since experience with the deep reservoir model has revealed
no problems in the present method of describing natural convection,
this phenomenon will not be discussed here. Details of its role in
distributing heat along the vertical axis of the reservoir can be
found in reference 15.
Mixing implies the transfer of materials or properties
within a system from points of high concentration to points of low
concentration, and vice versa. For a system which is experiencing
forced convection, it is a rule of experience that the time rate of
transport, F, of a property, S, through the system is (other things
being equal) proportional to the rate of change of concentration of
this property with distance, z. In equation form this rule is
- 30 -
-------
expressed as
F - -A g , 0)
where A is the coefficient of proportionality. The mixing process
as defined by equation 3 is often called "effective diffusion,"
"eddy diffusion," or the "diffusion analogy" since it is identical
in form to the equation describing the process of molecular diffu-
sion. The difference between the two processes, however, is that
for molecular diffusion, A is constant, while for turbulent transfer,
A is a function --of the dynamic character, or the turbulence level,
of the system. In general, A is a temporal and spatial variable,
and thus will be referred to as A(z,t). Equation 3 rewritten for
heat flow along the reservoir vertical axis is
H = -pcA(z.t) |J , (4)
where
H = heat flux, HL'V1;
p = density of water, ML~ ;
c = heat capacity of water, HM~ D~ ;
2 -1
A(z,t) = coefficient of eddy conductivity, L T
T = temperature, D;
z = elevation in the reservoir, L; and
t = time, T.
The proper or correct choice of A(.z,t) is a difficult task
because it represents the prototype response to many unknown and/or
poorly defined random inputs. However, experience with the coeffi-
cient (see Figure 4) has disclosed some of its general properties:
- 31 -
-------
TEMPERATURE, C
oo
ro
tO 15 20 25 \7
oc
20
40
60
80
100
120
5 10 15 20 2;
f
\
\
°o\A
°\
0\
A/
/ '
/
/
1
|
,
*r
f
\
V
V
r °
HUNGh
^^
'Y HORS
J*
E RESE
'RVOIR
EFFECTIVE DIFFUSION COEFFICIENT, cm2/sec
FIGURE 4. Effective Diffusion and Temperature Profiles for
Selected Impoundments (from Reference 15)
-------
I. Effective mixing in the epilimnion is usually
greatest near the surface but declines rapidly
with depth as the thermocIine is approached.
2. The effective diffusion coefficient is minimal
at, or very near the position of the maximum
temperature gradient, the so-called "thermocIine."
3. In the hypolimnion the effective diffusion coef-
ficient increases with depth in a more or less
erratic manner, reaching a maximum at about mid-
depth, thereafter decreasing as the bottom is
approached.
4. The value of A(z,t) is vastly greater than the
coefficient of thermal diffusivity for heat trans-
-7 2
fer by molecular conduction (K - 1.4 x 10 m /sec).
While it is impossible to give a detailed theoretical de-
scription for the observed behavior of A(z,t), the following arguments
are proposed as a possible explanation. A horizontal current with
an average velocity, u, at elevation, z, has a kinetic energy per unit
2
mass of KE = u /2. If it is assumed that the horizontal flow in a
reservoir is basically a boundary layer flow (except in the region
of the outlets and at points of tributary inflow), the vertical dis-
tribution of KE would then appear as shown in Figure 5-A. This amount
of energy is expended through turbulence; however, in a stratified
reservoir it is not all available for turbulent mixing as it is a
non-stratified impoundment, for the following reason: In a stably
stratified reservoir the vertical differences in density cause a
buoyant restoring force to be exerted on any parcel of fluid which is
displaced from its position of static equilibrium. Thus, of the
- 33 -
-------
A. KINETIC ENERGY
DISTRIBUTION IN
BOUNDARY LAYER FLOW
B. POTENTIAL ENERGY
DISTRIBUTION
C.NET FLOW ENERGY
DISTRIBUTION
NE=KE-PE-s*-
-j-—BE*.
CO
I
D.DISTRIBUTlON OF
WIND ENERGY
E.NET MIXING ENERGY
DISTRIBUTION
—WE stratified
A
-------
total KE available, only that energy which is in excess of that
required to overcome the buoyancy forces is available for mixing.
For a reservoir with a temperature profile of the form
shown in Figure 5-B, the vertical distribution of the work per unit
volume that is required to overcome the density stratification is
given by the curve PE which is shown in the same Figure. Thus, the
net amount of turbulent energy NE which is available for mixing is
the difference, KE - PE, with the restriction, NE >_ 0. A typical curve
of NE for the temperature profile given in Figure 5-B is shown in
Figure 5-C.
In addition to turbulence induced by the boundary layer
flow, mixing in the upper layers of the reservoir is contributed
to by the wind. Following arguments of Eckman as related by Hutchinson
(3), the energy input per unit volume, WE, should experience an expo-
nential decay with depth as shown in Figure 5-D. In the case of
a stratified reservoir, the same arguments apply for reduction of
the diffusion coefficient due to stratification as applied to the
boundary layer flow. As such, it is generally thought that the ef-
fects of wind mixing do not persist below the thermocline.
Figure 5-E, which is the vertical distribution of the
energy available for mixing, ME, follows from the linear combination
of NE for the boundary layer flow and WE resulting from wind mixing.
It is observed that the combined effect produces a form of ME which
is similar in its shape to the prototype forms of A(z,t) illustrated
in Figure 4. This is an expected result since the value of A(z,t),
at any elevation, z, is a direct measure of the turbulent energy
available for mixing at that elevation.
- 35 -
-------
While the preceding discussion has been devoted to a qualita-
tive description of the functional form of the eddy conductivity coef-
ficient, the pertinent energy processes can be expressed mathematically.
For the sake of completeness, these relationships are given below.
Details can be found in reference 17.
If the energies illustrated in Figure 5 are now thought of as
the rate of energy expenditure per unit time and per unit mass, then
KE becomes the equivalent rate per unit mass at which kinetic energy is
transferred from mean motion to turbulent motion and is expressed as
KE = AT(|f)2 (5)
where A is the eddy viscosity coefficient and u is the mean horizontal
velocity at elevation, z. The subscript,. T, has been added to distinguish
this coefficient from the eddy conductivity coefficient, A(z,t). Simi-
larly, defining PE as the average work done by turbulent motion against
gravity, or buoyancy, per unit mass and unit time leads to
PE - -9A(2.t)l|f f {6)
where g is the gravitational acceleration. Combining equations 5 and 6
gives the value of NE, the energy available per unit mass and per unit
time for the vertical transport of heat^ i.e. mixing;
NE = [AT(§Y)2 + gA(z,t) £f£ J ,
subject to NE > 0.
- 36 -
-------
If the effect of wind is neglected, the following observation
can be made. Since KE is the rate at which turbulent energy is supplied
to the system and PE is the rate at which turbulent energy is used to
overcome the gravitational stability of the system, then, in order to
maintain turbulence in the system, KE must be larger than PE. From
equations 5 and 6, it follows that
_q13fi A
9 p 3z_ < AT
,9jV A(z ,t) (8)
The dimensionless quantity on the left-hand side of equation 8 is called
the "gradient Richardson number," R. (18). The function, -^2- , is
I P oZ
commonly referred to as the "gravitational stability," E, of a stratified
water body. If this relationship is substituted into equation 8, there
results for the description of the Richardson number,
(3U/3Z)2
The efficacy of using the Richardson number as a measure of the
energy available for vertical mixing has been the subject of much dis-
cussion. The problem revolves around the ratio, A /A(z,t), which is, at
AT
best, difficult to define. If R. < -^7—ry , then mixing will occur;
1 f\(Z ,Z)
AT
conversely, if R- > A/7 -FT > turbulent mixing will not occur. Without
- 37 -
-------
going into further detail, suffice it to say that a small Richardson
number indicates maintenance or increase of turbulence, whereas a large
value indicates suppression or extinction of turbulence. In addition
to the problem of defining the critical value of R. (where R. = ,( T.^ ),
r I I *» \ ^ j ^ y
the lack of a theory describing the details of the reservoir's internal
flow structure precludes an accurate evaluation of the velocity gradient.
Hence, the evaluation of the Richardson number, itself, is a formidable
task.
Under the assumption that energy input by the wind is an
exponential function of the depth, WE is described simply as
WE = WEQ e"ad , (10)
where WE = energy input per unit volume per unit time at the
water surface;
a = decay constant; and
d = depth below the water surface.
The value of a is chosen so that WE is very small (say, 0.01 WE ) at
the depth of the thermocline. WE is a function of the wind stress at
the water surface. Some direction in the assessment of its magnitude
can be found in Hutchinson (3).
Finally, through the combination of equation 7 and 10, it is
possible to write an expression for the amount of turbulent energy in
the reservoir that is available for mixing. This expression is
- 38 -
-------
ME . HE0 e-" + [AT(|f)2 + gA(z,t) 1 ff 3 • <">
The examination of equation 11 and the prototype behavior of A(z,t) as
illustrated in Figure 4 reveals that while the behavior of a thermally
stratified reservoir seems to follow certain laws, it is apparent that
the overall complexity of real impoundments precludes an explicit, detailed
representation in mathematical form. Thus, certain simplifications must
be made. The simplifications made in defining the temporal and spacial
variation of the eddy conductivity coefficient, A(z,t), are described
in the following section.
EFFECT OF VARYING THE FUNCTIONAL REPRESENTATION OF THE INTERNAL MIXING
PROCESS
The functional representation of the vertical transport of heat
as a result of the internal mixing process was defined in the last section
by equation 4. It is apparent from this equation that any variation in
its functional form can be accomplished only by varying the functional
description of the eddy conductivity coefficient, A(z,t). Consequently,
this section is devoted to describing the results of using different
functional representations for A(z,t).
The functional representation of A(z,t) must, of necessity, be
empirical in nature since the eddy conductivity coefficient describes
the prototype response to many unknown and/or poorly defined random
influences. However, on the basis of the theory and prototype experience
described above, it is observed that the choice of a functional form to
represent the eddy conductivity coefficient must be based on the ability
- 39 -
-------
of the function to display the following characteristics:
I. A(z,t) is large in the epilimnion;
2. A(z,t) is small at the thermocline, or preferably,
A(z,t) should decrease as the stability, E, increases
. . . _ -I 9p . ,
(recalI E = — ~- ); and
P oZ
3. A(z,t) is large in the hypolimnion, being of the same
order of magnitude as its value in the epilimnion.
At this point, it is worth mentioning that Hungry Horse
Reservoir was an ideal prototype for studying internal mixing. Because
of its small discharge-volume ratio, only a small portion of the vertical
heat transfer occurred by advection. Consequently, "eddy diffusion"
played a major role in the vertical distribution of heat. This fact
can be readily observed in Figure 6 where the solar, vertical advection,
and diffusion fluxes within the reservoir are graphed as a function of
elevation for a typical day.
EXPONENTIAL REPRESENTATION
The first functional representation of A(z,t) to be tested was
the original function proposed in the Fish and Game Report (reference 15)
Summarily, this representation was:
A(z.t) = A0(t)
ACz,t) = A[zT,t) , ZT> z i
where AQ(t) - value of the eddy conductivity coefficient at the
water surface;
- 40 -
-------
-2
VERTICAL HEAT FLUX, k cal m sec.
TEMPERATURE, C
1090
1070
1050
1030
1010
u>
o>
"^ 990
o
— 970
,64
io2
10 3
LJ
950
1
1
L
ADVECTION"
FLUX ZERO
BELOW ELEV
9fi7
^
^
SO
\
^^
^
| AL
1
V
\
\
/
LAF
Vh
/
? FLUX .
:TION FLU
/
•q;-ADVECTW
X
/
/
VLJJ.
1UI IUIN rLUX
/
}
1UP\
DOV
/
/
/
/l
)IFF
VIM 1
/
/
/
'USION FL
UX
t
^
= 33l.3lcms 1 TRIBUTARY
= l40.74cms V
-0 = ll.4lcms [INFLOWS
-0 = 16 60cms
= 3.37cms
LOWER OUTLET
WITHDRAWAL ZONE
= 4.22cms
FIGURE 6. Contribution to Vertical Heat Flux by Individual Heat Transfer
Mechanisms, Hungry Horse Reservoir, 1965
-------
n = decay coefficient;
2T = elevation of the thermocline; and
z elevation of the water surface.
Examination of equation 12 reveals that the function as
defined in equation 12-A displays characteristic 1 and characteristic 2
above the thermocline. Equation 12-B, however, tends to be in violation
of characteristic 2 below the thermocline and does, in fact, violate
characteristic 3. Although the weakness of equation 12-B was realized
in the original derivation of the function, it was not thought to be a
significant weakness for the following reasons. In the hypoll runion,
the thermal gradients are very small; thus, according to equation 4,
the vertical rate of turbulent heat transport will also be small.
Furthermore, this rate would be significantly small, in comparison to
the rate of heat transfer by other mechanisms, so that the error
resulting from the use of equation 12-B would not be a significant
percentage of the total rate of heat transport.
For testing the functional representation of A(z,t) as hypo-
thesized in equation 12, Hungry Horse Reservoir proved to be an ideal
prototype. Examination of the observed temperature profiles for the
reservoir, as illustrated in Figure 3, reveals that for most of the
simulation period the thermocline (defined as the elevation at which
92T
—2" - 0 ) is either at, or very near, the water surface. Thus, the
9z
model was required to simulate a region,in which the functional repre-
sentation of A(z,t) was known to be v/eak, i.e., through most of the
simulation period, A(z,t) was described by equation 12-B.
- 42 -
-------
As originally proposed, the determination of A (t) in
equation 12-A would be made on the basis of experience and/or by de-
ducing its value from the prototype data. In addition, the value of
the eddy conductivity coefficient at the thermacline; A(zT>t), had to
be obtained in the same manner. Once these two values were established,
n could be determined by the ratio, A(zT,t)/A (t), and the elevation
z-p of the thermocline. In the initial simulations, however, A and
A(ZT) were assumed constant, i.e., time invariant, and various values
and combinations of values were evaluated for their ability to repro-
duce the prototype behavior.
The resuTts of the test simulations demonstrated that the weak-
ness involved in using equation 12-B was a drawback to its general use
in reservoir simulation. In particular, it was found that an accurate
representation of the temperature profile in the clinolimnion (the
region of the reservoir between the thermocline and the "knee" of the
temperature curve in the hypolimnion) was not possible with this
function. Without going into details concerning the individual test runs,
the problems involved with such a functional representation for A(z,t)
below the thermocline can be described as follows. First, if A(zy) is
chosen as the prototype value of A(z,t) at the thermocline, its value
is very small. The use of such a small value for the entire reservoir
below the thermocline over-suppresses the turbulent rate of heat
transport down through the reservoir; i.e. it "traps" the heat in the
surface layers. This results in very high computed surface tempera-
tures and temperatures that are too cold in the reservoir depths. A
typical case is illustrated in Figure 7 where the value of A(zT,t)
C O 1 '
was taken as 5.0 x 10" m sec" . Increasing the value A(z^) to a
value such that the computed surface temperatures were about the same
as the observed values allowed too much heat to pass through the
region of steep gradients but still did not allow enough heat to pass
through to the lower reservoir depths. A typical profile computed
by using this approach is shown in Figure 8.
- 43 -
-------
Hungry Horse Reservoir
^SIMUL A TION DA Y, 9 JUNE 1965
1050
CO
Ld
_J
LJ
~ 1000
20
TEMPERATURE, C
FIGURE 7. Typical Computed Temperature Profile with A(z,t) Described
by Equation 12 and A(z.pt) - 5.0 x 10"5 ~2
m sec
1100
1050
1000
LU
_)
LU
950
jjungry Morse Reservoir
SIMULATION DAY, 9 JUNE 1965
0
COMPUTED
OBSERVED
J_J_J.J_J_J
12 16 20
FIGURE 8.
TEMPERATURE, °C
Tyoical Computed Temperature Profile with A(z,t) Described
by Equation 12 and A(zT,t) = 1.0 x ICf4 m2 sec"1
- 44 -
-------
The conclusions drawn from these test simulations were:
1) the amount of heat transferred by deep reservoir turbulence can
be a significant portion of the total heat transfer; and 2) in the
region of steep thermal gradients below the thermocline, i.e. the
clinolimnion, A(z,t) must have a value substantially less than that
value used for the hypolimnion.
In reflecting upon the functional form of A(z,t) as proposed
by equation 12, it is observed that the function is primarily oriented
to the description of turbulence induced by wind action. This conclusion
becomes obvious when equation 12-A is compared with equation 10 and
Figure 5-D. A review of the circumstances under which equation 12 was
formulated confirms that such a formulation was not accidental. The
prototype reservoir upon which the original function was developed was
Fontana Reservoir in the TVA System. In this reservoir, a primary
factor in the development and persistence of the epilimnetic region is
wind mixing. Moreover, the reservoir itself has a moderate discharge
to volume ratio (2.3 yr" ). Thus, while wind mixing is one of the
primary movers of heat in the epilimnion, the high discharge to volume
ratio for the reservoir indicates that the prime mover of heat in the
hypolimnion is vertical advection. For such a case, the reasoning which
underlies the formulation of equation 12-B is sound.
In contrast to Fontana, Hungry Horse Reservoir is almost com-
pletely opposite in its characteristics. First of all, wind mixing (at
least during the simulation period in 1965) is insignificant. Examination
of the meteorological inputs for the simulation period reveals that the
development of the epilimnetic region resulted from convective mixing as
a result of surface cooling, not from wind mixing. With regard to
hydrologic and geometric characteristics, Hungry Horse is of the low
discharge-volume type having an average ratio of about 0.8 yr . With
such a small ratio it is not surprising that turbulent heat transfer in
- 45 -
-------
the hypolimnion should represent a significant portion of the total
heat transport.
In summary, the initial test simulations revealed that
equation 12-B was inadequate for defining the eddy conductivity coe-
efficient below the thermocline. The adequacy of equation 12-A in
describing A(z,t) above the thermocline could not be evaluated since
the epilimnion was formed by natural convective mixing rather than by
wind action.
STEP FUNCTION
On the basis of the experience derived from testing the
exponential representation of A(z,t), it was decided to investigate
the use of a step function for its adequacy in describing the vertical
variation of the eddy conductivity coefficient. The advantage of such
an approach was that the capability would exist to assign one value of
A(z,t) to the epilininetic region, another value to the region of the
thermocline, i.e. the region of steep thermal gradients, and a third
value to the hypolimnetic region. Such a function would display all
the characteristics required of the eddy conductivity function. The
generalized mathematical formulation for this type of functional
representation was proposed as follows:
A(z,t) = A£(t) z > z£
A(z,t) - AT(t) z£ > z > ZH (13
A(z,t) = AH(t) ZH > z (13
Schematic diagrams of this functional form are given in Figure 9 for the
case of the thermocline at the surface and for the case of the thermo-
cline below the surface.
- 46
-------
Key.: PROTOTYPE BEHAVIOR/——PROPOSED REPRESENTATION
A.THERMOCL1NE AT SURFACE B. THERMOCLINE BELOW SURFACE
FIGURE 9. Representation of A(z,t) with a Step Function
The functional form of the eddy conductivity function as
proposed in equation 13 essentially divides the reservoir along its
vertical axis into three zones with A(z,t) exhibiting certain distinct
characteristics in each zone. Re-examination of the above section on
MIXING PROPERTIES OF RESERVOIRS reveals these characteristics to be as
follows:
- 47 -
-------
ZONE DEFINITION OF CHARACTERISTIC FEATURES OF A(z,t)
ZONE
(1) z < ZH An(t) is primarily a function of the
intensity of the deep reservoir turbu-
lence; it is only slightly dependent
upon the gravitational stability of the
reservoir.
(2) z,, < z < ZE AT(t) is primarily a function of the
gravitational stability.
(3) Zp < z AF(t) is primarily dependent upon the
amount of wind mixing.
Although in general, Ap Ay and A,, are all time variant functions, they
were fixed at constant values in this set of test simulations since the
interest was more in how this type of functional representation would
shape the temperature profiles (or the vertical heat distribution) rather
than how precisely individual observed thermal profiles could be
matched.
One of the questions that arose in connection with the func-
tional representation as given by equation 13 was how to choose the
value of Zr and z,,. Since A(z,t) is functionally dependent upon the
stability of the system, especially in Zone 2, it was decided that zu
rl
and z^ should also be based on stability criteria. Through a series
of trial simulations, it was found that the best results were obtained
when these elevations were taken to be the elevations in the reservoir
where the stability was 4 x 10" meters .
- 48 -
-------
The functional representation of A(z,t) as a step function
was found to substantially improve the modeling capability over that
which was possible through the use of equation 12. However, since
there was an order of magnitude difference between the value of AT
and those of A,, and Ar-, a problem developed at the switching elevations,
Zj. and z^. The problem was specifically this: for any given
temperature gradient at the elevation zr or z.. the substitution of
t n
equation 13 into equation 4 indicates that the rate of heat transfer
on either side of these two elevations will differ by an order of
magnitude, since AT is an order of magnitude smaller than A,, or A£.
Consequently, immediately above zp there will be a "stacking" of
heat which creates large thermal gradients at that point; at the ele-
vation z,,, on the other hand, the large value of Au causes a "draining"
rl n
of the heat at that elevation, thus creating very low temperature
gradients. The result of this situation is to cause a series of
"bumps" in the computed temperature profile. While this is not
necessarily harmful to the performance of the model, it is undesirable
from an operational point of view. An example of the type of profile
obtained under this type of functional representation for the eddy
conductivity coefficient is illustrated in Figure 10. The diffusion
coefficients for this case were 1 x 10" m /sec for Ac and Au and
fi ? t H
5 x 10"° m /sec for AT-
- 49 -
-------
Hungry Horse Reservoir
-4 2 -
A£ = 1.0x10 m sec
AT= 5.0x10 m sec
A = 1.0 x 10 m sec
n
8 12 16 20
TEMPERATURE, C
FIGURE 10. Typical Computed Temperature Profile with A(z,t)
Described by Equation 13; Simulation Day is
1 September 1965
REPRESENTATION OF A(z.,t) BY A CONTINUOUS FUNCTION
The results of using the step function to represent the
spatial variation of the eddy conductivity coefficient were encouraging.
It appeared that, if the abru-pt step changes could be smoothed out of
the describing function, a good representation of the spatial behavior
of A(z,t) would be obtained.
Reviewing the mechanics of turbulent heat transfer in a
stratified water body reveals that a precise definition of the spatial
variation of A(z,t) is most important in Zone 2 because the temperature
gradients are steepest and the spatial variation of the eddy conductivity
- 50 -
-------
coefficient is greatest in this region. Thus, it was desirable to
represent A(z,t) in Zone 2 with a continuous function, A~-(z,t),
that would have a minimum value at the thermocline and would be
equal to Ar(t) and A,,(t) at the elevations zf and z,,, respectively.
Since A-j-(z,t) is primarily a function of the stability, it follows
that such a relationship should be sought in the prototype data.
The first step in looking for a continuous function of
Aj(z,t) was to use the observed prototype temperature profiles and
simulate the reservoir "backwards" to obtain the eddy conductivity
coefficients required to produce these profiles. The details of
this method of solution can be found in reference 15. Through the
employment of this methodology, the average spatial variation of A(z,t)
was found for the following time periods:
May 11 - May 26
May 26 - June 9
June 9 - June 23
June 23 - July 7
July 7 - July 21
July 21 - August 4
August 4 - August 18
August 18 - September 1
September 1 - September 16
September 16 - October 5
Linear plots of the average spatial variation of A(z,t) during these
time periods are shown in Figure 11. Note the scale change in the
A(z,t) axis at 10.
- 51 -
-------
1100
11 May to 26 May 26 May 9 June 23 June 7 July to 21 July
"-"• ~T"'to ~~Yto ~1~*~ "^
9 June 23June
1050-
1000-
950'
'to
7 July
02468 '°305070
02468 024 0246 02468
o
I-
>
UJ
_J
LJ
IIOO-,
1050 —
1000'
950-
21 July
~ to
4 Aug.
4 Aug. __!8Aug. to I Sept. ___ I Sept. to 16Sept
'to
I8 Aug.
\ 1 I
02468
024 02468 10
30
0 2 4 6 8 10
16 Sept.
'to
5 Oct.
~T—1
024
EDDY CONDUCTIVITY COEFFICIENT, m2 sec"1 x 10*
FIGURE 11. Average Eddy Conductivity Profiles for Hungry Horse Reservoir
1965
- 52 -
-------
It is observed that while these forms are similar to those shown in
Figure 4 below the thermocline, there is no significant increase in
A(z,t) above the thermocline. The reasons for this behavior have been
explained previously.
The search for a functional relationship between AT(z,t)
and the stability, E, culminated in plot of In A(z,t) versus In E.
These plots of the Hungry Horse data, for the time periods given
above3 are illustrated in Figure 12. It is quite evident from Figure 12
that A(z,t) is primarily a function of the stability when E is greater
than about 10" m" (for a stability of E = 10" m~ , the corresponding
temperature gradient is approximately 0.06 °C m~ ). Moreover, the
relationship in this region is of the form
A(z,t) - bE"a (14)
If the average functional variation indicated by the dotted line in
Figure 12 is used, the value associated with a_ is 0.7 and b_ is
1.48 x 10"8 m1'3 sec"1.
The use of equation 14 as a definition for Ay(z,t) in Zone 2
strengthens the deep reservoir model in several ways. First of all,
in addition to providing a continuous spatial description of the eddy
conductivity coefficient, it also yields a temporal description of
the coefficient in this zone since Aj(z,t) is a function of the
temporal and spatial variation of the stability, E. Second, by the
use of equation 14 the step changes in A(z,t) at the elevations ZH
and zr are eliminated since zu is now defined as the elevation at which
t n
- 53 -
-------
Legend: on MAY- 26 MAY 1955 o 9 JUNE - 23 JUNE 1965 A? JULY-SI JULY 1965
A26MAY-9JUNE 1965 Q 23 JUNE - 7 JULY I9S5
IU.U
8.0
6.0
4.0
2.0
1.0
.8
.6
.4
.2
.1
.08
.06
.04
.02
ni
\
i
('
•o
-A
Q
n
— "~~o
f
A A
-n-
- —
j
^
^^•BK
.£,
f
^
^M
s
AAAAA
,
D
n
..— — _
r-
*>
A
^-«
a
o
~"
A
a r
o
>
_ a
A
1
!
O
!• ^!W
^•V
A
3 L
a^ r.
^
Q
o co
X^
&
A
\
v^
\
c
x^
X
\
©
A ^"'zr'-^'v
^ o i J, /-^V'f ^
•~! '3 j-i r. "\" _.
0°°°
=..
^ -J^;>j LJ
"M D
1 JK)
o 0
— , —
J V
^
A
>
\
\
0
LJ
A
A
s v
^
0
t
Upper
-Envelope
>*/_/_
N
^
\
a
r
,
X
A
V
\
\
A
?^
, x
U(
0.7
'N|
"s
Lower
Enve!opex
S.
Js
^"
\
d
Ox
-v
\
N.
^
\
\
sn \
\ E
X
A
A
\
k
NS\
X\
S
O
x
~0
o
CVJ
u.
LL
LJ
O
a
\-
o
Z)
Q
12.
O
O
>
Q
Q
UJ
.01 .02 .04 .06.03 .! .2
STABILITY E(z,t),.m"1 x!06
.4 .6 .8 1.0 2.0 4.0 6.0 8.0 10 20 40 60 80100 200 400
FIGURE 12. Loq of Eddy Conductivity versus Log Stability—Hungry Horse Data
(Continued next page)
-------
Legend: 021 JULY-4AUG. isss ,o ISAUG. - ISEPT. 1965 A is SEPT.-5OCT. 1955
^4 AUG.- 13 AUG. 1965 D I SEPT.-16 SEPT. 1965
O
•x
~O
QJ
to
c\J
f-
o
ID
Q
12.
a
o
Q
O
LL!
IU.U
8.0
6.0
4.0
2.0
1.0
.8
•6
.4
.2
.1
.08
.OS
.04
.02
0 "
— — .
cf
—5 —
- — . —
A
^'3
D
f
Q&Q
» M^na
A
yO(
—
?c
«v
' mtiK'
n
0
,— _. —
51.
°°'
.
1^
LOO
1
- —
c
fl» ^K
c
1
••V 1
H
l\
X
„ _
\
\
\
'?D:
^^
N
\
'T3
a n
h
O^^v r*^1' "V
,
°0
".
\
a
^
SN
^
A
O
* V
Upper
/Enveiope
^7
N
b
\
a
n ch
X
X
A
A
'\
rv
\
\
n
X
0
c
.(97A
XN! 0
X
Lower
Enve!opex
\
ncn
X
0
o
A'
N
N
•N» c
%
U.
\
^
s
i
)\
\
X
o%0>
o
o
v
X
*
\
0s
x\
X
-o
\
.0! .02 .04 .06.03 .! .2
STABILITY E(z,t)srp:' x!05
FIGURE 12 (cont'd).
.4 .6 .3 !.0 2.0 4.0 6.0 8.0 !0
20
40 60 80100 200 400
Log of Eddy Conductivity versus Log Stabilitv--Hunqry
Horse Data
-------
¥^
and zr is the elevation at which EF =
AE(ty
Third, with the more precise definition of AT(z,t) the modeling
capabilities in the region of steep thermal gradients are significantly
improved. An examp-le of the application of this refined approach to
the definition of A(z,t) is illustrated in Figure 13. Note, again, the
scale change in the A(z,t) axis at 10.
EDDY CONDUCTIVITY COEFF.
A, m2sec' xlO4
1100
to
CD
"S I05°
E
0
h-
> 1000
LJ
_J
LU
950
)
Tr
5 10
~r~T""i i~" " i r~"| n i i i
1 f^
\ jf
\ V
— \_^ f
NO
^
C ,
««
— c
-
-
— -
— -
A(z,t) //T(z,t)
.
A 1
\ |
V 1
* v 1
1
^J.
SIMULATION
20
i i n
JH •-
DAY
1 SEPT. 1965
observed —
simulated —
*•-
^- — ""^
r"_T_T_.j_.) ...i . 1 i 1 i
0 4 8 12
TEMPERATURE,
.
dcaneiran
U.^
-------
Having formulated the description of AT(z,t) in Zone 2, there
remained the problem of defining Ac(t) in Zone 3 and Au(t) in Zone 1.
t n
Since wind mixing was not observed in the prototype data, no further
development of Ar(t) was possible; thus, the effort was directed toward
a temporal description of A,,(t).
The value of An(t) reflects the level of deep reservoir turbulence
or more precisely, the rate of dissipation of energy by turbulence.
According to Hutchinson (3), the hypolimnetic turbulence in lakes is due
almost entirely to the internal seiche current system. While this is
probably true during the summer months, turbulence is also supplied to
the hypolimnion in the fall and spring by the deceleration of the
tributary inflows which enter the hypolimnion. In addition, in operating
reservoirs, the energy contained in the transient momentum flows that
are established through operation of the outlets must be dissipated
through turbulence in the hypolimnion.
In searching for some relationship between A,,(t) and the above-
mentioned parameters, the following scheme was tried:
where (AQ ,) = average of the absolute values of the daily
UUL aVG
change in discharge through outlets below the
thermocline; and
(Q. ) , = average inflow minus the average outflow below
the thermocline; positive values or zero.
The time over which the flows were averaged was 14 to 15 days which
corresoonds to the time interval between successive measurements of the
- 57 -
-------
thermal profile in the reservoir. Figure 14-A is a linear plot of
observed values of AH(t) versus (AQout)net + (Qjn)ner Since more energy
is dissipated through deceleration of flows than through acceleration,
it was felt that a better fit might be obtained by the equation:
= fl
where (-AQ . ). = average daily decrease in discharge through
out ave • ~~~ ™
outlets below the thermocline, i.e.
(-AQoutW ^ all negative values of AQQUt ^
period
Equation 16 is illustrated in Figure 14-B. Comparison of the two plots
in Figure 14 indicates that equation 15 describes best the temporal
variation of A,,(t). From Figure 14-A the value of An(t) may be expressed
as
AH(t) - 0.0275? , (u)
where ? = [(AQout)ave + (Q1n)netl , m3 sec-l. and yt) has Un1ts Qf m2 sec-l
Equation 17 gives a general description of A,,(t) in Hungry Horse
Reservoir. The applicicability of this equation to deep reservoirs in
general, however, requires further testing. Although the functional
relationship appears to be sound, it is almost certain that the constant
0.0275 will not be the same from reservoir to reservoir. Its value
will most likely be primarily a function of the volume of the reservoir,
because for a given value of c, the turbulent energy dissipated per unit
- 58 -
-------
o
o
o
CO
CJ
O "~"
0-ls
14 16 18 20
],cms
A. RELATIONSHIP BY EQUATION 15
10-i
o
x
To
V)
CJ
AH(t)= 0.600?,
B. RELATIONSHIP
EQUATION 16
FIGURE 14. Relationship between Hypolimnetic Value of Eddy
Conductivity Coefficient and Deep Reservoir
Hydraulics
- 59 -
-------
volume (upon which A(z,t) is based) will be directly proportional to the
reservoir volume.
PROPOSED FUNCTIONAL FOM FOR THE EDDY CONDUCTIVITY COEFFICIENT
Briefly summarizing the results of the research effort reported
in this section, the most important finding was, by far, the log-log
relationship between A(z,t) and the stability E, as expressed by
equation 14, in the region of steep thermal gradients. Although this
formulation is based on the observed behavior of a single reservoir,
some credibility is afforded to its general applicability by the fact
that this same relationship was found to hold for Detroit Reservoir in
Oregon, which was analyzed under a separate contract with the U. S.
Army Corps of Engineers. Regarding the functional form of A,,(t) as given
in equation 17, the conclusions are not firm; thus the equation cannot
be regarded as a general function without further substantiation. The
evaluation of AE(t)5 of course, could not be accomplished due to the
lack of observed wind mixing.
On the basis of the existing theory and the results reported in
this section, the functional form of the eddy conductivity coefficient
is proposed as follows in equation 18. Some discussion of the proposed
form is given below.
/ \
Afztl-/^0'^""2' 7 < 7 Mo/n
r\\ t. 5 u j /Q _ i Z. r- ^ Z. g \ I OM rA /
A(z,t) = bE'a ZH < z < ZE , (18-B)
A(z,t) = c z < ZH (18-C)
- 60 -
-------
where zr and zu are the elevations at which E - 10~ meters" and
t n
E = (b/c) ' , respectively. The value of n is chosen so that
-n(z - ZF) _ b nn-6x -a
E - ^-(10 )
Equation 18-A is based on the theory of wind mixing in the
epilimnion. The value chosen for z^, although somewhat arbitrary at
this point, was chosen on the basis of the log-log plots of Figure 12,
the rationale being that the linear log-log relationship between
A(z,t) and E should cease to exist at about the same value of E in the
epilimnion as that observed in the hypolimnion. The proper choice of A
is purely speculative at this time; however, some insight into estimating
its value may be derived from research reported by Hutchinson (3).
Inspection of the linear portion of the eddy conductivity-
stability curves in Figure 12 indicates that, while the upper envelope
curve has a slope of -1.0 in this region, the average slope for all curves
is about -0.7. Consequently, the value proposed for _a is 0.7. It is
worth noting that this same value for the average slope was observed for
Detroit Reservoir.
The value of b^ in equation 18-B was purposely left unspecified
because there is some question as to v/hether this value might not be a
function of the hypolimnetic turbulence level. It would seem that the
greater the hypolimnetic turbulence, the larger the value of E at which
the log-log relationship between A(z,t) and E begins to exist. Hhile no
clear evidence of this fact was found in Hungry Horse, it must be borne in
mind that this reservoir has a very small discharge to volume ratio and
hence this variation would not be as apparent as in a reservoir with a larger
ratio. On the other hand, Detroit Reservoir, which has an average dis-
- 61
-------
charge-volume ratio, three times larger than that of Hungry Horse
_ Q "I O _ 1
(see Chapter II), was found to have a value of 2.0 x 10~ m ' sec
for b_. Thus, for the present, it is suggested that the value of b_ be
taken as 1.5 x 10"8 m1'3 sec"1.
For reasons previously given, A(z,t) in equation 18-C was
described as constant rather than by equation 17. In addition to these
previous arguments, consideration of the limited amount of time variation
in the hypolimnetic value of the eddy conductivity coefficient coupled
with the small temperature gradients in this region tend to negate the
need for a temporal and/or spatial description of the coefficient in deep
reservoirs. On the basis of the Hungry Horse analyses, it is suggested
that until more experience is gained with the hypolimnetic value of AM,
-42-1
that the value be assumed to be 2.5 x 10 m sec . This value is in
agreement with average values observed in deep lakes (3) and the oceans(17)
1"2 -1
A value of 0.3 x 10 ' m sec for AM* which was found in Detroit
n
Reservoir, lies within range of minimum observed values.
In summary, the following values are recommended for the
parameters in equation 18:
a = 0.7,
b = 1.5 x 10"8 m1'3 sec"1 ,
c = 2.5 x 10~4 m2 sec ~1 .
Finally, in the absence of a better definition, it is suggested
that the value of AQ be taken as equal to c and that n be taken as
zero, i.e., AE(Z) = c.
- 62 -
-------
MODEL SENSITIVITY TO THE DIFFUSION FUNCTION
The proposed functional form of the eddy conductivity coef-
ficient, as given in equation 18, contains four parameters, A , a, b,
and c. In order to acquire some feeling for the sensitivity of the
deep reservoir model to variations in these parameters, three test simu-
lations were performed with A(z,t) defined as follows:
1. by the upper envelope curve shown in Figure 12 with
a = 1.0,
b = 1.8 x 10"9 m sec"1 ,
c = 1.0 x 10"3 m2 sec"1 ;
2. by the average curve shown in Figure 12 with values of a,
b, and c as given above; and
3. by the lower envelope curves shown in Figure 12 with
a = 0.7,
b = 3.1 x 10"9 m1'3 sec"1 ,
c = 1.5 x 10"5 m2 sec"1 .
Typical temperature profiles produced in these three simula-
tions are shown in Figure 15. Examination of this Figure reveals
that while the average curve gives the best representation of the ob-
served profiles, use of the lower envelope curve results in computed profiles
that, for the most part, are within 1 °C of the observed values. The upper
envelope curve, on the other hand, produces the worst simulation; one
that is not acceptable within the reasonable limits of accuracy.
These simulations, while not providing conclusive evidence re-
garding the sensitivity of the model to individual parameter variation,
- 63 -
-------
/I
en
LEGEND:
TEMPERATURE
PREDICTED BY
1050
•LOWER ENVELOPE
-PROPOSED FUNCTION
•UPPER ENVELOPE
OBSERVED
TEMPERATURE
950
17 19 21 23
TEMPERATURE, °C
FIGURE 15,
Effect of Eddy Conductivity Parameter Variation on Simulated Reservoir
Temperature Profiles; Hungry Horse Reservoir, 1965
-------
are constructive in that they show that the model is not over-sensitive
to parameter variations within the ranges defined by the lower envelope
curve and the proposed curve. Within the parameter ranges defined by
the proposed curve and the upper envelope curve, however, parameter
variation should be approached with caution since the model appears to
be more sensitive to the parameter values within these ranges.
At the request of the FWPCA Study Team, one additional simu-
lation was run in order to demonstrate the difference between the
functional form of the eddy conductivity coefficient as proposed in
equation 18 and the functional form as proposed in the original model
development (see equation 12). One deviation from this original form was
made in that A(z,t) above the thermocline was taken as constant and
smaller than that value used below the thermocline. This modification
was made in order to retain sufficient heat within the surface layers of
the reservoir so that the simulated surface temperatures would be
approximately the same as those observed. At the suggestion of the
Study Team, the following values of A(z,t) were assigned:
_ c p _ 1
A(z) = 1.0 x 10" m sec above the thermocline; and
A(z) = 2.5 x 10 m sec" below the thermocline.
Typical profiles resulting from this simulation are illustrated
in Figure 16. In deference to the originally proposed function, there
are other combinations' of values for A(z) which will produce better
results. However, this simulation does illustrate the difficulties
encountered in using equation 12 to describe the functional form of
A(z,t). In particular, it demonstrates the difficulty in simulating the
clinolimnetic portion of the temperature profile as a result of assuming
A(z,t) constant in this region.
- 65 -
-------
1090
1
cr>
3 5 7 9 II 13 15 17 19
LJ
_1
LU
1030
970
LEGEND:
TEMPERATURE PROFILE
USING STEP FUNCTION
OBSERVED TEMPERATURE
•THERMOCLINE
13 15 17 19 21 23
TEMPERATURE, °C
FIGURE 16. Reservoir Simulation with A(z,t) Defined by Equation 12; Hungry
Horse Reservoir, 1965
-------
VERIFICATION OF THE DEEP RESERVOIR MODEL
Verification of the applicability of the deep reservoir model
to reservoirs in the Northwest, and a final test of the functional form
of the eddy conductivity as proposed in equation 18, were achieved
simultaneously by the simulation of Hungry Horse Reservoir for the
period 11 May to 20 October, 1965. The values used for the eddy con-
ductivity parameters a_, b_, and £ were the recommended values; A was
taken equal to c and n = o. The time increment between successive
computations was one day. The results of this simulation, which are
illustrated graphically in Figures 17, 18, and 19, show good agreement
between the simulated and observed data.
Figure 17-B shows a comprehensive summary of the behavior of
Hungry Horse Reservoir as observed and as predicted by the model during
the 173-day period. Isothermal lines, which are superimposed on an
elevation-time field, give a pictorial representation of the changing
thermal energy distribution throughout the simulation period. The primary
value of this plot for comparison purposes is to show the time and space
relationships between the observed reservoir behavior arid the model simula-
tion. The relationship between computed results and observed results at
any point in time is best illustrated by comparison of the thermal pro-
files which are shown in Figure 19.
The time-space relationship shown in Figure 17-B adequately
demonstrates the capability of the model to describe the same progression
of events as were observed in the reservoir. The primary difference
between the model results and the observed behavior is that the model
tends to enter the cooling cycle about one to two weeks early. The exact
- 67 -
-------
A. TOTAL DISCHARGE FLOW a DOWNSTREAM TEMP
o
o
300.00
250.00
200.00
150.00
100.00
50.00
0.00
DISCHARGE
OBSERVED TEMP.
TEMPERATURE
T
^N fi
V"~^ vl.Asr>Nl^ "' A,
#C1 .^|K?fe^4
u
JUN.
JUL
AUG.
SEP
OCT.
30.00
25.00
20.00
15.00
10.00
5.00
0.00
O
o
LU
cr
13
1-
<
or
LJ
Q_
^
UJ
B. TEMPERA TURE ISOTHERMS - RUN I DENT. NO. 5
cu
T" 20
20 20 I 18 16 14 12 10
- ~--:~c^--£r?^^"t"._.
3378.80 5
__ COMPUTED ISOTHERMS-
MEASURED ISOTHERMS-
JUN.
JUL.
FIGURE 17. Simulated and Observed Thermal Regime
for Hungry Horse Reservoir, 1965
- 68
-------
A. OUTLET AT SURFACE
250.00
200.00
150.00
100.00
50.00
0.00
300.00
250.00
CO
^ 200.00
O
^ 150.00
O
_J
^ 100.00
50.00
0.00
300.00
250.00
200.00
150.00
100.00
50.00
0.00
-r™" i r- ' — r" T " "
/TEMPERATURE
•xX/^-x
X
^-x -
/DISCHARGE
if
JUN. JUL. AUG. SEP OCT.
B. OUTLET AT ELEVATION 1013
^^-DISCHARGE
|Y
\ .TEMPERATURE ii
I /v A P^-'MJ,_
A _ \ i ii^^ " ''~*j ^n / • / 1
_f — A. » | ~\ ^ i('"^-sA i
^,_^iL^^.^i^L._..-^L,.™™lL^.^,™i-,™™:r
JUN. JUL. AUG. SEP OCT.
C. OUTLET AT ELEVATION 975
—
-
-
-
^^— — T E M P E R A T U R E
ts^^
^^— Dl SCH ARG E
JUN. JUL. AUG. SEP OCT.
25.00
20.00
15.00
10.00
5.00
0.00
30.00
25.00
20.00
15.00
100.00
5.00
0.00
30.00
25.00
20.00
15.00
10.00
5.00
0.00
O
o
cr
a:
UJ
CL
UJ
h-
FIGURE 18. Breakdown of Outlet Discharges and Computed Outlet
Temperature for Hungry Horse Reservoir, 1965
- 69 -
-------
MAY II, 1965
MAY 26, 1965
Tl.CO I.OC
u.w m.oo is.x IB.oo zo.ct
YDO j.co 4.co
JUNE 9, 1965
JUNE 23, 1965
TIHPE?^ILRE. DEC. C
14.DC 16.DO IB.00 ZC.K
1, 1965
Z.DC 4.ID fl.OO '.CO
FIGURE 19, Simulated and Observed Temperature Profiles for Hungry Horse Reser-
voir, 1965
(Continued next page)
- 70 -
-------
AUGUST 4, 1965
TEMPERflTUBE. OEG. C
lfl.CC 20.0C
SEPTEMBER 1,1965
SEPTEMBER 16,1965
1B.OO 20.00 Tl.OD 2-CO -J.OT
OCTOBER 5, 1965
OCTOBER ZO, 1965
FIGURE 19 (cont'd). Simulated and Observed Temperature Profiles for Hungry
Horse Reservoir, 1965
- 71 -
-------
reason for this cannot be stated due to the complexity of the interaction
between the reservoir and the meteorological factors. One reasonable
explanation, however, could be that the surface of the reservoir becomes
more turbid in the fall, thus trapping more short-wave solar radiation
in the surface layers than occurs earlier in the year. Such a condition
would warm the surface layers and consequently retard the rate of surface
cooling. However, even with the early cooling. Figure 19 shows that the
maximum variation in the computed and observed epilimnetic temperatures
is less than 2°C.
It should be noted that while the spatial variations between
the computed and observed 4°C and 6°C isotherms is quite large, the
actual difference between the observed and computed temperatures is
small (Figure 19 indicates the differences to be less than 0.6°C).
This apparently anomolous behavior is due to the very small thermal
gradients in the hypolimnion.
Figure 17-A supplies complimentary data to Figure 17-B in
that it gives the time history of the net thermal and hydraulic character-
istics of the reservoir discharge. The observed temperature line shown
in the Figure is the average daily stream temperature as recorded 1.5
miles below the dam. Except for the month of June and the first part of
July, the recorded temperatures are always lower than the computed tempera-
tures. Communication with persons familiar with the gage indicates that
it probably contains a negative bias. The high recorded temperatures
during June apparently resulted from the temperature probe being in slack
water during the low discharge period. Even though the absolute agree-
ment between the observed and computed temperatures cannot be compared,
- 72 -
-------
the observed temperatures are valuable in assessing the capability of
the model to produce the observed day-to-day behavioral pattern of the
downstream temperature. Examination of Figure 17-A indicates that the
behavioral pattern of the predicted temperature is in excellent agreement
with that observed, thus demonstrating the model's capability to simulate
the day-to-day variation in the downstream temperature.
Figure 18 provides a breakdown of the composite information
contained in Figure 17-A. In this simulation run, the method of flow
withdrawal from the reservoir was that used in the original derivation
of the model, i.e. the withdrawal depth is taken to be the thickness of
the element whose elevation corresponds to the elevation at the centerline
of the outlets (in this case, the element thickness was one meter).
Consequently, the temperatures shown in Figure 18 are the reservoir
temperatures at the elevation of the outlet. These particular plots will
be more useful in comparing the element withdrawal technique to the
selective withdrawal technique which is the subject of the following
chapter.
In summary, the verification run as described above adequately
demonstrates the validity of equation 18 in describing the internal
mixing process. Moreover, this simulation serves as verification that
the model is generally applicable for the simulation of deep reservoirs
in the Pacific Northwest.
- 73 -
-------
V. SELECTIVE WITHDRAWAL STUDIES
THEORY
When water is withdrawn from a stratified impoundment, it
does not necessarily follow that the entire fluid bulk will contribute
to the outflow as is the case for a non-stratified body. The theoretical
analysis of Yih (19) showed that there is a tendency for flow separation
when the dens i metric Froude numbe'r becomes less than I/TT. That flow
separation does, indeed, occur at low Froude numbers has been confirmed
experimentally by Debler (20).
For the flow situation shown in Figure 20, the densimetric
Froude number F is defined by Debler as
where
q = volumetric discharge per unit width of the channel;
d = height of the dividing streamline far back from the
outlet
po = density at the channel bottom;
0 = If ; and
g = gravitational constant.
- 75 -
-------
|
/
/
/
REG! OfJ
FLOWING
REGION
DIVIDING STREAMLINE \SLOPE =~P
DENSITY PROFILE \
d_P_
= dz
JL
D
FIGURE 20. Schematic Sketch of Flow Separation in
Stratified Flow
Through the performance of a set of experimental runs, Debler established
that the Froude number of the divided flow was a constant, with a value
of about 0.24. On the basis of this result, the withdrawal depth, d, can
then be defined as
d = 2.0
(20)
This equation is recommended by Brooks and Koh (21) for determination of
the withdrawal thickness close to -the dam (Brooks and Koh recommend
1.9 for the constant rather than 2.0). Farther back from the dam, they
recommend using the results of Koh's theory (reported in the same publica-
tion) which incorporates viscous broadening of the withdrawal layer for
the case of steady turbulent flow. However, since transients in the
withdrawal currents caused by fluctuating release rates are very persistent,
- 76 -
-------
they can significantly alter the steady flow solutions, especially far
back from the dam. Consequently, it was felt that the use of equation 20
would provide as good an approximation to the withdrawal zone as could
be obtained by application of the more sophisticated viscous broadening
theory.
In applying equation 20 to an operating reservoir, reference is
made to Figure 21 which shows a dam with three outlets. Each outlet
has a volumetric discharge Q. To convert these discharges to a discharge
per unit width, the reservoir width, w_, at the elevation of the outlet
is used. Equation 20 for this case is assumed to define the half-width
of the withdrawal layer. The value of p is taken as the density at the
centerline of the outlet and 3 is taken as the negative value of the
density at the centerline of the outlet.
TEMPERATURE PROFILE-*
A
FIGURE 21. Schematic Diagram of the Application of
Selective Withdrawal Theory to a Stratified
Impoundment
- 77 -
-------
By the use of these definitions for p and g, it is observed that the
p -, -1
ratio — = (—3) is the reciprocal of the stability, E, at the
3 po
centerline of the outlet. Thus equation 20 can be rewritten as
1/2
d = 2.0 [] (21)
Applying equation 21 to the outlet configuration shown in Figure 21.
the value of q is taken as follows:
Qs
surface outlet q = —
1/2 QA 1/2 QR
submerged outlets q -
u(z) - \4- U'
WA WB
To apportion the flow within the withdrawal zones, the inviscid
theory of Yin (22) is used in which the velocity within the flow field
is defined as
(22)
where U' is a constant called the "associated velocity." However,
considering the degree of density stratification encountered in prototype
reservoirs, the variation in the ratio — over the withdrawal depth is
P v
insignificant; hence, the velocity within the withdrawal zone is defined
simply as
u(z) = constant. (23)
- 78 -
-------
In applying equation 23 to the flow field, it is assumed that
the principle of superposition holds for regions in which the withdrawal
layers overlap. If the effect of channel side-slope is neglected, the
withdrawal layer configuration shown in Figure 21 has a spatial velocity
distribution given as
»M- S^T- zB + dB
-------
in Figure 22 for the case in which the computed streamline falls below
the reservoir bottom.
1.
•W
1
WITHDRAWAL LAYER COMPUTED
BY EQUATION 21
REVISED I
VELOCITY
DISTRIBUTION
REVISED WITHDRAWAL LAYER ^
~^7Z£~
'2^7
ORIGINAL
VELOCITY
DISTRIBUTION
FIGURE 22. Withdrawal Region for Case of Dividing
Streamline Falling Below Reservoir Bottom
A second situation, which can be encountered, is that the
dividing streamline, computed for an outlet on one side of the therrno-
cline, falls on the other side of the thermocline. An example of such a
case is shown in Figure 23. The objection to this condition is based on
the fact that the withdrawal depth is a function of stability which is
assumed constant over the withdrawal depth. The case illustrated in
Figure 23 is clearly in violation of this assumption. Since the withdrawal
depth varies inversely with the stability, the upper withdrawal region
should be smaller than the computed value. It is reasonable and con-
venient to limit the upper withdrawal zone to the elevation of the
- 80 -
-------
thermocline as shown in the Figure. A similar correction is made for
outlets in the epilimnion whose lower streamline falls below the thermo-
cline. The flow distribution for these cases is adjusted in a manner
analogous to that for the case of the dividing streamline crossing a
reservoir boundary.
TEMPERATURE PROFILE
REVISED
WITH DRAWAL
WITHDRAWAL
LAYER
FIGURE 23. Withdrawal Region for Case of
Computed Dividing Streamline
Crossing the Thermocline
- 81 -
-------
A special situation is encountred with the onset of fall
cooling when epilimnetic mixing by natural convection becomes pre-
dominant. The result of this mixing is to cause the reservoir in the
epilimnetic region to be isothermal as illustrated in Figure 24; hence,
the stability in this region is zero and Debler's criteria is no
longer applicable to outlets above the thermocline. To account for
this situation, the theory of Craya (23) as reported by Yin (22) has
been used.
o^x
UJ O
cco^E
T —
FIGURE 24. Effect of Convective Mixing on
the Thermal Profile
To use Craya's approach, it is assumed that the density
structure of the reservoir is a step function whose value is p above
the thermocline and p + Ao below the thermocline as illustrated in
Figure 25.
- 82 -
-------
THERMOCLINE-
P+AP
ASSUMED
DENSITY
STRUCTURE—-
OBSERVED
DENSITY
-STRUCTURE
P—e»-
FIGURE 25'. Assumed Reservoir Density
Structure for Application of
Craya's Criteria
Using the notations of this Figure, the critical epilininetic
discharge at which water will begin to be withdrawn from the hypolimnion
of the reservoir is
2h < dT ,
(25)
for submerged outlets, and
= 0.75
(26)
- 83 -
-------
for surface outlets. The velocity distribution in the epilimnion is
assumed to be uniform.
If q is greater than q for any outlet or combination of
\*
outlets, Debler's criteria (equation 21) is used to compute the depth of
withdrawal from the hypolimnion, taking q - qc for the discharge and
the value of E immediately below the thermocline for the stability.
The velocity distribution is then adjusted so that it is uniform through-
out the entire withdrawal depth. While there is no precedent for using
Debler's criteria for the excess flow, there does not appear to be any
existing theory which is more appropriate in this kind of situation.
Although some parts of the selective withdrawal theory related
above might be challenged on theoretical grounds, it must be recognized
that the approach to the incorporation of selective withdrawal theory
into the deep reservoir model has of necessity had a pragmatic orientation
due to the need to produce a functional tool which could provide required
answers when applied to prototype situations. Nevertheless , WRE believes
that the theory related above represents the best approach available at
the present time. Moreover, it points out the need for additional
research and development of a practical nature in the area of withdrawal
from stratified impoundments. The incorporation of the above theory, as
it stands, is regarded as being another step forward in refining the deep
reservoir model. The use of this theory should enhance the model's
capabilities particularly in simulating the individual discharge tempera-
tures from the reservoir outlets.
In order to test the validity of the proposed selective with-
drawal scheme, it was desirable to compare computed withdrawal tempera-
tures with observed temperatures. The only observed withdrawal temperatures
- 84 -
-------
on Hungry Horse were the cooling water temperatures for the turbines and
the accuracy of these temperatures was not known. Consequently, Fontana
Reservoir was chosen to test the theory since the observed data was more
reliable. The data for the computations was taken from reference 24.
The comparison of the computed turbine intake temperatures with the observed
turbine cooling water temperatures for Fontana Reservoir is shown in
Figure 26 (the cooling water temperatures had been correlated and
adjusted to the reservoir temperatures during a period when the reservoir
was isothermal). It appears that the computed outflow values contain a
small amount of bias, tending on the average to be about 1°F too cool.
Very rarely are they biased more than 2°F. At the present time the reason
for this bias is not known; however, it is felt that the predictions are
close enough to the observed temperatures to adequately satisfy the
modeling requirements.
As a further test, a few trial computations were made on Hungry
Horse Reservoir and Detroit Reservoir. The results shown in Table 2
indicate that the prediction is satisfactory.
SIMULATION OF HUNGRY HORSE RESERVOIR USING SELECTIVE WITHDRAWAL THEORY
Having incorporated the selective withdrawal theory into the
deep reservoir model, the simulation run illustrated in Figures 27 and
28 was executed. The only difference between this run and the verifica-
tion run in Chapter IV is that the verification run used the original
method of element withdrawal while the run illustrated here had the
selective withdrawal theory incorporated. To facilitate the comparison
of these two simulations, Figures 17 and 18 have been overlaid on Figures
27 and 28, respectively,
- 85 -
-------
UJ
46 48
50 52 54 56 58 60 62 64
COOLING WATER TEMPERATURE,°F
66 68 70
FIGURE 26. Comparison of Outflow Temperatures Computed from Debler's
Criteria with Measured Cooling Water Temperatures (Fontana
Reservoir)
- 86 -
-------
TABLE 2
HUNGRY HORSE RESERVOIR
TURBINE INTAKE TEMPERATURES
DATE
21 July 1965
18 August 1965
COMPUTED, °C
4.4
4,3
COOLING WATER TEMPERATURE, °C
4.5
5.0
DETROIT RESERVOIR
TAILRACE TEMPERATURES
DATE
20 May 1964
7 November 1958
10 November 1958
TURBINE INTAKE
TEMPERATURE
COMPUTED, °F
42
54
52
SPILLWAY TEMPERATURE
COMPUTED,°F
47
NO DISCHARGE
NO DISCHARGE
TAILRACE TEMP.5°F
COMPUTED OBSERVED
44
54
52
43
54
53
- 87 -
-------
Examination of the results indicates that the effects of in-
cluding the selective withdrawal theory into the model are something
less than dramatic. The most obvious result is a minor alteration
of the thermal regime of the reservoir in that the vertical rate of heat
transfer is slightly increased, i.e., the isotherms for the selective
withdrawal case are somewhat deeper than in the element withdrawal case.
In addition, the onset of the convective cooling cycle is slightly
retarded.
The complexity of the reservoir mechanics precludes any simple
explanation for this observed behavior. It can only be hypothesized
that the selective withdrawal scheme induces larger thermal gradients
into the reservoir with the result that the rate of turbulent heat
transfer is increased.
It was expected that the use of selective withdrawal would
cause a decrease in the computed net outflow temperature during the
period of spillway operation, since some of the cooler water below the
surface would be included in the spillway discharge; however, this was
not the case. In fact, the selective withdrawal scheme produced slightly
higher outflow temperatures than did the element withdrawal scheme. The
reason for this is hypothesized as follows. During the period of spill-
way operation, the element withdrawal scheme takes all of the surface
spill from the surface element of the reservoir, thus tending to expose
the cooler underlying water. The selective theory, on the other hand,
causes a portion of the spill to be drawn from elements below the surface;
consequently, not as much of the warm surface water is removed, and the
reservoir's surface has a tendency to remain warmer. Unless the spillway
- 88 -
-------
o
o
u_
300.00
250.00
200.00
150.00
100.00
50.00
0.00
A. TOTAL DISCHARGE FLOW a DOWNSTREAM TEMP.
DISCHARGE
TEMPERATURE
ELEMENT
WITHDRAWAL
TEMPERATURE,
SELECTIVE
WITHDRAWAL .
JUN.
JUL
AUG.
SEP
OCT.
30.00
25.00
20.00
15.00
10.00
5.00
0.00
o
o
UJ
o:
13
h-
<
or
LjJ
Q_
•^
UJ
B. TEMPERATURE ISOTHERMS - RUN I DENT. NO. 5
CD
E
O
UJ
_i
UJ
1090.00
1075.00
1060.00
1045.00
1030.00
1015.00
1000.00
985.00
970.00
955.00
940.00
1416 18 20 20
*~ .-^-—T—^ 1—•—cr ] i i\—i
SELECTIVE WITHDRAWAL -*\\*-ELEMENT WITHDRAWAL
PREDICTIONS
PREDICTIONS
\
3576.00
3526.70
3477.40
3428.10
3378.80
3329.50
3280.20 >
UJ
JUN.
JUL.
AUG.
SEP.
OCT.
3230.90
3181.60
3132.30
3083.00
UJ
FIGURE 27.
Simulated Thermal Regime of Hungry Horse
Reservoir, 1965, Using Selective Withdrawal
Theory
- 89 -
-------
A. OUTLET AT SURFACE
5UU.UU
250.00
200.00
150.00
100.00
50.00
Or\r\
T "| | | |
TEMPERATURE,-^
SELECTIVE >.
WITHDRAWAL /r\J\^.
X\ f f»-" r "*vk*»
^ s ^^
TEMPERATUREr"^^ VV
ELEMENT "^i
WITHDRAWAL ^=*--
DISCHARGE— ^^^
, , .r-^-T ^--^7- — " — — ^
^>v.\j^
25.00
20.00
15.00
10.00
5.00
n r>n
JUN. JUL. AUG. SEP. OCT
B. OUTLET AT ELEVATION IOI3
300.00
250.00
CO
^ 200.00
o
^ 150.00
O
^ 100.00
c.r\ on
^j\j \j\j
0 00
1 1 1 1 i
_
_ —
nicrwARrr K /TEMPERATURE,
DISCHARGE-^N /ELEMENT
\\ 1 WITHDRAWAL •
TEMPERATURE, i \ 1 i M
- SELECTIVE 1 \ \ li 'i'
WITHDRAV/AL-x j \ \ ; lif-.Aip <
r I__XA— — -X l4 - j- '^/r —
iV « A I / 1 "^' 'l 1
30.00
25.00
20.00
15.00
100.00
s nn
•~>.\j\j
ri r\ri
JUN. JUL. AUG. SEP OCT.
C. OUTLET AT ELEVATION 975
•ij\J \J . \J \_/
250.00
200.00
150.00
100.00
50.00
0 00
III II
— —
^- TEMPERATURE,
^^ SELECTIVE
/ WITHDRAWAL
/ TEMPERATURE,
/ /ELEMENT
£ / WITHDRAWAL
/c: x^rr--^^^
^ — DISCHARGE
(, — ^ c_.j — v | | |
^J^J. WW
25.00
20.00
15.00
10.00
5 OO
^«* . v V/
r\ r\r\
JUN. JUL. AUG. SEP OCT.
u
o
^
LU
cr
^
o:
UJ
CL
LlJ
(—
FIGURE 28. Simulated Discharge Temperatures for
Outlets as Con.: jtad with Selective Wi
Hungry Horse Reservoir, 1955
- 90 -
Theory;
-------
discharge is sufficiently high, the withdrawal layer will not be deep
enough to result in a new withdrawal temperature that is less than that
of the element withdrawal case.
The examination of Figure 29 reveals one other observation that
is worth discussing. During the period of operation of the discharge
tube at elevation 975, the selective withdrawal case produces outflow
temperatures significantly higher than does the element withdrawal
\
scheme. The reason for this is that the stability gradient at this point
is sufficiently low to cause withdrawal from the entire hypolimnion;
consequently, the outflow temperature for the selective withdrawal case
represents the mean temperature of the entire reservoir below the thermo-
cline, whereas the outflow temperature for the element withdrawal scheme
represents only the temperature of the reservoir at elevation 975.
In summary, while the simulation run incorporating the
selective withdrawal theory did not produce the dramatic results expected,
it is felt that this was due to the physical nature of the reservoir
and its mode of operation, not to a failure of the selective withdrawal
theory. This conclusion becomes more clear in the following chapter
where two test cases are simulated using the net outflow temperature
from the reservoir as the criteria for the operation of the outlets.
- 91 -
-------
VI. MULTIPLE OUTLET STUDIES
As a demonstration of the deep reservoir model's versatility
in assessing the consequences of reservoir regulation for downstream
temperature control, two test cases were simulated on Hungry Horse
Reservoir. These tests were conceived mutually by the FWPCA Study
Team and WRE and are as follows.
TEST CASE I
Using the existing outlet configuration in Hungry Horse Reser-
voir, operate the outlets on a daily basis in such a manner as to pro-
duce a prescribed net outflow temperature.
TEST CASE II
HypotheticaIIy, add a second turbine intake to the reservoir.
On a daily basis, operate the spillway in the same manner as it was
historically operated in 1965 and distribute the remainder of the out-
flow between the two turbine intakes in such a manner as to produce a
prescribed outflow temperature.
- 93 -
-------
For both test cases, the total reservoir outflow was taken as
the historical outflow in 1965. The downstream object temperature for
the two cases was taken as the observed temperature of the major reser-
voir tributary, i.e., the temperature of the South Fork of the Flathead
River above Twin Creek for the year 1965. The results of these two test
simulations are reported below.
TEST CASE I
Test Case I was not the product of whimsical imagination.
Rather, it was designed around the fact that, since the construction of
Hungry Horse Dam, the reservoir release temperatures have been so cold
that the native fish population below the dam has been eliminated. By
operating the reservoir outlets in such a manner as to produce a release
temperature equal to that of the natural stream, it would be assured that
the fish population could be re-established. In addition, the analysis
of this type of operation would reveal the hydropower tradeoff required
to maintain the downstream fishery.
The simulation results of Test Case I are illustrated in Figures
29 and 30. In order to compare this simulation with the historical case,
the results of the selective withdrawal simulation (Figures 27 and 28)
have been overlaid. Inspection of Figure 29-A reveals that excellent
agreement was obtained between the reservoir release temperature and
the object temperature. This result demonstrates that it is possible to
operate the reservoir in a manner such that the required downstream
temperature control can be achieved. One exception, however, must be taken
to this general statement. While the test case shows that spillway
- 94 -
-------
o
£
o
300.00
250.00
200.00
150.00
100.00
50.00
0.00
A. TOTAL DISCHARGE FLOW a DOWNSTREAM TEMP
OBJECT
TEMP.
I
-HISTORICAL TEMP
-TEMPERATURE
DISCHARGE
JUN.
JUL
AUG.
SEP.
OCT.
30.00
25.00
20.00
15.00
10.00
5.00
0.00
O
o
UJ
a:
z>
<
rr
LU
a_
^.
UJ
B. TEMPERATURE ISOTHERMS - RUN I DENT. NO. 6
O
h-
<
UJ
_J
UJ
JUN.
JUL.
AUG.
SEP
OCT.
FIGURE 29.
Simulated Thermal Regime of Hungry Horse Reservoir
1965, with Modified Reservoir Release Scheme
- 95 -
-------
operation was required in the early part of the summer, this could not
have been possible physically, since the water surface was below the
spillway crest prior to 24 June. Consequently, in practice, the object
temperature would not have been met during the period 11 May to 24 June
but rather, the outflow temperature would have been the same as that recorded
historically. This exception does not detract, however, from the general
result.
Examination of Figure 29-B indicates that changing the operation
schedule of the outlets had no significant effect on the thermal regime
of the reservoir except to raise slightly the reservoir isotherms toward
the middle and end of the simulation period.
To evaluate the magnitude of the power tradeoff required to main-
tain the object temperature, reference is made to Figure 30 which illus-
trates the release rate and temperatures of the individual outlets. The
outlet at elevation 975 is not shown since it was not operated in this
simulation. If the amount of power production is taken as a linear
function of the turbine discharge rate, Figure 30 indicates that the
hydropower production on a day-to-day basis has been reduced approxi-
mately 40 percent from its historical value over the period from 7 July
to 27 August; thereafter the percentage loss decreases until by 13
September the loss is insignificant. The loss of revenue from hydro-
power, as measured in terms of turbine discharge rates, has, so to
speak, gone over the spillway. These dollars lost must be weighed
against the dollars (and intangibles) gained through re-establishment
of the fishery below the dam, and other benefits of downstream tem-
perature control.
- 96 -
-------
300.00
250.00 -
200.00 -
A.OUTLET AT SURFACE
HISTORICAL TEMPERATURE
HISTORICAL DISCHARGE
,-_ — i A,
-' :--~~- --'
30.00
- 25.00
-20.00
- 15.00
- 10.00
AUG.
SEP.
OCT.
-5.00
0.00
o
o
O
B. OUTLET A T ELEVA TION IOI3
300.00
250.00
200.00
150.00
100.00
50.00
0.00
-HISTORICAL DISCHARGE
-DISCHARGE
-TEMPERATURE
(Historical Values
'•i
!!
JUN.
JUL.
AUG.
SEP.
OCT.
UJ
tr
I
o:
UJ
CL
30.00 S
Ld
I-
25.00
20.00
15.00
10.00
5.00
0.00
FIGURE 30. Discharges and Discharge Temperatures for Individual
Outlets with Modified Reservoir Release Scheme; Hungry
Horse Reservoir, 1965
- 97 -
-------
The examination of Figure 30-A reveals that, at the beginning of
tfxe high spill period in July, the spill temperature for Test Case I was
significantly lower than the historical value, even though the thermal
regime of the reservoir was identical during this time for both cases
as evidenced by Figure 29-B. This observed phenomenon is the result of
applying selective withdrawal theory to the surface outlets and pro-
vides the culmination of the discussion in Chapter V which alluded to
this Figure. In this case, the withdrawal layer was of sufficient
depth (six meters) to produce a net outflow temperature significantly
lower than the observed surface water temperature.
TEST CASE II
The second test case was hypothetical in that it was assumed that
a second turbine intake was available for power generation. The objective
of this test case was to determine if it was possible to distribute the
1965 historical power demand between the two sets of turbines and still
achieve the release temperature objective. Consequently, the hypothetical
intake was located at elevation 1055 which is the highest elevation at
which it could be placed and remain submerged over the simulation period.
The results of this simulation are illustrated in Figures 31 and 32.
Again, as in Test Case I, the 1965 historical simulation with the selec-
tive withdrawal scheme has been overlaid for comparison purposes.
Inspection of Figure 31-A reveals that the object temperature could
not be achieved during June, July, and August. While the release temp-
eratures during this period were one to two degrees Centigrade higher
than the historical release temperature, they fell from one to six degrees
Centigrade below the object temperature. Hence, it is concluded that
1965 power demand and the object temperature could not have been stmul-
- 98 -
-------
taneously satisfied with two fixed turbine intake locations. If the
power demands were to be met, it appears that a set of multiple intakes
would have been required in order to meet the object temperature.
Figure 31-B indicates that the inclusion of a second turbine in-
take did not significantly alter the thermal regime of the reservoir.
The breakdown of temperatures and discharges through the indi-
vidual outlets is illustrated in Figure 32. It is apparent from
Figures 32-B and 32-C that the reason why the object temperatures could
not be met during the summer months was that the reservoir temperature
at elevation 1055 was not high enough. Moreover, the only time of the
year when the lower turbine intake was required was in the fall when the
upper portions of the reservoir were warmer than the object temperature.
The two test cases described above have been demonstrative of
the types of analysis which can be performed for re-evaluating the
operating rules for existing reservoirs and for aiding in the design
of outlet placement in proposed reservoirs. Such analyses (and there
are many that could be performed depending upon the particular problem
which may be encountered) would not be possible without a tool like
the deep-reservoir model. In fact, this is the real value of the model. Its
utility lies not in its ability to reproduce historical data but rather in
its potential for planning and management through its prediction capabilities
- 99 -
-------
A. TOTAL DISCHARGE FLOW a DOWNSTREAM TEMP.
HISTORICAL TEMP.
TEMPERATURE
OBJECT TEMP.
0.00
JUN.
JUL
AUG.
SEP.
OCT.
0.00
o
o
LJ
cc
Z>
1-
<
or
UJ
o.
^
UJ
B. TEMPERATURE ISOTHERMS - RUN I DENT. NO. 7
CO
L_
CD
O)
e
o
i-
2
UJ
_J
UJ
1965 RESERVOIR REGIME- A
Selective Withdrawal Predictions^
JUN.
FIGURE 31.
JUL.
AUG.
SEP
OCT.
Simulated Thermal Regime of Hungry Horse
Reservoir, 1955, with Addition of Additional
Turbine Intake
- 100 -
-------
A. OUTLET AT SURFACE
250.00
200.00
150.00
100.00
50.00
Onr^
1 1 1 1 1
TEMPERATURE,
~X//V/^\ 1
v)
v_ _
DISCHARGE--^
i i r ' ' ^'i ' "^K
25.00
20.00
15.00
10.00
5.00
n nn
JUN. JUL. AUG. SEP. OCT.
B. OUTLET AT ELEVATION 1055
300.00
250.00
if)
^ 200.00
O
..
^ 150.00
O
1
^ 100.00
50.00
n nn
i i i r~ T"
-
^-DISCHARGE
S
\\
i\ TEMPERATURE^
\ \
! \ \
__ : — ' ;" f • . n
r '/•• / ^•^. i •!'. A
~~ '^ / *^ » j • '"A* y • \ A . ^
V \ / ^^ ^^ " V'*^ i/ i /- i • ~^ \
30.00
25.00
20.00
15.00
100.00
5.00
OOO
O
o
UJ
or
^
IT
UJ
DL
UJ
1-
JUN.
JUL.
AUG.
SEP.
OCT.
C.OUTLET AT ELEVATION 1013
5UU.UU
250.00
200.00
150.00
100.00
50.00
c\ r\r\
-
HISTORICAL DISCHARGE— ^^
DISCHARGE -~^^\ ^\
"HISTORICAL \ TEMPERATURE-, \A
TEMPERATURE-^^^ \ /, j{
"^\ r % ~^' ^\>j ^L
JW. WW
25.00
20.00
15 00
10.00
5/-\/~\
.UU
0.00
JUN.
JUL.
AUG.
SEP
OCT.
FIGURE 32. Discharges and Discharge Temperatures for Individual
Outlets with Addition of Additional Turbine Intake;
Hungry Horse Reservoir, 1965
- 101 -
-------
VII. SUMMARY AND CONCLUSIONS
GENERAL
The major accomplishments emanating from the efforts expended
in Part I of the research were the incorporation of selective with-
drawal theory into the model, as related in Chapter V, and the refined
functional description of the eddy conductivity coefficient as given
by equation 18. In addition, refining the method for introducing inflows
into the model, and the incorporating the minor model refinements as
described in Chapter III, have enhanced the operational sophistication
of the solution technique. It is felt that these improvements represent
a major step forward in generalizing the model's capability to simulate
deep reservoirs.
The model, at its present level of development, has essentially
eliminated the need for reliance on prototype thermal profiles for evalu-
ating the empirical coefficients in the eddy conductivity function;
consequently, its potential as a prediction and design tool for proposed
reservoirs has been substantially increased. The two example test cases
in Chapter VI served to demonstrate this potential in design as well as
to indicate the value of the model for re-evaluating the operating rules
of existing reservoirs. Thus, the deep reservoir model, as developed to
this point, meets the requirements imposed by objectives 1 and 2 stated
- 103 -
-------
in the INTRODUCTION, and is considered to be adequate for the types of
analyses contemplated by the FWPCA Study Team. The supplement to this
report supplies pertinent information regarding the operation of the
model and the interpretation of its output.
LIMITATIONS
While the deep reservoir model has been significantly generalized
and upgraded from its originally derived state, it is not without limi-
tations and these limitations must be omnipresent in the mind of the
user if the results of simulation runs are to be properly interpreted.
The primary limitation of the model is that it is one-dimensional. While
the solution describes a three-dimensional body of water, keeping complete
inventory of all water and energy in the system, it permits the in-
ternal transfer of heat and mass to occur only along the vertical axis of
the reservoir. Such a solution implies that the reservoir properties are
constant, or uniform, throughout a horizontal plane. Although this
implication is reasonable for the majority of deep reservoir problems
of interest, exceptions will be encountered. On the basis of the exper-
ience gained through the several applications of the model, these excep-
tions are expected to occur when the discharge to volume ratio becomes
more than about 3 yr" , or when the reservoir geometry is atypical, e.g.
and extremely long reservoir, such as Lake Roosevelt, or a severly
branched reservoir with,significant inflows in the branches. On terms of
the densimetric Froude number as defined by equation 2 in Chapter II, which
simultaneously accounts for the discharge volume ratio and the reservoir
geometry, the applicability of the model should be questioned if F is
on the order of 0.1 or greater.
- 104
-------
A second limitation is associated with the selective withdrawal
theory. The theory, as incorporated into the model, assumes that the
reservoir outlets are slots which extend laterally across the face of
the dam. The error of this assumption is not known at present; however,
the tests on Fontana, Hungry Horse, and Detroit reservoirs, as related
in Chapter V, tend to indicate that it is not a serious problem.
Finally, the functional form of the eddy conductivity coeffi-
cient must be used with caution until more experience is gained through
the application of the model to additional reservoirs. Since time and
money constraints limited the developmental investigations to Hungry
Horse Reservoir, it was impossible to form generalized functional de-
scriptions for three of the empirical parameters.
The first of these parameters is A (t), the value of the eddy
conductivity at the water surface as a result of wind mixing. This
value is probably a function of the wind speed, wind direction, and
reservoir surface geometry. Its value could not be investigated since
no wind mixing was observed on the re'servoir.
The second parameter is the value b in equation 18-B. While a
value was obtained for Hungry Horse Reservoir and is recommended for
use at the present time, it is possible that this value could be a
function of the turbulence level of the hypolimnion, particularly in
reservoirs with larger discharge to volume ratios. Although it is
felt that this parameter requires further study, it is believed that
the use of the recommended value will not causa misleading conclusions
to be drawn in general modeling applications.
- 105 -
-------
The third parameter which was not generalized was the value of
AH(t), the hypolinmetic value of the eddy conductivity. While this
parameter was shown to vary with the inflow-outflow relationship in
the hypolimnion, the conclusions were too weak to generalize without
further study. However, as in the case of the parameter b, the recom-
mended value should suffice as a reasonable approximation until addi-
tional investigations on different reservoirs can be conducted.
RECOMMENDATIONS
Notwithstanding the limitations described above, WRE believes
that the current version of the deep reservoir model constitutes the
most generally applicable model available at the present time for pre-
dicting the thermal regime of stratified reservoirs. However, some
amount of additional refinement would greatly enhance its prediction
capabilities. This refinement is specifically concerned with the gener-
alized definition of the three empirical parameters A (t),b, and
Anft), contained in the eddy conductivity function. To this end, it is
recommended that a set of reservoirs be analyzed in order to derive gener-
alized descriptions for A (t),b, and Au(t). The reservoir set should
0 n
display the following characteristics:
I. Differing geometrical shapes;
2. Various discharge to volume ratios;
- 105 -
-------
3. Wind mixing should be evident in some of
the reservoirs; and
4. Each reservoir must have the required complement
of data, including thermal profiles, for effecting
a "backward" simulation in order to determine ob-
served eddy conductivity coefficients for the
reservei r.
The incorporation into the deep reservoir model of generalized expressions
for these three parameters would constitute the ultimate development of the
model consistent with its one-dimensional limitation.
- 107 -
-------
PART II
SIMULATION OF WEAKLY STRATIFIED RESERVOIRS
-------
VIII. INTRODUCTION
The weakly-stratified reservoir is characterized pre-
dominately by isotherms which are tilted along the longitudinal
axis of the reservoir. This tilting results primarily because
of the low residence time of the reservoir coupled with a large
length-depth ratio for the reservoir, i.e., because of the large
Froude number (see Chapter III). The predominant hydraulic
phenomenon which affects the degree of tilting is the interflow
within the reservoir which possesses sufficient inertia to disrupt
the density structure of the reservoir from its static equilibrium
state.
As noted in Chapter II, the description of the weakly-
stratified reservoir requires that the variation of the reservoir
properties be specified in two dimensions; namely in the vertical
direction, and along the longitudinal axis of the reservoir. The
approach taken herein has been to divide the reservoir longitudinally
into several segments and to assume that the heat and mass transfer
phenomena within each segment can be described by the deep reservoir
model. Thus, the weakly-stratified reservoir is envisioned as a set
of deep reservoirs connected in series. Each "reservoir" is coupled
- 109 -
-------
to the one above and below by the requirements for heat and mass
continuity within the reservoir as a whole. By solving the
vertical temperature distribution in each segment and comparing
the profiles of adjacent segments, it is possible to describe the
thermal properties of the reservoir in the tv/o required directions,
The idea of longitudinal segmenting is not unique to this
model. In particular, Jaske (25), in his recently published model
of Lake Roosevelt, has used a segmenting scheme that is nearly id-
entical to the one proposed here. Without going into detail, the
basic differences between his model and the model presented here
are:
I. the method of handling heat transfer across the
ai r-water " i nterface,
2. the description of the heat and mass transfer
between adjacent segments,
3. the formulation of the mechanisms of internal
heat transfer in the vertical direction (Jaske
hypothesizes molecular diffusion of heat while
the WRE model employs a turbulent transfer
mechani sm), and
4. the numerical solution scheme.
In addition, Jaske's model, in its present state of development,
cannot handle a reservoir with fluctuating water surfaces whereas
this poses no problem in WRE's model.
- 110 -
-------
In the chapters which follow, the details of segmenting
the reservoir and of interfacing adjacent segments with respect to
heat and mass transfer are given. The prototype application of
the model is demonstrated by a simulation of Lake Roosevelt for
the year 1967.
- Ill -
-------
IX. SEGMENTING THE RESERVOIR
NOMENCLATURE
In the following text, reference will often be made to
the words segment and element. Segment is used to denote a
certain reach in the reservoir. Since the heat transfer character-
istics of each segment will be represented by the deep reservoir
model, it is necessary to divide each segment vertically into a number
of elements. These elements are conceptually identical to the elements
of the deep reservoir model. Figure 33 illustrates schematically a
reservoir divided into segments and elements.
SEGMENT,
ELEMENT
FIGURE 33 Schematic Illustration of a Reservoir Divided
Into Segments and Elements
- 113 -
-------
PROCEDURE
The purpose of segmenting the weakly-stratified reservoir
is to break it into a set of units each of which can be considered
to behave essentially like a deep reservoir, i.e., the isotherms
within the segment can be considered to be approximately horizontal.
While the decisions regarding the number of segments into which the
reservoir is to be divided and the location of the physical boundaries
of each segment are primarily subjective, there are at least four
criteria upon which these decisions should be based. These criteria
are:
I. reservoir geometry;
2. residence time within the segments;
3. points of tributary inflow; and
4. desired longitudinal resolution of the isotherms.
Proper attention to these factors will eliminate much of the guesswork
from the segmenting process.
GEOMETRICAL CONSIDERATIONS
It is important for two reasons that a suitable number of
reservoir cross sections be available for examination. First, the
location of the endpoints of the reservoir segments depends upon the
channel geometry; and, secondly, these cross sections will be used later
to determine the volume-stage relationship for each segment. Ideally,
the sections should be taken at points in the reservoir where there
are significant changes in channel geometry, e.g., at narrow points,
at wide points, and in bends.
- 114 -
-------
The best modeling will result if the ends of the segments
are in straight sections of the channel and are of small cross
section since the horizontal flow will be defined best through
these sections. The broader sections lying within the segment
will be characterized by lower horizontal velocities and are more
suitable for representation by the deep reservoir model.
RESIDENCE
The residence time of the reservoir flow within a segment
(the volume displacement time) is a very important consideration
because, if the residence time is less than the simulation time
step, a violation of the necessary conditions for a valid heat budget
for the segment can result. Consequently, some estimate must be
made of the shortest expected residence time within each segment.
In general, this event will occur during the peak discharge period.
This period of high flows usually occurs before the onset of the
stratification period; hence, it can be assumed that the reservoir
is in motion over its entire depth and the residence time can be
computed as the ratio of the segment volume to the maximum expected
daily value of discharge. It should be noted however, that during
the stratification period the flow does not generally extend over
the entire reservoir depth, but rather is confined to some fraction
of the depth. Thus, while the total flow may not be large, if it is
sufficiently concentrated within some vertical range of depths, the
resulting velocities could be large enough to cause the flow to
travel completely through a segment during a simulation time step,
- 115 -
-------
even though the residence time is sufficient for the peak flow period.
If this situation is encountered during the simulation runs, the
segment lengths can be extended or the simulation time step shortened.
In most cases, the latter procedure is the most expedient remedy.
TRIBUTARY INFLOWS
Consideration of the points of tributary inflow in segmenting
the model is not critical; however, some consideration of them in the
initial segmenting of the model should produce better results. At
certain times of the year, an individual tributary, as a result of
its temperature or density, may create an interflow in the reservoir
at an elevation independent of that of the main body of the flow.
This phenomenon is best handled in the model if the tributary is
introduced at about the middle of the reach because if it does not
enter the reservoir in the region of the mainstream interflow, it may
very well experience an internal lateral spreading in the upstream
direction as well as downstream.
LONGITUDINAL RESOLUTION OF ISOTHERMS
A fourth consideration in segmenting the reservoir is the
resolution that is desired in describing the longitudinal variation
in the elevation of the isotherms. Since the deep reservoir model
is applied to each segment, there results one temperature at each
elevation for each segment. Consequently, the more segments there
are in the reservoir, the more resolution is obtained in determining
- 116 -
-------
the position of the isotherms. On the other hand, if the segments
are too small, the representation of the segment by the deep reser-
voir model becomes less accurate; hence, resolution within the seg-
ment decreases. The choice of segment size at this point is subjective,
however, on the basis of experience with the deep reservoir model,
it is felt that the best longitudinal resolution can be obtained,
consistent with adequate internal resolution within the segment if
the minimum residence times are set up as about one day. This, of
course, implies a simulation time step of about one day, the interval
that was previously recommended in Part I.
To summarize the segmenting procedure, the most important
physical considerations are:
I. placement of the segment ends in narrow straight
portions of the channel; and
2. assurance that the residence time of any segment
does not exceed the simulation time step interval.
In addition to the guidelines prescribed above, however, technical
sagacity is an important asset to the segmenting process, since many
of the required decisions are subjective.
- 117 -
-------
INTERFACING ADJACENT RESERVOIR SEGMENTS
Interfacing as discussed in this chapter is the process
of transferring heat and mass between adjacent reservoir segments.
In contrast to the reservoir segmenting process, the interfacing
is executed internally in the weakly stratified reservoir model
and thus, requires no action on the part of the individual. The
discussion which follows describes the method of interfacing as
programmed into the model.
Since the deep reservoir model, as applied to an individual
segment, is one-dimensional along the vertical axis, a horizontal flow
distribution is not explicitly defined within the segments. Rather,
it is implicitly defined by specifying the flow distribution across
the upstream and downstream interfaces. Consequently, it is the inter-
facing criteria which determines the horizontal flow pattern through
the reservoir.
In the absence of a fundamental description of unsteady
flow in a thermally stratified impoundment, it becomes necessary to
surmise the true nature of the flow on the basis of existing results
from steady state theory and experiment. Through such conjecture,
it is possible to arrive at several plausible schemes for describing
the horizontal flow distribution within-the reservoir. The method
described here is one of these schemes. While it is by no means a
complete answer to the problem, it is believed that the method does
provide an adequate description of the horizontal flow through a
weakly stratified reservoir.
- 119 -
-------
It has been tacitly assumed that the horizontal flow
through the reservoir may be characterized as an interflow, i.e.,
a flow that is vertically confined between two elevations in the
reservoir. The difference in elevations which defines the thickness
of the interflow is generally less than the reservoir's depth. This
assumption is based upon work done by Jaske (26) on Lake Roosevelt
and upon theoretical and experimental investigations as reported by
Long (27) and Yih (22). Since the location of the interflow and
its thickness are governed by momentum forces, buoyancy forces, and
the degree of stratification, the interfacing criteria must also
reflect these factors. In addition, heat and mass continuity must
also be maintained across the interface.
MOMENTUM AND BUOYANCY CONSIDERATIONS
The location of the interflow on the vertical axis of the
reservoir is governed jointly by the momentum of the flow and by
buoyancy forces. The momentum force tends to propel the fluid along
horizontal lines or planes; however, as the interflow travels through
the reservoir on this plane, it will experience vertical buoyancy
force as a result of the tilted isotherms. The net path which is as-
sumed by the interflow will be determined by the resultant of these two
forces.
In an effort to take the buoyancy and momentum forces into
account in the reservoir model, it has been assumed that the interflow
within a reservoir segment is contained between two horizontal planes,
thereby conserving momentum within the segment. At the interface
- 120 -
-------
between adjacent segments, however, there is an abrupt change in the
elevation of the interflow. Upon crossing the interface of the up-
stream segment, it is assumed that the interflow centers itself about
the elevation z' in the downstream segment where the temperature is
equal to the arithmetic mean temperature of the outflow from the up-
stream segment. The change in elevation at the interface is a step-
wise attempt to account for vertical movement of the interflow as a
result of buoyancy forces.
INTERFLOW THICKNESS
The thickness of the interflow is known to be functionally
dependent upon the discharge and the degree of stratification; hence,
it is a function of the densimetric Froude number. In the absence of
a better definition, Debler's criteria as discussed in Chapter V of
Part I has been used to evaluate this thickness. Using equation 21,
the interflow thickness D can be evaluated as
D
= 2d = 2(2.04) •
1/2
(27-A)
or
D = 2.88
w /gE
1/2
(27-B)
where Q is the total discharge across the upstream interface of the
segment, and w and E are the reservoir width and stability, respectively,
at the elevation z1 as defined in the preceding paragraph. With the
flow spatially centered at z1, the upper and lower elevations bounding
- 121 -
-------
the flow are in general, z' +^D, and z1 - T>-D, respectively. Excep-
tions occur when z1 + i D is higher than the water surface, or when
z' - i D falls below the reservoir bottom. In these cases, D remains
the same but the interflow layer is shifted up or down as required to
contain it within the vertical reservoir boundaries.
MASS CONTINUITY
Mass continuity within the reservoir is maintained as follows:
a mass balance is first performed on the reservoir as a whole and the
water surface elevations are determined. It is assumed that the water
surface is horizontal throughout the reservoir. Given the water sur-
face elevation for each time step, it is then possible to do a mass
balance around each segment beginning with the segment farthest upstream
and working downstream toward the dam. Given the inflow to the farthest
upstream segment as the mainstem river flow, the storage equation is
applied and the outflow from that segment is computed. This outflow
becomes the inflow to the next segment. The solution is repeated
for successive segments until the last segment is reached. Segments
with tributaries must include the tributary inflow as part of the total
inflow. If the overall mass balance has been done correctly, the out-
flow computed for the segment abutting the dam will be identically
equal to the reservoir outflow.
- 122
-------
HEAT CONTINUITY
Since the thickness of the interflow generally will be
different on each side of the interface, the temperature distri-
butions will also be different. Approximating the temperature dis-
tribution within the interflow as a linear function of the elevation,
the temperature distribution of the inflow into the downstream seg-
ment can be given as
T.1 = m(z.' - zo') + b (28)
where the primes indicate properties associated with the downstream
segment, z ' is the elevation at the bottom of the interflow zone,
and b is defined as T , the temperature of the bottom element within
the interflow zone in the upstream segment.
The value of m must be determined in such a manner that heat
continuity is maintained across adjacent segments, i.e.,
E pc Q ' T ' = I DC Qn. T (29)
n n i n uj j
where QT.' is the 'inflow to the i^" element of the downstream segment,
th
CL. is the outflow from the j segment in the upstream element, and
n is the number of elements in the interflow zone. As before, the primes
indicate properties associated with the downstream segment. If Q is
the total discharge across the interface of the adjacent segments, the
velocity, u, which has been assumed constant over the interflow zone,
- 123 -
-------
can be defined as
u = £ = .SL = , and (30-A)
where L is the length of the segment, A is the mean cross-sectional
area of the interflow zone within a segment, V is the volume of the
segment occupied by the interflow zone, and V. is volume of the i
element in the interflow zone. Solving the right-hand equality of
equation 30 for Q~. and Qj.' gives
V. .
Qn- = Q -fr xnA (31-A)
Uj V j anu
Qn' = Q 4r (31-B)
Finally substituting equations 28 and 31 into equation 29, there results
the following definition of the slope m:
m ~~ (32)
Substituting equation 32 for m in equation 30 and taking the value of
b as T gives a temperature distribution to the interflow entering the
upstream face of a reservoir segment that assures the maintenance
of heat continuity across the segment interfaces.
- 124 -
-------
The interfacing criteria as described above is summarized
schematically in Figure 34-A for a reservoir with three segments.
The assumed relationship between the properties of the segmented
model and the prototype is illustrated in Figure 34-B. Note that
the interflow zone in Segment 3 is centered vertically at the ele-
vation where the temperature is equal to that of the mainstem.
Note also that the outflow from Segment 1 is governed by the selective
withdrawal criteria as described in Chapter V of Part I.
Tributary inflows to the segments are accounted for in
the heat budget in the following way. If the elevation at which
the inflow enters the segment falls within the interflow zone, the
flow is added' to the total flow at 'that elevation. If the elevation
falls outside the region of interflow, the tributary is then considered
as an independent interflow and handled in the same manner as
described above. If,at any segment,the two interflow layers coincide,
the flows are then combined at that segment, and treated as a single
interflow through the remainder of the reservoir.
- 125 -
-------
A. MODEL REPRESENTATION
rv>
CTl
SEGMENT I
SEGMENT
SEGMENT.
T =
/, fQ = Mainstream
1 R River Discharge
" = Mainstream
River Temp.
| = ARITHMETIC MEAN
TEMP OF THE OUTFLOW
FROM SEGMENT j
-J,.
B. PROTOTYPE
Figure 34. Relationship of the Interflow Properties of the Model to
Those of the Prototype Reservoir
-------
XI APPLICATION TO LAKE ROOSEVELT
The weakly-stratified reservoir model, as developed in
the preceding three chapters, was based on the knowledge that presently
exists with regard to the thermal behavior of such water bodies. The
test of its validity required that an appropriate prototype reservoir
be chosen and simulated. Lake Roosevelt, behind Grand Coulee Dam,
was chosen as the test case for the following reasons: 1) the amount
of data and its accuracy are the best available; and 2) since this
reservoir is the most important weakly-stratified reservoir in the
Columbia River System, it is imperative that the FWPCA Study Team
has a working model of this reservoir as a part of its Columbia
River Thermal Effects Study. The details and results of this
simulation are given below.
GENERAL INFORMATION
Lake Roosevelt is located in the mainstem of the Columbia
River in Northeastern Washington (see Figure 1). At full pool, this
reservoir backs up the Columbia River from Grand Coulee Dam to the
- 127 -
-------
Canadian Border, a distance of nearly 150 miles. The impoundment,
which has a usable storage capacity of 5.23 million acre-feet and
a total capacity of 9.56 million acre-feet, is used for power
development and irrigation. Irrigation water is diverted from the
reservoir at Grand Coulee Dam and pumped through a lift of about
280 feet for a distance of two miles into Banks Lake. From there
it is distributed through a system of canals to the Bureau of
Reclamation's Columbia Basin Project.
As was the case for Hungry Horse Reservoir, time and money
constraints required that the period of simulation be limited to
a single year. The year 1967 was recommended by the FWPCA Study
Team as the period for which the data base was most reliable. The
simulation period was further reduced to a length of 76 days since
thermal profiles were taken in the reservoir only between 11 July
and 25 September of that year. This period is considered adequate,
however, since it covers the time from the onset of stratification
within the reservoir to a time that is well into the cooling phase
of the annual thermal cycle.
PHYSICAL AND THERMAL DATA
A plan view of Lake Roosevelt is shown in Figure 35. From
this F.igure, it is seen that there were-some 92 cross-sections avail-
able in the reservoir proper. These sections, which were composited
- 128 -
-------
COULEE
DAM
C TEMPERATURE STATION
FRANKLIN D. ROOSEVELT LAKE
Figure 35. Plan View of Lake Roosevelt
- 129 -
-------
by the FWPCA Study Team from a set of sounding maps, were taken
at each reservoir bend and at each point of significant change in
channel geometry.
Gaged inflows'to the reservoir are the Columbia River
mainstem, Kettle River, Colville River (entering at Kettle Falls),
and the Spokane River. Ungaged tributaries contribute less than
0.05 percent to the total flow. Additional physical and thermal data
were summarized in Table 3. Unfortunately, there were no tempera-
ture data on the reservoir outflows; consequently, it was not possible
to compare the computed outflow temperatures to observed values.
METEOROLOGICAL DATA
In the absence of meteorological observations at Lake
Roosevelt, it was necessary to transpose the data from Spokane,
Washington. The use of these point values in the model was equi-
valent to assuming uniform meteorological conditions to exist over
the entire length of the reservoir. While neither the transposition
nor the assumption of uniform meteorology is untenable for average
daily values of solar radiation and cloud cover, there is some
question as to the validity of transposing air temperatures, humidity
data, and windspeeds from Spokane to Lake Roosevelt let alone assuming
them to be uniformly constant over the entire 150 mile reservoir length,
Moreover, the windspeeds taken at Spokane were measured at a distance
of 100 feet above the ground. Consequently, it was necessary to
- 130 -
-------
TABLE 3
PHYSICAL AND THERMAL INFORMATION
Lake Roosevelt
PHYSICAL DATA
Location (Grand Coulee Dam): Latitude 47°57', Longitude 118°59'
Elevation:
393 meters (full pool elevation)
Intake Elevations:
Outlet Tubes
Penstock Intakes
Diversion Intake
285 meters
316 meters
346 meters
318 meters
364 meters
Spillway
Type - 504 meter crest with drum gates;
located in middle half of darn.
Elevation of Crest 384 meters
Mean Gaged Inflow:
Columbia River
Kettle River
Colville River
Spokane River
2785 CMS
82 CMS
8 CMS
230 CMS
Average Reservoir Discharge _
to Volume Ratio: 8.30 yr
SURFACE AREA, rn*x !0~8
o.o
£400
1.0
2.0
CD
35°
< 300
UJ
LLJ 250
SURFA
RESERVOIR VOLUME
. * BOTTOM OF RESERVOIR
L_J_,J__.1__1
3.0
0.0 2.0 4.0 6.0 8.0 10.0 12.0
RESERVOIR VOLUME, m3x I0~9
- 131 -
-------
TABLE 3
(Continued)
TEMPERATURE DATA
Reservoir:
Temperature profiles for each of the
ten stations shown in Figure 35 at
about 7 day intervals.
Tributaries
Columbia River, mean daily temperature
at International Boundary.
Kettle River, synthesized mean daily
temperature from FWPCA.
Colville River, assumedat same
temperature as Kettle River.
Spokane River, mean daily temperature
of the river at Long Lake about 30
miles upstream from Lake Roosevelt.
- 132 -
-------
adjust the evaporation coefficient to a value such that the computed
evaporation rates for the reservoir were equal to the average values
of evaporation observed in the areas surrounding Lake Roosevelt. This
adjusted value was taken as 1.3 x 10 mps/mb/mps. The solar radia-
tion extinction coefficient within the reservoir was estimated at
0.691 nf which gives 99.9 percent extinction at a depth of 10 meters.
The extinction depth was based on secchi disk measurements of 10
meters that were taken in the reservoir.
SEGMENTING THE RESERVOIR
Following the guidelines' established in Chapter IX,
Lake Roosevelt was divided into six segments as illustrated in
Figures 36 and 37. The downstream faces of the segments were
located in narrow and, as often as was practical, in straight sec-
tions of the reservoir channel. The compromise in placing the
interfaces in straight portions of the channel was necessitated by
residence time considerations.
Residence times for each segment are given in Table 4.
It is observed that the minimum residence time in Segment 6 is
less than one day, which is not particularly desirable, since it
is less than the recommended simulation time step. There are two
solutions to this situation as previously stated in Chapter IX:
1) increase the reach length; or 2) decrease the simulation time
step. In this instance, however, neither of these steps was taken
for the following reasons. First, as seen in Figure 36, Segment 6
- 133 -
-------
SEGMENT 2 SEGMENT
Figure 36. Plan View of Lake Roosevelt as-Segmented for Simulation
Model
- 134 -
-------
oo
i
400 r
SPILLWAY
DIVERSION
INTAKE
w
L.
o
£ TURBINE Tzpc L
£ INTAKE ^r^2_J
O
OUTLET
TUBE
^ OUTLET
5 TUBE 97,
_i 275
LU
FROM DAM, KILOMETERS
Figure 37. Reservoir Segments in Profile View, Lake Roosevelt
-------
is already quite long (- 27 miles); to Increase its volume sufficiently
to give a minimum residence time of one day would require the reach
length to be extended another 50 percent. Such a reach would be too
long to be realistically modeled. Second, the peak discharge period
occurs about the middle of June which is prior to the onset of strati-
fication. Hence, the small energy budget violations which occur
during this period are not serious. By the end of June when stratifi-
cation begins to set in, the reservoir discharge has decreased to a
value such that the residence times are greater than one day.
TABLE 4
RESIDENCE TIMES FOR INDIVIDUAL SEGMENTS
IN LAKE ROOSEVELT
SEGMENT
average
1 7.4
2 8.1
3 5.0
4 8.7
5 7.6
6 3.2
RESIDENCE TIME,
maximum
26.0
28.6
17.6
30.5
26.8
11.1
days
minimum
1.9
2.1
1.3
2.3
2.0
0.8
Since the reservoir above station 57 is essentially completely
mixed, it was not explicitly included in the simulation. Rather, the
temperature changes through this upper reach were determined by the
- 136
-------
FWPCA Study Team through the application of their river-run model
between the Canadian border and station 57. The temperatures
thus predicted at station 57 were used as the input temperatures
to Segment 6.
Referring again to Figure 36, note that the Sanpoil and
Spokane arms of the reservoir have not been included as part of
their appropriate segments. The reason for excluding these arms
is that the model simply includes their properties as being a part
of the properties of the main channel. This, of course, alters the
characteristics of the interflow through these segments. If it is
felt necessary to include such arms in the simulation, they should
be considered as separate reservoir segments and should be interfaced
to their appropriate segments in the main channel. Such a refine-
ment was not considered necessary in this simulation since previous
studies on the reservoir have indicated that these arms do not signi-
ficantly affect the thermal structure of the main body of the reser-
voir. Consequently, the Sanpoil arm is completely unaccounted for in
the simulation. The Spokane arm is also neglected; however, the
Spokane River is introduced as a tributary to Segment 3 using the
temperatures and discharges recorded at Long Lake.
- 137 -
-------
SIMULATION RESULTS
For each of the reservoir segments shown in Figures 36 and
37, the volume-stage and area-stage relationships were determined from
cross-section data. The thermal regime of the Lake Roosevelt was
then simulated over the 77 day period from 11 July through 25 September
1967. As with the deep reservoir model, the performance of the weakly-
stratified reservoir model was judged on the basis of its ability to
reproduce the observed temperature profiles in the reservoir. Typical
results for four days during the simulation period are illustrated in
Figure 38 which shows the computed and observed profiles for e-ach
reservoir segment, the elevation of the isotherms along the longitudinal
axis of the reservoir, and the position of the flowing layer within the
reservoir. Note that there were two observed profiles in Segments 2,
4, and 5.
GENERAL OBSERVATIONS
A comparison of the observed and computed temperature profiles
in Figure 38 indicates that the model simulates, quite well, the
general progression of the thermal events in Lake Roosevelt, and that
it does, in fact, produce tilted isotherms. The only outstanding
discrepancy between the simulated results and those observed in the
prototype is that the reservoir tends to cool too rapidly toward the
end of the simulation period (see Figure 38-D). Even so, over the 77
day simulation period, the maximum deviation of the computed temperature
from those observed was 2.6 °C; the average deviation for the period
was less than 1.3 °C.
- 138 -
-------
400
3682 CMS
16.t °C "*
271 CMS
15.8 °C
1986 CMS ,325
• 14.9 °C
O
1-
LJ
_J
LL!
300 I-
275 iM^^tt3®m
c-{° ^^&^w^^^
^uz^^mzM^Az.
0 25
200
400 r
SEGMENT I
50 75 !00 125 150
DISTANCE FROM DAM, KILOMETERS
SEGMENT 2 SEGMENT J SEGMENT 4 SEGMENT 5 SEGMENT 6
/
14 16 18
14 16 18 20
i i i i i
,13 15 17 19
J!3 i5 17 19 2!
13 !5 17 19 2!
LL!
13 15 17 19 2!
TEMPERATURE , °C
OBSERVED
COMPUTED
JULY 24,1967
Figure 38-A. Typical Results from Simulation of the Thermal Regime
of Lake Roosevelt, 1967
-------
o
I
400
790 CMS
22.4 °C
SEGMENT I
SEGMENT 2
75 !00 125 150 !75 200
DISTANCE FROiYi DAM, KILOMETERS
SEGMENT 3 SEGMENT 4 SEGMENT 5 SEGMENT 6
400 r
CO
17 19 2!
15 17 19 21
14 16 18 20
I I I I I JI4 16 18 20 22
14 16 18 20 22
14 IS 18 20 22
TEMPERATURE, °C
OBSERVED
COMPUTED
AUGUST 14 ,1967
Figure 30-B. Typical Results from Simulation of the Thermal Regime
of Lake Roosevelt, 1967
-------
400
1
315 CMS
19.7 °C
CD
', 2282 CMS
' 17.9 °C
o
K
LU
LU
SEGMENT I
25 50
SEGMENT 2
150
!75
200
DISTANCE FROM DAM, KILOMETERS
SEGMENT 3 SEGMENT 4 SEGMENT 5 SEGMENT 6
400
15 17 19
j //\ i i i
15 17 19 21
19
J_JLJ_J_LJI4 16 18 20 22
14 16 18 20 22
275
15 17 19 21 23
TEMPERATURE, °C
OBSERVED
COMPUTED
AUGUST 29,1967
Figure 38-C. Typical Results from Sfmulati'on of the Thermal Regime
of Lake Roosevelt, 1967
-------
ro
!
400 r-
,-r
118.7 CMS
17.3 °C
2103 CMSJ
_~ 17.3 °C
O
LU
_]
LU
0 25
50
SEGMENT I
SEGMENT 2
75
DISTANCE FROM DAM, KILOMETERS
SEGMENT 3 SEGMENT 4 SEGMENT 5 SEGMENT 6
400
375
£ 350
cu
£
z~ 325
o
< 300
LU
_J
LU OTP;
r
/
/
/
/
/
1
1
1
1
1 I 1 1 1
-
-
4 16
/
ii
ii
!'
n
i/
n
1
\\ ! 1 I II
18 20 22
-
I 1 1
/
/
/
/
1 1 I ! 1
, J
/ / L // L"
/ / if ill
If \ 14 16
// L ! I ,' 1 1 1 1
// 14 16 18 20
' i i i i
V
1 J
18
4 16 18 20
4 16 18 20 22
OBSERVED
COMPUTED
14 16 18 20 22
TEMPERATURE, °C
SEPTEMBER 25 ,1967
Figure 38-D, Typical Results from Simulation of the Thermal Regime
of Lake Roosevelt? 1967
-------
The behavior of the reservoir interflow during the simulation
period is interesting. Referring again to Figure 38, it is observed
that during July, when the temperature gradients are small, the flowing
layer extends over most of the reservoir depth. As the intensity of
the reservoir stratification increases, the bottom of the flowing
layer moves progressively higher in the reservoir until it assumes
the position indicated on 14 August. Thereafter the reservoir begins
to cool and the bottom of the flowing layer begins to move downward
as observed on 29 August. Finally in late September the intensity of
the stratification has become sufficiently weak so that the flowing
layer once more extends over the total reservoir depth. This is a
rather important phenomenon because it indicates that during the period
when the reservoir surface layers are the warmest, this surface water
is being drawn off and discharged through the outlets. Consequently,
the release temperatures from the reservoir would be expected to be
significantly warmer than the mean temperature in a horizontal plane
at the elevation of the reservoir outlet. According to Jaske (26)
this phenomenon was, in fact, observed during the 1963 operation of
the Columbia River Cooling Program undertaken by the Hanford Operation.
Finally, it is noted that the Spokane River (which, except
for the mainstem, is the largest tributary inflow to the reservoir)
did not appear as a separate interflow at any time during the simulati
period. Segment 3 was purposely made small in order to account for
this phenomenon should it occur; however, the river temperature was
such that the tributary always entered the reservoir at an elevation
contained within the primary interflow zone.
on
- 143 -
-------
SPECIFIC OBSERVATIONS
In comparing the computed and observed profiles, one should
be aware that the computed profiles are constructed from the average
daily values of temperatures on horizontal planes within the segment
while the observed profiles are composed of instantaneous point values.
As a result of the unsteady nature of the reservoir on an hour to hour
basis (outlet operation, internal waves, lateral temperature variations,
diurnal surface heating and cooling, etc.), the instantaneous observed
profiles may be expected to deviate from their average daily values
by 0.5 °C in the reservoir depths and by 1 to 2 °C at the surface.
While these deviations are not an important consideration when com-
paring profiles at some distance below the water surface, they are
quite significant when comparing surface temperatures. Since the
reservoir profile measurements were taken during the day, it is reason-
able to expect that the observed surface temperatures are somewhat
higher than the average daily surface temperature. Allowing for this
amount of variability in the observed surface temperatures5 it is
surprising how closely the computed and observed values compare, es-
pecially when it is recalled that point meteorological data, measured
at Spokane, was assumed to apply over the entire reservoir.
The examination of the computed and observed profiles for
the several reservoir segments revealed that there was often a
temperature difference of 1.5 to 2.0 °C in the profiles near the
bottom of Segments 1, 2, and 3. Several examples of this behavior
can be seen in Figure 38. It is now thought that this deviation re-
sulted from neglecting the effect of the Spokane arm of the reservoir
in the model simulation. As the reservoir interflow warms, there
is a significant quantity of cold water stored in this arm that flows
- 144 -
-------
out through the reservoir mainstem and is replaced by the warmer water
of the interflow. Neglecting to account for this cold water outflow
is equivalent to short circuiting the main reservoir interflow through
Segment 3. Consequently, the simulated rate of heating downstream is
greater than that observed in the prototype. Unfortunately, there
was insufficient time to thoroughly investigate this problem; however,
some studies of a very preliminary nature indicated that if the volume
properties of the Spokane arm are lumped together with those of Seg-
ment 3, the simulation of the lower portions of the temperature pro-
files in Segments 1, 2 and 3 would be improved.
Another feature of the simulation which requires discussion
is the accelerated cooling in the fall of the year. This behavior,
which can be observed in Figure 38-D, actually starts sometime between
5 and 12 September, and is thought to occur for one or more of the
following reasons: 1) the simulated rate of heat loss through the
air-water interface is too large; 2) the present method of describing
the velocity and temperature distribution of the interflow is
incorrect; or 3) the storage effect of the Spokane arm is not taken
into account. Again time restricted the amount of effort which could
be devoted to this problem; however, on the basis of one of the test
runs, some insight was gained on possible ways to decrease the simulated
rate of cooling.
In order to acquire a better understanding of why the simulated
rate of reservoir cooling was greater than that observed in the proto-
type, a test run was executed in which the convective mixing (mixing
that results when negative temperature gradients form in the reservoir)
- 145 -
-------
was suppressed. Figure 39 shows a typical shape of the resulting
profiles. It is readily apparent that heat is being lost through the
air-water interface; hence, it was most logical to suspect that this
simulated heat loss rate was too large; i.e., the meteorological data
was incorrect. However, a check of the reliability of the data by the
FWPCA Study Team resulted in nothing suspicious. Moreover, the success
with which the surface temperatures were modeled during the warming
period inspired confidence in the reliability of the data during the
cooling period. Thus, the possibility of bad meteorological data
was discarded as a consideration.
400
390
.380
Q>
1 370
-
2 360
r~~
>
y 350
LJ
340
•a^n
- OBSERVED PROFILES-,
"N. 4
\ i j
COMPUTED S\\\
PROFILE \l|
yi
|i
)
i
i
— 1 ;•
W
ft
'!*
'
ri ••
n '
i /
i i i
x-
SEGMENT 5
Sepf. 12, 1967
i
15 16 17 18 19
TEMPERATURE,°C
FIGURE 39. Typical Shape of Computed Temperature Pro-
file during Cooling Period when Convective
Mixing is Suppressed
- 146 -
-------
A second, and more probable, explanation of the cooling
phenomena is the manner in which the temperature is distributed
through the interflow as it crosses the segment interfaces. It was
assumed that this distribution is linear. However, if it is not
linear, the effect of the linearization is to displace heat from the
reservoir surface to a lower depth as seen in Figure'40. The net
result is to induce aritificial surface cooling. Rather than linearize
the temperature profile the actual profile should be used; however, as
noted in Chapter 'X, this is a more difficult task and requires additional
bookkeeping which was not required in the linearization scheme.
o
h-
>
LiJ
_J
LJ
ACTUAL
DISTRIBUTION
DISTRIBUTION
TEMPERATURE
FIGURE 40. Effect of Linearization on the
Temperature Distribution of the
Interflow
- 147 -
-------
A third possible explanation for the observed cooling
involves consideration of the horizontal velocity distribution within
the interflow. As presently developed, it is assumed that the
velocity distribution is uniform in the interflow zone and zero out-
side of this zone. If the simulated reservoir velocities are larger
near the surface than in the prototype, an excessive amount of
surface heat will be carried out of the reservoir by longitudinal
advection. This additional heat loss, coupled with the loss of surface
heat to the atmosphere, would produce an accelerated cooling rate. A
recent report by Jaske (25) contains actual measurements of the velocity
profile in Lake Roosevelt during 1967. These measurements indicate
that the horizontal velocity in the reservoir near the surface is
significantly lower than that observed deeper in the interflow zone.
The incorporation of such a velocity distribution into the model would
significantly decrease the cooling rate. Unfortunately, this data was
not available at the time when the reservoir model was being developed;
consequently, the uniform velocity distribution was assumed due to
the lack of a more suitable definition.
Finally, the short-circuiting through Segment 3, which results
from neglecting the Spokane arm, is thought to contribute to the
accelerated cooling in Segments 1, 2 and 3. The effect of neglecting
this arm in the fall is just the opposite of that which results from
neglecting it in the spring. As the Columbia River cools, there is a
significant amount of warmer water in the Spokane arm that is replaced
by the cooler Columbia River flow and flows out through the reservoir
mainstem. Neglecting this storage effect causes the cooler water to
arrive at the lower reaches of the reservoir too early in the fall,
- 148 -
-------
thereby causing Segments 1, 2, and 3 to cool more rapidly in the model
than in the prototype. Unfortunately, there was insufficient time to
investigate this phenomena, even in a preliminary way.
REMARKS
On the basis of the results from the 1967 test simulation,
it is believed that, on the average, the weakly-stratified reservoir
model can simulate the temperature regime of Lake Roosevelt within
1.5 °C of the true values; however, deviations as high as 2.6 °C
might be expected during fall cooling. As such, the model is considered
to be adequate for the types of analyses contemplated in the FWPCA
Columbia River Thermal Effects Study except possibly during the period
of cooling when its adequacy is judged as marginal.
If the accuracy of the model during the fall cooling portion
of annual temperature cycle is judged acceptable, the model can be used
for simulating the annual thermal cycle of the reservoir within the
confidence limits specified above. Moreover, several years may be simulated
in succession without concern over the possibility of carryover error
in the heat budget from year to year since the Columbia River discharge
is sufficient to completely displace the old reservoir water between
the end of the fall cooling cycle and the onset of warming the follow-
ing spring.
If efforts are made to improve the model during the cooling
period, it is suggested that the description of the temperature profile
of the interflow at the interface between adjacent segments be refined
- 149 -
-------
first, since of all the probable factors contributing to the acceler-
ated cooling in the model, it is thought that linearizing the temperature
profile of the interflow is probably the most significant. Secondly,
the effects of the Spokane arm should be investigated in more detail.
It is thought that the effects of this arm can be adequately
represented if its volume-stage characteristics are simply lumped
together with those of Segment 3. Finally, the velocity distribution
of the interflow is considered to have the least influence on the rate
of cooling. Consequently, it is suggested that this refinement be given
the lowest priority in further development of the model.
- 150 -
-------
XII SUMMARY AND CONCLUSIONS -- PART II
GENERAL
The weakly-stratified reservoir model as developed and
tested in the preceding four chapters was based on the assumption
that a weakly-stratified reservoir can be visualized as a series of
deep reservoirs. Adjacent "reservoirs" are coupled to one another by
the requirement for conservation of heat and mass in the reservoir as
a whole. The longitudinal flow in the reservoir was considered to be
a density current the thickness of which was determined by Debler's
criteria. The vertical placement of the interflow was established
on the basis of its mean temperature and the elevation at which that
temperature existed in the reservoir.
A test of the model on Lake Roosevelt over a 77 day period
from 7 July through 25 September 1967 showed that, in general, the
model simulated with reasonable accuracy the thermal regime of the
reservoir; however, some difficulties were encountered in simulating
the reservoir during the period of cooling in September.
- 151 -
-------
On the basis of the test case related in Chapter XI, it
is concluded that the weakly stratified reservoir model, as developed
herein, is adequate for the types of analyses contemplated in the
FWPCA Columbia River Thermal Effects Study, except during the period
of cooling in the fall of the year. During this period, the adequacy
of the model is considered marginal. Preliminary studies indicate,
however, that a refinement of the temperature distribution in the
interflow zone should significantly improve the model's capability
to simulate the fall cooling period.
On the basis of evidence gathered through the test runs on
Lake Roosevelt, it was concluded that the Spokane arm of the reservoir
does affect the thermal regime of the reservoir and that a better
representation of the temperature profiles near the reservoir bottom
in Segments 1, 2, and 3 could be obtained if the storage effect of
the Spokane arm was accounted for in the model. In addition, inclusion
of the storage effect should improve the simulation during fall cooling
although the extent of the improvement is not known.
RECOMMENDATIONS
At this juncture, it might be well to point out that the con-
clusions drawn concerning the applicability and limitations of the
weakly-stratified reservoir model have been derived from observation
of the results of a simulation of Lake Roosevelt over a single 77 day
period. While it is believed that these observations and the resulting
- 152 -
-------
conclusions are basically valid, it is strongly recommended that the
first priority in FWPCA's use of the model be the simulation of Lake
Roosevelt for different years (preferably two or three) during which
internal temperature measurements are available. These simulations
would serve to establish additional confidence in the model and to
confirm its weaknesses as pointed out above.
Assuming the weaknesses noted in Chapter XI are confirmed
and that it is desired to improve upon the accuracy of the model in
its present form, the following steps are recommended:
I. Include the Spokane arm in the reservoir by simply
lumping its volume-stage characteristics with those
of Segment 3. This step is recommended first since it
is the simplest to do. If simple inclusion of the volume-
stage relationship improves the simulation in Segments
I, 2, and 3, it should be retained. If not, it should
be dropped and further investigation deferred until
after the refinement of the temperature distribution
of the interflow.
2. The description of the temperature distribution of
the interflow as it enters the upstream face of a seg-
ment should be reformulated to reflect more precisely
its actual distribution. It is believed that this
refinement would, at least, make the accuracy of the
model during fall cooling consistent with that over
the remainder of the year; in addition, it is also expected
that the simulation would be improved over the entire
year.
- 153
-------
3. The velocity distribution of the interflow zone should
be reformulated on the basis of the information that
has recently become available on Lake Roosevelt (25).
This refinement, however, is presently rated as a low
priority since it is not thought that a better description
of velocity distribution would significantly improve the
mode ITs accuracy.
CAVEAT
Since the weakly-stratified reservoir model was developed
from general concepts and theory, it should be applicable to any reser-
voir falling within the weakly-stratified reservoir class as defined
in Chapter II. However, the only test which has been made of its valid-
ity is the simulation of Lake Roosevelt as reported herein. Hence, its
general applicability remains to be proven through the simulation of sev-
eral other reservoirs of the weakly-stratified class. Subsequent to the
establishment of its general applicability the model could then be used
with confidence as a planning and analysis tool for proposed reservoirs
and for reservoirs having insufficient thermal data for verification
purposes.
- 154 -
-------
XIII LIST OF REFERENCES
1. "Water Temperature; Influences, Effects and Control," Proc.,
12th Pacific Northwest Symposium on Water Pollution
Control, U.S.P.H.S., Corvallis, Oregon, 7 November 1963.
2. Klein, L., River Pollution II, Causes and Effects, Butterworths,
London, 1962.
3. Hutchinson, E. G., A Treatise on Limnology, Vol. 1, Geography,
Physics, and Chemistry, John Wiley and Sons, Inc.,
New York, 1957.
4. Wurtz, C. B.j and Renn, C. E. , "Water Temperature and Aquatic
Life," Department of Sanitary Engineering and Water
Resources, The Johns Hopkins University, June, 1965.
5. McKee, J. E., and Wolf, H. W., "Water Quality Criteria," State
Water Quality Control Board, Sacramento, California,
Publication No. 3-A, 1963.
6. Laberge, R. H. , "Thermal Discharges," Water and Sewage Works,
Vol. 106, December, 1959, p. 536.
7. Hynes, H. B. N., The Biology of Polluted Waters, Liverpool
University Press, Liverpool, 1960.
8. Collins, W. D., "Temperature of Water Available for Industrial
Use," U. S. Geological Survey Water Supply Paper 520(f),
U. S. Govt. Printing Office, Washington, D. C.
9. Sylvester, R. 0., "Effects of Water Uses and Impoundments on
Water Temperature," Proc., 12th Pacific Northwest
Symposium on Water Pollution Research, Corvallis,
Oregon, 7 November 1963, p. 6.
- 155 -
-------
10. Churchhill, M. A., "Effects of Density Currents in Reservoirs
on Water Quality," Journal of Water and Sewage Works,
30 November 1965.
11. "A Ten Year Hydro-Thermal Power Program for the Pacific North-
west," Bonneville Power Administration, Portland, Oregon,
January 1969.
12. "Summary Report on Nuclear Power Plant Siting in the Pacific
Northwest for Bonneville Power Administration," Battelle
Memorial Institute, Pacific Northwest Laboratories,
Richland, Washington, 1967.
13. "The Treaty with Canada for Joint Development of the Columbia
River," Bonneville Power Administration, Portland,
Oregon, 1 December 1965.
14. "A Proposal for Formulation of a General Mathematical Model
for the Pre-diction of Thermal Energy Changes in
Impoundments," Presented to the Federal Water Pollution
Control Administration Thermal Effects Project, by
Water Resources Engineers, Inc., Walnut Creek, California,
7 May 1968.
15. "Prediction of Thermal Energy Distribution in Streams and
Reservoirs," Prepared for the Department of Fish and Game,
State of California, by Water Resources Engineers, Inc.,
Walnut Creek, California, Revised Edition, 30 August 1968.
16. "Applications of Mathematical Models for Prediction of the
Thermal and Quality Behavior of Lake Washington," Prepared
for the Washington Pollution Control Commission by Water
Resources Engineers, Inc., Seattle, Washington, November,
1968.
17. Neuman, G.3 and Pierson, W. J., Jr., Principles of Physical
Oceanography, Prentice-Hall, Inc., Englewood Cliffs,
New Jersey, 1966.
18. Richardson, L. F., "Some Measurements of Atmospheric Turbulence",
Phil. Trans. Roy. Soc.3 A, Volume 221, London, 1920.
19. Yih, C-S., "On the Flow of a Stratified Fluid," Proc., 3rd
U. S. National Congress of Applied Mechanics, 1958.
20. Debler, W. R., "Stratified Flow into a Line Sink," Journal of
the Engineering Mechanics Division,, ASCE, Vol. 85, EMS,
July 1959.
- 156 -
-------
21. Brooks, N. H., and Koh, R. C. Y., "Selective Withdrawal from
Density Stratified Reservoirs/'Preprint from Proceedings
of ASCE Specialty Research Conference on "Current Research
into the Effects of Reservoirs on Water Quality," Portland,
Oregon, 22-24 January 1968,
22. Yin, C.-S., Dynamics of Non-Homogeneous Fluids^ The Macmillan
Co., New York, 1965.
23. Craya, A., "Theoretical Research on the Flow of Non-Homogeneous
Fluids," La Eouille Blanche, 1949.
24. Duncan, W., Harleman, R. F,, and Elder, R. A., "Internal Density
Currents Created by Withdrawal from a Stratified Reservoir,"
A Cooperative Study, Tennessee Valley Authority and U. S.
Corps of Engineers, Norn's, Tennessee, February 1962.
25. Jaske, R. T., A Three Dimensional Study of Parameters Related
to- the Current Distribution in Lake Roosevelt, Battelle
Memorial Institute, Pacific Northwest Laboratory,
Richland, Washington, April 23, 1969.
26. Jaske, R. T. and Snyder, G. R,, "Density Flow Regime of
Franklin D. Roosevelt Lake," Journal of the Sanitary
Engineering Division^ ASCE, Vol. 93, SA3, June 1967.
27. Long, R. R., "Velocity Concentrations in Stratified Fluids,"
Journal of the Hydraulics Division^ ASCE, Vol. 88, HY1,
January 1962.
- 157 -
-------
APPENDIX A
MATHEMATICAL MODELS FOR THE
PREDICTION OF THERMAL ENERGY
CHANGES IN IMPOUNDMENTS
GRAPHICAL OUTPUT ROUTINES
FOR OUTPUT FROM PROGRAM TSIP
A - 1
-------
INTRODUCTION
In an effort to enhance the interpretation of the output
from program TSIP, two graphical output routines have been included as
a part of the computer software package. The first of these routines
generates, on the computer system's line printer, a plot of segment
temperature vs. depth. This display is produced for each day that printed
output is requested. The second output routine, which is a complete
and separate program, is programmed to produce pen plots of segment
temperature vs. depth for specified days, as well as plots of the
temperature and flow of the reservoir system's downstream releases.
The pen plots are written for the Calcomp plotting system, and one
must have access to such equipment in order to use these routines.
The details and requirements for use of these routines are given on
the following pages.
A - 2
-------
LINE PRINTER PLOT
The printer plot will produce a temperature vs. depth
plot for each segment for each day that printed output is requested
from program TSIP. These plots are produced automatically, and
no changes in program input are required for their production.
An example of the format of this plot is shown on the next page,
and a FORTRAN listing of the plot producing routines is included
on subsequent pages. At the present time this plotting routine
is included as an integral part of program TSIP. If machine
storage requirements become critical to the implementation of
TSIP, these routines can be removed without any adverse effects
on the actual thermal simulation procedure. The only reference
to these routines is the call to SUBROUTINE CURVE from the output
section of SUBROUTINE SUBB.
A - 3
-------
TEMPERATURE VERSUS DEPTH FOR JULIAN O.AY 21Q
METERS
ABOVE
BOTTOM
200.000 I
I
T
I
I
I
I
I
I
I
160.000 -
I
I
T
I
I
I
I
I
I
120.000 -
I
I
I
I
I
I
I
I
I
80.000 -
I
I
I
I
I
I
I
I
I
************************
********************
*********
******
* ** *
** *
* *
* *
**
* *
* *
* *
*
**
*
I
I
I
I
I
I.
I
I
I
G .r
**
*
EXAMPLE PLOT OUTPUT
LINE PRINTER PLOT
PROGRAM TSIP
.000
.*-_! i-
P..O 10.0
12.0 lt.0 lfi.0
WATER TFUPERATURF * OFGRECS CFNTIGRAnr
18.0
20.0
22.0
1
26.0
-------
FORTRAN PROGRAM LISTINGS
Printer Plots for Program TSIP
Listings for Subroutines
1. CURVE
2. PINE
3. PPLOT
4. SCALE
5. LABD
A - 5
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
SU8ROUTINE CUR VE (X »Y »NPT »NCV tNPLOT)
DIMENSION X ( 203t 1) »Y (203 t 1 ) tNPTd )
CO MV ON /LAB/ TITLE( 18 )f XLAB( 11 >.YLAB (£)
1 tH09 IZ (20) »VERT (6)
**** CURVE IS THE ENTRY TO A GENERALISED WINTER PLOT ROUTINE. THE
ROUTINE PLOTS SEQUENTIALLY PAIRED VALUES TAKEN FROM THE X AND
Y ARRAYS. THE SCALING VALUES FOR ROTH ARRAYS ARE STORED IN THE
LAST TWO ARRAY LOCATIONS IN THE SAME MANNFR AS CALCOMP SCALING.
THE ARGUEMENTS IN THE CALLING SEQUENCE APE DEFINED BELOW...
X...THE ARRAY CONTAINING THE X AXIS COORDINATES OF THE
POINTS TO BE PLOTTED.
Y...THE ARRAY CONTAINING THE Y AXIS COORDINATES OF THE
POINTS TO BE PLOTTED.
NPT...THE NUMBER OF POINTS TO BE PLOTTED.
NCV...THIS VALUE IS ALWAYS ONE (1).
NPLOT...USED FOR PLOT IDENTIFICATION THIS VALUE
ORINTED ABOVE FACH PLOT FOR EACH CALL TO
IS
CURVE.
NPTS-NPT (1 )*NCV
AXLEN^IO.
CALL SCALE(X .AXLEN»NPTSt 1 )
AXLENZ5.
CALL SCALF ( Y » AXLEN iNPTS* 1 )
XM TN=X (NPTS+ 11 1)
OELTX=X ( NPTS+2 tl )
XL A3 (1 )rXM!N
DO ?60 l-l .10
260 XL A" (1+1 )=XL AB(I ) + DELTX
SET UP X ANP Y SCALES
FORM X LABELS AND FACTORS
r
c
c
c
c
c
FORM Y LABrLS AND FACTORS
YMIN=Y(NPTS+ltl)
DELTYrY(NPTS+2 »1)
YL Aq (6 )=YMIN
DO 270 lrl,5
270 YL AB (6-1 J =YLAB (7-1 )
YSCAL-5a./(YLAB(l)-YMIN)
INITIALIZE PLOT OUTLINE
A - 6
-------
CALL PPLCT(Q»0tNCDtNPLOT)
K r 1
C
C HRflW I N EACH CU9VE
C
00 450 L-l»NCV
IF
-------
SUBROUTINE PINF.(Xl*YltX?tY2fNSYMtNCTr
AXA-X1
AXB-X2
AYAzYl
AYB-Y2
N-l
IF( ABS (AXB-AXA ) . LT .A8S (AYB-AYAM GO TO 290
C
C SET PARAMETERS FOR X DIRECTION
C
IF (AXB.GT.AX A) GO TO 2^5
AX ArX2
A X 8 - X 1
AYftrY2
AYRrYl
245 CONTINUE
IXA-AXA+.5
IX8=AX8*.5
lYArAYA+.5
IYB:AYB*.5
250 CONTINUE
IF (TXA .LT.O. OR .TXA .GT. 100) GO TO 260
IF (IYA.LT.O.OR .IYA.GT.50) GO TO 260
CALL PPLOTdXA »IYA»NSYM»NCT )
250 CONTINUE
IXAr IXA-H
YAr(N*(AYB-AYAM/(AXB-AXA)
IYArAYA+YA+0.5
N = N+1
IF(IXA.LE.IXB) GO TO 250
GO TO 400
C
C SET PARAMETERS FOR Y DIRECTION
C
290 CONTINUE
IF(AYB.GT,AYA) GO TO 295
A Y B - Y 1
AYA-Y2
AXB-X1
AXA-X2
295 CONTINUE
IXA=AXA+.5
IXB-AXB+ .5
IYA=AYA+.5
IYB=AYB+.5
300 CONTINUE
IFdXA.LT.O.OR .TXA .GT. 100) GO TO 310
IF(I YA.LT.O.OR.IYA.GT.50) GO TO 310
CALL PPLOT(IXAtlYAfNSYM.NCT)
310 CONTINUE
IYArIYA+1
XA=(N*{AX8-AXA))/(AYB-AYA)
IXArXA+AXA+0.5
A - 8
-------
N-NH-1
IF(TYA-IV8I 30D»320.fOO
320 IXA = 1X8
GO TO 300
«»00 RETURN
END
A - 9
-------
SUBROUTINE PPLOT IIX» IV ? KtNCT )
DIMENSION A(51.101> »SYMI 9)
COMMON /LAS/ TITLE < 18 ) »XLA<3 (1 1 ) .YLAfMfi >
1fHO^TZ(2G)»VERT(6)
DATA SYM / 4H****»4H++++» 4H*'**. 4HXXXX* 4H....» 4H2222.
1 4H , 4HIIII, 4H /
IF(K-gS) 200.220.230
200 A(51-IY.IX+1)-SYM{K>
RETURN
?2P CONTINUE
1-0
WRITE16.103) TITLE.MCT
DO 225 11=1.6
1 = 1+1
WRITE(B.lOl) YLA3(II).(A(T.J).J=l?ini)
IF(II.EO.G) GO TO 228
:DO 224 JJ-1.9
1=1 + 1
TFd.NE.28) GO TO 221
WRTTE(6tl06) VERT(5)»vrRT(6).(fl(I»J)»Jrl,101)
GO TO 224
271 IFCI.NE.?4) GO TO 222
WRTTE(B.infe) VERT(l».VERT(2).(j5(I,J)rJ = l ?1 01 )
GO TO 224
222 IFCI.NE.2G) GO TO 221
WRTTE(6.in6) VERT(3)»VERTC4),(A(ItJ),J=l,101)
GO TO 224
2?3 wRirE(6,ina) (A(i,j).j = i,101>
224 CONTINUE
225 CONTINUE
228 CONTINUE
WRITE(£.102) XLAR
WRITE
101 FORMAT(F17.3.1X.101A1)
102 FORMAT(^20.1 .1 OF 10-1 )
103 FORMATUHl.20Xtl8A4.I6/)
105 FORMAT
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
SUBROUTINE SCALE (ARRAY»AXLENtN"TStINC)
DIMENSION ARRAY(NPTS ) »INT(5 )
DATA INT»/2»4t5fS»10/
INCT-IARS(INC)
250
255
260
2&5
275
C
C
c
c
c
c
280
300
320
SCAM F 09 MIX AND M IN
AMAX zARRAY( 1 )
AMINrARRAY( 1 )
DO 250 N = l »NPTSf INCT
IF(AMAX.LT.ARRAY (N ) )
IF(AMIN.GT.ARRAY IN) )
CONTINUE:
IF( AMAX - AMIN ) 275»25
A M A X - A R R A Y ( N )
AMINrARRAY(N)
>?75
RFSFT Mft X AMD MIN FOR 7FPO PflMGF
IFt AMIN ) 265» t*00t 2&0
AMIN - 0.0
AMAX - 2.0 * AMAX
GO TO 275
AMAX - 0.0
AMIN = 2.0 * AMIN
CONTINUE
RATE-(AMAX-AMIN)/AXLEN
A-LOG10
-------
IF (AMIN.LT .0.) KrK-1
C
C CHECK FOR ^AX VALUE IN RANGf
C
IF (AMAX.GT .(K*AXLEN)*RANGE) GO TO 330
I:NPTS*INCT+1
ARRAY(I)rK*RANGE
I-I*INCT
ARRAY(I)-RANGE
RETURN
C
C TF OUTSIDE RANGE RESET L A NO
C
330 L = L+1
IFCL.LT. 11 ) GO TO 280
L-?
3 It? GO TO 230
C
C SET UP NEGATIVE STEPS
C
IF(AMAX.GT.Q.) KrK+1
IF (AMIN.LT. ( K-«-AXLEN>*R ANGE) GO TO 330
I = INCT*NPTS-«-l
ARRAY(I)-K*RANGE
ARRAY(I»=-RANGE
RETURN
400 WRITElGf100)
100 FORMATf // lOXf 'RANGE AND SCAL.F ARE ZERO ON PLOT ATTEMPT* )
RETURN
END
A - 12
-------
BLOCK DATA
COMMON/LAB/ TTTLE(18)fXLA8(ll)»YL4B(F)
1»HORIZ(20)»VERT(6)
DATA VERT/4HMETE .UHRS . <* H ABOV t <4HE , <* H3 OT T» HH OM /
DATA TITLE/8*4H tMH TEM i UHPFR A 14 HTURE »4 H VERitHSUS »tHDEPT»
A 4HH FOttHR JUttHLIAN* 4H DAY /
DATA HORIZ/tHWATE» 4HR TF »<4HMPER »<4HA TUP »t HE * . tHDE GR 14 HEES t
A»»HCENTf 4HIGR At 4HDE »10*MH /
END
A - 13
-------
CALCOMP PEN PLOT
The pen plotting routine graphs certain requested por-
tions of the output from program TSIP. In order to use this routine,
an output file (normally a magnetic tape) must be written during the
execution of program TSIP. The production of this file is a user
defined option, and is accomplished by placing a positive value in
columns 61-70 of card number 1 of Card Group II in the input to TSIP.
The file thus produced, together with the small amount of card input
explained below, constitutes the total input for this program. Examples
of the plots produced by the program and its FORTRAN listing are shown
subsequent to the following explanation of plot's input and output.
A - 14
-------
INPUT
The card input to this plotting routine consists of two
or more cards as explained below:
CARD NO.
CARD COLUMNS
DATA DESCRIPTIONS
1-80
(free field)
Comment to be printed on line
printer.
1-5
6-10
Logical Unit Number upon which
the tape written by TSIP is to
be mounted for input to the
plot routine.
Logical Unit upon which Calcomp
tape is to be written on output,
11-15
The total number of days for
which segment profile plots
are to be produced.
3,4,5*
1-80
(16 fields of
5 columns each)
A list of the days for which
profile plots are desired.
Up to 48 days may be requested.
* Use only as many cards as
necessary.
A - 15
-------
OUTPUT
The output from this routine are pen plots of the follow-
ing items.
1. Temperature vs. depth for each segment for each day
specified.
2. Temperature vs. time for total downstream release.
3. Discharge vs. time for total downstream release.
4. Temperature vs. time for the discharge from each
individual outlet.
5. Discharge vs. time for each individual outlet.
A - 16
-------
o
o
o
CM.
CO
TEMPERRTURE PLOT F0R DRY 200
o
o
o
CO.
CM
O
O
O
•*.
CM
O
O
UJ
• o
o
I •
o
.CO.
LUO
,CM.
D
O
O.
00
D
O
0_
O
o
J45.00 50-00 55-00 60-00 65-00
TEMPERRTURE IN DEGREES F
EXAMPLE PEN PLOT OUTPUT
Segment Temperature vs. Elevation
A - 17
70-00
75-00
-------
o
o
D0WNSTRERM DISCHRRGE
120-00
140.00
160-00
180-00 200.00 220-00 240-00
TIME IN JULIRN DRYS
EXAMPLE PEN PLOT OUTPUT
Downstream Temperature vs. Time
260-00
280.00
300.00
-------
o
o
0-,
o
o
•
o
oo
LUQJ
LU
CD
•o.
(£>
LU
»— O
CE°
LULO
Q_
LU
O
o
OJ
o
o
DQWNSTRERM TEMPERRTURE
'120.00 140.00 160-00 180-00 200.00 220.00 240.00
TIME IN JULIRN DRYS
EXAMPLE PEN PLOT OUTPUT
Downstream Discharge vs. Time
260.00
280.00
300.00
-------
o
o
ro
o
DISCHRRGE FROM OUTLET NUMBER 1
"120-00
140-00
160-00
180-00 200-00 220-00 240-00
TIME IN JULIflN DRYS
EXAMPLE PEN PLOT OUTPUT
Individual Outlet Temperature vs. Time
260-00
280-00
300
-------
o
o
o_
o
o
o
oo
cng
LU .
LUOJ
LU
o
o
—-o.
(O
LU
ct:
ID
1—0
or°
Q_
SI
LU
V-o
o
OJ
o
o
o
CO.
TEMPERRTURE OF OUTLET NUMBER 1
120-00 140-00 160-00 180.00 200-00 220-00 240-00
TIME IN JULIRN DRYS
EXAMPLE PEN PLOT OUTPUT
Individual Outlet Discharge vs. Time
260-00
280-00
300-00
-------
FORTRAN PROGRAM LISTING
CALCOMP PEN PLOT FOR OUTPUT FROM PROGRAM TSIP
A - 22
-------
C ********* CALCOMP PEN PLOT FOR OUTPUT FROM PRO-AM TSIP
r
DIMENSION ISUFdOaO). IDOUTC50). X(3£3). Y{3&*), QOTC368.5).
A TOUT(368i5)f DST(368)» 050(368). MFAD(?H)» T03J(368)
C **** READ AND WRITE CARD IMFORMATTON
REAO(5»501) HEAD
WRITE(G»601) HEAD
REAO(5r503) ITAPEt NTP. NSD
WRITE(G»603) TTAPEt NTP. NSD
IF ( NSD .LE. 0 ) GO TO 103
READ(5»503) ( IDOUT(J)t J r 1, NSO )
WRITE(6f6Q5) I TDOUT(J)t J - 1. NSD )
C **** SETUP CALCOMP ORIGIN
103 CALL PLOTSJIBUF. 1000.NTP )
CALL PLOT< 0.0.-30.0.-3)
CALL PLOT( 0. Ot 1. 5.-3 )
CALL SYMBOL( 0. 0. G.0.0. 21 .25HPLOT FOR Fypcfl S TU OY TEAM
A 90.0.25 )
CALL PLOT( 4.0. 0. 0.-3 )
C **** SETUP AND PLOT REQUESTED TEMPERATURE PROFILES
105 REAO(ITAPE) IDAY. LDAY. SDZ» NOUTSt NSEGt NLrTS
WRITE<6»607) IDAY» LDAY. SDZ» NOUTS. NSEC=. NL^TS
Y(1) - 3.2808 * 0.5 * SDZ
TA - 3.2808 * SOZ
IDAY r IOAY - 1
DO 125 L - IDAY. LDAY
REAO(ITAPE) ID. NP TS t ( X«J). J - 1. M PT S )
IF ( NSD .LE, 0 > GO TO 125
DO 107 K - 1. NSO
IF ( IDOUT(K) .EQ . 10 ) GO TO 111
107 CONTINUE
GO TO 125
111 DO 115 J - 2. NPTS
Y(J) r Y(J-l) + TA
115 CONTINUE
DO 1 19 J - 1 » NPTS
X ( J) - 1.8 * X (J) + 32.0
119 CONTINUE
CALL SCALE( Y( 1 ) . 8.0.NPTS .1 )
CALL SCALECXd > .6.0. NPTS .1 )
CALL AXIS( 0. 0. 0. 0. 2UHTE^PER ATURE IN DEGREES F» - 24 t 6, 0 » 0. 0 »
A XfNPTS+1)»X(NPTS+2) )
CALL AXIS(0.0»0.Ot 17HELEVATION IN F EE T , 17,8.0.90.0 ,
A YtNPTS+1 )fY(NPTS + 2) )
CALL L INE(X. Yt NP TS » 1 tO»0)
CALL S.YMBOLC 1.04 .8.00. 0 . 11 .2UHTEMP ER A TU RE PLOT FOR DAY.O.O»?4
FPN = ID
A - 23
-------
CALL NUMBER(4.5«4.R.,aat0.m»FPN»0.a.-l)
CALL PLOT ( lO.OtO .0 » - 3 )
125 CONTINUE
REAO QTAPE .END^l 27) TA
127 IF ( NSEG . NE . NLFTS ) GO TO 105
WRITE(6»609»
**** INPUT AND CONVERT DOWNSTREAM VALUES
NPTS - LDAY - IDAY + 1
REAOCITAPEJ ( < QOTCJ»K)»K - 1» NOUTSJt J - 1» NPTS J
REAOIITAPE5 ( C TOUTJJ»K)» K - 1» NOUTS)» J - It NPTS )
DO 1««5 J - 1 » NPTS
DSQCJJ - 0.0
T A ~ 0.0
DO 111 K ~ It NO UTS
QOTfJtKJ - 35.3115 * QOT(JtK)
TOUTUrK) r 1.8 * TOUT(J»KJ + 32.0
DSO(J) = OSQ(J) + QOTCJ.KS
TA r TA + QOT(J,K» * TOUTJJ»K)
IHI CONTINUE
DST(J) r TA / DSQKJ)
145 CONTINUE
**** GENERATE THE TIME AXIS
X ( 1) - 0.5
DO If 9 K - 2 » NPTS
XI K> - X (K-l ) + 1.0
IHS CONTINUE
TA - IDAY - 1
CALL SCALE J X.9.0.NPTS* 11
**** PLOT THE FLOWS AND DISCHARGES
DO 155 K - 1 , NO UTS
TOUT(NPTS+l.K) - 30.0
TOUT(NPTS+2*K9 - 10.0
CALL SCALEJQOT II »K ) t 6.0.NPTS tl )
155 CONTINUE
CALL SCALE IDSQ (1 )» 6.0»NPTS tl 1
MAX - NOUTS + 1
DO 167 L : It MAX
CALL AXIS( 0.0» 0. 0» 19HTIME IN JULIAN D 5 YS t- 19 .9 .25 » 0 . 0 »
fl Tfl »X(NPTS + 2) J
CALL AXIS(0,0» 0. 0, 2«* HTEM PER A TURE IN DEGREES F( 24 . 6. 01 90 . 0 t 3 0, Of
A 10.0 8
IF I L .GT. 1 5 GO TO 159
**** PLOT TOTAL DOWNSTREAM VALUES
CALL L INE(X»DST» NP TS» 1 .O.OJ
CALL SYMBOLS 3, 08.6.00tO. ]«* » 22HDO WNS TPE AM TEMPERATURE tO.0.22 )
A - 24
-------
CALL PLOT! 13 .25t 0, Ot-3>
CALL AXISC O.Of 0. Ot 19HTIMF IN JULIAN DAYSf-19 i 9 .25100 .0f
A TAt X(NPTS-»-2) )
CALL AXIS! 0.0? 0. Ot 16HDISCHARGE IN C FS . ] 6 *S . 0 f9 0 . 0 »
A DSQC NPTS + 1 ) f DS OCNPTS+2) )
CALL LINEC XfOSQ fNPTSt If Of 0 )
CALL SYMBOH 3. 22f 6.00f 0. 14 f 20HnnWNSTPEAM 0 ISCHARGE»Q.0f20 )
GO TO 163
**** °LOT DOWNSTREAM COMPONENTS
159 K = L - 1
FPN r K
CALL L INE( X t TOUT (1 tK )t NPTS tl t Ot 0)
CALL SYMBOLC 2. 52 f6 .OOf 0. 14 f 28HTEMPERATURE OF OUTLET N!U MBER »0 .T ?28 )
CALL NUM8ER(6.58.6.00f0.1<4tFPNf0.af-l)
CALL PLOT ( 13 ,25f T,Of-3 J
CALL AXIS( 0. Of 0. Of 19HTTME IN JULIAN 0 A YS »- 19 f 9 .25 f 0 . 0 f
A TAfX(MPTS+2) )
CALL AXISC O.Of 0. Of 16HDISCHARGE IN C FS f 16 f 6 . Of 90 .0 f
A GOT C NPTS + 1 fK ) » 30 T ( NP TS+2 f K ) >
CALL L INE( Xf OOT( If K) f NPTSf 1 fO fO )
CALL SYMBOLC 2. 52 fG.OO. 0. 14 f 28HDISCH AP^E FROM OUTLET NU MR ER f 0 .0 ,28 )
CALL NU^BER(6.58f6.00f0.14fFPNf0.nf-l)
163 CALL PLOT( 1 3 .25f O.Of-3)
WRTTE(6f613) L
167 CONTINUE
REWIND ITAPE
STOP
**** INPUT FORMAT STATEMENTS
501 FORMAT! 20A«4 )
503 FORPATC 1615 )
**** OUTPUT FORMAT STATEMENTS
601 FORMAT( 1H1 / 35Xf 55HCALCOMP PLOTTING ROUTINE FOR THERMAL SIMUL4T
1ION RESULTS // 40Xf 20A4 )
603 FORMAT( // 20Xt 10HINPUT TAPE 110 / 13Xf 12HTALCOMP TAPE 110 /
A 17Xf 13HPROFILE PLOTS 110 )
605 FORMATf /// 23Xf 17HPROFTLE PLOT DAYS //( 140 ) »
607 FORMAT( // 19Xf 11HTAPE HEADER// 21Xf 9HFTRST DAY 110 /
A 21Xf 9HFINAL DAY 110 / 23Xf 7HELEMENT F1Q.2 /
R 23Xf 7HOUTLETS 110 /
C 20Xf IOHSEGMENT NO 110 / 20Xf 10HOUTLET SEGMENT 110 )
609 FORMAT( //// 40Xf 22HPROFILE PHASE COMPLETE J
613 FORMAT( // UOXf 15HCOMPONENT PHASE 12f 12H IS COMPLETE )
END
A - 25
-------
1
5
Accession Number
~ Subject Field & Group
02H, 05B
SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Organization
Water Resources Engineers. Incornorat-prI
Walnut Creek, California
Title
MATHEMATICAL MODELS FOR THE PREDICTION OF THERMAL ENERGY CHANGES IN IMPOUNDMENTS,
10
Authorfs)
Orlob, Gerald T.,
Roesner, Larry A., and
Norton, William R.
16
Project Designation
FWPCA Contract No. 14-12-422, 16130 EXT
21
Note
22
Citation
Final Report, FWPCA Contract No. 14-12-422, FWQA Report No. 16130 EXT 12/69
23
Descriptors (Starred First)
*Reservoirs, * Impoundments, *Water temperature, ^Thermal stratification,
*Mathematical models, *Heat budget, Computer programs, Columbia River,
Stratification, Diffusion, Diffusivity
25
Identifiers (Starred First)
*Temperature prediction, *Selective withdrawal, Lake Roosevelt
071 Abstract^- mathematical model developed previously by Water Resources Engineers, Inc., was
refined and extended to produce a more generalized model capable of describing the thermal
regime of a stratified reservoir under various hydrologic,meteorologic,and hydraulic conditions.
Two important refinements of the original model were made: (1) A generalized empirical ex-
pression for the eddy conductivity coefficient was derived, and (2) The influence of reservoir
stability on the vertical extension of the withdrawal zone for reservoir outlets and the zone
of influence for inflows was accounted for.
The refined model accounts for external heat exchange by the usual heat budget components.
Internal heat transfer is accomplished by the penetration of short-wave radiation, eddy dif-
fusion, and vertical advection. The model was tested and verified on Hungry Horse Reservoir
and is applicable on reservoirs for which the assumption of horizontal isotherms is valid. A
method for testing this assumption is given.
The model described above was extended to simulate weakly stratified reservoirs (tilted iso-
therms) . The weakly stratified reservoir is represented as a set of smaller reservoirs (or
segments) that are coupled to one another by heat and mass continuity. Each small reservoir is
described by the basic model,and the tilted isotherms are produced by connecting the tempera-
ture profiles at the longitudinal midpoints of the segments. A test of the segmented model on
Lake, Roosevelt gave satisfactory resultsjhowever,its general applicability remains to be
rais report was submitted in fulfillment of project 16130 EXT,Contract No. 14-12-422. under
the sponsorship of the Federal Water Quality Administration.
Abstractor
Larry A. Roesner
Institution
Water Resources Engineers, Incorporated
WR:I02 (REV, JULY'1969)
WRSI C
SEND TO: WATER RESOURCES SCIENTIFIC INFORMATION CENTER
U.S. DEPARTMENT OF THE INTERIOR
WASHINGTON. D. C. 20240
* GPO: 1969—359-339
------- |