WATER POLLUTION CONTROL RESEARCH SERIES • 16130 DJH 04/71
TEMPERATURE PREDICTION
IN
STRATIFIED WATER:
MATHEMATICAL MODEL-USER'S MANUAL
U.S. ENVIRONMENTAL PROTECTION AGENCY
-------
WATER POLLUTION CONTROL RESEARCH SERIES
The Water Pollution Control Research Series describes the
results and progress in the control and abatement of pollution
in our Nation's waters. They provide a central source of
information on the research, development and demonstration
activities in the Environmental Protection Agency, through
inhouse research and grants and contracts with Federal,
State, and local agencies, research institutions, and
industrial organizations.
Inquiries pertaining to Water Pollution Control Research
Reports should be directed to the Chief, Publications Branch
(Water), Research Information Division, R&M, Environmental
Protection Agency, Washington, D.C. 20^60.
-------
TEMPERATURE PREDICTION IN STRATIFIED WATER:
MATHEMATICAL MODEL-USER'S MANUAL
(Supplement to Report 16130DJH01/71)
by
RALPH M. PARSONS LABORATORY
FOR WATER RESOURCES AND HYDRODYNAMICS
Department of Civil Engineering
Massachusetts Institute of Technology
Cambridge, Massachusetts 02139
for
ENVIRONMENTAL PROTECTION AGENCY
Research Grant No. 16130 DJH
April 1971
For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402 - Price $1.25
-------
EPA Review Notice
This report has been reviewed by the Environmental Protection
Agency and approved for publication. Approval does not
signify that the contents necessarily reflect the views and
policies of the Environmental Protection Agency nor does
mention of trade names or commercial products constitute
endorsement or recommendation for use.
ii
-------
ABSTRACT
The annual cycle of temperature changes in a lake or reservoir may be quite
complex, but predictions of these changes are necessary if proper control
of water quality is to be achieved. Many lakes and reservoirs exhibit hori-
zontal homogeneity and thus a time-dependent, one-dimensional model which
describes the temperature variation in the vertical direction is adequate.
A discretized mathematical model has been developed based on the absorption
and transmission of solar radiation, convection due to surface cooling and
advection due to inflows and outflows. The mathematical model contains pro-
vision for simultaneous or intermittent withdrawal from multi-level outlets
and time of travel for inflows within the reservoir. Heat transport by tur-
bulent diffusion in the hypolimnion is neglected.
Good agreement has been obtained between predicted and measured temperatures
in both the laboratory and the field. Field verification consisted of the
simulation of the thermal structure of Fontana reservoir during a nine-month
period. Criteria for the applicability of the model are given. The mathematical
model is a predictive one, since the required data is that which would normally
be available before the construction of a reservoir.
Emphasis has been placed on a detailed explanation of the physical basis for
the mathematical model and on the computer program inasmuch as the report is
intended primarily as a user's manual.
This report was submitted in fulfillment of Research Grant No. 16130 DJH
between the Water Quality Office, Environmental Protection Agency and the
Massachusetts Institute of Technology.
Key words: reservoir temperature distribution; thermal stratification
in reservoirs; reservoir water quality.
-111-
-------
TABLE OF CONTENTS
Page
ABSTRACT iii
TABLE OF CONTENTS v
LIST OF FIGURES viii
FORWARD x
!_. INTRODUCTION 1
1.0 Thermal Stratification in Reservoirs 1
1.1 Proposed Model 3
2. DESCRIPTION OF PHYSICAL BASIS OF MODEL 4
2.0 Introduction 4
2.1 Schematization of Reservoir 4
2.2 Assumptions 4
2.3 Important Model Characteristics 8
2.3.1 Variable Area 8
2.3.2 Direct Absorption 8
2.3.3 Convection 10
2.3.4 Inflow and Outflow 10
2.3.5 Calculation of Inflow Velocity Distribution 10
2.3.6 Calculation of Outflow Velocity Distribution 12
2.3.6.1 Laboratory Case -,»
2.3.6.2 Field Case 15
2.3.6.3 Outlet Geometry 16
2.3.6.4 Multiple Outlets 16
2.3.7 Vertical Advective Velocity 17
2.3.8 Travel Time for Inflows 17
2.3.9 Variable Surface Level 21
2.3.10 Entrance Mixing 21
2.4.1 Derivation of Basic Equation 27
2.4.2 Initial and Boundary Conditions 28
-v-
-------
3. DESCRIPTION OF NUMERICAL MODEL
3.0
3.1
3.2
3.3
3.4
Solution of Heat Transport Equation
3.0.1 Choices of Numerical Scheme
3.0.2 Description of Numerical Scheme
3.0.3 Limitations on Explicit Scheme
Formulation of Explicit Finite Difference Scheme
3.1.1 Internal Element
3.1.2 Surface Element
3.1.3 Bottom Element
Convective Mixing
Energy Check
Numerical Dispersion
4. DETERMINATION OF PARAMETERS
4.1
4.2
4.3
4.4
4.5
4.6
Net Solar Radiation
Net Longwave Radiation
Evaporation and Conduction Losses
4.3.1 Field Evaporation
4.3.2 Laboratory Evaporation
4.3.3 Negative Evaporation
Required Input Data
Selection of Time and Distance Steps
4.5.1 Variable Time Step
4.5.2 Lower Limit on Ay in Field Reservoii
Possible Program Changes
4.6.1 Cut-Off Gradient
4.6.2 Entrance Mixing
5. CRITERIA FOR APPLICABILITY OF MODEL
5.1
5.2
Wind Mixing Criterion for Lakes and Reservoirs
Flow Through Criterion for Reservoirs
6. VERIFICATION OF MATHEMATICAL MODEL FOR FIELD RESERVOIR
6.1
6.2
6.3
6.4
Description of Reservoir
Input Data
Results
Verification from Other Sources
Page
29
29
29
29
30
30
30
34
37
39
40
40
43
43
43
44
45
46
47
47
47
48
49
49
49
49
50
50
51
53
53
53
54
54
-VI-
-------
Page
7. THE COMPUTER PROGRAM 61
7.1 Main Program - Flow Chart 61
7.2 Function FLXOUT - Surface Losses 69
7.3 Subroutine TOUT - Outflow Temperatures 71
7.4 Subroutine SPEED - Velocities 72
7.5 Subroutine AVER - Convective Mixing 74
7.6 Remaining Subprogram 75
7.7 Input Variables 75
7.8 Program Listing 81
7.9 Input Data 108
ACKNOWLEDGEMENT 119
BIBLIOGRAPHY 120
DEFINITION OF NOTATION 122
-V11-
-------
Figure
List of Figures
Page
1.1 Seasonal Variation of Temperature Profiles in Fontana ^
Reservoir, 1966
2.1 Schematization of Reservoir
2.2 Control Volume
2.3 Effect of Variable Area on Temperature Predictions in a
Clear Lake
2.4 Transmission of Solar Radiation Vs. Depth 11
2.5 Dye Concentration Profiles in a Tributary of Fontana
Reservoir
2.6 Outflow Temperatures. Alternate Withdrawal from Multiple
Outlets
2.7 Outflow Temperatures. Simultaneous Withdrawal from Multiple
Outlets
2.8 Calculation of Vertical Velocity 17
2.9 Effect of Entrance Mixing on Outflow Temperatures. No
Insolation. Constant Inflow, Outflow 23
2.10 Effect of Entrance Mixing on Outflow Temperatures.
Variable Insolation. Constant Inflow, Outflow 24
2.11 Effect of Entrance Mixing on Outflow Temperatures.
Variable Insolation. Variable Inflow, Outflow 25
2.12 Effect of Entrance Mixing on Outflow Temperatures.
Fontana Reservoir, 1966 26
3.1 Internal Element 31
3.2 Surface Element 34
3.3 Bottom Element 37
3.4 Effect of Numerical Dispersion on Outflow Temperatures 42
5.1 Effect of Wind on Thermocline 51
6.1 Map of Fontana Reservoir and Watershed 55
6.2 Measured and Predicted Outflow Temperatures Through
Fontana Dam 56
-viii-
-------
Figure Pa§e
6.3 Measured and Predicted Temperature Profiles, Fontana
Reservoir, June 22, 1966 "
6.4 Measured and Predicted Temperature Profiles, Fontana
Reservoir, July 20, 1966 58
6.5 Measured and Predicted Temperature Profiles, Fontana
Reservoir, September 15, 1966
6.6 Measured and Predicted Temperature Profiles, Fontana
Reservoir, November 10, 1966 °U
-IX-
-------
FOREWORD
This is the fourth report issued in conjunction with a continuing research
program on thermal stratification and water quality in lakes and reservoirs.
The previous reports are as follows:
1. Dake, J.M.K. and D.R.F. Harleman, "An Analytical and Experimental Investi-
gation of Thermal Stratification in Lakes and Ponds", M.I.T. Hydrodynamics
Laboratory Technical Report No. 99, September 1966. (Portions of this
report have also been published by the same authors under the title:
"Thermal Stratification in Lakes: Analytical and Laboratory Studies",
Water Resources Research, Vol. 5, No. 2, April 1969, pp. 484-495.)
2. Huber, W.C. and D.R.F. Harleman, "Laboratory and Analytical Studies of
the Thermal Stratification of Reservoirs", M.I.T. Hydrodynamics Laboratory
Technical Report No. 112, October 1968.
3. Markofsky, M. and D.R.F. Harleman, "A Predictive Model for Thermal Stratification
and Water Quality in Reservoirs", M.I.T., Department of Civil Engineering,
Ralph M. Parsons Laboratory for Water Resources and Hydrodynamics, Technical
Report No. 134, January 1971. CAlso published in Water Pollution Control Research
Series No. 16130 DJH 01/71 by Environmental Protection Agency, W.Q.O., Wash. B.C.)
This report is intended as a user's manual for the transient temperature distribu-
tion model developed in the second report (Technical Report No. 112). The computer
program presented herein supersedes that given in Technical Report No. 112. The
numerical scheme has been changed from an implicit to an explicit scheme. Other
changes include the option of accounting for the travel time of inflows within
the reservoir, simultaneous or intermittent withdrawal from multi-level outlets
in the reservoir and the continuous variation of the water surface elevation.
The emphasis in this report is placed on an explanation of the computer program
rather than on the development and testing of the theory.
H V *•
-------
1. Introduction
1-0 Thermal Stratification in Reservoirs
Thermal stratification occurs in almost all lakes and reservoir impoundments.
In shallow "run of the river" reservoirs the stratification may be relatively weak
and in certain seasons the isotherms may be tilted in the downstream direction. In
deep lakes and in reservoirs in which the storage volume is large compared to the
annual through-flow, the isotherms are horizontal during most of the year and strong
stratification is generally developed during the late summer and autumn seasons. The
mathematical model developed in this study is concerned with the latter situations
in which the water temperature is a function of depth and time.
The primary causes of thermal stratification are the low thermal conductivity
of water, the limited penetration of radiant heat and light, and the fact that stream
inflows in late spring and early summer tend to be warmer than the reservoir surface
waters. These warm inflows spread out over the reservoir surface. Furthermore,
virtually all heat, apart from advected heat, enters the reservoir through the sur-
face in the form of radiant energy. A large percentage is absorbed in the first
few meters and thus the water near the surface is heated more quickly than that in
the lower layers. This warm water tends to remain at the surface, absorbs more
heat, and thus a stable condition tends to be set up. However, evaporation will
always cool the surface layer, setting up convection currents. Surface cooling and
hence convection will be enhanced by back radiation and conduction losses, espec-
ially at night. Wind stresses on the water surface will cause mixing whenever a
neutral or unstable density gradient is set up by surface cooling. These processes
of heating, cooling, and wind action lead to the development of a warm, freely cir-
culating, turbulent upper region, called the epilimnion. This overlays and to a
great extent insulates a colder, relatively undisturbed region called the hypolimnion.
The generality of thermal stratification in temperature climates has been established
since the end of the 19th century (Hutchinson Q.3 ) )•
A typical annual thermal cycle in a reservoir is illustrated in Figure 1.1
which shows the nearly isothermal condition in early spring, the development of
thermal stratification in spring and summer, and the return to the initial condition
in winter. The destruction of thermal stratification is accompanied by vertical
mixing throughout the reservoir, and is usually referred to as the overturn.
Under thermally stratified conditions, with the hypolimnion insulated from
the atmosphere by the epilimnion, renewal of oxygen in the lower layers cannot take
-1-
-------
i
ro
i
0)
p
a)
>
0)
cfl
-------
place. This can lead to anaerobic conditions and low quality water. Anaerobic
decomposition may produce undesirable tastes and odors, and occasionally toxic
effects. During the overturn the mixing of these bottom waters with the rest of
the reservoir may pollute all the water for a short period. Furthermore, release
of this poor quality water could cause a deterioration of water quality downstream
of the impoundment. For these reasons, knowledge of thermal profiles in impounded
waters is essential for water quality control, and prediction of the profiles is
necessary for proper design of outlet works, and a systematic approach to water
quality management.
1.1 Proposed Model
The emphasis has been placed on the development of a predictive model which
uses only such data as is available before the reservoir exists. In the case of
existing lakes or reservoirs, the model can be used to predict future behavior
under various meteorological conditions or due to man-made changes in inflow and/or
outflow conditions.
A discretized mathematical model has been developed, based on the absorption
and transmission of solar radiation, convection due to surface cooling, and on a
realistic treatment of inflows and outflows. The initial verification of the mathe-
matical model was carried out in the laboratory. Details of the laboratory physical
model are given by Huber and Harleman (12). The mathematical model was then extended
to the field case and verified using data from Fontana Reservoir. The model has
also been verified by other investigators using data from Lake Norman in North Caro-
lina (23). The thermal behavior of reservoirs has been simulated over a complete
year allowing for inflow and outflow, and calculating surface heat losses daily as
a function of meteorological data, and the surface temperature generated by the mathe-
matical model. Good agreement has been obtained for both laboratory and field cases.
The mathematical model has been extended to include withdrawal from multi-level
reservoir outlets.
The annual thermal cycle of lakes, with little inflow or outflow, has also
been simulated using a simplified form of the model. This option is available in
the present computer program. The option of using the mathematical model for a
laboratory reservoir has been retained in the computer program. Wherever specific
allowance is made for the laboratory case, this will be pointed out so as to avoid
confusion.
-3-
-------
2. Description of Physical Basis of Model
2.0 Introduction
A comprehensive temperature prediction model must account for internal
radiation absorption, boundary heat sources and sinks, and heat transport by
advection, convection and diffusion. The geometrical simplifications and assump-
tions leading to the development of a comprehensive mathematical model will be
presented. The important characteristics of the model are discussed.and the equa-
tion governing the temperature distribution T(y,t) in a stratified reservoir will
be derived.
2.1 Schematization of Reservoirs
The reservoir is schematized, as shown in Figures 2.1 a,b,c, by considering
it as a series of horizontal elements similar to that in Figure 2.2. The element
at elevation y has a horizontal area A(y) and thickness Ay. A portion of the river
inflow enters the element at the upstream end, and a portion of the outflow through
the dam leaves the element at the downstream end. Heat also enters the element
through the horizontal surface by transmission of radiation, by vertical advection
and by diffusion. The equation governing the temperature distribution will be form-
ulated by considering the conservation of mass and heat within a typical element.
However, before this can be done it is necessary to make certain simplifying assump-
tions, and to consider the mechanisms of radiation absorption, heat advection, con-
vection and diffusion.
2.2 Assumptions
2.2.1 The two principal assumptions involved in the model are:
(a) The isotherms in a stratified reservoir are horizontal, and hence thermal
gradients exist in the vertical direction only. This assumption has the effect of
reducing the problem from a three-dimensional to a one-dimensional, variable area
problem. Observations both in the field and in the laboratory have shown that hori-
zontal isotherms occur in many lakes and reservoirs. However, this assumption limits
the applicability of the model, and certain criteria should be satisfied before using
it. These criteria are discussed in Chapter 5.
(b) Heat transport by turbulent mixing is accounted for only in the epil-
imnion region and only during times at which the temperature induced density pro-
file is unstable. Previous mathematical models for thermal stratification in lakes
and reservoirs have usually included turbulent diffusion coefficients for heat as
important numerical parameters. In general, these diffusion coefficients are functions
-4-
-------
(a) Three-Dimensional View
(b) Control Volume Slice
out,
(c) Side Elevation
Figure 2.1 Schematization of Reservoir
-5-
-------
Internal Radiation Absorption
Heat Source
Inflow Heat Source
Diffusive Heat Source
Outflow
Heat Sink
Advective Heat Source
Control Volume Illustrating Heat Conservation
Figure 2.2 Control Volume
-6-
-------
of depth and time, and since these functional relations cannot be specified a
priori, such mathematical models tend to lose their predictive value. For any
given model it is always possible to find turbulent diffusion coefficients which
will match field data. However, it is difficult to determine to what extent these
coefficients represent diffusion or merely the effect of simplifications or inaccur-
acies in the basic formulation of the mathematical model.
The temperature prediction model developed in this study does not require
assigned values for turbulent diffusion coefficients. Whenever the temperature pro-
file in the epilimnion develops an unstable density gradient, vertical mixing is in-
duced to produce a surface layer of uniform temperature. Thus, even though turbulence
may exist in the surface layer due to wind shear and wave motion, the heat transfer
will be dominated by convection currents and surface cooling effects. These currents
tend to eliminate near surface temperature gradients and nullify the role of gradient
driven turbulent diffusion.
In the hypolimnion region vertical temperature gradients are small and diff-
usive heat transport will not be significant even if turbulence does exist. In the
thermocline region, the density stratification will tend to inhibit turbulence, al-
though it will not necessarily remove it altogether. It should be noted, however,
that calculations of thermal diffusivity in the thermocline region of lakes often
give values of the same order as the molecular diffusivity. Examples are given
by Hutchinson (13), Orlob (19) and Sweers (26)- (See Table 1) These values are
obtained by assuming that all heat transport, not accounted for by other methods,
takes place by turbulent diffusion. Since vertical advective transport is not con-
sidered in any of the above cases, these values are probably too high.
Table 1
Molecular Diffusivity of Heat = 0.0014 cm /sec
Lake
Sodon
Linsley Pond
Mendota
Castle Lake
Ontario
Eddy Diffusivity at
2
Thermocline (cm /sec)
0.007
0.0033
0.025
0.02
0.02 - 0.07 (seasonal
average)
-7-
Reference
Hutchinson (13)
Hutchinson (13)
Hutchinson (13)
Orlob (19)
Sweers (26)
-------
The approach adopted in this study is to neglect turbulent diffusion as
a first approximation, and to take all other known forms of heat transport into
account as accurately as possible. If marked discrepancies occur, which cannot
be explained by other factors, e.g., mixing of the inflow as it enters the reservoir,
then allowance should be made for turbulent diffusive transport. However, in both
laboratory and field reservoirs good agreement has been obtained between predicted
and measured temperatures, and thus it appears that this approach is justified.
Since turbulent diffusion is neglected, it would seem only reasonable to
neglect molecular diffusion as well. This process is included in the general model
for three reasons: 1) Molecular diffusion may be significant in the laboratory case;
2) If accurate values of vertical turbulent diffusivity become available, they may
be included in the model at a later date; 3) A numerical solution to the heat trans-
port equation will be presented, and numerical schemes, regardless of type, behave
better when diffusion is present. It should be mentioned here that numerical schemes
involving advective terms introduce errors which have the same effect as an additional
diffusive term. This effect is discussed in Section 3.4.
2.2.2 Other assumptions, less fundamental than (a) and (b) are:
(c) Solar radiation is assumed to be transmitted in the vertical direction
only. This assumption is not exact but, due to refraction at the water surface and
the fact that when the solar angle is small a large proportion of the short wave radia-
tion comes from the sky; this should not lead to large errors.
(d) The sides and bottom of the reservoir are assumed to be insulated, the
only heat crossing the boundaries (apart from the surface) is via inflow and outflow.
(e) The density (p), specific heat (c) and the coefficient of molecular diffus-
ivity (a) are assumed constant in all heat budget and heat transport calculations.
(f) Solar radiation energy, transmitted by the water and intercepted by the
reservoir sides, is assumed to be distributed uniformly over the cross section at the
depth of interception. This effect may be quite significant in lakes and reservoirs,
which have both clear water and a rapid change of area with depth near the surface
(see Figure 2.3).
2.3 Important Model Characteristics
2.3.1 Variable Area
Variation of area with depth is taken into account. This affects both the
vertical advection of heat, and the distribution of the solar radiation energy.
2.3.2 Direct Absorption
Solar radiation is absorbed directly in the body of the fluid, as well as
-8-
-------
s—'
J3
Temp. °C
10 12 14 16 18 20 22
2
4
10
12
14
16
18
20
22
24
-60-
days
—120
days
B
a
1-D constant area
1-D variable area
Castle Lake
Depth 35 m.
Surface Area 50 acr.
Area at 5 m. depth
27 acres
T = 4°C (end of
°March)
4 Measurements
May 25 1963
0 Measurements
July 27 1963
Data from Bachmann & Goldman (2 )
Figure 2.3 Effect of Variable Area on Temperature
Predictions in a Clear Lake
-S-
-------
at the surface. Transmission of radiation at elevation y is given (as per Dake
and Harleman (5 ) ) by
- 3) e~n(ys ~ y)
where
y = water surface elevation
s
$ = net incident solar radiation (gross-reflected)
o
g = fraction of absorbed at the surface (,-0.4 - 0.5)
n = extinction coefficient of reservoir water, n varies for different
reservoirs and may also vary in time and space. This flexibility is included in the
model. See Figure 2.4 for examples of measured and theoretical curves.
2.3.3 Convection
Convection in the epdlimnion is accounted for by allowing mixing to take place
whenever the temperature gradient (8T/9y) is negative (y is elevation, positive up-
wards). Since surface losses (due to evaporation, back radiation, etc.) are generally
greater than surface absorbed energy (due to solar radiation, g , and atmospheric
long wave radiation) this leads to a mixed zane at the surface in most cases, al-
though in spring and early summer this mixed zone may be rather shallow.
2.3.4 Inflow and Outflow
To include inflow or outflow in heat transport calculations, two separate
items of information are required, namely the level of entry or discharge, and the
distribution of flow about that point.
It is assumed that inflow enters the reservoir water column at the level at
which its density matches that in the water column. If entrance mixing (see Section
2.3.10) is allowed to take place, the mixed inflow temperature is then taken as
the criterion. Outflow is assumed to be centered about the level of the outlet.
This assumes a linear temperature gradient in the vicinity of the outlet, which may
not be the case. However, Elder and Wunderlich (8 ) point out that even for non-
linear gradients near the outlet, the flow is fairly symmetrical about the outlet
centerline.
2.3.5 Calculations of Inflow Velocity Distribution
The fact that warm inflows flow directly over the surface, and dense inflows
(whether cold or sediment laden) flow along the bottom is accepted in the literature
-10-
-------
-------
on reservoir flows (Howard (n), Goda (10) ). However, the evidence for flows
at intermediate depth is sparse. Elder and Wunderlich (8 ) give the results of
dye tests in Fontana Reservoir (see Figure 2.5). These dye profiles also show
that the inflow velocity profile may be approximated by a Gaussian curve. It is
assumed that 0
where
-u (y) = u (t) e
(2.2)
max
u. (y) = inflow velocity at elevation y
(t) = maximum value of the inflow velocity at time t
max
y. (t) = elevation of inflow at time t
'in
= inflow standard deviation
a.
Little is known about the value of a.- In this model it is held constant
and related to the depth of the entering stream.
u. (t) is found by equating the sum of the inflows into each layer to
max
the total inflow, i.e.,
Q(t) =
B(y) u(y) dy
(2.3)
where
y = surface elevation
s
y = bottom elevation
b
B(y) = width of the reservoir at elevation y, and can be included in the
input data, or as in this case calculated from the horizontal cross-sectional areas
and reservoir lengths, i.e., B(y) = A(y)/L(y). This assumes a rectangular basin.
2-3.6 Calculation of Outflow Velocity Distribution
The outflow velocity distribution is treated in a manner similar to that of
-12-
-------
8-18-66
8-16-66
t-
LU
Lul
LJ
_J
UJ
LU
LJ
L-
LJ
_J
LU
MOO
1500
423'
-" TEMPERATURE F
DYE CONCENTRATION PPB
CONFLUENCE
UITLE TENNESSEE
NANTAHAIA
76
9-16-66
163885
78
80
S-7-66
1643.77
82 80
82
8-31-66
1447.43
P
84 80
84
8-25-66
U49.49
75'--
75'-
1600
1500
1400
68
7068
70 72 70 72 74
LITTLE TENNESSEE RIVER MILES
76 74
76
78
FIGURE 2-5. DYE CONCENTRATION PROFILES IN TRIBUTARY OF FONTANA RESERVOIR
(After Elder and Wunderlich, 1968)
-13-
-------
the inflow, with one important difference; namely, the standard deviation is
calculated using the results of either Koh (16) or Kao (14). The problem of
withdrawal from a stratified reservoir is treated rather extensively in the litera-
ture. The solutions have been summarized by Brooks and Koh (4 ). Huber and Harleman
(12) also have a rather lengthy discussion of Koh and Kao's work.
2.3.6.1 Laboratory Case
Koh's solution(16)for a steady two-dimensional case involving both
viscosity and molecular diffusion is applicable. The thickness of the withdrawal
layer & is given by
(2-4)
where
(eg/av)
x = horizontal distance from the outlet
e = normalized density gradient, ^
P 9y
a = molecular diffusion coefficient
v = kinematic viscosity
For given values of x, g, a and \>, Eq. (2.4) can be put in the form
cons t
6 "" (2.5)
For the laboratory case, using standard values g, v, D,
2
g = 980 cm/sec
2
v = 0.01 cm /sec
a = 0.0014 cm /sec
and x chosen at the mid-point of a horizontal line between the outlet and the
reservoir bottom (x = 240 cm)
K TO'1/6
6 = 2.2 e (2.6)
-14-
-------
Huber and Harleman (12) note that the operating conditions in the laboratory reser-
voir were in the range of Koh's experiments, and that the calculated values of 6
agreed well with those observed during experimental runs. Equation 2.6 was used
in the model verification involving the laboratory reservoir.
2.3.6.2 Field Case
Kao (14) obtained a solution for the case of a dif fusionless flow towards
a line sink at the bottom of the end of a uniform channel. Huber and Harleman (12)
have modified Kao's solution for the unbounded case, and substituted a Gaussian velo-
city profile instead of the uniform one proposed by Kao. They obtained the following
formula
2
.8 f-9-)
\ge/
6 = 4.8 -- (2.7)
\ge/
where q = outflow/unit width, the width being taken as the average width of the
reservoir at the elevation of the outlet. This is very close to the formula pro-
posed by Elder and Wunderlich ( 9 ) , based on some field measurements in Fontana
Reservoir. They suggest
r 2 \
§ = 5.0 I-3-) (2.8)
* ge '
Equation 2.7 has been used in the verification of the model using field
data. The units used for the field case were meters, kgms, days, and for these
units Equation 2.7 can be written as
1/2
S = 0.0092 *-j£ (2.9)
where
5 = thickness of withdrawal layer (m)
2
q = outflow/unit width (m /day)
e. = normalized density gradient (m )
The outflow standard deviation, a , is calculated on the basis that 95% of
the outflow comes from the calculated withdrawal layer. Thus
% - 1756
-15-
-------
and
u (y) - u Ct) a 2°o C2.1D
o o
max
where
u (t) = velocity at y = y = maximum velocity
o J J •'out
max
y = elevation of reservoir outlet centerline.
'out
In the above calculations, the density gradient is determined from the
temperature gradient at the outlet. Early in the warming period, this gradient
may be very small, and may lead to withdrawal over the full depth of the reservoir.
However, in many cases a warmed surface layer may have already developed and it is
unrealistic to expect this lighter water to be drawn down to any extent. To avoid
this, a cut-off gradient has been specified, which limits the thickness of the with-
drawal layer. A further refinement could be included which would lead to an asym-
metrical withdrawal layer, the asymmetry being a function of the temperature profile
on each side of the outlet. This asymmetry was observed in the laboratory reservoir.
However, when this refinement was incorporated in the model, no significant improvement
was noticed and it has been deleted from the final program.
2.3.6.3 Outlet Geometry
The actual outlet geometry is not usually significant. In the model the
outlet appears as a line sink, while in practice it is usually rectangular. However,
measurements in Fontana (7) have shown that at a relatively short distance upstream
of the dam, withdrawal takes place over the full width of the reservoir. Similarly,
as long as the calculated withdrawal layer thickness is large compared to the outlet
height, then the latter may be neglected. If, however, the outlet height is similar
to, or larger than, the withdrawal layer thickness, then a minimum thickness, say
1-2 times the outlet height, should be postulated.
2.3.6.4 Multiple Outlets
In many reservoirs, withdrawal can take place from several levels. This
facility has been included in the mathematical model. The number of outlets, the
level of each outlet and the individual outflow rates, are specified in the input
data. Withdrawal may take place from one or all outlets at any time. The velocity
field of each outlet is calculated as in Section 2.3.6.2, and superimposed on one
-16-
-------
another. No field verification has been carried out for the multiple outlet case,
but the comparison with laboratory results is quite promising. See Figures 2.6, 2.7.
2.3.7 Vertical Advective Velocity
The vertical advective velocity, V, is obtained by requiring that the con-
tinuity condition is satisfied at each level. Starting at the bottom element of
the reservoir, and setting the inflow across the bottom equal to zero, each element
was considered in turn and the vertical velocity calculated so that continuity was
exactly satisfied. See Figure 2.8.
n+1
rv
n+1
n
A 1
n,n-l
V
n
n-1
Figure 2.8 Calculation of Vertical Velocity
A
n,n+l
( VA . + B Ay (u. - u ) } (2.12)
V n n,n-l n J i o /
x n n '
2.3.8 Travel Time for Inflows
In previous field and laboratory studies it was observed that the measured
values of outlet temperature lagged the predicted values. The most obvious cause
is that some time elapses between the time the inflow enters the reservoir and the
-17-
-------
i
^-j
00
Inflow - Outflow Rates
(cm /min)
= Qout
- Outflow
' (elev. 25 cm)
Outflow
elev. (-30.0 cm)
Insolation
2
(cal/cm -min) -
Time (mins)
Measured Outflow -
Temperatures
Predicted Outflow
Temperatures
0 30 60 90 120 150 180 210 240 270 300 330 360 390 420
Time (minutes)
Figure 2.6 Outflow Temperatures. Alternate Withdrawal from Multiple Outlets
-------
to
i
34
33
32
31
30
29
28
QJ 07
S-t Z/
4-1
Cfl
S 26
I 25
24
23
22
21
20
Outflow
El. 20.0cm)
Inflow
/ (Elev. -30,0 cm)
Variable Insolation
Constant Inflow-Outflow Rates
Q^. = 11355 cm3/min
3
= 0 =0 = 3785 cm /min
Outflow'
(Elev. -77.5 cm)
1.0
0.5
0
Insolation (cal/~
cm -min)
0 120 240 360
Time (mins)
—— Measured Outflow-
Temperature
Predicted Outflc*
Temperature
0 30 60 90 120 150 180 210 240 270 300 360 390
Time (minutes)
Figure 2.7 Outflow Temperatures. Simultaneous Withdrawal from Multiple Outlets
-------
time its effect is observed in the measured temperature profile. This effect is
enhanced by the fact that temperature profiles are usually taken near the outlet.
This time lapse has been called the travel, or lag time, and is calculated as follows.
The inflow entering on day N is calculated and labelled QQIN(N). The time lag is
calculated by estimating the velocity of a uniform underflow, or density current,
travelling down the upstream slope of the reservoir. The equations for the velocity
and depth of the density current are given by
, v 1/2 \2/3
u =
and
\23
qj (2-13)
In \I/3
d = 1.92 (V-) (2-14>
\g s/
where
g' = g Ap/p
q = mixed inflow rate per unit width (see Section 2.3.10)
s = bottom slope of reservoir
v = kinematic viscosity
d = thickness of density current
These formulae are obtained using Keulegan's (15) approach, under the
assumption of laminar flow. The velocities and thicknesses predicted for Fontana
were similar to those measured in the field by Elder and Wunderlich (8 ), and hence
it seems reasonable to use these formulae for both the field and laboratory cases.
The density difference, Ap, and the travel distance are calculated as follows.
The entering inflow QQIN(N), temperature TTIN(N) is allowed to mix with the surface
waters. On the basis of the mixed temperature, the final position of the inflow
in the water column is established. The downslope distance to this level, and the
average density difference between the mixed inflow and the surrounding water are
then easily established. Once the time lag LAGTIM(N) is known we have
QIN[N + LAGTIM(N)] = QQIN(N)
TIN[N + LAGTIM(N)] = TTIN(N)
-20-
-------
i.e., an inflow entering the reservoir during time period N is not taken into
account in the heat transport equations until time period JN + LAGTIM(N)J. The
possibility of several days inflow affecting the temperature profiles at the same
time is allowed for.
Alternatively, lag times were taken directly from laboratory observations
of dye traces. A significant improvement in the predicted results followed, showing
that discrepancies between predicted and measured occurrence times can be at least
partly accounted for by this method. In the final analysis, however, it was decided
that the improvement in results was more than, offset by the non-predictive nature
of the modifications. The facility has been kept in the model as part of the main
program, but is usually bypassed.
2.3.9 Variable Surface Level
The surface level is calculated as a function of the initial surface level
and the cumulative inflow and outflow. This was done to avoid the accumulation of
roundoff errors. The continuous variation in the surface level is accounted for by
allowing the thickness of the surface layer to vary between 0.25 Ay and 1.25 Ay
where Ay is the layer thickness. The reason behind the lower limit on the thickness
will become clearer when the overall schematization of the model, and the method of
heat transport into the surface layer is discussed. Briefly, however, an extremely
thin layer (~Ay/100) can lead to the calculation of very high surface temperatures,
and unrealistic long wave radiation and evaporation losses. Within the range spec-
ified above, the model is not sensitive to the layer thickness.
In this model, there is no attempt to include the direct rainfall and eva-
poration into the calculation of surface level. This could be easily incorporated,
if desired.
In the calculation of the surface level, an initial approximation is used,
in that the surface element is considered to have vertical instead of sloping sides.
A correction for the surface element thickness is later included to correct this
approximation. The maximum correction observed was less than 1% of Ay. This pro-
cedure was used to simplify calculations, as the surface area, and hence the volume
of the surface element is a function of the surface element thickness.
2.3.10 Entrance Mixing
When a river enters a reservoir it will entrain some of the reservoir water.
The amount of this entrainment is one of the biggest unknowns in the present knowledge
-21-
-------
of reservoir behavior. It was found that changes in the. amount of entrance mixing
had a very significant effect on the computed temperature profiles and outflow
temperatures. In fact it is possible to bracket the measured laboratory results
with quite reasonable changes in the entrance mixing, as shown in Figures 2.9, 2.10,
2.11 and it therefore seems unreasonable to be concerned with heat transport due
to eddy diffusion until more is known about the entrance mixing effect.
Direct measurements of entrainment in the laboratory reservoir gave the
following results.
Inflow at surface (i.e. warm water) 50% entrainment
Inflow at intermediate depth 200% entrainment
Use of these values in the mathematical model for the laboratory case led
to significant improvements in the correlation between predicted and measured results.
See Figures 2.9 - 2.11.
Modification of the entrance to the laboratory flume reduced entrance mixing
to between 10% and 50%, and this indicates that the geometry of the river entering
the reservoir will probably be very significant.
No field measurements of entrance mixing were available, and hence the
choice of an entrance mixing coefficient is somewhat arbitrary. Work in this labora-
tory on buoyant surface jets (28) indicates that an entrance mixing of 100% is real-
istic. This value gave excellent results (see Figure 2.12) but field studies should
be undertaken.
Entrance mixing is simulated in the mathematical model by withdrawing a
specified amount of water from a selected depth, d , and mixing this water with the
inflow. The amount of water entrained per unit of inflow, is called the mixing
ratio r. The mixed inflow rate Q'. is therefore given by
Q'in = (1 + r) ^in (2.15)
where
Q. = stream inflow rate.
For the Fontana case, with an entrance mixing of 100%, r is equal to unity,
and the mixed inflow rate is twice the stream inflow rate.
-22-
-------
ro
OJ
i
32
30
28
0)
J-l
n)
M
0)
I" 26
4)
H
24
22
T
Laboratory Reservoir
No Insolation
Constant Inflow-Outflow
Measured
• Predicted, r = 0
m
— Predicted, r_ = 0.5,2.0
m
J.
J.
50 100 150 , 200 250 300
Time(minutes)
Figure 2.9 Effect of Entrance Mixing on Outflow Temperatures
350
-------
I
ro
01
M
3
4-1
CO
>-i
0)
OJ
H
30
28
26
24
22
20
18
Insolation
(cal/cm^-min)
T
Laboratory Reservoir
Variable Insolation
Constant Inflow, Outflow"
120 240
Time (minutes)
360
Predicted, r = 0
m
Predicted r = 0.5, 2.0
m
Outflow
I.
_L
60
100
25CT
150 2UD
Time (minutes)
Figure 2.10 Effect of Entrance Mixing on Outflow Temperature
300
350
-------
32 -
ro
ITI
Outflow rate (cm /min)
Inflow rate
Predicted, r = 0.0
Predicted, r = 0.5,2.0
m
0 90 180 270 360
Time (minutes)
Laboratory
Insolation
(Cal/cm^-min)
Insolation
Inflow. Outflow
120 180 240 300
(Time - minutes)
I
22
150
200
250
50 100
Time (minutes)
Figure 2.11 Effect of Entrance Mixing on Outflow Temperature
300
350
400
-------
ro
Crt
i
o
•U
a)
0)
H
4-1
O
24
22
20
18
16
14
12
10
8
Measured
Predicted, r = 0
' m
Predicted, r =1.0
' m
Outflow Temperatures
Through Fontana Dam
1966
Range of Variability
of Measured Values
Jan I Feb I Mar I Apr I May I June I July I Aug I Sept I Oct I Nov I Dec I
50
100
150 200
Time (days)
250
300
350
Figure 2-1Z. Effect of Entrance Mixing on Outflow Temperatures in Fontana Reservoir
-------
The mixed inflow temperature T' is
in
rQ. T + Q T, rT +
ln in £
where T is the average temperature of the water over the depth d and T is the
temperature of the inflowing stream. The outflow velocity from the surface layer
due to the entrained flow is
m TO.
Thus entrance mixing affects both the inflow and outflow velocity profiles
and hence the vertical advection terms. This simple method of allowing for entrance
mixing is probably reasonable for warm water inflows, but for cold inflows it is
probably not realistic to assume that the mixing water comes only from the surface
layer. It is assumed here that d is equal to the depth of the entering stream.
However, the model is not sensitive to moderate variations in d .
m
2.4.1 Derivation of the_Basic Heat Transport Equation
The reservoir is schematized as shown in Figures 2.1 a,b,c. Using the
assumptions and model characteristics discussed in Sections 2.1 - 2.3, a basic
heat transport equation can be obtained from consideration of heat flow through a
control volume bounded by the reservoir sides (Fig. 2.2).
We obtain
8T(y)_ a J_ (A(y) 3T(y)\ _1_ J_ /A/vU/v) \ _ J, J
Jt ~ A(y) 3y V 3y / ~ pcA(y) 3y \A(y)*(y) j A(y) 9y
1
A(y) V"iVJ" "VJ" ^i "o
B(y) T ~ U(y) B(y) T(y))
where
T(y) = temperature at elevation y
V(y) = vertical velocity at elevation y
u_(y) = inflow velocity at elevation y
-27-
-------
u (y) = outflow velocity at elevation y
1 = inflow temperature
A(y) = area at elevation (y)
t = time
a = molecular diffusivity
y = elevation (positive upwards)
2.4.2 Initial and Boundary Conditions
Equation 2.18 requires one initial and two boundary conditions. At the
beginning of spring (t = 0), the reservoir is assumed to be in an isothermal state,
and this provides the initial condition
T = T for all y at t = 0 (2.19)
Conservation of heat at the water surface and the reservoir bottom give
the two boundary conditions.
The detailed formulation of these boundary conditions is given in Sections
3.1.2 and 3.1.3.
-28-
-------
3. Description of the Mathematical Model
3.0 Solution of the Heat Transport Equation
It was not possible to solve Equation 2.18 analytically, subject to the
prescribed initial and boundary conditions, and thus a numerical solution was
sought.
3.0.1 Choice of Numerical Scheme
Three basic finite difference methods are applicable to the convective-
diffusion equation under consideration: (1) an explicit or forward difference
scheme, (2) an implicit or backward difference scheme, (3) a combination of the
first two which results in an implicit method also. In general, explicit schemes
have the advantage in ease of formulation, and often in reduction in computer time.
However, they are subject to certain stability requirements, which can restrict the
size of the time step used in the solution. Implicit methods are unconditionally
stable for any time step, but require more complicated methods of solution.
It will be shown, Section 3.0.3, that the limitations imposed by an explicit
scheme are not unduly restrictive for the reservoir model. Furthermore, this type
of scheme has the advantage of keeping separate the various terms in the heat trans-
port equation, which facilitates alterations and corrections to the mathematical
model.
It is recognized that the explicit scheme lacks the accuracy of a mixed
scheme, such as that proposed by Stone and Brian (24) and used by Huber 0-2). This
lack of accuracy, however, appears only in the treatment of the high frequency har-
monics, and since the temperature profiles tend to be relatively smooth, these har-
monics should not be too important. Furthermore, results of the explicit scheme have
been compared with those from the Stone and Brian scheme, and the differences were
insignificant.
3.0.2 Description of Numerical Scheme
The general mechanism of an explicit finite difference scheme is to find the
value of a variable (e.g. temperature), as a function of space (e.g. depth), at a
time step (n + 1), when the spatial distribution of the variable at time step (n)
is known. The usual approach is to put previously derived differential equations
such as Equation (2.18) in a finite difference form and to obtain approximate solu-
tions for specific input data. Under some conditions, notably when large gradients
exist, either in the dependent variable (e.g. temperature) or in coefficients (e.g.
velocity), many finite difference formulations may lead to significant errors in
-29-
-------
heat or mass conservation. These errors may be avoided if the problem is form-
ulated in terms of a control volume, such as in Figure 2.2. This approach is
used here.
3.0.3 Limitations on Explicit Scheme
The use of an explicit finite difference scheme entails some limitations if
numerical stability is to be maintained. The stability criteria for this case may
be written
< (3.1)
" 2
(Ay)
V < 1 (3.2)
Ay
where
AT = time increment
Ay = distance increment
D = diffusion coefficient
V = magnitude of velocity in y direction
As long as turbulent diffusion is neglected, Equation (3.1) is not at all
restrictive. However, Equation (3.2) can lead to a rather small At because of
inflow conditions which usually occur only a small percentage of the time. This
was remedied by use of a variable time step At. See Section 4.5.
3.1 Formulation of Explicit Finite Difference Scheme
This scheme is formulated by applying the heat balance approach to a typical
reservoir control volume or element. The elements used were horizontal sections of
the reservoir, bounded by the reservoir sides (see Figure 2.2). An internal element
will be treated first, followed by surface and bottom elements.
3.1.1 Internal Element
The jth internal element for both the laboratory and field case is rectangular
in plan, average length L average width B., average area A. and height Ay. In
_. .. . ... ^ J IT
_3T
3n
normal to the side. In the laboratory case, however, heat losses through the sides
the field case the sides are considered to be insulated, i.e., — = 0 where n =
3n
-30-
-------
are accounted for.
The rate of temperature change within the jth element (see Figure 3.1) is
equal to the difference between the heat flow rate into and out from the element.
The heat flow is made up of five terms, the direct absorption term, the diffusion
term, the vertical advection term, the horizontal advection term, and the side
loss term (laboratory case only).
— j + 1
Figure 3.1 Internal Element
3.1.1.1 Direct Absorption
The temperature change due to the direct absorption of transmitted radia-
tion is given by
1 f
AT, = T~T~ ' U--LI • A-J_I • ~ • • i A- ^
1 pcA.Ay V J+1,J J+ljJ J^J""1 J»J
J ^
.At
(3.3)
The transmission of radiation has been discussed in Section 2.3.2.
-31-
-------
3.1.1.2 Diffusion Term
"2 ' A^ I -*%-*- Vl,J--V^ V-'>- 4t (3'4>
\
. Y
,3-J/
a = molecular diffusion coefficient.
3.1.1.3 Vertical Advection Term
The vertical advection terms are of particular interest in this problem
due to the configuration of the vertical velocity field. Depending on the elevation
of the inflow, vertical velocities may be all positive (upwards), all negative, or
positive at some elevations and negative at others. Thus the formulation of the
vertical advection term is not straightforward. Four separate combinations of
vertical velocities are possible and each combination requires a separate formula-
tion, or serious errors result (Bella (3) ).
Taking upwards as the positive direction and noting that V. is defined at
the bottom of the jth element we have
(a) V,,
Aip _ •*- (VTA VTA. I A t- ("3 "^ \
3 A A-(7 I V^-L-!_1A-i -?_1 V-iJ_1 •"•-!•"•-5.L.1 .! ) • At (,J-->)
(b)
AT = f V T A — T7 T A IA«- / o £\
3 A.AV I 11 i.1-1 vi+lii+TAi+l.i J0, V <0
AT
3 - A^7 I Vi-lVi-1 ' VlVlVl < Kt (3.7)
-32-
-------
(d) V.< 0, V. >0
^7— ( V.T.A. . . - V....T.A.,.. , \At (3.8)
•jAy V^ j j j,j-l j+1 3 3+1 »J )
AT3 * A
3.1.1.4 Horizontal Advection Term
AT/ = -T-T- (u- T- ~ u T-) • B-, Ay .At (3.9)
4 A.Ay 1.1 o.i j
3 J 3
3.1.1.5 Side Heat Loss Term - Laboratory Case Only
Huber and Harleman (12) assume that the outside wall of the plexiglas
flume is at the same temperature as the inside and that radiation is the prime mech-
anism for heat loss. Both the reservoir sides and the laboratory surroundings are
assumed to radiate almost as black bodies, leading to a net radiation loss rate for
the sides of the jth element.
b = 0.97 k fT.4 - T 4) (3.10)
pm \ j a /
where
T = room temperature
T., T are in °K
J a
k = Stefan Boltzmann constant
This side heat loss term does not appear to have an important effect on
the calculated temperature profiles, and Huber's treatment is deemed sufficient.
The side heat loss term is
AT, = r-r- [<(> (L. + B.).2Ay].At (3.11)
5 pcA.Ay Tmx J J
AT =__i— [2 (L. + B.)] .At (3.12)
-33-
-------
3.1.1.6 Total Temperature Change
For the case of the field reservoir the total temperature change
in time increment At, is given by
AT
(3.13)
For the special case of the laboratory reservoir the extra side heat loss term
must be added to the above Equation (3.13).
3.1.2 Surface Element
The surface boundary condition is obtained in a similar manner to the equa-
tion for the internal element, namely, by writing a heat balance equation for a
specified control volume. In this case, however, the thickness of the control volume
varies, and extra heat sources and sinks have to be considered. (See Figure 3.2)
jm
LU I Jm-l
Figure 3.2 Surface Element
-34-
-------
The new surface level is calculated before any heat flow is considered
Once Ay is known, velocities are calculated such that
S Hi
u. B. Ay + V A = u B AY (3 14)
im > 'sur vjin jm,jm-l V jm sur ^'^}
An average surface area S is also calculated such that
' area '
Sarea = (AysUr
3.1.2.1 Direct Absorption Term
AT = — [ (j, A - ((,. . A. . .. 1 .At (3.16)
1 pcSarea Ysur \ ySUr ysur J^^1""1 jm,;jm-l
where
3.1.2.2 Diffusion Term
/ (T. -T. ). \
( •1" •^-1 A. . J .
\ AY jm,jm-l/
AT2 = - S Z? I A,"' A, , j .At (3.18)
area sur
where
a = molecular diffusion coefficient.
3.1.2.3 Vertical and Horizontal Advection Terms
V >0
A rn -^ f |\7T A -4-11 RAV T
3,4 S Ay V jm jm-1 jm,jm-l i. jm sur i
' area sur \ J J J 'J im J
u B. Ay T. ^ .At (3.19)
o. im sur jm '
jm
-35-
-------
Substituting from Equation (3.14)
S y jm-l jm- i .
area sur -J -> >-J J J
V. < 0
Jm
" (V'mA-m •m-lT'm + Ui B'mTiZ
area sur x ' Jm
- u B. T. Ay A -At (3-21)
o. im im sur I
jm J J I
Substituting from Equation (3.14) the V. term disappears and
— (u B. (T. - T. )).At (3.22)
* jm X m 7
area
3.1.2.4 jiide Loss Term - Laboratory Reservoir Only
This term is independent of layer thickness, and is therefore given by
Equation (3.12)
3.1.2.5 Surface Absorbed Solar Radiation Term
In Section 2.3.2 it was pointed out that a fraction, g, of the solar
radiation can be considered to be absorbed at the surface
AT6 = ^S Ky— I ^o Aysur; ' At (3-23)
area sur '
3.1.2.6 Surface Heat Loss Term
All heat losses from the field reservoir, apart from advective losses,
take place from the reservoir surface. The total heat loss rate per unit area,
can be written as
ra
-36-
-------
where
cj) = rate of heat loss due to evaporation
(j> = rate of heat loss due to conduction
= rate of heat loss due to net longwave radiation
IT 3.
See Chapter 4 for further discussion of these parameters.
The surface heat loss term is given by
AT.
1
pcS Ay
area sur
ysur
• At
(3.25)
3.1.2.7 Total Temperature Change
For the case of the field reservoir, the total temperature change in the
surface element AT , in time increment At, is given by
s
AT = (AT + AT0 + AT . + AT, + AT )
S ± / -J j 4 D /
(3.26)
j=jm
For the laboratory case the side loss term (AT ) must also be included.
j=jm
3.1.3 Bottom Element
Ay/2
• _ O
— 5=1
Figure 3.3 Bottom Element
-37-
-------
The bottom element differs from the internal elements in two ways. Firstly,
it is half as thick, so that the point j = 1 coincides with the bottom. Secondly,
there is no heat flow through the lower interface. The average area is defined as
(A2jl+A1)/2 (3-27)
3.1.3.1 Direct Absorption Term
AT1 -
3.1.3.2 Diffusion Term
/ T — T \
i— ( a ~ — ' A_ , ) .
^- \ ^ 2'1/
X
AT = — — a ~ — ' A_ ) .At (3.29)
2
3.1.3.3 Vertical and Horizontal Advection Terms
AT., = ——
3.4 A t-
But V is defined by
V2A2,1 = Blf
AT3,4 - ' VBl Ti-Tl 'At (3.32)
-38-
-------
For V < 0
AT = — I 11 R y IT — T \ _ \7 A IT
AJ-o A A .— /o I u_. JJi o IJ-. IT I — VIAO i I •'-'-
.,12\i LI 22,1 \
-------
where y . = the elevation of the bottom of the mixed layer
•* Tn~i ~v
mix
T = T (t) = mixed temperature at time t
This procedure is carried out in subroutine AVER. The surface layer is
examined. If the temperature is lower than that of the next layer, the top two
layers are mixed. The mixed temperature is compared with that of the third layer,
and if necessary, the top three layers are mixed. This process continues until a
stable profile is obtained.
3.3 Energy Check
Since the surface losses are a function of the surface temperature, it is
possible for large errors to occur in the overall energy balance without this being
evident from the calculated profiles, since any error tending to increase or decrease
the surface temperature, causes an increase or decrease, respectively, in the heat
loss from the surface, and thus a compensating effect is introduced. For this reason
a check is made on the energy balance for each time period. The most sensitive check
is the ratio between the increase in stored energy and the net energy inflow. How-
ever, when the net energy inflow approaches zero, this ratio can be misleading, and
for this reason, a second ratio is given, which compares the total energy stored
plus the energy outflow, with the initial energy plus the energy inflow.
3.4 Numerical Dispersion
One of the drawbacks of numerical schemes involving advection terms, is
that inaccuracies tend to be introduced. These take the form of an additional diffu-
sion or dispersion effect, and are called numerical dispersion.. A detailed treatment
of this effect is given by Bella (3 ). For the one-dimensional variable area case,
the numerical dispersion coefficient can be written as
Dn = -i- (Ay - V^ ^- . At J (3.37)
j
where
A. = horizontal cross section area at center of jth element
A. . , = horizontal cross section area at interface between jth and
(j-l)th element, see Figure 3.1.
-40-
-------
For most reservoir cases the one-dimensional constant area form is suffi-
cient i.e.
D = -J- (Ay - V. At) (3.38)
nj ^ 3
By differentiating Equation (3.38) with respect to V, it can be shown that D
max
occurs for V = Ay/2At. The maximum numerical dispersion is therfore given by
.2
(3.39)
tJi_l L.
max
For the field case At = 1 day. Ay = 2 meters, this results in D equal to 40 times
max
the molecular diffusion coefficient, which is still a relatively small value. For
the laboratory case, Ay = 2 cm, At = 1 minute, this results in D equal to ~7 times
max
the molecular value.
Since this model uses a very small diffusion coefficient, it is not possible
to minimize the effect of numerical dispersion by replacing a with (a - D ) in the
diffusion term. However, it seems reasonable to use the opposite procedure, namely
replacing a by (a + D ), thus doubling the numerical dispersion, to see if this
error has a significant effect. From Figure 3.4 it can be seen that the effect of
doubling the error is a minor one, and hence it seems likely that numerical disper-
sion is not significant in this case.
-41-
-------
24
0)
M
D
22
20
18
16
14
• 2 12
? i
§ 10
H
O
H
O
2
0
Jan
Measured
Predicted, r =1.0
> m
Predicted, r = 1.0
' m
Numerical Error Doubled
Range of. Variability
of Measured Values
Outflow Temperatures
Through Fontana Reservoir"
1966
Feb I Mar I Apr I May ' June I July I Aug I Sept I Oct ' Nov I Dec '
50
100
250
300
150 200
Time (days)
Figure 3.4 Effect of Numerical Dispersion on Outflow Temperatures
350
-------
4. Determination of Parameters
The principal heat budget parameters, used in the model, are net solar
radiation, net longwave radiation, evaporation and conduction. These parameters
are the same as those used by Huber, and for further details the reader should
consult (12). A short discussion of each parameter is given below.
4.1 Net Solar Radiation (cj> )
Several formulae are available for obtaining the solar radiation $ in terms
2 s
of the solar constant (1-98 cal/cm /min) , solar angle, cloud cover and other atmos-
pheric factors. These formulae were examined by Andersen (U.S.G.S. - 1954) as part
of the Lake Hefner studies, and he concluded that indirect methods of evaluating
solar radiation are accurate only within about 15%, and therefore measured values
of $ should be used whenever available. A reasonable value for albedo is 7% (l ).
S
Values of § =
y
The net longwave or back radiation is defined as the difference between
the longwave radiation emitted by the water, and that emitted by the atmosphere,
and absorbed by the water.
4.2.1 Longwave Radiation from Water
The longwave radiation emitted by the water, , is known accurately (within
+ 1/2% (1) ), and is given by
d, = 0.97 k T 4 (4.1)
Yrw w
where T = water surface temperature in °K
w -624
k = Stefan Bolzmann constant = 1.171 10 kcal/m - °K - day
The values of d> are calculated by the mathematical model as a function
Yrw
of the surface temperature(T °K) generated witnin the model.
4.2.2 Atmospheric Longwave Radiation cj>
3.
If available, measurements of atmospheric radiation should be used, and values
of may be read in as input data.
3,
-43-
-------
If measurements are not available, reasonable values of a may be calculated
as a function of air temperature and cloud cover. The best approach appears to be
that of Wunderlich, who uses Swinbank's (27) clear sky formula, modified for cloud-
iness.
* = 0.937 . 10~5 k T 6 (1..0 + 0.17 c2) <4'2>
a a
where
c = fraction of sky covered by clouds
T = air temperature °K, 2 meters above the water surface.
3.
Thus for the field case, assuming 37, reflection of by the water surface
3>~
$ = 0.97k [T 4 - 0.937 T & (1.0 + 0.17 c2)] (4.3)
The option of calculating the net longwave back radiation
-------
= relative humidity
The heat loss is given by
where L = latent heat of vaporization
= 595.9 - 0.54 T
T = surface temperature in °C
S
(4.7)
Conduction losses are usually related to the mass flux (E ) by the Bowen ratio R
m
where
(T - T )
(4.8)
where
Thus,
N = constant
> = p(a + bW) (e
c s
f
+ N
_
s
•u)
(4.9)
Note - T is in °C
a
4.3.1 Field Evaporation
The model contains the option of using the Rohwer (21) or the Kohler (17)
formula, both of which are based on Equation 4.9. In the field case used to verify
the model (Fontana Reservoir) best results were obtained with the Rohwer formula,
and thus this formula is recommended.
4.3.1.1 Rohwer Formula
2
The value of (cj> + ) is in kcal/m - day and the units are as follows;
3 e c
p in kg/m ,
a = 0.000308 m/day-mm Hg,
b = 0.000185 sec/day-mm Hg,
W in m/sec (measured six inches above the surface),
e in mm Hg,
S
-45-
-------
e in mm Hg ,
a
4> Is expressed as a fraction,
L in kcal/kg,
c in kcal/kg- °"C,
N =269,1 kcal-mm Hg/kg-°C
No specification of the elevation for the measurement of T& is given.
Rohwer's formula can thus be expressed in the form of Equation 4.9 as
Rohwer - Field / 269.1(1^) \
*e + ^c = (0.000308 + 0.000185W) p (eg- i(;ea> f L+cTg + - ^ j - J (
N s a
4.3.1.2 Kohler Formula
This field evaporation formula is due to Kohler (1954) as given by
2
Wunderlich. The quantity $ + 4> will have units of kcal/m -day where
... 3
p in kg/m ,
a = 0,
b = 0.000135 sec/day-mb,
W in m/sec (not less than 0.05 m/sec),
W is measured two meters above the surface,
e and e in mb (milli-bars) ,
S cl
is expressed as a fraction,
L in kcal/kg,
c in kcal/kg-°C
N = 372 kcal-mb/kg-°C
T is measured two meters above the surface.
a
The Kohler formula is expressed in the form of Equation 4.9 as
/ 372(T - T )\
fl + cTg + ° ) )
N s v '
Kohler - Field / 372(T - T )
- 0.00
-------
value of the constant (a) was found to depend on the mode of surface heating.
Values of (a) were obtained from laboratory measurements and applied in Rohwer's
formula.
4.3.3 Negative Evaporation
In some cases the vapor pressure gradient given by the term (e - i|>e )
S 3.
is negative which leads to a negative evaporation loss. While it is known that this
phenomenon occurs in nature, in the formation of dew, nevertheless little is known
about the process when it occurs over a water surface, and the constants in Equations
4.10, 4.11, apply only in the case of a positive vapor pressure gradient. For this
reason, when the evaporation is found to be negative, it is put equal to zero. The
conduction term is kept, however, although since it was initially obtained as a func-
tion of evaporation, the accuracy of the formula used is open to question.
4.4 Required Input Data
For best results the following parameters should be measured in the field:
1. Solar radiation
2. Atmospheric radiation
3. Air temperatures
4. Relative humidities
5. Wind speeds
6. Streamflow rate
7. Streamflow temperature
If necessary, solar radiation and atmospheric radiation can be calculated,
and in this case the cloud cover should be known. Inflow rates and temperatures
can be obtained from Streamflow records. The extinction coefficient of the water (n )
can be obtained either from Secchi disk, or underwater photometer readings taken
in the lake, reservoir, (or river, if the reservoir is not yet built). However some
change in n may be expected when the reservoir is built. It should be noted that the
required input data is such that it would be available if a reservoir were planned.
4.5 Selection of Time and Distance Steps
An important step in the use of the model is the choice of the distance incre
ment Ay, and the time step At. Both Ay and At are restricted by the two stability
criteria
' D -^-j < \ (4.12)
(Ayr
-47-
-------
and
v At < (4.13)
Ay -
where
D = diffusion coefficient
V = vertical velocity
Other factors such as the input data, and depth of reservoir may also be
important. A reasonable approach would be to choose Ay first, based on the depth of
the reservoir. A minimum of twenty (20) increments is recommended. Approximately
fifty (50) increments have been used in the model verification runs. Too small a
Ay leads to very costly runs and the possibility of instabilities near the surface.
See Section 4.5.4. Once Ay is chosen, At may be obtained from the stability
UlclX
criteria. A value of V for use in Equation (4.13) may be taken as
max
V
max A
out
where
Q = outlet discharge
A = horizontal area of reservoir at outlet
out
Once At is known, a reasonable value of At, based on the input data, may be chosen.
in 3.x
For a field reservoir, a At of one day is a reasonable choice.
4.5.1 Variable Time Step
Under certain inflow conditions, in particular high inflow rates and low
inflow temperatures, V may be considerably larger than that given by Equation (4.14)
TTlcLX
and thus the stability criterion (4.13) may be violated. These inflow conditions may
exist for a very small percentage of the annual cycle, and it is unrealistic to base
the choice of At on them. This is remedied by reducing At during the critical period.
For each time period the vertical velocity profile is scanned and the maximum velocity
(in either direction) is selected. The criterion (4.13) is now considered and if
violated, a new At is selected such that the smallest integer number of new At's is
equal to the previous time interval. In the cases considered the extra time involved
-48-
-------
was insignificant, but in some reservoirs with a high inflow/capacity ratio this
could lead to very expensive runs. However, a basic assumption of the model, namely
that isotherms are planar and horizontal, restricts the model to reservoirs with a
low inflow/capacity ratio.
4.5.2 Lower Limit on Ay in Field Reservoir
It has been found that if the surface layer is much less than 0.5 meters,
oscillations may occur in the surface temperature under some conditions. This is
due to the fact that when there is a net positive heat flux into a thin surface
layer, large temperature increases result. This is followed in the next time step
by unrealistically large surface heat losses, which lower the surface temperature,
thus increasing the likelihood of a large positive heat flux in the next time step.
Large oscillations in the surface temperature may result. The usual case of a net
negative heat flux into the surface layer is not affected by the layer thickness
as mixing with the lower layers occurs, and the layer thickness is accounted for
here. To eliminate these oscillations, a lower limit of 0.5 meters is applied to
the surface layer thickness.
4.6 Possible Program Changes
Apart from the redefining of the minimum surface layer thickness to allow
for a smaller Ay, other changes may be found necessary when the program is applied
to different field cases.
4.6.1 Cut-Off Gradient
As discussed in Section 2.3.6.2, to avoid the withdrawal of warm surface
waters when this would be physically unrealistic, the width of the withdrawal layer
is limited by means of a cut-off gradient. The choice of the limiting gradient is
somewhat arbitrary, and if unrealistic outflow temperatures are obtained, this value
may have to be changed.
4.6.2 Entrance Mixing
The present method of allowing for dilution of the streamflow as it enters
the reservoir is considered unsatisfactory, and should be changed as soon as informa-
tion about this process becomes available.
-49-
-------
5. Criteria for Applicability of Model
The assumption of horizontal isotherms is basic to the mathematical model
developed in this study. The following criteria may be used to determine whether
this assumption is valid.
5.1 Wind Mixing Criterion for Lakes and Reservoirs
Shallow lakes and reservoirs in exposed conditions tend to have tilted
isotherms during periods of high wind velocity. Since this effect is only a temp-
orary one, it can probably be neglected, except in the case where the thermocline
reaches the surface. In this case, the hypolimnion is mixed with the epilimnion
by direct wind action, and Mortimer (18) considers that it is this mechanism which
accounts for the increases in depth of the epilimnion after a windy spell. The
inclination of the thermocline (here defined as the bottom of the surface mixed
layer) can be obtained as follows.
Sverdrup (25) gives the inclination (i ) of the surface due to wind effects
s
as
-7 TT2
i - 4 x 10 ' . Z- (5.1)
s h
where
h = average depth of water (m)
W = wind speed (m/sec) (use mean of the maximum wind speed)
The inclination of the thermocline i is given by
i = i -^-
t s Ap
The maximum displacement d of the thermocline by the wind is therefore
where
Lfc = average length of lake in the thermocline region, in the direction
of the wind (m) .
-50-
-------
w
Figure 5.1 Effect of Wind on Thermocline
The criterion for use of the model is that the mixed depth (h ), calculated
m
in the model for midsummer, should be greater than the thermocline displacement (d )
for mean maximum wind conditions, i.e.,
4.10 7W2
m
Ap
(5.3)
Ap = density difference between epilimnion and hypolimnion
If the inequality (5.3) is satisfied, the mathematical model is applicable
in the case of a lake. For a reservoir, however, the following criterion must also
be satisfied.
5.2 Flow Through Criterion for Reservoirs
When the flow through the reservoir is large in comparison to the reservoir
capacity, the best criterion is that of Orlob (20), who introduces a densimetric
Froude number which is a function of the flowthrough rate, the reservoir geometry
-51-
-------
and a density difference. The criterion is
-
-------
6. Verification of Mathematical Model for Field Reservoirs
The mathematical model was verified for the field case using data from
Fontana Reservoir. During 1966 the TVA Engineering Laboratory at Norris, Tennessee,
conducted an extensive field study of the reservoir in which all pertinent hydrolo-
gical and meteorological parameters were measured. These data are among the best
presently available with which to compare predicted and measured temperature dis-
tributions for the field reservoir situation.
6.1 Description of Reservoir
Fontana Reservoir, shown in Figure 7.1, is located on the Little Tennessee
River in western North Carolina. It is approximately twenty-nine (29) miles long
and is fed by three rives, the Little Tennessee, the Tuckasegee, and Nantahala rivers
as well as many smaller streams. The reservoir is highly stratified during the summer.
6.2 Input Data
All data used in the program are listed in Section 7.9 in the form of computer
input. Certain data were available on an hourly basis and other data on a daily basis.
The program was run with time steps of one day and all hourly data were reduced to
daily averages.
Air temperatures, relative humidities and wind speeds were available on an
hourly basis. The latter were measured at a reservoir shore location and transformed
to mid-lake values by an empirical correlation provided by the TVA Engineering Labora-
tory.
Inflow data was available on a daily basis for the three tributaries listed
above and for runoff from the watersheds bordering the north and south shorelines.
The sum of these sources is used as the reservoir inflow, and the inflow temperature
was taken as the weighted average of the inflow temperatures from each of the five
sources. The outflow rate and temperature were available on a daily basis.
The solar radiation values were computed by Huber (12) and correlated for
a limited number of cases with measured values obtained by TVA. There is some
doubt as to the accuracy of these values. See Huber (12) for details. Values of
atmospheric radiation were computed by Equation 4.2 for each hour and averaged on
a daily basis.
Geometric data, namely areas and lengths of the reservoir were available at
fifty foot intervals. These values are used as input data, and values at intermediate
elevations are found by linear interpolation.
The values of the absorption coefficient, n> and surface absorption fraction,
-53-
-------
3, were obtained from photometer measurements in Fontana Reservoir (29}•
6.3 Results
The comparison of predicted versus measured values for outflow temperature
as a function of time is given in Figure 6.2. Predicted versus measured temperature
profiles are given in Figures 6.3 - 6.6.
The measured temperature data was taken during the 1966 study (29).
The agreement between predicted and measured values for both outflow temperature
and temperature profiles is very good. As noted previously, the effect of entrance
mixing is significant, at least towards the end of the yearly cycle.
6.4 Verification from Other Sources
A modified version of this model has been checked against data from Lake
Norman in North Carolina with good results (23). The modifications were necessary
because of a thermal power station which uses Lake Norman as a source of cooling
water.
A simplified version of the model has been checked against data from Lake
Tahoe, and Castle Lake, California, and from Manton Reservoir in Northern Australia
(22). Good agreement was obtained in all cases, but the quality of the available
meteorological and hydrological data was much lower than for Fontana Reservoir, and
these results are not presented here.
-54-
-------
River
Fontana
Dam
Nantahala River
Little Tennessee River
Figure 6.1 Map of Fontana Reservoir and Watershed
-55-
-------
24
22
20
18
16
u 14
§ 12
a; 10
C«
oi
2
0
T
Outflow Temperatures
Through Fontana Dam
1966
Measured
Predicted
Range of Variability
of Measured Values
Jan I Feb I Mar ' Apr ' May I June I July I Aug I Sept I Oct I Nov I Dec
50
100
150
200
250
300
350
Time (days)
Figure 6.2 Measured and Predicted Outflow Temperatures Through Fontana Dam
-------
I
tn
0)
4-)
0)
E
C
o
•H
4-1
ca
cu
520
510
500
490
480
470
460
450
440
430
420
410
Fontana Reservoir
June 22, 1966
outlet C.L.
Measured
Predicted r =1.0
m
I
I
I
I
I
I
02 4 68 10 12 14 16 18 20 22 24 26 28 30
Temperature (°C)
Figure 6.3 Measured and Predicted Temperature Profiles
-------
Ul
CO
I
0)
4-1
0)
e
c
o
•H
+J
tfl
W
510
500
490
480
470
460
450
440
430
420
410
I I I T
Fontana Reservoir
July 20, 1966
outlet C.L.
Measured
Predicted
I
I
I
I
I
4 6 8 10 12 14 16 18 20 22 24 26 28
Temperature (°C)
Figure 6.4 Measured and Predicted Temperature Profiles
30
-------
or
10
^
0)
510
500
490
480
470
460
2J 450
0)
440
430
420
410
I I I
0 2
Fontana Reservoir
September 15, 1966
Measured
Predicted, r = 1.0
m
Outlet C.L.
I
I
8 10 12 14 16 18
Temperature (°C)
I
20 22 24 26 28 30
Figure 6.5 Measured and Predicted Temperature Profiles
-------
en
o
i
510
500
490
480
« 470
M
cu
§ 460
c
o
•H
0)
iH
W
450
440
430
420
410
Fontana Reservoir
November 10, 1966
^ Outlet C.L.
Measured
Predicted
I I
I
10 12 14 16 18 20 22 24 26 28 30
Temperature (°C)
Figure 6.6 Measured and Predicted Tem|>eratuxe Profiles
-------
7. The Computer Program
The basic logic of the computer program is given by the following flow
charts. However, little attempt is made to reproduce details, particularly
where these have already been discussed.
7.1 Main Program
This program contains the input and output commands, initializes variables,
and carries out the heat flow calculations.
-61-
-------
[ Start J
Read Input
Data
Calc.
1 (Laboratory,
Calc. Labora-
tory Long Wave
Radiation
(Field)
Read Extra
iInput Data
J
Calc. Elevations, Len-
gths, Widths, Areas at
Center of each Volume
Element
Calculate Elevations,
Lengths, Widths, Areas,
at Center of Each Vol-
ume Element
Calc. Max.
Surface Level
Write Out
Specified
Data
Initialize
Variables
[Continue
-62-
-------
Calc. Initial Surface
Area, Thickness of Sur-
face Element, Average
Area of Surface Element
Initialize Variables
for Energy Check, Calc
Initial Stored Energy
(EO)
Calc. JP (used
in output formal
Set All Velo
cities Equal
to Zero, a
1.0 °
(Start of
Iteration Loop
of Program
1
t =
N =
t +
N +
At
1
(Inflows and Outflows Not Equal to Zero)
No Entrance
1 Mixing
Calc. Tempera-
ture of Mixing
Water
-63-
-------
Calc. Temperature of
Mixed Inflow
Calc. Inflow Level
Travel Time
Included
1 No Travel Time
Calculation of
Travel Time At
Calculation of Inflow,
Temperature at Time
(t + At)
Write Velocity of Inflow
in Reservoir,
Qln(t). Att, T.(t)
Relocate Level of Inflo
ontinue
-64-
-------
Surface Level Constant with Time
Calculate Surface Level
Qm), Thickness of Sur-
face Element Ay
sur
Calculate Surface
Area, Average Area of
Surface Element
Check Accuracy of Surface
Level, and Correct Thickness
of Surface Element
Ay = Ay + Ay
sur 'sur J
Set Temperatures for New
Surface Layer Equal to
Temperature of Old Surface
Layer, if Necessary
Compute JP
Continue
-65-
-------
Call SUBROUTINE SPEED
Calc. At' so that „ At' n
v ——< i
I = integer
f Start of Heat
I Transport
V Calculations
M = 0
M = M + 1
Calc. AT for Internal Layer
for all J
a) Direct Absorption
b) Vertical Advection
(4 possibilities)
c) Horizontal Advection
d) Diffusion
-66-
-------
Calc. AT for Surface
Layer
Call FUNCTION FLXOUT
for Calc. of Surface
Heat Losses
Calculate AT
for Bottom
Half Layer
T = T. + AT.
J2 ^ 3
(Field - No^v
Side Losses) 1*
Laboratory
1 Side Losses
Calculate Side
Losses AT
T = T - AT
Sj
Check Reasonableness of
Results. If T. > 100°C,
jm
then Set t = t _
stop
Call SUBROUTINE AVER
-67-
-------
Calculate Outflow Tempera-
tures for each Outlet.
Call SUBROUTINE TOUT
Energy Check
(Yes, Print OvJ
is Required)
No. (Print Out Not Required)
Write Specified
Output
Return to Start of
Main Iterative Do-
Loop
-68-
-------
7.2 Function FLXOUT (N)
Calculation of Surface Losses L due to Evaporation, Conduction, and
Radiation. If Evaporation Losses are Negative, Assume Zero.
f Start J
Calculate Air Temperature
(T ) at Time t, Relative
3.
Humidity O) at Time t
Water Surface
Temperature
T = T
s jm
Calculate Latent
Heat of Vaporization
Calculate Vapor
Pressure at Water
Surface and in Air
Field 2
1 Laboratory
Calculate Evaporation Plus
Conduction Losses, Radia-
tion Losses, Sum Losses
Calculate Wind Speed
for Time t
[ Return J
Continue]
-69-
-------
Calc. Cloud Cover at Time
t. Calc. Atmospheric
Radiation from Eqn. 4.3
Calc. Evaporation
Using Kohler Formula
Calc. Total Losses
Calc. Atmospheric
Radiation at Time T,
from Read in Values
Calculate Evaporation
Using Rohwer Formula
Calc. Total Losses
+4> +
L Yra ye yc
9
( Return J
1 '
I Return )
-70-
-------
7.3 Subroutine TOUT (HEATOT, FLOWOT)
Compute Outflow Temperature
I Start )
I = Outlet Number
Sum Outflow to Outlet (I)
Multiplied by Temperature
for each Layer (HEATOT)
Sum Outflow to Out-
let (I) for each
Layer (FLOWOT)
< P
f Return J
-71-
-------
7.4 Subroutine SPEED(N). Calculate Velocities
( Start J
Compute Normalized
Inflow Velocity
Profile
Compute Maximum Inflow
Velocity u , Hence u
max
Calculate Temperature
Gradient at Outlet (I)
(DERIV)
Calculate Density Grad-
ient
Calculate Withdrawal
Layer Thickness, 5,
Using Kao's Modified
Formula
Calculate
€>
For Very Small Gradients, Put
Withdrawal Layer Thickness =
to Twice the Distance Between
the Outlet and a Specified
Temp. Gradient. Calculate 0,
If Specified Temp. Gradient °
Does not Exist, Put a =
100 Ay. °
Calculate Withdrawal
Layer Thickness, 5,
Using Koh's Formula
-72-
Continue)
-------
Compute Normalized
Velocity Profile for
Outlet (I)
Calculate Max. Out-
flow Velocity uo (I)
TT / -w \ TT13.X
Hence u (I)
j
Repeat for all (I)
Sum Outflow Velocities at
each Level, Including Those
Due to Entrance Mixing,
Thus Obtaining u
j
Calculate Vertical
Velocities from
Continuity
i'
I Return J
-73-
-------
7.5 SUBROUTINE AVER
Convective Mixing in Region of Instabilities. Complete Profile Checked
Allow Convective Mix-
ing to Occur Until In-
stability Disappears
j = j - 1
Calculate New Tempera-
tures for Bottom Mixed
Region
Calculate New Tempera-
tures for Surface Mixec
Region
Calculate Mixing Depth
h
m
-74-
-------
7-6 Remaining Subprograms
FUNCTION TTIN(N), FUNCTION FLXIN(N), FUNCTION QQIN(N) and FUNCTION QOUT(N,I)
all calculate the appropriate values of Inflow temperature, solar radiation, inflow
rate and outflow rate at time t = NAt. These values are obtained by linear inter-
polation from the input data.
7.7 Input Variables
Sample input to the program is provided in Section 7.9 to clarify the
input formats. Note that surface elevations, apart from an initial elevation, are
not essential, and can be ignored, if necessary.
In the following, input variables are grouped according to the cards on
which they appear.
Where options are available, asterisks are put in front of the option recom-
mended for the field reservoir. Where constants are used, recommended values of the
constant for the field reservoir are given.
Card 1, FORMAT 20A4
WH = Alphanumeric variable used to print a title at
beginning of output. Anything printed on this
card will appear as the first line of output.
Card 2, FORMAT 20A4
WH = Alphanumeric variable used to list units used in
computation prior to output at each time step.
Card 3, FORMAT 1615
JM = Initial number of grid points = number of the
surface grid point.
** KATRAD = 1 Atmospheric radiation is read in.
= 2 Atmospheric radiation is calculated in the program.
KSUR = 1 for a constant surface elevation.
** =2 for a variable surface elevation.
KOH = 1 for use of Koh's Equation 2.6 or 2.7, for computing
the withdrawal thickness.
** =2 for use of Kao's Equation 2.10.
** KQ = 1 for computations with inflow and outflow.
= 2 for computation with no inflow or outflow.
-75-
-------
KLOSS
NPRINT
KAREA
KMIX
MIXED
KTRAV
Card 4T
YSUR
DY
DT
TSTOP
TZERO
EVPCON
Card 5.
SPREAD
SIGMAI
ETA
BETA
RHO
HCAP
DELCON
= 1 for Rohwer laboratory evaporation formula.
= 2 for Kohler field evaporation formula
(Equation 4.11).
= 3 for Rohwer field evaporation formula (Equation 4.10).
= Number of time steps between print outs of calculations.
= 1 for laboratory reservoir calculations.
= 2 for calculations for any other reservoir.
= 1 far no entrance mixing.
= 2 to include entrance mixing.
= Number of grid spaces in surface layer for
entrance mixing =4.
= 1 travel time of water within reservoir is neglected.
= 2 "travel time of water within reservoir is accounted for.
FORMAT 8F10.5
= Surface elevation at beginning of calculations.
= Distance increment Ay.
= Time step, At.
= Time at which program ceases calculations.
= Initial isothermal reservoir temperature.
= Constant, a, in evaporation formulas of Chapter 4 for
KLOSS=1 or 2. For KLOSS=3, EVPCON = 0.01.*
FORMAT 8F10.5
= Number of outflow standard deviations> a » equal to half
the withdrawal thickness (see discussion of Equation 2.10)
= 1.96.
= Inflow standard deviation, a,. Equation 2.2.
= Radiation absorption coefficient, n, Equation 2.1.
= Fraction of solar radiation absorbed at the water
surface, g, Equation 2.1.
= Water density, p, = 997 kg m
= Water specific heat, c, = 0.998 kcal/kg.
= Half the value of the constant of Equation 2.5 used
to predict the withdrawal thickness, 5. For KOH=2, DELCON=0.00461.*
-3
-76-
-------
RMIX = Mixing ratio, r .(Equation 2.15)= 1.0
m
Card 6, FORMAT 1615
NTI = Number of inflow temperatures to be read in.
NTA = Number of air temperatures to be read in.
NSIGH = Number of relative humidities to be read in.
NFIN = Number of insolation values to be read in.
NSURF = Number of surface elevations to be read in.
NDD = Number of values of the diffusion coefficient
to be read in = 2.
NQI = Number of inflow rates to be read in.
NQO = Number of outflow rates to be read in.
NOUT = Number of outlets.
Card 7, FORMAT 8F10.5
DTTI = Time interval between input values of TI.
DTTA = Time interval between input values of TA.
DTSIGH = Time interval between input values of SIGH.
DTFIN = Time interval between input values of FIN.
DSURF = Time interval between input values of SURF.
DTDD = Time interval between input values of DD(= TSTOP for
constant diffusion coefficient.)
DTQI = Time interval between input values of QI.
DTQO = Time interval between input values of QO.
Card Group 8, FORMAT 8F10.5
TI = Values of inflow temperatures, T .
Card Group 9, FORMAT 8F10.5
TA = Values of air temperature, T .
cl
Card Group 10, FORMAT 8F10.5
SIGH = Values of relative humidities, if), in decimal form.
Card Group 11, FORMAT 8F10.5
FIN = Values of insolation, $ .
Card Group 12, FORMAT 8F10.5
SURF = Values of surface elevations, y .
o
Card Group 13, FORMAT 8F10.5
2
DD = Values of diffusion coefficients, D = 0.0125 m /day.
-77-
-------
Card Group 14, FORMAT 8F10.5
QI = Values of inflow rates, Q..
Card Group 15, FORMAT 8F10.5
QO = Values of outflow rates, Q .
Card Group 16, FORMAT 1615
LOUT = Numbers of grid points corresponding ta outlet elevations.
Card Group 17, FORMAT 8E10.5
ELOUT = Outlet elevations.
Card Group 18, FORMAT 3F12.2
SLOPE = Slope of reservoir bottom =0.0 if KTRAV = 1.
10 2
GRAV = Acceleration due to gravity = 7.315 x 10 m/iay .
2
VISCOS = Kinematic viscosity = 0.0864 m /day*
Card Group 19, FORMAT 8F10.5
THICK1 = Observed thickness of inflow layer when traveling
along the reservoir bottom =0.0 if ETRAV = 1.
THICK2 = Observed thickness of inflow layer when traveling
horizontally = 0.0 if KTRAV = 1.
The following parameters are read in when KAREA = 2. (i.e. for field reservoir)
Card 20, FORMAT 1615
NAA = Number of areas to be read in.
NXXL = Number of lengths to be read in.
NWIND = Number of wind values to be read in.
NATRAD = Number of atmospheric radiation values to be read in-..
JMP = Number of grid points for which program variables should
be initialized. (This should be the maximum value of JM
expected to occur in the calculations.)
Card 21, FORMAT 8F10.5
DAA = Vertical distance interval between input values of AA.
DXXL = Vertical distance interval between input values of XXL.
DTWIND = Time interval between input values of WIND.
DATRAD = Time interval between input values of ATRAD.
AAB = Elevation of first (lowest) value of AA.
XXLB = Elevation of first (lowest) value of XXL.
Card Group 22, FORMAT 8F10.5
AA = Values of horizontal cross-sectional areas, A.
-78-
-------
Card Group 23, FORMAT 8F10.5
XXL = Values of reservoir lengths, L.
Card Group 24, FORMAT 8F10.5
WIND = Values of wind speeds, W.
Read in Card Group 25 when KATRAD = 1
Card Group 25, FORMAT 8F 10.5
ATRAD = Values of atmospheric radiation, .
Si
Read in Card Groups 26, 27 when KATRAD = 2.
Card Group 26, FORMAT F10.5, 15
DCLOUD = Time interval between input values of CLOUD.
NCLOUD = Number of cloud percentages to be read in.
Card Group 27, FORMAT 8FI0.5
CLOUD = Percentage of sky covered by clouds.
-79-
-------
oo
$JOB RYAN,KP=29.TIME=5,PAGES=99,RUN=CHECK
C RESERVOIR STRATIFICATION PROGRAM* 1971.
^ MAIM PROGRAM
rOMMDN T(K2,2) ,EL(K2),XL( 1?2),A(K'2) , TI ( 366 ) , TA( 366) , SIGH (366)
CCMMON FIN (36 A), WIND (366 ) , DO ( 366 ) , Q T ( 366 ) , 00 ( 366, 5 )
UOMAX(5) ,UIMAX( 1 ) , DTTI , DTTA , DTSIGH ,3TF IN, DTWI NO, DTDD, DTQI
DTQO, JM , JOUT , J IN, KSUR , KOH , KO , KLOSS, YSUR, YOUT , OT, OY
TSTOP. E VPCON. SPREAD, SIGMA I, S I GMAO
EVAP,RAD,TAIR,PSI,DERIV,HAFDEL,EPSIL
YBOT, BETA , DELCON, V ( 1 "'2 , 1 ) , UI ( If- 2,1),DTT
RHO,HCAP,KMIX,RMIX,JMIXB,MIXED,QMIX,KAREA,DATRAD,ATRAO(366)
AR,WINDY,B<102) ,S( 1'2).EX(102) , EXO( 1D2 ), U0( 102,1 )
OIN(4iSi)) ,TIN(46r; ),OOMIX(i^2 ) , M IXH, DM I X, EX I ( ID 2 ) , OX( 10 21
NOUT,LOUT(5),ELOUT(5),TOUTC( 5) ,UOT(102,5)
SURF (366) , GRAV . SLOPE , VISCOS, L AGT IM( 366 ), EO , E 1 , E2, E3, ET
THICK1 , THICK2, KTRAV, OYSUR , AYSUR
DCLOIin,CLOUD(366),KATRAD. I I I
DIMENSION WH(20) ,AA(1 ^2) , XXL (.1C 2)
C READ IN ALL DATA FOR PROGRAM,
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
-si
•
oo
o
T)
C
rt
IP
i-i
8
00
H
READ (5,9r' )
WRITE (6, 90' )
READ (5,9r!:)
READ (5,9,1)
1MIXED.KTRAV
READ (5.9S.2)
READ (5,9f?)
READ ( 5,9M )
READ (5.9T2)
READ (5,9U?)
READ (5.9^2)
READ (5,9-r??)
READ (5,91<2)
READ (5.9C2)
READ (5,902)
READ (5.902)
DO ^1'.
1=1,2' )
,I=1,2C)
(WH(I ) ,
( WH( I )
(WH(I),I=1,2C)
JM,KATRAD,KSUR,KOH
KO,KLOSS,NPRINT,KAREA,KM IX,
YSUR,DY,DT,TSTDP,TZERO,EVPCON
SPREAD,SIGMAI,ETA,BETA,RHO,HCAP,DELCON,RMIX
N TI,NTA,MS IGH,NFIN,NSURF,NDO,NO I,NOO,NOUT
DTTI,DTTA,DTSIGH,DTFIN,DSURF,OTDD,DTOI,OTOO
(TI(I ) ,I = 1,NTI)
(TA( I ), 1 = 1,NTA)
(SIGH( I ),I = 1,NSIGH)
(FINd ),I = 1,NFIN)
(SURF(I),1=1,NSURF)
-------
B'K'l READ(5,9r>2)( 00 IN, I ),N=1,NQO)
RE AD ( 5,9"! I (LOUT (I ) , I=1,NOUT)
READf 5,9f 2 H FLOUT (I »,I=1,NOUT)
REAO<5,9'<3) SLOPE, GRAV,VISCOS
READ C5,9f2) THICK1,THICK2
YBQT=ELOUT<1 )-OY*FLOAT( LOUTC 1 ) -1 )
GO TO (4,2> , KAREA
C PFAO TM DATA FOR OTHFR THAN LABORATORY RESERVOIR IF INDICATED.
2 REAr?(5,9Cl> NAA,NXXL,NWIND,NJATRAD,JMP
READ (5,9r2) D AA, DXXL,OTWIND ,DATRAD, AAB, XXLB
READ (5,9^2) ( AA( I) , 1=1 , NAA )
RE AD (5, 9C2)(XXL( I ).! = !, MX XL)
READ (5,9f2) ( W IND( I ) , I -1 ,NW I ND)
GO TO (9.1CM.KATRAD
9 READ (5,9^2) ( ATR AD( I ) , 1=1 » MATRAD )
, GO TO 11
~ 1"? READ(5,9f;l) OCLQUD, MCLOUD
1 READ(5,9'r 2)(CLOUD(I ) ,I = 1,NCLOUD)
11 DO 3 1=1 , JMP
RA = (EL» I)-AAB) /DAA
L = RA
Ad) = AA(L + 1 )+( RA-FLOAT(L) ) * ( AA{ L + 2 ) -AA(
RA = (EL( I )-XXLB)/DXXL
L = RA
XL( T ) = XXL(L+l )+(RA-FLOAT(LI ) *( XXL ( L+2 )-XXL ( L
3 B( I )=A( I ) /XL( I )
GO TO 5
4 JMP = JM^-IFIX ( (33«?-YSUR.)/DY+r».5l
C CALCULATION OF LONG WAVE RADIATION IN LABORATORY
AR = •1,7B88E-1?> «< TA( 1 ) +273. 16 ) **4
DO 8 1=1, JMP
B(T) = 3C948
EL(I) = YBQT+DY*FLQAT(I-1)
IF (EL( I )-22«4) 6,7,7
6 XL^< EL ( I ) +87, 0 )
-------
GO TO 8
XL (I) = J.093,
A(I ) = XL ( I)
CONTINUE
WRITE
WRITE
WRITE(6,9r5)(LOUT(II,EL OUT(I),I=1,NQUTJ
WRITE
WRITE
WRITE
WRITE
WRITE
WRITE
INITIALIZE
DO P5;l
6. 9v
(WH(I ), 1 = 1,2'
JM,YSUR,RHO
6,9^6) DY,YBGT,ETA
6,9f7) OT,TZERO,BETA
6,9^8) KTRAV,SIGMAI,HCAP
6.9^91 TSTOP, SPREAD
6,91*M KATRAD.KSURtKDH , KO , KLQSS, KARE A, E VPCON f DELCON, KMIX
(6,923) MIXED, RMIX
MANY VARIABLES.
N=l,36'!
oo
CO
I
85:) CONTINUE
DO 851 1=1,JMP
T(T,1) = TZERO
DTT=DT
JM-MIXEO
FT = 0,^
RAD = nB-r
EVAP = "><
TAIR = -\
EPSIL =
HAFDEL =
JIN = JM
JXM=JM
N = ••>
JMIXB =
-------
DYSUR=YSUR-EL(JM
IF(YSUR-EL(JM)) 858,858,859
85fi AYSUR = A(JM)-(DY/2»D-DYSUR)*-( A(JM)-A( JM-1) ) /DY
GO TO 86"-
859 AYSUR = A(JM»+(DYSUR-QY/?<»f)M A
SARFA1=SARFA
DYSUR1=DYSUR
JM1=JM
Z INITIALIZE VARIABLES USED IN ENERGY CHECK
E'=A( 11/2, r+SAREA*DYSUR/DY
JMM = JM-1.
DO 13 1=2, JMM
13 Er: = E;1-»-A( I)
IF(JM-6^) 15,15,16
15 JP = JM
GO TO 17
16 JP=60
17 GO TO (?0,18) , KO
18 UIMAX(l) =y9l
UOMAX(l) = f.,"j
TS=0.1
DO 19 J=1,JM
SIGMAG = 1.2
C STATEMENT 2* IS BEGINNING OF MAIN ITERATION LODP OF PROGRAM.
20 N=N+1
ET=ET+DT
GO TO (21,5D3),KO
21 GO TO (24,22 ), KMIX
C MIX INFLOW WATER IF INDICATED,,
-------
27 TP=",T>
DO 23 J=JMIXR,JM
23 TP = TP+TU. 1)
TP = TP/FLOAT( MIXED+1 )
TS= ( TTT N ( N ) + TP*R M I X ) / ( 1 . r +R M I X )
GO TO 25
24 TS=TTIMN)
25 CONTINUE
C LDCATE ACTUAL LEVEL OF DAYS INPUT
DO 27 1 = 1 , JM
J = JM+l-J
IF
-------
CO
XLAG=SLDIST/VELF+XL( JINJ/HVELF
87? LAGTIM
QIN(ML»=QIN(ML )+QOIN
-------
36 CONTINUF
DO 38 M=1,JM
IF(ABS(Q!0)-SUM) 39,39.38
38 SUM=SUM+A< JM1-MJ*DY
37 YSUR=EL(JMl)-KM-Co5)^DY4-(OIO-SUMI/A<
GO TO 4"
39 YSUP=EL( JM1 )-(M-Do5 )#DY-M OIQ+SUM
40 DYS=YSUR-EL< JM1J+DY/2.0
IF(DYS) A1,42,A2
4? M=IFIX(OYS/DY)
GO TO 43
41 M=IFIX(DYS/OY)-1
43 JM=JM1+M
DYSUR = YSUR-EL( JM ) -i-DY/2,(;
C CALCULATE MEASURED SURFACE LEVEL
R = ET/OSURF
L=R
RR=R-FinAT(L )
SURMES=SURF(L)+RR^( SURF (L+l) -SURF (L) )
C CALCULATE SURFACE AREA
IF( YSUR-EL (JM) ) 58,58,59
58 AYSUR = A( JMH(OY/2«n-OYSUR)*(A( JM»-A( JM-1) )/DY
GO TO 61
59 AYSUR=A( JM) + (DYSUR-DY/29f))*( A( JM+1 )-A( JM) )/OY
61 SAREA=( ^YSUR + ( A(JM)+A( JM-1) )/2»';)/2«0
C CHECK ACCtfRACY OF SURFACE LEVEL
JMM=JM-1
IF(JM-JMI) 511,511,512
512 00 513 J=JM1,JMM,
513 SUMV=SUMV+A( J):*OY
511 S!JMV = SUMV+SAREA*OYSUR-SAREAI*DYSUR1
GO TO 515
510 JMM1=JM1-1
00 514 J=JM,JMM1
-------
514 SUMV=SUMV+A( JH;DY
SUMV=-( SUMV+SAREA1*'DYSUR1-SAREA*DYSUR)
515 ERRGR=SUMV-OI3
OYCQR=ERROR/SAREA
OYSUR=DY$UR-OYCOR
IF(DYSUR-0Q5) 5t6,5D6,5U7
506 DYSUR=DY$UR+DY
JM=JM-1
. 5f;7 MM=JM-JJM
IF (MM) 44,44,5:}
5" DO 51 1=1,MM
J=JM+I-I
51 T(J,l)=T(JJM,l)
44 T(JM,1>=T(JJM,1)
JMIXR=JM-M!XED
, IF
501 CONTINUE
VM=OY/DTT
IF(VVV-VM) 5^3,504,504
5«H DT=DY/VVV
IDT=OTT/DT+1
DT=DTT/IDT
GO TO 5v5
503 IDT=1
5-5 DO 79 M=l, IDT
-------
C HEAT TRANSPORT CALCULATIONS
OH 1114 J=?,JMM
C DIRECT ABSORPTION TERM
ARJ1=(A< JI+AU+1 ) )/2.T
ARJ2=(A< J)+A( J-l) ) /?.!}
DELTA=M.O-RETAI*FLXIN I/A( J)/DY
C HORIZONTAL ADVFCTION TERM
1162 DELTC = (UI( J, 1 ) - TS-IJCM J , 1 )*:T< J, 1 ) ) *B( J ) *DY/A ( J ) /DY
C DIFFUSI3N TERM
OFLTD=OD(1 )*-•( (T(J + 1,1)-T( J.I) ) /DY? ARJl-( T ( J t 1 » -T ( J-l , 1 ) )/DY*'ARJ2)
1/A( J) /DY
DELT=(OELTA-i-OELTB+OELTC+DELTD)*DT
1114 T( J,2)=T( J.I )*DELT
C CALCULATIONS FOR SURFACE LAYER
IF ( V(JM,1 ) JH63, 1164, 1164
1164 DELTJM=nT*« 1, '/—BETA) ^FLXIN( N ) * { A YSUR-EXP (-ETA^OYSUR )*
1 ( A< JM)+A( JM-1 ) )/29n)/SAREA/DYSUR/HCAP/RHO+
1V( JM,U«CT( JM-l,l)-T(JM,l) J*(A( JM)+A(JM'l) ) /2*i5/S AREA/DYSUR
1+UIt JM.1)*(TS-T
-------
1-OD(1)MT( JM,1)-T( JM-1,1) )/DY«MA(JM)+AUM-l) )/2.G/SAREA/DYSUR
1+(RETA^FLXIN(N)-FLXOUT(N))*AYSUR/RHG/HCAP/DYSUR/SAREA)
GO TO 1165
1163 DFLTJM=nT*({l,,0-BETA)*FLXIN(N)-T(1,1) )*(A(2l+Aa) )/20Q/OY) / A( 1 )/DY*2«D
RO TO 1168
1166 DELTl=DTi--( (1. ^-RETA) *FLX INC N )*EXP I-ETA*( YSUR-EL ( 1 )-DY/2,«) )*
1 (A(2J+A(U)/Zal/RHO/HCAP
2+UT(l,l)*R(1)«OY/2«3*(TS-T{1,1))-V(2,1)*
-------
C SUB AVER MIXES SURFACE LAYERS IN THE EVEMT OF A SURFACE INSTABILITY.
6^ IF
-------
GO TO (85,89), KO
85 F = 2»^-HAFDEL
WRITE (6,918) EPSIL,F,SIGMAO
DO 88 I=1,NOUT
89 WRITE (6,917) ELOUT(I),OOUT(N,I),TOUTC
DO 9« 1=1,in
93 WRITE (6,921) 93 LL = 7Q
7* 94 DO 95 1 = 61,LL
95 WRITE (6,921) ( J , EL( J ) , T (J, 1 ) , J=I , JM, 10 )
lf?'1 IF (ET-TSTOP) 2a,l»l
1 CONTINUE
9GO FORMAT (2C«A4)
9f?l FORMAT (1615)
9P2 FORMAT (8 PI?.5)
903 FORMAT(3F12. 2)
9!H FORMAT (« NUMBER OF GRID POINTS=•I 3,17X,'SURFACE ELEVATION='F7«2,
118X,«OENSITY='E12.5)
955 FORMAT (' OUTLET LEVEL=»I3, 26X,'OUTLET ELEVATION='F802)
°36 FORMAT (• D/='F6.2.33X,'BOTTOM ELE VATIQNJ='F8o 2,1 8X, • ETA= • F6. 3 )
907 FORMAT.(' DT=«F6.2,33X,'INITIAL TEMPERATURE='F6»2,17X,'BETA=•F5«2)
9GB FORMAT (» KTRAV='I5 ,31X,'INFLOW STD» DEV.='F6.2,21X,
I'HEAT CAPACITY='F8a5)
9D9 FORMAT (' STOP AT TIME=«F7a2,22X,'OUTFLOW SPREAD CnNSTo=«F5,2)
910 FORMAT (» KATRAD='I 2,15X,'KSUR=•I 2,13X,'KOH ='I 2,15Xf«KQ=«I2,16X,
l'KLOSS=«I2,13X,'KAREA=«I2/« EVAPORATION CONSTANT=»E11.4,1DX,
2 'CONST IN EON FOR OUTFLOW DELTA=»F5.2,7X,•KMIX=«12)
-------
911 FORMAT ( ""N = 'I3,'v ABOUT TO ENTER SUBROUTINE AVER,*)
912 FORMAT ( ' ELAPSED TI VE=•F7,2,22X,«ACTUAL SURFACE ELEVATION*1
IF7,2.1IX,•INFLOW TEMPERATURE=•F6«2 )
913 F3RMAT <« NO, OF TIME STEPS= ' 14 ,21 X , • SURFAC E ELEVATION USED= •
1F9.2.11X.'AIR TEMPERATURE='F6.2>
914 FDRMAT (• N0a OF GRID POINTS= • 13 ,23X , • EL EVAT ION OF INFLOW='F7,2,
116X,'RELATIVE HUM ID IT Y= « F5» 2 )
935 FORMAT (' LEVEL OF INFLOW=•T3,23X,'EVAPORATI ON FLUX=«E12«5,14X,
1'INSOLATTON FLUX=•El?.51
916 FORMAT (• HEAT LOSS FLUX=•E12»5,15X,«RADIATI ON FLUX=»E12.5,
116X»IINFLOW RATE='Fllol)
917 FORMAT (• OUTLET ELEVATION=•FID,5,15X,• OUTFLOW RATE=•Flr,1,19X,
I • OUTFLOW TEMPERATURE='Fft02lC')
918 FORMAT (« FPSILQN=«El1.4,23X,•WITHDRAWAL THICKNESS=«F7.2t15X,
1'nilTFLOW STD, OEV0='F6.2)
925 FORMAT (/7(• J ELEV TEMP(C) '))
vo 9?1 FORMAT{7{ I 3, F6« l«F6a 2 t 3X1 )
V 9?? FORMAT (/' TIME FROM BEGINNING OF CALCULATIONS FOR THIS DATA SET=«
113,' MINUTES,'i3,« SECONDS,,'>
923 FORMAT {• NO. GRID SPACES IN MIXED LAYER='I 3,8X,«MI X ING RATIO*'
1F5.2*
924 FORMAT (' TEMP OF MIXING LAVER='F6«2,15X,'MIXED INFLOW TEMP=«
1F6.2.19X,'TOTAL ENERGY RATIO='F9.5)
925 FORMAT (« MIXING OEPTH=•F5.2,24X,•ATMOSPHERIC RADIATION='E12a5,
1 9X,«WIND SPEED='F5»2)
926 FHRMATC ENERGY RATIO=«F7.3, 5X.'ENERGY ST3RED='E12.5, 5X,'ENERGY
1INFLOW='E1295, 3X,'ENERGY OUTFLOW=«E12»5)
CALL EXIT
END
-------
FUNCTION FLXOUT(N)
C CALCULATION OF SURFACE LOSSES DUE TO EVAPORATION, CONDUCTION, AND RADIATION,
COMMON T(102,2),EL(1(»2),XL<102),A<1D2),TI(366),TA(366),SIGH(366)
COMMON FIN<366),WIND(366),OD(366),01(366),00(366,5 I
COMMON UOMAX(5),UIMAxm,DTTI,DTTA,DTSIGH,OTFIN,DTWINO,DTDD,DTQI
COMMON DTOO,JM ,JOUT,JIN,KSUR,KOH,KQ,KLOSS,YSUR,YOUT,DT,DY
COMMON TSTOP,EVPCON,SPREAD,SIGMAI,SIGMAO
COMMON EVAP,RAO,TAIR,PSI,OERIV,HAFDEL,EPSIL
COMMON YBOT,BETA,OELCQN.VC102,1 I,UI(1G2,1),DTT
COMMON RHO,HCAP,KMIX,RMIX,JMIXB,MIXED,QMIX,KAREA,DATRAD,ATRAO<366)
COMMON AR,WINDY,B(1^2),SUQ2),EX<102),EXO(1C2),UQ< ID2,1)
COMMON QIN(46a),TIN(46f5) ,OOMIX<1!)2) ,MI XH,DMI X,EX I ( 102) ,OX (102 )
COMMON NOUT,LQUT(5),ELOUT(5),TOUTC<5) ,UOT(102,5)
COMMON SUPF(366),GRAV,SLOPE,VISCOS,LAGTIM<366),EO,El,E2,E3,ET
COMMON THIC<1,THICK2,KTRAV,OYSUR,AYSUR
COMMON DClCUD,CLOUD(366),KATRAO
C KLOSS = 1 FOR LABORATORY USING ROHWER FORMULA,,
C 2 FOR FIELD USING KOHLER FORMULA,,
C 3 FOR FIELD USING ROHWER FORMULA*,
ET=DTT*FLOAT(N)
R = ET/DTTA
L = R
RR = R-FLOAT(L)
TAIR=TA(L)+RR«(TA(L-H)-TA(L) )
R = ET/DTSIGH
L = R
RR = R-FLOAT(L)
PSI=SIGH
-------
1
2
*( TAIR+273. 16) **«•
M 273., 16+TS )*#4-AR
FIR
20
15 CHI = RHO
IF(FVAP)1,1,2
EVAP=0.t
EVAP=EVAP+CONOUC
UNITS OF RADIATION ARE CAL/CM-CM-MIN.
AR = 0.7888E-K'
RAO = 3.7888E-13
w = r?.a
FLXOUT=EVAP+RAD
RETURN
FIELD DATA, WIND SPEED IS IN M/SEC.
R = ET/DTWIND
I = R
W=WINO(L)+(R-FLOAT(L) > * ( WIND ( L-H ) -W I NDl L) )
WINDY=W
C CALTULATinN OF ATMOSPHERIC RADIATION
C UNITS OF RADIATION ARE KC AL/M-M-DAY.
GO TO (7,8),KATRAD
7 R=FT/DATRAD
L = R
AR=ATRAD(U + (R-FLOAT< L) )*( ATRAD{
RAD = 1,13587E-6*(TS+273,16)^*4-AR
GO TO 9
8 R=ET/DCLOUD
L=R
CC = CLOUD (L )-KR-FLOAT(L» )*( CLOUD ( L + l ) -CLOUD ( LI )
RAD=l»13587F-6*-( ( TS+273, 16 J ^*4--'J. 937 E- 5*( TAIR+273,, 16)**6*( loO+
)-ATRAD( L ) )
9 GO TO (15.25.301* KLOSS
C CALCULATION OF FIELD FVAPORATION USING KOHLER FORMULA
C VAPOR PRESSURES IN MB,
25 DE = DE/3- 75 "i "5 62
FVAP = H*DE+HCAP*DE*TS
EVAP = EVPCON^RHO^W**EVAP
-------
0*(TS-TAIR)
IF(EVAP»3t3,4
3 EVAP=0«C
4 EVAP=EVAP+CQNDUC
FLXQUT=EVAP+RAD
RETURN
C CALCUIATION OF FIELD EVAPORATION USING ROHWER FORMULA*
30 CHI = RHO#(H*DE+TS*HCAP*DE)
EVAP = CHI*FW !!EVPCON
CONDUC=RHO^EVPCON*269al«(TS-TAIRMFW
IF(EVAP)5,5,6
6 EVAP=EVAP-«-CnNDUC
FLXOUT =EVAP-»-RAD
RETURN
1 FND
tO
(7k
I
-------
SUPROUTINE TOUT(HEATOT,FLOWOT )
C COMPUTE WEIGHTED AVERAGE OF OUTFLOW TEMPERATURE,
COMMON T(K2,?) ,EL(K?) ,XL( 102) ,A(1?2) ,TI (366) , TA ( 366) t SIGH (366)
COMMON FIN(366 ) .W IND( 366 ) .00(366 ). 01 (366) ,00(366,5)
COMMON UOMAX( 5) ,UIMAX(1 ) ,OTTI , DTTA , OTSI GH, DTF IN, DTW I NO , DTDD, DTQI
COMMON DTOO, JM.JOUT. J IN. KSUR , KOH , KQ, KLGSS, YSUR, YOUT,DT,DY
COMMON TSTOP, EVPCON, SPREAD, SIGMA I, S I GMAO
COMMON EV AP, R AO, TAIR.P SI. OER IV. HAFOEL.EPS I L
COMMON YBQT, BETA. DELCON.VU1* 2,1 ).UI( 102.D.DTT
COMMON RHO.HCAP,KMIX.RMIX.JMIXB,MIXEO,QMIX,KAREA,DATRAD,ATRAD(366)
COMMON AR,WINDY,B(lf,2),S( 1^2 ) , FX( 1C 2 ) , EXO( 1 <. 2 ) ,UO( 102,1 )
COMMON OIN(^6»,TIN(46r ), COM IX (132 ) , M IXH, DM I X, EX I ( ID 2) , OX( 132 )
COMMON NOUT,LOUT(5) ,ELOUT(5) ,TOUTC(5) ,UOT(182,5)
COMMON SURF(36fc),GRAV,SLOPE,VISCOS,LAGTIM(366),E;?,El,E2,E3,ET
COMMON THICK1 ,THICK2 , KTRAV, DYSUR , AYSUR
COMMON OCL QUO, CLOU D( 366) .KATRAD, III
I
to
HEATQT=T(JM,1)*B( JM)*UOT( JM,I
FLOWOT = B(JM)*UOT(JM, I
JMM=JM-1
00 2 J=2.JMM
HEATOT = HEATOT-»-UOT( J, I)*B( J
FLOWOT=FLOWOT+UOT( J,I )*B( J)*DY
IF(FLOWOT.EO,0«!i) FLOWOT=1.0
CONTINUE
RETURN
FND
*UOT
)*DY/2.0
)*UOT(1,I
-------
SUBROUTINE SPEED(N)
C COMPUTATION OF VERTICAL AND SOURCE AMD SINK VELOCITIES.
C ALSO, COMPUTATION OF WITHDRAWAL THICKNESS,,
C SOURCE AND SINK VELOCITIES ARE ASSUMED TO HAVE GAUSSIAN DISTRIBUTION,
COMMON T<102,2),EL(ir:2),XL{lD2),A(102) t T I ( 366 ) , TA( 366) , SIGH (366)
COMMON FIN (366) ,W IND( 366 ) , DD ( 366 ) ,01(366) ,00(366,5)
COMMON UQMAX(5),UIMAX(1) ,DTTI , DTTA , OTSIGH, L)TF IN, DTWI ND,DTDD, DTOI
COMMON DTOO, JM,JOUT , J I N,KSUR, KQH,KQ, KLOSS » YSUR t YQUT , DT, DY
COMMON TSTOP, EVPCQN, SPREAD, SIGMAI .SIGMAO
COMMON EVAP,R AD, TAIR, PS I , OER IV , HAFDEL , EPS IL
COMMON YBOT,BETA,DELCON,V(1)2,1 ) , U I ( 1 02, 1 ) , DTT
COMMON RHO,HCAP,KMIX,PMIX.JMIXB,MIXED,QMIX,KAREAfDATRAD,ATRAD(366)
COMMON ARfWINDYtRU*21*Sf l?2)«EX*OY/2oO+EXI(JM)*B(JM)*DYSUR
JMM=JM-1
DO 3 1=2, J MM
3 VOLIN=VOLIN+EXI( J)*'B( J)*DY
-------
UIMAX( 1)=QIN( N)/VOLIN
GO TO <8,7),KMIX
7 UIMAXC1 )=UIMAX(1)*(1,0+RMIX)
8 DO 6 J=1,JM
6 UI (J,lUUIMAX(l)*EXnJ)
C COMPUTE OUTFLOW VELOCITIES
DO 10 LT=1.NOUT
JOUT=LOUT(LT)
C COMPUTE WITHDRAWAL THICKNESS.
C NOTE THAT ONLY HALF THE WITHDRAWAL THICKNESS IS COMPUTED
DERIV = (T{JOUT+1,1)-T< JOUT-1,1) ) /2oO/DY
IFCDERIV-ft.nin > 11,11,15
11 JOUTl=JOUT+2
C CUTOFF DUE TD SHARP CHANGE IN DENSITY GRADIENT
00 I?. J = JOUT1,JMM
IF( (T( J+1,1)-T(J,1) )/DY-0rf5)12,13f 13
12 CONTINUE
GO TO 19
13 HA FOE 1= FLOAT < J-JOUT)*DY
S I GMAO=HAFDEL/ SPREAD
19 JOUT2=JOUT-2
DO 21 I=1,JOUT2
J=JOUT2+2-I
IFC(T(J,1)-T(J-1,1) J/OY-O.OS) 21*21,22
21 CONTINUE
GO TO 14
22 HAFD1=FLOAT( JOUT-J )*DY
SIGM1=HAFD1/SPREAD
IF( SIGM1. oLT, STGMAO) SIGMAO=SIGM1
GO TO 14-
APPROXIMATING FORMULA USED FOP DENSITY IS RHO=1* 0-0« OCir:On663*( T-4, 0 J**2
15 FPSIL= ?»n*(T( JOUT,l)-4«C)/U5inG000-(T( JOUT , 1 ) -4e 0 ) **2 )*DERI V
GO TO (17,16) ,KOH
CALCULATION OF WITHDRAWAL THICKNESS USING KAO FORMULA.
16 OPUW = OOUT(N,1 T )/R( JOUT)
-------
HAFDEL = DELCON*SQRT(OPUVM/EPSIL**0«25
GO TO 18
CALCULATION OF WITHDRAWAL THICKNESS USING KOH FORMULA,
17 HAFDEL = DELCON/EPS IL**D. 1666667
18 SIGMAO = HAFDEL/SPREAD
IF(SIGMAQ) 20,29,14
2^ SIGMAO=l0r-
14 CONTINUE
COMPUTE EXPa FACTOR
DO 10H 1=1, JM
S( I )=
-------
31 UO(J.U=0«'
DO 35 LT=1,NOUT
35 UnU,l»=UO(J,l UUOT( J.LT)
36 CONTINUE
C COMPUTE VERTICAl ADVECTTVE VELOCITY
Vf2tl)=(lJIflf 1 )-UO(l,11 )* R(1)*DY/(AU H-A<2))
I im V — I M .L 1
vl P^ A "" vl ™ * 1
OH 5^0 J=3,JMX
o
i
V(J,1>=(V(,J-1,1) '(A(J-2 )+A( J-l)
l«nY)/(A( J)*A( J-l) l*=2.r.v
5?~ CONTINUE
RETURN
FNO
-------
SUBROUTINE AVEP(N)
C PERFORMS COMVECTIVF MIXING OF SURFACE LAYERS,,
:OMMON
COMMON!
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
COMMON
o
IV)
i
2) ,ELiio2) ,XL(1Q2) ,AU32) ,Ti<366),TA(3&&),siGH(3&6)
FINI (366 ) , W IND( 366 ), DO ( 366 ) ,01(366) ,00(366,5)
UOMAX(5),UIMAX(1),DTTI,DTT A,DTSIGH,OTFIN,DTWIND,DTDO,DT01
DTQD,JM,JOUT,JIN,KSUR,KOH,KQ,KLOSS,YSUR,YOUT,DT,DY
TSTOP,EVPCON,SPREAD,SIGMA I,SIGMAQ
EVAP,RAD,TAIR,PSI,OERIVtHAFDEL,FPSIL
YBOT,BETA,DFLCON,V(152,1),UI(102,1),DTT
RHO,HCAP,KMIX,RMIX,JMIXB,MIXEDtQMIX,KAREA,DATRAO,ATRAD(366)
AR , WI NDY , B (1D2 ) , S (10 2 ) , EX (132 ) , EXO (10 2 ) , UO ( 132 , 1 )
QIN(46'1),TIN(460) ,QQMIX(102) , MI XH,DMI X,EX I ( 102 ) ,OX( 102 )
NOUT,LOUT(5),ELOUT(5),TOUTCt5),UOT(132,5)
SURF(366),GRAV,SLOPE,VISCOS,LAGTIM(366),ED,E1,52,E3,ET
THICK1,THICK2.KTRAV,DYSUR,AYSUR
DCLOUD,CLOUD(366) .KATRAD
JMM=JM-1
DO 5 1=1,JMM
JJ=J-1
IFCKJ.l )-T< JJ,1) ) 6,7,7
6 CONTINUE
TF(J-2) 8,8,9
3 T( 2 . 1) = (T(211)*A<2)+T(1,1 ) *
GO TO 7
9 DO I-1) K = l , JJ
KJJ=KJ-1
TF(JM-KJ) 2,2,3
2 AREA=( AYSUR-KA(JM)+A( JM-1) )/2,
GO TO ^
3 AREA=A(KJ)
4 AV1 = AV1+T( KJ,1)*AREA
-------
o
CA>
I
AV2=AV2+AREA
TAV=AV1/AV2
IF(TAV-T(KJJtlM 1
1»> CONTINUE
20 IF(J,EO«JM) MIXH=K
DMIX=( MIXH-1 )
DD 31 L=KJ.J
CONTINUE
RETURN
END
-------
FUNCTION TTIN(N)
C COMPUTE INFLOW TEMPERATURE FROM READ IN VALUESo
C LINEAR INTERPOLATION BETWEEN READ IN VALUES.
COMMON T(l"2*2) tEL(l?2 )*XL(lf,2lt A(1Q2) » TI ( 366 ) ,TA( 366) , SIGH (366)
COMMON FIN(366),WIND(366),DD(366),QI(366),00(366,5)
COMMON UOMAX(5),UIMAX(1),DTTI,DTTA,DTSIGH,DTFIN,DTWINO,DTDD,DTQI
COMMON DTQO,JM,JOUT,JIN,KSUR.KOH,KQ,KLQSS,YSUR,YOUT,DT,DY
COMMON TSTOP,EVPCON,SPREAD,SIGMAI,SIGMAO
COMMON EVAP,RAD,TAIR,PSI,DERIV,HAFDEL,EPSIL
COMMON YBOT,BETA,DELCON,V(132,1),UI(1D2,1),DTT
COMMON RHO,HCAP,KMTX,RMIX,JMIXB,MIXED,QMIX,KAREA,DATRAD,ATRAD(366)
COMMON AR,WINDY,B(19?),S( K> 2),EX(It 2),EXO(102),UO(102,1)
COMMON QIN(460),TIN(460),QQMIX(1^2),MIXH,DM I X,EX I(132),OX(102)
COMMON NOUT,LOUT(5),ELOUT(5),TOUTC(5),UOT(132,5)
COMMON SURF(366),GRAV,SLOPE,VISCOS,LAGTIM(366),ED,El,E2,E3,ET
, COMMON THICK1.THICK2.KTRAV
Q ET=DTT^FLOAT(N)
-f R = ET/DTTI
L = R
RR = R-FLOAT(L)
TTIN=TI(L)+RR*(TI
-------
FUNCTION FLXIN(N)
C COMPUTE INCCMINC, SOLAR RADIATION FROM READ IN VALUES.
C RFAD IN VALUES TREATED AS A STEP FUNCTION.
COMMON T(lr2,2).EL(lD2J,XL(l')2),A(102 ),TI<366),TA(366),SIGH(366)
COMMON FIN(366),WIND{366).OD(366),QI(366),QO(366,5)
COMMON IJQMAX(5)«UIMAX<1 ), DTTI, DTT A, DTSI GH, DTF IN, DTWI ND ,DTDD , DTQI
COMMON DTQO,JM,JOUT,JIN,KSUR,KOH,KO,KLOSS,YSUR, YOUT,DT»DY
COMMON TSTOP,EVPCHN,SPREAD,SIGMA I,SIGMAQ
COMMON FVAP,RAD,TAIR,PSI,DER IV t HAFDEL, EPS IL
COMMON YBOT, BETA, DELCON,V(1'» 2,1) ,UI (1>>2 » 1) , DTT
COMMON RHO,HCAP,KMIX,RMIX,JMIXB,MIXED,QMIX,KAREAtDATRAD,ATRAD(366)
COMMON AR,WINDY,B(1'121,S <1?2),EX(102),EXO(152),UO(102, 1)
COMMON OIN(461),TIN(460),OOMIX( H'2) ,MIXH,DMIX,EXI(1C 2 I,OX(1J2»
COMMON NOUT,LOUT(5),ELOUT(5),TOUTC(5),UOT(102,5)
COMM1N SURF(366) ,GRAV,SLOPE,VISCOS,LACTIM(366),EG.El,E2,E3,ET
_^ COMMON THICK1,THICK2,KTRAV
"" ET=OTT*FLOAT(N)
R = ET/DTFIN
L = R
FLXIN = FIN(L)
RETURN
END
o
en
i
-------
FUNCTION OOIN(N)
C COMPUTE INFLOW RATE FROM READ IN VALUES.
C READ IN VALUES TREATED AS A STEP FUNCTION,
C f'MMON T(102,2),EL(1C2)« XL(13 2)t A(1G2 I•TI< 366)iTA(366),S1GH< 366)
COMMON F I \M366), WIND (3 66 ),DD<366) ,QI (366) ,00(366,5)
COMMON UOMAX<5),UIMAX(1)»DTTI,DTTA,DTSIGH,OTFIN,DTWIND,DTDD,DTQI
COMMON DTQO.JM,JOUT,JIN,KSUR,KQH,KQ,KLQSS,YSUR,YC}UT.DT,OY
C3MMON TSTOP.EVPCON,SPREAD,SIGMAI,SIGMAD
COMMON EVAP,RAD,TAIR,PS I, DERIV, HAFDEL, EPS IL
COMMON YBOT, BETA,DELCON,VUi>2tl) , UK 10211) ,OTT
COMMON RHO.HCAP.KM!X.RMlX*JM!XB,MIXEDtOMIX,KAREA,DATRADtATRADO66)
C OMMON AR,WINDY,B(1J2),S(1"2),EX(It 2),EXO
-------
FUNCTION OOUT(N,I)
t COMPUTE OUTFLOW RATE FROM READ IN VALUES*
C READ IN VALUES TREATED AS A STEP FUNCTION,
COMMON T(1-2,2),EL(102).XL(102),A(1"2),TI(366),TA(366),SIGH(366)
COMMON FINn66),WIND<366),DD(366),OI(366),QQ(366,5)
COMMON UOMAX(5),UIMAX(1),DTTI,DTTA,DTSIGH,DTFIN,DTWIND,DTDD,DTOI
COMM1N DTQO,JM,JOUT,JIN,KSUR,KOH,KO,KLOSS,YSUR,YOUT,DT,DY
COMMON TSTOP,EVPCON,SPREAD,SIGMA I,SIGMAO
COMMON EVAP,RAD,TAIR,PSI,DERIV,HAFDEL,EPSIL
COMMON YQOT, BETA,DELCON,V(1D2,1).UI1102,1),DTT
COMMON RHQ,HCAP,KMIX,RMIX,JMIXB,MIX ED,OMIX,KAREA,DATRAD,ATRAD(366)
COMMON AR,WTNOY,B(1^2),5(102),EX{102),EXO(102),UO(172,1 )
COMMON QIN(46"M,TIN(46rn,QOMIX(n2),MIXH,OMIX,EXI(lQ2) ,OX(i:2)
COMMON NOUT,LOUT(5),ELOUT(5),TOUTC(5),UOT(1^2,5)
COMMON SURF(366),GRAV,SLOPE,VISCOS,LAGTIM(366),ED,E1,E2,E3,ET
COMMON THICK1,THICK2,KTRAV
ET = DTT*-FLOAT(N)
R = ET/DTOO
L = R
00UT = 00(L,I)
RETURN
END
-------
*FNTRY
1FIELD DATA FOR FONTANA RESERVOIR FOR MARCH 1 TO DECEMBER 31, 1966.
-UIL UNITS IN METERS. DAYS. KILOCALORIES, KILOGRAMS, AND DEGREES CENTIGRADE,
o
47 1
492»49
1,96 4.
3>"6 3D 6
1 = 5
7. 58
6n ;"j 6
If; .84
9,, 4r
9n91
9«95
13, 87
15»68
13,26
12«59
I5o8:j
16,36
17a98
19,73
19, 39
19,4-7
2H.43
23,22
2:% 54
18,**4
18,49
19*55
19,30
18,79
18,77
17, 21
15«75
13, 69
2 2
?oT
C
356 3r: 6
I.1)
7,24
6,67
ir»83
8, 52
i%23
1" , 41
13, 55
15,00
13,41
12,28
16,46
17, 13
18, O'j
19325
19,28
19*29
2U06
21.49
21*25
18*45
18,71
2C'*B6
18, 66
18,96
19»OD
16»43
15,14
14, 29
1 3
l.J
.75 C.
3 '56 2
1»"
6a93 -
7.61
11*41
8,25
1C'. 93
11,63
12, 81
14. 55
13*93
13,28
16*97
17. C- 6
18.37
18. 53
19,34
18a76
20.75
21* 64
21.19
1 8a 44
18,48
19,96
18,48
1.9,14
18.76
15,36
15,89
14844
13 2
305.0
5ft 997,
356 306
1.0
8927
8,29
11,39
8,27
9,75
12,29
13.31
14079
14,38
14,30
17,25
16.2'?
18. 65
18,91
19.87
18,81
22, 16
22, 3H
21.53
18, 27
17, 83
19,, 89
18,62
18.84
18, 18
16,03
15a95
13.94
2 4
6*7
0 0.
1
1,0
8.35
9613
HB94
80 50
9047
12.32
13997
14993
15.14
13.81
16026
15.68
18,46
19a42
19»67
20'e31
23a63
22.84
21.33
18,31
17«62
19*7"
19.12
19.17
17,16
16.28
16036
13.97
1
D.01
998 9,
306=0
6.34
9,96
11.36
8,73
9.55
12066
14.25
14, 3D
15.11
13,86
15.73
15,31
18,00
18007
20.57
20«81
24.24
21.37
21-.*?8
18.48
17.74
19.54
19.58
19*08
17,37
15.85
16.32
13.89
00461 1,00
1.4
5.62
1007!5
12011
9,83
9845
12^45
15.29
13,59
13.81
14.16
15,81
15-. 78
18,58
18a42
19.39
2C*55
23a73
20.65
2U« 63
18,48
18914
19,18
18.84
18.26
17, 3H
16.00
1 5, 20
13,98
i.C
5,51
10,97
11,68
9*93
9.48
12,92
15.42
12.86
12»83
14.51
158 51
16,91
2")B45
18.88
19,75
20..59
23,74
2D.48
19.01
18.20
18,91
19.86
18.44
18,88
17,28
16*10
13,84
14.28
vO
OT
03
M
I
rt
-------
o
10
14,12
11.23
12,11
8.. A3
12,22
9, 72
8,55
6,21
6*82
6,45
4*69
6,^34
l«"9rr
11,451
2,285
7S706
3, 872
14,503
15.288
12,76r*
13.811;
18, 252
190289
19* 189
19,782
1 9,^87
2% 961
2", 7( 5
21,89^
?•% 2:7-1
19..991
21,359
23,398
18*333
19,851
18,398
13, 78
1' .82
11,81
8, 81
12,14
8,65
6,68
7,78
6.47
6,63
4,82
5,248
4*^79
IT, 478
1,821
9a162
8B892
14032Q
15,481
14,871
14«753
15.99C,
16«884
i8«24r-
17*334
19oll9
2f*126
2: «616
21,837
2^ *539
2U,521
2' ,835
24, ?36
1 8. 157
19,873
17.52"1
13, 84
1C* 84
11™ 49
K.< 1
lln6<-:
8,61
6.29
8«84
6,24
5»82
8.752
5«834
1%325
2.808
7a!X 9
14, 036
15. 785
17,253
16,124
16»r77
18,244
15, 868
2-eOf;l
18, 962
20.^22
19.251
21e225
21, 199
22, 264
20. 838
2% 423
22. 82H-
18, 172
2G.279
180C19
14«25
11, 11
IB. 93
l! ,92
1C, 7!
8.53
6, 33
9388
6,46
4. 08
11. 746
7.473
10,118
1.612
4» 49^
14,995
16.072
15.828
17,406
14.543
19.78D
11. 156
2C-.13P
18, 722
20, 839
21.147
22.387
23.723
22.816
20.266
2U372
21. 81 C
17. 669
20,783
16, 330
14.28
11.58
11.38
11.14
9.91
8.04
6» 46
9019
6.49
3.76
C0281
9.3^6
ID -.673
2.176
4.0 18
l'o,94-8
156 669
16. 199
16, 865
12c711
15.585
1^.545
2^,709
18.485
21» 466
22.369
24.D86
21.759
23*756
2f-,037
2ia924
22, "92
19.054
20.953
16,854
130'J4
12, 46
12.C3
11.87
9. 64
8.55
6o27
8,24
6.59
4.12
-3,432
12.868
15*065
6*157
60440
9a377
15.946
16,398
1303r'9
19,122
179954
12,878
20,640
16,408
20.298
20,963
24.779
22.718
23.990
2^.121
22*145
22.846
2<}9 546
18.379
18.640
13.02
12.79
10.81
11.14
K«*14
9.47
5.68
7.22
6.46
5.05
-2.56f;
11.472
13a138
6.616
5.253
7.864
16,830
14.063
9, 9 5 6
17.264
17.391
15.439
19,662
18.199
21,017
22,457
23.719
2?. 908
22.134
20,964
22,389
21.757
18.785
18.892
15.742
11.74
12.56
9. 2C
12.52
1J027
90 )8
5.26
7.^2
6.31
5.38
-1,432
lloJ 03
4» 11 7
7,283
6.247
9,398
16.717
12,950
12,2)4
18.814
17,497
170855
20,2F>6
18.628
21,295
21,289
24,490
21.449
2-% 767
25.333
229619
2Q.887
19.387
18.478
17.Q59
-------
15,840
13,559
10, 174
9*513
80749
a.25»-
*,428
13.C76
6,94^
-^2<">7
5,543
•1.497
5,695
-3. 691
L5, 687
f%726
9.512
%655
0,480
0«592
f?a 629
1.899
%620
(?, 849
n,&93
0.811
Q , 7-' - 4
0, 751
0, 740
0.878
0,851
3,911
0. 741
C, 818
?'}„ 876
£•,855
14.965
14.2L-8
13.469
9*866
8,107
7,339
5,3D4
1^,286
3*493
-1.522
9»435
-0.226
5,164
1,718
fc.7^7
C.674
0,730
r.a852
Ct.,618
0*568
0.757
ft«9C 7
^643
«.95t>
C«818
"«856
•::-.775
C.9C5
C«771
0«886
C«827
;Ba47
0.707
5.816
il»873
n,,832
16.749
18,583
13. 655
12.215
7B193
8. 6C.2
If.*. 52 3
10.962
4.233
1.251
13.2.33
1.092
^4* r;45
%899
D.66C
5.661
C.65*
0,835
0.681
0»755
t. 857
')»670
0.789
C. 861
0.671
0. 77 1
D»744
Cc741
0.951
%799
i>»888
0.7D3
D.832
0. 888
Q.875
18.1D1
17.673
8. 635
14.964
10,727
8.783
10,444
6e45!>
4,929
0,795
16.696
1,279
-4.841
0.80ft
0,688
0.544
0.699
0. 546
0.769
0.852
0,924
0.647
3 . 73 9
0, 74C;
3.666
9.806
3.825
D.765
0. 860
0,811
0,827
0. 779
Go85!J!
0. 92P
0*912
16,357
16*044
9,879
17.684
14o354
9.233
10. 746
5,769
5. 464
2.252
1D.667
2,979
-0.781
C.845
P. 767
0.696
0*658
0.583
0,793
0.922
0, 960
0.738
Oo945
OB 882
0«668
t«858
9, 798
D»778
0.819
D. 772
009DO
0, 764
3.825
T/.931
00 882
14.567
15.379
11.651
14.084
14.866
10,010
12,747
6«126
7. 803
0,443
3,067
3,808
Da409
0, 989
0.73G
0.669
0.569
0.597
0. 795
C»796
0«929
D.767
0*781
0.871
0.683
9.806
SB 906
0,825
Do848
0.778
0.761
i?«775
&.8S9
3',886
08837
15.396
10.724
13.383
6.452
13,859
2.928
14.179
9.373
11,906
-0.513
2.058
3.705
4»788
0. 992
0.943
0,797
0. 609
0.785
0. 715
0.699
0.888
C«. 562
0.915
0.948
3,739
0.751
0.774
0*848
0.833
0.828 .
3,801
(1.920
Oe858
5. 875
0.917
13.470
8.225
16,243
10.092
Ilol52
-1.490
12,162
10.375
ID, 372
1,178
1.567
3,997
-1,212
!J. 807
0.841
0.637
0.652
0.567
0,653
0.750
13.684
5.565
De806
0,950
08730
0,697
5,753
G.846
0.874
6. 847
0.743
0,912
0.878
5.853
0.807
-------
X818
CU 855
"•i » 8 ? 9
D« 87 D
'."», 873
"..,? *> B * >
0,751
",846
'"* « 8 5 5
0. 794
% 873
••>, 790
' -• ,917
0.812
""-, 854
0,848
3641, 3C6
423" a 89"
388^,664
5<~ 68, 949
4949^656
5826, 9C 2
5734. 128
24r 30 943
7C.44, 164
7f'48,488
7^1 5* 25~
5559* 433
7924,957
6625, 542
7283, 171
3S86-, 595
5713, 582
4666.128
8765, 7r3
Co 876
C, 8 54
0«837
0,945
U 978
r s 8 6 2
0,800
•'"•« 7?'"f
rr«8f 5
'\ 7 3 1
• ,932
'.,891
1 r, "' Vi '•?
r. •• •> 9 9 1
{,892
r.:.979
3136o 242
4277,375
4A..7-632
31-17.281
522U113
5365*847
3 f 8 4, 714
3279*037
6846,023
3290,461
60 38 o 2 46
6334,777
692nfll95
3973, f)46
7374,636
3^89. 93C-
5293,636
6914,285
7922,957
1,834
G. 859
t. 837
f/.96D
0.915
n.785
^f. 748
f',769
^o 816
'». 816
C.851
".878
-------
ro
i
73 52 . 761
3576<,f'29
7157, 363
7f !-(6» 195
569?0394
7688.058
2635.14L*
6453,46''
4987a582
6796, T'l 1
7783. «13
1898* 651
4913,996
4285, 332
3541-, 468
1331,164
1725. H85
2561. 5?0
1824, 208
1"23, 124
492*493
498,287
5^,582
5'n»771
501.298
500ffl856
5^792
501.119
506,218
568,382
539.522
51D»555
511* 101
511.308
510, 863
510, f 52
8427.151
&r?7. f^7
7423«613
3184.957
6765.558
6r?r-5,988
262^.025
2491.281
6679, IT 5
6163.191
4513»?&1
f 133*597
4644*339
1572.299
1755.187
1316,963
1213,383
2 8?' » 626
1^62.705
K18.83C'
492a928
498fl64^
5r»r.63i
5 11 „ 6 46
5r 1.451
5*r.673
5-G.591
5r; 1.478
5;)6,632
508.687
5r;9.65")
51P. 769
511.C98
511*253
51 C. 769
5"9«936
4519,984
6161, 382
2982..816
6643. 335
7256,316
5875.125
3766, 197
4740 Bn6&
6f'59. 664
6183. 347
4534e425
5162. 824
3883« 798
3816. 518
35C6.288
2519^199
12P1.373
3129,675
2111«588
493*355
498o948
5^0. 9P 8
5f, 1.524
501, 487
500.737
50 C-. 442
501.886
537.Q25
5D8i 958
509.7611
510,945
5 11 «, 040
511.186
510.677
5?) 9. 939
5979,816
3859»921
2974»871
2902.418
5962.2K;
415G0406
55128746
4966. 089
6013«75f'
4702, 656
5D81.917
5026, 875
3 627.. 60 5
2723.992
3187. 323
362?»78&
219?% 775
1427, 263
2846a838
495ffl 12C
499, 238
501,122
501.399
51)1,438
505,737
500.341
502.493
537.349
5H9.125
509.766
51UD52
5110n76
511,082
51C.57C
510.034
6725o464
5848.3D8
676D.871
2931.037
6051,398
2694.489
4632.312
6499*613
5967.421
2139.C97
19799243
4445o445
1662,454
348C, 351
3431, 094
2712a 884
1178.839
1706, 212
1778,112
496,147
499, 5fl!6
501.256
5UU387
501,317
500,619
503,649
503,383
507*623
509,183
509. 900
511,082
511.189
511,015
510.543
510.0*3
4487.339
3224.782
7236.601
3847.717
7635.425
5788.476
3936.839
4273*535
592D.7C3
4582..531
1959.093
2977.059
1644.110
233&.179
1375,295
3135,385
1168.D23
2733a 521
1768.080
496a 845
499,790
501.481
501.365
501,228
500. 667
5^3,829
504.334
537.897
509.312
509,924
511.064
511.399
511,076
519,521
51^.022
67 2?. 61 7
6898.820
4727.816
3647.945
6418.589
3130.831
7065.171
3538,553
2255,749
5443,476
1938.940
1778.650
1625.914
1486.542
2689,957
3103.937
23S6.Q22
216C.58G
1032.739
4970354
499*991
501,554
501.375
501,061
500,719
531.042
505.124
507,986
509.406
510,107
511.028
511,357
511,156
5111,339
509.973
3054.296
6777«746
7908,640
3789,301
7223.468
6886,033
6485.597
5879,285
5331,382
2079,610
2338,230
4176«777
1607,875
1469,986
1345.624
1237.375
1147,303
2702*813
2868,579
497*848
500,277
501,746
501,311
500,911
5D0«,868
501,^94
505.739
508.175
509,507
510.229
511,003
511,375
51U012
5 It ,198
509.936
-------
5r9C97!
5"0
5^9
939
125
571
co
5^6,977
513,161
532,911
5T 1,941
51-"?, 646
499, 323
498, 6ir-
496.818
494,568
494.5;r>4
495,912
, 495,928
Ii 497,491
499,951
497,4'"'
495,3V.
495, 361
493,514
492,343
%' 1245
1 46794^ 2, f<
137C~'82?,0
1 '33 6781, C
76-R849, •'•.
6851409,n
5969643,."
8196027, ''•
16636713.
1^226685.
8318355 =
12844521.
5^9.836
5f 9,162
5^8, 333
506.245
502«764
5Mi P32
5-;; ,469
409,396
498,418
494,3^3
495.:'95
496, :H 4
497.199
495,^29
495.^74
493^282
492»264
-UT1245
1211^548, «
12942384,:
5r, 2.67^
51 1. 67'"
5i ' -,289
499*323
498.3''<2
496.(;ni
494.If 2
495, 553
495.729
495,979
498.415
499,421
497.(68
494,751
494.867
49 3.111
5"9. 823
5:9,028
5;j8. 153
5"'5, 885
5 .3,578
51,13. IC3
5"2»578
5^1,518
^99,323
498.119
495.632
493.803
495.9Q9
495.742
495,983
4980787
499.146
496, 799
494,641
494.678
492,834
5 J9, 686
5J8.946
5;':7, 937
5'.5. 572
5~,3.511
5^3.212
502,603
5*?la 341
499.902
4990 241
497,979
495, 355
4933 547
496,159
495.557
495.94'
499,18?
498.835
496.580
495.193
494,446
492.551
509.586
518*863
5D7.693
5f,>5.279
503.301
503.206
532.457
5 01,143
499.854
499.149
497.650
495,144
493,431
496,120
495.385
495,912
499,445
498.372
496,33a
495,550
494.206
492,264
509.491
5'} 8. 799
507. 501
504.947
5D3.28J
503.130
502,280
5IH.021
499,686
498,951
497.372
494.919
493.16S
496.092
495,723
496.229
499.671
497.994
496.034
495.650
493.983
492.163
509.357
508.684
507.26G
504.617
593.225
503.030
5C2.136
SC'0.829
499.521
498,735
497.193
494.794
493.105
496,050
495.894
496.757
499,823
497.671
495.635
495.611
493,754
492.380
9761R340)
6728-381,0
64834?4aC
6189835,')
123552!»6. ^
15413425,"
l.->948424. ^
7776561.^
1"J92123«"1
1.474^618.
12122780 .
8734274,
1*>226685,
9847465,
66;,5752,
6728^82.
12783357,
13823151^
If 764931.
6373328,
8685342,
'"45212688s,
9859698.
7645548.
9272520.
8049233,
8012533.
0 6972739,
017615328.
f11988219.
C 9125725.
? 7596616.
0 7388658.
028991888,
0 9847465,
0 7241862.
C 837952C,
0 7046136,
t> 851408U
Olf'887261.
G26545328,
0 8636411*
C 7486519,
3 722963C,
019939568.
010911726,
S 8171561,
C 8134862«
•3 722167^.
9 8318355e
J 9541643,
028624912.
01388726C.
C 9541643.
0 7780109,
0 7462054.
^16636713,
Oil? 177753.
0 8196027,
0 7890205,
0 6923806.
D 69238^7,
0 7951369,
325566672.
D10838329,
0 9052328.
0 9248055.
0 6728081.
^14679452. C
010520274.0
01C52&274.0
0 7951369,0
0 7241862.0
«•> 6581287.0
3 8196026,0
Di:)226685,0
0 9^52328.0
011743562.G
C 7119534.0
-------
5896?46»n 614i1903a
55>4794
5?;' 9 "^67
-,
a
55f-4794,
5 ft? ? 848
4293738
429373R
6B74876
5578191
8269423
631216?
499 m 2
433">438
5162273
4 5 >" I 69 7
6361-^94
6; 18573
29736544
1
1
2
1
1
">9484?5
i
a
*t
*i
•»
o
i>
T»
"»
9
1
0
*
,
a
3896549,
[551232
746? 553
9296986
6 8 53 41 f
1621233
719293T
241637S
1223287
1223287
17126D3
8 563?- 15
8J"737f:f
8318358
7339727
7829042
8563^1.
,
s
D
fl
3
,
a
o
»
.
&
9
**
,
«
,
f- 66T5752
0 5186738
C 5749451)
r- 5932944
0 5357998
?• 33*2875
" 545586?
r 5945177
*• 7963602
J 6!''675f'6
>> 4833149
"i 433f»438
0- 52112C4
r- 4893149
•: 6?63232
U 55*4793
F 16RP, 136C
3 959^576
C,1113191R
C1671tnil
f 8196^26
C- 9345917
e> 7192931
H lf^53251'7
n 7657 7 8^
P11254247
0 11-0959
^ 2935891
C 9663974
^13333838
0 2446575
111743564
U12232879
r. 2446575
6 73397?
,
»
A
«
a
,
a
.,
,
•n
13
*
<9
if>
O
1
W
O
-a
0 63733279
D 6226533,
0 5f'.f;-3244,
- 76G8849,,
r= 5125574,
D 41714r8«
0 3376273»
'' 5272368,
0 8391753a
C; 7608849B
0 74131 22o
••' 4771)82?.
n 34986f-l.
01t>D3f1957.
0 5627122*
j 5822847*
1 5565958,
C123552D6,
0 8954466,
- 96762D4,
•"'149241ir-.
S 8991163,,
0 88G767U
f? 7767877.
0 9125726»
C P6241770
r.
0 1223287*
0 3180548,
5 978630,
01321151D,
^ 905233fto
1J 4893151«
010764934,
*) 1467945.
7 n
*">
^
r
.-
r>
?j
L.-
o
01
D
0
*?
f-
r-i
r.
r
0
5
r?
0
69238D7,
5945178*
5492560,
7339725,
4146943.
3767724a
3865588.
5198971,
3015782*
8^9232.
6f-91971,
43f;597i.
344967! ,
12Q5315,
6752547.
55G4793,
5431396«
9969794,
62387660
9431548.
013333836,
••..*'
0
0
•J
n
!?
0
3
01
01
'-.•••
8832137.
87j98D8»
8318355.
8538547a
6544588,
9786 3U
1957260,
97863f'»
2722194.
C764934,
8D737D00
G10D3D961.
n
D
D,
733972.
0 9125726=,
Q 7975834,
0 515004C:0
C 44^3835,
H 3731024,
0 44f?3834»
C 411D246,
C 8269424,
015046440,
C 8807671*
f? 4770820,
D 4648492,
P 77f6712»
C 8514083,
0 7^46137B
' 5284602,
0 5C'f"'3245-,
C 7095t67,
C- 672808C.
t 91134930
('123C6274,
0 6801479,
D 8440685»
G27279296,
C- 7315259.
C. 6165369,
3 1957260,
0 19572608
3 318C548,
0 9C5233Ge
01(5275618,
C 11743564,
0 1957260,
0 0.
7 4893153
C
V
D
V-
o
?
0
0
fi:
D
0
D
78 2 9041 o
6434492*
466T7240
5504794^
3792191',
4966547,
417141D.
8281&57»
10826095^
7767876«
4342 67C9
4991^12^
'j
0
9
A
'")
0
0
0
0
;j
0
0
51D6426D2.0
r?
D
Cl
Cr
0
0
0
•'^
0
D
D
L
^
*}
**,
a
D
0
0
0
A
;*»•
1
6899342«
60797399
51378CI7*
7486521^
8073697*
81348639
9419313,
11804725*
7254096,
7D95066.
1963376S«
7168464^
61653690
1467945S
195726G,
1467945^
8563015,
9296988.
66C5755*
4159179,,
0.
733972.
0
v)
0
0
*>
0;
y
0'
i S
D
D
''l>i
•^
•o
;j
0
0
(J
.n
D
0
7
6972738.
5039944.
5321300.
5431396.
3559765.
4257U40.
47218900
5884012.
9321451.
6831478.
4991012,
5688285.
6752548.
6287698.
9052329.
5015478.
6544588.
8933001.
023184240.
Q18JJ55712.
010948425.
7474288*
6238766.
014924110.
7706711.
9174657.
2691233.
3865589.
6605755.
8318358.
010642605.
5382467.
1223287.
0.
7584385.
0 5015478. 0
0 5076643,0
D 6140903.0
O 5773917.0
0 411D245.0
0 433D437.0
0 8477384.0
0 5529259oG
D 9431548, «
0 7229630.0
D 4819752*0
0 4379370,0
0
0 540693J.O
u 8440683.0
0 6581287.0
C' 99^8631,0
011865891.0
318 71 62 88*0
020257632.0
015397945.0
0 6617985*2
0 6728581,9
012905685,0
0 7278560,0
018716288*0
0 1345616.0
& 2446575. €
0 5627124*43
0 9296988.0
9 1364263 5. Q
0 2691233.0
a 6361097*0
0 183493.2
0 3302877.0
-------
2691233,
685^412,
97863! ,
22^1918,
8f;737;)f .
1'J 523276.
13525276.
44S3836,
8196029,
7829C42,
1 5535757.
1688136?.
9')5233Jk.
97863'>3o
10397947.
, 9541646.
II 1C 397947.
YI 7829^42,
11743564,
1467945,
13?r r 825,
12.232879,
611643,
4159179,
15413428,
13945482.
17CC-3696.
1223287°.
12477537.,
22
"s 744657.6
r 2446575* ^
3 1 467945 „<'•>
r 6116439,'?
r. 8 44? 6 86.-'
f> 8318358.0
C 9786 3r? Or;
f 636K97,':
'? 97863r3.n
• 34?52H6,0
016392^58,"
"16636716»^
•) RR'" 7673,0
r 1^''3">961« )
.; 81*'737e^.O
T 92969R80:"5
"> 4159179nH
• i'^764934^ '}
(.11988222, "-
r, 489315»1
ri?232879.H
n 85S3fl.5«0
9 '.'' n "t
"1565fl.J85-,r?
-14924113, ;"
^1419^14!?. -3
r'17ff 3696,0
^12232879,'')
,nl2477537»"
733972,7
2446575, (
1712603. o
8(. 737? *-.?•>
8563015, '.>
8196C29. j
5627124, 3
F563f'*l 5« G
33)2877, .j
110, .9591.0!
1737C688, )
156 58085. T
9541646, 1
9541646.":
9541646.0
8807673, G
1" 5 2'>2 76, }
8563015. •>
1^764934.0
.**H f i
n 5 201 2 76.1
9786303,3
4648494.^
15658C'85.r:
12966852,3
13945482»0
14679455.0
12232879,0
443»'
r'- a
r
»
9
15.24
, "; 73156608
r* ^0"!
9 3"?6 3;)6
1 5 a 2 4
C')C« .T86
58
In ^
3914521.
685T412.
3180548,
6116439.
97863t3.
88t 7673.
2935891.
4159179,
77'* 6714.
010030961.
)16881360.
U516877D.
7829042.
8807673.
9541646.
i 11254249.
'30961.
•> 122 3 2 879,
?114190140.
12232879.
r-j
5871782.
86853^4.
195726%
£15658^85.
314924113.
97863C3.
012232879.
012232879,
0 6361C97,
0 2446575.
r 5871782.
r 7339727.
1.11 520276,
P 6361C97.
0 526C'138a
C 7217398»
011254249=
016881360,
t.1516877?.
0 7095070,
ff 4159179,
011> 275618.
U11621235,
C 9786303»
C11254249,
011254249,
r. '"- 3
211988222,
C 9419317.
C 489315,
517126016,
t'1321151Co
5137809.
66D5755.
7951371.
1345616.
3669863.
4893151,
4770823.
',1' 7C9507D.
C' 7584385.
^11865893,
316514387,
jl68813639
} 8563015.
5 88^7673,
^1^764934.
H11498906.
e 9296988.
!J14679455.
'10397947.
C 8807673,
j 8318358.,
131D52"'276.
1 9541646.
1 3547535.
C12232879,
C12232879,
C 92969883
012232879.
01211u550.
.? 5627124.
3 4648494.
f) 7829042.
0 8318358.
J 1957260.
0114989^6.
0 7095C70.
J 7U95D70.
1 6361097.
) 11743564.
016881360.
fJ 9^52330.
'! 9296988,
01J275618.
0 8073700.
•^11449975.
D11254249.
C15902743.
010642605.
311988222.
0 9296988.
31^764934.
3 7339727,
? 3669863.
01737,3688.
•513945482,
51023D961.
012232879.
01211J55-J.
C 4893151.^
? 7339727. C
0 7829042.0
0 4159179,0
0 9908632.^
•3115:9591. C
0 7584385.0
0 8563015."'
?.< 7829042 »v
0 14 193 140. C-
CIL'275618. 0
OICC'30961.0
Q 9296988.^
0 9786303,0
019275618.0
011254249,0
012232879,0
0 9'?52330.D
010275618.0
112232879. 0
J137C3825.0
a 4893151.0
0 55&4795.0
014801784,0
314924113.0
^12232879."
012232879,0
012232879.0
1.0
396C24
396=24
-------
783279,9 169968*", 4249199. 7243872, 10643231, 14487744. 21286463. 3QD27672.
16u77,34 23480,33 28211,85 34552.62 41338.27 43259.17
en
i
1 770,, 2 8
45737,56
8,193
1,618
2«653
6,6^3
9,852
"> , 68'*
3. 187
^. 788
2, 281
% 533
2,899
1, 586
0.979
1,848
I,f58
UC43
2,386
1. 359
1*464
1,257
%623
1,345
?,25C
"U691
% 741
r>.729
1.626
1. 315
2*737
3»R89
1, 394
1.780
1^863,07
2.*r.s
H. 992
3«775
1*949
3.328
3«?38
20649
1,355
1.862
f - * 1 5 2
f ,721
1,550
0.887
* e 2 5 1
0.8^2
0.639
1,435
1» 5i 6
t « 8 9 1
0.849
1.211
U215
<%842
P.646
D.592
?»229
f o 2 4 1
1,378
1 <,nOO
2a956
1«261
2«, 11"
1»921
1.C35
6.938
6.13*
5.;, 95
4.248
U969
1.731
1 . 942
4. 1 2 0
1,734
3,277
1.967
l.T-81
9. 534
G, 418
^•,929
U. 909
1.181
1.628
0, 824
1.028
1, 151
0.819'
0.328
C, 121
%439
2,961
1. 324
1,661
1.455
2. 028
3.429
2.853
4,893
5.487
6.488
5.726
D.756
0. 43 C
1,666
2, C54
2.097
4.391
1,293
1,867
G, 586
1.129
1.3?9
1,157
1,291
?9885
1.176
0.717
Oo 87ii
1.134
0.456
1.130
0.635
1,881
2. 234
f',962
2,289
D0 763
6.9H
1,644
1, 839
3.184
3a 390
5a434
3=322
0,721
1,2? 5
0.259
&. 844
29774
0,875
1,564
1*413
1.561
1,676
0.727
C. 909
2«238
5.724
U313
3.654
1.564
0. 346
1.802
4,610
0,439
1.934
0.555
1»28*
0. 500
7,768
C'a889
29 6*3
8»166
*3515
3,313
2»*55
0.972
4,386
2.256
1.101
U630
2.457
t-,382
1.089
C3163
1,735
2.642
1*662
0«8*3
0.462
1, 69>0
0. *97
1,488
2a *00
2.178
1.279
D, 65!'j
2,138
0.956
2.372
0,*71
5, *66
'),642
3* 13*
6.91*
2.67*
2.67*
*a!26
2.521
2.721
C.*5*
v»2l7
0.903
2.107
1,029
1.268
2,760
1*550
1.517
1.027
n.777
1.483
0.733
0.390
1.65*
0.7*5
2.115
3.110
0,328
1.951
1.208
5,*38
1.273
3» 086
3.361
U.I 44
6,889
8.37*
1,958
2.821
*, 52*
2.672
1.5*1
D,292
1*006
1.478
1.094
D.721
2»133
1,329
2,137
0,482
0,789
1.681
2, 076
0.71*
0.92*
1.705
1.302
1,121
*,608
0.515
0.381
7.725
0.713
-------
1, 655
1. 729
1 *% 2 H 3
•'"•.332
3.984
rU 59 u
•j.592
5359, 863
4751, 516
6C59.133
4844, 613
5563,699
5G 6SS 1^2
6647, 836
7476.668
6121, 16*
6784, 8>?9
7229, 883
7679, 789
6979.^82
7456,379
7371, 131
8397.141
7948, 527
8385, 234
7111,285
7484, 746
84563773
8*51, 590
7172.234
7612.547
6904, 562
7492.859
6446,547
6f-89«, 328
5683.039
1.234
I «414
12»457
i'J, 569
1.179
0*966
0.745
5688*148
5168a273
5946,841
5268,578
6222.957
6249a648
7204=191
7381,379
6517.695
7338*683
73^3.213
70 2 4« 7 58
7283.1^9
7296,387
7^24. VI 2
8?36,129
77 8 2. 41 ^
7933,77:.'
743^93^
7312*430
7911«742
8111,734
7778,113
750r-*578
73340477
742 6« 094
7223*926
6254.922
5787.1 59
1.435
1.010
6.32fl
1, 388
0, 346
6, 726
6561, 992
5329,066
5 86 w« 2 54
4951*293
6089.039
6 99 1). 77?
7469,715
7621.145
72P2.363
717'i. 375
7506.473
6 53 C. 402
7931,371
7184.183
7124. 5HO
8159,238
7603o 152
8112,770
7484. 81'5
8139,949
79 li). 773
8672»945
7246, 906
762.1. 980
7337,852
7384, 02C
7535,242
6492« 477
6088.355
0.875
fi,3^3
0,787
3. 839
D0423
6, 212
6726.355
5536S418
5818,348
4768,641
5304. 258
7127925r
7473, 594
7414. 383
7320.734
6572,730
7809,570
5955,070
8C76.301
7413. 969
7408, 723
8148. 437
8224.410
8157, 992
7619B 887
7657, 551
8286, 695
8511*437
7823.852
7872a359
7252,727
7509.039
7V' 9 3. 656
5860.102
6770, 586
•f>. 7D3
0. 415
4. 82 5
3. T""' 2
0.32^
3.251
4716.891
6547S 926
6395, 285
49- lo 27:;
5186. 629
66573328
73520336
7334.086
7075.352
6812.236
7539, 875
58u4s23D
8137. 871
7594» 6D2
7457. 172
8150, 547
8191,645
78993 594
8073.980
7533.656
8054.516
7937oT 74
8'?-33. 828
7780.195
7731, 4!>6
7437.426
6696,, 422
6032^348
7805. 184
!?. 960
f;9474
2^617
4,136
lo t>07
G . 3 8 8
4418.652
6811,750
68773 48'"^
5 651 . 8'"- 5
5403.684
64U9.637
7115.891
7499,074
6282.3-87
7898.969
7616.016
61499875
79C'6o 078
7555.770
79C7.836
8283o 824
8113o930
7615.902
7767, 598
8^23.031
8446,711
8026,980
812rf,164
7038a 840
7522«S94
6908,430
7397,156
6270.711
669Da336
1.631
^.381
1.822
r.84i
3,952
Io767
4471 0 2 58
6803.648
68 -H. 520
5523.246
5479. D39
5656.371
6959.629
7236.093
5659.367
7567.121
767C.324
6748,168
7472. C 82
7352.023
8425,734
8697. 852
8767.758
7695.6C9
8149.699
7916,973
79434oQ47
8064,906
7847.305
7214,680
7469.223
6482.023
6460. 887
7118,867
5657.578
10459
2.876
00319
1.35 5
1.731
7*885
4479.535
6515.793
5295,723
5532, 844
5328,777
5769. 660
7336,344
6135.754
6227,539
7732.457
7732, 340
6984a 199
7157.809
7167,137
8442,891
7919,339
859D.437
7329,574
7942.430
8248,227
8006.719
7314. 211
7812.555
7055.676
7047.305
6458,586
5851, 910
6726,227
6612.137
-------
6? 12, 4P.-2
6428,273
469 '•„ 746
*455n687
5780,593
5342a922
61?.?«422
5075.309
5729834G
4974*98-?
*STOP
5745*551
573<5e641
56TU93!";
6696,343
5732,281
5226,633
66?< «^82
4^25,215
6^06a?2^
5551«!"?66
5771,695
5842.367
6322ol5?
62130141
5407.797
51^3.184
7186»574
4758,594
4583* 145
6197,914
5840.879
6285,121
5854B039
5575,625
47t'?6, 953
7554.059
54-31,308
4203. noc
7323,012
5997,723
6800, 262
5677,031
5514*562
534&.715
6732, 836
5605.648
4903, 75D
7402,449
6459.695
7113,973
59470832
64Q40789
4976,117
5731,562
5335.344
5178.578
7231.988
57?. 6.117
7264.023
6686.859
6588,035
4839.344
5346.113
5526.586
5878.484
6738.133
4758.484
7015.832
6690,953
6673.324
55S0.930
5550.945
5338.875
4487*641
00
I
-------
ACKNOWLEDGEMENT
This investigation was supported, in part, by the Water Quality
Office, Environmental Protection Agency, under Research Grant No. 16130 DJH
entitled "Thermal Stratification and Reservoir Water Quality". The project
officer was Mr. Frank Rainwater, F.W.Q.A. Pacific Northwest Water Laboratory,
Corvallis, Oregon.
A portion of the laboratory experiments were carried on in 1967 by
Dr. Wayne C. Huber, then a graduate research assistant in the Hydrodynamics
Laboratory. Computations were made at the M. I.T. Information Processing
Service Center. Our thanks to Miss Kathleen Emperor who typed the report.
-119-
-------
Bibliography
1. Anderson, E. R., "Energy Budget Studies", part of "Water Loss Investigations -
Lake Hefner Studies", Technical Report, U.S.G.S., Professional Paper 269,
1954.
2. Bachmann, R. W. and Goldman, C. R., "Hypolimnetic Heating in Castle Lake,
California", Limnology and Oceanography, Vol. 10, No. 2, April 1965.
3. Bella, D. A. and Dobbins, W. E., "Difference Modelling of Stream Pollution",
Proc. ASCE, SA5, October 1968.
4. Brooks, N. H. and Koh, R. C. Y., "Selective Withdrawal from Density Strat-
ified Reservoirs", ASCE Specialty Conference on Current Research into the
Effects of Reservoirs on Water Quality, Portland, Oregon, January 1968.
5. Dake, J. M. K. and Harleman, D. R. F., "An Analytical and Experimental
Investigation of Thermal Stratification in Lakes and Ponds", Technical
Report No. 99, M.I.T. Hydrodynamics Laboratory, September 1966.
6. Dake, J. M. K. and Harleman, D. R. F., "Thermal Stratification in Lakes,
Analytical and Laboratory Studies", Water Resources Research, Vol. 5,
No. 2, April 1969.
7. Elder, R. A. and S. Vigander, "Capability and Potential of USAGE's Current
Meter for Ultra Low Velocity Measurements", Isotopes in Hydrology, Proc.
Symposium, November 1966, International Atomic Energy Agency, Vienna, Austria,
1967.
8. Elder, R. A. and Wunderlich, W. 0., "Evaluation of Fontana Reservoir Field
Measurements", ASCE Specialty Conference on Current Research into the
Effects of Reservoirs on Water Quality, Portland, Oregon, January 1968.
9. Elder, R. A. and Wunderlich, W. 0., "The Prediction of Withdrawal Layer
Thickness in Density Stratified Reservoirs", TVA Division of Water Control
Planning, Engineering Laboratory Report No. 0-6781, March 1969.
10. Goda, T., "Density Currents in an Impounding Reservoir", Eighth Congress,
I.A.H.R., 1959.
11. Howard, C. S., "Density Currents in Lake Mead", Proceedings, Minnesota
International Hydraulics Convention, September 1953.
12. Huber, W. C. and Harleman, D. R. F., "Laboratory and Analytical Studies
of Thermal Stratification of Reservoirs", M.I.T. Hydrodynamics Laboratory
Technical Report No. 112, October 1968.
13. Hutchinson, G. E., "A Treatise on Limnology Vol. 1", Wiley & Sons, 1957.
14. Kao, T. W., "The Phenomenon of Block in Stratified Flow", J.G.R., Vol. 70,
No. 4, February 1965.
-120-
-------
15. Keulegan, G. H. , "Laminar Flow at the Interface of Two Liquids", J.
Research for the National Bureau of Standards, Vol. 32, June 1944.
16. Koh, R. c. Y., "Viscous Stratified Flow Towards a Line Sink", W. M. Keck
Laboratory, Report KH-R-6, Caltech, 1964.
17. Kohler, M. A., "Lake and Pan Evaporation" in "Water Loss Investigations,
Lake Hefner Studies", Technical Report, U.S.G.S., Professional Paper
269, 1954.
18. Mortimer, C. H. , "Water Movements in Lakes During Summer Stratification",
Phil. Trans, of Royal Society of London, Vol. 236, 1952.
19. Orlob, G. T. and Selna, L. G., "Progress Report on Development of a Mathe-
matical Model for Prediction of Temperatures in Deep Reservoirs - Phase 3:
Castle Lake Investigation", Water Resources Engineers, Inc., Lafayette,
California, January 1967.
20. Orlob, G. T. , "Mathematical Models for Prediction of Thermal Energy Changes
in Impoundments", Final Report to FWQA by Water Resources Engineers, Inc.
December 1969.
21. Rohwer, C., "Evaporation from Free Water Surfaces", U.S. Department of
Agriculture, Technical Bulletin No. 271, December 1931.
22. Ryan, P. J., "Thermal Stratification and Overturn in Lakes and Reservoirs",
M. Eng. Sc. Thesis, University of Melbourne, June 1968.
23. Slotta, L., Personal Communication, August 1970.
24. Stone, H. L. and Brian, P. L. , "Numerical Solution of Convective Transport
Problems", Journ. A.I. Ch. E., Vol. 9, No. 5, 1963.
25. Sverdrup, H. U. , "Oceanography for Meteorologists", George Allen and Unwin
Ltd., London, 1945.
26. Sweers, H. E. , "Vertical Diffusivity Coefficient in a Thermocline", Limnology
and Oceanography, Vol. 15, No. 2, March 1970.
27. Swinbank, W. C. , "Long-Wave Radiation from Clear Skies", Quart. Jour, of
the Royal Met. Soc. of London, Vol. 89, July 1963.
28. Stolzenbach, K. D. and Harleman, D. R. F., "An Analytical and Experimental
Investigation of Surface Discharges of Heated Water", Technical Report No.
135, Ralph M. Parsons Laboratory for Water Resources and Hydrodynamics,
M.I.T., February 1971.
29. Tennessee Valley Authority, "Fontana Reservoir 1966 Field Data, Part I,
Temperature and Flow Data", Report No. 17-91, TVA, Division of Water
Control Planning, Engineering Laboratory, Norris, Tennessee, March 1969.
-121-
-------
Definition of Notation
Representative units of variables are given in m, kgm, days, kcal and C.
Small Roman Letters
a constant in evaporation formula (m/day - millibar)
b constant in evaporation formula (millibar )
c specific heat (kcal/kgrnSc) = 0.998
d depth of layer for entrance mixing (m)
d displacement of thermocline due to wind (m)
e saturated water vapor pressure at air temperature (millibars)
3.
e saturated water vapor pressure at temperature of water surface (millibars
s
2 10
g gravitational acceleration m/day = 7.315 x 10
h mean depth of reservoir (m)
h depth of mixed layer (m)
i subscript denoting inflow
i slope of water surface
S
i slope of thermocline
j number of grid point
jm number of surface grid point
2
q flow per unit width(m /day)
x entrance mixing ratio
s slope of reservoir bottom
t time (day)
u. inflow velocity (m/day)
u outflow velocity (m/day)
x horizontal distance from outlet (m)
y elevation (m)
-122-
-------
elevation of reservoir bottom (m)
elevation of bottom of mixed layer (m)
elevation of outlet (m)
y elevation of water surface (m)
o
Roman Capitals
2
A horizontal cross-sectional area (m )
2
A^ average area of bottom element (m )
2
A surface area (m )
B width of reservoir (m)
2
D coefficient of thermal diffusivity (m /day)
2
D coefficient of numerical dispersion (m /day)
2
E evaporation mass flux (kgm/m -day)
F densimetrie Froude number
2
H heat flux (kcal/m -day)
L length of reservoir (m)
L length of reservoir at elevation of thermocline (m)
L latent heat of vaporization (kcal/kgm)
N ratio of heat conduction coefficient to mass transfer coefficient
(kcal - millibar/kg°C)
3
Q reservoir flow through rate m /sec
3
Q. inflow m /day
3
Q outflow rate m /day
out
R Bowen ratio
2
S average area of surface element (m )
area
T temperature (°C)
T air temperature (°C) (unless otherwise specified)
a
-123-
-------
T. inflow temperature
T1. mixed inflow temperature (°C)
T temperature of mixing water (°C)
T temperature of surface water (°C)
s
T temperature of surface water (°K)
w
V vertical advective velocity (m/day)
¥ volume of reservoir (m )
W wind speed (m/sec)
Greek
2
coefficient of molecular conductivity (m /day) = 0.0125
3 fraction of solar radiation absorbed at the water surface 0.4 - 0.5
6 thickness of withdrawal layer (m)
A increment
e normalized density gradient (m )
n radiation absorption or extinction coefficient (m )
2
v kinematic viscosity (m /day) = 0.0864
3
p density (kgm/m ) = 997
3
p reference density (kgm/m ) = 1000
a. standard deviation of inflow velocity distribution
a standard deviation of outflow velocity distribution
2
heat flux (kcal/m -day)
2
longwave atmospheric radiation (kcal/m -day)
3.
2
heat loss flux due to evaporation (kcal/m -day)
2
<|>. short wave radiation at element j (kcal/m -day)
-124-
-------
2
j>L total surface heat loss flux (kcal/m -day)
2
j) long wave radiation from sides of laboratory flume (kcal/m -day)
2
net short wave radiation at water surface (kcal/m -day)
2
shortwave radiation reaching water surface (kcal/m -day)
s
2
shortwave radiation transmitted through water surface (kcal/m -day)
Yysur
ij) relative humidity
-125-
-------
1
Accession Number
w
5
rj Subject Field & Group
05G
SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Urbanization
Cambridge, Massachusetts 02139"
Ti(/e
TEMPERATURE PREDICTION IN STRATIFIED WATER: MATHEMATICAL MODEL-USER'S MANUAL
(Supplement to Report No. 16130DJH01/71)
10
Authorfs)
Patrick J. Ryan
Donald R.F. Harleman
16
Project Designation
Grant 16130 DJH 04/71
21 Note
22
Citation
Environmental Protection Agency, Water Pollution Control Research Series, 16130 DJH
01/71, pp. 125, January 1971, 6.6 fig., 29 references.
23
Descriptors (Starred First)
*Reservoir temperature distribution; *Thermal stratification in reservoirs;
*Reservoir water quality.
25
Identifiers (Starred First)
mmmmmw*^m&mM^^
27
Abstract
The annual cycle of temperature changes in a lake or reservoir may be quite complex,
but predictions of these changes are necessary if proper control of water quality is to be
achieved. Many lakes and reservoirs exhibit horizontal homogeneity and thus a time-dependent,
one-dimensional model which describes the temperature variation in the vertical direction is
adequate. A discretized mathematical model has been developed based on the absorption and
transmission of solar radiation, convection due to surface cooling and advection due to in-
flows and outflows. The mathematical model contains provision for simultaneous or inter-
mittent withdrawal from multi-level outlets and time of travel for inflows within the reser-
voir. Heat transport by turbulent diffusion in the hypolimnion is neglected. Good agreement
has been obtained between predicted and measured temperatures in both the laboratory and the
field. Field verification consisted of the simulation of the thermal structure of Fontana
reservoir during a nine-month period. Criteria for the applicability of the model are given.
The mathematical model is a predictive one, since the required data is that which would
normally be available before the construction of a reservoir. Emphasis has been placed on
a detailed explanation of the physical basis for the mathematical model and on the computer
program inasmuch as the report is intended primarily as a user's manual.
)naTcTR.F. Harleman
Institution
Massachusetts
Instit.ut.p of Technology
ICUMENT. TO: WATER RESOURCE*>*C
WR:I02 (REV. JULY 1969)
WRSIC
SEND WITH COPY OF DOCUMENT, TO: WATER R ESOU R C E*>*C I E N T I F I C INFORMATION CENTER
U.S. DEPARTMENT OF THE INTERIOR
WASHINGTON, D. C. 20240
* GPO: 1970-389-930
------- |