University of Minnesota
       St. Anthony Falls Hydraulic Laboratory
             Project Report No. 329


Water Temperature Characteristics of Lakes

        Subjected to Climate Change
                      by

                 Midhat Hondzo

                     and

                 Heinz G. Stefan
                  Prepared for

   ENVIRONMENTAL RESEARCH LABORATORY
  U.S. ENVIRONMENTAL PROTECTION AGENCY
                Duluth, Minnesota
                  August 1992
              Minneapolis, Minnesota

-------
                                 Abstract
      A  deterministic,  one  dimensional,  unsteady  lake  water  temperature
model was  modified  and validated to simulate  the seasonal  (spring  to  fall)
temperature  stratification structure over a  wide range  of lake morphometries,
trophic  and  meteorological   conditions.      Model  coefficients   related   to
hypolimnetic   eddy   diffusivity,   light  attenuation,   wind   sheltering,   and
convective  heat  transfer  were  generalized  using  theoretical  and  empirical
extensions.

      Propagation of  uncertainty in the lake  temperature  model  was  studied
using a  vector state-space method.   The  output  uncertainty was  defined  as
the  result of  deviations of meteorological  variables from  their  mean values.
Surface water  temperatures were affected  by  uncertain  meteorological  forcing.
Air  temperature and dew point temperature fluctuations  had significant effects
on  lake  temperature  uncertainty.  The method  presents a  useful  alternative
for   studying  long-term  averages  and variability  of the  water  temperature
structure  in  lakes due to variable meteorological forcing.

      The lake water  temperature  model  was  linked to  a daily meteorological
data base to simulate daily water temperature in  several specific lakes as  well
as  27 lake  classes  characteristic for  the north  central  US.    Case studies  of
lake  water  temperature  and  stratification  response  to  variable  climate were
made in a  particularly  warm  year (1988)  and a  more  normal one (1971).   A
regional   analysis  was   conducted  for  27  lake   classes  over  a  period  of
twenty-five years (1955-1979).   Output from a  global  climate model (GISS)
was used to modify the  meteorological data base to  account  for a doubling  of
atmospheric  COz-   The simulations  predict that  after  climate  change:  1)
epilimnetic  water temperatures  will be higher but  will  increase  less  than air
temperature,  2)  hypolimnetic  temperatures  in   seasonally  stratified  dimictic
lakes will be largely  unchanged  and  in  some  cases lower than at present,  3)
evaporative water loss will be  increased by  as much as  300  mm  for the open
water season,  4) onset  of stratification will  occur  earlier  and  overturn  will
occur later  in  the season,  and 5) overall lake stability will become greater  in
spring and summer.
The  University of Minnesota is  committed to the policy that all persons shall  have  equal
access to its  programs, facilities, and employment  without regard  to race  religion,  color,
sex,  national origin, handicap, age,  or veteran status.

-------
                                  Preface
      This  study  addresses  the  question  of  how lake  water  temperatures
respond to climate  and climate changes.   The study is conducted  by model
simulation.    The  chapters  of this  study  are  a  collection  of  papers  or
manuscripts previously published or submitted  for publication in professional
journals.  Each chapter has  its own abstract and conclusions.  Each chapter
-of this study deals with a subquestion of the  problem.
                                      11

-------
                       Table of Contents
                                                          Page No.
    Abstract                                                       i
    Preface                                                        ii
    List of Figures                                                 v
    List of Tables                                                ix
    Acknowledgements                                             xi

1.   INTRODUCTION AND LITERATURE REVIEW                 1

    1.1  Introduction                                               1
    1.2  Previous  temperature  prediction model                        2
        1.2.1  Model formulation                                    2
        1.2.2  Model coefficients                                    5
    1.3  Overview of study                                          6

2.   REGIONAL LAKE WATER TEMPERATURE
    SIMULATION MODEL DEVELOPMENT       .                 8

    2.1  Introduction                                             '  8
    2.2  Model generalization                                        9
        2.2.1  Hypolimnetic diffusivity closure                        9
        2.2.2  Attenuation coefficient                               13
        2.2.3  Wind sheltering coefficient                            13
        2.2.4  Wind function  coefficient                             16
    2.3  Water temperature model validation after
        generalization of hypolimnetic eddy diffusivity                19
    2.4  Numerical uncertainty of model after
        hypolimnetic closure                                       29
    2.5  Accuracy of the  regional model after
        implementation  of all changes                               30
    2.6  Conclusions                                               36

3.   PROPAGATION OF UNCERTAINTY DUE TO
    VARIABLE METEOROLOGICAL FORCING
    IN LAKE TEMPERATURE, MODELS                          37

    3.1  Introduction                                              37
    3.2  Numerical model                                          38
    3.3  First and second moment development                       39
    3.4  Lake Calhoun - application                                43
        3.4.1  First  moment  analysis                               46
        3.4.2  Second moment analysis                              46
    3.5  Conclusions                                               54
                              111

-------
     4.   CASE STUDIES OF LAKE TEMPERATURE AND
         STRATIFICATION RESPONSE TO WARMER CLIMATE       55

         4.1  Introduction                                              55
         4.2  Method of lake temperature modeling                       56
         4.3  Model validation                                          59
         4.4  Results and discussion                                     59
             4.4.1  Thermal energy  budget                              59
             4.4.2  Equilibrium temperatures                            64
             4.4.3  Vertical mixing/Onset of stratification                69
             4.4.4  Water temperatures                                 69
         4.5  Conclusions                                              70

     5.   WATER TEMPERATURE CHARACTERISTICS OF MINNESOTA
         LAKES  SUBJECTED TO CLIMATE CHANGE                 74

         5.1  Introduction                                              74
         5.2  Method of lake temperature modeling                       76
         5.3  Climate conditions simulated                               76
         5.4  Regional lake  characteristics                               80
         5.5  Simulated lake water temperature  regimes for
             historical  and  future weather        . .  .                    83
             5.5.1  Water temperatures                                 83
             5.5.2  Thermal energy  flexes                               89
             5.5.3  Vertical mixing/stratification/stability                 95
         5.6  Conclusions                                             102

     6.   SUMMARY                                                 103

     7.   REFERENCES                                              105

APPENDIX A.      Vertical diffusion in  small stratified lake:
                   Data and error analysis                            A.I
                   A.I   Introduction
                   A.2   Study site
                   A.3   Vertical eddy  diffusivity
                   A.4   Sediment  heat storage
                   A.5   Water temperature observation
                   A.6   Vertical eddy  diffusivity estimates            A-17
                   A.7   Error  analysis
                   A.8   Conclusions

APPENDIX B.      Temperature equation discretization                 A-24

APPENDDC C.      Temperature equation linearization                 A-25

APPENDIX D.      Cross-term evaluations                            A.28

APPENDIX E.      Regional lake water  temperature simulation model   A-30
                   E.I   Lake input data file
                   E.2   Example  input data  file
                   E.3   Meteorological data file
                   E.4   Example  meteorological data file              E.4
                   E.5   Program  listing


                                   iv

-------
                               List  of Figures
Figure 1.1     Schematic diagram of  source and sink  terms in the heat budget
              model.
Figure 2.1     Hypolimnetic eddy diffusivity dependence on lake  surface area.
Figure 2.2     Hypolimnetic eddy diffusivity  forcing parameter  (a)  dependence
              on lake surface area.
