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                                                                    Pae

 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  = 4C (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     2o                                 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  JJ
           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. > 100C,
                 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 )-224) 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/2D-DYSUR)*-( A(JM)-A( JM-1) ) /DY
       GO TO 86"-
  859  AYSUR = A(JM+(DYSUR-QY/?
       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/2n-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';)/20
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
    5H 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,UCT( 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) )/DYMA(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*2D
       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/23*(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,ll
       1  CONTINUE
     9GO  FORMAT (2CA4)
     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='F72,
       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='F62,17X,'BETA=F52)
     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,15XfKQ=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=F62 )
    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=E125,14X,
       1'INSOLATTON  FLUX=El?.51
    916 FORMAT  (  HEAT LOSS  FLUX=E125,15X,RADIATI ON FLUX=E12.5,
       116XIINFLOW  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 lF6a 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='F62,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='F52)
    926 FHRMATC  ENERGY RATIO=F7.3, 5X.'ENERGY  ST3RED='E12.5, 5X,'ENERGY
       1INFLOW='E1295, 3X,'ENERGY OUTFLOW=E125)
        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=l13587F-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(EVAP3t3,4
        3  EVAP=0C
        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)-4C)/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**025
    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)
       lnY)/(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,EOJM) 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,DTDY
       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 ITI< 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
49249
1,96 4.
3>"6 3D 6
1 = 5
7. 58
6n ;"j 6
If; .84
9,, 4r
9n91
995
13, 87
1568
13,26
1259
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
1575
13, 69
2 2
?oT
C
356 3r: 6
I.1)
7,24
6,67
ir83
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
19OD
1643
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
1967
20'e31
23a63
22.84
21.33
18,31
1762
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
2081
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
1283
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
14753
15.99C,
16884
i824r-
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
884
6,24
582

8.752
5834
1%325
2.808
7a!X 9
14, 036
15. 785
17,253
16,124
16r77
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
1425
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
0592
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
9435
-0.226
5,164
1,718
fc.7^7
C.674
0,730
r.a852
Ct.,618
0*568
0.757
ft9C 7
^643
.95t>
C818
"856
::-.775
C.9C5
C771
0886
C827
;Ba47
0.707
5.816
il873
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
0755
t. 857
')670
0.789
C. 861
0.671
0. 77 1
D744
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
0668
t858
9, 798
D778
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
6126
7. 803
0,443
3,067
3,808
Da409

0, 989
0.73G
0.669
0.569
0.597
0. 795
C796
0929
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
4788

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
0837
0,945
U 978
r s 8 6 2
0,800
'" 7?'"f
rr8f 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'n771
501.298
500ffl856
5^792
501.119
506,218
568,382
539.522
51D555
511* 101
511.308
510, 863
510, f 52
8427.151
&r?7. f^7
7423613
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^
5rr.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"9936
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
2111588

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
3859921
2974871
2902.418
5962.2K;
415G0406
55128746
4966. 089
601375f'
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
6777746
7908,640
3789,301
7223.468
6886,033
6485.597
5879,285
5331,382
2079,610
2338,230
4176777
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
              502764
              5Mi P32
              5-;; ,469
              409,396
              498,418
              494,3^3
              495.:'95
               496, :H 4
               497.199
               495,^29
               495.^74
               493^282
               492264
           -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"2578
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?46n 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 433f438
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, 4893151
ri?232879.H
n 85S3fl.50
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
139454820
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
2653
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
3R89
1, 394
1.780
1^863,07

2.*r.s
H. 992
3775
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
1261
2, 11"
1921
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
2238
5.724
U313
3.654
1.564
0. 346
1.802
4,610
0,439
1.934
0.555
128*
0. 500
7,768
C'a889
29 6*3
8166
*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
08*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*
v2l7
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
2133
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
12457
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
7911742
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
8672945
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
7522S94
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
7186574
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/kgC)
                                      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

-------