Figure 2.3     Maximum hypolimnetic  eddy  diffusivity (at  N2=7.5*10'5 sec"2)
              dependence on lake  surface area.
Figure 2.4     Relationship  between  total attenuation coefficient  and  Secchi
              disk depth.
Figure 2.5     Wind sheltering coefficient  dependence  on lake surface area.
Figure 2.6     Wind function coefficient dependence on lake  surface  area.
Figure 2.7     Lake "wind speed measurements.
Figure 2.8     Ecoregions and spatial distribution of selected lakes.
Figure 2.9     Cumulative    distributions   (%)   of   lake   parameters   in
              Minnesota.Lakes  selected  for  model  validation  are  shown  by
              symbols.
Figure 2.10   Lake  Calhoun water temperature profiles.
Figure 2.11   Square Lake  water temperature profiles.
Figure 2.12   Waconia Lake water temperature profiles.
Figure 2.13   Thrush  Lake water  temperature  profiles.
Figure 2.14   Williams Lake water temperature profiles.
Figure 2.15   Standard  deviations  of,  estimated   lake   water  temperature
              uncertainties.
Figure 2.16   Standard   deviations  of   estimated   lake   water  temperature
              uncertainties.
Figure 2.17   Standard   deviations  of   estimated   lake   water  temperature
              uncertainties.

-------
Figure 2.18   Simulated  temperature (isotherm)  structure  in  Thrush  Lake.
              Top  shows  results  from  validated  model  and  bottom shows
              results from regional model.

Figure 3.1     Schematic  illustration  of  the  lake  temperature  perturbation
              system.

Figure 3.2     Meteorological variables at  Minneapolis/St. Paul.   Daily means
              and standard deviations for the  period  1955-1979.

Figure 3.3     Estimated   long-term  average   epilimnion   and  hypolimnion
              temperatures.

Figure 3.4     Long-term  average isotherms in Lake Calhoun.

Figure 3.5     Standard   deviations  of   estimated   epilimnion  temperature
              uncertainties.  Contributions by several meteorological  variables
              and totals are shown.
                  x
Figure 3.6     Standard  deviations   of   estimated deep  water  temperature
              uncertainties.  Contributions by several meteorological  variables
              and totals are shown.

Figure 3.7     Long-term   average  temperature  profiles  plus  or  minus  one
              standard deviation in Lake  Calhoun.

Figure 3.8     Epilimnion   temperature  long-term  average plus  or  minus  one
              standard deviation.

Figure 4.1     Lake Elmo water temperature profiles.

Figure 4.2     Lake Holland water temperature profiles.

Figure 4.3     Lake Calhoun water temperature profiles  in 1971.

Figure 4.4     Cumulative  evaporative losses (simulated).

Figure 4.5     Mean monthly equilibrium temperatures (simulated).

Figure 4.6     Mixed layer depths (simulated).

Figure 4.7     Simulated epilimnion  temperatures.

Figure 4.8     Simulated hypolimnion temperatures.

Figure 5.1     Regional boundaries  and   geographic  distribution  of  lakes in
              MLFD database.

Figure 5.2     Geographical  location  of  the   closest  GISS  grid  points   for
              Minneapolis/St. Paul and Duluth.
                                     VI

-------
Figure  5.3     Climate parameters  a Minneapolis/St.  Paul  in the  past  and
              under a 2xCC>2 (GISS) climate  scenario.

Figure  5.4     Geographic distribution of lakes  according to key  parameters:
              Secchi depth, maximum depth,  and surface area.

Figure  5.5     Cumulative distributions  (%)  of  key parameters  in Minnesota
              lakes (from MLFD database).

Figure  5.6     Horizontal area vs. depth relationship for lakes.  Area
              and depth are normalized.

Figure  5.7     Simulated  weekly epilimnion temperatures.

Figure  5.8     Simulated  weekly hypolimnion temperatures.

Figure  5.9     Examples   giving   range  of  epilimnetic  and   hypolimnetic
              temperatures  over a  25 year period (95% confidence  interval).

Figure  5.10    Simulated  temperature (isotherm) structure in (a)  three medium
              deep (13 m maximum depth) lakes of  large (10 km2), medium
              (1.7  km2)  and  small  (0.2 km2) surface  area,  (b) three medium
              size lakes  (1.7 km2  surface area)  with maximum depths of 4 m,
              13 m,  and 24 m.  Isotherm  bands  are in increments of 2C.
              Simulated  water  temperatures  are past 1955-79  climate  (top)
              and 2XC02 (GISS)  climate scenario (bottom).

Figure  5.11    Examples  of individual surface heat flux components.

Figure  5.12    Simulated  cumulative net heat flux.

Figure  5.13    Simulated  cumulative evaporative losses.

Figure  5.14    Simulated  weekly mixed layer depth.

Figure  5.15    Simulated  lake numbers as a function of lake depth  and trophic
              status.

Figure  A.I    Ryan Lake bathymetry.

Figure  A.2    Temperatures  recorded in  lake sediments  and overlying  water,
              1990.

Figure  A.3    Calculated and measured temperatures in lake sediments,  1990.

Figure  A.4    Hypolimnetic  lake   water  temperatures   recorded  at  2  min
              interval in Ryan Lake, May 7 to August 9, 1989.

Figure  A.5    Meteorological  conditions  (solar  radiation, wind speed and  air
              temperatures)  during  the period of investigation.
                                    Vll

-------
Figure A.6    Seasonal  lake temperature structure. Isotherms  (C)  shown  in
              this figure are derived from  measurements  at  20 minute and  1
              m  depth intervals.

Figure A.7    Heat  flux  through  the  sediment-water   interface   and  solar
              shortwave  radiation  received  at  the  4m  depth,   May  7  to
              August 19, 1989.

Figure A.8    Vertical turbulent diffusion coefficient time  series in  Ryan  Lake.

Figure A.9    Calculated vertical eddy  diffusion  coefficients for time  intervals
              of  1, 5, 10, and  15 days,  with  and without sediment  heat  flux.

Figure A.10   Mean values  plus  or  minus two standard deviations  of eddy
              diffusion  coefficient as a function of  depth.

Figure A.11   Estimated  eddy  diffusion errors  (2ffv  )  for  different  sampling
                                                  J^-Z
              intervals.

Figure A.12   Estimated   eddy  diffusion  errors  20w  (cmV1)    space-time

              tradeoff,  7 m depth July  19.
                                     Vlll

-------
                                List of Tables
Table 1.1     List  of  calibration  coefficients  with  ranges used  in  previous
              simulations.

Table 2.1     Quantitative  measure  of  the success  of  simulations-Validated
              model.

Table 2.2     Coefficients for calibration of water temperature model.

Table 2.3     Coefficients for uncertainty analysis.

Table 2.4     Quantitative measure  of the  success of the simulations-Regional
              model.

Table 3.1     Correlation  coefficients of  daily  meteorological  variables  for
              Minneapolis/St.  Paul,  1955-1979.
                              -\

Table 4.1     Lake data.

Table 4.2     Mean monthly meteorological data.

Table 4.3     Differences   (C)  in   simulated   mean  daily  epilimnetic  and
              hypolimnetic  temperatures for  different  starting  dates of  the
              model  (April 1 reference).

Table 4.4     Monthly averages of  daily heat balance components  (1000  kcal
              m-2 dayi).

Table 4.5     Cumulative heat balance components  (1000 kcal nr2).

Table 4.6     Net  cumulative  heat input  (content)  per  meter  of  average
              depth  (1000 kcal m-i).

Table 5.1     Weather parameters  changes  projected  by  the  2xCC>2  climate
              model  output for Minneapolis/St.  Paul.

Table 5.2     Lake classification.

Table 5.3     Morphometric  regression  coefficients  in  the  area  vs.  depth
              relationship.

Table 5.4     Maximum  temperatures of southern Minnesota lakes.

Table 5.5     Seasonal   stratification  characteristics  of  southern   Minnesota
              lakes.
                                      IX

-------
Table A.I     Regression coefficients for K2 = o(N2)7.
Table A.2     Errors of the eddy diffusivity estimation.

-------
                         Acknowledgements


     We are grateful  to the U.S. Environmental Protection  Agency, Office of
Program  Planning and  Evaluation,  Washington, D.  C.  and  Environmental
Research  LaboratoryDuluth, Minnesota; and the International Student Work
Opportunity Program  (ISWOP), University of Minnesota  for  support  of  this
study.

     Staff  members   of  the  Environmental  Research Laboratory,  Duluth,
provided  information,  constructive  comments, and  suggestions for the study.
Special thanks goes to John Eaton,  Kenneth Hokanson,  Howard McCormick,
and Brian Goodno.   Finally we  wish  to thank  those who provided field data
for  the  study,  especially  Dennis  Schupp  and David   Wright  (Minnesota
Department  of Natural  Resources), Richard  Osgood  (Metropolitan  Council),
Donald   Baker,   David   Rushee  (Soil  Science  Department   University  of
Minnesota),  Roy  Janne,  Dennis Joseph (Center for  Atmospheric  Research,
Boulder), and  Tom Winter  (U.S. Geological Survey).

     Last but not least,  we extend our gratitude to members of the SAFHL
research staff for  assistance  and support.

-------
              1. Introduction and Literature Review
1.1  Introduction

     The concentrations  of some gases (CC>2, H20, N20,  CE^)  have  been
increasing  in the  atmosphere  (Bolin  and  Doos,  1986;   NRG,  1982;  1983;
Houghton et al.,  1990).    These  commonly  called  "greenhouse  gases"  are
absorbing and reradiating  energy  at both long  and short  wavelengths.  As  a
consequence,  greenhouse  gases  are  able  to  affect   global  climate  possibly
resulting in  global mean warming of the earth's terrestrial  and aquatic surface
and  the lower  atmosphere  (Bolin and  Doos, 1986; NRG,  1982;  1983;  Wanner
and  Siegenthaler, 1988; Waggoner, 1990).

     Special  attention  has been  paid  to  the  increase  of carbon  dioxide
because it  is estimated that about half  of  the temperature change is  due to
the  increase of atmospheric  CC-2  alone.    Mathematical  models  of global
climate  change  lead  to  the  conclusion  that   the increase  in  mean global
equilibrium surface temperature for a doubling  of CO2 is most likely to be in
the range of 1.5 to 5.5 C  (Bolin and  Doos, 1986, Waggoner, 1990).  One of
the  uncertainties  is  due to  the  transfer  of increased  heat  into  the oceans
(NRG,  1982;  1983, Waggoner,  1990).   Surely due to  their  high  heat capacity,
oceans will act  as a sink for heat and  delay the warming.

     The question which we want to address in this  report is  how  freshwater
lake temperatures respond  to atmospheric conditions.   Changes  in  lake  water
temperatures and temperature  stratification dynamics may  have  a profound
effect on lake ecosystems  (Meisner et  al.,  1987; Coutant,  1990;  Magnuson et
al.  1990; Chang et  al.,  1992).   Dissolved  oxygen, nutrient cycling, biological
productivity,   and  fisheries  may   be  severely   affected  through  temperature
changes.

     Considerable effort has gone into global climatological modeling with the
objective to  specify future climatic conditions in a  world with high greenhouse
gases.,    Some  models use  statistical analysis  of past  climatological  data in
order to provide  analogies for  future  climatological  changes.    Unfortunately,
all  causes of past climate  changes are not  fully understood (Bolin  and Doos,
1986;  Waggoner,  1990),  and  predictions   of   future  climates   are  difficult,
especially on a  regional basis.   Nevertheless simulated  climate  conditions are
and  will be  used in  numerous effect 'studies.    Another approach  to  finding
both  climatic trends  and  their effects is to examine long-term records.   In
few  lakes,  e.g.  in  the  experimental lake  area (ELA) in Ontario,  Canada,
weekly  or biweekly  vertical profiles of water quality  and biological parameters
have been  collected  over periods of 20 or more years  and  these records reveal
e.g.  rising  average surface  water  temperatures,  shorter ice cover periods,  etc.
(Schindler  et al., 1990).    A data record  of more than  100 years  for  Lake
Mendota was analyzed by  Robertson (1989)  and Magnuson  et  al. (1990).

-------
      To  make  generalizations  to lakes of different geometries and latitudes,
and  to extrapolate to  possible  future  climates,  numerical simulation models
(McConnick,  1990; Robertson  and  Ragotzkie,  1990) are  useful.    Herein  the
use  of such  a  model  is  demonstrated by  application  to  morphometrically
different lakes with sparse  data sets.  The  lakes are located near  45  northern
latitude and 93 western longitude in the northcentral United States.


1.2   Previous Temperature Prediction  Model

      A one-dimensional lake water quality model,  which has been  successfully
applied to simulate hydrothermal processes  in different lakes and  for a variety
of meteorological conditions (Stefan and Ford, 1975; Stefan et al., 1980a; Ford
and   Stefan,  1980) was used   in this study.   The  model was previously
expanded  to   include  suspended  sediment  (Stefan  et  al.,   1982V  light
attenuation (Stefan et al.,  1983),  phytoplankton  growth and nutrient dynamics
(Riley and Stefan, 1987).     Only  the hydrothermal  part of the  model  was
applied in this  study.
        ,        /
1.2.1     Model formulation

      In  the  model the lake  is  described  by a system  of horizontal  layers,
each  of  which  is  well  mixed.    Vertical transport  of heat is described by  a
diffusion   equation in  which   the  vertical  diffusion   coefficient  Kz(z)   is
incorporated in  a conservation  equation of the form:


                      A fr = I 
-------
CO
               Ha
          atmospheric
            radiation
        Hsn
He
He
Hbr
incoming   reflected evaporation   convection  back
  solar    fraction                          radiation
radiation
                                                   t           r      t
                                            surface absorption
                                            absorption  with depth
                                        fraction  transmited
                                           to next layer
                                             absorption with depth
                                        fraction transmited
                                           to next layer
          Fig.  1.1    Schematic diagram of source and sink terms in the heat budget model.

-------
where Hs is incoming  solar  radiation  (kcal  m^day"1),  r is  the  reflection
coefficient computed   as  a  function  of  the  angle  of  incidence  and  the
concentration of suspended  sediment  in the surface layer  (Dhamotharan, 1979;
Stefan et al., 1982).  /? is  the surface absorption  factor  (Dake and Harleman,
1969).  The  attenuation of  solar  radiation with depth follows  Beer's law:

                         H(i) =  H(i-l) exp(-M  Az)                    (1.4)
                                   gn
where  Hsn(i-l) is  solar  radiation at  the  top of a  horizontal layer  of  water
        __ n   ^*
(kcal m  day  ), Hsn(i) is solar radiation  at the  bottom  of a layer,  A is

thickness  of  a  layer (m), fj, is the extinction coefficient (m~ )

                         fj. =  \i  +  u -SS + /i.Chla                    (1.5)
                         r    *w    rsa       ;ch                        v   '


where  /^ is  the extinction coefficient of lake water (m~ ),  ^ss is the specific
extinction  coefficient due  to suspended sediment  (1 m^mg"1); SS  is  suspended

inorganic sediment concentration (mg  I"1); \L& is  the  extinction coefficient  due
to   chlorophyll   (m2   g"1Chla)(Bannister>   1974),   Chla   is   chlorophyll-a
concentration (g m'3).

                                H  = a e  T4                            (1.6)
                                 a       a   a                            v   '
where  a is  Stefan-Boltzmann constant, Ta is absolute temperature (K),  ea is
atmospheric  emissivity (Idso and  Jackson,  1969).  Back radiation  Hbr follows
the  same  formulation  (6),  but  the  emissivity   is  fixed  at  0.975,  and
atmospheric  temperature is replaced by water surface temperature Ts  .

     Aerodynamic bulk formulae were used to  calculate surface  wind shear r,
latent  heat flux H , and the sensible heat flux H   across  the  water surface
                   c'                              c
(Keijman,  1974;  Ford and Stefan, 1980; Strub and Powell,  1987; Sadhyram et
al., 1988):
                                                     Uf                  (1-7)

           H  =  p c   8 ' u'  = p c C u, 6. =  p c f(U )(T  -  T )       (1.8)
             c    ra p            'a p s *  *    'a p v a'v s     a'       v   '
            H  = p L   q ' u'  =  p-L C^u. =  p L f(U )(q  - q )        (1.9)
             e    ra v  n         ra v C1* *     a v x  a/vT     a'        v   '
where T is the surface wind  stress,  pa is  the  density  of  the  air, u'  and  u;'
are. turbulent fluctuations  of  velocity in horizontal and vertical direction; the
overbar  represents a  time average;  u,   is a  velocity  scale, Ua is the  wind

speed above the water  surface,  Cd is the momentum or drag  coefficient  (Wu,
1969),
0'  is turbulent fluctuation in temperature,  6^  is a temperature scale, Cs and

Cf are  heat transfer  and vapor transfer  coefficients, respectively, and  together

-------
with u. are expressed as a  function of wind  speed, f(Ua), (Ford, 1976), Ts is

water surface  temperature,  Ta  air temperature  above the water  surface, Lv is
'.atent heat  of vaporization, q'  is the specific humidity fluctuation, q^ is  the

specific  humidity  scale,  qa is the specific humidity above the water surface, qs
is   the  specific  humidity  at  saturation  pressure   at  the  water   surface
xemperature.

     Turbulent  kinetic  energy supplied  by  wind  shear  and  available  for
possible entrainment at  the interface  was estimated  (Ford  and  Stefan, 1980)
bv
                          TKE = Wstr
dA
                                                                       (1.10)
where  As  is  lake surface  area  (m2), Ut  is  shear velocity  in the  water (m

day-i), and Wstr is the wind sheltering coefficient.

     The  model distributes the  surface heat  input  in  the water column using
turbulent diffusion (Eq.  1)  in response to  wind and  natural  convection (Ford
and  Stefan,  1980).   The numerical  model  is applied in  daily timesteps using
mean daily values for the meteorological  variables.   Initial conditions,  model
set-up parameters, and daily  meteorological variables  average air  temperature
(Ta), dew  point  temperature (Td),  precipitation  (P),  wind  speed  (Ua),  and
solar radiation (Hs) have to be  provided to use  the model.
1.2.2
          Model  coefficients
      Model  calibration  coefficients  needed  for  simulations  of  lake  water
temperatures  are given in  Table  1.1.   These  coefficients are  kept at  their
initially specified value throughout  the entire period of the simulation.

Table  1.1  List  of  calibration  coefficients  with  ranges  used  in  previous
            simulations.
Coefficient
                         Symbol     Units
       Range of values
                                                   Previous      Literature
                                                   Simulations   Values
Radiation extinction
by water

Radiation extinction
by chlorophyll

Wind sheltering

Wind function
coefficient
 Maximum hypolimnetic
 eddy diffusivity          Kzn,ax
                                  (mi)
                                  (m2 g-i Chla)
   0.4-0.65
0.02-2.0
                         W
                           str
   8.65-16.0      0.2-31.5

   0.1-0.9        0.1-1.0


   20.0-30.0      20.0-30.0
                                  (m2 day1)
   0.1-2.0
0.086-8.64

-------
      Radiation  extinction  coefficients  by  water  (^w)  and  chlorophyll
specify the rate of attenuation of short-wave radiation energy  as  it penetrates
through  the  water  column.    Both coefficients  vary  as  a  function of  the
wavelength.   Usually  these coefficients are reported by  a  single mean  spectral
value  for  a  given lake.   Smith   and  Baker (1981)  measured   a  range  of
0.02-2.0  (m*1) for  /iw as a function of the wavelength.   Values  of fa in  the
range  of 0.68   0.35  (m-i)  have  been  reported by Megard  et al.  (1979).
Chlorophyll extinction coefficient  is species dependent.  Values in  the range of
the  0.2-31.4  (m2  g'1 Chla)  with  a mean spectral  value  of  16.0 have been
reported  by  Bannister   (1979;  1974)  for  IL< while  Megard  et al.  (1979)
reported  values  of 22   5  (m2 g'1 Chla)  for the  photosynthetically active
radiation.

      The  wind   sheltering   coefficient  (Wstr)  determines  the   fraction   of
turbulent kinetic  energy from  the  wind  applied at  the lake  surface  and
available for mixing.  The coefficient can  range from 0.1  to 1.0 depending  on
the  size  of the  lake  and the terrain surrounding the  lake.   The  coefficient
defines the "active" portion   of  the lake  surface  area  on  which wind  shear
stress contributes to the turbulent kinetic energy.

      The wind function coefficient  is  defined for  the neutral boundary  layer
above  the lake  surface.   This  condition  occurs  for the  case  of  negligible
atmospheric stratification.  The  wind speed  function used  is  linear with  the
wind speed

                                f(Ua) = cUa                          (l.H)

where  c  is defined as  a wind function  coefficient.   The  atmosphere above
natural water bodies is  often  nearly neutrally  stable.  A significant amount of
experimental  and  theoretical  research  has  been  done  in  regard   to  wind
function  coefficient estimation (e.g.  Dake,  1972,  Ford,  1976,  Stefan et  al.,
1980b, Adams et  al.,  1990).   Different  ranges of coefficients were reported
depending  on measurement location of the windspeed Ua relative to the lake
surface.   Herein the  wind function coefficient is taken  to be  in  the  range
20-30  if wind speed is  in mi h'1, vapor  pressure in mbar, and  heat flux in
kcal m^day1.

      Maximum hypolimnetic   eddy   diffusivity  is the  threshold value  for  the
turbulent diffusion under negligible  stratification.   In modeling this condition
is  assumed to be  satisfied  by  small  stability frequency e.g. N2 =  7.5*10'5
sec"2.  Maximum hypolimnetic eddy diffusivity ranges from 8.64  m2 day'1  for
large lakes (Lewis, 1983) to 0.086 m2 day1 for small  lakes  (Appendix  A).


1.3   Overview of Study

      The goal of this study is to develop an  understanding of how  freshwater
lake temperatures respond to  atmospheric conditions.

      Chapter  2   presents   the  regional   lake   water  temperature  model
development  and  validation.   The lake water temperature model, which  was
originally  developed  for particular  lakes  and   particular  years  has   been
generalized to a  wide range of lake  classes and meteorological conditions.

-------
     Chapter  3 presents  a first order analysis  of  uncertainty propagation in
lake  temperature  models.     The  source  of  the  uncertainty  is  variable
meteorological  forcing which  enters  the lake temperature  equations through
the source term and boundary conditions.   The  analysis  presents  a  useful
alternative for  the  study  of long-term  averages  and variability of temperature
structures in  lakes  due to variable meteorological forcing.

     Chapter  4  presents a  lake water  temperature model application  to a
particularly warm year (1988)  and a normal year (1971) for comparison.   The
comparison is made for  morphometrically different  lakes  located in  the north
central US.  The  analysis  was  a first step  in  quantifying  potential thermal
changes in inland lakes due to climate change.

     Chapter  5  presents a  lake water  temperature model application  to a
representative  range  of  lakes  in Minnesota  for past  climate  and  a  future
climate scenario associated with doubling of atmospheric  CO2-   Emphasis  was
on long  term behavior and a  wide  range  of lake  morphometries and trophic
levels.    The   base  weather  period  (or reference)  was  for  the years  from
1955-1979.   For future  climate scenario  the daily weather parameters  were
perturbed by the 2XC02 GISS climate model output.   The simulation results
showed how  water temperatures in  different freshwater  lakes  responded to
changed  atmospheric  conditions in a region.

     Chapter 6 summarizes the results of the study.

-------
       2.  Regional Lake Water Temperature Simulation

                         Model Development


     A lake water temperature model  was developed to simulate the seasonal
(spring to  fall) temperature  stratification  structure  over a wide  range  of lake
morphometries,  trophic  and meteorological  conditions.    Model  coefficients
related  to  hypolimnetic eddy diffusivity, light attenuation,  wind  sheltering,
and  convective heat transfer, were generalized  using theoretical  and  empirical
model  extensions.   The  new relationships  differentiate  lakes  on a regional
rather  than  individual  basis.     First   order  uncertainty   analysis  showed
moderate   sensitivity  of  simulated   lake   water   temperatures  to   model
coefficients.   The proposed regional  numerical model  which  can  be  used
without calibration has an average 1.1* C  root mean square error, and  93%  of
measured  lake water  temperatures variability  is explained by the  numerical
simulations,  over wide  ranges   of  lake  morphometries,  trophic levels,  and
meteorological conditions.


2.1  Introduction

     Changes  in  meteorological  variables   in   the  future  "greenhouse"
atmospheric conditions are usually specified through the global climate change
models  output  on a regional   rather  than  a local  scale.    Usually water
temperature  data are  only   available  for a  few  lakes,  not necessarily  for
"typical" lakes in order  to   calibrate  lake water temperature  model  and  to
validate predictions.   Some  coefficients such  as eddy  diffusion coefficients  or
turbulence  closure coefficients used in  lake water temperature models  are not
universal  due  to their  dependence  on   stratification,  and  the  longer  than
subdaily time step of the simulations (Aldama  et al.,, 1989).

     The  purpose of this   chapter  is  to describe how  a  lake temperature
model, which was described  in  Chapter  1, and which  was initially  developed
for particular lakes  and particular meteorological years,  could be generalized
to a wide  range  of lake  classes  and meteorological conditions.   To do  this,
new   functional  relationships had  to  be  introduced  for   the  calibration
coefficients which differentiate lakes  on  a regional rather than  an  individual
basis.   The generalized model  can  than be applied  to  lakes  for  which  no
measurements exist.   Fortunately -it  can be  demonstrated that the  regional
model  makes  prediction  almost  with  the  same  order  of  accuracy  as  the
validated previous calibrated to  particular lakes.  Therefore regional and long
term lake  temperature structure modeling rather than  short  time  behavior of
particular lakes can be accomplished  with same confidence.

-------
2.2   Model Generalization

      In  order  to apply the lake water  temperature numerical model to lakes
for which there  are no  measurements, the model has to be  generalized.   This
was   accomplished  by  introducing   functional  relationships  for  the  model
coefficients  which  are valid for lakes  on a  regional  rather  than  individual
basis.

2.2.1  Hypolimnetic diffusivity closure                                .    .

      Although the hypolimnion  is isolated from the surface  (epilimnetic) layer
by  the thermocline and  its associated density  gradient, strong and  sporadical
local  mixing  events  have  been observed in  the  hypolimnion  (Jassby  and
Powell,  1975;  Imberger,  1985;  Imberger  and  Patterson,  1989).    Heat flux
between  water   and  lake  sediments  was  found  to  be  important  in  eddy
diffusivity  estimation  for  inland  shallow  (10  m maximum  depth)   lakes,
representative for  north  central  United  States (Appendix A).   Hypolimnetic
eddy  diffusivity  dependence on  stratification strength  measured by buoyancy
frequency has been pointed  out  consistently (Quay et al.,  1980; Gargett, 1984;
Gargett  and  Holloway  1984;  Colman  and  Amstrong,  1987;  Appendix  A).
Stability frequency is  related to  hypolimnetic eddy diffusivity by :


                                K2 = a (N2)7                           (2.1)

where stability  frequency N2=-(dp/cte)(g//?),  in which  p is  density of water,
and  g is acceleration  of gravity, 7  is determined by  the mode of  turbulence
production (narrow or broad band  internal waves,  local shear etc.), and a is
determined  by   the  general level  of  turbulence.    For  most  inland  lakes,
coefficient  7 ranges from 0.4  to 0.6 (Jassby and Powell, 1975; Quay  et  al.,
1980; Gerhard  et al.,  1990;  Ellis  and Stefan, 1991, Appendix A).

      Hypolimnetic  eddy  diffusivity  estimations  in five  northern  Minnesota
lakes  follow  Eq.  2.1  as shown  in  Fig.  2.1.   Lakes  were  selected from  the
regional  prospective  i.e. with  different  surface areas  and  maximum  depths.
Dimensionless  analysis  (Ward,   1977)  suggests that  lake  surface  area,  can
provide  the horizontal scale for the  vertical eddy diffusivity estimation.   The
vertical  scale (lake depth) is implicitly built into the stability frequency.  The
a  coefficient in  Eq.  2.1  can be  interpreted  as a  measure of turbulence level
and  is  plotted  as a  function  of lake  surface area in Fig. 2.2.    A general
relationship applicable to lakes  on  a regional  scale  was therefore  summarized
as:
                    Kz = 8.17  x 10-4,(Area)0'56  (N2)-'                (2.2)

where Area is  lake surface area  (km2), N2 is  in sec'2,  and Kz is in cm2 sec'1.

      Maximum   vertical   hypolimnetic   eddy   diffusivity   Kzmax   was  also
correlated  with  lake  surface area because turbulent mixing in non-stratified
lakes  depends strongly on kinetic wind energy supplied, which in turn depends
on  lake  surface  area.   Maximum  hypolimnetic  eddy  diffusivity  versus lake
surface area for  eight different  lakes is  plotted in  Fig.  2.3.   Data are from
Jassby and Powell (1975),  Ward (1977), Lewis  (1983), and from this study.

-------
     10-
               Ann Lake  (0.47  km2)
 o
 (U
 E
 o
 o
 d>
 c/i
(N

 E
 o
     10-
     10-J^
10-M

   10
                 % :  
                          *
     .

Big Carnelion Lake ( 1 .89 km")
                                       Square  Lake (0.85 km2)
                                                       Big  Marine  Lake (74  krna)
          -
                                              10-2
                                                        10
                                                                          -3
10'2
                        N2(sec-2)                              N2(sec-2)



             Fig.   2.1   Hypolimnetic eddy diffusivity dependence  on lake surface  area.

-------
             = 8.17X10-4(Area)-56
           Ryan  Lake  (0.06 km2)
           Ann Lake (0.47 km2)
           Square Lake  (0.85 km2)
           Big Cornelian Lake ( 1 .89 km")
           Big Marine Lake (7.4 km)
           Lake Mcllwanie  (26.3 km)
10-*-
                             Area (km2)
Fig.   2.2   Hypolimnetic eddy diffusivity forcing parameter (a)  dependence on
           lake surface area.

-------
 U
 (D
 (f)
CM
 U
  N
10-=
       IU '
         -2_
       10
         -J_
              o
              D
i i i iii|	1ii i 11ii|	1ii i 111ii	r
Ryan  Lake  (0.06 km2)
Castle Lake (0.2 km2)
Ann Lake (0.47  km2)
Square Lake  (0.85 km2)
Big Cornelian Lake (1.89km2)
Big Marine  Lake (7.4  km2)
Lake  Mcilwanie  (26.3  km2)
Lake  Valencia (350 km2)
Ryan  Lake (winter)
K7mox = 0.048(Area)-56
  10-2
                       10
                                 -'
                                    Area
                                                    101
102
            Fig.  2.3   Maximum  hy poll nine lie eddy diffusivity (N2=7.5*10'5 sec"2)
                       dependence on lake surface area.

-------
          Attenuation coefficient

     The specific radiation  attenuation  coefficients  for  water and chlorophyll
 -TO replaced by the total  attenuation  coefficient.   This  was done following
. '..( parsimonious principle i.e. the fewer  coefficients  in the  model, the less
.-.ncLTtain the  model estimate.  In addition uncertainty  analysis showed that
 hlorophyll-a   made  a   minor   contribution   to  lake   water   temperature
uncertainty.   A  relationship  between total attenuation coefficient \i (m'1) and
y-cchi   depth   zsd  (m)  was  obtained   from  measurements in  50  lakes  in
Minnesota (Osgood,  1990)  and is plotted in Fig.  2.4.

                               It =  1.84 (Zsd)-i                          (2.3)

I'he form of this relationship has been found  to  be valid  in inland  waters  in
St-neral  (Idso and Gilbert,  1974) and in the ocean (Poole and Atkins,  1929).

L' 2.3      Wind sheltering  coefficient

     The wind  sheltering  coefficient is a function of lake surface area  (fetch).
The turbulent  kinetic energy computation (Eq.  1.10) uses  a  wind  speed and
direction  taken, from  off-site  weather  station  at  10 m  elevation and  adjusts
that wind speed for  fetch over  the lake  in the direction  of  the wind.   As
wind  speed typically increases  with fetch,  the   calculated downstream  wind
speed  is  an  estimate of  the  maximum  wind   speed  on  the  lake  surface.
Typically fetch  on a lake is  reduced  by wind sheltering the upwind  side  of
the lake  where  the  wind makes  a transition  from  a landbound  turbulent
velocity  profile  to the open  water.   This  was explained by Ford and Stefan
(1980).    The reduction in fetch  or surface area sheltered from direct  wind
access  by trees or buildings  along the  shoreline  will be more  significant  for
small lake than a large one  because a) a relatively  larger portion of the total
lake surface  area will be wind  sheltered  b)  the  downwind  maximum  wind
speed does not  grow linearly  with  fetch and will on a large lake be near  the
real wind speed over a large  portion of the lake  surface area, and c) wind
gusts will be less effective over  a small  lake surface than a large one  because
of  spatial averaging.   Also  lake  morphometry,  i.e. distribution  of area with
depth  will be  a factor in the  translation  of wind energy into mixing.   A
maximum wind speed at the  downwind  end  of a large lake will  also  be more
representative  for  a  large  lake  than   a  small   one,  especially  if  the  lake
morphometry is taken into consideration.

     For all  these  reasons  a very  strong dependence of the wind  sheltering
coefficient  (Wstr) on lake  surface  area can  be expected.     A  functional
relationship was obtained  by plotting  the  wind  sheltering  coefficient  obtained
by  calibration  in  several  previous numerical model  simulations  (Fig.  2.5).
Biweekly temperature profile measurements  in ten  lakes and  throughout  the
summer  season  were used to optimize the wind  sheltering  coefficients plotted
in  Fig.  2.5. ' The empirical relationship is

                        Wstr =  1-0 - exp(-0.3*Area)                    (2.4)

where  Area is  the lake surface area in km2.
                                     13

-------
              84*(secchi  depth)-1
          .,,,,,,.,
        0.3     0.5    0.8     1.0     1.2    1.5
                    (Secchi  Depth  m  )"'


Fig.  2.4   Relationsliip  between total attenuation coefficient and Sechi disk
           depth.

-------
 c
 OL)
 O
 0>
 O
O
 O>
 C
  102-
  10'-
  10-:
 CD
L   10--:
 0>
00

?   10-H
 io-3-
    10-'
            Wstr= 1,0-EXP(-0.3*Area)
10-2
10-
                                          rnrj

                                            10
                       Surface  Area (km2)
10'
                                                                    E  Williams (0.35)
                                                                    X  Riley (1.2)
                                                                    ^  Waconia (10.0)
                                                                    O  Square (0.8)
                                                                    A  Thrush (0.07)
                                                                    A  Greenwood  (7.70)
                                                                    ffl  Fish (1.16)
                                                                    D  Elmo (1.23)
                                                                      Cedar (3.30)
                                                                    O  Calhoun (1.71)
102
Fig.   2.5   Wind sheltering coefficient dependence on lake surface area.

-------
      The  result  in  Fig.  2.5 seems to  indicate  that  the modeling of  wind
mixing in  lakes, especially small ones, depends  more on a correct  amount of
energy  supplied  than on  a energy  dissipation.  This  is  a new  insight  which
appears to result from this study.

2.2.4     Wind function coefficient

      Wind function  coefficients (c) defined in  Eq.  1.11  enters  into  the heat
transfer  relationships  (1.8) and (1.9),  and depends  also  on lake surface area
(fetch)  as  was  found by  Harbeck  (1962),  Sweers  (1976),  and summarized  by
Adams et  al. (1990).   Harbeck (1962) analyzed data from several lakes  of
different sizes and pointed out that evaporation rates in  small and  large  lakes
-might be the same.   The fetch dependence  is  introduced  mainly due  to the
wind speed increase  over  the water.   As  air flows from land  to a  smoother
water surface,  at a constant height above  the water (e.g.  10 m), its velocity
increases with  fetch.   In  this numerical  model off-lake  wind speeds  measured
at permanent weather stations  are are used, but  they  are adjusted for lake
fetch  (Ford  and Stefan,  1980).   Nevertheless, some residual   wind function
coefficient  dependence on  lake  fetch  is  shown  in Table 2.1.    A  functional
relationship  was obtained  by  plotting  the  wind  function  coefficient  from
several previous numerical model  simulations against  lake surface  area  (Fig.
2.6).   The estimated relationship is

                              c =  24+ln(Area)                          (2.5)

where Area is  again  in km2.  This relationship  shows  only a week  dependence
of c  on lake surface  area, and can be viewed  as  a minor  adjustment.  The
need for this adjustment  can be  explained by  examining the wind  boundary
layer development over the  surface of small and large  lakes  (see  Fig.  2.7).
Wind speed  increases  with  fetch (distance from the  leeward shore)  but
non-linearly.   In our model  wind speed  is  taken  from  an  off-lake weather
station and  a  maximum   wind speed  at  the  downwind  end of the lake  as
shown  in  Fig.  2.7   (top)  is  computed  for  the  use  in the  heat  transfer
equations (1.8)  and  (1.9).   This calculated wind  speed  is an overestimate of
the  area!  average  wind   speed  over  the lake  surface.    Because  of the
non-linearity of wind speed with  distance  the overestimate is more severe  for
small lakes than for  large lakes.   Therefore  the wind function  coefficient has
to be smaller for smaller  lakes in  order to compensate for the wind velocity
overestimate.    If on  the  other hand,  wind speeds are  measured on the lake
(middle  of the  lake)  as shown  in  Fig. 2.7 (bottom) the situation is reversed.
In   that   case   the   wind  measurements  on  a  small  lake   are  severely
underestimated  relative to the  area!  average than  on a big lake.   For this
reason the wind function  coefficient has  to decrease with fetch  (surface  area)
to compensate for this  non-representativeness of the wind speeds measured in
midlake.   This  decreasing  trend of wind function coefficient with lake surface
area  was   found  and reported   by  Harbeck   (1962),   Sweers (1976)  and
summarized  by  Adams  et   al. (1990).    In  addition   Adams called  upon
increases  in  relative  humidity  with fetch over  cooling  ponds   to  justify the
decrease in wind function  coefficient with fetch.

      It  is concluded  from all  of the  above that the wind function coefficient
can  increase  or  decrease  with lake  surface  area  depending on the location
where wind speed is measured.
                                      16

-------
    40-
    35-
 G>
'u
 o  25^
    20-
    15-
      10-2
Wcoef = 24  + In(Area)
                                  .O   O  
                             X
       IO-1           10            10'
           Surface Area  (km2)
                                                                \ ~
                                                       X  Williams  (0.35)
                                                       X  Riley (1.20)
                                                       4  Woconia  (10.0)
                                                       O  Square (0.8)
                                                       A  Thrush (0.07)
                                                       A  Greenwood (7.70)
                                                       ffl  Fish (1.16)
                                                       D  Elmo (1.23)
                                                         Cedar (3.3)
                                                       O  Calhoun  (1.71)
      Fig.   2.6   Wind function coefficient  dependence on  lake surface  area.

-------
             offlake wind measurement
 wind
            on lake wind measurement
wind
       Fig.   2.7   Lake  wind speed  measurements.
                        18

-------
2.3   Water Temperature Model Validation After Generalization of
     Hypolimnetic Eddy Diffusivity

     The  model  was   first   modified   by  adding  the  hypolimnetic  eddy
diffusivity  closure (Eq.  2.2).   The  number  of calibration coefficients  (Table
1.1)  was thereby reduced  from five  to  four.   The  modified  numerical model
than had  to  be validated with  water  temperature  measurements in  several
selected  lakes  over  a  period of  several  years.    Representative  lakes  in
Minnesota  were  selected  through an analysis  of  the state's  extensive  data
bases.   Differences between waterbodies  in adjacent  ecoregions were found  too
small to justify  further subdivisions on this basis.   The state was  divided into
a northern part,  roughly coinciding with three  ecoregions,  and   a  southern
part,  roughly coinciding with three other ecoregions (Fig.  2.8)   which  also
extended into Wisconsin,  Iowa,   and South  Dakota.   "Representative"  lake
meant  either having values of lake  surface area, maximum  depth, and secchi
depth  near the  median  as  identified in  a  state  report  by  the Minnesota
Pollution Agency  (Heiskary  et al,  1988) or  being  near the far  ends  of  the
respective frequency  distributions  for ecoregions.   Selected representative lakes
with  their  position  on  the  cumulative  frequency  distribution   curves  for
northern and  southern  Minnesota  are given  in Fig. 2.9.   Lakes  covered  the
entire   range  of maximum  depths  (shallow-medium-deep),  surface  area
(smaU-medium-large),  and trophic status (eutrophic-mesotrophic-oligotrophic).
Geographical distribution of these lakes in Minnesota is given in Fig.  2.8.

     To  validate  'the   model   numerical   simulations  were   started   with
isothermal  conditions (4  *C)  on  March  1 and continued  in daily timesteps
until November  30.   Ice  goes out  of Minnesota lakes sometime between  the
end  of March and beginning  of  May.   Dates  of  spring  overturn vary  with
latitude and year.   To  allow  for  these variable conditions, a 4"C isothermal
condition  was maintained in  the lake  water  temperature  simulations until
simulated  water  temperatures began  to rise above  4*C.    This  method"
permitted  the model  to find  its  own date of  spring  overturn  (4C) and  the
simulated summer heating cycle started  from  that  date.

     Daily  meteorological  data   files  were  assembled  from Minneapolis/St.
Paul, and  Duluth, for  southern and northern  Minnesota respectively.

     A quantitative measure  of  the success of the simulations for  the  nine
representative lakes is given in Table 2.1.  Different gauges  of the simulation
success are  defined as: (a) volume weighted temperature averages
             p                     P
                 V > Tsi               V i
                                                                        (2-6)
                                     19

-------
          NorVhern  Minnesota Wetlands
              ^""s^^^^   i




          N 0 RTH   "\_XNa- GREENWOOD LAKE
                                   *&

                                  A

                                      THRUSH LAKE
                 rrgSv,

                  
-------
         !00
           10-'

         100
       10-'       10       10'
             AREA (km2)
         TOO-
                       MAXIMUM DEPTH (m)
          80-
   60-
i   40-
          20-
           0-r
                                      O
                                      o
                                                     D
                                                     O
                                                     A
                                                     D
                                                     O
                                                     A
                                                     A
                                                     A
                                                     
                                                     m
                                                     D
                                                     O
              468
          SECCHI DEPTH (m)
                                                 10
                                            Waconia (10.0)
                                            Greenwood  (7.2)
                                            Cedar (3.3)
                                            Williams (0.35)
                                            Calhoun (1.7)
                                            Elmo (1.2)
                                            Fish (1.1)
                                            Square (0.8)
                                            Thrush (0.07)
                                            Elmo (42.0)
                                            Greenwood  (34.0)
                                            Calhoun (24.0)
                                            Square (20.0)
                                            Waconia (11.0)
                                            Williams (10.0)
                                            Thrush (14.0).
                                            Fish (8.0)
                                            Cedar (4.7)
                                            Square (7.0)
                                            Greenwood  (6.5)
                                            Elmo (3.0)
                                            Calhoun (2.5)
                                            Williams (2.1)
                                            Thrush (7.3)
                                            Waconia (1.9)
                                            Cedar (1.2)
                                            Fish (1.0)
Fig.   2.9
Cumulative  distributions  (%)  of lake  parameters  in  Minnesota.
Lakes selected for model validation are shown by  symbols.
                                    21

-------
(b) temperature root mean square errors

       P                                 P

                    0.5                r  2rfVil

                                  Ejf=
                                                   I*
                                                              0.5
                                                                       (2.7)
 and (c) r2 i.e.  portion  of  the temperature  measurements  explained  by  the
 simulations (Riley  and Stefan,  1987).  In the above  equations subscripts  i, s,
 and m  refer to  the counting  index, simulated, and measured temperature
 respectively.  Vj is the  water  volume of a layer in  the  stratified lake.   The
 above   parameters  are  estimated  by  summing  over lake  depths.    Overall
 seasonal  average  parameters  are   reported  in  Table  2.1.    Examples  of
 simulated  and measured  vertical lake water  temperature  profiles  are given in
 Figs.  2.10,  2.11,  2.12,  2.13,   and  2.14.    The model  simulates  onset  of
 stratification,  mixed layer depth and water temperatures well.
      Table 2.1    Quantitative measure of the success of the
                   simulations-Validated model.
Lake
Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Year
1971
1984
1988
1987
1985
1985
1986
1986
1984
Tm
(oC)
14.37
20.64
13.94
24.40
14.37
20.14
11.80
11.97
17.26
Ts
(oC)
14.52
20.86
14.09
24.13
14.52
20.12
11.97
11.91
16.37
(oC)
0.86
0.94
1.77
0.80
0.86
0.78
0.89
0.90
1.08
E2
(oC)
0.79
0.99
1.80
0.82
0.79
0.73
0.79
0.91
1.07
.*
0.97
0.93
0.92
0.90
0.97
0.92
0.93
0.96
0.97
Number of
field data.
136
20
214
32
136
43
46
114
110
 Average

'T--
 Ts-
 E,-
 E2-
             16.54    16.49    0.97
0.96
0.94
95
     measured volume weighted average temperature
     simulated volume weighted average temperature
     temperature root mean square  error
     volume weighted temperature root mean square error
     portion of the measured water temperature variability  explained by
     the simulations
                                     22

-------
                            simulated
                            measured

                              s/20/7i
                                                                 e/oa/7i
                              6/39/71
                                                                 7/15/71
&
o  .J
    1
e/02/71
                                   S/2V1
I
LJ
o
  16-
                              tO/ll/71
                                                                 TO/25/71
           10      IS      20      25
               TEMPERATURE (*C)
                10      IS      20      25
                    TEMPERATURE (-C)
     Fig.    2.10    Lake  Calhoun water temperature profiles.
                              23

-------
I  10

a
   15-
                           measured

                          simulated

                                V15/85
                                                                   5/O9/SS
   10-
a
Ul
o
   15-
  20-
                                6/11/85
                                                                   7/12/85
  .5-

a.
UJ
o


   15-
  20-
                                S/07/SS
                                                                   9/CS/85
          S     10     IS    20    25


               TEMPERATURE (-C)
                                             5     10     IS    20    25    20


                                                   TEMPERATURE (*C)
      Fig.    2.11    Square  Lake  water temperature profiles.
                              24

-------
to
Cn
                              X

                              CL
                              LJ
                              O
                              I

                              CL
                              UJ
                              Q
u-
1-
2-
3-
4-
5-
6-
7-
8-
9-
10-
n_
1-
2-
3-
4-
5-
6-
7-

fl-
9-
10-
11-
 measured
i i _i
simulotea

5/17/85
'
1
i


T
.

 - -
-
.
-
-


-







8/20/85




1 * * 1 " *


 -


,






\ *

7/22/85 







l
l
l
i
i
(
/

















9/18/85
/*
f .
1

                                               10     15      20     25
                                               TEMPERATURE (C)
10     15     20     25   ,  30

TFMPERAfURE (C)
                                                Fig.   2.12   Waconia Lake  water temperature profiles.

-------
1  6
Q.  8
UJ
a
  10

  12

  1*
                           measured
                       	  simulated
                                                                 6/17/86
/
  10-

  12
                               7/15/86
O-  6-
I
CL
t. i
o
  10-

  12-
                               9/'6/56
                                                                  10/'/B6
               10    15     20
               TIMPERATURL CC)
                                            5     10     15    2:
                                                  TEMPERATURE cc;
      Fig.    2.13    Thrush Lake water temperature profiles.
                                 26

-------
 e *-
   8-
   10
                              S/17/84
   2-




I :
*~r

x
   8-
   10
                              6/lt/St
                                                                  7/03/8<
   J-
   8-
                              8/OS/S*
                                                                  8/26/84
 *

x

CL
UJ
O
                              S/IJ/8*
                                                                  11/02/8*
                10     15    20


                TEMPERATURE (*C)
                                 25     0      5
10    IS    20


TEMPERATURE (*C)
                                                                     25
       Fig.    2.14   Williams Lake water temperature profiles.
                               27

-------
     Volume  weighted and unweighted root  mean  square  error was less than
1C  for  all lakes  (Table  2.2) except the deepest (Lake Elmo has a maximum
depth  40  m).    This   is  mostly  due  to  small  differences  in  predicted
thermocline depth for the deepest simulated lake.   Difference between two
estimated root mean  square errors  (Ej and  3)  indicate the vertical position
of the  maximum  simulation  error.    If  E2 is greater  than  Ej,  than  the
difference  between  measured  and simulated  lake  water  temperatures  are
greater  in  the surface layers because  2  values  are  volume  weighted  and EI
values  are  not.
  Table 2.2  Coefficients for calibration of water temperature model
Lake




Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Year




1971/
1984
1988
1987
1985
1985
1986
1986
1984
Max.
depth

Hmax
(m)
24.0
4.70
41.8
8.20
21.0
11.0
34.0
14.0
10.0
Surface
area

As
(km2)
1.71
3.30
1.23
1.16
0.85
10.0
7.70
0.07
0.35
Wind
funct.
coeff.
c
H
24
24
26
26
24
27
29
24
22
Wind
shelt.
coeff.
Wstr
(-)
0.40
0.60
0.50
0.50
0.10
0.90
0.80
0.01
0.20
Attenuation
coefficient

Mw Hch
(m-i)(m2g-iChl-a)
0.65 8.65
0.65 8.65
0.65 8.65
1.00 8.65
0.50 8.65
0.65 8.65
0.65 8.65
0.65 8.65
0.65 8.65
Chl-a



(g m-3)
4-371
6-1302
3-83
18-484
1-47
11-348
1-35
2-l
3-79
Field data given by:
                  i
                  2
                  3>7>8
                  4'5
                  6
                  9
Shapiro and Pfannkuch, 1973
Osgood, 1984
Osgood 1989
Minnesota Pollution Control  Agency, 1988
Wright et al., 1988
Winter, 1980
      The average root mean square error for all lakes was  10C, and 94% (r2
= 0.94)  of water  temperature measurements  variability  was explained  by the
numerical model (Table 2.2).

      Model  coefficients  used in  the  simulations  are  given  in Table 2.2.
These coefficients  give  minimum  values  of  root  mean  square error, and
highest  value  of  r2   between  measurements  and  simulated  lake   water
temperatures.

      In  the  following sections the modified model with  the  hypolimnetic eddy
diffusivity  closure  as  described in  this section will be  referred  to  as  the
validated model.
                                     28

-------
2 4   Numerical Uncertainty  of Model After Hypoliranetic Closure

     Uncertainty in the lake water temperature  simulations was  considered in
.'rrns of all model coefficients except  maximum hypolimnetic eddy  diffusivity
15 specified in Table  1.1.    To  first-order  the  uncertainty  in lake  water
u-mperature depends on the  uncertainty in the  model coefficients, and on the
sensitivity of the lake water  temperatures  to  changes in  the coefficients:

     PT = E{[T  - TJ[T - T]tr) s


          E{[T(u) + ^ (u - u) - T][T(u)  + ^ (u - u)  - T]tr) =
                     dn                       dn

              (u -  n)}[~ (u -  u)]tr}=
             dn          dn
           cm                         on
           OM     on

where  P   is  the (m x  m)  covariance  matrix  of  the simulated  lake  water
temperatures,  m  is the total number  of  discritized lake control  volumes, E{.}
is  the  mathematical  expectation,  T is the mean lake  water temperature, (tr
is  the  transpose,  u is the vector of the n  coefficients,  Pu  is  the  (n  x  n
                                        frr
covariance  matrix of system coefficients,  -  is the (m x n) sensitivity matrix
of partial  derivatives of  the  lake water temperatures  with  respect  to  the
coefficients.    Sensitivity  matrix  is estimated using  the  influence  coefficient
method (Willis and Yeh,  1987).

     Data for the system coefficients covariance matrix  are given  in  Table
2.3.    These  values  were chosen  to  be in  the  range  of theoretical  and
simulated values (Tables  1.1 and 2.2),  and  to  have coefficients  of variation
(standard deviation/mean)  equal  to  0.3.    This  value  is  chosen  because
first-order uncertainty analysis  could  be  questionable when  the  coefficient of
variation of a nonlinear function increases above  0.3.
                 Table 2.3 Coefficients for uncertainty analysis

Coefficient           Lake Calhoun      Williams Lake       Cedar Lake
                   mean    st. dev.    mean     st. dev.   mean     st: dev.
/__ .o
/iw lm lj
/A:h (m2
\vstr
c
1
g-tChla)


0.65
8.65
0.60
24.0
0.20
2.65
0.18
7.20
0.65
8.65
0.20
20.0
0.20
2.65
0.06
6.00
0.65
8.65
0.60
24.0
0.20
2.65
0.18
7.20
                                     29

-------
     Three  lakes  are  selected  for  the lake  water  temperature  uncertainty
estimation.  Lake  Calhoun is a eutrophic,  deep (24 m  maximum depth) lake,
Williams Lake  is  oligotrophic, and has maximum depth close  to  the median
depth of  3002  lakes in  Minnesota  (Fig.  2.9),  and  Cedar Lake is  a  highly
eutrophic shallow (4.7 m  maximum depth)  lake.

     Standard  deviations  of smoothed  simulated  epilimnion and  volume
weighted average  hypolimnion temperatures  are  given  in Figures  2.15, 2.16,
and  2.17.    Although  high  variability  in  model coefficients  was  imposed,
maximum standard  deviation in epilimnion  temperatures  was less  than 1C,
and  less   than   1.5 C  for   the   hypolimnion  temperatures.      Epilimnion
temperatures are  most  sensitive  to the wind function coefficient for  all three
lakes.   In  the  shallow  and  well  mixed  Cedar  Lake  the  wind  function
coefficient  is  the  only  one  that  significantly  contributes to  lake   water
temperature uncertainty.   The lowest  variability of lake  water temperature
uncertainty  is  associated   with  radiation   attenuation   by  phytoplankton
(Chlorophyll-a).    Variability  in  water  attenuation   and  wind   sheltering
contribute less  t>  uncertainty in epilimnion lake water  temperatures  than  the
wind  function  coefficient.    Volume  weighted  hypolimnion  temperatures
displayed  higher  uncertainty than  epilimnion  temperatures.   For  Williams
Lake  and  Lake Calhoun, all three  coefficients  i.e.  water attenuation, wind
sheltering  and  wind  function coefficient  significantly contributed to  the lake
water  temperature   uncertainty.     Schindler  (1988)  pointed   out  that   in
oligotrophic  lakes  dissolved  organic  carbon  is  one   of  the  major  light
attenuating  factors.


2.5  Accuracy of the Regional Model After Implementation of all Changes

     The  Number  of calibration coefficients was reduced from four to zero.
Functional relationships substituted into the  model in Equations 2.2, 2.3,  2.4,
and  2.5.     The   model  output   was  compared  with  water  temperature
measurements  in  nine selected representative lakes.   Simulations  started with
isothermal  conditions  (4eC)  on  March 1  and  progressed  in  daily  time steps
until November 30.    Quantitative  measure of the success of the  simulations
and differences  between the regional model and the validated model of  section
2.3 and 2.5 are given  in Table 2.4.  The average weighted and  unweighted
root  mean  square error  was 1.1  C  (16.5  C average  measured  lake  water
temperature).    Ninty  three percent of  measured  lake  water  temperature
variability   was explained  by  the  numerical  simulations  (r2=0.93).     The
regional model has  in average 0.15 C higher  temperature root mean  square
error.

     One  example  of the  daily  simulated isotherms  for  the regional  and
validated model (section 2.3)  is given in Fig. 2.18.   Both  models simulate
onset of  stratification,   mixed  layer depth  and  water  temperatures  in  a
virtually identical  way.
                                     30

-------
                         LAKE CALHOUN
       - 1
wind  sheltering
wind  function
water attenuation
chlorophylla  attenuation
                                                    EPILIMNION
     2.0-
.-J
     0.0-
wind  sheltering
wind  function
water attenuation
chlorophylla  attenuation
                                                    HYPOL1MNION
           MAR   APR   MAY   JUN'   JUL   AUG   SEP   OCT   NOV
  Fig.   2.15   Standard deviations of estimated lake water  temperature
               uncertainties.
                                    31

-------
                             WILLIAMS LAKE
     2.5-
     2.0-
o:
UJ
Q_
     I.O-i
               wind sheltering
               wind function
               water attenuation
          	 chlorophyll-a attenuation
     0.0
2.0-
               wind sheltering
               wind function
               water attenuation
               chlorophyll a attenuation
                                               EPILIMNION
                                               HYPOLIMNION
     0.0   .-Ti --IT
           MAR   APR   MAY   JUN   JUL   AUG   SEP   OCT   NOV
 Fig.   2.16   Standard deviations of estimated lake water  temperature
               uncertainties.
                                  32

-------
                         CEDAR LAKE
     2.5-
^   2.0H
     1.5-
n
UJ

     1.0-
     0.5-
     0.0-
cm
UJ
CL
2
UJ
t
     2.0-
     1.5-
     1.0-
Ld
<   0.5H
     0.0-
          	 wind sheltering
          	 wind function
          	 water attenuation
          	chlorophylla  attenuation
               wind  sheltering
               wind  function
               water attenuation
               chlorophylla attenuation
                                                   EPILIMNION
            |  ...... T  ,. ,  I j  i .  . -j I  :  I jir- , j i  :  I ,
                                                   HYPOLIMNION
          MAR   APR   MAY   JUN  - JUL    AUG   SEP   OCT   NOV
     Fig.   2.17   Standard deviations of estimated lake water  temperature
                  uncertainties.
                                 33

-------
                          SIMULATION YEAR -  1986
         14-r
              MAR   APR   MAY   JUN   JUL   AUG   SEP   OCT  NOV
Fig.   2.18   Simulated temperature  (isotherm)  structure in Thrush Lake.  Top
            shows  results  from  validated  model  and bottom shows  results
            from regional model.
                                   34

-------
                      Table 2.4  Quantitative measure of the success of the simulations - Regional model
CO
tn
Lake

Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Average
Year

1971
1984
1988
1987
1985
1985
1986
1986
1984

Regional model
Tm
CC)
14.37
20.64
13.94
24.40
14.37
, 20.14
11.80
11.91
17.26
16.54
T8
(C)
14.44
20.68
14.31
23.90
14.90
20.09
12.61
12.54
16.57
16.67
E,
CC)
1.02
1.07
1.83
0.87
1.24
0.68
1.24
0.95
1.26
1.13
Ej
(C)
0.89
1.15
1.93
0.89
1.03
0.68
0.99
0.97
1.25
1.10
H
0.96
0.91
0.90
0.89
0.95
0.94
0.92
0.95
0.95
0.93
regional model -
AT8
(C)
-0.08
-0.18
0.22
-0.23
0.38
-0.03
0.64
0.63
0.20
0.17
\
AEt
0.16
0.13
0.06
0.07
0.38
-0.10
0.35
0.05
0.18
, 0.14
Differences
- validated model
AE2
(c)
0.10
0.16
0.13
0.07
0.24
-0.05
0.20
0.06
0.18
0.12
w
-1
-2
-2
-1
-2
2
-1
-1
-1
-1

-------
2.6  Conclusions

     A  lake  specific  water  temperature  model   was  generalized  for  the
application  to  a  wide  range  of  lake classes  and  meteorological  conditions.
Functional relationships  which  differentiate lakes on a regional  rather than on
an individual basis were developed.

     Hypolimnetic  eddy  diffusivity  was  estimated  as  a  function  of lake
surface  area, and  stability  frequency.   Equation  2.2 extends  Ward's  (1977)
analysis  to  a   wider  range  of  lake  geometries.    Although  the  proposed
relationship is  a significant  simplification of the turbulent diffusion  processes
taking  place in the  hypolimnion, it  was  found to  be useful  in the seasonal
lake water temperature modeling.

     Total  attenuation coefficient  was estimated as a function of Secchi depth
(Fig.  2.4).   Secchi depth  is  chosen because  it  can be  measured  easily  and
values  are commonly available.

     Wind  sheltering and wind function  coefficient increase with surface area
(fetch) of the lake (Figs. 2.5 and 2.6).  The wind  function coefficient increase
is  very  likely an additional  adjustment of the wind  velocity coming from land
over the lake surface.

     Uncertainty   analysis   revealed  moderate sensitivity  of  simulated lake
water temperatures to the  variability of  individual model coefficients.   This
could be  due  to   the  high  thermal  inertia of the  water especially for  the
seasonal   lake   water  temperature   modeling.      Nevertheless   epilimnion
temperatures showed 1C  standard  deviations  due  to  the  wind  function
coefficient  variability.  Water  attenuation, wind function  and wind sheltering
coefficients  equally contribute to  the  hypolimnetic  temperatures variability in
an oligotrophic lake.

     The proposed model has  practical  application in lake water  temperature
modeling,  especially  in  lakes  where  measurements are  not  available.    The
regional model  simulates  onset  of stratification, mixed layer depth, and water
temperatures well.   Average temperature mean square  error was  1.1 C,  and
93% of  measured  lake  water  temperature  variability  was  explained  by  the
numerical simulations over  a wide range of lake classes and trophic levels.
                                      36

-------
        3.   Propagation of Uncertainty  Due to Variable

     Meteorological Forcing  in Lake Temperature Models


     Propagation of uncertainty in  lake temperature models is studied  using  a
-ctor state-space method.   The output uncertainty is defined as  the result of
.rviations  of  the  meteorological  variables  from  their  mean values.    The
ir.ilysis  is applied to systems with correlated and uncorrelated meteorological
> .\nables.    Surface  water   temperatures  are strongly  affected  by  uncertain
~:<:teorological   forcing.     Air  temperatures  and  dew  point  temperature
>.:ctuations have significant effect on lake temperature uncertainty.   Ignoring
  rrelation  in   meteorological  variables underestimates  uncertainties  in  lake
.'rnperature  estimates.   Long-term average  water temperature  structure  in
akos can be estimated by  computer model simulation for just one  year when
-suits  from  a  statistical  analysis  of  meteorological  variables  are  used  as
nput.   The analysis presents a useful alternative for  the study  of long-term
ivcrages  and  variability of water  temperature   structures  in lakes  due  to
variable  meteorological  forcing.


31   Introduction

     It  was  shown in Chapters  1 and  2  that vertical  water  temperature
profiles  in  lakes are related to meteorological  variables  by heat  transport
yquations which apply basic conservation principles.  Atmospheric conditions
.ire the  driving force for heat transfer through a lake water  surface.   Surface
water temperatures  of lakes are  primarily  related to the meteorological  forcing
and secondarily to lake morphometries  (Ford  and  Stefan,  1980a).

     Observed   meteorological  variables  used   in  lake  water  temperature
modeling (Harleman and Hurley,  1976; Ford  and  Stefan,  1980b)  such as solar
radiation, air and dew  point  temperature, and  wind speed are usually single
realizations   of  the weather  process  for   a  particular  year.    For  lake
management  purposes   and  decision  analysis  we  are  interested  in  mean
temperatures as well as expected  ranges of water temperature variation  due to
the weather  variability over a  longer period  of time.    Deterministic  lake
water temperature  models  cannot   provide  such information  from a single
model simulation  for  a particular  year.    The  stochastic  alternative  is  to
consider meteorological variables as  random variables with estimated statistical
properties  in terms of  first  and second  moments,  and  correlation  structure.
First and second moment of lake temperatures can then be  predicted  from  a
single mode  simulation.

     Lake   water  temperature  models   are  nonlinear   dynamic   systems.
Approximation  techniques,  for  obtaining  the second  moment  of a  dynamic
system output  from the moments of its input have been employed in the area
of groundwater hydrology  (Dettinger  and  Wilson, 1981;  Mclaughlin,  1985;
Townley and  Wilson,  1985;  Protopapas  and  Bras, 1990; McKinney,  1990).


                                    37

-------
Generally,  three  techniques  are available  i.e.  (1)  Monte Carlo, (2)  derived
distribution, and  (3)  perturbation approach techniques.   Monte  Carlo  methods
have  been  proposed in  lake  water  quality  modeling of   phytoplankton,
herbivores, nitrate, and  available phosphorus  (Scavia et  al.,  1981;  US Army
Corps  of  Engineers,  1986;   Canale  and  Effler,  1989).   Although  simple,
limitations of  this   approach  have  been  related  to  the large  number  of
simulations.    In  addition,  the  prescribed  probability  distribution  for  the
coefficients could  change in time-varying  systems.   The  derived distribution
approach is not applicable because of the complex relationship between inputs
(meteorological variables] and outputs  (lake temperatures).  The perturbation
approach  utilizes  generally  two  methods: time domain  (state-space)  methods
of the Taylor series  expansion type,  and spectral (frequency domain) methods.
As  pointed out  by  Protopapas and Bras  (1990),  state  space methods are
advantageous for  the  time variable boundary conditions.

     In  this   Chapter  we   employed  the  perturbation  vector  state-space
approach  to  propagate  uncertainty  of  meteorological  input  variables into  a
lake temperature  model.   This  study  follows  the work of Protopapas  (1988)
who  used   the   state-space   approach  for   uncertainty   propagation   of
meteorological inputs  in  a soil/plant  model.

     The question we want to address in  this Chapter  is how  to predict the
lake temperature  uncertainty due  to  the  variability  of the  meteorological
forcing in  time.   This analysis  quantifies  contribution  of each  meteorological
variable to temperature uncertainty separately.   Secondly, we  will demonstrate
that a long-term  average thermal v structure in a lake can be obtained  without
running a water temperature model for several years of meteorological data.


3.2  Numerical Model

     In  this study  a one-dimensional  lake water  temperature  model, which
has been  previously  described in Chapter  1, was  used.   Lake  temperature  is
represented  by a  nonlinear  partial  differential  equation (1.1).    Nonlinearity
comes  through  the   boundary  conditions and  hypolimnetic  eddy  diffusivity.
Analytical   solution  of   this   equation   is  possible   only   under   certain
approximations (Dake and  Harleman,  1969).   Equation  (1.1)  is  discretized
numerically (Appendix  B)  using an  implicit  control  volume  method.    This
leads  to  a system of  equations in the form

                     Ac(K(k),G) T(k+l) = T(k) +  H(k)                (3.1)

where  Ac  is  a  system  (mxm)  tridiagonal    matrix,  m is the number  of
discretized control volumes,  T(k-fl)  is a (mxl)  vector  with lake temperatures
at  time  step  k+1,  Kfk)  is  a  (mxl) vector  with   lake  eddy  diffusivity
parameters; note  that K(k)  = f(T(k),  Ws(k)), Ws is a wind  speed, H(k)  is a
(mxl)  vector function with  source term  parameters, and G is a (mxl) vector
with lake geometry parameters.  Boundary conditions  are treated through the
source term.  The control  volume approach, satisfies the heat  balance in the
computational domain regardless of the number  of discretized control  volumes
(Patankar, 1988).
                                     38

-------
     The  numerical model is  applied  in daily time  steps  using  mean  daily
-iJut's for the meteorological variables.   The required meteorological variables
i.-r  solar  radiation, air  temperature, dew point temperature, wind speed and
.::'Xtion.  Initial conditions, model setup parameters,  have to be  provided  to
.u-  the  model.

     Taylor  series  expansion   is  commonly   used  for  the linearization  of
runctional  relations  around  nominal   values.    The  function  and its   first
Derivative  must  be defined at  the nominal  point.   Expanding  equation  (3.1)
-  a Taylor's series around the nominal value  and keeping  first  order terms,
ovts a linear perturbation temperature equation.


                Ac(k) T'(k+l) = B(k)  T'(k)  + F(k) C'(k)            (3.2)

Nominal  (mean) values  and first order  derivatives evaluated at  these values
x-'i.-  denoted by  circumflex.  Perturbations of  the  water temperatures T'(k),
-.r.d meteorological variables C'(k)  are  denoted  by primes are defined as:
                 /                                                        .

               T'(k) = T (k) - T (k),  C'(k)  = C(k) - C(k)           (3.3)


     The  tridiagonal matrix Ac(k)  is  equivalent to the  matrix Ac(k) of the
                                           A         A
icterministic  temperature model.   Matrices  B(k) and F(k)  require evaluation
 :'  the  first  order  derivatives  of  all terms  in  equation  (3.1) which  contains
..\ke temperature  and meteorological  variables  at  time  step k  respectively.
                                     A      A           A
Details  about  entries  in  matrices  Ac(k),  B(k),   and  F(k)  are  given  in
Appendix  C.   Terms  with the  same  perturbation  parameter  are collected
before entries into the matrices.  Equation (3.2) can be rewritten  as
                   T'(k+l)  = 0(k) T'(k)  + #k) C'(k)                (3.4)

where ^(k) is  transition matrix 0(k)=Ac~1(k) B(k), and ^k)=Ac~1(k)F(k).

     The  first term  in  equation  (3.4)  describes  unforced dynamics of  the
system while  the  second term describes  the  variation  of the meteorological
forcing  function.     A   schematic  illustration  of  the  lake  temperature
perturbation system is given  in Fig.  3.1.   Air temperature  (Ta),  dew  point
temperature (Td),  solar  radiation  (Hs),  and  wind  speed  (Ws)  are forcing
meteorological  functions.   A transition matrix  ^(k)  connects  the  state  of  the
system between time  steps.


3.3   First  and Second Moment Development

     Taking   the  expected  value  of  equation (3.1)  yields  the  first  order
estimate of the mean of water temperature

              E [T(k+l)] =  T(k+l)  =  Ac-i(k) (T(k) + H(k))          (3.5)
Notice that  the first  order estimate of the mean water temperature  is exactly
the value  obtained  through the deterministic  approach.
                                     39

-------
Ws(k)
 Fig. 3.1    Schematic illustration of the lake temperature perturbation
            system.
                                 40

-------
     A recursive, solution of equation (34) is
     T'(k) =  #k,0)T'(0)  +     #k,n+l) #n) C'(n);    T'(0)=0
                            n=0
     T'(k) = S   #k,n+l) V
-------
    ST/T,(k) =     ^(k,n+l)^n)Muc(k)k)V
-------
     Lake  temperature covariance propagation  is  calculated  in  the following
       (l)  Set  isothermal  (4C ) initial steady  state conditions  for  lake water
 *.;.<:ratures   in  spring,   initialize   covariance   matrix   of  meteorological
,-:v.ubations; (2) read meteorological variables,  mean and perturbation values,
 :  the  next time  step;  (31  using mean  values  for  meteorological  variables,
  mpute  first  moment of  lake  temperature  profile for  the next  time  step
                                                              A
' r-iion (3.5), store matrix Ac(k); (4) compute matrices B(k),  F(k), i.e.  first
 -' in  derivatives with respect to  the perturbed  meteorological  variables  and
T.:mated  lake  temperatures (Appendix C); (5) calculate matrix N(k)  for the
  -dated  case  (Equation 3.12);   (6)  compute transition  matrix  0(k),  and
>ii.  (7)  calculate  additional term  P(k)  (Equation  3.13)  for  the  correlated
 ii-c. (8) propagate  and  store temperature covariance  matrix ,   ,  for the

vtt time  step  (correlated case Equation 3.14, uncorrelated Equation 3.10);
 -   store transition  matrix <(k),  and ^(k), if correlated case,  store in addition
 k),  and S(k);  (10) go to  step 2  if  last day of simulation is not reached.


: *   Lake Calhoun  Application

     The  test lake,  Lake  Calhoun, is  a  temperature zone dimictic lake.   The
 lie is eutrophic with maximum depth of  24 m, and  surface  area of 1.7 km2.
Vu-teorological  data used  are from   the   Minneapolis  St. Paul International
i.-rport,  located 10 km from the  lake.

     Daily  meteorological  data  time series  (1955-1979),  averaged  over  25
 cars, are given in Fig.  3.2.   Long  term means of solar radiation,  dew point
.-mperature,  and air temperature  display  typical  seasonal cycles.   Means  are
increasing till the end of summer and decreasing towards  fall.  Perturbations
inandard  deviations) for  meteorological forcing  variables are  also obtained  by
direct data processing.  They  describe weather  variability over  25  years for a
particular  day.    Standard  deviations were higher  in spring   and  fall than  in
lummer (Fig. 3.2).

     The  time series for  each meteorological variable is reduced to  a residual
S'-ries  by removing  periodic means and standard deviations  as pointed out  by
l-jchardson (1981).   The dependence among the  meteorological  variables  was
described  by  calculating   cross  correlation coefficients  of  the  residual  time
cries.   The serial correlation coefficients for  time lags up to  3  days are given
:n  Table  3.1.   The  serial correlation coefficients for the one  day lag  were
significant for air temperature (0.69)  and dew point  temperature  (0.66).   A
significant cross correlation coefficient (0.8)  was  calculated for  zero time  lag
ithe same day) between  dew  point  temperature and air temperature.   Other
meteorological variables were  uncorrelated for the  same day.

     The first  order  estimate of the  mean  and  covariance temperatures  is
constrained to parameter  perturbations within only the linear region about the
model  trajectories.    Linear  approximation could  be  questionable   when  the
coefficient  of variation  for  the  parameter  of a  highly nonlinear function
increases above  0.3  (Gardner  et al.,  1981).   Average  coefficients of variation
for  input  meteorological   variables   are:  air   temperature  0.13,   dew  point
temperature 0.17, wind  speed 0.33,  and  solar  radiation  0.37.   Although  the
                                     43

-------
fN
 I

 E
 o
 O
 h-
 <
 O
a
3
o
GO
UJ
a:
LJ
a
2
 O
 a


 UJ
 Q
700



600-


500


too



300-



200-



100-


  n


 35


 30-


 25-


 20


 15-


 10


  5-


  0-


 T5-


-10
           I I 1  I I  ' *  'I '

           overage

           standard deviation
                            i  " :;i i-   -
         M'R    MAV   JUN    JUL    AUG    SIP    OC1
            ovmoge

            standard devialion
         '  i ''~r       i      i  - i i       i      i  
         APR    MAY   JUN    JUL  . AUG   SEP   OCT
                                                                                                          (
                                                           AI'H    MAY    JUN    JUL   AUC    SL'P   OCI
                                                            ..".:'',V.;i.                          .1. > '!.-.
                                                          <  '( V ^ '' (' : ^'!"\ft,-.l.i..;: ;V /Vi-V.A'-'^^A:"'^ '
-------
         Table 3.1      Correlation  coefficients of  daily  meteorolo/;icnl  viiriahlrtt  f..i
                         Minneapolis-St. Paul,  1955-1979.

Solar
Radiation
Air
Temperature
Dew Point
Temperature
Wind
Speed
Solar
Radiation
Time Lag (days)
01 2
1.00 0.39 0.14
0.18 x 0.17 0.11
-0.25 -0.14 -0.06
-0.16 -0.05 -0.04
Air
Temperature
Time Lag (days)
0 1 2

1.00 0.69 0.38
0.80 0.54 0.26
0.11 0.08 0.07
Dew Point
Temperature x
Time Lag (days)
0 1 2

0.80 0.58 0.29
1.00 0.66 0.33
0.09 0.07 0.06
Wind
Speed
Time Lag (days)
0 1 2



1.00 0.38 0.18
Cn

-------
 solar  radiation had the  highest  variability  note  that it is linearly  related
 through the source term to the  water  temperature equation (Equations  1.1,
 1.2, 1.3, and  1.4)

 3.4.1     First moment analysis

      The  nonlinear lake temperature model was  used for the first  moment
 temperature estimation.   Model  setup parameters  which  are  basically  related
 to  lake geometry  have been estimated by comparing  model  simulations with
 measurements  (Chapter  2).   The  standard  error  between measurements  and
 simulations was about 1.0 C.   The  error  is  mostly associated  with small
 differences between measured and  predicted thermocline depths.

      Long  term average temperature structure in  Lake Calhoun was obtained
 using  two  different methods.  In the first  method,  the  lake temperature model
 computed  the vertical  temperature  structure in the lake for each  of  twenty
 five years  (1955-1979), separately using  daily values for  meteorological data.
 The results  of these twenty  five years  of simulated  lake temperatures were
 statistically Analyzed in terms of mean  (Teav)  and  standard  deviation (oeav)
 for  the particular  day.   In the  second  method,  twenty five years of daily
 meteorological data were first  statistically analyzed to  provide daily  means
 and standard  deviations.    This  averaged meteorological  year was  used in a
 single   simulation   run  to obtain  the average daily  water temperature (Tav)
 throughout  a season.

      Epilimnetic  and  hypolimnetic  lake  temperatures  obtained  by  these  two
 methods are compared in Fig.  3.3.   Epilimnion  temperature  is defined  as  the
 temperature of the upper isothermal  (mixed)  layer.  Hypolimnetic temperature
 is a volume  weighted  average temperature below  the  upper  isothermal layer
 down   to  the  lake bottom (Equation  2.5).    Nearly identical  temperature
 distributions  were obtained by  the two  methods.   Maximum  difference  was
 less than  1C at  any time  of the season.  Isotherms obtained by the  two
 methods are compared for the entire period  of  simulation in  Fig.  3.4.   Onset
 of stratification and mixed layer depths can  be  seen  to be nearly identical.

 3.4.2     Second  moment analysis

      Uncertainty in the lake temperatures is measured by the variance of the
 model  output.   Temperature  covariance  propagation was calculated by using
, the  proposed  vector state-space  perturbation model.  Two cases  were
 considered:    (1)   uncorrelated   meteorological      variables,   (2)   correlated
 meteorological  variables.    "Uncorrelated" means  that   daily  meteorological
 variables were independent of each  other at any time.    "Correlated"  means
 that  a correlation  between air, and  dew  point  temperature  at  zero and   one
 day time  lag existed.    These  two meteorological  variables  were  considered
 because they  had  significant correlation, and  as.  will be  shown   later,  this
 resulted in a significant contribution to lake temperature uncertainty.

      The  Long-term average  temperature structure in  the lake  was calculated
 using  the   second  method  in the first  moment  analysis.   Perturbations for
 meteorological forcing  variables  are given in  Fig. 3.2.
                                      46

-------
    35-
     0-
O
o
    3CH
             'eav 1955-1979
            T,
eav 1955-1979

av
     0	r
                                  EPILIMNION
                                               KYPOLIMNION
         APR    MAY     JUN     JUL    AUG    SEP     OCT
  Fig. 3.3    Estimated long-term average epilimnion and hypolimnion

            temperatures.
                             47

-------
24	^
      APR    MAY    JUN     JUL    AUG    SEP    OCT
    Fig. 3.4   Long-term average isotherms in Lake Calhoun.
                         48

-------
      Standard   deviations   of   simulated   epilimnion   and    hypolimnion
 '.rmperatures  are  given in Figs.  3.5 and  3.6, respectively.   Contributions  by
 -".-rturbations  of individual  meteorological  variables  perturbations  as  well  as
 .he  total  contribution of  all  perturbation  variables  were  calculated  with
  or related  and uncorrelated input variables at  a daily timestep.   Air and dew
 >--:nt temperature  contributed the most to the temperature uncertainty, while
 viiar radiation and wind  speed had smaller effects.   Furthermore,  the overall
 uncertainty in water  temperature was found  to  be larger in the case of the
 Correlated  daily  process  than in the  uncorrelated  one.    Uncertainty  in lake
 ater  temperature  varies with  time,  since sources of  uncertainty  vary with
 iirae.     These   sources   are,   the  sensitivity   of  lake   temperatures   to
 meteorological variables as well as  the  amount  of the  uncertainty  concerning
 these variables. At  the beginning  of the simulated period uncertainty  was  set
 ;.o  zero  since  initial  conditions  were  considered  perfectly known.    Isothermal
 initial conditions of 4"C (after ice thaw) April   1  are appropriate for the 45
.latitude.   Although isothermal water  at 4'C  may not  exactly exist on April
 i, thermal inertia of  the water makes  summer predictions insensitive to  initial
  nnditions  if  a starting date at  or before "ice-out" is  chosen  (Chapter  2).
 I'hree periods can/be distinguished  in Fig. 3.5  : a steep rise in  temperature
 .mcertainty  in spring,  more  or  less  constant  uncertainty after onset  of
 * i ratification  in  summer,  and  decreasing  uncertainty  in  fall   when  lake
 .cmperature  is driven towards isothermal conditions.  Temperature  uncertainty
 :s decreasing  in fall when observed meteorological variables and estimated lake
 *ater  temperatures are both decreasing.  First  order derivatives with respect
 u  the  lake  temperatures and meteorological  variables  are evaluated  at these
 bserved and  estimated  values respectively.   Thus, they have less weight in
 uncertainty propagation.

      Uncertainty propagation for deep hypolimnetic  temperature  (1m  above
 lake sediments)  is  given  in  Fig.  3.6.    In spring and fall, during  well-mixed
 conditions  (overturn periods), standard  deviations of 0.4  C (correlated case)
 and 0.3 C   (uncorrelated case)  are calculated.   During stable  stratification,
 uncertainty  was  not  significant.    This is a result  of  the fact   that  Lake
 Calhoun has  no  significant  continous  point  inflows (tributaries).    .Summer
 temperature  in the hypolimnion  was  determined by  mixing  events  in spring,
 and remained almost  constant  throughout  the fall  overturn (Ford and Stefan,
 1980a).   In  lakes  with  point inflows  this would not   be  the  case,  due  to
 plunging flow  phenomena.

      Vertical  profiles of  the first moment lake  temperatures, plus or  minus
 one  standard  deviation interval,  are shown in Fig.  3.7.   Spring  (April)  and
 fall  (October) indicated periods when  uncertainty propagates throughout  the
 entire lake depth.   These are the periods of weak stratification  or well  mixed
 conditions.    Uncertainty was  decreasing  with   depth.    After  the onset  of
 stratification   estimated  uncertainty  was  much  more  significant  for   the
 epilimnetic layer  than  for the hypolimnetic layer.   For the same period of
 time, deep water had insignificant lake  temperature uncertainty.

      The  first moment epilimnion  temperature  estimates plus  or  minus  one
 standard deviation obtained by  two  different approaches for the  1955-1979
 period  are compared  in  Fig.  3.8.    In the  first  approach  the deterministic
 water temperature model" was  run  for 25  years using  daily  meteorological
                                      49

-------
      STANDARD  DEVIATION  OF SIMULATED EPILIMNION TEMPERATURES
     6.0-
O   5.0-
o
LU
rr
Z)
     4.0-
 rr
.UJ   3.0-
 o_
 UJ
 "~   -2.0-
 rr
 UJ
 S   1.0H
     0.0-
          	 all inputs
          	air temperature  ,
          	 dew point temperature
          	 solar radiation
          	wind speed
                                  UNCORRELATED DAILY PROCESS -
O   5.0-
o

UJ
 cr
 LJ
 Q_
 LJ
             all inputs
             air temperature
            -- dew point temperature
             solar radiation
             wind speed
                                  CORRELATED DAILY  PROCESS
     0.0
            APR     MAY
                           JUN     JUL
AUG
.SEP
OCT
    Fig. 3.5    Standard deviations of estimated epilimnion  temperature
              uncertainties.   Contributions by several  meteorological
              variables and totals are  shown.
                               50

-------
    STANDARD DEVIATION OF SIMULATED DEEP  HYPOLIMNION  TEMPERATURE
     O.c
LU
Ct
ID

<
C
UJ
CL
cz
UJ
               all inputs
               air temperature
               dew point temperature
               solar radiation
                              UNCORRELATED DAILY PROCESS
          	  wind  speed
O.-H
0.2-
     O.C
O
o
               all inputs
                              CORRELATED DAILY PROCESS
       	  air temperature
       	  dew point temperature
       	  solar radiation
          	 wind  speed
           APR
                                                       OCT
 Fig.  3.5    Standard deviations  of estimated, deep water temperature
            uncertainties.   Contributions by several meteorological
            variablesand totals are shown.
                                 51

-------
         4-





         8-





        12-





        1S-





        20-
        24
                         4/10
                                                        6/08
         4.





         a-





        12-



         ,'

        16-





        20-





        24
                         S/M
      I
      t^
      Q.
      U)
      o
 4-





 s-





12-





16-





10-
        24'
                                                                B/24
12-





16-





20-





34
                                                                 10/28
              5   10   15
                              CC)
                                          10    IS   23    25   X


                                           TEMPERATURE fC)
Fig. 3.7    Long-term average temperature profiles plus  or minus  one

            standard  deviation  in Lake Calhoun.
                                   52

-------
Cn
CO
                 40-
             ^  35-1
             O
             o
             "~"  30-
             Ld
             cr
             P  25-
             r
             <
             o:
             UJ  20-
             CL
                 10-


                  5-


                  0
Teov (1955-1979)-f(~)^"eav
                                                              EPILIMNION
                      I  I  T  I ^|  |  i  i|  I  T^ 1^I  I ^  i  1  j^ t  1   I  |  I  F  I T  T

                      APR     MAY    JUN    JUL    AUG     SEP    OCT
                     Fig. 3.8   Epilimnion^temperature long-term average plus or minus one
                             standard deviation.

-------
data.     Long  term  average  (Tgav)  and  standard  deviations  (aeav)  were
estimated from the  simulated lake  water  temperatures  over  the  25  year
period.   In  the  second approach the  long  term  average (Tav)  temperature
structure in  the  lake was  estimated using the  method  described in Section
5.1.   Water  temperature variability (<7av) was estimated  using the  proposed
perturbation   model.     Results  shown  are  for  correlated   meteorological
perturbation  variables.   The  maximum difference  was less than 2C for the
range of 23 C variability.


3.5  Conclusions

     A first   order  analysis of uncertainty  propagation  in lake  temperature
modeling  has  been  made.    The  source  of  the  uncertainty  is   variable
meteorological forcing which  enters  the lake  temperature equations  through
the source term  and boundary conditions.   The  analysis presents  a  useful
alternative  for the study of long-term  averages  and variability  of  temperature
structures in  lakes due to variable meteorological forcing.

     The analysis applied  herein  can  be  applied  to  systems with  correlated
and uncorrelated -meteorological parameters.   The main findings are:

     (1) Long-term  average temperature  structure  in lakes can be  estimated
by  using the  results of  a  statistical analysis  of long-term  meteorological
variables as input  in  a computer model simulation  for just  one  year.

     (2) Air  temperature  and  dew  point  temperature have significant effect
on lake temperature  uncertainty.

     (3) Epilimnetic  temperature  uncertainty  has  three  distinct  periods   :
steep  rising uncertainty in  spring, steady uncertainty in summer, and falling
uncertainty in  fall.   The maximum  standard deviation of 4C  of epilimnetic
temperature  uncertainty was  estimated  in the  summer  for  the 25 year  a
period.

     (4) Hypolimnetic  temperatures were  not strongly affected by  uncertain
meteorological forcing.  Standard  deviations  of  less than  1C were  estimated
in spring and fall  during the overturn periods.

     (5)  Ignoring  the   correlation   of  air   and  dew   point   temperatures
underestimates  uncertainties in lake temperature  estimates.   Accounting for
correlations gives  better  agreement  with lake water temperatures  obtained by
25 years of estimated lake temperatures.
                                     54

-------
              4. Case Studies of Lake Temperature

        and Stratification Response to Warmer Climate


     The  impact of climatic warming on lakes will  most likely have serious
implications  for water resources  and water quality.   Rather than using model
predictions of greenhouse warming, this chapter looks at the changes in heat
balance  and  temperature profiles in a particularly  warm year (1988) compared
to a more normal one (1971).  The comparisons are  made for three   different
morphometrically  different lakes located  at  45 north  latitude  and 93  west
longitude  (North Central USA)  and for the  summer  period  (April 1  to
October 31).  Water temperatures are daily values  simulated  with  a  model
driven  by  daily  'weather  parameters and   verified  against several  sets  of
measurements.   The results  show  that in the warmer  year epilimnetic water
temperatures were  higher,  evaporative  water  loss  increased,   and  summer
stratification occurred earlier in the season.


4.1   Introduction

     A  validated one-dimensional lake water temperature model, which  has
been described in Chapters 1 and  2, was used to  study the  changes in a lake
as a result  of  different  weather conditions.   In  this  chapter use of such a
model  is  demonstrated  by   application  to   three  different  morphometrically
lakes with sparse data sets.   The  lakes are located near 45* northern latitude
and 93 western longitude in northcentral United  States.  The lakes  are Lake
Calhoun,  Lake  Elmo,   and  Holland  Lake in  the   Minneapolis/St.  Paul
metropolitan area.

     In the summer of 1988,  the northcentral region  of the  United  States
experienced  very  dry and hot  weather and  this  was selected to represent a
"warm climate"  in  this  study, while for  "normal" conditions, the  year 1971
was chosen.   1988 tied  for the warmest year in the  100-year global  record of
instrumentally  recorded  air  temperatures  (Kerr, 1989).   Uncertainty  analysis
of the effects  of variable meteorological forcing  on  lake  temperature models
indicates  that  air   temperature  has  the  most  significant effect   on  lake
temperature  uncertainty (Henderson-Sellers,   1988;  Chapter  3).    1971  was
normal in the sense that mean  air temperature from May to September  was
only 0.2* C  below  the normal from 1941 to 1970.   The  effects  of  the 1988
(warmer)  and  the  1971  (normal)   summer climate  on  temperatures  and
stratification in the three lakes are reported herein.
                                    55

-------
4.2  Method of Lake Temperature Modeling

     The test lakes,  Lake Calhoun,  Lake Holland, and Lake Elmo,  are three
temperature zone  dimictic lakes.   Water temperature data  were collected in
Lake  Calhoun in  1971  (Shapiro and Pfannkuch,  1973) and used  to validate
the model  for normal weather conditions.  For warmer conditions (1988) the
model  was validated  with data from Lake Elmo  and Lake  Holland (Osgood,
1989).   The  terrain  in which the lakes and  weather stations  are  located  is
flat and  quite  uniform with  respect to land use  (residential and park  land).
Morphometric characteristics,  Secchi-depths  and  chlorophyll-a measurements
for all three lakes in  1984 (Osgood, 1984) are given in Table 4.1.  Lake Elmo
surface area is equal  to the  median value of  970 statistically  analyzed  lakes
in the North  Central Harwood Forests  ecoregion  in Minnesota (Heiskary and
Wilson, 1988).   Lake Calhoun and  Holland  Lake have  a larger and smaller
surface area than  the median, respectively.   All three lakes  were classified  as
eutrophic.  Secchi depths  and chlorophyll-a  were  close to the  median  values
of the lakes in the ecoregion.

                         . Table 4.1  Lake data
Lake
Calhoun
Elmo
Holland
Mean
depth
[m]
10
13.4
4.6
Max
depth
[m]
24.0
41.8
18.8
Surface
area
[km*]
1.71
1.23
0.14
Volume
[10m3]
17.1
16.5
0.65
Secchi
Depth
M
2.5
2.8
2.2
Typical
Summer, Chla
[gm-3]
20
8
28
     Meteorological  data used  are from the  Minneapolis-St. Paul Internationa]
Airport located  5 to  18  miles from the studied lakes.   The  meteorological
data file contains  measured daily  values  of average air  temperature (Ta), dew
point temperature (T.'
duration have been  continuously  recorded  on Lake  Mendota (Wisconsin,  C
latitude) since  1855.  The mean  date  of ice  thaw was April 5 with  11 cU*-
standard deviation  (Robertson, 1989).    Model sensitivity  to  the  date  of :>
initial   isothermal   conditions  is  summarized  in   Table 4.3.     Epiiiranr.;
temperatures are  very well  simulated throughout  the entire  summer  per.--*
                                     56

-------
               Table 4.2  Mean Monthly Meteorological Data
      Max  Min   Aver. Diff. from
                        Normal*
                                     Td     P     W          HS
                                    [C]    [mm]  [ms'1]     [cal cm^d'1]

U'R
MAY
;r.\
-TL
U.'G
 Y. P
CT

14.9
19.4
27.4
26.5
27.5
22.9
15.7

1.7
6.5
16.5
14.3
14.2
11.3
5.8


8.3
13
21
20
20
17
/io
.0
.9
.4
.9
.1
.8
Year
0.6
-1.3
2.0
-2.3
-0.8
1.2
0.5
1971
-1
2
15
12
12
10
6

.7
.8
.0
.8
.2
.0
.7

0.9
2.6
3.0
3.2
1.4
2.3
4.6

5
4
4
4
3
4
4

.1
.6
.0
.1
.8
.1
.5

411
482 .
450
563
479
338
192
MEAN(MAY to SEPT)
   ,   24.7   12.6   18.7
  (MAY to SEPT)
       3.15   3.45
3.26
       -0.2
1.60
         10.6
       2.5   4.1
4.20    0.63  0.26
MEAN(MAY - SEPT)
      28.2   15.0   21.6      2.7       11.8     1.7    4.8

a (MAY - SEPT)
       3.42    3.23  3.29     1.20       3.03    1.16   0.23
462
72.7

APR 15.1
MAY 25.7
JUN 30.5
JUL 32.3
AUG 29.3
SEP 22.8
OCT 12.5

2.0
11.4
16.6
18.8
17.3
10.9
0.8

8.5
18.5
23.5
25.6
23.3
16.9
6.7
Year
0.8
4.2
3.6
2.9
1.6
1.0
-3.6
1988
-2.4
7.1
12.3
14.5
15.3
9.8
-0.6

1.3
1.4
0.2
0.9
3.5
2.4
0.7

4.6
5.1
4.8
4.4
4.8
4.8
4.7

469
584
654
610
497
331
284
                                           535
                                           114
*Normal is the  30-year average from 1941 to 1970
a  standard deviation
                                  57

-------
Or
OO
      Table  4.3  Differences (C) in  simulated  mean  daily epilimnetic and hypolimnetic temperatures for different
      starting dates of the model (April  1  reference)
MAY
JUN
JUL
Epilimuiou
AUG SEP OCT SEASON
MAY
JUN\
Hypolimnion
JUL AUG SEP
OCT
SEASON
Year 1971
MAR 1 -0.05
MAR 12-0.09
MAR 22-0.09
APR 10 0.05
APR 20 0.91
-0.07
-0.08
0.07
0.01 v
-0.03
0.00
0.00
0.01
0.00
0.00
0.00
0.00
-0.01
0.00
0.00
0.00
0.00
-0.05
0.00
0.02
0.00
-0.01
-0.03
0.00
0.08
-0.02
-0.03
-0.02
0.01
0.16
-0.19
-0.29
-0.18
0.08
1.77
-0.24
-0.34
-0.18
0.03
1.59
-0.28
-0.39
-0.18
0.01
1.46
-0.30
-0.41
-0.19
0.00
1.35
-0.32
-0.43
-0.19
0.00
1.26
-0.33
-0.45
-0.19
0.01
1.20
-0.28
-0.39
-0.18
0.02
1.44
Year 1988
MAR 1 0.14
MAR 12 0.07
MAR 22 0.07
APR 10 0.70
APR 20 1.60
0.00
0.00
0.00
0.03
0.08
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
-0.02
-0.03
0.01
0.00
0.00
-0.04
-0.07
0.03
0.01
0.01
-0.06
-0.15
0.03
0.01
0.01
0.10
0.24
0.48
0.01
0.29
1.42
3.14
0.49
0.13
0.30
1.46
2.93
0.51
0.10
0.32
1.44
2.79
0.52
0.10
0.33
1.43
2.77
0.52
0.10
0.33
1.42
2.67
0.44
0.08
0.28
1.40
2.59
0.49
0.10
0.31
1.43
2.82

-------
regardless of  the  starting  date of the  model.   Surface water  temperatures
''catch up"  in time.   Hypolimnetic  summer water  temperatures are good  as
long as  the  model  is  started  before  seasonal  stratification  sets  in.   Better
results are obtained if temperature is not  allowed  to drop  below  4C  after
start of  the  simulation.   Although  isothermal  water at 4C  may not exactly
exist  on  April  1, thermal inertia of the  water  make  summer  predictions
insensitive  to initial  conditions  if a  starting  date  at or before  "ice-out"  is
chosen.
4.3  Model Validation

     The model was  validated  with  water temperatures  measured in Lake
Elmo and Holland  Lake in 1988.   Eight examples  of measured and  calculated
water temperature  profiles for these  lakes  are  given  in  Figs.  4.1  and 4.2.
Actually 16  profiles were measured in each lake.   Simulations  started with
isothermal conditions (4C) on April 1 and progressed in daily timesteps until
October 31.  Model coefficients were kept at their  initially  specified  value
throughout this period.    The model  simulates onset of stratification, mixed
layer   depth   and   water  temperatures   well.     Standard   error   between
measurements and  simulations was  2.0 *C and 1.5 *C  for  Elmo  and Holland,
respectively.    This  is  mostly  due  to  small  differences  in  the   predicted
thermocline  depth.   A model validation for Lake Calhoun was made for 1971.
Measured and  calculated  water  temperature profiles  are .given  in   Fig. 4.3.
Comparison  shows  that  the  onset of  stratification, mixed  layer depth and
temperature were well predicted.   Standard  error was 1.4C.


4.4  Results and Discussion

4.4.1      Thermal  energy  budget

     Mean monthly heat balance terms for 1971 and 1988  are given in Table
4.4.   Short  wave  solar  radiation  (Hsn)  and longwave  atmospheric  radiation
(Ha)  increase  the  water  temperature,   while  evaporation  (He),   and back
radiation (Hbr) cool the water.  Conductive heat transfer (Hc) can either heat
or cool the  water.   These  five mechanisms,  mainly responsible for  the  net
heat energy input to the water,  changed from  month to month and  from year
(1971)  to year (1988).   Solar radiation (Hsn)  and atmospheric radiation (Ha)
are  only  given  once   because   they  are  the   same  for  all  three  lakes.
Cumulative  heat balance terms for both simulated periods  are given in Table
4.5.

     Under  warmer conditions (1988) more solar radiation  reached the lake
surfaces.  The cumulative difference at the end  of the  simulation period was

5000  kcal  m'2.    The  additional  available  solar  radiation  increased  the
surfacewater  temperature and   stability  (defined  as   a   density  difference
between adjacent layers) of the  water column  (Spigel et al.,  1986)  as  will  be
shown.
                                     59

-------
                  TEMPERATURE (C)
^  40 h
c.
Ul
o
      40 h
      Fig. 4.1     Lake Elmo  water temperature profiles.
                          60

-------
                 TEMPERATURE  (C)
   10  -
   20  -
             o  MEASURED
             CALOAATED
  20  -
Fig. 4.2    Lake Holland water temperature profile.
                        61

-------
                    TEMPERATURE  (C)
      o (-
     10  -
    20  -
OJ
Q
                O MEASURED
                CALCULATED
    10 -
    20  -
Fig.  4.3     Lake  Calhoun water temperature profiles in 197.1.
                         62

-------
     Table  4.4  Monthly averages of daily heat balance components [1000 kcal
H8n
Lake Calhoun
Ha Hbr He
Hc Hn
Lake Elmo
Hbr He HC
Hn
Lake Holland
Hbr He HC
Hn
                                                             Year 1971
APR
MAY

JUN
JUL
AUG
SEP

OCT
MEAN

3.89
4.60
i
4.27
5.37
4.54
3.18

1.83
(MAY to
4.39
5.76 -6.92
6.27 -7.64

7.44 -8.58
7.11 -8.85
7.20 -S.71
6.87 -8.36
\
6.18 -7.69
SEPT)
6.98 -8.43
-1.04
-1.86

-2.00
-3.30
-2.68
-2.14

-1.39

2.39
0.34
-0.08

0.07
-0.46
-0.20
-0.25

-0.42

-0.18
2.03
1.28

1.20
-0.13
0.15
0.70

-1.49

0.36
-6.88 -1.04 0.44
-7.49 -1.65 0.13

-8.44 -1.74 0.25
-8.82 -3.47 -0.46
-8.67 -2.76 -0.17
-8.35 -2.26 -0.25

-7.71 -1.55 -0.49

-8.35 -2.38 -0.10
X 2.17
1.85

1.77
-0.27
0.13
-0.81

-1.75

0.53
-7.03
-7.76

-8.63
-8.83
-fl.71
-8.30

-7.49

-8.45
-1.37
-2.22

-2.24
-3.38
-2.79
-2.10

-1.09

-2.55
0.17
-0.25

0.01
-0.45
-0.21
-0.19

-0.15

-0.22
1.41
0.64

0.84
-0.19
0.04
-0.54

-0.73

0.16
o>
co
                                                             Year 1988
APR
MAY
JUN
JUL
AUG .
SEP
OCT
MEAN

4.46
5.57
6.27
5.83
4.72
3.11
2.69
(MAY to
5.10
5.76 -7.08
6.84 -8.06
7.52 -8.99
7.84 . -9.17
7.56 -9.00
6.79 -8.23
5.54 -7.50
SEPT)
7.31 -8.69
-1.29
-2.49
-^J.40
-4.06
-3.74
-2.32
-1.97

-3.40
0.16
0.32
-0.18
0.03
-0.21
-0.27
-0.79

-0.06
2.01
2.18
0.22
0.47
-0.67
-0.92
-2.03
 
0.26
-6.96 -1.15 0.39
-7.84 -2.00 0.71
-8.87 -4.29 -0.04
-9.11 -4.14. 0.11
-8.99 -3.98 -0.21
-.24 -2.55 -0.31
-7.39 -1.86 -0.47

-8.61 -3.39 0.05
2.49
3.28
0.59
0.52
-0.89
-1.20
-1.44

0.46
-7.25
-S.21
-9.01
-9.14
-8.95
-8.16
-7.30

-8.69
-1.71
-3.07
-4.63
-4.15
-3.70
-2.21
-1.65

-3.55
-0.12 .
0.10
-0.20
0.05
-0.15
-0.18
-0.48

-0.08
1.13
1.23
-0.06
0.43
-0.52
-0.65
-1.21

0.09

-------
                                    air
                                   was
     Atmospheric  long  wave  radiation  and  back radiation from the  water
surface  are  proportional to the fourth power  of absolute temperatures.  Both
were  higher  under  wanner  conditions.    Higher  back  radiation  was   an
indication   of  higher   surface   water  temperatures   under  increased   air
temperatures and solar  radiation.

     Cumulative  evaporative  losses  resulting from  the  average 2.9'C
temperature increase  are plotted in Fig.  4.4.  Cumulative evaporative loss

higher  by  about  180,000  kcal m'2  for  the  1988  season  compared  to  1971.
This  translates into  an  additional  water  loss  of  about  0.3   m  in  1988
compared to 1971.  This loss occurred in each of  the three lakes despite their
differences  in  size  and depth.   Increased  evaporation  not  only represents  an
additional water loss  but  also contributes  to  increased natural  convection  due
to surface cooling.

     Conductive heat  transfer through  the lake  surface made only  a  small
contribution to the heat budget.  The cumulative  conductive  heat  input was
not significantly different  during  the  two  years,  but  the onset of cooling  by
convection  was delayed until August in 1988.
              /
     Net heat fluxes  on  a  monthly  time  scale are  shown in  Table 4.4,  and
on a cumulative basis in Table 4.5.   Cumulative net heat  flux (Hn)  from the
atmosphere  to the  water increased from  April to June  in 1971  and from April
to July in  1988, and  then began to decrease indicating that the lakes  received
heat for a  longer  period  in 1988  than in 1971.   The net cumulative  heat
input is  also  a  measure of heat content relative  to April 1.   The maxima of
the net  cumulative heat  input were  only  slightly  different in  1971 and  1988
(see  Table  4.5), but very  different among the three lakes because of the effect
of depth especially  surface mixed  layer depth.  Normalized  values  with respect
to depth are given  in Table 4.6.  The trend is from  higher to  lower values as
the depth  increases.   This  reflects the  thickness of the surface  mixed layer
depth relative to the total lake depth.

4.4.2     Equilibrium temperatures

     Equilibrium temperature is  defined as that  water temperature  at  which
the  net rate of  heat exchange  through  the   water  surface is  zero  and
continually  changes  in response  to   the  meteorological  conditions.    Mean
monthly equilibrium  temperatures for Lake Calhoun are  shown  in  Fig.  4.5.
These  values  were obtained  by  a separate calculation  setting the  net  heat
transfer  rate  Hn equal to zero.   Calculations  were carried out for  the  entire
year (12 months) to see  how the  dates of the  0C  crossings  and  hence the
date  of ice  formation might  shift   from  year   to  year.    Under  warmer
conditions equilibrium temperature was higher from  March to  August.  From
August  up to the ice formation in  November no difference between  the  colder
and  the warmer  year  was noticed  probably  because  the fall of  1988  was
cooler than in 1971 (see Table 4.2).  The 0*C crossings in Fig. 6 occurred at
about the same time in 1988  and 1971  indicating  that dates of ice formation
and  ice  thaw  were not significantly affected by the  heat in  July  and August.
There  could  be a larger  change in  ice  thaw  and  freeze-over  dates  if  air
temperatures  were  changed yeararound, not  only in summer  as  in  this  case
study.
64

-------
         Table 4.5   Cumulative heat  balance components   [1000 kcal m"2]
Lake Calhoun
Hgn Ha Han
He
Hc
Hn
Lake Elmo
He Hn
Lake Holland
He Hn
                                                    Year 1971
APR

MAY
JUN
JUL
AUG
SEP
OCT
117

259
387
554
694
790
846
173

367
590
810
1034
1240
1431 
-208

-444
-702
-976
-1246
-1497
-1735
-31

-89
-149
-251
-334
-398
-441
10

8
10
-4
-11
-18
-31
61

101
137
132
137
116
70
-31
1
-82
-135
-242
-328
-396
-444
65

122
175 ,
167
171
147
93
-41
\
-110
-177
-282
-369
-431
-465
42

62
87
81
82
66
44
01
                                                    Year 1988
APR
MAY
JUN
JUL
AUG
SEP
OCT
134
306
495
675
822
915
998
173
385
610
863
1088
1291
1463
-212
-402
-732
-1016
-1295
-1542
-1775
-39
-110
-248
-374
-490
-559
-620
5
15
9
10
4
-4
-29
60
128
134
149
128
101
38
-35
-97
-225
-354
-477
-553
-611
75
176
194
210
183
147
102
-51
-M7
-286
-414
-529
-595
-647
34
72
70
84
68
48
11

-------
                SIMULATION PERIOD  4/01-10/31
  800000-

  700000-

5 600000-

  500000-
   100000


       0


^ 700000-I


J5 600000-1
*

g 500000 -|

S      J
5 400OOO-1
o
$
S 300000-
LU
>
3 200000-
3

0 100000-


       D


^ 700000-
"2

5 600000-
2.
to
cl 500000-
J
LJ
>
^ 400000-
cr
O
CL
 30C-000-
LJ

3 200000-
3
2
O
0 100000-


       0
          1971  losses (kcol/m2)
          1988 losses (kcol/m2)
          1971  losses (m)
          1988 losses (m)
           LAKE CALHOUN
          ""["   "*   |    '   I   _
          1971 losses (kcol/mz)
          1988 losses (kcol/m2)
          1971 losses (m)
          1988 losses (m)
           ~\	'  T	'	f  ~
          1971 losses (kcol/m^)
          1988 losses (kcol/m2)
          1971 losses (m)
          1988 losses (m)
                                         LAKE HOLLAND
                                                                   1.2
                                                                  p.1
                                                                  j-i.o
                                                             -0.9 ~

                                                            E-0.8  
                                                            E-o.s  >

                                                            r"  I
                                                             rO.3  3
                                                             -0.2
                                                            =-0.1
                                    1-1.1
                                                             r-8  $
                                                             t     o
                                                             r-0.6  0
                                                             -0.4
                                                             -0.3
                                                             -.,
                                                             E-o.i
                                                              o.o
          APR
                      MAY
JUN    JUL    AUG   SEP   OCT
Fig.  4.4      Cumulative evaporative  losses (simulated).
                                 66

-------
                            30-
o>
                        O
                        o
bJ
0



i
Ul
Q.
2
UJ
                            20-
                            15-
                            10-
                             5-
                            -5-
                           -10-
                                     1988-89

                                     1971-72
                             FEB  MAR APR MAY JUN  JUL AUG  SEP  OCT  NOV  DEC  JAN
                                Fig. 4.5    Mean monthly equilibrium temperatures (simulated).

-------
Table 4.6   Net cumulative heat input  (content) per meter of average depth
(1000 kcal m'1)
Lake Holland
: (4-6 m)
Lake Calhoun
(10 m)
Lake Elmo
(13.4 m)
                             Year  1971
APR
MAY
JAN
JUL
AUG
SEP
OCT
11
16
22
20
21
17
11
6
10
14
13
14
12
7
5
10
14
13
14
12
7
                             Year 1988
APR
MAY
JUN
JUL
AUG
SEP
OCT
8
18
18
21
17 
12
3
6
13
13
15
13
10
4
6
14
16
17
15
12
8
                                     68

-------
4.4.3     Vertical  mixing/onset of stratification

     Surface  mixed  layer depths  are  shown  in  Fig. 4.6.   The mixed layer
depth is  defined as the thickness of the upper isothermal layer.  Large mixed
layer depths  at the beginning  and at the end of  the simulated period  indicate
spring  and fall overturns.  After ice-out in spring,  mixing  depths  were high,
i.e. temperature was  uniformly distributed  throughout  the entire lake.  That
was also the  justification for selecting April as the  initial time for numerical
simulations.

     In  summer  mixed  layers  were deeper in two  of the  three lakes under
warmer  (1988) conditions.    Increased  net heat  flux  to the lake caused  a
slightly earlier onset of stratification.  The simulated onset  of stratification is
first  observed  in  the smallest  lake (Holland Lake).  Lake Calhoun and Lake
Elmo  started  to  stratify  later  and showed similar  mixing events on  a daily
timescale.    Vertical  mixing  is   caused  bv  wind  and  natural   convection.
Surface mixed layers were  deeper in Lake Elmo, mainly because  more wind
energy  was available for entrainment at the thennocline due to the longer
fetch (greater  surface area) of the lake.   Natural convection is mainly driven
by net surface  (evaporative, conductive) heat loss.   Under  warmer  conditions
evaporative loss was  much higher, ana  1 to 2  m  deeper mixed  layers were
probably produced in this way.

     Fall  overturn occurred  earlier after  the  wanner  1988 summer  because
lower  fall   temperatures  produced   stronger  cooling  and   surface  water
instabilities,  i.e. thermals and  convective  negatively buoyant  (cold) currents
earlier  (Horsch  et  al.,1988).  In the presence of convective cooling, less
turbulent kinetic  energy,  supplied by the wind, is needed for the deepening of
the thermodine.

4.4.4      Water temperatures

     Daily epilimnetic  temperatures at a  depth  of 1.5 m are shown  in  Fig.
4.7.    Although   lakes   have   different  morphometries,  similar temperature
patterns  were  observed.   This is  in agreement with  field measurements made
by Ford and  Stefan  (1980) in 1974  and  1975.   In  both 1971 and 1988  the
surface  temperatures of  the  three  lakes  exhibited  similarities  and  parallel
trends  which  are  predominantly related  to weather  phenomena  and  only
secondarily  to  lake  morphometry  (Ford  and Stefan,  1980).   From  April
through  August epilimnion water temperatures were higher in 1988 (average
water  temperature  increase    3C  compared  to   3'C  in  air temperature
change)  and  lower after  the lake  staried cooling.

     Daily hypolimnetic  temperatures  are  shown in  Fig.  4.8.   Values  are at
depths  well  below   the  thennocline,  and  water   temperatures   are  nearly
isothermal below  that  depth.   Lake Calhoun and  Lake  Elmo received
additional  heat during the spring turnover  periods (Fig. 3.6) when the climate
was  warmer  (1988).   Average hypolimnion temperature  was  0.6  and 1.4* C
higher  in  1988 in Lake  Calhoun and  Lake Elmo, respectively, than in 1971.
Lake  Holland experienced  an opposite trend.   The lake  started  to  stratify
earlier  too, but  due to  the  increased  stability and small  lake  surface area,
wind mixing  throughout  the  entire lake in spring  under warmer  conditions
                                      69

-------
was  weaker.    Average  hypolimnion  temperature  was  1.2C  lower  under
warmer  conditions.    Once  a  stable  stratification   was  established,   the
hypolimnetic temperature was almost constant throughout the summer  for all
three lakes.                                          -

     As is  typical  for  dimictic  lakes  in  temperate  regions,  the  summer
temperature  in  the hypolimnion was determined  by mixing  events in  spring
and  remained almost  constant  throughout the rest of the simulation  period.
Lake  Elmo,  although  twice  as   deep  as  Lake  Holland,  had  a  higher
hypolimnetic temperature.   Point inflows in these lakes were not significant,
and  the  hypolimnetic  temperature  difference  is  therefore   related   to  the
differences  in spring mixing dynamics, which through wind fetch,  is related to
the  surface areas  of the  lakes.  Greater wind shear stresses and hence wind
energy  inputs  are usually  associated with  larger  lake  surface  area  (longer
fetch).


4.5  Conclusions

     A validated  onedimensional  and unsteady lake water quality  model  can
be used to  study the changes in  a lake  as a result  of different  weather
conditions  including global  warming.   The analysis  described herein  is  a first
step  in quantifying  potential  thermal changes in  inland lakes due to climate
change.  Water temperatures in three lakes  in a sensitive latitude have been
simulated with weather from  two very  different summers.  Mean lake  depths
were 4.8, 10, and  13.4 m.

     The main  findings are:

     (1) Simulated  epilimnetic water  temperatures  responded  strongly  to
atmospheric changes.
     (2) Simulated hypolimnetic  temperatures  responded less   strongly   and
inconsistently (plus or minus) to atmospheric changes.   They were determined
by mixing  events  in spring, and lake morphometries.
     (3) Simulated evaporative heat losses increased about 40 percent  in  the
wanner summer.  Evaporative water  losses increased by about 300 mm out of
800  mm or about  40 percent.
     (4) Dates of ice  formation in fall  seemed  only  weakly  affected  by  the
hot  midsummer weather.   Dates shifted by  a few days.  This  may  not be
typical  because of the cool  fall.
     (5) Simulated conductive heat  transfer had  a negligible  effect  on  heat
budget  changes.
     (6) Higher  atmospheric radiation  due  to  higher  air  temperature  was
compensated by higher  backradiation  from the water.
     (7) Simulated surface mixed- layer  depths increased slightly (by  1  to  2
m) in  the  warmer summer, probably  due to  stronger convective mixing.
     (8) Simulated stratification onset  occurred slightly earlier in the warmer
year.
                                     70

-------
                SIMULATION PERIOD 4/01-10/31
  45

  40-

  35-j

  30-j

I 25-|
r
j 20-


  15~

    3
   5-;
    ]
   o4-
    j
                     LAKE CALHOUN
     1971
     1988
  40-
  35 -i
& 20-(


  15-
                     LAKE ELMO
                                             1971
                     tAKE HOLLAND
  40-
	 1971
	 196=
  jo-:
2 2*-3
& 20-J
         ;PR   K/.-V    JUN  . JUL   AUG   SEP   CCT
      Fig. 4.6     Mixed  layer depths (simulated).
                               71

-------
          SIMULATION PERIOPD 4/01-10/31








g

UJ
*-^i*' A
' ^"' A/V\

j^\ MI/ ^/V\/\AA A . 
v/-  \ :
 , i  i  i  i i i  t  i
	 " 1971 LAKE ELMO
--- 1 QOQ
(j7OO -
 DIFFERENCE
i* /*-* '.' ' *"* * -
,/'*.* ..;*"* V'"\ _....--"*t
/ * """"" * " \ 
'X ;"v' ~"*'\
- ' ./ "l<1*..
/,>. ;*'. ; '' V-^. ..,
* ." V '""" *


^ . f "**" A 

"' """" / u^ yA ~
Pv Afv y\/^ \. .A
W" ^ \ -
-
  i  i  i  t  t i
~- 	 1971 LAKE HOLLAND "
	 1988
	 DIFFERENCE  - , /-, >. __ ,,
/'^^V'-'-'V' '"">.'.
A-; /U ^ ''\"/'\-
** ' '  ' *** tir*
,;/l/^ /^ : 'v.'
     APR    WAY   JUN   JUL   AUG   SEP
Fig.  4.7     Simulated epilimnion temperatures.
                      72

-------
                SIMULATION PERIOPO 4/01-10/3!
  12-
 *-
2
     ....._._ 1971
     	 1988
     	 DIFFERENCE (1988-1971)
                                   LAKE CALHOUN \
         i  _...
  -2-
	- 1971
	 1988
	 DIFFERENCE
   2H
                                        LAKE ELMO
  -2'
          1971
          1988
          DIFFERENCE
   6-4
                                   LAKE HOLLAND \

  -2
         APR    MAY   JUN   JUL   AUG    SEP   OCT
    Fig. 4.8     Simulated hypolimnion  temperatures.
                                73

-------
     5.  Water Temperature Characteristics of Minnesota

               Lakes Subjected  to Climate Change


     A deterministic,  validated, one-dimensional,  unsteady-state  lake water
quality  model  was linked to  a daily  weather  data base to  simulate  dailv
water  temperature  profiles in  lakes  over  a period of twenty-five  (1955-79)
years.   27 classes of lakes which are characteristic for the north-central  US
were investigated.   Output from a global  climate  model (GISS) was used  to
modify the weather data  base  to account for a doubling of atmospheric  CO^-
The  simulations  predict  that  after  climate change  epilimnetic temperatures
will  be  higher   but   increase  less  than  air  temperature,  hypolimnetic
temperatures  in seasonally stratified  dimictic lakes will be largely  unchanged
or even lower than/at present, evaporative  water loss will  be  increased by  as
much as 300  mm  for the season, onset of stratification will occur earlier and
overturn later in  the  season, and overall lake stability will become greater  in
spring  and summer.


5.1  Introduction        v

     This Chapter deals with the question  of how climate change may affect
thermal aquatic habitat  in  lakes.  A regional  perspective is  taken, and the
scope is  to estimate temperature  changes  in lakes of different  morphometric
and  trophic characteristics in a region.  Southern Minnesota is chosen as  an
example  because  an  extensive  lake  database is  available  (ERLD/MNDNR,
1990).    The  geographic  boundaries  of Southern  Minnesota  are  defined  in
Figure  5.1.

     Herein a  dynamic and  validated regional lake water  temperature model
(Chapter 2) will  be applied to a representative  range of lakes in a  region for
past  climate  and  one  future climate  scenario.    Rather  than  analyzing
particular years and  lakes,  emphasis is on long  term behavior and  a wide
range  of  lake  morphometries  and trophic  levels.   In  this study  the  base
period  (or comparable  reference) was  from  1955  -  1979.  For the same period
of time  weather  parameters  were  perturbed by  the  2XC02  GISS  (Goddard
Institute for  Space Studies)  climate  model  output.   The  regional  impact  of
these climates  on  different  lake  classes  in  southern  Minnesota  is reported
herein.  The  simulated water temperatures,  past  and  future, will be  presented,
interpreted and related to the  lake characteristics  and climate characteristics.
The  results will  show  how  water  temperatures  in different freshwater  lakes
respond to changed atmospheric conditions  in a region.

     Lake levels  will  be largely controlled  by  the  water budget including
evaporation and runoff.  The response of watershed (surface) runoff to climate
change is the subject of other  investigations not included herein. Lake depths
                                    74

-------
Fig.  5.1      Regional  boundaries and-geographic distribution of lakes in MLFD
             database.
                                   75

-------
will  therefore be treated herein as either  invariant or be lowered  to account
for increased evaporative water  losses, where applicable.   Changes  in  the
watershed  may  affect nutrient loadings and  hence primary  productivity  and
transparency of  the water.   Such secondary effects, also were not investigated,
but a  sensitivity analysis indicates that water temperature predictions  for the
types  of lakes   studied  herein  are   usually  not  sensitive  to  transparency
(Chapter 2).


5.2  Method of Lake Temperature Modeling

     The numerical model  is applied in  daily  timesteps  using  mean daily
values  for the meteorological variables.  The  required  weather  parameters are
solar  radiation,  air  temperature,  dew point  temperature, wind speed, wind
direction,   and    precipitation.      Initial  conditions,   lake   morphometry
(area-depth-volume), and Secchi depth have to be provided  to  use the model.
Simulations  were made from spring overturn to fall overturn.   Since the date
of spring overturn  is  unknown,  the  initial  conditions were- set  at  4"C  on
March 1, and  water, temperature  was not allowed to drop  below  4 C (well
mixed  conditions).   Although isothermal water at 4'C may  not exactly  exist
on March 1, the isothermal 4" condition  continues until  the model  simulates
warmer temperatures and the onset of stratification.   The summer predictions
are  thus    made   quasi-independent   of  initial   conditions   and  match
measurements well  (Chapter  3).  The model is one-dimensional in  depth  and
unsteady, i.e.  it  simulates   water temperature   distributions over   depth  in
response  to  time  variable  weather.   Vertical water  temperature simulations
are made over  an entire season (March 1  to November 30)  and in time steps
of one  day.  The  calculated  daily  water temperature profiles  are analyzed
statistically  and presented graphically.

     The regional  water temperature  simulation  model was  validated  against
data from  nine Minnesota lakes  for   several  years  (Chapter 2).  The model
simulates onset  of  stratification,  mixed layer depth,  and water  temperature
well.   Root mean  square error is 1.2C,  and 93%  of measured  lake water
temperatures variability  is explained  by the numerical simulations,  over  wide
range of lake morphometries and trophic levels.


5.3  Climate Conditions  Simulated

     Meteorological data from the Minneapolis-St.  Paul  International  Airport
(93.13*  longitude,  44.53* latitude) were used.   The meteorological data file
assembled  contains  measured daily values of  average air  temperature, de
point  temperature,  precipitation, wind s,peed, and  solar radiation from 1955 u
1979 (March -  November).   The  period from  1955 to 79 was  chosen becauw
it  is long enough  to give  a representative average  of base  conditions befcrr
climate warming.    In the  1980s  warmer  than  average air temperatures
observed  (Jones et  al.,  1986; Kerr, 1989),  and  therefore  this  period
excluded.  Sources  of climate data were as follows:
                                     76

-------
     Climate  scenarios  were  selected  following  EPA  guidelines  on  global
dimate  change  effect studies  (Robinson  and  Finkelstein,  1990).    Climate
projections  by  four different  models  (GISS, GFDL,  OSU,  UMKO)  for the
doubling  of atmospheric  C02 were provided  by NOAA (1990).   The  monthly
dimate projections  by the four models  are different from each other and their
explicit effects on  water  temperature dynamics can be studied for each  model
separately.   In this study only the GISS projections for  the  grid point closest
to Minneapolis/St. Paul were  used  (Table 5.1),  as  suggested  by EPA for
effect studies.   The geographical  location of this grid point is given in Figure
5.2.   A  comparison  of the  mean monthly weather parameter values  (for
Minneapolis/St.  Paul) projected  by  the four models shows that  the  GISS
projections  are not substantially different from  GFDL and  OSU, except for
wind speeds in November.   No adjustments  were made to those wind speeds,
however,  for a lack of a rational basis and because  late fall  winds do not
affect the  summer water temperature  dynamics.   No interpolations  between
grid  points were made, following  explicit EPA recommendations.


Table 5.1  Weather parameters changes  projected by the 2XCO2
               climate model output for Minneapolis/St.  Paul.
MONTH
JAN
FEB
MAR
APR
MAY
JUN
JUL
AUG
SEP
OCT
NOV
DEC
AIR. TEMP
(diff.'Q*
6.20
5.50
5.20
5.05
2.63
3.71
2.15
3.79
7.02
3.73
6.14
5.85
SOL. RAD.
(Ratio)t
0.92
1.04
0.98
1.03
1.00
0.99
0.98
1.04
1.04
1.12
1.03
0.99
WIND S.
(Ratio)l
0.92
1.12
0.47
0.69
0.67
0.85
0.93
1.00
1.07
2.23
5.00
0.77
REAL. HUM.
(Ratio) t
1.16
1.01
1.13
1.00
1.09
1.01
0.93
1.02
0.90
0.95
1.00
0.98
PRECIP.
(Ratio)t
1.17
1.03
1.28
1.03
1.12
1.08
1.10
0.98
0.70
0.88
0.99
1.24
* Difference  = 2XCOz GISS - PAST
t Ratio = 2XCO2 GISS/PAST
     The uncertainty  of the climate predictions  is  not  the  subject of  this
paper.   It  is understood that relative humidity and wind speeds are not  well
predicted at the local scale  by global climate models.  Fortunately, .
uncertainty  analysis of the effects of variable  meteorological forcing  on  lake
temperature models indicates that  air  temperature has the  most significant
effect in lake temperature  uncertainty  (Henderson-Sellers,  1988;  Chapter  3),
and  that parameter is  better predicted than others.

     Seasonal  distributions  of  the 25-year  average  of  observed  weather
parameters (which were used as model inputs) are shown in Figure 5.3.   Past
climate and  the 2XCOa GISS scenario were  used  as  inputs to the  water
temperature simulations.


                                    77

-------
oo
                                 MINNEAPOLIS/SLPAUL
                                9GISS GRID POINT MINNEAPOLIS/St.PAUL
                                 DULUTH                            >
                                BBGISS GRID POINT DULUTH
                                      north dokoto   o
                                                 o
                                                 V)
                                                 0)

                                                 c  -  wscons
                                    km
                           Fig. 5.2     Geographical location of the closest GISS grid points for
                                      Minneapolis/St.  Paul and Duluth.

-------
, I
                                                                                                       -
                                                                                      MAR  APR  MAY  JUN   JUL  AUC   SEP  OCT  NOV
                                                         	 MSTP (1955-1979)

                                                         	 2XC02 CLIMATE SCENARIO (CISS - model)

                                                                                      MAR  APR  MAY   JUN   JUL  AUC   SEP   OCT   NOV
                                                                                                                                   30
                                                                                                                                       o

                                                                                                                                    20  8
-15




 10




5
                                                                                                      ui
                                                                                                      ui
                                                                                                      a.
                                                                                                                                       a
                                 Fig. 5.3        Climate  parameters  a  Minneapolis/St.  Paul  in  the past  and under

                                                  a 2xC02 (GISS) climate  scenario.

-------
5.4  Regional  Lake Characteristics

     Regional  classification of  lakes  was  approached in  a variety  of ways.
The ecoregion  approach was  considered first, but found to  give too  detailed a
picture.  The entire  state was considered  as a regional entity but rejected as
too  large because of  the  diversity  of  climate.    Dividing the  state  into  a
northern and southern region  was considered appropriate  and not as arbitrary
as  might   seem   because   there  is  a   significant  gradient  in   geological,
topographic,  hydrological dimatolpgical  and ecological  parameters  across the
mid-section  of  the state  (Baker et  al.,   1985, Heiskary  et al.,  1987).   The
southern and northern  regions are about equal in 'size (Fig. 5.1).

     The  Minnesota  Lakes  Fisheries  Database,   MFLD  (ERLD/MNDNR,
1990),  which contains  lake  survey data for  3002  Minnesota  lakes,  is  for the
southern region.   The  MLFD database includes 22  physical variables and fish
species.  Nine  primary variables explain 80 percent of the  variability between
lakes.   These  nine variables  include  surface  area,  volume, maximum  depth,
alkalirJty, secchi depth, lake shape, shoreline complexity,  percent  littoral  area,
and length of growing season.   For regional classification  of  the  lakes  in this
study,  the possible thermal  structure  (i.e.  whether  lakes are stratified or not)
and trophic  status are of  primary  concern.   Observations  in  the northern
hemisphere  show  that  onset and maintenance  of stratification in lakes  is
dependent on surface area and maximum  depth (Gorham  and Boyce, 1989) as
well &s climatological forcing i.e. solar radiation  and wind (Ford and  Stefan,
1980).    Lake  trophic  status  contributes  to solar  radiation  attenuation  and
oxygea  balance.   Trophic status  was assessed  by  using a Secchi depth  scale
(Heisiary and  Wilson,  1988)   related  to  Carlson's  Trophic  State   Index
(CarUon, 1977).  Secchi depth information  was  available in the  MLFD.

     A statistical  analysis  of southern and northern Minnesota  Lakes  in the
MLFD  in terms of surface  area,  maximum depth and Secchi depth  was made.
The ceographic distribution  of different classes of lakes in  Minnesota is  given
in  Fipre 5.4.   Cumulative frequency distributions  shown  in Figure 5.5  were
used **o subdivide all lakes into  three ranges of surface area, maximum depth
and fecchi depth, as shown in Table  5.2.   These represent 27 classes of  lakes
in e^'h of the  two regions  of the state.         "

                         Table 5.2 Lake classification
     Key                        -Cumulative                    Description
 Pa;anter            Range        Frequency           Class      Value


     >m2)              < 0.4     Lower 30%          0.2        Small
                   0.4-5       Central ^60%         1.7        Medium
                       > 5       Upper 10%          10        Large


         Depth         < 5       Lower 30%           4        Shallow
   )                  5-20      Central  60%         13        Medium
                       > 20      Upper 10%          24        Deep


  cci- Depth            < 1.8     Lower 20-50%       1.2        Eutrophic
  (z                1.8 - 4.5      Central  20-50%      2.5        Mesotrophic
                       >4.5      Upper 0-10%        4.5        Oligotrophic
                                      80

-------
      J\
                                                      I  I I  I  I  I  I  I  I  i

                                                        SECCM DEPTH > 4i METERS
                                                        AREA > 5.0 SO. KU.     cinLn
Fig.  5.4      Geographic  distribution  of  lakes   according  to key  parameters:
             Secchi depth,  maximum depth,  and surface area.
                                    81

-------
                  100
                    10-*
                             10"
                                       10V
                                   AREA (km2)
                                                10'
                  100
                                MAXIMUM DEPTH (m)
                  100
                            2468
                                SECCHI DEPTH (m)
                                                          10
Fig.  5.5      Cumulative  distributions  (%)  of  key  parameters  in  Minnesc1.;
             lakes (from MLFD database).
                                      82

-------
     A  representative value  for  surface  area,  maximum depth  and Secchi
depth was chosen in each lake class as input to the model simulation.  Those
values  are shown  under the heading "class" in Table  5.2.

     Representative  area-depth relationships  for  three different  lake  classes
(by surface area)  were obtained from  35 lakes which covered the  entire range
of distributions in a set of 122 lakes (Figure 5.6).

     After  areas  are  expressed as  fractions of  surface area  and  depths are
expressed as fractions of maximum depth,  an equation of the form

                       Area = a    exp(b- Depth) -f  c                    (5.1)

is  fitted  to the  data  and  subsequently  used  in  the  simulation  as  a
representative area-depth  relationship.    Coefficients  a,  b; c,  calculated  by
regression  analysis are given  in Table  5.3.   This  procedure is equivalent to
self-similarity of depth-area relationships within  a given class.
       .
Table  5.3  Morphometric  regression  coefficients
in the area vs.  depth relationship.
Area
Small
Medium
Large
1.19
1.14
1.14
-1.76
-2.10
-2.91
-0.20
-0.15
-0.08
      Lake  basin  shape was assumed circular  for  the purpose  of  wind  fetch
calculation.   The water  temperature  simulation  results  were  shown to  be
insensitive  to  these  assumptions  of  morphometric  self-similarity  and  basin
shape.


5.5   Simulated lake water temperature regimes for
      historical and future weather

5.5.1     Water temperatures

      Simulations   of  daily  water  temperature  profiles  from  March   1  to
November  30  (275 days)  in each  year  from  1955  to  1979  (25  years)  were
made  for  each of the 27  lake classes.    In addition  to  lake morphometric
input, i.e.  surface  area, maximum depth  and  depth-area  relationship,  these
simulations used  actually  recorded daily  values  of weather  parameters, i.e.
solar  radiation,  air  temperature,  dew  point  temperature,  wind  speed,  and
precipitation for each day  simulated.   A massive weather-database had to be
developed  prior to the  simulations.   The calculated output  of 185,625 vertical
water temperature profiles,  each consisting  of  24  water  temperature values,
provided  base line information  on lake characteristics  during  a period of the
past  when  little climate change occurred.
                                      83

-------
                            SMALL  LAKES
                                       G-O Corver
                                       - Crystol
                                           Deep
                                       G-Q Forquhor
                                       E-ffl Islond
                                           McCorrons
                                           McDonough
                                       i-A Long-N
                                           Lucy
                                       O-O Porkers
                                       * Reitz
                                           Sucker
                                           Twin-Middle
                                           Wolsfeld
                                                    Aubum Eost
                                                C (D Auburn West
                                                    Bold Eagle
                                                Or-EJ Boss
                                                -83 Bovorio
                                                    Big Cornelian
                                                    George
                                                    Colhoun
                                                    Crystol
                                                    Elmo
                                                    Harriet
                                                    Eagle
                            LARGE  LAKES
                                       3-O Woconia
                                       Z O Big Marine  ~
                                           Coon
                                       3-Q Whits Bear  ~
                                       9-3 Greenwood  :
               0.1   0.2  0.3  0.4.   0.5   0,6   0.7  0.8   0.9   1.0
       0.0
0.0
Fig. 5.6
      Horizontal  area vs. depth  relationship for lakes.
      and  depth  are normalized.
Area
                             84

-------
     To simulate  possible future  water  temperature  regimes,  the  monthlv
corrections  specified by the 2XC02 GISS  model  scenario were applied to the
weather data base and the simulations  were repeated.

     From these simulated water temperature data bases under historical and
future  climates,  each  consisting of  4,455,000 water  temperature  values, the
following characteristics were extracted.

     Epilimnetic water temperatures  were defined as  water temperatures at
1.0 m  below the water surface  regardless whether maximum-depth is  4  m, 13
m or  24 m, respectively.  The seasonal  course of epilimnetic  temperatures,
averaged weekly  over  25  years  is  shown in  Figure 5.7  for both past climate
and  the 2XC02  GISS climate  scenario.   The difference between  the two  is
also  shown  in  Figure  5.7;  the  associated  air temperature  increments due to
climate change were presented  in  Table 5.1.  The largest  change in weekly
water  temperature  change in response to climate change, is on the order of 6
to 7"C, and occurs in spring (April),  the  minimum  is  on the order of 0 to
2C  and occurs  either  in  fall (October  and November),  or  in July.

     The GISS  scenario  gives  a  seasonal  surface  water temperature pattern
different from that for the past.   The cooling phase, for example, commences
later and has stronger water temperature gradients.  Maximum  weekly surface
water  temperatures and  the time  of their  occurrence are  given in Table 5.4.
The  highest  surface water temperatures, 27.4* C  ( 0.1 C) were calculated  for
the shallow  lakes and  the lowest, 26.2 C (  0.1 C) for the deep lakes.   With
climate change  the predicted rise in the  seasonal surface water  temperature
maxima is  1.9 to 2.2C,  which  is small compared to air temperature changes
in Table 5.1.  The occurrence of the maximum water surface  temperatures  is
shifted by 11 to  20 days  towards the fall with the climate change.

     Surface water temperatures are fairly independent of lake  morphometry
within the range of lakes investigated.   Extreme values in lakes  of  different
geometry  vary   by no   more  than  4*C  on any  given day.    Maximum
differentials   occur  in  spring  and  fall.   From June through  September, i.e.
during  the  period  of  seasonal  water  temperature stratification, surface  water
temperatures in  lakes  of  different  morphometic  characteristics   (depth  and
area)  are very  similar (within  1.0 C).   In  very large  lakes (e.g. the  North
American Great  Lakes)  the  significantly  greater  water volumes  and  mixed
layer depths cause a substantial lag in heating  and  cooling leading to water
temperature  differences larger than 4C.

     Weekly averages of 25 years of simulated  hypolimnetic temperatures are
shown in Figure 5.8.   Values   are  1, m  above  the lake  bottom  (maximum
depth).  Hypolimnetic temperature  responses to  climate  change  show  wider
variability  than  epilimnetic responses.    In  shallow  (polymictic)  lakes, the
hypolimnetic  and  epilimnetic  water  temperature  rise  is very   similar  in
magnitude   and  time  of  occurrence.    In   deep small  lakes  hypolimnetic
temperatures are as much as 3.5C colder after  climate  change  than before.
Hypolimnetic wanning  during the summer is dependent on vertical turbulent
diffusion and therefore wind fetch and  hence surface area.   Dependence of
hypolimnetic temperatures on lake morphometry is very evident in Figure 5.8.
The seasonal  pattern  of hypolimnetic water temperatures  was  altered  by
climate change most  significantly in shallow  lakes.  All others showed typical
seasonal wanning patterns in response to vertical diffusion.


                                     85

-------
00
0>
     35



     30



     25



I     tO



     15







      i



      0



     30





   uj
   a: 20-
                                            CPIIIMHIOH ItMPERAIURCS

                                        SIUULAI10H P'HIOU 1951 - 19''J (I'ASI)
                                                              EPIUUNION TCUPCRATUttCS

                                                         }XC02 CUUMC JCCKWIO (CISS - model)
  CPICIMNION TtMPCRAIUNCS


DIFFERENCt (CISS modl - PST)
                                   SUAll LAKES
                                   MEDIUM LAKES
                                   (XBCC LAltrS
                                                                             SMALL LAKES
                                                                             MEDIUM LAKES
                                                                            -r-
                                                                             LARCC IAKCS
                                                                                                                                     I	I
                                                                                                                                                   i '    i    I '

                                                                                                                                                  SMALL LAKES
                                                                                                                                                   i '  '  I ' '  1
                                                                                                                                                  MEDIUM LAKES
                                                                                                                                                   LARGE LAKES
                                                                                                                                                                -1
                                                                                                                                                                4
                                 MAR APR MAY  JUN JUL  AUG SEP OCT NOV  MAR APR MAY JUH JUL AUG SEP OCT NOV MAR APR MAY JUN  JUL AUG SEP  OCT NOV
                                                                                                                                                                 -I
                                                        Fig.  5.7        Simulated weekly  epilimnion  temperatures.

-------
                                               HVPOL/MNtON
oo
35

30


25


10

IS


10


 5

 0

30




20-


15-

ID-


 S'

 0


30




20


15


10

 5


 0
                                                       1955 - 1979 (PASf)
                                                        MtPOtlUNION IEWPCRAIURCS
                                                  }XC02 CUUA1C SCCMKBIO (CISS - inodel)
  KTPOUVINION TEUPCPATURCS

OVTERENCE (CISS rmxltl - PAST)
                                     SUAU. UKCS
                                     UCOIUU UKES
                                     URGE LAKCS
                                                                               SKUIL LAXCS
                                                                                                                 \
                                                                               MEDIUM LAKES
                                                                               URGE IAXES
                                    *&
                                                   fj-v neiijM wigciropnie
                                          	- thaltoo Mlrophte	

                                    MAR  APR MAY JUN JUL  AUG SEP OCT NOV MAR APR MAY  JUN JUL  AUG SEP OCT  NOV MAR APR MAY JUN  JUL AUG SEP  OCT NOV
                                                       Fig.  5.8        Simulated  weekly  hypolimnion  temperatures.

-------
     Table  5.4  Maximum temperatures of southern Minnesota lakes
OO
oo
PAST 1955-1979
Maximum
Depth
m
SHALLOW
(4.0)







MEDIUM
(13.0)







DEEP
(24.0)







day = Julian
Area

km2
SMALL
(0-2)

MEDIUM
(1.70)

LARGE
(10.0)

SMALL
(0.2)

MEDIUM
(1.70)

LARGE
(10.0)

SMALL
(0.2)

MEDIUM
(1.70)

LARGE
(10.0)

day when
Trophic
Level

eutrophic
mesotrophic
oligotrophc
eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophc
> eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophc
eutrophic
meeotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
Epilimnion

ec
27.5
27.4
27.3
27.4
27.4
27.3
27.4
27.4
27.3
20.0
26.5
26.6
26.4
26.4
20.5
26.5
26.5
26.6
26.4
26.3
26.1
26.2
26.2
26.1
26.1
26.1
26.1

day
203
203
203
203
203
203
203
203
203
203
206
203
204
207
207
203
206
207
206
204
207
206
206
206
206
206
207
Hypolimnion

C day
24.9 206
26.8 204
27.0 203
26.2 204
27.0 203
27.1 203
26.5 203
26.9 203
27.1 203
11.9 278
12.8 277
17.5 261
18.7 254
19.9 252
23.0 233
24.0 220
24.6 218
25.5 211
7.3 308
7.4 308
7.8 308
11.6 294
11.8 293
12.6 291
18.2 261
18.4 263
19.4 259
GISS-2XCO2
Epilimnion

C
29.4
29.4
29.3
29.4
29.5
29.4
29.5
29.6
29.5
28.7
28.7
28.7
28.5
28.6
28.7
28.6
28.7
28.7
28.5
28.3
28.10
28.2
28.1
28.1
28.1
28.1
28.2

day
217
217
217
217
217
217
217
217
217
217
218
218
218
218
218
223
218
218
217
218
220
218
223
223
223
223
218
Hypolimnion

c
26.2
28.3
29.2
27.5
\ 29.1
29.4
28.3
29.1
29.4
12.6
13.0
17.5
18.2
20.3
25.3
26.0
26.6
27.3
10.3
10.4
10.6
12.8
12.9
13.3
18.4
18.7
20.1

day
229
218
217
205
181
217
181
181
217
289
284
276
274
271
248
233
224
218
305
305
305
291
291
291
276
276
273
DIFFERENCE
Epilimnion

c
1.9
2.0
2.0
2.0
2.1
2.1
2.1
2.2
2.2
2.1
2.2
2.1
2.1
2.2
. 2.2
2.1
2.2
2.1
2.1
2.0
2.0
2.0
1.9
2.0
2.0
2.0
2.1
(GISS-PAST)
Hypolimnion

c
1.3
1.5
2.2
1.3
2.1
2.3
1.8
2.2
2.3
0.7
0.2
0.0
-0.5
0.4
2.3
2.0
2.0
1.8
3.0
3.0
2.8
1.2
1.1
0.7
0.2
0.3
0.7
maximum temperature occur

-------
     The highest hypolimnetic  water temperatures  (27.1'C) were  calculated
for  shallow  oligotrophic lakes which are typically  polymictic or well-mixed for
the entire simulation period.  The lowest  maximum hypolimnetic temperatures
(7.3 C) occurred in small  and deep  eutrophic  lakes.   Climate change  raised
by  0  to 3*C the maximum  hypolimnetic  water temperature or lowered it  by
as  much as 3.5"C,  depending  on the particular  stratification dynamics of a
lake.

     In addition to  long-term changes of water temperatures (Figures 5.7 and
5.8) variations from year to  year  are also of interest.   Unfortunately weather
parameters  for  the GISS  climate scenario  were only  given  as  long term
monthly  averages.   Therefore  variability  on  an  annual basis  could not  be
explored  for  the  GISS  scenario.    On  the  other  hand,  annual  weather
information  was available for the 1955-79  period, and  therefore could be used
to  give the range of simulated  daily water temperatures.    Bands  of  water
temperatures  within  the  95%  confidence  interval are  shown  in Figure 5.9.
The spread  is significant  and on the order  of   3 to 5 C  around the  mean.
This  range  is  about  twice  as  wide  as  that  due  to differences  in lake
morphometry  (Figures 5.7   and  5.8).     This  is  in  agreement  with  field
measurements by  Ford and  Stefan (1980) and has  some bearing on habitat.
Examples of water temperature structures  in  typical  lakes are given  in Figure
5.10.

5.5.2     Thermal energy fluxes

     The water temperatures discussed  above are, of course, the result of net
heat energy input or losses  through the water surface, and vertical
distributions of that heat within the lake.   For better understanding of the
water  temperatures, it is therefore  appropriate  to consider,  at  least, briefly
heat fluxes  and stratification dynamics.   Simulated net heat flux through the
water  surface is plotted  in  Figure  5.12  for past and  future  (GISS) climate
conditions.

     Five  heat  transfer  processes  are responsible  for  heat input  into  the
water:     short  wave solar  radiation,   long  wave  atmospheric   radiation,
conductive  heat transfer,  evaporation, and back radiation.   Short  wave  solar
radiation  and atmospheric  radiation increase  the water  temperature,  while
evaporation  and back  radiation cool  the  water.   Conductive heat transfer can
either  heat  or cool the water.  All five fluxes together  comprise  net  heat flux
at the water  surface.

     Individual  daily  heat  fluxes  vary  dramatically  with  weather  as   is
illustrated in  Figure 5.11.  To keep-track of the  extraordinary dynamics  and
to explain  them  would take  more space  than available  here, and may not be
particularly  fruitful.    As  a  summary,  the  cumulative net  heat fluxes  are
presented in  Figure  5.12  for  past  and future (GISS) climate.   The  difference
between  the two is  also shown in Figure  5.12.  Lakes with  large surface areas
will receive  more net heat  input  (up  to  30%)  than smaller  ones, and in
extremely small lakes the  difference  is even negative, meaning less  heat  will
be  transferred through the water  surface  and stored in the lake!  All net heat
fluxes  are per unit surface area of a lake,  not total values.
                                     89

-------
               EPIJJUNION TEMPERATURES EUTROPHIC LAKES
                     PAST 1955-1979
HYPOUUMON TEUPERATURES EUTROPHlC LAKES
      PAST IS55-1979
                                               depth dcp (24 m)
                                               oreo small (0.2 km
                                               depth deep (24 m)
                                               oreo targe (10 km*)
         MAR APS MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JUN' JUL AUG SEP OCT NOV

Fig. 5.9       Examples    giving    range    of    epilimnetic    and    hypolimne:
               temperatures  over  a 25 year  period  (95%  confidence interval).
                                          90

-------

to
                          Fig.S.lOa
Simulated temperature  (iiotherml structure in three  medium deep (13 m maximum depth) lakes  of
large (10 km), medium (1.7 km') and  small (0.2  km') surface area. Isotherm  bands  are in increments
of 2'C.  Simulated water  temperatures are for past  climate (1955-79)  (lop)  and the 2XCO}  GISS
climate icenario (bottom).

-------
Fig.S.lOb   Simulated  temperature  (isotherm)  structure  in  three  medium  area  (1.7  kmj)  lakes  of  different
           maximum depths:  shallow (4 m), medium (13 m) and deep (24  ro). Isotherm bands are in  increment
           of 2'C  Simulated water temperatures are for  put  climate  (1056 70)  (tup)  and the 2XCO] G1SS
           climate scenario (bottom).

-------

-------
CO
                       600000
                       400000
                                       CUUULA1IVC NO MEAT flXIX

                                    SIMULATION PRIOO 1955 - 1979 (PAST)
                                                   CUMULATIVE NET MEAT fUJX

                                               2XC02 CLIMATE SCENARIO (GISS - model)
  CUMULATIVE NET HEAT rwx
DlFFtREMCE (CISS model - PAST)
                                                                                                                                                          190000
                                                                                                                                                          100000
MAR  APR MAY JUN JUL  AUG SEP OCT  NO^MAR APR  MAY JUN  JUL AUG  SEP OCT NOV
                                                                                                                                             I	t ' t -100000
                                                                                                                MAR APR MAY JUN JUL AUG SEP OCT  NOV
                                                               Fig. 5.12      Simulated cumulative  net  heat flux.

-------
      Back  radiation and evaporation are the  main  processes by  which lakes
lose heat in  the  summer.   Evaporative  losses were  found to be  significantly
increased after climate  change  (GISS).    In  all lakes,  regardless  of depth,
surface area  and  trophic status, the  computed evaporation water losses  were
uniformly 0.30 m (   0.01  m) higher (Figure  5.13).   In  other  words, lake
water budgets will be put  under significant  stress.  This  increased evaporation
also  explains why  the  water  temperature  increases  after  climate  change
remains  at  a  relatively  moderate 2*C, when  air temperature increases by  an
average  seasonal  simulation  (March  1,-  November  30)  value  of  4.4  "C.
Evaporative  cooling  is  a key to   the  understanding  of  the  temperature
responses to changed climate.

-5.5.3     Vertical mixing/Stratification/Stability

      Vertical   mixing   and  stratification   affect   lake  water   temperature
dynamics.   A surface  mixed  layer  depth is defined  here as the  thickness of
the  isothermal  layer from  the water  surface downward.   Surface  mixed layer
depths are  calculated daily by  the wind  mixing algorithm in the model  and
averaged over  a  week  (Figure  5.14).   Mixed layer  depths  at  the beginning
and  before  the  end  of simulation  are  equal  to  the total lake  depth  and
indicate  spring and fall overturns.  The most  shallow mixed layer  depths were
calculated for small, deep, eutrophic lakes based on the classification in Table
5.2.   Vertical  mixing is  caused by  wind  and  natural  convection.    Surface
mixed layer depths were the  shallowest for  small lakes because of short fetch.
Smaller  wind stresses  and hence  wind  energy inputs axe  usually  associated
with smaller  lake surface area  (shorter  fetch).   In these lakes  the  smallest
amount  of  turbulent  kinetic  energy is  available  for  entrainment   of  the
thermocline.    Wind  energy   required  for  entrainment  of layers   at  the
thermocline is inversely proportional  to the  stability (defined  as  a  density
difference between  adjacent  layers)  of the water  column and  depth of the
mixed  layer.    The  lowest   hypolimnetic  temperatures   and   the   highest
temperature  (density)  gradients were calculated for  small,  deep, eutrophic
lakes.  That  was  the reason for  the smallest mixed layer depths calculated  for
these  lakes.

      For the same  morphometric  lake characteristics, oligotrophic  lakes  had
deeper surface mixed  layers than eutrophic  lakes because of higher penetration
depth of irradiance.

      Climate change will  impose  higher  positive net heat fluxes  at  the lake
surface earlier in the season  than  in the past.  That  causes an earlier onset
of   stratification.     This  is  in  agreement  with   a  conclusion  derived  by
Robertson (1989) from  field data for Lake  Mendota.   In  the period from the
onset  of stratification  until September, mixed layer  depths were  projected in
the average  1.2 m smaller than in the past.   From  the end  of September,
mixed layer  depths were deeper after climate  change,  mainly due to  stronger
natural  convection and  higher winds caused  by climate change.  In spring1 and
summer evaporative  losses were  also increased  by  climate change   but  no
significant persistent cooling occurred  because  of net  heat input  from radiation
and  convection.  The earlier onset  of stratification  in spring and the mixed
layer  depth increase in  fall were also found by  Schindler et al. (1990) in  his
analysis  of  observations  in  the  EL A.     In  the  EL A mixed layer depths
increased due to transparency  increase  and increased winds due  to   reduced
forest cover resulting from increased incidence  of forest fires.


                                      95

-------
        SIMULATION PERIOD 18SS - l7 (PAST)
                                          ZXCOj CUMATE SCEN*0 (CISS - nyxl.l)
   MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY J'JN JUl AUG SEP OCT NOV
Fig.  5.13      Simulated  cumulative  evaporative losses.
                                96

-------
               MIXED LAYER DEPTHS

         2XC02 CLIMATE SCENARIO (GISS - model)
    5-
   10-
   15-
   20-
    5-
  10-
ui  15-
Q
   20-
    5-
   10-
   15-
   20-
   25
                                MEDIUM LAKES
                                LARGE LAKES
      MAR APR MAY JUN JUL AUG SEP OCT NOV
          MIXED LAVER DEPTHS

   SIMULATION PERIOD 1955 - 1979 (PAST)
    MEDIUM LAKES
                                                                        " i     i 
         X-X medium 0ltgoIropHic\

                  i eutropnlc

         O~O shallow otlgotrophic

                  eu trophic
APR MAY JUN JUL AUG SE? OCT NOV
                Fig.  5.14      Simulated  weekly  mixed layer depth.
                                             97

-------
     The stabilizing  effect of the  density  stratification  and the  destabilizing
effect of the  wind  can be quantified  using a  Lake  number  (Imberger  and
Patterson, 1989):

                             gSt(l  -Zt/zm)
                       Ln  = 	m	                   (5.2)
                             Po U, AQ '   (1 -  zg/zm)


where g is  acceleration due  to gravity  (m s~2), zt  is height from the lake
bottom  to the center of the thermocline (m), zffl is  maximum lake depth (mj,
zg is the height  of the center of volume of lake, A0 is lake surface  area (m2),

po is hypolimnion  density  (kg  m"3),  St  is  the  stability  of  the  lake  (kg  m;

Hutchinson,  1957), u^  is  surface shear  velocity (m  s"1).   Estimates  for  the

different elements in  the  Lake  number  are obtained from  daily  lake  water
temperatures  simulations,  daily  meteorological  data,  and  lake  geometry.
Larger Lake number values indicate stronger stratification and  higher  stability
i.e. forces introduced  by the wind stress will have minor effect.   Lake  number
dependence on lake/area, depth, and  trophic status,  for different lake  classes
is  given in  Figure  5.15.   Stability is  higher  for oligotrophic  lakes  than
eutrophic  lakes.    Oligotrophic  lakes  had  deeper   thermoclines  and  required
greater  wind force in  order to overturn  the  density structure  of  the  water
column.   Climatic change caused  higher  lake numbers,  i.e.  more   stable
stratification among the same lake classes.

     Seasonal   stratification  is  defined   herein   as  the  condition   when
temperature  difference between  surface and deep water  temperature is  greater
than 1C.   Although  1*C  is an arbitrary criterion, it  is useful  to  identify a
possible stratification  shift  with climate change.   With the above definition,
stratification  characteristics for  southern  Minnesota  lakes are  given in  Table
5.5.   A seasonal stratification ratio (SSR)  is defined as  the total  number  of
days when stratification stronger than 1C  exists, divided by the period from
the  earliest  to  latest  date of stratification.    A   SSR  ratio less  than   1.0
indicates a  polymictic, typically  shallow  or   a medium-depth  large  lakes.
Other lake categories  were dimictic since  the seasonal  stratification ratio  was
1.0.   In other  words, once seasonal  stratification  was  established, it  lasted
until fall overturn.

     Climate  change  advanced the  onset  of  seasonal  stratification  in   the
average  by 50 days for shallow lakes, and 34 days  for deep  and  medium deep
lakes.   Length of stratification  was prolonged by 60  days for  shallow  and by
40 days for  deep and  medium  deep lakes.
                                     98

-------
                           i- .otrphic PAST<1955-1979)
                           CHD shollo. oSgotrophi* PAST (1955-197?   B:ALL LAKES
                               ihollow eutrophic 2XC02CtSS
                       80- X-K stxiuow digotropnic 2XCO2SS        ,,
                                           10      15      CD     25     30
Fig.  5.15     Simulated  lake
               status.
ce numbers as a  function of lake  depth and  trophic
                                               99

-------
    Table 5.5  Seasonal stratification  characteristics of southern  Minnesota lakes
o
o
LAKE CHARACTERISTICS
MAXIMUMSURFACE
DEPTH AREA
m
SHALLOW
(4.0)







MEDIUM
(13.0)







DEEP
'?! 0)







km2
SMALL
(0.2)

MEDIUM
(1.7)

LARGE
(10.0)

SMALL
(0.2)

MEDIUM
(1.7)

LARGE
(10.0)

SMALL
(0.2)

MriMUM
'

* ** **


TROPHIC
STATUS

Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
OliRotrophic
Kutrophic
"'i '. ;!: :<
' <
f .....-. ^ ,
*

BSS
day
118
134
0
132
0
0
133
0
0
100
100
101
105
106
106
106
106
124
101
101
101
104
HJ1
* ( i *
".A
*

PAST 1955-1979
ESS
day
269
241
0
244
0
0
240
0
0
293
290
268
262
256
241
241
233
210
312
313
312
295
295
'!'!'.'
?*
VI i

LSS
day
152
108
0
113
0
0
108
0
0
194
191
168
158
151
136
136
128
87
212
213
212
192
192
1KH
l')9
1 i
i i
SSR

0.89
0.12

0.63


0.19


1.00
1.00
1.00
1.00
1.00
0.99
0.86
0.87
0.99
1.00
1.00
1.00
1.00
1.00
1. 00
0.99
\ r ^^

MAXSD MINSD BSS
m
1.7
1.9

1.6


1.8


5.0
4.5
7.0
4.9
4.9
6.1
3.5
4.5
5.5
9.0
10.0
13.0
10.0
10.0
12.0
8.0
<) <

m
0.1
0.2

0.1


0.1


0.2
0.2
0.2
0.4
0.4
0.4
0.4
1.0
1.0
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0 .1

day
68
85
0
76
116
0
85
116
0
68
69
70
69
69
70
70
72
73
69
70
70
70
. 71
72
72
71
7?
GlSS-2xC02
ESS
day
271
246.
ox
255
138
0
255
137
0
288
287
276
274
270
251
250
250
247
301
302
302
290
290
290
275
275
273
LSS
day
204
162
0
160
23
0
171
22
0
221
219
207
206
202
182
181
179
175
233
233
233
221
220
219
204
205
202
SSR

0.98
0.54

0.84
0.22

0.54
0.14

1.00
1.00
1.00
1.00
1.00
1.00
0.99
0.97
0.90
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
GISS
MAXSD MINSD BSS
m
1.7
2.1

1.9
1.3

1.3
0.9

5.3
6.7
9.1
5.9
4.0
5.3
2.9
4.3
5.3
5.8
6.7
8.6
7.7
8.6
10.6
11.5
13.4
9.6
m
0.1
0.1

0.1
0.4

0.1
0.3

0.2
0.1
0.1
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
day
-50
-49

-56


-48


-32
-31
-31
-36
-37
-36
-36
-34
-51
-32
-31
-31
-34
-33
-33
-34
-35
-34
- PAST
ESS
day
2
5

11


15


-6
-3
8
12
14
10
9
17
37
-11
-11
-10
-6
-5
-2
11
11
13
LSS
day
52
54

67


63


27
28
39
48
51
46
45
51 
88
21
20
21
29
28
31
45
46
47
SSR

0.09
0.42

0.22


0.35


0.00
0.00
0.00
0.00
0.00
0.02
0.13
0.10
-0.10
0.00
0.00
0.00
0.00
0.00
0.00
0.01
0.00
0.00

-------
HSS      Beginning seasonal stratification, i.e. first Julian day  when  difference
         between surface and  deep water temperature is greater than 1"C.

".SS      End seasonal stratification,  i.e.  last  Julian  day  when  difference
         between surface and  deep temperature is less than  1C.

LSS      Length of seasonal stratification (ESS-BSS)+1(

SSR      Seasonal   stratification  ratio,  i.e.  total   number   of  days  when
         difference  between surface and deep  water temperature  is  greater
         than 1 C  divided by LSS

MAXSD  Maximum  stratification  depth, MINSD  - Minimum  stratification
         depth
                                    101

-------
5.5  Conclusions

     A regional simulation study  was  conducted  for 27  classes  of  lakes in
Minnesota.'  Lakes  were  classified according to area,  maximum  depth,  and
trophic level.   A  validated, one-dimensional,  unsteady  lake water  quality
model  was  linked  to  global  climate  model output  in  order   to  quantify
potential  thermal changes in inland  lakes'  due  to climate  change.    Water
temperatures were simulated  on  a daily  time  base for past weather conditions,
1955-1979  and the 2xC02 GISS  model climate scenario.

     The main findings are as follows:

     (1)  Simulated epilimnetic  temperatures were  predominantly  related to
weather  and  secondarily to  lake  morphometry.   Weekly  average  epilimnetic
temperatures  were  raised  by climate  change  for  all lake  classes.     The
seasonally  averaged  water  temperature rise was  3'C, compared to 4.4C  air
temperature increase caused by the climate change.  The largest differences in
water  temperatures  occurred  in  April and  September,  and  were  7.2 C  and
4.9 C,  respectively.  The  seasonal daily  maximum of epilimnetic temperatures
rose only  about 2/C with climate  change.
     (2)  Hypolimnetic  temperatures  were  predominantly  related  to   lake
morphometry  and mixing events in spring, and only secondarily to weather in
summer.    The highest  temperatures  were calculated for  large,  shallow,
eutrophic  lakes.  After  climate change,  hypolimnetic water temperatures  were
as follows:  shallow lakes,  warmer  by  an average 3.1'C; deep lakes, cooler by
an  average  1.1 C;  small-area,  medium  depth  lakes,  cooler by  1.7C;  and
large-area  medium-depth lakes, warmer by 2.0* C.
     (3) Simulated  evaporative  heat  and  water losses  increased by about 30
percent for  the 2xCOz GISS  climate  scenario.    Evaporative water  losses
increased by about 300  mm,  making the total water loss 1200 mm.
     (4) Net heat  flux  at  the  lake  surface  increased  with  changed climatic
conditions.   The largest difference in  calculated  cumulative  net heat storage
between  past  and future climate was  100,000 kcal  m*2 and  occurred  in  April
and  September with climate change.
     (5)  Simulated  mixed layer depths  decreased  about 1  m in  the  spring
and  summer, and  increased in the fall.
     (6) With climate change, lakes stratify earlier, and overturn  later  in thr
season.  Length of the stratification period was increased by  40 to  60 days.
     (7) Climate  change caused  greater lake stability in spring and  summer
In fall lakes were driven faster towards isothermal conditions.
                                     102

-------
                               6.  Summary


     JLs  a  result of the  research described here,  a  better  understanding  of
v. *  freshwater  inland  lakes  respond to variable atmospheric  conditions has
'XT'  trained.

      Chapter  2  describes  how a specific lake water  temperature model was
r-cer-j-Iised to simulate the seasonal  (spring to fall) temperature stratification
 TCT   a  wide  range  of  lake  morphometries   and   meteorological conditions.
Mcc=L coefficients related  to  hypohmnetic  eddy  diffusivity,  light attenuation,
ice  sheltering and  convective heat  transfer were generalized using theoretical
xr.d   empirical   model   extensions.      The  proposed  regional  lake,  water
ictnp-exature model simulates   the onset  of  stratification,  mixed  layer depth,
ir.d Crater temperatures well.

      Hypolimnetic  eddy   diffusivity  was estimated   as  a  function   of  lake
surface area and stability  frequency.   Although the  proposed relationship is  a
simplification  of  the  turbulent diffusion  processes   taking  place   in the
hypcKmnion, it  was  found to be useful in  seasonal  lake  water temperature
modeling.    Heat exchange between water  and lake  sediments,  a  process
commonly  neglected  in  previous work,  was found  to  be  important  for the
analysis of vertical hypolimnetic  eddy diffusivity  (Appendix A).   Estimates  of
hype-limnetic eddy diffusivity  made without  sedimentary heat flux were up  to
one third  smaller than those  made  with the  heat  flux.   Effects of errors  in
temperature measurements and sediment heat  flux estimates on  the estimated
vertical  eddy diffusivity were  evaluated as well.

      Chapter  3 describes  a first order  analysis of  uncertainty  propagation  in
lake temperature modeling.   The output uncertainty  is defined as the  result
of deviations  of  the  meteorological  variables  from  their mean  values.   The
analysis   was   applied    to   systems    with   correlated   and   uncorrelated
meteorological variables.    Surface water temperatures  are strongly affected by
uncertain  meteorological  forcing.   Air temperature and dew point temperature
fluctuations  have  a  significant effect on  lake   temperature  uncertainty.
Long-term average water  temperature structure in  lakes can be estimated by
computer  model  simulations  for just  one  year  when  the  results from the
statistical   analysis of  meteorological , variables  are  used as  input.     This
analysis presents  a useful  alternative for the study  of long-term averages and
the   variability  of  temperature   structures   in   lakes   due   to    variable
meteorological  forcing.     In   addition,   the  analysis   revealed  the  separate
contribution of each meteorological variable to water temperature uncertainty.

      The  analysis  described  in   Chapter 4  was a  first  step  in quantifying
potential thermal changes  in inland lakes due to  climate change.  Rather  than
using  global climate change predictions,  this analysis looked  at  the changes  in
heat  balance  and. temperature  profiles  in  a  particularly warm year  (1988)
                                     103

-------
compared to  a "normal" year  (1971).   A comparison  was  made for three
morphometrically  different  lakes  located  in north  central  US.    Simulated
water temperatures  were daily values driven by daily weather parameters  and
verified against several  sets of measurements.   The results show that in the
warmer year,  epilimnetic water temperatures were higher;  evaporative  water
loss increased;  and summer  stratification occurred  earlier in the season.

     Rather than analyzing particular years and  particular lakes, emphasis in
Chapter  5 is  on long term  behavior  and a wide  range of lake morphometries
and trophic levels.  The regional lake water temperature model was linked to
a daily meteorological data  base to simulate daily  water  temperature profiles
over a  period  of  twenty-five  (1955-1979)  years.   Twenty seven  classes of
lakes  which  are  characteristic  of the  north-central  US  were  investigated.
Output  from  a global climate model  (GISS) was  used to modify the  weather
data base to  account  for the  doubling  of  atmospheric  C02-   The  simulations
predict  that after climate change  epilimnetic temperatures  will be higher  but
increase  less  than  air  temperature;  hypolimnetic temperatures  in  seasonally
stratified  dimictic  lakes will  be  largely  unchanged or  even  lower  than  at
present; evaporative water loss will be increased  by as much as 300  mm for
the season,  onset  67 stratification  will occur earlier and. overturn later in the
season; and overall lake stability will  become greater in spring and summer.
                                     104

-------
                                References
 Adams,  E.E.,  Scott,  S.A.  and  Ho,  E.  K.  (1987).   "Vertical  diffusion in a
       stratified cooling lake,  Journal of Hydraulic  Engineering,  ASCE, Vol.
       113(3), pp. 293-307.

 Adams,  E.E.,  D.J.  Cosier, and K.R. Helfrich   (1990).   Evaporation from
       heated water  bodies: Predicting combined forced plus  free convection,
       Water Resources Research, Vol. 26(3), pp. 425^135.

 Aldama,  A.A,  D.R.F.  Harleman,  and  E.E.  Adams  (1989).    Hypolimnetic
       mixing in  a  weakly stratified lake,  Water Resources  Research, Vol.
       25(5), pp. 1014-1024.

 Andricevic,  R.  (1990).   Cost-effective network -design  for groundwater flow
       monitoring,  Stochastic Hydrology and  Hydraulics, Springer International,
       Berlin, Germany, 4(1), pp. 27-41.

 Baker,  D.G., Kuehriast, E.L.,  and  Zandlo,  J.A  (1985).  Climate of Minnesota,
       Part   15  - Normal  temperatures   (1951-1980)  and their  application,
       Agricultural Experimental Station University of Minnesota, ADSB-2777.

 Bannister,   T.T.  (1974).    Production   equations  in  terms  of  chlorophyll
       concentration,  quantum  yield, and upper  limit to production, Limnology
       and Oceanography, 19(1), pp. 1-12.

 Bannister,  T.T.  (1979).    Quantitative description  of  steady  state,  nutrient-
       saturated   algal  growth,   including   adaptation,   Limnology   and
       Oceanography,  Vol. 24(1), pp. 76-96.

Birge,  E.A.,  Juday, C. and March,  H.W. (1927).  The temperature of the
       Bottom  deposits of  Lake  Mendota;  A chapter in  the heat  exchange  of
       the lake," Transactions  of the  Wisconsin  Academy of Sciences,  Arts and
       Letters, Vol 23,  pp.  187-231.

Bloss,   S.,   and  Harleman,  D.R.F.  (1979).  "Effect  of  wind-mixing on the
       thennocline  formation in  lakes  and reservoirs,"  Tech.  Rep.  249,  Ralph
       M. Parsons  Lab. for Water Resour. and  Hydrodyn;, Dep. of Civ.  Eng.,
       Mass.  Inst. of  Technol.

Bolin,  B. and  B.R.  Doos  (1986).    The  Greenhouse  Effect,  Climate  Change
       and Ecosystems,  John Wiley  and Sons, New York, p. 541.

Canale,  R. " P.,  and  S.  W.   Effler (1989).  Stochastic  phosphorus model  for
       Onondaga Lake,  Water  Research., 23(8),  1009-1016.
                                    105

-------
 Carlson,  R.E.(1977).   'A  trophic  state  index  for   lakes',   Limnology  and
        Oceanography, Vol.22, pp. 361-369.

 Carslaw,  H.S.  and  Jaeger, J.C. (1959).  Conduction  of heat in solids,  Oxford
        University Press, N.Y.

 Chang, L.H.,  S.F.  Railsback,  and R.T.  Brown (1992).   Use of reservoir water
        quality model to simulate global  climate change  effects  on  fish habitat,
        Climatic  Change, Vol. 20, pp. 277-296.

 Colman,  J.A.,   and  Armstrong,   D.E.   (1987).    Vertical  eddy  diffusivity

        determined with 222Rn  in  the benthic  boundary  layer  of ice-covered
        lakes,  Limnology and Oceanography,  32(3), pp.  577-590.

 Coutant,  C.C.  (1990).   Temperature-oxygen habitat for freshwater and coastal
        striped  bass  in  a  changing  climate,   Transactions  of the  American
        Fisheries Society, Vol.. 119(2), pp.  240-253.

 Dake,   J.M.K.  and  D.R.F.   Harleman  (1969).    Thermal  stratification   in
        lakes:    Analytical  and  laboratory  studies,  Water Resources  Research,
        Vol. 5(2), pp. 404^95.

 Dake,  J.M.K.   (1972).     Evaporative cooling  of  a  body  of  water,   Water
        Resources Research, Vol. 8(4),  pp. 1087-1091.

 Dettinger, M. D., and  J.  L. Wilson (1981). First  order analysis of uncertainty
        in numerical models of  groundwater flow,  1, Mathematical  development,
      , Water Resour. Res., 17(1), 149-161.

 Dhamotharan,  S.,   (1979).    A mathematical  model   for  temperature  and
        turbidity   stratification  dynamics  in shallow  reservoirs, Ph.D.   thesis,
        University of Minnesota,  p. 319.

 ERLD/MNDNR:  (1990). Minnesota  Department of Natural Resources Fisheries
        Division  Lake Data Base, expanded by US  Environmental  Protection
        Agency Environmental Research Laboratory  Duluth.

 Ellis,  C.  and  E.G.  Stefan  (1991).   Water  temperature dynamics  and heat
       transfer  beneath  the ice cover  of  a lake,  Limnology  and Oceanography,
       Vol. 36(2), pp. 324-335.
''                                    .                   '
 Ford,  D.E.  (1976).   Water temperature dynamics of  dimictic lakes: Analysis
       and prediction using integral  energy  concepts, Ph.D. thesis,  University  of
       Minnesota, p. 432.                                             -.    '

 Ford,   D.E.  and E.G.   Stefan  (1980a).    Thermal predictions  using integral
       energy  model, Jour. Eydraulics  Division,  ASCE,  Vol.  106, EY1, pp.
       39-55.

 Ford,   D.  E.,  and  E.   G. Stefan  (1980b). Stratification variability in  three
       morphometrically different  lakes  under  identical  meteorological forcing,
        Water Resour. Bull., AWRA, 16(2),  243-247.
                                     106

-------
 Gardner,  R.H.,  R.V.,  O'Neill,  J.B.,  Mankin,   J.H.,  Carney  (1981).     A
       comparison of sensitivity analysis and  error analysis  based  on a stream
       ecosystem  model, Ecol Modeling, 12, 173-190.

 Gargett,   A.E.  (1984).    Vertical   eddy  diffusivity   in  the  ocean  interior,
       Journal of Marine Research,42(2), pp.  359-393.

 Gargett,  A.E. and Holloway,  G.  (1984).  Dissipation  and diffusion by internal
       wave breaking,  Journal of Marine Research,  42(1), pp. 15-27.

 Gerhard,  H.,  J.  Imberger, and  M. Schimmele   (1990).    Vertical  mixing  in
       Uberlinger See, western part  of  Lake Constance, Aquatic Sciences,  Vol.
       53(3),  pp.  256-268.

 Gorham,   E.,  and  Boyce,  P.M.  (1989).  Influence  of lake  surface  area  and
       depth  upon  thermal  stratification  and  the   depth  of  the  summer
       thermocline,  J.  Great Lakes Res., Vol.  15, pp.  233-245

 Gu, R. and H. G.  Stefan (1990).  Year-round simulation of cold climate lakes,
       Cold Regions Science and Technology 18, pp. 147-160.

 Harbeck,  G.E. Jr.  (1962).  A practical field  technique for measuring  reservoir
       evaporation utilizing mass-transfer  theory, U.S.  Geol. Surv.  Prof. Pap.,
       272E, pp.  101-105.

 Harleman, D.R.F.,  and K.A. Hurley (1976).  Simulation of the vertical thermal
       structure of lakes under transient meteorological conditions,  Dynamics of
       Stratification  and  of Stratified  Flows  in  Large Lakes: Proceedings  of
       Workshop  of  the  Committee  on  Lake  Dynamics  (Winsdor,  Ontario:
       International Joint Commission, Regional Office), pp. 79-96.

 Heiskary,  S.A.,  Wilson, C.B., and   Larsen, D.P.  (1987).  Analysis  of regional
       patterns in lake water quality: Using ecoregions  for lake management  in
       Minnesota, Lake and Reservoir Management,  Vol 3,  pp.  337-344.

 Heiskary,  S.A.  and   C.B.  Wilson   (1988).    Minnesota  lake  water  quality
       assessment report, Minnesota  Pollution Control Agency,  St.  Paul, p.  49.

 Henderson-Sellers,  B.  (1988).  Sensitivity  of  thermal  stratification  models
       to  changing  boundary  conditions,  Applied  Mathematical Modelling,  Vol
       12,  pp. 31-43.

 Horsch,  G.  M.   and  H. G.  Stefan  (1988).    Cbnvective  circulation  littoral
       water due  to surface cooling, Limnology and Oceanography,  Vol.  33(5),
       pp.  1068-1083.

Eoughton,  J.  T.  et al. editors  (1990).   Climatic Change:   The IPCC Scientific
       Assessment.  Cambridge University  Press.

Idso,   S.B.  and   R.  D.  Jackson   (1969).    Thermal   radiation  from  the
       atmosphere, Jour,  of Geophysical Research, 74(23), pp. 5397-5403.
                                     107

-------
 Idso,  S.B.  and  R.G. Gilbert (1974).    On  the  universality  of the Poole  and
       Atkins. Secchi disk-light extinction equation,  J.  Appl. Ecol,  Vol.  11, pp.
       339-401.

 Imberger,  J.  (1985).  The  diurnal  mixed  layer,  Limnology and  Oceanography,
       Vol 30(4), pp.  737-770.

 Imberger,  J., and  Patterson,  J.C. (1989).  Physical  Limnology,  in  Advances
       in  Applied  Mechanics,  edited  by  J.W.  Hutchinson   and  T.Y.  Wu,
       Academic Press, Vol 27, pp. 303-475.

 Imboden,  D.M.,  Lemrnin,  U., Joller,  T., and Schurter,  M.  (1983).    Mixing
       processes  in  lakes:  mechanisms  and  ecological relevance,  Schweiz.   J.
       Hydrol., 45(1),  pp.  11-44.

 Jassby,  A.   and Powell,  T. "(1975).    Vertical  patterns of  eddy  diffusion
       during  stratification   in   Castle  lake,   California,   Limnology   and
       Oceanography, 20(4), pp. 530-543.

 Jones, P.D.,  Wigley, T.M.L., and Wright,  P.B.  (1986).  Global  temperature
       variations between 1861 and 1984, Nature, Vol 332(31), pp. 430-434.

 Keijman, J.  Q.   (1974).    The estimation of energy  balance  of a lake  from
       simple weather data, Boundary-Layer  Meteorology, Vol. 7,  399-407.
                    *
 Kerr,  R.A.  (1989).  1988  ties for  warmest  year,   Science,  Vol  243(17),  pp.
       891-892.

 Lewis,  W.M. (1982).    Vertical  eddy  diffusion  in  a  large  tropical  lake,
       Limnology ana  Oceanography, Vol. 27(1), pp. 161-163.

 Lewis, W. M.,  Jr., (1983). Temperature,  heat,  and mixing  in Lake  Valencia..
       Venezuela, Limnology and  Oceanography, Vol.  28(2),  pp.273-286.

 Likens,  G.E., and  Johnson, N.M. (1969). Measurement  and  analysis  of th
       annual heat budget  for the sediments  in two Wisconsin lakes.  Limnolo?,
       and Oceanography, Vol. 14(1),  pp. 115-135.

 Magnuson, J.J,  J.B.  Meisner, and D.K.  Hill  (1990).   Potential changes  ::
       thermal habitat  of  Great  Lakes  fish after global  climate   warmim
       Transactions   of  the  American  Fisheries   Society,  Vol.   119(2),  p:
       254-264.

McCormick, M. J. (1990).   Potential changes in  thermal structure  and cycle :.-
       Lake Michigan due  to global  warming,  Transactions  of the   Amer*^
       Fisheries Society,  Vol. 19(2), pp.  183-194.

McKinney,  D.  C.  (1990).  Predicting  groundwater   contamination: uncertain-
       analysis and network design, Ph.D. thesis,  Cornell University, 256 pp

McLaughlin,   D.   (1985). - A  distributed parameter  state  space  approach  v
       evaluating  the  accuracy of  the  groundwater predictions,  Ph.D.  th
-------
     a.   R.O.,  W.S.  Combs,  Jr.,  P.D.  Smith,  and  A.S.  Knoll   (19791
      Attenuation of light and  daily integral rates  of  photosynthesis  attained
      ty  planktonic  algae,  Limnology  and  Oceanography,  vol.  /4^bj,  pp.
      1338-1050.

      r.  J.D.,  J.L.  Goddier, H.A. Regier,  B.J. Shuter,  and  W.J.  Christie
      '1987).    An  assessment  of  the  effects  of climate  warming on  Great
      I^akes  Basin fishes,  Journal of Great Lakes  Research, Vol.  13(dJ, pp.
      340-352.

       l  Research  Council (1982).   Carbon Dioxide/Climate  Review  Panel.
      Carbon  Dioxide   and  Climate:     A  Second  Assessment.     National
      Academy  Press, Washington, B.C.,  p. 72.

       l   Research  Council  (1983).     Changing   climate:   report  of  the
      Carbon  Dioxide   Assessment Committee.    National  Academy  Fress,
      Washington, D. C.,  p. 496.

sCAA.  (1990).  National.  Center for  Atmospheric  Research - data  support
       section/Personal communication Roy Jenne and  Dennis Joseph,  Boulder,
       Colorado.

SrSder.   U.P.,   Schindler,  P.W.,  Wirz,  U.E.  and  Imboden,  D.M.  (1983).
       Chemical and  geochemical studies of Lake  Biel H. A chemical  approach
       to lake mixing. Schweiz. Z. Hydrol, Vol.  45(1), pp. 45-61.
Os-ood,  R.A. (1984).   A  1984 study of the water quality (T431???{^tai1
   "   area lakes,  Metropolitan Council,  St. Paul, MN,  Publ. No.  10-84-172, p.
       40.

Qs^ood,  R.A  (1989).   An evaluation  of  the  effects  of  watershed  treatment
   ~   systems on  summertime phosphorus concentration in metropohtan area
       lakes,  Metropolitan  Council, St. Paul, MN, Publ. No.  590-89-062, p. 85.

Osgood,  R.A.  (1990).    Personal  communication, Metropolitan Council, St.
       Paul,  Minnesota.

Patankar,  S.  V. (1988). Numerical heat transfer and fluid flow,  McGraw-BjH,
       New York, N.Y., 197 pp..

 Poole,  H.H.  and  W.R.G.  Atkins  (1929).   Photo-electric measurements of
        submarine illumination  throughout the year, J.  Mar. Biol. Assoc.  U.K.,
        Vol. 16,  pp. 297-324.

 Priscu,   J.C.,  Spigel,  R.H.,  Gibs,  M.M., and  Downes,   M.T   (1986).   A
        numerical   analysis    of   hypoh'mnetic   nitrogen  . and   phosphorus
        transformations in Lake Rotoiti, New Zealand: A geothermally mfluenced
        lake,  Limnology and Oceanography, Vol.  31(4), pp. 812-8.J1.

  Protopapas,  A.  L. (1988). Stochastic hydrologic analysis  of soil-crop-dimate
        interactions,  Ph.D. thesis,  Parsons  Lab., Dep.  of Civ.  Eng.,   Mass.
        Inst.of Technol.,  Cambridge, 239  pp.
                                      109

-------
 Protopapas,  A.L.,  and  R.L.   Bras  (1990).  Uncertainty  propagation  with
       numerical models  for flow and transport in the unsaturated zone,  Water
       Resour. Res., 26 (10), 2463-2474,  1990.

 Quay,  P.D.,  W.S.  Broecker,   R.H.  Hesslein,   and  D.W.  Schindler  (1980).
       Vertical  diffusion  rates determined by tritium tracer experiments in the
       thermocline and hypolimnion of two lakes, Limnology  and. Oceanography,
       Vol. 25(2),  pp. 201-218.

 Richardson,  C.  W.  (1981).   Stochastic  simulation   of  daily  precipitation,
       temperature,  and  solar  radiation,  Water  Resour. Res.,  17(1), 182-190,
       1981.

 Riley, M.J. and H.G. Stefan (1987).  Dynamic  lake  water  quality  simulation
       model  "Minlake",  University of Minnesota, St.  Anthony  Falls Hydraulic
       Laboratory,  Report No. 263, p. 140.

 Robertson,  D.M. (1989).   The use of lake water  temperature and ice cover as
       a climatic  indicators,  Ph.D. thesis, University of Wisconsin-Madison, p.
       330.

 Robertson,  P.M.,   and  R.A.  Ragotzkie  (1990).    Changes  in  the  thermal
       structure of  moderate to  large size lakes  in  response to changes in air
       temperature, Aguatic  Sciences,  Vol. 52(4), pp.  362-379.

 Robinson,  P.J.,  and  Finkelstein, P.L.  (1990)  Strategies for  the development
       of climate scenarios for impact assessment:  Phase I final report, USEPA
       Atmospheric    Research    and    Exposure   Assessment   Laboratory
       EPA/600/53-90/026.

 Sadhyram,  Y.,  P.  Vethamony,   A.  Suryanarayana,  G. N. Swamy,  and J. 5
       Sastry  (1988).   Heat  energy  exchange  over a  large water  body und*-:
       stable atmospheric conditions,  Boundary-Layer Meteorology, Vol. 44,  ~,
       171-180.

 Scavia D., W. F. Powers, R.  P. Canale, and J.  L. Moody (1981).  Compare =
       of   first   order   error   analysis   and   Monte   Carlo   simulation  -.-
       time-dependent lake  eutrophication models, Water Resour. Res., 17u
       1051-1059.

 Schindler, D.W., K.G. Beaty, E.J. Fee, D.R. Cruikshank, E.R.  DeBruyn, !  -
       Findley, G.A. Londsey, J.A. Sherer, M. Stainton, M.A. Turner  (IS*
       Effects  of  climatic  warming  on  lakes  of the  Central  Boreal  Fc-r*.
       Science, Vol. 250(16), pp.  967-970.

Shapiro,  J.  and H.  Pfannkuch   (1973).   The  Minneapolis  chain  of lakn
     .  study of urban drainage and its effects 1971-1973, Limnological n
       Center, University  of Minnesota, Minneapolis, Report No.  9.

Smith,  C.R.  and  K.S.  Baker   (1981).   Optical properties  of  the
      natural waters, Applied Optics, Vol. 20(2),  pp.  177-184.
                                     110

-------
         H H..  J.  Imberger,  and  K.N.  Rayner  (1986).   Modeling  the diurnal
        ^iicd  layer,  Limnology and Oceanography, 31, pp. 533-556.

        ;  R.E., and  Armstrong, D.E.  (1983).   Lake  mixing and  its relationship
        ui  epilimnetic  phosphorus in  Shagawa  lake,  Minnesota,  Can.  J.  Fish.
        .\
-------
Wdander,  P.  (1968). Theoretical  forms  for the vertical exchange coefficients in
       a stratified fluid with  application  to  lakes  and  seas,  Acta  Regiae
       Societatis  Scientiamm  et  Litterarum Gothoburgensis, Geophysica 1, pp.
       1-26.

Wu,  J.  (1969).  Wind  stress  and  surface  roughness  at  air-sea   interface,
       Jour. Geophysical Research, 74, 2, 444-455.

Ward,  P.R.B.  (1977).     Diffusion  in  Lake  Hypolimnia,  Proceedings   17th
       Congress,    International    Association    for    Hydraulic     Research,
       Baden-Baden,  IAHR, pp.  103-110.

Willis,  R.  and  W.W-G.   Yeh  (1987).    Groundwater  Systems  Planning and
       Management, Prentice-Hall,' Englewood  Cliffs, NJ 07632, p. 415.

Winter, T. (1980).   Personal communication, U.S. Geological  Survey,  Denver
       Federal  Center, Lakewood,  CO 80225.

Wright, D.} R. Lawrenz, W.  Popp, and M. Danks  (1988).   Acid precipitation
       mitigation program,  Progress  report Thrush Lake  Minnesota,  Minnesota
       Department  of   Natural   Resources   Division  of   Fish   and  Wildlife
       Ecological  Services Section, p. 177.
                                    112

-------
Appendices

-------
                                APPENDIX A

                         Vertical  diffusion in a small
                    stratified lake: Data and error analysis


      Water temperature  profiles  were measured at  2 minute  intervals in a
 stratified temperate  lake  with  a  surface  area of 0.06 km2 and  a maximum
 depth of 10  m from  May 7  to  August  9, 1989.   The  data were  used  to
 calculate the  vertical eddy diffusion coefficient  (Kz)  in the hypolimnion.  The
 depth is representative of a large  number  of lakes in the north central  United
 States.  (KE)  was calculated  over time intervals of  1  to  15  days and  varied
 from  10"3  to  10"1  cm2s~l.    A  numerical  model  was  developed  for heat
 conduction  in  the  sediments, and  heat flux  between water  and sediments was
 incorporated into the relationship  from  which eddy  diffusivity was estimated.
 Heat flux between water and  lake sediments, a term  commonly neglected, was
 found to be  important  in the Kz estimation.   Kz  values  were  related  to
 stratification  stability  as  measured by  the  Brunt-Vaisala  frequency  N  using

 Welander's  expression  of  the  form  Kz =  o(N2)T  Values of a were on  the
 order  of 10'4  and  7 varied from -  0.36 to -  0.45  when Kz was  given  in
 cm2s'1 and  N  is in  s"1.   An  error  analysis  was  conducted and the effects  of
 different  choices  of  sampling  intervals  in time  and depth  on the  eddy
 diffusivity estimates  were  evaluated.  The  longest time interval (15 days) and
 the smallest depth increment (1 m) used in  this  study  were found  to  give  the
 best Kz  estimation.
A.1   Introduction

      Density  stratification  due to  vertical  temperature  gradients  inhibits
vertical  mixing  in  lakes  and  reservoirs,  and mixing in  turn  affects  the
distribution of phytoplankton, nutrients, and  other  water quality  constituents.
Quantifying  turbulent transport phenomena is one  of the major challenges in
lake and  reservoir hydro thermal and water quality analysis.   Specification of
vertical turbulent (eddy) diffusion coefficients  in one-dimensional  water  quality
models, which are  often used for decision-making, is particularly  difficult.

      In  the  analysis of vertical  turbulent  mixing  by one-dimensional  wind
energy models, the depth of the surface  mixed (epilimnetic)  layer  is calculated
by  an integral  entrainment  model while  the  vertical  transport  in  the
hypolimnion  is  taken into  account by  a diffusion equation (Stefan and  Ford
1975;   Bloss   and  Harleman  1979;  Ford and  Stefan  1980).    Although  the
hypolimnion  is isolated  from the epilimnetic layer by the thermocline and its
associated  density  gradient,  strong  and  sporadical local  mixing  events  have
been  observed in the hypolimnion (Jassby and Powell  1975;  Imberger 1985;
Imberger  and Patterson  1989).   Such mixing events  can  originate  from
oscillating boundary layers  induced by seiche  motions on the bottom of lakes,
                                     A-l

-------
 internal wave interaction and breakdown, shear instabilities  in the thermocline
 (billows),  epilimnetic turbulent kinetic  energy  leakage to the hypolimnion  and
 double  diffusion  processes.     Scales  for  such  events  range   from  the
 Kolmorgorov scale to the  lake  basin  scale.   Eddy  diffusion  dependence  on
 stratification strength as measured by  buoyancy frequency has  been pointed
 out consistently  (Colman and  Armstrong  1987; Gargett  1984;  Gargett  and
 Holloway  1984;  Imboden et al.  1983;  Jassby and  Powell  1975; Quay  et  al.
 1980;  Welander 1968).

      Direct measurements of vertical turbulent diffusion in  lakes are not easy
 because of the 3-D nature of the diffusion  field, and  the spatial  and temporal
 scales.   To estimate  diffusion values, one  can rely  on measurements of water
 temperatures or  concentrations  of natural tracers present in  the water.  Water
 temperature measurement is the most  commonly used  method because  of its
 simplicity;  however,  a  careful  assessment  of all  external  and  internal  heat
 sources is  required.

      The  purpose of this study  was to estimate vertical eddy  diffusion from
 water temperature measurements  in a  typical inland  shallow lake.   Sediment
 heat  exchange, commonly neglected along with error analysis, is  also included
 in  the estimation.    Lastly,  criteria for measurement  intervals  in  space  and
 time that  minimize  the error  in eddy diffusivity  estimation  are  proposed.
 The  latter  uses  principles  which  are  also used  in  groundwater  monitoring
 network design (Andricevic 1990).


 A.2  Study Site

      Ryan Lake,  located in Minneapolis,  Minnesota,  has  a  surface  area of

 0.06 km2,  mean depth of 5  meters and maximum depth of  10.5 m  (Fig. A.I).
 The  lake,  located  in a flat  terrain,  suburban  residential  area,  is  highly
 eutrophic,  and regularly  experiences winterkill of fish.   The maximum  depth
 of  10 m is equal  to the  median  maximum  depth of 779 statistically analyzed
 lakes in Minnesota.   The depth  of Ryan  Lake can be considered  as  typical
 for the north central  United  States.

      Lake  water  temperatures  were measured  every  two  minutes at 1  m
intervals from the lake  surface  to the 10  m  depth.   Every 20  minutes,  the
 previous ten  measurements  were  averaged  and  recorded.   The  measurement
 scheme  was  adopted   to   reduce  the   "high  frequency"  electronic   and
measurement  noise, while retaining  the  fluctuations expected at  timescales of
hours  and  days.    The  probes  are rubber coated  thermisters  with a time
constant of 10 seconds.    They were  calibrated in  a  water  bath prior  to
installation.   Absolute accuracy  (values' measured  by two  adjacent  probes at
known temperatures)  was   0.05 C (95%   confidence  interval),  while relative
accuracy (the difference  between  successive  measurements by the same probe)
was O.OrC.  A Campbell Scientific CR10  datalogger installed on a  small raft
recorded the water temperatures.  Hypolimnetic data for the period  from May
7  to August 9,  1989,   were selected  for  analysis because this period was
characterized  by  stable  seasonal  lake  stratification.   In 1990  measurements
were extended to  sediment temperatures using  probes identical to those in  the
water.  The sediments were  soft,  organic material  and poorly consolidated, as
                                    A-2

-------
                           10m max *  I I

                               9m 
                       MEASUREMENT
                           LOCATION
Fig.  A.I    Ryan Lake bathymetry.
                   A-3

-------
indicated by  the  ease  with  which  the thermister  probe  support  rod  was
installed.
A.3  Vertical Eddy Diffusivity

      Many studies  have assumed lake  basins to  be  closed  systems  and have
estimated  an  average vertical eddy  diffusion  coefficient  over the whole  basin
(Adams  et al.  1987;  Gargett  1984;  Imboden et al.  1983; Jassby  and Powell
1975; Lewis 1983; Nyffeler et al.  1983;  Priscu et  al.  1986; Quay et  al. 1980).
One  of  the methods for such  estimation  is through the  budgets  of . scalar
quantities  such  as temperature (Gargett 1984).

      The  one-dimensional, unsteady heat  transport equation applied  along  the
vertical axis of  a water  column is:
The flux-gradient method for  the computation of Kz reduces this  equation  to
the form


            K,   [ f? ]"' { It  /  T(Qd< + T|S - / Sdz }        (A.2)
     It is assumed that  there  is  no vertical  advection  (w  =  0)  of water
anywhere.  There is, however, a conductive heat flux H  , from the sediment

into the water at the  lake  bottom (z =  0)  and a radiation  (penetrative)  heat
flux H ..   A  heat  balance for the  water  column  between z  =  0  and  z
       sol
therefore leads to a replacement  for equation (A.2).
                                                                       (A-3)
                               3t

Vertical kinematic thermal  eddy diffusivity can be explicitly expressed  as:
                               z                      z
                             /
                       _
                       dt
                K2 = 	-	2	2	           (A.4

                                     Jz

where  T(z,t)  =  measured  water  temperature  distribution,  z  =  upwar-
coordinate starting at the lake  bottom, t =  time, w = vertical component c:
velocity,  and S  = internal  source term. At the  time scale of 1 to 15  days ^
which  Kz is computed,  and without  significant inflow  or outflow to  or iio
the lake,  net vertical velocity w is customarily small enough to be  neglectec
Short term  effects during storms and turnovers  will show up implicitly in  th-
value of  Kz.  The source  term  "S" in  Eq. A.4 accounts for solar short wav-
                                    A-4

-------
radiation absorbed  in the water  and heat  flux through  the  water column to
or from the  sediments at the bottom.   For  shallow  inland  lakes, the source
term can be  particularly significant.

     As  pointed   out  by  Gargett  (1984),   the   budget  method   has   two
advantages.    First,  few additional  assumptions are  needed  to  estimate Kz
from Equation A.4.  Second, time averaging is implicit in the estimate of Kz.
With  the  exception of the surface mixed layer, turbulence in lakes occurs in
patches and   intermittently.   Turbulent  "bursts"   involve  small  volumes of
water  (tens  of  cubic meters) and last on the order  of minutes  (Imberger
1985).   Therefore,  time averaging for such  systems appears to be essential to
capture the long-term behavior.

     Eddy diffusion  dependence  on  buoyancy frequency  was pointed out  by
Welander  (1968) and others.  Welander  derived an expression relating Kz to
the square of the  Brunt-Vaisalla frequency (N)  as Kz  = a (N2)^,  where  N2
         &    p ^  density of water  and  g  = acceleration of gravity.   If
turbulence  is generated by  the dissipation of energy from large-scale motions,
7 = -1.0;  otherwise, if it is generated by shear flow, 7 = -0.5.   Welander's
analysis  was  very  informative, but  it  was  based  on  several  assumptions:
steady state,  no  boundary  effects,  and  a  linear  dependence of  density  on
temperature.   Such assumptions are only  marginally valid for  lakes.   The
results to be presented herein will show that Welander's theory fits lake data
reasonably  well.


A.4  Sediment Heat Storage

     Few previous analyses  include heat flux to or from the sediments in the
eddy  diffusivity estimation  for  summer  conditions.   A  notable exception is
Stauffer  and  Armstrong's  (1983)  study  of Shagawa  Lake's western  basin
(maximum  depth  14 m).   In principle, sediment  heat flux is related  to the
water  temperature gradient at  the  sediment/water  interface (Nyffeler  1983);
however, unknown turbulent  heat  transfer coefficients relating the  flux to the
gradient   as  well  as  exceedingly  small   temperature   gradients  in  the
near-sediment  water limit  the usefulness  of the  relationship.    Relying  on
measurements and computer simulations,  Priscu  et al.  (1986)  assumed that
the heat  flux from  the sediments  to the  water  was  constant.   This was
physically  justified for the  geothermally  influenced  lake which they  studied.
In the more general situation, conductive heat flux through the sediments is
variable in  depth  as well as with time  (Birge et al.  1927; Likens and Johnson
1969).

     In  this  study, a numerical  model  was developed  to  simulate sediment
heating or  cooling by  the  overlying water.   A one-dimensional,  unsteady heat
conduction  equation  was  applied  since  conduction  into   and  out  of  the
sediments  is  essentially a  1-D  process.    The   unsteady heat   conduction
equation for  the  sediments  is a  simplified "version of  Eq. A.I.  Vertical
velocity, w, is zero because  there is no  advection and  S  = 0 because there
are no internal heat  sources  or sinks  in the sediments.   Kz in  Eq.  A.I is
replaced by Kzs =  sediment thermal  diffusivity,  and T  is  replaced by  Ts =
                                     A-5

-------
sediment temperature.   So  5Ts/dt  = Kzs(52Ts/&2)  js  the  heat conduction
equation applied to  the sediments.   The partial  differential equation  was
discretized  using a control  volume  method  (Patankar 1988) and  solved  by a
tridiagonal  matrix  algorithm.   The boundary  conditions are:   (1) measured
water temperatures at the water/sediment interface and  (2) no flux (adiabatic
boundary)  at  zj  =  6 m  depth below  the sediment  surface.    It could be
shown  by unsteady heat transfer analysis  of a  semi-infinite slab  that  seasonal
heat storage  did not penetrate significantly beyond  6 m depth in an annual
cycle.   Heat  flux (Hsed) through the sediment/water  interface is  calculated as
the rate of change  in  sediment heat storage given by  integration  of computed
sediment temperature profiles T(z,t):                         -
                                          Zd

                      Hsed  = Ps CPS "If  / Ts(2'1)  dz                  (A-5>
                                         0

where  ps  bulk sediment  density, cps is  sediment  specific heat  and  (ps cps)
is  bulk specific heat of the  sediments per  unit volume.

     Time series of measured sediment and overlying water temperatures  are
given in Fig. A.2  down  to depths of 1.5  m into the sediment.   No probes
were placed at  any greater depths.   These temperatures were recorded  from
April  3  to July  9,  1990.    In the  early  part of this season,  temperature
gradients, and  hence  heat   fluxes   into the sediments,  are at  a maximum.
High fluctuations  of water   and sediment  surface temperatures correspond to
the spring  overturn.   Water  temperatures  at 0.5 m above  the sediments and
at the sediment  surface plotted in Fig.  A.2 were almost identical, indicating
the  presence  of significant  turbulent  mixing  in  the water  boundary  layer.
The differences  between  20  minute readings  at  the  two  elevations  had an
average of 0.00845" C  and  a standard  deviation of  0.088C.   The square of
the correlation coefficient (R2) was  0.98.

     The  unsteady  sediment  heat  conduction  model  was  calibrated  for

thermal diffusivity.  A sediment thermal  diffusivity of 0.0022 cm2sec"1 applied
uniformly over  depth  (Gu  and  Stefan   1990)  simulated  measured  sedimen*.
temperatures  well  (Fig. A.3).   The  maximum  difference between calculated
and  measured temperatures  was 0.2"C  in  a range of 1.7 "C.  Accuracy of the
measurements was  0.05'C.   The maximum  discrepancy was  observed at  0.5 ~
below  the sediment surface in April. For  the rest of the period the difference
was  less than 0.1 C.   Since  sediment thermal  properties do  not  change wi:::
season, the calibrated  values should hold for year-round  simulation.

     Using  the calibrated model and measured deep  water  temperatures  frcr:
the 1989 data as  the  boundary condition  at the water sediment  interface,  ti=
sediment temperatures for  the period of lake  eddy diffusivity  studies  (May
to  August  9, 1989) could  be simulated.   Heat flux  from  the" sediments w^j
calculated using the simulated  sediment temperatures  in equation (2).  Likerj

and  Johnson  (1969) used values of ps=l-2 g cm"3 and cps=0.8 cal g'^C"1 i"
soft  bottom  sediments,  giving (ps  CnS) =  0.96  cal cm'3  'C'1.   A  water v
solids  ratio of 3:1  corresponded to  the  calibrated thermal diffusivity  (Carsu

and Jaeger  1959),  and hence values of 0.8 cal  g'^C"1  for specific  heat and 1  '


                                    A-6

-------
     7
     6-
O
o

 d>
 L_
 Z3
-t~'
 D
 v_
 CD
 QL

 E
 
-------
                                                      calculated
                                                      measured
                                                      measured
                                                      measured
April 2      April 10      April 18      April 26       May 4       May 12
              \lruhu-il and measured  temperatures in lake sediments,  1990.

-------
                                    accuracy
equation (2) is  estimated to be 2  kcal nrMay1.
A.5   Water temperature observations

      Measured  hypolimnetic  water temperatures based  on 20 minute averages
are plotted in Fig. A.4.  Depths  are below the water surface which fluctuated
by  less than 0.1  m.  There  were no spatial variations of water temperatures
to speak  of, other  than in the vertical  direction.   For  the entire period of
record, stratification was stable.   Above  6 m depth, water temperatures were
influenced by the deepening  thennocline.  'The rise in  water temperature at
the  4  m  depth  in July  was  due  to  epilimnetic  wanning associated with
increased  solar  radiation and  deepening of the mixed layer.

      Daily  time  series  of on  site  incident  solar radiation,  air temperature,
wind  speed, 'and  epilimnetic  water temperature  are  given  in   Fig. A.5.
Following  standard  weather bureau procedure, solar  radiation is a  daily total,
wind speed  is an average of  three-hourly  readings, and air temperature is the
mean  of a daily  maximum and  minimum  reading.   The strong  rise of suriace
water temperature from May 7 to 17 coincides with high solar radiation and
low  wind.   Water temperature fluctuations  from May  17  to  June  23  are the
result of  high  fluctuations  in  solar radiation  and  air  temperatures.   From
June  12  to the  end of the  observation  period, high  solar  radiation  and air
temperatures increased  water  surface  temperatures.

      The  entire  stratification  dynamics  are put  in  perspective in Fig. A.6.
The  isotherms  were developed from the  water temperature records  such as
plotted  in  Fig.  A.4.    The  window  of  data  analyzed  for  vertical  eddy
diffusivities  is shown in Fig.  A.6.


A.6  Vertical eddy  diffusion estimates

      With  temperature T(z,t)  given  by  discrete  measured  values   Ti,j  at
increments Az  and At  in  depth and time,  respectively,  Equation  (A.4) nust
be  discretized numerically  to yield  the  eddy diffusion  coefficient estimate: in
the  form:
where i  and N are the bottom and  topmost  layer of the  lake,  respectively,
ATt  = water temperature  difference over a time interval At at a fixed depth
z,  Hsoi  =  solar  radiation  heat  flux  at  depth  z,  Hsed  = water/sediaent
interface  heat flux, and AT  = temperature difference over a vertical distance
Az  averaged over  the time  interval  At,  p  =  density  of water and c?  =
specific heat of water.
                                     A-9

-------
'
o
                          May 7     Moy ?l    June A   June Ul    July 2     July  16    July .50  Aucjust  tJ
                          Fig.  A.4     Hypolimnetic  lake water  temperatures recorded  at  2 min interval
                                       in Ryan Lake, May 7  to  August 9, 1989.

-------
       800-
                  1 m depth water temperature
                  oir  temperature	_^
                                                              wind  speed
                                                              _j    ,_i   j
8.0
                                                                           -6.0
                                                                           -4.0
                                                                           -2.0
-0.0
                                                                                 in

                                                                                 2
          May 7     Mo/ 21     June 4    June 16    July 2    July 16    July JO   August 13
Fig. A.5      Meteorological   conditions  (solar  radiation,  vnnd  epetd   and   air
               tcmpcratnrcm) during lh  period of in ventilation.

-------
to
                E
                o.
                
-------
     Heat   flux  through  the  water-sediment   interface  was  calculated   by

Equation (A.5).  The average heat flux was 7.0  kcal  m2 day1,  and the actual
time series is  shown  in Fig. A.7.  The heat flux was from the  water into  the
sediment,  i.e.  the  sediment  acted as a  heat sink  throughout the  period of
investigation (May  7 - August 9).

     Internal  solar  radiation  absorption  was calculated for each  depth from
measured  radiation  at the  lake surface  and  an  attenuation  coefficient (Eq.
1.4).   Bi-weekly Secchi  depths varied  from 0.8  m to  1.25  m  during  the
period  of  analysis  (May  7 to  August  9,  1989) with  a  mean of  1.0   m.
Relationship  between  total   attenuation  coefficient  and  Sechi   disk   depth
translates  a Secchi  depth  of 1.0 m into  an attenuation coefficient of  1.8  m*1
which  was used  throughout  the analysis  (Fig. 2.4).   Solar radiation adsorbed
at the 4  m depth  amounted  to  about one  third of the sediment  heat flux
(see  Fig. A.7)  and was  less at greater depth.   In general, if an internal heat
source  due to solar radiation exists but  is  neglected,  the  eddy  diffusion
coefficient  will  be  overestimated.   Although not shown in Fig.  A.7, solar
radiation  and  water to sediment heat  flux  had  different signs  in Equation
(A.6) because absorbed solar radiation is  an input to the water and heat flux
to the  sediments is  a loss from the water during the period  of  study.

     Vertical  eddy  diffusion  coefficients  calculated  for  sampling  intervals of
five  days  and depth increments  of  1 m (Fig.  A.8)  show  decreasing  values
with time  as  seasonal stratification progresses.   High variability in space  and
time is apparent.    Vertical  eddy diffusivity coefficient  values  ranged from
approximately  0.001 to 0.1  cmVi with an  average near  0.01 cmV1.   The
highest eddy   diffusivity  was  found   at  the  greatest  depth  (near  the  lake
bottom)  while the 4 m depth had the lowest values.   Eddy  diffusivity near
the lake sediments  is produced mainly by  the interaction  of  internal waves
with the lake  bottom resulting in internal breaking  waves and by turbulence
induced by bottom  shear during internal  seiche  motion.   All of  these  are
contributing to an intensely mixed lake bottom boundary layer, as previously
shown  by  the temperature records in Fig.  A.2.   One result of this  mixing is
the decrease in stratification intensity with greater depth in the  hypolimnion.
There  is  also  a  positive  feedback  because  shear-induced  turbulence  is
dampened  by density stratification conditions.

     Eddy  diffusion coefficient estimates versus  static  density stability (N2)
for  different   sampling   periods  with   and  without  consideration   of  the
sedimentary heat source term  are given  in Fig. A.9.  A least squares  linear
regression  was used  to  estimate coefficients a  and 7 in  the   relationship  Kz

= Q(N2)7  .   As expected,  an inverse  relationship between Kz  and N2  was
observea.    When the sediment heat  flux term  was  omitted, eddy diffusivity
was  underestimated.  The error was  up  to one  third of the estimated  values.
It is noteworthy that a stronger dependence  of  Kz on  N2 was  observed when
the sedimentary  heat source term was considered (Fig. A.9).
                                  ' A-13

-------
I-1
>.
Q
fN
CJ
 X
 Z)
 _J
 U_
 Id
 n:
                  water-sediment  heat flux

                  solar radiation heat  flux  4m  depth

                  water temperature above sediment     y
                  water temperature above sediment  ../""
F-5.4

 r5.3

 r5.2

 1-5.1

 r5,0

 r4.9

 I-4.8

 j-4.7

 j-4.6

 r4.5

 r4.4

 r4.3

 r4.2

 -4.1
                                                                                    O
                                                                                    o
UJ
C
ID
h-
<
cn
Ld
Q_
^
LJ
                                                                              4.0
         May 7    May 21    June  4   June 18   July 2    July  16   July 30  August  13
        Fig. A.7     Heat  flux  through  the   sediment-water   interface  and   solar
                    shortwave radiation received at the 4m  depth, May 7 to August

-------
                              10
                                                      RYAN LAKE

                                          TURBULENT DIFFUSION COEFFICIENT TIME SERIES
ii
en
                         I
                         u
                         O)


                       "E
                        u
                              ID"
IO'J
                            10'J
                            10-
                              Apr 30
              	7 M

              	 6 M

              	 5 M

              	 4 M DEPTH

                        i i [ i 'i i  i i i i i i j i i  r ' i r i i i i [  i i i i i i i i
            May 20      Jun 9     Jun 29      Jul 19      Aug 8
                         Fig.  A.8     Vertical turbulent diffusion coefficient time series in Ryan Lake,

-------
 o
 0>
 00
 CM


 O
 o
 4>
 CO
 u
 
                 R'=0.32 . T=1 doy 
                 SEDIMENT HEAT NOT
                         K =1.15X10""(N2)-041
                 RJ=0.67 . 7=5 doys*     ,

                 SEDIMENT HEAT NOT CONSIDERED
                 R=0.87 , T=10 doys

                 SEDIMENT HEAT NOT CONSIDERED
                                                     rrj   i .-..
-------
        A.I   Regression coefficients for Kz =
Sampling
Interval
 (days)
1
5
10
 15
Hsed * 0
a 7
9.14*10-5
9.77*10-5
1.12*10-<
1.37*10-<
-0.45*0.017
-0.45*0.022
-0.44*0.018
-0.43*0.018
Hsed = 0
a 7
2.20*10-5
1.15*10"*
1.47*10-4
1.87*10'*
-0.360.041
-0.41*0.028
-0.40*0.03
-0.38*0.02
     Kz is a bulk estimate of the diffusivity over a  time interval rather  than
 i.n estimate of  an instantaneous value.  Variability of the eddy diffusivity was
 .he highest for a Campling interval  of one day.    Different regression  lines
 could  be fitted  to  the  one  day results  without  changing  the  regression
 coefficient R2.    By increasing the  sampling interval, the effects of  variable
 meteorological  conditions  and  mixing  events  were averaged  out  over longer
 and longer periods.   With a longer sampling interval, the linear  regression  fit
 was better.   Regression coefficients  with standard errors are  summarized  in
 Table A.I.
A. 7  Error Analysis

     Uncertainties  in  the  estimated  Kz  values  are  introduced  by  water
temperature  measurement   errors,  model   parameter   values,   and  model
formulation.   Magnitudes  of errors  for  key variables  in Equation  (A.3)  are
listed in Table A.2.   The first  three of these errors are measurement errors.
The  last value is based on  uncertainty in  estimates of  specific heat and  soil
temperature measurements.   2.0  kcalm*2  day1  is  about  30%  of  the average
net heat flux value of 7 kcal nr2dayl.

  Table A.2 Errors
         Variable
Symbol
  Error
  Temporal temperature
      differential ATt (measured)
  Depth temperature
      differential ATZ (measured)
  Depth increment Az (measured)

  Source heat flux AH (calculated)
   Az

   c
   s
0.01 (C)


0.05 (-C.)

0.01 (m)

 2.0 (kcalnr'day1)
                                    A-17

-------
     To assess the estimation error of eddy diffusivity, a small perturbation e
of the variables in  Table A. 2  is  introduced  into  equation  (A. 6).   Temporal
and  depth temperature differentials,  depth increments, and  source heat fluxes
are considered  random variables  and represented by  the  mean  value plus  a
perturbation  term.   The  perturbation  terms  have zero  mean,  and  standard
deviations equal to  one-half the  values given in  Table A. 2.  Mean  (denoted
by overbar) and perturbation .(denoted by e) were expressed as:


                              Kz  = K7  + ekz           '               (A.7)

                             ATt = ATT +  et                         (A.8)

                              Az  = Az + e                            (A.9)
                   /.  '          S = 5 +  eg    -                       (A.ll)


with the statistical properties

                           E(e) = 0 (mean of e)                      (A. 12)


                         E(f2) =  
-------
VAR(KZ)=

                                      A2iA2j
          e2e2
          tAz
       4Az2
       AT2
               AT
                     4Az2

                     AT2
                       z
                     1   N
                    i-  S
                    At J = 1
                                        AZ
                                TjAZj
                                             4 AZ
                   Az
AZ<
                 AT
                     4Az
                           Az
                                         N
                                         E
At
          16
                N
        AtAz3  ia,
                       AZ
                                                                (A.15)
where  i or j designate  the  depth  under consideration,  Az = depth increment,
ATj  or   ATj  =  temporal  temperature  difference  at   the  depth  i  or  j,
respectively.   N  =  total number of measuring points below  the flux  surface
under  consideration, ATZ = temperature difference over  the  depth increment
Az averaged over the  time interval At, S is source term. Other variables are
given in  the  main text  (Table A.2).

     Vertical profiles  of the mean  eddy diffusion coefficient plus or  minus
two  standard deviation intervals  are given  in Fig. A. 10.   These  profiles are
for five day  sampling  intervals  and depth increments  of  1  m.   The error in
Kz estimation  increased with  depth mainly due to  the  smaller temperature
gradient  near the lake bottom.

     The dependence of calculated KZ  values on the frequency and spacing of
the water temperature measurements is illustrated  in Figs. A.11  and  A.12.   If
sampling  intervals  exceeded four  days, the  error in estimated Kz values did
not change significantly, regardless of depth.  Sampling intervals  of three days
or less increased  the error.

     The tradeoff between  depth  and time  intervals with regard to  the error
in Kz  estimation is illustrated  in  Fig. A.12 by isolines  of  equal 1a\z  values.
Three  regions  can be  distinguished  on  the graph:   (1) up to  a  sampling
interval  of three days  the  error  was a function  of the  time increment only.
The bigger  errors  correspond  to  the smaller sampling  time increments. (2)
From  5  to  10  days sampling interval,  errors were  a function of both Az and
At. That  was the region where  trade-off between space and time was possible
in  order to  obtain  the   same  estimation error.  In  general,  errors  were.
                                   A-19

-------
    3-


    4-


    5-


1  6-
x

a  7-
u
o

    U


    9J


   10
Jun 9,  1989
                       -2(T
    0.000  0.010  0^620  0.030   0.040  0.050 'o.060
    4



    5



1  6

x
h
Q.  7.
UJ
O

    8



    9



   10
                    K^cm'sec"1)
                                                             Jun 29, 1989
 I  '
mean

+ 2
-------
    O.O70
    0.016-

 '
b
(N
    0.012-
    0.008-
    0.004-
     0.000
                       i  '  i '  I  '  l  '  l
                                 m  depth
                               5 m
Storting date: Jul 19,  19Q9  	 6

                               7 m
                           	8 m
           1   2   3   4   5   6   7   8   9101112131415
                         TIME  INCREMENT (DAYS)
        Fig. A.ll   Estimated  eddy diffusion  errors  (2crK ) for  different  sampling

                   intervals.

-------

to
                   3.00
              -CN   2.60
                   2.20
Z
UJ
2
LJ
cr
o
     1.80  -]
              ^-
              o_
              UJ
              Q
      .40
                   1.00
                            88888888  88
                          1    2   3   4   5   6   7   8   9   10  11.12   13  14  15
                                        TIME  INCREMENT  (DAYS)
                         Fig. A.12    Estimated eddy diffusion errors 2aK (cmV1) space-time tradeoff,

                                    7 m depth Jul 19.

-------
decreasing  for  smaller Az and  larger At.   (3)  For sampling  intervals larger
than 10  days  the estimation error became  primarily a function of Az, i.e.  the
more measuring points used in a profile,  the smaller the error  in  Kz.


A. 7  Conclusions

     The vertical turbulent eddy  diffusion  coefficients in the hypolimnion of a
temperature  stratified  temperate  lake  with  a  depth  typical of  the north
central United States were evaluated  from water  temperatures  measured at 2
minute  intervals  from  May 7  to  August  9,  1989.    Kz  values  typically
increased  by  a factor  of 10 between 4 m depth  and  9 m  depth.  Eddy
diffusion  coefficients  Kz  were  computed  for periods  from  1  to  15  days  and
varied  from 0.001 to  0.1  cm^/s for the 1-day intervals and from 0.002 to 0.04
cm2/s  for  15-day intervals.    Kz  also varied  with  stratification  stability
measured by the Brunt-Vaisala  frequency  N.  The relationship  Kz =
produced  the  best  fit  when  a = 0.00014  and 7  = - 0.43  for  periods of 15

days, where Kz is  in  cmV1 and N in s~l  .   As the time step  was  shortened
to one  day, the  fit became  poorer  and  7  values  changed slightly (Fig.  A. 9).
The   value  7  =  -  0.43  fits  Welander's  (1969) model  for  shear  driven
thermocline erosion.  The a  value is related  to lake size (Fig.  2.2).

      The water temperatures measured  and  recorded every 20 minutes at  the
sediment /water  interface  and  at 0.5  m above  showed a mean difference of
only 0.008' C and nearly  identical responses  in time (Fig. A.2)  over  a  period
of three months.   This is indicative of a well-mixed boundary layer  near  the
lake  bed.

      Heat exchange between  water  and  lake  sediments  was  found  to  be
important to the analysis of vertical  thermal diffusivity.  A numerical  model
was  used to estimate sedimentary  heat  flux for  incorporation into  the  eddy
diffusivity  estimation.   A mean value of  sedimentary  heat  fluxes during  the
summer period (May - August)  was 7 kcal  nr2dayl.   Estimates of Kz  made
without sedimentary heat flux were up  to  one third smaller than those  made
with  that heat flux.   Heat input by solar  radiation through the water surface
also  influences  the  estimates  of Kz,  but  this  influence diminishes with depth.

      Effects  of errors in  temperature measurements  and sediment  heat  flux
estimations  on  calculated  vertical eddy  diffusion   coefficients  were evaluated.
Estimation  errors   were  much larger   near  the   lake   bottom   (in   the
hypolimnion)  than  in the  thermocline  region  (Fig.   A. 10).    The  smallest
estimation errors of the eddy  diffusivity were obtained for sampling  intervals
of 15 days and  depth  increments of 1.0 m.   At  the 7  m depth,  the  error was
about  0.0011 cm2/s relative  to  a value  on  the order  of 0.010  or 11 percent
(see  Figs. A.8  and A. 11).  The  error doubled when the depth increment was
trippled to 3.0 m or when the sampling  interval was reduced from 15 days to
5 days  (Fig. A. 12).  At the shorter sampling interval  the error  was,  however,
insensitive to  the depth  increment.   When the sampling interval was reduced
to 1.5  days, the estimation  error increased  to 0.008 cm2/s  or nearly 80% of
the estimated value calculated at the 7  m depth.   Thus the  recommendation
is to  sample  at   longer  time   intervals (several   days)  and  at  finer  depth
resolution (order of 1  m).
                                    A-23

-------
                               APPENDIX B

                     Temperature equation discretization


     Temperature  equation  (l.ll  is  discretized  using  the  control  volume
method.   For intermediate control volumes (i = 2, m-1).

            A   .    k        k-fl
                       ~
            Ai
where  At is time step,  Az  is control volume width (constant)  , k  stands  for
time, i stands  for  control  volume location.  Source term  H, is described  by
equation 1.4.  System  matrix of the deterministic model Ac(k) contains terms
on the left  hand side  of equation (B.I).

     For  the surface  control  volume  (i=l) equation (B.I) will differ:   eddy
diffusivity  K 1-0.5 is  zero, the  first  entry  in the matrix is  term multiplied  by

T$+1, source terms are equations (1.2), (1.3), (1.6), (1.7), (1.8)  and (1.9). For
the  bottom  control volume (i=m) insulation  is  imposed  by  setting Ki+o-s

equal zero.  Diagonal  entry  in the  matrix is term  multiplied by Tik+1.

     Eddy  diffusivities  at   the  control   volume   interfaces  are  defined  as
harmonic  mean           .
                                       2  K.K.  .
                            Ki-o.5  = - 1-^~                        (B-2)
                                    A-24

-------
                               APPENDIX C

                     Temperature equation linearization
     Matrix Ac(k) has the same  entries as the system matrix Ac(k) given in
Appendix  B.   Matrix  B(k)  is  tridiagonal,  and  contains  derivatives  with
respect  to  lake  water  temperature at  time  step  k. For  the intermediate
control volumes :
        A-  *r   ^--n^  -k+i     5K- nK  -k+i
         i0.5       iO.o  rp    _      iQ.5  rp
         "     I  ~" XY.     ^- :  i      " i~     * :
         A.

        A...!-   dK.  ...  .k-t-l     <9K.     -k+l -,    k
         10.5 I      1-0.5  rp           1-0.5  rp       T"   -1.
        	1  	. .     1 .  .      ,     i. .    \  i .   T

           i        i                 i
         (A. .,  _ r  l/IV.   . .k+1     UJV. , _ -  .K
          1 + 0.5       1+0.5 rp	1+0.5 rp

           A   I  ffTk       i+1     dTk       '
A   L  5T        "*"     5T           J    At
  i        i                 i
        A.  -cr   5K.  n, .k+1     5K..,.k+l
         1 + 0.5]      1+0.5 ,p    _   ^ 10.5 rr,     I  rp.
         ^   L   ark      i+l     pr-k           I    i+l

where  ft = i if i =  l else ft = 0

     For the surface control volume matrix  B(k) will  slightly differ.   First,
terms  which  are  multiplied by 'Tj are the  first  entry (bn).  Secondly,  eddy
diffusivity  K j-o-s is  equal  to zero.  Thirdly,  additional terms  which take into
account  boundary conditions have to be added (Equations 1.6, 1.8, and 1.9).
These  equations  have to  be differentiated and evaluated  with respect  to the
water temperature in the first  control volume.
                            k        k         k
                      r  dHbr     5HC       dHe i  Az     k
                Had =    -    + -     +  -r   - T             (c.i)
      For the  bottom  control  volume  matrix B(k)  will  also slightly  differ.
Diagonal  term (bmm) is the one which  is multiplied  by Tj.   Eddy  diffusivity
Kio-s is  zero.
                                    A-25

-------
     Eddy  diffusivity  vector  K(z,k)  contains  epilimnion  and  hypolimnion
diffusivities.  Epilimnetic diffusivities are a function of wind speed,  thus  their
partial derivative with respect  to lake temperature  is zero  entry.  This is not
the case for the hypolimnetic eddy diffusivity.
           i                                              0P(T)  g
K = a  (P)  where N2 is Brunt-Vaisala frequency N' =(	)  -^
                                                          dz     p

                    dK . ..          K2         dK. 8N.
                    	LM = 2 	I	i -1               (C.2)
                      5T.        (K._1 + K.)2   5N. 0T.
             A
     Matrix  F(k)  contains terms which require  evaluation   of first  order
derivatives with respect  to  uncertain meteorological  variables.   Equations  1.3,
1.4, 1.6,  1.8, and 1.9  have to be  differentiated  and evaluated  at  the nominal
values  of  those variables.   Entries in the matrix are grouped  with  respect to
the perturbed meteorological parameters.  Matrix  dimensions are m x  m+3
                        [ F1(kj;
                                         ^
where air perturbation temperature vector F (k) is:
                                             Az     k
                      Hcl = (+      ' ~T/a                 (c-3)
dew  point temperature perturbation vector  F (k) is:
                                     
                                5He 0ea    Az    k
                       Hc2 = (  TT 7^ )  T'di
wind speed perturbation (mxm) matrix F,(k) is:


           ..r   5K  ._. .k+l     5K  .  .. .k+i -,      k
           21| 	Li T	i-Q.5 T     |  w ,
              L   flw.k       '-1     flw>     Xi   J  Wsi-l
Ai-^.
  Ai
Ai.o.5
  A.
                 5K  -_^^ -k+1     3K _/> -k+!  1      k
                     i-Q.5 rp	iQ.5 IT.       W '    -4-
                  n       J>   A  ~LI     ** "^.     1     I  * * S *    T
                 0Wk      1~1     5W        '    J     1
                                    A-26

-------
                                                       &HCS ) ws;
                      -T,., -    .  ruj T;    |  +  &HC, I ws;  +
        A.
                                    1-4S  -k+1 ^      k
       A.

where
                                                 T
                                                   w /
                                                 J    si
                                             Az    k
                     HC3  = (   +     )  - Ws                (C.5)
     First  an4  last  control volume have  interface diffusivities  K.     and

K.    equal to zero,  and the additional  term HC3 is present only in the first
  l"~^J.5
control volume.

                                       A
     Solar  radiation perturbation vector  F4(k) is

                                dHsn     Az    k
                        HC4 = (  j- ) - H's.                   (C.6)
                                 9H      wC
                                   A-27

-------
                                APPENDIX D

                            Cross-*enn evaluations
     Air   and  dew  point  temperature  have  significant  correlation.  The
covariance matrix between successive days for these two  parameters has been
evaluated according  to Protopapas (1988):
                    Cov [ C'(n)  C'(k) ]  = S(n) Mc S(k)T               (D.I)
where
  S(n) =
<7taln)
0
0
0
0
<7td(n
0
0
0
) 
0
0
0
0
0
0
                Aa(n,k)       Ptatd(n,k)      0         0

                              />td(n,k)       0         0
   LC 
                00             00

                00             00
  S(k) =
0u(k)
0
0
0
0
tftd(k
0
0
0
) 
0
0
0
0
0
0
where  Aa  is air temperature standard deviation,  k) is  air  temperature correlation between  two  day
Ptd(n>k) is  dew point temperature correlation,  ptatd(n,k) is  air temperature
dew point  temperature correlation.

     Replacing index  n  by index k  leads  to  the  form  of the  covari;u:<-
matrix  M(k,k) for the zero time lag.   In addition, the diagonal terms an- ^
equal to one  in the correlation matrix  Mc.  They are equal  to the  standt-
deviations  of  air  temperature  ffat(k),  dew  point temperature  ortd(k),  *--
speed o"ws(k),  and  solar  radiation <7sr(k)  in  the  matrices S(k), and S(n).
                                     A-28

-------
      If  meteorological perturbations are not cross-correlated covariance matrix
Muc(k,k) is

                     
-------
               APPENDIX E





Regional lake water temperature simulation model
                   A-30

-------
                        LAKE INPUT DATA FILE
TITLE
SEKI
NDAYS
NPRNT
DZLL DZUL BETA EMISS WCOEF WSTR
WCHANL WLAKE DBL ST S FT
ELCB ALPHA BW

Initial conditions^
Secchi depth reading
Number of days for output
Dates for output
Heat budget and mixing coefficient
Initial stage & Outflow channel
MBOT NM NPRINT INFLOW DY MONTH ISTART MYEAR
Z(1)...Z(MBOT)
T(1)...T(MBOT)
Field data section

NFNPRFLE'
NFLD
DEPTH(1)...DEPTH(NF)
FD ATA( 1 )...FD ATA(NF)
Number of depths & parameters
Parameter code (1)
Depths
Temperatures
NDAYS
                                      A-31

-------
                           EXAMPLE INPUT DATA. FILE

 LAKE CALHOUN 1971  ** from APRIL through december **
2
9
520 608 629 719 802 824 913  1011  1028
0.15 1.00 .4 .97 24.5 0.4
100. 100. 200 224.0 .001 .035
205 18 100
24 8 1 0 91 4 1 1971
0.5  1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5
12.5 13.5 14.5 15.5  16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5
4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4, 4. 4. 4.
11 1
1
0.   2.  4.  6.  7.  8.  10. 12T15.  20. 23.
13.1 12.8 12.5 12.5  12.5 10.4 8.6 7.9 7.4 7.  6.9
14 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12.  15. 20. 25.
20.4 20.3 20.1 19.8  15.1 13.6 113 11.7 11. 9.2  8.1 7.3 7.1 6.9
12 1
1
0. 1. 2. 3. 4. 5. 6. 8. 10. 12. 15. 20.
24.4 24.6 24.6 24.3  23.5 19. 15.1 11.5  10. 8.8 7.2 7.
141
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12.  15. 20. 24.
23.1 23.1 23.1 23. 22.9 22.8 17.6 12.4  11.4 10.2 9. 8.2 7.8 7.7
14 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12.  15. 20. 24.
21.  21. 20.8 20.8 20.8 20.6  19.4 14.7 11.3 9.8 9. 8.2 7.5 7.4
11 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12.  15. 20.
22.7 22.7 22.7 22.7  22.6 22.6 20.5 15.7 12. 10.2 9.2 8.1 7.7
15 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 12. 15. 20. 25.
21.7 21.7 21.5 21.5  21.5 21.3 21.1 18.5 13.3  10.6 10. 9. 7:9 7.6 7.5
13 1
1
0. 1. 2. 4. 6. 8. 10. 12. 14. 16.  18. 20. 22.
14.5 14.5 14.3 14.3  14.2 14.2 14.2 11.1 10.1  9.3 8.6 8.4 8.3
71                        '
1
0. 12. 13. 14. 15. 20. 23.
117 117 10.5 9.0 8.4 7.9 7.8
                                                 A-32

-------
                       LAKE INPUT DATA VARIABLES
SEKI
NDAYS
NPRNT
DZLL
DZUL
BETA
EMISS
WCOEF
WSTR
WCHANL
WLAKE
DEL
ST
S
FT
ELCB
ALPHA
BW
MBOT
NM
NPRINT
INFLOW
DY
MONTH
ISTART
MYEAR
Z
T
NF
NPRFLE
NFLD

DEPTH
FDATA
Secchi depth reading (m)
Number of particular dates selected for output
Dates selected for output
Minimum layer thickness (m)
Maximum layer thickness (m)
Surface absorption factor for solar radiation
Eraissivity of water
Wind coefficient for convective heat loss
Wind sheltering coefficient
Width of the inlet channel (m)
Width of lake  (m)
Elevation of the bed in the bottom layer (m)
Stage of the lake on the first day of simulation (m)
Bed slope at inlet channel
Roughness coefficient at inlet channel
Elevation of the bottom of the outlet channel (m)
Side  slope of outlet channel
Bottom width  of outlet channel (m)
Total number  of layers in the initial conditions
Number of months to be simulated
Interval between days for tabular output
Number of inflow and outflow sources
Julian day of first day of simulation
First month of simulation
Day of the month that simulation starts
Year of simulation
Depths in the  initial conditions  (m)
Temperatures  in the initial conditions (C)
Number of depths in field data profile
Number of profiles in the field data
Parameter code to match field data profiles to state variables
1 = Temperature (C)
Field data depths (m)
Field data temperatures (C)
                                             A-33

-------
                 METEOROLOGICAL DATA FILE
MONTH KDAYS MYEAR
TAIR TDEW WIND DRCT RAD CR PR
KDAYS

where

MONTH
KDAYS
MYEAR
TAIR
TDEW
WIND
DRCT
CR
PR
Month
Total number of days in the month
Year
Average air temperature (oF)
Dew point temperature (oF)
Average wind speed (mph)
Wind direction
Percent of sunshine
Precipitation (inches)
           EXAMPLE METEOROLOGICAL DATA FILE
4
32
21.5
26
32.5
34
39.5
53
55
49
54
54
50
39
40
55
61
55
52.5
30
22
10
14
12
14
23
29
32
27
28
34
25
20
23
32
46
36
37
1971
22.4
19.3
13.9
7.6
3.9
6.7
9.3
13.5
13
15.3
13.5
6.8
9.7
4.3
13.9
12.9
6.6
11.7

320
320
310
280
310
30
80
50
290
80
340
280
280
50
100
40
190
140

203
324
514
568
558
489
533
398
553
527
432
373
389
471
531
347
565
76

9
8
93
100
100
100
100
94
100
100
62
- 78
73
79
99
79
97
7

13
6
0
0
0
0
0
0
3
0
0
0
1
0
0
2
0
11
           KDAYS
                                  A-34

-------
      1U4W789
      SLARGE
         PROGRAM REGIONAL

      C
      C     THIS PROGRAM IS MODIFIED VERSION OF THE MINLAKE
      C     MODEL (JULEY * STEFAN, 1987. UNIVERSITY OF MINNESOTA
      C     SAFHL, PROJECT REPORT * U3\ THE COMPUTER CODE
      C     IS ADDAPTED FOR THE REGIONAL LAKE WATER TEMPERATURE
      C     SIMULATION IN A LAKE ATTATCH USER SUBROUTINE DURING
      C     UNKING.
      C
         COMMON/MTHD/TAIR(31).TDEW(31).RAD(31),CR(31).WIND(31).
         + PR(31).DRCT(3I)
         CX)MMON/RESULT.T2('W),CHLATOT()).PA2(),PTSUM(W).BODy40).
          COMMON/FLOW/HMK(41).OE(40),FVCHLA{5),PE(5.41)
          COMMON/VObJZMAX.DZ(-).Z<40)A(l).V(W),TV(W)J\TOP(41)J)BL
          COMMON/SURSDZ{>).SZ<>),LAY(0).AVGt(<.60).SVOL(0)
          COMMON/CHOICE-MODEL.NrrRO.IPRNT(6XNDAYS_STRNT(30VNCLASS.PLT(30)
          COMMON/WATER.BETA.EMISSJCKljaClHKMAX.WCOEF.'WSTR
          COMVON/CHANE^'WCHANUELCaALPHA.BW.WLAKE
.          COMMON/STEPS:DZU_DZUUMBOT.NM.NPRINT.MDAY.MONTKILAY.DY
          OOMMON/STAT,SUMXY(10).SUMX(10)5UMY(10)J{SQ(10).
          +YSOX10).RSOC10).RMS(10).RELM(10),MTHREL(10).MDAYREUlO).ZREmO).
          * ZRELM( 10).RS( 10XREU lOJ.MTH RMS( 10).MDA YRMS( 10).ZRS( 10).2RMS( 10)
          COMMONnNFi.OU,OIN(5).TlN(S),PAlN(S).BODIN(5).DOIN(5).CiN;S).
          +CDIN(5)_XNH!S(5;.KNOIS(5).CHLA1N(3J)
          COMMONOURCRM(3.JO).PROD(40).XMR(3.*)),PRODSUM(-!0)
          COMMON/FIELD/ !FLAG(10).FLDATA(ia50).DEPTHC50),NFLD{10)
          COMMON/FILE' DIN.MET.FLO.TAPE&TAPEUREC
          COMMONmn.' TTTLE
          DIMENSION B(10)5TATVAR(60VICX(-))
          CHARACTER'04 17TLE
          CHARACTER-16 DlN.METJLO.TAPE&TAPEl
          CHARACrER'I TS
-------
 505 FORMATC REQUEST CHANGE OF TITLE (Y/N) T.SXJ)
 506 FORMATC ENTER NEW TITLE .-)
    WRTrec.505)
    READC.XA)1) XS
   IF(XSEQ.-r .OR. XS.EQ.V) THEN
    WR1TEC.SOIS)
    READC.-(A)-) TTTLE
   ENDIF  '
   WRJTE(8,999) FF
   WRrTE(8,1900)
 999 FORMAT{LXA1)
 1001 FORMATC ENTER INPUT DATA FILE NAME :'.10X/)
 lOO: FORMATC ENTER RLE SAME FOR TABULAR OUTPUT:  '-0
 1000 FORMATC ENTER METEOROLOGICAL DATA FILE NAME : Vj
 1003 FORMATf ENTER INFLOW DATA FILE NAME :'.9XJ>
 1900 FORMATC//.8XO
C
    CALLSTART(STAFT.lSTART.INFLO'.MYEAR.IRUN.ILSEn)
C
C"'" Call LAKE routine u> change or >dJ input quanliba 

   ZMAX=ST-DBL
   CALL SETUP

  11 IF(DZ(I).GT.D2UI_AND.MBOT.LT.W! THEN
    CALL SPUT(ULAY)
    GO TO 11
   ENDIF
   l-I+l
   IF(LGT.MBOT.OR.I.GT.40) GOTO 12
   GOTOn
  12 CALL SET2(MBOT)
   CALL VOLUME(MBOT)
   CALL AREA
   CALL ABOUND
   CALLTVOL(MBOT)
C_DETERMINATION OF INITIAL MLXED LA'i'ER DEPTH
   1LAY-1
   DO7I-1.MBOT-1
   IF(T2(I).GT.T2(I + 1)+CONST) GO TO S
  1 ILAYILAY1
  8 DMLX=Z(lLA^-)+DZ(ILAY)-aS
   DO 9 1=1.10
    RLM(l)a
  9  RMS(I)=aO
   MP.O
   IPRNT(1)1
   KDAYS=0
   IPRNT(5)=IPRXT(5)-1
   MDAY=0
   CALLPTABLE
   1PRNT(5) = IPRNT(5) * 1
C~ Sun simubiion for each montb
   EDA Y= 345.
   YEARREAL(MYE.\R)
   IF(AMOD(YEAR.4.0).EO.aO) EDAY-.Vid
   ETSUM=a
   HT5UM-a
C
   DO100JLNM
   CALL MTHDATA(MONTH.KDAYSJrfYEAR)
   EDAY-345.
   YEAR=REAUMYEAR)
   IF(AMOD(YEAR.4.0).EO.a) EDAY3oi
C-START SIMULATION FOR EACH DAY
   DO 200 MDAY-ISTART.KDAYS
   NFLOWolNFLOW
   CALL LAKE
-------
      c
      C-DETERMINATION OF WIND MIXING ORDER
      C_HEAT IS ABSORBED FIRST, THEN WATER COLUMN IS MIXED BY THE WIND
      C
         CALLHEBUG(ILAY.TMIXQNET,HS.HA,HBR,KE.HC.
         +TAIR(MDAY),TDEW(MDAY),CR(MDAY).RAD(MDAY).WIND(MDAY),
         HPAN.PCOEFF.SEKI)
         CALL CONMlX(ILAY.TMrXMBOT)
      C-.CALCULATION OF EVAP. IN TERMS OF VOLUME
      C-CALCULATES LATENT HEAT OF VAPORIZATION ALV
         HEDHE/((597Jl*531'n(l))'RHOCn(l),C2(l),CD2(l)))
         HEVHED*ATOP(l)
         CALL WINM1X(RKE,TMDULAY.MBOT)
         DMIXZ(lLAY).EQ.NPRNT(NDAYS))THEN
      C  Access and output field data
           CALL FDATA(NF)
           WRTTEC.3020)
           READ(*,99) XS
       99   FORMAT(A1)
      3020  FOR-MATC///: DO YOU WANT TO MEW GRAPHICAL RESULTS: ( Y/N)'.
         +    3)
      C_.Ca[l plouing routines
           IF(XS.EQ/Y' .OR. XS.EQ.VT THEN
            CALL SUBPLOT(NF.MYEAR)
           ENDIF
           ENDIF
         IPRNT{1)=0
       200 CONTINUE
         ISTART=1
       100 CONTINUE
      C_CotBputc 2nd list statistics!
      C_   l}relatjve am) absolute mnimmn deviacioRf ben>cen
      C    model and Geld data and day of occurrence
      C   2) slope of regression of Held dau on simulation results
      C-   3) regression coefticiem
      C   4) standard error of the regression
         WRJTE(8,3000)
         DO 101 1=1.10
          ff(XSO(I).GT. 0.) B(I)=SUMXY(I)/XSQ(I)
          IF(IFLAG(I).GT.2) THEN
           SUMXY(l)=(YSQ(l)-SUMXY(l)'SUMXY(!)/XSQ(l))/(IFLAG(l)-2)
           RSCKI)-L-SUMXY(I)'
         +   IFLAGCO-OFLAGt^iydFLAG
           SUMXY(I)=SQRT(SUMXY(I))
          ENDIF
       101 CONTINUE
         WRTTE,7X'rT,X,IBOD-.7X.-DO1.
         +X.TSS'.XTDS-.6X'N03'.OC,'NH4'y)
       3013 FORMAT(LX,'MAXIMUM RELATIVE ERROR ^MOXlOtSXF-t.O))
       3014 FOR.MAT(1XTJATE OF MAX REL, ERR,M.10(4XI2.r.I2))
       3015 FORMATCIX^MAXIMUM ABSOLUTE ERRORM4X10(2XF7J))
                             A-37

-------
    301* FORMAT(LVDAT OF MAX ABS. ERR.M.10(XSQOO),
       +YSQ(10).RSQ{10).RMS<10).RLM(10)_MTHREmO).MDAYREUIO).ZREmO).
       +ZRELM(10).RS(10XREU10).MTHR.MS{10XMDAYRMS(IO)JZRS(10).ZRMS(10)
       COMMON/VOUZMAXDZ<40).Z<40XA(40).V(4).TV(40)ATOP(41).DBL
       COMMON/STEPSCZlJ.DZin.MSarNM.NPRirn'.MDAy.MONTHILAY.DY
       COMMON/CHOICEMODEL.NJTRO.IPRNT(<).NDAYS,NPR,ST(>0),Naj^SS,PIjOT(30)
       COMMON/FIELD/ IFLAOCIO).FLDATA(:a50).DEPTH(50).NFLD(10)
       DIMENSION COMP(40.10)
'      EQUIVALENCE (T2(1).COMP( 1.1))
       CHARACTER'16 DIN.MET.F1X).TAPES.TAPE1
       DO303IetlO
        RS(1)=0
    303  RELd)-0
       DO304je!.20
        DO304I = UO
    304  FU)ATA(U)=0.0
       READ(?.*)SF.NPRFLE
       NDAYS-NDAYS*!
       IFCNF.GT.O) THEN
        READ (7.')(NFUX1).1=1.NPRFLE)
        READ(7.-)(DETrH(l).I-= LNF)
        DO30SI-J.NPRFLE
         READ(7.-)(FLDATAf NFLD(I)J)J - 1.NF)
    305   CONTINUE
        U,l
    C Locate simulation values corresponding to sampled
    Cconsu'tuenu and depth of field data
        DO 310 KS= 1.NF
         L=LL
         DO 315 LL
-------
    WRrrE(8>3015) (RS(l).l-UO)
    WRrTE(8,3016) (ZRS(I).I=U10)
C_Siore dau on plot file (Upe&PLT)
    IFOPRNT(5).GT.O) THEN
      WRTTE(1.REC=IREC) REAL(NF)
      IREC-IREC+1
      WRTTE(tREC=IREC) REAUNPRJFLE)
      [REOIREC + 1
      DO5001-LNF
       WRTTEdREColREC) DEPTH(t)
 500     IREC-IRE01
      DOJ011-1.NPRFLE
       WRTTEaREC-IREC) REAL(NFLD{I))
 SOI     IREC=IREC+1
      DOS0212-1.NPRFLE
      DO5021-UNF
       WRrTE(l.REC-!REC) FLDATA(NFLD(t2).l)
 502     IRECIREC+1
      ff( NFLD(l).NE.l) THEN
        X-0.0
      ELSE
C-Mixcd bycr depib in Held diu ukcn ai dT/dZ- 1.0
       DO 503 J2=1NF
        X=(FLDATA(U:-D-FLDATA(U2)V
   +                     (DEPTH(J2)-DEPTH(JM))
        IF(XGT.LO) GOTO 504
 503      CONTINUE
 504     X=(DEPTH(J2)*DEPTH(J2-1))'0.5
       IF(JZGENF)X-ZMA.X
      END1F
      WR1TE(1.REC=IREO X
      IREC=iRECl
    ENDIF
   ELSE
    NF=0
    IF(IPRNT(5).GT.O) THEN
      WRrTE(l.REC=lREC) REAUNF)
      KEC=IREC+1
    ENDIF
   ENDIF
 29 FORMAT(1XIXI5)
 3009 FORMAT(///5X,'DATE: \if f ;aj/SX.
   +SUM OF ERRORS BETWEEN MODEL AND FIELD DATA1)
   +-LAKE TEMPERATURE STATISTICS'/)
 3013 FORMATfUfMAXIMUM RELATIVE ERROR (%)MOX.10(5XF4.0))
 3014 FORMATCIXUEPTH OF MAX REL. ERR. (M)'.9X10(5XF4.1))
 3015 FORMATWMAXIMUM ABSOLUTE ERROR'.14X10(2XF7J))
 30W FORMAT(1XT)EPTH OF MAXABS. ERR. (M)'.9X10(5XF4,1V//)
 3020 FORMAT(/42XTEMP'.iX1CHLA>.7XTP-.7XTXn
   RETURN
   END
 C
   SUBROUTINE PTABLE
 C"***
 C"'" Send simulation raulu to the tabular
 C" output file (upeSDAT)
 C*****
   COMMON/RESULT/ T2<40).CHLATOT()).PA2<40).PTSi;M(40).BOD2<40),
   +DSO2(40).C2(40),CD2(40)JCNO2{40)JCNH2(40).CHLA2(3.40),
   +PC2(3,40)JCNC2(3.40).T20(40),SI2(40)
   COMMON/VOL/ZMAXDZ440).Z(40) A(40).V(40).TV(40) j\TOP(41 )JDBL
   COMMON/STEPS/'DZLL.DZUL.MBOT.NM.NPRINTJUDAY.MONTHILAY.DY
   COMMON,'CHOlCE'MODEL.NrrRO,IPRNT().NDAYS,NPRNT(30)J
-------
     IF(NrTRO.GT.O) THEN
     ELSE
     ENDIF
     WRITE(a3006)
     IFCNITRO.Crr.O) THEN
     ELSE
     EVD;F
   ENDIF   .
 1200 FORMAT(///5X'ZOOPLANCTON PARAMETERSV-SX^v
   + 'CONG (*.3X'GRAZING(MG/M3)-.3X
   * ' DAYDEPTH(M)V.9Xn.0.9XF7.1.9XF&4.1IXF4.1.'/)
 1201 FORMAT(Fii.}                       ..   . ..
 399 FORMAT(//5XlNmAL CONDmONSVSX18(Vy) .'
 3000 FORMAT(//5X1NFORMATION ON LAKE QUAUTYVSX27r>V)
 3001 FORMAT^X'ZfM)  T(C) SS(PPM) TDS(PPM) CHLA(PPM) PC(PPM)',
   + ' P(PPM) PT(PPM) BODCMG/L) DO(MG/L) DZ(M)'.
   + ' AREA(M2)   VOUWJIT
 3002FORMAT(5X3(F5.2.2X).lXFS.J,4X4<2XFSXF<2.
 3006 FORMATCWX-BIOLOGtCAL PARAMETERSV5X21C-VJX'Z(M)-.X
   +-CHLOROPHYLL'.XTOTAL-.6X1NERNAL PM1X1NTERNAL SV16X
   +-V.XT.6X-3'JX1CHLA-.4XT.X7-.X'31.6XT.XT.
 3007 FORMAT(5XFJ.2.2X10(F5J.2X))
 3008 FORMAT(//.5X' TEMPERATURE PROFILESV.5X'Z(M)-.
p****
C**"** Prtxiucc profile of slate variables and
** field data (if niilable)
^
   CHARACTER'1 XS.CH1.CH2,CH3.CH4
   CHARACTER'S; TITLEO)
   CHARACTER'4 MNTH(15
   COMMON/SOURCE/RM(3,40).PROD(40)JCMR(3.40)J'RODSUM()
   COMMON/STEPS'DZLL.DZUUMBOT.NM.NPRINT^IDAY.MONTHILAY.DY
   COMMON/RESULT/ T2(40).CHLATOT(40).PA2(40).FrSUM(40).BOD^-IO),
   +DS02(40).C2(40).CD2(40-^XNO2(40).XNH2(40XCHLA2(3..
   COMMON/F1ELD/1FLAG(10).FLDATA(10.50).DEPTH{50)>.'FLD(1I')
   COMMON/\'OL.'ZMAXDZ(40).Z())A(')).V(40Vn'(-<0)j\TOP(41).DBt
   DIMENSION ZD(23).FD(23).\TM(43).ZV(43).VAR().IO)
   DIMENSION FCT(10).LEN(10).1SCAL(10).ICX(4)
   EQUIVALENCE (T2(1XVAR(1.1.(2D(1).PROD(1)).(FD(1).PRODSUM(I)).
   EQUIVALENCE (CHUCX(1)MCH1ICK(2)).(CH3.ICX{3)).(CH4.IC.\(4))
   DATA lCX/16*lB.16*5B.iw3il6-4A
   DATA BCAL/UMJ-U'-I/
   DATA Fcr/u*iooo,4'i.ri(x
   DATA MNTWJAN '.'FEB '.TvlAR '.'APR '.TklAY VJUS VJUL '.'AUG '.
   + 'SEP '.XXTT ,'NOV '.T)EC V
   DATA TTTLE /TEMPERATURE (QV
   DATA LEN /1^27.2&25.23.<-24.2S/
 1000 WRTTEC.109) CHLCH2.CH3.CH4
C-IUl and idea desired suie variable for plotting
   DO100I=1.1
 100  U-RITE{M01) l.TrTLE
-------
 110 CONTINUE
   12-0
C_if Geld data b nibble, locale Held data corresponding
C-to selected state variables
   IF(NF.CT.O) THEN
    DO 111 1-tNF
    IF(FLDATA(ICLI).GT.O) THEN
     12-12+1
   .  FD02)=FU>ATA(ICU)
     ZD(I2)--DEPTH(I)
     IF(X,LT.FD(12)) THEN
      X-FD(12)
      IND 1
     END1F
    ENDIF
 111  CONTINUE
    ENDIF
    NPEN-1
    NPEN1-1
    IPORT-99
    MODL=W
    ZV(MBOT+1)-0.
    ZD(I2+1)=0.
    VTM(MBOT+lj=0.
    FD(I2+1)=0.
    FCTOR.-0.7I
 C_begin plotting sequence
 990 CAUL PLOTS(O.IPORT,MODL)
    CA1JLSIMPLX
    CALL FACTOR(FCTOR)
    CALL NEWPEN(NPEN)
 C-dflenninc maximum i &. y values either in simulated
 C atate variables or in fieid data
    IF(IND.GT.O) THEN
      13=MBOT+1
      CALL SCALE(VTM,10,13. 1)
      DX=VTM(MBOT+3)
    ELSE
      13-12+1
      CALL SCALE(FD.10.,13.1)
      DX=FD(U+3)
    ENDIF
    lF(ZV(MBOT)rr.ZD(I2)) THEN
      D-MBOT+1
      CALL SCALE(ZV.i,I3.-l)
      YSC ZV(MBOT+3)
     ELSE
      D=I2+1
      CALL SCALECZD.4,D.-1)
      YSC ZEKI2+3)
     ENDIF
     ZV(MBOT+2)YSC
     ZD(I2+2)-YSC
     VTM(MBOT+2)-DX
     FD(I2+2)-DX
     DAX1S=DX*FCT(IC1)
     CALL STAXIS(i.I5.2,15.ISCAL(ICl})
  C_draw x-axis
     CALL AXlS(U,TTTLEaCl)aN(ICl).-10,a.O,DAXJS)
     CALL STAXB(.2,25_2_li,l)
  C-dnw y-axis
     CALL AX1S(U1,T)EPTH  (M)MO.-^9a.YD.-YSC)
     XMDAY=MDAY
     XMYEAR-MYEAR.
  C_prim title to diagram
     CALL SYMBOL(4J J..25.MNTH(MONTH),0,4)
     CALL NUMBER(999,999_2S.XMDAY.O.-1)
     CALL SYMBOU9M,999_2S.-. '.0. 2)
     CALL NUMBER(W,999_25.XMYEAR_0,-1)
     CALL NEWPEN(NPENl)
     CALL PLOT (U7,-3)
  Cplot simulated profiles as a line
     CALL UNE(VTM.ZV,MBOT.L
-------
     WRTTELO ) : 'A)
      GOTO 990
     ENDff
     GOTO 1000
  1001 RETURN                          '    .
     END
     SUBROUTINE ABOUND
 C
 c
 c
 C*"
 c
 c
     COMMON.'STEPS.'DZLUDZUl-MBOT.NM.NPRINT.MDAY.MONTH.Il-AYDY
     COMMOS.VOL/2MAXDZ<40).Z( 40)^(40). V(40).TV(40),ATOP(41).DBL
     DUM-tt
     DO10011.MBOT
/    ZDUM=ZMAX-DUM
     DUM=DUM+DZ(I)
     CALL LAKE(ZDUMw\DUM0.1)
  100 ATOP(1)=ADUM
     ATOP(MBOT+l)-a
     RETURN
     END
     SUBROUTINE AREA
 c*****
 C""" Compute ihe area through the middle of each layer
 C***** using the depth-area relationship in LAKE
 C*****
     COMMO.VSrEPS^>2U.D2UL.MBOT.NM.NPRJNT.MDAY.MONTHJlJ\Yjy
     COMMONA'OlJ2MAX.DZ<40).Z(40)j\C40).V(0).TV(JO).ATOP(41).D3L
     DUM>a
     DO1001.1.MBOT
     ZDUM-ZMAX-DUM-DZ(!)/1
   CALL LAKE(ZDUMADUM,0,1)
 100 A(I)ADUM
   RETURN
   END
    SUBROUTINE COEF(MODELMBOT.NCLASS)
^
C*""** Compute some coefficients used in the constant
C"**** volume and finite difference solutions
C*****
     COMMON/COEFF/ DUM^40).DUM3(40)
     COMMON/VOUZ.MAXDZ(40).Z(40)JV(40).V(40).TV(40).ATOP(41).DBL
     COMMON/FLOW/ HMK(J!).OE(40).FVCHLAr5).PE<5.41)
    DOlOOi.lMBOT-1
 100   DUM^D^DUMI'ATOPd +
   KK2
   IF(MODELGT.l) KK=1
   DO 300 K=KK.NCLASS+1
     DOaWi-iMBOT
     XFVCHLA(Kr(DZ(I-l)*D?(l))'J/HMK(I)
C - PE-tLO-.rABSfX))".^ -
     AO.L0..1-ABS(X)
     A1=AO-AO
     PEfK.I)=AI'ArA(VX
   PE(K.l)-aO
 200  PE(KJHBOT+1)0.0
   RETURN
   END
   SUBROUTINE CONMLX(ILAY.TMIX.MBOT)

C***** Remo\e deniity insubilities by mhang umublc
C""* \rrcn dcwnuranl and mcrjjngunth lower layen.
^*
   COMMON RESULT/ 17(40).CHLATOT{40).PA2(40).PTSUM(40).BOD^40).
   + DSO).C2(40).CD2(40).XNO2(40)JCNH2(40).CHLA2(3.40).
   +PO(3,40).XNC2(3.40).TIO(40).SI^40)
   COMMON VOL/ZMAXDZ(40).Z<40)JV(40).V(40).TV(40)J\TOP(41).DBL
                       A-42

-------
           DIMENSION RHOT(40)
           DO100I-UMBOT
         : RHOT(I)-RHO(T2(l).a.O.)
         6 1FLAG-0
           1-0
           M-MBOT-l
         1 1-1*1
           IF(LEO.M) GO TO 5 .
           IF(RHOT(l).LE,RHOT(ltl))GOTO 1
           IFLAG-1
           IB-I
           TVDUM-TI)'V(I)
           VDUM-V(I)
           TDUM=TVDUM/VDUM
           RHODUM=RHO(TDUM.a.O.)
           J=!-l
         3 J=J+1
           IF(RHODUM.LE.RHOT(J + 1)) GO TO :
           1FLAG.1
            VDUMVDUM+V(J + 1)
            TDUM-TVDUM.'VDU.V!
            RHODUM = RHO(TDUM.O..a;
x           IF(J.EO.M) OO TO I
            GO TO 3
         2  IE-J
            !F(J.EO.M)1E"IE + 1
            IF(iaEO.IE) GO TO 4
            T2(K)=TDUM
         y RHOT(K) = RHODUM
         4  IF(I.SE.M)GOTO1
         5  IF(IFLAG.EQ.1)GOTO6
        C-DETERMINE MLXED LAYER DEPTH...
            DO700I = 1.MBOT-1
            IF(CnWOU2.MAXDZ<40)i(40)J\(40).V(40).TVC40)w\TOP(41).DBL
            COMMONi'RESUUT/T2('W),CHLA.TOT(40).PA2(40).fTSUM(40).BOD2(40),
           +PO(3.40)JCNC(3.40).Ta)(40)^C(40)
           COMMONOiOIC/MODEL,NrrRO,]PRNT(6).NDAYS.NPRNT(30),NCUASS.PLOT(30)
           ff(V(I)J-E.O.) THEN
              !F(LEQ.M8OT) THEN
                D-MBOT
                MBOTMBOT-t
                IF(V(II-l).LE,aO) THEN
                 II-IM
                 GOTO:
                ENDIF
                RETURN
              ENDIF
              V(Itl)-V(I)V(l + l)
              DZZD2(0
              KK=!
              MBOT=MBOT-1      >
            ELSE
             U=I
             !F(LEQ.MBOT) 11=1-1
             KK1! + 1
             VC=V(IO+V(KK)
             VCOMB=UVC
             CD2(n)=(CD^SI)*V(l!)+CDZ(KK)'V(KK))*VCOMB
             DOS5K=1.NCLASS
         55    CHLA2(K.II)(CHLAKK.n)-V(II)+CHLA2(K.KK)%V(KK))>VCOMB
                               A-43

-------
         IF(MODELEOJ) THEN
         DOS7K-1.NCLASS
   57     PC2(O)(PC2(K.BrV(7I
         ENDIF
        PA2(H)=CPA2(II)'V(ir)+PA2(KK)'\'(KK))'VCOMB
        BOD^I1)={BODZ(II)*V(II)+BOD^KK)*V(KK))'VCOMB
        DSO2(II)=(DS02(IU'V(II)+DSO2(KK)'V(KK)}'VCOMB
         IF(NITRO.E0.1) THEN
          DOSK=1.NCLASS
   56      XNC2(K.n)=CXNC2(K.II)>V(Il)+XNCKK.KK)tV(KK))'VCOMB
         ENDIF
        v(ii)=vc
        DZ(1I)=DZJBOT
      T2(K)-T2(K+1)
    .  0).
      * DSO2()).a(40).CDI(40)_XNO^))_XNH2(t.CHLA^3.<0).
      COMMON/TEMP/PARIO(4).PCDUM(3.40)JCVHD(VXNOD(40).
      + CHLADUM(3,40)JCNCD(3,40).PADUM(-W)51D(40)
      COMMON/STEPaDZU.DZUl-MBOT.NM.NPRlNT.MDAY.MONTH.ILAY.DY
      COMMON/FLOW/HMK(41XOE()JVCHLA(5).PE(5.41)
      COMMON/WATER/BETA.EMISSJCK.tJCKZHKMAX.WCOEF.WSTR
      COMMON/SOLV/ AK(*0).Bm40),CK<40).DtC(<0)
      DIMENSION CK40)
   C.-.CALCULATION OF THE HEAT ABSORPTION FROM METEOROLOGICAL PARAMETERS
   C_IN A COLUMN OF WATER
   C-.CALCULATION OF HEAT FLUXES INTO THE WATER BODY
      CALL FLXIN(HS,HA,TAIR,RAD.CR.C2)
      CALLFLXOUT(TS.HBR.HE,HC,TA1R.TDEW.W1ND.IPAN.PCOEFF.WCOEF)
      HOOUT-HBRHE+HC
   C..CALCULATJON OF EXTINCTION COEFF. (ETA) AS A FUNCTION OF SUSPENDED
   C-SED1MENT CONCENTRATION
       ETA=l.W(l/SEia)
   C.CALCULATION OF HEAT ABSORBED IN EACH LAYER
        HO=(L-BETA)'HS
      EX=EXP(-ETA'DZ(1))
                          A-44

-------
 C-CONVERSION FACTOR OF 1000 USED FOR DENSfTY-HEAT CAPACITY OF WATER.
    HQ-HO/EX
 C-CALCULATE THE SOURCE TERM Q FOR EACH LAYER
    D010U2.MBOT
     EX-EXP(-ETA'DZ(I))
     ODZ.:ir(DZ
-------
   HE=FCN'(ESA-EA)
   OOTQ30
 20  HE=15-lTD'PCOEFF*59a'ia
C-.CALCULATES BACK RADIATION
 30  HBR=*(U71E-*0.7TSK.TSKTSKTSK)
C_CALCULATES CONDUCTIVE LOSS USING BOWENS RATIO
   HC=FCN-0.617W(TSK-TAK)
   RETURN
   END "
   SUBROUTINE SETAMKCW.HKMAX11AY.MBOT)
C"
 ' Compute venial diffusion coefficient in each layer.
 ** Diffusion coefficient berween layers as tbe bannonic
 1 mean of tbe diffusion coeffidenu in adjacent layer*.
 
DIMENSION AMK(41)
COMMON/FLOW/HMK(41).OE(40).FVCHLA(5).PE{5.41)
COMMON/VOU2.MAXDZ(40).Z<40)A(40).V(40).TV(40)ATOP(41).DBL
COMMON/RESULT/T2(40),CHLATOT(40).PA2(40).PTSUM(40).BOD2(40X
C"
C"
c-
C"
   +PC2a0)_XNC2(3.40).T20(40).SI2<40)
   DIMENSION PSON(40)
   AREA-(ATOP(iy(10"0)
   ALFA=ai7*O.OOOr(AREA)"0.5
CVenical diffusion coefficient in tbe mixed layer
C computed as a function of tbe wind speed
   DKM2S.*WL3
   IF(DKM.GT.HKMAX) THEN
   SON =0.000075
   DKM.(ALFA/((SON)"0.43)r&
   ENDIF
C
   D0100I=ULAY
 100 CONTINUE
C
C 929 FORMATtlXHYPOUMNION1)
C_RMK TEMPORARILY USED TO STORE DENSITIES.
C_diffu>ion coefficient below tbe miied layer computed as
C_a function of the square of tbe Brunt-Vaula frequency (SON)
C_ SON = (G/RHO) ' d(RHO)/dZ
   DO 200 1=ILAY.MBOT
 200 HMK(I)-RHO(T2(I).C2(1).CDI(I))
C   WRTrE(99.)
   DO 300 lcILAY + l.MBOT-1
   SON = ABS(HMK(I-1)-HM1C(I + 1))/((Z
-------
20  CONTINUE
   DO 1WK-LKDAYS
   READ(9/)TA]R(K).TDSW(K).WIND(K).DRCT(K).RAD(K).CR(K),PR(K)
   PRdC)-PR(KyiOa
104 CONTINUE
   DO 100Kl,KDAYS
    TAIR(K)=(TAIR(K)-31)'0.5S56
    TDEW(K)-(TDEW(KV31)-aS55<>
100  CR(K)-(100.-CR(K))W.01
C
C MAKE CORRECTIONS FOR CLIMATE MODEL OUTPUT
C IF NEEDED. THESE ARE MONTHLY ADJUSTMENTS.
C
C   DO1287K=1.KDAYS
C   IF(MONTH.EQ.3) THEN
C   TAIR(K)=TA!R(K)+4.SO
C   TDEW(K)=TDEW(K)7.25
C   WIND(K)=WIND(K)-aS2
C   RAD(K)=RAD(K)>L04
C   PR(K)PR(K)*1-02
C   ENDIF
C   !F(MONTH.EO.4) THEN
C   TAIR(K)=TAiR(K1 + 497
C   TDEW(K)=TDEW(Kl-i.40
C   WlNDPR(Kri.!4
C   ENDIF
C   IF(MONTKEO.ll) THEN
C   TAIR(K)=TAIR(K)+5.93
C   TDEW(K)=TDEW(K)5.51
C  WIND(K)WIND(K)'1.0S
C   RAD(K) = RAD(K)'0.95
C   PR(K) = PR(K)'L1
C   ENDIF
C 12S7 CONTINUE
 1001 FORMAT(//.SX.70f')^.MX'PROGRAM ABORTED.V.10.X.
   + 'METEOROLOGICAL DATA FILE DOES NOTJ.15X.
   + MATCH YEAR OF SIMULAT.ONV.'.5X.70(~))
                      A-47

-------
     RETURN
     END
     SUBROUTINE PMETE(HS.RAD.HA,WIND.HBR,P.HETA1R,HC.TDEW.HED>.HE.TA1RJ1CTDEW.HED.HEV,
     +QNET.DMIXZEUPRSECCHI
  C
  MOO FORMAT(/.JX26HMETEOROLOGICAL INFORMATION./.5X26(1H-).
     +/.TX17HNET SOLAR RAD. = .4XF9.2.13H KCAL/M'M/DAY.
     f Z5X13HSOLAR RAD. - , 9XF6.1.8H LANGLEY,
     +/.TX16HNET ATM. RAD. - 4X.F9.113H KCAUM'MT>AY,
     +2SXWHWIND VELOOTY - , 7XF5.17H M.P.H,
     +/.7XUHBACKRAD.  . 9XF9.113H KCAL/M'M/DAY,
     +2SXlHPREaPITATION - . 8XFJ.13H METER(S).DAY.
     +/.7X19HEVAPORATTVE LOSS  ,IXF.9Z13H KCALM'M/DAY,
     +25X12HA1R TEMP. - .12XF5.Z10H DEGREES C
     +/.7X20HCONDENSATION LOSS  . 1XF9.Z1JH KCAUM'MOAY.
.     +25X16HDEW PT. TEMP.  . SXF5Z10H DEGREES C
     +/.7X20HEV. HEAT TRANSFER  ^XF7.<.10H IWDAY OR. .Eltt^.TH M"3.T)
     +2HAYJ.7X19HNET HEAT FLUX IN - . 2XF9.Z9H KCAL/M'M.
     +/.7X1 MKED LAYER DEPTH (M)'.ZXF5.i
     +36X1EUPHOTIC DEPTH (M) '.2XF5.J
     +/8X'SECCHI DEPTH (M) ='.7XF5.2)
     RETURN
     END
     SUBROUTINE SETUP
  C"
C
C
C
C
       1 Determine the initial thicbnos. volume. >nd Area of lvn
       ' and the total vc4ume of above each layer from the depths given
      " in the input dau Hie.
     COMMON/STEPS'DZLL.DZUL.MBOT.NM.NPRJNT.MDAY.MOSTH.ILAYI>Y
     DZ(MBOT)ZMAX-(Z(MBOD +Z(MBOT-1))' J
     Z
-------
   SUBROUTINE SOLVE(VARiMBOT)
C*****
C*"** Tri-dbgonal matrix solving routine
C*****
   COMMON/SOLV/ AK(40).BK(40)1CK(40).DK(40)
   DIMENSION VAR2(40).TX(40)
   DO 60 I-2.MBOT
     TTAK(I)/BK(M)
     BK(1)-BK(I)-CK(I-1)TT
 0    DK(I)=DK(I)-DK(1-1)TT
C"  "BACK SUBSTITUTION""""
   TX(MBOT)=DK DZUL) into iwo or oore
C***** layers of equal thickness. All slate variables are (he ume in
C***** each new layer. Volume is adjusted later.

   COMMON/VOL/ZMAX.DZ(40).Z(*))>\(*)),V(40).TV(40}ATOP(J1).DBL
   COMMON/CHOICE/MODEL.NITRO.IPRNT().NDAYS.NPR.NTC50).NCLASS.PLOr.:-0:
   COMMON/RESULT/ T2(40).CHLATOT(40).PA2(40).rTSUM()l.BODajl<').
   + DSOy40).C2(40).CD2(')),XNO2(40).XNH:().CHLAi-3.JOV
   PC2(3.40),XNC2(3,'10).T20(40),S12(40)
   COMMON/STEPS/DZLL.DZULMBOT.NM.NPRINT.MDAy.MOSTH.!LAY.DY
   DO100K=I.MBOT
   n=MBOT+I-K
   CD2.40).T20(40),SI2(-W)
    COMMONrtEMP/PARJO(4),PCDUM(3,40),XNHD(T(30).VCLASS,PLOT1'30)
    COMMON/CHANEUWCHANL.ELCBjUJ'HA.BW.WLAKE
    COM.MON/WATER^ETA.EMISSJCKLXK2,HKMAXWCOEF>STR
    COMMON/STEPS/DZLL,DZUI_.MBOT.NM.NPRJNT.MDAY.MONTH.ILAY.DY
    COM.MONA'OUZMAXDZ<40).Z(-tO)^(40).V(40).TV(40)^TOP(-
-------
     NITRO=0
     DO 200 1-1.30
   200 NPRNT(I)=0
  C*	INPUT MODEL OPTIONS AND INITIAL CONDITIONS "
     READ(7.') SEKI
      WR]TE('.1006)
      READ(V(A)-) X
      IF(X.EQ.-r .OR. XEQ.y) THEN
      1PRKT(2)=1
      TAPE1=TAPE8
      Tl(LEN8+2) = -O-
      Tl(LEN8+3)-ir
      OPEN'(1FILE=TAPEISTATUS=-NEW)
     ENDIF
     IPRNT(4)=0
   1000 FORMATC PLOT FILE TO BE CREATED (Y/N) ?    'A
   1001 FORMAT(A1)
   1004 FORMATC ENTER UP TO 10 DEPTHS TO BE SAVED V.
     +' END WITH A CHARACTER (.*.*__X>:  V)
   1006 FORMATC OUTFLOW FILE TO BE CREATED (Y/N) 1  V)
     READ(7,') NDAYS
     IF(NDAYSGT.O) READ(7,')(NPRNT(1).1-1,NDAYS)
     READft*) DZLLDZULBETA.EMISS,WCOEF,WSTR
     READ(7,') WCHANL.WLAKEDBLST.S.FT
.    READ(7,') ELCB.ALPHA.BW
    READ(7,') MBOTJ4M.NPR1NT.INFLOW.DY.MONTH1START.MYEAR
     READ(7,-) (Zfl).ll.MBOT)
     READ{7,*) (T2(I),I=1.MBOT)
   1500 FORMAT(/, 5X
     +/JX41HMIN1MUM THICKNESS OF EACH LAYER (DZLL) - J5ZH METERCS)
     +yXBRR'JX-FVBODVJX
     +-(iyDAY)'.2(2X-(GM/M2)-)JX-(MS)ViX5(nj.2X))
      RETURN
      END
      SUBROUTINE STATS(FLDATAJCXIFLAG.DEPTH,I)
  C*****
  C***** Compute lUiiuia >nd tuiaitical quaniiua with
  C***** Y daipialinj field daa tad X daiputinf model raulu.
  C""***
      COMMON/STAT/SUMXY(10)5UMX(10)5UMY(10)_XSO(10).
     +YSQ(10).RSQ(10).RMS(10).RELM(IO).MTHREL(10).MDAYREL(10).ZREL(10).
     +ZRELM(10).RS(10).REL(10).MTHRMS(10),MDAYRMS(10).ZRS(10).ZRMS(10)
      COMMON/STEPS/DZLL.DZUL.MBOT.NM.NPRINT.MDAY.MONTH,ILAY.DY
      SUMXY(I)=SUMXY(I)+FLDATA'XX
      SUMXO)-SUMX(I)*XX
      SUMY(I)=SUMYa)+FLDATA
      XSO(I)=XSO(I) XX'XX
      YSO(I) - YSQ(I) * FLO ATA' FLDATA
      XX=XX-FLDATA
      X2=XX'PLDATA'10a
                        A-60

-------
   X3ABS(XX)
   !F(XlGT_ABS(REm))) THEN
    RL(I)X2
    ZREHO-DEPTH
   ENDIF
   ff(XlGTJ\BS(RELM(l))) THEN
    RELM(I)-X2
    MTHREm)=MONTH
    MDAYREL(t)=MDAY
    ZRELO)-DEPTH
   ENDIF
   ff(X3.GT.ABS(RS(l))) THEN
    RS(I)-XX
    ZRS(r)"DEPTH
   ENDIF
   tF(X3.GTj\BS(R-MS(I))) THEN
    RMSMMONAfOL/2MAXDZ(),Z(-W).A(40).V(40).TV{40)^TOP(41).DBL
   AZ=tt
   AVOL=a
   DO100I=1.MBOT
   II=MBOT+1-I
   AVOL=AVOL+V(II)
   CALL LAKECZDUM^VOL.0.4)
   DZ(II)=ZDUM-AZ
   A2=AZ+DZ(II)
 100 CONTINUE
   RETURN
   END
   SUBROUTINE TVOUMBOT)
C**"**
C***** Determine the volume of water above a layer
C*****
   COMMON/VOUZMAX.DZ(40).Z(-iOyv(40).V(40).TV(40)^TOP(41).DBL
   SUM=tt
   D01001-1.M80T
   SUM=SUM+V(I)
   TVO)=SUM .
 100 CONTINUE
    RETURN
   END
   SUBROUTINE W1NEN(T.V,WIND)
C_
C-CA1.CULATIO.N OF THE SHEAR VELOCITY AND THE WIND SHEAR STRESS
C-CONVERSION OF WIND SPEED FROM M.P.R TO M/S
C-DENSITY OF WATER AND AIR ASSUMED TO BE 1000 AND L177 KC/M3
C_
   CALL LAKE(FTCH.aao,:)
   ZB-ALOG(FTCH)'a8- 1.0718
C-CALCUUAT10N OF WIND SHEAR STRESS
   CZ=SQRT(W)'.0005
   ff(W.GE35.) C2=.0026
C-ASSIGNMENT OF CHARACTERISTIC SURFACE VELOCITY
c-usiNG CALCULATION OF SHEAR VELOCITIES
   V=.03J3>SORT(CZ)'W
   RETURN
   END
   SUBROUTINE WINMLX(ETS,IL.MBOT)
C_
C.CAJXULATE THE AMOUNT OF ENTRAPMENT RESULTING FROM WIND MIXING.
C-USE THE DEPTH OF CENTER OF MASS OF MIXED LAYER TO DETERMINE THE
C-POTENTUL ENERGY THAT MUST BE OVERCOME BY THE KJNECTIC ENERGY
C-OF THE WIND FOR ENTRAPMENT TO OCCUR.
C_
   CMMONA'OUZMAXDZ(),Z<).A(40).V().TV()XrOP(41).DBL
   COMMON/RESULT/ T2(*)).CHLATOT('(0).PA2()).PTSUM('W).BOD2(40),
                         A-51

-------
   +DS02(),CI(40),CD2(40),XN02(40)JCNH2(40).CHLAI(3.).
   + PC2(3,40)XNC2(3.40).T20(40)312<40)
   !F(ILGEMBOT) GO TO 35
   SUM2=0.
   DOI01-1.1L
   ILAY-l
   SUM1=SUM1+RV
 10  SUM2=SUM2+RVZ(I)
   1 = ILAY           ,  .  .
 20  DCM=SUMZ
-------
               LAKE SPECIFIC SUBROUTINE
Area computation section
Depth-area functions of the form AREA=f(ZMAX-Z), written as DUM=f(ZD).
AREA(m:), Z(m).

Fetch computation section
The longest distance across the lake surface area in the wind direction (m).

Volume computation section
Volume-depth functions of the form VOLUME(below depth Z)=f(ZMAX-Z), written as
ZD=f(DUM). VOLUME (m3), Z(m).
                              A-53

-------
            EXAMPLE  LAKE  SPECIFIC  SUBROUTINE
   SUBROUTINE LAKE(ZD.DUM.NFLOW.ID)
   COMMOWMTHD/rAJR(3I).TDEW(3I).RAD(31).CR(3I).WtND(31).
   + PR(31).DRCT(31)
   COMMON/RESULT/ T2(40).CHLATOT(40).PA2(40).PTSUM(40).BOD2(40),
   COMMON/SOURCE/RM(3.40),PROD(40).XMR(3.40).PRODSUM(40)
   COMMONVFL.OW/HMK(4I).pE(4b).FVCHLA(5).PE(5.4I)
   CO.MMON;TEMP/PAR1O(4).PCDUM(3,40)XNHD(<0)_XNODC-10).
   + CHLADU.V{3.40),XNCD(3.40).PADUM(40).SID(40)
   COMMONAX3L^MAXDZ(40).Z<40)^(40).V(40).TV(40).ATOP(41).DBL
   COMMONJSUBJSDZ<60).SZ().LAY(40).AVGI(4.40).SVOL<>)
   COMMON,CHOlCE/MODEl-NrniO,IPR>!T().KDAYS.NPRNT(30).NCLASS.PLT(30)
   COMMONAVATER^ETA.EMISSJ6-ASURF'(ZD+l.)"l.'<209
   RETURN
C* ......... FETCH COMPUTATION SECTION
200 ZD=crt(O.l<59>ASURF)
   RETURN
C* ......... VOLUME COMPUTATION SECTION
300  CONTiNUE
1234 ASUR?=:.?1E
   CONS=0.00>J'ASURF
   RETURN
C" COMPUTE DEPTH FRO.V VOLUME .....
400  CONS=0.0054'ASURF
   DUM1=DUM'CONS
   ZD=(DUM1"0. 406)6). 1.
   RETURN
C ..... WRITE ON THE SCREEN. DAY. MONTH 
 500 WR1TEC.S01) MONTH.MDAY
 501 FORM.-kT(2X' monthM3.' dj> -.i2)
C ..... WRITE EPILIMNION  i HYPOUMNION TEV?ERATL!RES 
   DO IIIH I=1.MBOT
   IF (LEO.:) THEN
   WRITE (2L8TI) MONTH.MDAY.Z(1).T2(I)
   ENDIF
   IF (LEO^2) THEN
   WRITE (n.871) MONTH.MDAY.Z(I).T2(I)
   ENDIF                                    '
 871 FORM.AT(1X.I4.1XI-1.2XF9.2.3X.F9.3)
11111 CONTiNUE
   RETURN
 00 RETURN
 700 RETURN
 800 RETURN
 900 RETURN
 1000 RETURN
 1100 RETURN
 1200 CONTINUE
 1300 CONTINUE
   RETURN
   END
                                               A-54

-------