(   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
hypolirnnetic   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  COa-   The simulations  predict  that  after  climate  change:  1)
epilimnetic  water temperatures will be  higher but  will  increase  less than  air
temperature,  2)  hypolirnnetic  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
                              in

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

         WATER TEMPERATURE CHARACTERISTICS  OF  MINNESOTA
         LAKES SUBJECTED TO CLIMATE CHANGE                 74
         5.1  Introduction
         5.2  Method of lake temperature modeling
         5.3  Climate conditions simulated
         5.4  Regional lake characteristics
         5.5  Simulated lake water temperature regimes for
             historical and future weather
             5.5.1  Water temperatures
             5.5.2  Thermal energy flexes
             5.5.3  Vertical mixing/stratification/stability
         5.6  Conclusions

     6.   SUMMARY

     7.   REFERENCES
APPENDIX A.
APPENDIX B.

APPENDIX C.

APPENDIX D.

APPENDIX E.
Vertical diffusion in small  stratified lake:
Data and error analysis
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.7   Error analysis
A.8   Conclusions

Temperature equation  discretization

Temperature equation  linearization

Cross-term evaluations
                                                   74
                                                   76
                                                   76
                                                   80

                                                   83
                                                   83
                                                   89
                                                   95
                                                  102

                                                  103

                                                  105
                                                                     A.I
A-17



A-24

A-25

A.28

A-30
Regional lake water temperature  simulation model
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

-------
                               last 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 axe shown.

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 2°C.
              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.
                     m
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 (2av )  for  different  sampling
                                                  JVz
              intervals.

Figure A.12   Estimated   eddy   diffusion  errors   2av (cmV1)   space-time
                                                     -K-z
              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  m"2).

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  2xCOs  climate
              model  output  for Minneapoiis/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  Kz =  o(N2)T



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  Laboratory—Duluth, 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.
                                   XI

-------
                      s
              1.  Introduction and  Literature Review
1 1    Introduction

      The  concentrations  of some gases  (CO2,  H20,  N20,  CH4)  have  been
increasing  in  the  atmosphere  (Bolin  and  Doos,  1986;  NRC,  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  C02  alone.   Mathematical  models  of global
climate  change lead  to  the  conclusion  that   the increase  in  mean global
equilibrium surface temperature for  a doubling  of CO 2 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
(McCormick,  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 nqrthcentral 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., 19SOa; 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:
where  T(z,t)  is  water temperature  as a function of depth  (z)  and time  (t),
A(z) is the horizontal area of the lake as  a  function of depth, H(z,t) is  the
internal  distribution  of  heat  sources due  to radiation  absorption inside  the
water  column, pw is  the  water density, and  cp is the specific heat  of water.

     The vertical temperature profile in  the lake is computed from a balance
between  incoming heat from solar and longwave radiation and the outflow of
heat through  convection, evaporation,  and  back radiation.    The  net  increase
in  heat  results  in  an  increase  in  water  temperature.   The  heat  balance
equation (see  also Fig. 1.1) is given  by

                     H=H   +H+H+H+H,                 (1-2)
                      n      sn      a     c      e     br

                                                         _9   _ 1
where  Hn is net heat input at  the water surface  (kcal m "day  ), Hsn is  net
solar  (short  wave)  radiation, Ha is  atmospheric long wave radiation,  Hc is
conductive loss (sensible  heat),  He is evaporative loss  (latent lieat),  and  Hbr
is   back  radiation.    The  heat  budget  components  in  equation  (1.2)  are
computed as follows:

                                            - /J)H                      (1-3)
                                                  §

-------
    Ha

atmospheric
  radiation
        Hsn
He
He
Hbr
incoming   reflected evaporation   convection   back
  solar     fraction                          radiation
radiation
                                         t          t          t
t
'

fe>. surface absorpt

	 ^ absorption with
fraction transmited
to next layer
- j
j 	 ^ absorption with
i
i
j fraction transmited
on
depth


depth

                                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:
where Hsn(i-l) is  solar  radiation at  the  top of a  horizontal layer of  water
        _ n   _ -I
(kcal m  day  ), Hsn(i) is  solar radiation  at the  bottom  of  a layer,  A  is

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

                         H = p  +  p -SS + /x.Chla                    (1.5)
                         r    'w    'ss       'ch                        ^   '


where /^ is the extinction coefficient of lake  water (m~~ ),  pss is  the  specific

extinction  coefficient due to suspended sediment  (1 m^mg"1); SS is suspended

inorganic sediment concentration (mg  I"1); /ich is  the  extinction coefficient due

to   chlorophyll   (m2   g"1Chla)(Bannister,   1974),   Chla   is   chlorophyll-a
concentration (g  nr3).

                                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 T,
latent heat flux  H  . and  the sensible heat  flux  H   across  the  water  surface
                  e                              c
(Keijman,  1974;  Ford and  Stefan, 1980;  Strub and  Powell,  1987; Sadhyram et
al., 1988):
                      T = p   U'W'  = p Ml = p C,  U2                 (1.7)
                                        r  *                            v   '
                                                   ,
                            a            ra *     a  d   a
           H  = p c   6 ' u'  = p c C u, 9. = p c f(U )(T  - T )       (1.8)
             c    an            ra D s *  *   ra D v  a'x q     a'       v   '
                                                             - qa)       (i-9)

where  T is the surface  wind  stress,  pa  is the  density of the air,  u'  and uj'
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,  C
-------
*iih u. axe expressed as a function of wind speed, f(Ua),  (Ford,  1976), Ts  is

cater surface  temperature,  Ta air temperature above  the water surface, Lv  is
.lu-nt  heat  of vaporization, q' is tiie specific  humidity  fluctuation, q^ is the

 peafic  humidity scale, qa  is  the specific humidity above the water surface,  qs
-s   the  specific   humidity  at  saturation  pressure  at   the  water   surface
•-emperature.

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

day1), 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
Radiation extinction
by water
Radiation extinction
by chlorophyll
Symbol Units
V* (m-1)
/zch (m* g'1 Chla)
Range of values
Previous Literature
Simulations Values
0.4-0.65 0.02-2.0
8.65-16.0 0.2-31.5
Wind sheltering          Wstr    (-)

Wind function
coefficient                c        (-)

Maximum hypolimnetic
eddy diffusivity          Kzmax   (m2  day1)
                                                  0.1-0.9       0.1-1.0


                                                  20.0-30.0     20.0-30.0
                                                  0.1-2.0
                                                       0.086-S.64

-------
      Radiation  extinction  coefficients  by  water  (/%)   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 fa as  a  function of the wavelength.   Values  of /^ in  the
range  of 0.68 ± 0.35  (m'1)  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  /ich  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) =  C  Ua                          (1.11)

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 nr2dayl.

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


                                Kz = a (N2)7                           (2.1)

where stability  frequency N2=-(5p/5z)(g/p)J  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  sheax  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-* (Area)°'56 (N2)-o-43               (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   Kzaax   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   hypoHmnetic  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.

-------
                Ann Lake  (0,47 km2)
 o
 «U
 to
 O
     10-3-
 o
 0>
 
n
 E
 o
M
     ID'1-:
     10
       -2_
     10-
         o-
Cornelian Lake  i 1.89  km')

 Hr4          10"3         10
                                   Square Lake (0.85  km'

                                   KZ=1.0X1Q-J(N2)-°4J
                                                   ®   Big Marine Lake (1A  km'
                         N'(soc-2)                               N3(sec-2)


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

-------
            = 8.17X10~'t(Area)0
           Ryan Lake (0.06 km2)
       4-   Ann Lake (0.47 km
       O   Square Lake (0.85 km2)
       W   Big Cornelian Lake ( 1.89 km")
       A   Big Marine  Lake  (7.4 km2)
       A   Lake Mcllwanie (26.3 km2)
                             Area  (km2)
Fig.   2.2   Hypolimnetic eddy diffusivity forcing parameter (a) dependence on
           lake surface area.

-------
 U
 0)
 CO
CM
 U
  N
       102u
       10'-
10°-
        U '
D

A
A
X
X
T
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)
       :0.048(Area)°-56
                                         —rTTrrn	1—i
                                    Area
            Fig.   2.3   Maximum hypolinmetic eddy diffusivity (N2=7.5*10'5 sec"2)
                       dependence on lake surface area.

-------
  '.'  'l      Attenuation coefficient

      The specific radiation attenuation  coefficients for water  and chlorophyll
 *••:••  replaced  by the total attenuation  coefficient.   This  was  done following
  ,'  parsimonious principle i.e. the fewer  coefficients in the model, the less
 .::  main the model estimate.  In addition uncertainty analysis  showed that
 :.,  .rophyll-a  made  a   minor   contribution  to  lake   water  temperature
 ,n  rriaanty.   A  relationship  between total attenuation coefficient [i  (nr1) and
 ^
-------
                          [—i—r-~i—r~]—r
          fj,= 1.84*(secchi depth)
       T—i—|  I—r—i—i—|—i—i—i—i—|—i—i—i—i—|—i—i—i—i—j—i—i—i—i—|—
         0.3     0.5    0.8    1.0     1.2     1.5
                     (Secchi  Depth m
Fig.   2.4    Relationship between total attenuation  coefGcient and  Sechi disk
            depth.

-------
 C
 0)
'u
 0)
 o
O
 0)
•4— •
13
T)
C
 1 O1 -
10-'
      o-*-
                 Wstr= 1.0-EXP(-0.3+Area)
                    mi,
                     10-2
                                        X
                                       10°
                            Surface  Area  (km2)
10'
I02
                 X   Williams (0.35)
                 X   Riley (1.2)
                 *   Waconia (10.0)
                 O   Square (0.8)
                 A   Thrush (0.07)
                 A   Greenwood (7.70)
                 m   Fish (1.16)
                 D   Elmo (1.23)
                 •   Cedar (3.30)
                 O   Calhoun (1.71)
    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 frojn 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-
 y
'<-»—
"
 o
o
    20-
    15-
      10
— 2
      Wcoef = 24  + In(Area)
o-1            10°
    Surface Area  (km2)
10'
                                                              X   Williams  (0,35)
                                                              X   Riley (1.20)
                                                              4   Waconia  (10.0)
                                                              O   Square (0.8)
                                                              A   Thrush (0.07)
                                                              A   Greenwood (7.70)
                                                              BB   Fish (1.16)
                                                              D   Elmo (1.23)
                                                              •   Cedar (3.3)
                                                              O   Calhoun  (1.71)
                                                          O2
      Fig.   2.6    Wind function coefficient dependence on lake surface area.

-------
       «     off—lake  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
 liffusivity  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
(small-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  (4°C) 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
                 ViTSi           2^ViTmi

     Ts  = _iZi	   Tm = -i^	                          (2.6)
              X  v-
                                     19

-------
          Northern  Minnesota Wetlands
                                       GREENWOOD  LAKE
                                      THRUSH  LAKE
         WILLIAMS LAKE
               Northern Lanes and  Forests
                                           Twin Cities
                                          Metropolitan Area
           North Cent-'aFHardwood Forests
    -H,'Western Corn Belt Plai;
               FISH LAKE
                     Driftless Area
  0.
100.
Fig.   2.8    Ecoregions and spatial  distribution of selected lakes.
                                 20

-------
           10-
       1Q-1       10°       10'
             AREA (km2)
         100
                       MAXIMUM DEPTH (m)
           0~
                        SECCHI DEPTH (m)
                                                     A
                                                     n
                                                     o
                                                     A
                                                     A
                                                     A
                                                     D
                                                     O
                                                     A
                                                     A
                                                     A
                                                     •
                                                     m
                                                     D
                                                     •
                                                     ©
                                                     O
                                            Waconia  (10.0)
                                            Greenwood (7.2)
                                            Cedar (3.3)
                                            Williams (0.35)
                                            Caihoun  (1.7)
                                            Elmo (1.2)
                                            Fish  (1.1)
                                            Square (0.8)
                                            Thrush (0.07)
                                            Elmo (42.0)
                                            Greenwood (34.0)
                                            Caihoun (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)
                                            Caihoun (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

            ;i-Tmi)2-|0.5                 r IrfVi

                                   £2=
                                                    (Tsi-
                                                   X"
                                                                       (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.   Vi 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
Ei

(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
r2 Number of

(-)
0.97
0.93
0.92
0.90
0.97
0.92
0.93
0.96
0.97
field

136
20
214
32
136
43
46
114
110
data










Average


Tm -

T  —


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

                          5/20/71
                                                            6/08/71
                          6/29/71
                                                            7/15/71
                         8/02/71
                         10/11/71
                                                            10/25/71
      10     15     20     25
           TEMPERATURE  (*C)
10      15      20
    TEMPERATURE CC
Fig.    2.10   Lake Calhoun water temperature profiles.
                         23

-------
X  ,0-j

a.
UJ
a



   15-
                            measured

                            simulated


                                V15/85
                                                                   5/04/85
                                6/11/85
                                                                   7/12/85
X  10-

s:
UJ
a


   15-
S/07/S5
                                   9/05/S5
                10     15    20    2S    0


               TEMPERATURE (-C)
                   10    15    20


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

-------
N3
Cn
                              X
                               ^
                              UJ
                              O
                              a.
                              UJ
                              Q
u~
1-1
2-
3-
4-
6-
7-
8-
9
10-
1-
2-
3-
4-
5-
6-
7-

H
9-
10-
1 1 -
• measured
	 simulated

5/17/85
1

•
•
T







8/20/85




.
•
•


















-i








, ;
•
«
•
•
•
•

•
•
•




7/22/85

















                                                10     15      20

                                                TEMPERATURE (°C)
                                                                                                           9/18/85
25
                    10      15     20

                    TFMPFRATURE (°C)
                                        25
                                              30
                                                 Fig.   2.12   Waconia Lake water  temperature profiles.

-------
   «
z
a.  e-
UJ
o
   10-

   12-
   14'
                        •   measured
                       	 simulated
                                5/U/86
                                                                  6/17/66
   2-
Q
  10-

  12-
  14-
                               7/15/86
   8-

  10-

  12-

  14
9/16/86
          S     10    15     20     25
               TEMPERATURE (-C)
             5     10    is
                  TEMPERATURE
      Fig.    2.13    Thrush Lake  water temperature profiles.
                                 26

-------
£  *.
                             5/17/8*    j
               S/JO/S*
  10-f-
                             6/1Z/B4
£  4-
                             8/08/S4
                                                               S/26/8<
                             9/!3/8<
                                                               11/02/84   1
               10     15    20
              TEMPERATURE 1'C)
10    15     20    25
TEMPERATURE {'C>
      Fig.    2.14    Williams Lake water  temperature profiles.
                             27

-------
      Volume  weighted and unweighted root  mean  square error was less than
1°C 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 (Ei and  £3)  indicate the vertical position
of  the  maximum simulation  error.    If  Ej  is greater  than EI,  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

Umax
(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
(-)
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 /^ch
(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^84
l^l7
11-348
1-35
2-4 6
3-79
Field data given by:
                  i
                  2
                  3'7»8
                  4'5
                  6
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  1°C,  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 Hypolirnnetic Closure

     Uncertainty in the lake water, temperature simulations  was considered  in
'.-•rras of all  model coefficients  except  maximum hypolimnetic  eddy diffusivity
i,-  specified  in Table 1.1.   To first-order  the  uncertainty  in  lake water
:< niperature  depends  on the uncertainty in  the model  coefficients,  and on the
-•nsitivity of the lake water temperatures to changes in the coefficients:

     PT  =  E{[T - T][T - T]tr] «

          E{[T(u) +  — (u - u) - T][T(u) + — (u  - u) - T]tr} =
                      du                       dn

              — (u - *)][—  (*  - *)]*}=
              on          on
             - E{(u - u)][(u - u)p }
           dn                         dn
             • Pu                                                       (2.8)
           dn     dn

* here  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
   the transpose,  u is the  vector  of the n  coefficients, Pu  is  the (n  x n
                                         orp
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
Mw (m-i)
p,± (m2 g-iChla)
W'str
c
Lake Calhoun
mean st. dev.
0.65 0.20
8.65 2.65
0.60 0.18
24.0 7.20
Williams Lake
mean st. dev.
0.65 0.20
8.65 2.65
0.20 0.06
20.0 6.00
Cedar Lake
mean st. dev.
0.65 0.20
8.65 2.65
0.60 0.18
24.0 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 1° C,
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 to 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
oHgotrophic  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  (4° C) 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
               wind sheltering
               wind function
               water attenuation
               chlorophyll —a  attenuation
EPILJMNION
LJ
               wind  sneltenng
               wind  function
               water attenuation
HYPOLJMNION
          	chlorophyll — a  attenuation
    0.5-
    0.0
          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.0J
              wind sheltering
              wind function
              water attenuation
              chlorophyll —a attenuation
   0.0
              wind sheltering
              wind function
              water attenuation
              chlorophyll —a attenuation
                                                    EPIUMNION
                                                    HYPOLI.MNION
   0.0
         MAR   APR   MAY   JUN   JUL   AUG   SE?   OCT   NOV
Fig.    2.16   Standard deviations of estimated  lake  water temperature
             uncertainties.
                                 32

-------
                         CEDAR LAKE
o
o
     2.5-
     2.0-
     1.5-
     1.0-
     0.5-
     0.0-H
2.0-
     0.0
           wind sheltering
           wind function
           water attenuation
           chlorophyll—a attenuation
wind  sheltering
wind  function
water attenuation
chlorophyll—a  attenuation
                                                    EPiLIMNION
                                                   HYPOLIMN1ON
          MAR   APR   MAY   J'JN  . JUL    AUG   SEP   OCT   NOV
     Fig.   2.17   Standard deviations of estimated lake water  temperature
                  uncertainties.
                                 33

-------
                           SIMULATION YEAR -  1986
              MAR   APR   MAY  JUN   JUL   AUG   SEP   OC
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

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


Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Average
Year


1971
1984
1988
1987
1985
1985
1986
1986
1984

Regional model
Tra
CC)
14.37
20.64
13.94
24.40
14.37
20.14
11.80
11.91
17.26
16.54
Ts
CC)
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
E2
(°C)
0.89
1.15
1.93
0.89
1.03
0.68
0.99
0.97
1.25
1.10
r2
(-)
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
AE,
CC)
0.16
0.13
0.06
0.07
0.38
-€.10
0.35
0.05
0.18
0.14
Differences
- validated model
AE2
CC)
0.10
0.16
0.13
0.07
0.24
-0.05
0.20
0.06
0.18
0.12
Ar2
(%)
-1
_2
-2
-1
-2
2
-1
-1
-1
-i

-------
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   epiKmnion
temperatures  showed  1° C  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
 '-".or state-space  method.  The output uncertainty is  defined as the result of
 irviations of  the  meteorological  variables  from  their  mean  values.    The
 .,-,aiysis  is applied to systems with  correlated  and uncorrelated  meteorological
 ir.abies.    Surface  water  temperatures  are  strongly  affected  by  uncertain
 •*,' -teorological   forcing.     Air  temperatures  and  dew  point   temperature
  filiations  have significant effect on  lake temperature uncertainty.   Ignoring
  rrelation  in   meteorological  variables  underestimates  uncertainties in  lake
  rr.perature  estimates.    Long-term  average  water temperature structure  in
 .K"S can be estimated by computer model  simulation  for just one year when
 -suits  from a  statistical analysis  of meteorological  variables   are used  as
 -.put.  The analysis  presents a useful  alternative  for the study  of  long-term
 .v.-rages  and  variability  of water  temperature  structures  in  lakes  due  to
 inable  meteorological forcing.


3 1   Introduction

     It  was shown  in  Chapters  1  and 2  that  vertical  water  temperature
; roGles  in   lakes are related to  meteorological variables by heat  transport
•--1 nations 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
ind 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
rnanagement  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 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+l) 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
«»iues  for the meteorological  variables.   The required meteorological variables
-'••   solar  radiation, air  temperature,  dew point temperature, wind speed and
 ..'•-ction.  Initial conditions,  model setup parameters,  have to be  provided  to
 •«  the model.

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


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

*< rnmal  (mean) values  and  first order  derivatives  evaluated at  these values
.:<•  denoted by  circumflex.   Perturbations  of  the  water temperatures T'(k),
.• : 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

:>'terministic  temperature model.   Matrices  B(k) and  F(k) require evaluation
 :"  the  first  order  derivatives of  all  terms  in equation  (3.1) which  contains
:.ike  temperature  and meteorological variables  at   time  step k  respectively.

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  
-------
Ws(k)
 Fig.  3.1    Schematic illustration of the lake temperature perturbation
            system.
                                 40

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

-------
     A recursive, solution of equation (3.4) is
     T'(k) = #k,0)T'(0) +   "  #!,n+l)  #n)  C'(n);    T'(0)=0
                             n=0
     T'(k) = S   <£(k,n+l) $n) C'(n)     (3.6)
              n=0

 rutial  conditions  are  assumed  to  be  known  with  certainty  i.e. T'(cn=0,
:tk.s)=^(k-l)^(k-2)...^(s), and  $(k,k) is an  identity  matrix.   Equation  (3.6)
:-.iys  that  temperature  perturbation at time  step k is a  linear combination of
;he  meteorological forcing perturbations C'(k), C'(k-l),...,C'(l).

     The  first  order  estimate  of  the  mean  and  covariance  temperature
;*-rturbations are obtainable  from equation (3.6) (Protopapas, 1988).   In the
::fference equation form
     T'(k+l) =  #k)  T'(k)                                            (3.7)


     £T/T/(k+l) = E [(T'(k+l)  - f '

          E   T'k  - T'kT'k  -  f
          E  [(T'(k) -T'(k))  C'(k)T]
          E  C'kT'k  -   '
          E  [C'(k) C'(k)T]  ^
-------
    ET/T/(k) = IT ^(k>n+l)^n)Muc(k)k)^n)T^(k,n+l)T;  ET/T/(0) =  0,


or in  difference .equation form


          ET/T/(k+l) = #k) ET/T/(k) #k)T+ ^k)Muc(k,k)^k)T     (3,10)


where  the  covariance  matrix  Muc(k,k)  has  diagonal  terms equal  to  the
variances of the perturbed meteorological variables  (Appendix C).

     If weather  perturbations  are  correlated  between  successive  days, cross
terms (third  and  fourth) in equation (3.9) have to  be evaluated.  If we define
a disturbance covariance matrix as


                M(n,k) = E  [ C'(n) C'(k)T]  = S(n)McS(k)T
then
                                  k 1
              E [T'(k)  C'(k)T] = E~ <£(k,n+lMn)S(n)McS(k)T        (3.11)
                                  n=0

where Mc is  a correlation matrix between successive days,  and S(n), Sfk)  are
"standard deviations" of the covariance matrix M(n,k).  If  we define N(k)  as:


     N(k) =  S"1  $k,n+l)$n)S(n)

          =  J(k-l)N(k-l)  -f ^(k-l)S(k-l)                            (3.12)

additional cross terms can be  written in difference  equation form:


     Pt(k) = <£(k)N(k)McS(k)TV
-------
     Lake  temperature covariance propagation is  calculated  in  the following
s.  ;i   (l)  Set  isothermal  (4°C ) initial  steady state conditions for  lake  water
 ••. ;«-ratures   in   spring,   initialize   covariance  matrix   of  meteorological
 «•".•.;: bat ions; (2) read meteorological Variables, mean and  perturbation values,
    the  next time  step;  (3)  using mean  values for  meteorological  variables,
   -,;me  first  moment  of  lake  temperature  profile for  the  next  time  step
                                                          *     A
-  , ;.i',ion (3.5), store matrix Ac(k); (4) compute matrices B(k), F(k), i.e. first
 •:- r  derivatives with respect to  the  perturbed meteorological  variables and
•-.-..mated  lake  temperatures (Appendix C); (5)  calculate matrix N(k)  for  the
  : Mated  case  (Equation  3.12);   (6)  compute transition matrix  $(k), and
-.«»,:   (7)  calculate  additional term  P(k)  (Equation  3.13) for the  correlated
 i.v  (8) propagate and store temperature covariance  matrix £  ,   ,  for  the

 -i:  time  step  (correlated case  Equation  3.14, uncorrelated  Equation  3.10);
   -.lore transition  matrix $(k),  and ^(k), if  correlated case,  store in addition
 •  i   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
 if is eutrophic with maximum depth of  24  m, and  surface  area of 1.7 km2.
'.' • J,eorological  data  used  are from  the  Minneapolis   St. Paul International
i. :-,>.)rt,  located 10  km from the  lake.

     Daily  meteorological   data  time  series   (1955-1979),  averaged over  25
>'irs, 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
 standard  deviations) for meteorological forcing variables are  also obtained  by
hrcct data processing.  They  describe  weather variability  over 25  years  for a
; articular  day.    Standard  deviations were  higher  in  spring   and  fall than  in
summer (Fig. 3.2).

     The  time  series for each meteorological  variable is reduced to  a residual
•••Ties  by removing  periodic means and  standard deviations as pointed  out  by
k;:hardson (1981).    The dependence among  the meteorological  variables was
i> scribed  by calculating   cross  correlation coefficients  of  the residual  time
•Ties.  The serial correlation coefficients for time lags up  to  3 days are  given
;r.  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
 the  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

-------

c%r
i
E
CJ

"o
o
Z
O
5
Q
•C
rii
cc
5
O

700-
600-


500-


•100

300-


200

100-
ri
            overage

            standard deviation
         Af'R    MAT    JUN   JUL    AUO
                                           SEP
    35
UJ
a:
a:
LU
a
O
CL
u
Q
30:


25


20


15-


10


 5-.


 0-
   -10-
            standard deviation
                                                 OCI
         APR    MAY    JUN   JUL    AUG
                                          -T—,-T-

                                           SEP
                                                               averagff

                                                               standard deviation
                                                                                                      35




                                                                                                      30



                                                                                                     -25




                                                                                                     ;2Q



                                                                                                     ;I5



                                                                                                     -10
                                                            .  , .  . , ,  ... ,  ,i i  | .  .  , ,  .  , ,  ,  , ,  .  | ,

                                                            APR    MAY    JUN   JUL    AUG   SLP   OCI
                                                                                                              CJ
                                                                                                              bJ
                                                                                                              0.
                                                                                                              2
                                                                                                              u
                                                                                                              (—

                                                                                                              O-

                                                                                                              <
                                                               startdord deviation
                                                 OCT
                                                           A                                             I
                                                           '••'.M.-.r-.rfS*.;'' .                           ,(,  •,   '•-.
                                                            <  v v •' • • r '••'-•••^'••'"\ff,:;._..;-. rv./-vi^-/.-/-^'v'';'<;;'" '•••':'
                                                            APR    MAY    JUN    JUL    AUG    SEP   OCI
                                                                                                          •10
                                                                                                          -6
                                                                                                          -4
                                                                                                              Q
                                                                                                              L.J
                                                                                                              UJ
                                                                                                              Q.
                                                                                                              Q



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

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

Solar
Radiation
Air
Temperature
Dew Point
Temperature
Wind
Speed
Solar
Radiation
Time Lag (days)
0 1 2
1.00 0.39 0.14
0.18 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
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 (<7eav)
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 1*C 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

-------
APR
MAY
JUL
AUG
SEP
                                                        OCT
Fig. 3.3   Estimated  long-term average epilimmon and hypolimnion
          temperatures.
                            47

-------
 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
u nperatures are given in  Figs. 3.5 and 3.6, respectively.   Contributions by
 ••rturbations of  individual meteorological  variables  perturbations  as well as
 :.!•  total  contribution  of all perturbation  variables  were  calculated  with
  rrelated and uncorrelated input  variables at a daily timestep.  Air and  dew
,. ;nt  temperature contributed the  most  to the  temperature uncertainty, while
- iir radiation  and wind  speed had  smaller effects.  Furthermore, the overall
:n certainty  in  water temperature  was found  to be larger in the case of the
 undated daily process than  in  the uncorrelated one.   Uncertainty in  lake
• aier temperature  varies  with time, since sources of  uncertainty  vary with
'.une.     These  sources   are,   the  sensitivity  of   lake   temperatures  to
"leteorological variables as  well as the  amount  of the  uncertainty  concerning
these  variables. At  the beginning  of the simulated period uncertainty was  set
•-•) zero  since  initial  conditions were considered  perfectly known.   Isothermal
r.nial conditions  of  4°C (after ice thaw) April   1 are  appropriate for the 45°
.uitude.   Although  isothermal  water at 4°C may not  exactly exist on April
   thermal  inertia of the water makes summer predictions insensitive to  initial
  .nditions  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
incertainty  in  spring,  more or  less   constant  uncertainty  after  onset of
 t ratification in  summer,  and  decreasing  uncertainty  in   fall  when   lake
.- mperature is  driven towards isothermal conditions.  Temperature uncertainty
s  decreasing in fall  when observed meteorological  variables and estimated  lake
•vater temperatures  are both  decreasing.  First  order derivatives  with respect
:  > the lake temperatures and meteorological  variables  are evaluated  at  these
if;served  and estimated values  respectively.   Thus, they have less  weight in
uncertainty  propagation.

     Uncertainty propagation for  deep  hypolimnetic  temperature (1m  above
iake  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,
,ind 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

Ld
^   4.0-1
Ld
Q_
2
Ld
o;
Ld
t—
        -1
      .0
     2.0-
      .0-
               all  inputs
               air temperature
               dew point temperature
               solar radiation
                                  UNCORRELATED  DAILY PROCESS  -
          	wind speed
O
o
LU
cr
P
<
ry
uJ
Q_
2
LJ
fy
              all inputs
              air temperature
              dew point temperature
              solar radiation
                                  CORRELATED DAILY PROCESS
          	  wind speed
     0.0-t
           APR     MAY
                           JUN
Ui
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 S!MULATED DEEP HYPOLIMNiON TEMPERATURE
               all inputs
                                  UNCORRELATED DAILY PROCESS
o

LJ
             -  air temperature
             •-  dew point temperature
             -  solar radiation
          	  wind speed

-------
      I 12-
      51






/ ~Z*a
; To>-<7
/ «/10






•" jf

r
1 =T:*a :
,„- -
i 6/08
"
•
;
     J.

     X 12-
                               6/29
        16-


        20-
                 10   IS    20   2

                  TEMPERATURE fC)
Fig.  3.7   Long-term average temperature profiles plus or minus one
           standard deviation in Lake Calhoun.
                                52

-------
Cn
CO
             O
             o
             LU
             o_
             cr
40



35-



30-;



25-



20-



15-;



10-



 5-



 0-
                           T
                            eav
eav
          EPILiMNION
                                      av
                       APR    MAY     JUN     JUL    AUG     SEP     OCT
                     Fig. 3.8   Epilirnnion temperature long-term average plus or minus one

                              standard deviation.

-------
data.     Long  term  average  (Teav)  and   standard  deviations  (<7eav)  were
estimated  from  the  simulated  lake water  temperatures  over  tie  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  (o"av)  was  estimated  using the proposed
perturbation   model.     Results   shown  are   for   correlated  meteorological
perturbation  variables.   The maximum  difference was  less than  2"C  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 4" C 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  1°C  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
instrurnentally  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 ia  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
[10«m3]
17.1
16.5
0.65
Secchi
Depth
[m]
2.5
2.8
2.2
Typical
Summer, Chla
[g m-3]
20
8
28
      Meteorological  data used are from the Minneapolis-St. Paul International
Airport located 5 to  18 miles from the  studied lakes.   The meteorologies.
data  file contains  measured daily  values  of average air temperature (Ta), de**
point temperature (T^), precipitation  (P), wind  speed  (Ua) and solar  radiation
(Hs).   Mean and  standard deviations  (S.D.)  for  those  parameters  average-:
over  the simulation  period, from  May through  September, are given  in Tab.*
4.2.   Mean  summer  air  temperature  in 1988 (^1.6 °C) was 2.9°C  higher tha:
in 1971 (18.7 °C).  May to  September is  the  main period  of  interest.   Mea_*
April air  temperature  was about the same in  1971  and 1988, but  Octo'tx-
1988  was  much colder than  normal.   Wind,  the most important   extern.
hydrodynamic  force  causing mixing  in the lake, had  similar  values  for bcu
periods in  terms of mean and  standard  deviation.   Mean solar  radiation  w^
18%  higher  in  1988  than in 1971.

      The model  assumes  isothermal  initial  conditions  of 4"C on  April
This  is appropriate  for the 45'  latitude.   Dates of ice  formation,  thaw  i:
duration have  been continuously  recorded  on  Lake Mendota  (Wisconsin,  C
latitude) since  1855.  The  mean  date of ice  thaw was April  5  with  11 di*
standard deviation  (Robertson, 1989).   Model  sensitivity to  the  date  of "•>
initial  isothermal   conditions   is  summarized  in  Table 4.3.    Epilimr."*-
temperatures are  very  well simulated throughout  the entire  summer  ptr •
                                     56

-------
               Table 4.2   Mean Monthly Meteorological Data


iril
M \Y
,.N
: ' ' L
U'G
• r'.P
;:T
Max

14.9
19.4
27.4
26.5
27.5
22.9
15.7
'!EAN(MAY

24.7
Ta
[°C]
Min

1.7
6.5
16.5
14.3
14.2
11.3
5.8
Aver. Diff. from
Normal*

8.3
13.0
21.9
20.4
20.9
17.1
10.8

0
-1
2
-2
-0
1
0
Year
.6
.3
.0
.3
.8
.2
.5
Td
[•C]
1971
-1
2
15
12
12

.7
.8
.0
.8
.2
10.0
6
.7
P
[mm]

0.9
2.6
3.0
3.2
1.4
2.3
4.6
W
[ms-1]

5.
4.
4.
4.
3.
4.
4.

1
6
0
1
8
1
5
HS
[cal cm^d"1]

411
482
450
563
479
338
192
to SEPT)
12.6
18.7
-0
.2
10.6
2.5
4.
1
462
- (MAY to SEPT)


APR
MAY
JUN
JUL
AUG
SEP
OCT
3.15

15.1
25.7
30.5
32.3
29.3
22.8
12.5
3.45

2.0
11.4
16.6
18.8
17.3
10.9
0.8
3.26

8.5
18.5
23.5
25.6
23.3
16.9
6.7
1

0
4
3
2
1
1
-3
.60
Year
.8
.2
.6
.9
.6
.0
.6
4
1988
-2
7
12
14
15
9
-0
.20

.4
.1
.3
.5
.3
.8
.6
0.63

1.3
1.4
0.2
0.9
3.5
2.4
0.7
0.

4.
26

6
5.1
4.
4.
4.
4.
4.
8
4
8
8
7
72.7

469
584
654
610
497
331
284
MEAN(MAY - SEPT)
      28.2   15.0  21.6     2.7

c (MAY -  SEPT)
       3.42   3.23  3.29     1.20
11.8
1.7    4.8
 3.03    1.16  0.23
535
                  114
*Normal  is the 30-year average from 1941 to 1970
a =  standard deviation
                                   57

-------
Cn
00
Table 4.3 Differences (°C1 in simulated mean daily epilirnnetic and hypolimnetic temperatures for different
starting dates of the model (April 1 reference)
MAY
JUN
JUL
Epilimniou
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
-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  4" C  after
start  of  the simulation.   Although isothermal water at 4°C may not  exactly
exist  on  April  1,  thermal inertia  of the  water  make   summer  predictions
msensitive  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 (4° C) 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.4°C.


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 (Ht,r) 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
surface—water  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 -
c.
UJ
Q
    40 -
    Fig. 4.1     Lake Elmo water  temperature profiles.
                         60

-------
                   TEMPERATURE  (°C)
CL
UJ
c
  Fig. 4.2     Lake Holland water temperature profile
                          61

-------
                    TEMPERATURE (°C)
5


X
I—
Q.
UJ
D
Fig.  4.3     Lake Calhoun  water temperature  profiles in 1971.
                         62

-------
     Table  4.4  Monthly averages  of daily heat balance  components [1000 kcal
Lake
Hsn Ha
Calhoun
Hbr
He
Hc
Hn
Lake
Hbr H
Elmo
H0
Hn
Hbr
Lake
I
Holland
Ie 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 -8.71
6.87 -8.36
6.18 -7.69
SEPT)
6.98 -S.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 -4).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
2.17
1.85

1.77
-0.27
0.13
-0.81
-1.75

0.53
-7.03
-7.76

-8.63
-8.83
-8.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
OS
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 -1.29 0.16 2.01 -6.96 -1.15 0.39 2.49 -7.25 -1.71 -0.12 1.13
6.84 -8.06 -2.49 0.32 2.18 -7.84 -2.00 0.71 3.28 -8.21 -3.07 0.10 1.23
7.52 -8.99 -4.40 -0.18 0.22 -8.87 -4.29 -0.04 0.5!) -9.01 -4.63 -0.20 -0.06
7.84 .-9.17 -4.06 0.03 0.47 -9.11 -4.14 0.11 0.52 -9.14 -4.15 0.05 0.43
7.56 -9.00 -3.74 -0.21 -0.67 -8.99 -3.98 -0.21 -0.89 -8.95 -3.70 -0.15 -0.52
6.79 -8.23 -2.32 -0.27 -0.92 -8.24 -2.55 -0.31 -1.20 -8.16 -2.21 -0.18 -0.65
5.54 -7.50 -1.97 -0.79 -2.03 -7.39 -1.86 -0.47 -1.44 -7.30 -1.65 -0.48 -1.21
SEPT)
7.31 -8.69 -3.40 -0.06 0.26 -8.61 -3.39 0.05 0.46 -8.69 -3.55 -0.08 0.09

-------
      Atmospheric  long wave radiation  and  back  radiation  from the  water
surface  axe  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° G  air
temperature increase are plotted in Fig.  4.4.  Cumulative evaporative  loss was

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  0°C  crossings and  hence the
date  of  ice  formation might  shift  from  year  to   year.    Under  wanner
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 wanner  year was   noticed  probably  because  the fall of  1988  was
cooler than  in  1971 (see Table 4.2).  The O'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 year-around, not  only in  summer  as  in  this  case
study.
                                     64

-------
         Table 4.5   Cumulative heat balance  components  [1000  kcal m
-21
Hsn
Lake Calhoun
Ha Han
He
He
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
-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
o>
en
                                                     Year 1988
APR
MAY
JUN
JUL
AUG
SEP
OCT
134
306
495
675
822
915
998
173
385
610
853
1088
1291
1463
-212
-462
-732
-1016
-1295
-1542
-1775
-39
-116
-248
-374
-490
-559
-€20
5
15
9
10
4
-4
-29
60
128
134
149
128
101
38
-35
-97
-225
-354
-477
-553
-«11
75
176
194
210
183
147
102
-51
-147
-286
-414
-529
-595
-647
34
72
70
84
68
48
11

-------
                      SIMULATION PERIOD 4/01-10/31
   800000-
   70COOO-
   70COOO-
5 60COOO-
^.
t/)
o 500000-
_J
U}
>
J 4000OO-
3
G.
3J
a
3
o
  300000-
  200000-
   10C-000-
               1971 losses (kcol/m2)
               1988 Josses (kcol/m2)
               1971 losses (m)
               1988 losses (m)
                                         LAKE CALHOUN
                i       t       I
               1971 losses (iccol/m2)
               1988 losses (kcol/m^)
               1971 losses (m)
               1988 losses (m)
  70C-OOC-
  60C001
               1971 losses (kccl/m2)
               1938 losses (kcci/rn2)
               1971 losses (m)
               198S losses (m)
                                              HOLLAND
F
-1.0
^0.9 ~
                                                                 "°-8 y,
                                                                      .0
                                                                  0-7 ~!
                                                                 -0.5 2
                                                                 -0.4  5
                                                                      2
                                                                 -0.3  3
                                                                 -0.2

                                                                 -0.1
              A»R
                     MAY
                            JUN
     Fig. 4.4     Cumulative  evaporative losses  (simulated).
                                     66

-------
o>
                        O
                        o
                        LJ
                            30-
                            25-
                            20-
                            15-
                            10-
                        UJ
                        2    5H
                        UJ
                            -5-
                           -10	r
                                     1988-89
                                     1971-72
                                  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 mev ~ 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
                                   m
     Surface  mixed layer depths  are shovra 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  by wind  and   natural  convection.
Surface mixed layers were  deeper in  Lake Elmo, mainly  because more wind
energy  was  available for  entrainment  at the therrnocline 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, and  1  to 2  m  deeper mixed layers  were
probably produced in this way.

     Fall  overturn occurred  earlier  after  the  warmer 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 thermocline.

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  morphonaetries,   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  ~  3°C   compared  to  3°C  in  air  temperature
change) and lower after  the lake started cooling.

     Daily  hypolimnetic temperatures are  shown in Fig.  4.8.   Values are at
depths   well   below  the thermocline,  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.2° C  lower  under
warmer  conditions.     Once  a   stable  stratification   was  established,  the
hypolimnetic temperature  was almost constant throughout  the  summer for all
three lakes.
                    s
     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  one—dimensional 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
warmer 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 —
  1

40-
    H
  30-

§25-]
X   -4
| 20-j
                     LAKE CALHOUN
1971
1988
                                                .(.'  -i
                                            •mi     q
    1
  15-
    3
    j
               Ifl
                    LAKE ELMO
                                          1971
                                          198S
  40-
    5
  35-
                    LAKE HOLLAND
                                         1971
                                         198=
                     JUN   JUL   AUG   SEP    CC'
      Fig. 4.6    Mixed  layer depths (simulated).
                              71

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






p'
O3 —
30-

25-

20-

	 1971
	 1988
	 DIFFERENCE
(1988-1971)

• /
i
i ' 1 ( i ' i ' i ' i
LAKE CALHOUN -j
H
i
.t'\,,jt" \tf~: .;..•': - \ ~]
!' /'s-'v Vv \ f'~-~ ? .../-A _ j
A./ '"*-' ^ ' \ • J
; v '""\ j
   -5-
 30-
3
<
a.
5
    0-
 -5-


-10


 30-


 25-


 20-
           1971
           1988
           DIFFERENCE
           1971
           1938
           DIFFERENCE
                                       LAKE ELMO
                                                  \
LAKE HOLLAND '
          APR   MAY   JUN    JUL    AUG   SEP
     Fig. 4.7     Simulated epilimnion temperatures.
                             72

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

   2-
     •-••	  1971
     	  1988
     	  DIFFERENCE (1988-1971)
                                      LAKE CALHOUN -
      	  1971
      	  1988
          DIFFERENCE
                                      LAKE ELMO
  -2-
                                     —I	1	j	•—
                                      LAKE HOLLAND
10-
   6-
I  *'

          1971
          1988
          DIFFERENCE
         APR   MAY   JUN   JUL   A'JG   SEP   OCT
    Fig. 4.8     Simulated hypolimnion temperatures.
                                73

-------
                     s
     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 daily-
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 CO3.
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

     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  2XCOj  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

-------
     therefore be  treated herein as either  invariant or be  lowered  to account
 fjr  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.2°C,  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 Airpcrt
(93.13°  longitude,  44.53°  latitude) were  used.   The meteorological  data  £:r
assembled  contains  measured  daily values of  average  air  temperature, de*
point temperature,  precipitation, wind  speed, and solar radiation from 1955  t
1979 (March -  November).  The  period from  1955 to 79 was chosen beca^
it  is long  enough  to  give  a representative average  of  base  conditions  befcrr
climate wanning.   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  foEowing  EPA  guidelines  on  global
 i:mate  change  effect  studies  (Robinson  and  Finkelstein,  1990).    Climate
 projections by  four different models (GISS,  GFDL,  OSU, UMKO) for  the
 loubling of atmospheric COj were provided by  NOAA  (1990).  The monthly
 innate  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.'C)'
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) t
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 = 2XC02 GISS - PAST
t Ratio =  2XC02 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  2XCC-2  GISS  scenario  were used  as inputs  to the  water
temperature simulations.

                                    77

-------
                                • MINNLAPULIS/Sl.PAUL
                                ©GiSS GRID POINT MINNEAPOLIS/St.PAUL
                                • DULUTH
                                  IGISS GRID POINT DULUTH
oo
                          Fig. 5.2     Geographical location of the  closest  GISS  grid points for
                                     Minneapolis/St. Paul and Duluth.

-------
-4
co
                              o:
                              O
                              oo
700-



600-'



500:



400



300-



200:



100-



  0-


100-
                                      	 MSTP (1955-1979)

                                      	 2XC02 CLIMATE SCENARIO (GISS - model)
                                      MAR  APR  MAY  JUN  JUL  AUG  SEP  OCT  NOV

                                  80-
                              9   60-
                              ID
                              X
                              LU
                              UJ
                              a:
                                  40-
 20-
     	 MSTP (1955-1979)

     	 2xC02 CLIMATE SCENARIO (GISS - model)

     A
•- MSTP (1955 1979)

--•• ZXCCb CLIMATE SCENARIO (CISS - model)
                                                 UJ
                                                 K
                                                 ID


                                                 £
                                                 LU
                                                 CL
                                                 s
                                                 LU
                                                 I—

                                                 CK

                                                 <
                                                                                        -p^-T-r-r-r-l-. I '1 I •» I I I <  I I I  I ) I r I I I I  I I I I  I I

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

   2xco2 CLIMATE SCENARIO (Giss - model)
                                                                                              /. ti;; •
                                                                                              im
                                                                                       M*R  APR  MAY  JUM   JUL   AUG   SEP   OCT   NOV
                                             -10


                                             30
                                                                                                                                    -25
                                                                                                                                     20
                                                                                                   -15 Q
                                                                                                       ui
                                                                                                       tu
                                                                                                       a.
                                                                                                    10 to
                                                                                                                                    -5
                                  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 climatological  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,
ajid  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
oxvgei  balance.   Trophic status was assessed by using  a Secchi depth scale
 Eeistary and  Wilson,  1988)   related  to  Carlson's  Trophic  State   Index
        , 1977).   Secchi depth information was available in the  MLFD.
     A  statistical  analysis of southern  and northern Minnesota Lakes  in  the
MLF2 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
an 3  5ecchi 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
Lais - 5
< 5
5-20
> 20
< 1.8
1.8 - 4.5
>4.5
- Cumulative
Frequency
Lower 30%
Central 60%
Upper 10%
Lower 30%
Central 60%
Upper 10%
Lower 20-50%
Central 20-50%
Upper 0-10%
Class
0.2
1.7
10
4
13
24
1.2
2.5
4.5
Description
Value
Small
Medium
Large
Shallow
Medium
Deep
Eutrophic
Mesotrophic
Oligotrophic
                                      80

-------
        J\
    1 1 •-'-.•• Jj"' --•
    - '.I-'..-3 ^jVa.
                                  ' 5<.3ECCH:DEPTH<.4.5MET=?S
                                                             SECCHICe-TH>45METERS :
                                        _-

                                           TV/
                                   <= DEPTH <» 20 McTERS
                                                             AfiEA>5.0SQKU.
Fig. 5.4       Geographic  distribution   of  lakes   according  to  key   parameters:
               Secchi depth,  maximum  depth, and surface  area.
                                       81

-------
                  100
                                MAXIMUM DEPTH (m)
                  100
                            2468

                                SECCHI DEPTH (m)
                                                         10
Fig.  5.5      Cumulative  distributions  (%)  of  key  parameters   in  Mirmesc
           •  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

-------
   UJ
   cr
   <
   Ld

   <
   UJ
   rv
1.0

0.9

0.8'

0.7

0.6-

0.5'

0.4-

0.3-

0.2-

0.1-

0
                            SMALL LAKES
G-O Carver
S-© Crystal
    Deep
GHD Farquhar
E-ffl Island
    McDonough
t-A Long-N
    Lucy
O-O Parkers
    Re'itz
    Sucker
    Twin-Middle
    Wolsfeld
           Area=1.l9-exp(-1.76»Depth)-0.20
                                             Auburn East
                                             Auburn West
                                             Bald Eagle
                                         5-B Boss
                                         S-ffl Bavoria
                                             Big Carnelion
                                             George
                                             Caihoun
                                             Crystal
                                         C— O Elmo
                                             Harriet
                                             Eagle
                            MEDIUM  LAKES c
           Area= 1.14»exp(-2.1 *Depth) -0.15
                             ARGE LAKLS   3-0 Waconia
                                                •—-€» Big Marine  ~
                                                    Coon
                                                3-3 White Bear
                                                EH3 Greenwood
           Areo = 1.14«exp(-2.91»Dep
         0.0   0.1   0.2   0.3   0.4  0.5  0.6   0.7   0.8   0.9   1.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  monthly
corrections  specified by the 2XCO2 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 5
to 7°C, and occurs in spring (April), the minimum  is  on the  order of 0 to
2° C  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.2°C,  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.   Fjctreme 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 4° C.

     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.5° C  colder after  climate change than before.
Eypolimnetic warming  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 warming patterns in response to vertical diffusion.


                                     85

-------
00
o>
                                            EPlUWNiON ICMP


                                       SiUUl>HOH I'HIUU  19^5
   30-



   25-



   20-



   15



   10



   5-







   30



g 25

ul
a: 20

\—


a 1S
a.

S 10
I—


   5



   0



   30






   20



   15



   10



   5



   0

                                  MEDIUM LAKES
                                  LARCt LAKFS
                                 ****

         EPIUMNION TEMPERATURES

    3XCO2 CUMATT StrrVfliO (CISS - mod-M)




SMALL LAKES
                                                                            MEDIUM LAKES
                                                                            lARCt LAKES
                                                                                                   tPHIMNON ICWPtHAIURES


                                                                                                 DIFrERCNCE (CISS mod«l - PA5T)
                                                                                                                                                MEDIUM LAKES
                                                                                                                                                LARGE LAKES
                                                                                                                                                                a
                                                                                                                                                                2
                                                                                                                                                             -a
                                                                                                                                                              -i
                                 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.7        Simulated  weekly  epilimnion  temperatures.

-------
oo
                                           WPOUMNtON TtMPEfUIURCS
                                       SIMULATION PRlOO 1955 - 1979 (PAST)
     HYPOLIMNION IEMPER»TURES
2XCOj CUMA1E SCt>WW (OSS - modtl)
  KTPOUUNION rEMP

DIFFIRENCE (CISS mod.l - PAST)
                                                                                                                                                     — 5
                                 MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY  JUN JUL AUG SFP  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
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
Epilimnion

"C
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
maximum temperature

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
occur
Hypolimnion

°G
24.9
26.8
27.0
26.2
27.0
27.1
26.5
26.9
27.1
11.9
12.8
17.5
18.7
19.9
23.0
24.0
24.6
25.5
7.3
7.4
7.8
11.6
11.8
12.6
18.2
18.4
19.4


day
206
204
203
204
203
203
203
203
203
278
277
261
254
252
233
220
218
211
308
308
308
294
293
291
261
263
259

GISS-2XCO2
Epilimnion Hypolimnion

"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 °C
217 26.2
217 28.3
217 29.2
217 27.5
217 29.1
217 29.4
217 28.3
217 29.1
217 29.4
217 12.6
218 13.0
218 17.5
218 18.2
218 20.3
218 25.3
223 26.0
218 26.6
218 27.3
217 10.3
218 10.4
220 10.6
218 12.8
223 12.9
223 13.3
223 18.4
223 18.7
218 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 (GISS-PAST)
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

Hypolimnion

°C
1.3
1.5
2.2
1.3
2.1
f.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


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

-------
               EPt/UNON TTUP£RATURES EUTROPHIC LAKES

                     f*AST 1955-1979
NION TEMPERATURES EUTROPMIC LAKES
 PAST 1955-19?9
                                              depth medium (13

                                              ores smoii (0.2 wm
         MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JUN J'JL AUG SEP OCT NOV
Fig. 5.9       Examples    giving    range   of    epilimnetic    and   hypolimne;
               temperatures over a  25 year period  (95%  confidence  interval).
                                          90

-------
CO
                          Fig.S.lOa
Simulated temperature  (isotherm'J structure in three medium deep  (13  m  maximum  depth)  lakes of
large (10 km2), medium (1.7 km')  and small (02 km2) surface area. Isotherm hands are in increments
of 2'C.  Simulated water  temperatures  are for  past climate (1955-70)  (top)  and the  VXCOj O1SS
climate scenario (bottom).

-------
CO
10
                             Fie.S.lOb   Simulated  temperature  (isotherm)  structure  in  three  medium  area  (1.7  km3)  lakes  of  different
                                        maximum depths,  shallow (•! m),  medium  (13 m) and deep  (2<1  in).  Isotherm  bands are in  increment
                                        of 2'C.  Simulated water temperatures  are  for  past  climate  (1065 70)  (top) and  the SXCO? G1SS
                                        climate scenario (bottom).

-------
                                                                                 MEDIUM ABC* DtCP CUmOPMIC IAXCS


                                                                                      SIMULATION YfAS 1970
to
CO
                                                                                                  	AIMOSFMRtC IONS WAV!'
                                                         APR  MAY  JUN   JUL  AUO  SEP   OCT   NOV
                                                                                                  APR  Mst   JUH  JUL  AUO   S£f>  CXil  HOV
                                                              t ' ' ' I ' '  ' i ' ' ' i
                                                         APR  MAt  JUN   JUL  MiG  SCP   OCT  NOV
                                                                                                  APR  HAY   JUN  JUL  AUG   S£P  OCT  NOV
                                                        APR  MAY  JUM   JUL  AIJG
                                                                                                  APR  MAY   JUM  JUl   AUG   SEP  OCT  NOV
                                             Fig.  5.11       Examples  of  individual  surface  heat  flux  components.

-------
CO
                                       cuMuunvt NFT HOT nt;x
                                    SIMULATION PRIOO 'ois - 1979 (PAST)
                                                                                 CUUULA1M Nrt HEAT O.UX
          CUMULATIVE NET HtAT FLUX
        OKTEREHCE (GISS model - PAST)
                       200000-
                       400000-
                       200000-
                       400000-
                                                                                                                                                           150000
                                                                                                                                                          -t 00000
                                                                             2XCO, CUUAIE SCENARIO (GISS - model)
                                                                                                                                               *
                              MAR APR MAY JUN JUL
MAR APR MAY JUN  JUL AUG SEP OCT NOV
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 aU  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 are  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 spring  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  ELA.    In the  ELA  mixed  layer  depths
increased due to  transparency increase and  increased  winds  due to reduced
forest cover resulting from increased incidence of forest fires.


                                     95

-------
                    1QS5 - 1979 (PAST)
   MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JJN JU. AUG SEP OCT NOV
Fig. 5.13      Simulated cumulative evaporative losses.
                            96

-------
           SIMULATION PERIOD lass - 1979 (PAST)
                                            2XCOj CUWATE SCENARIO (CISS - rrxxj«l)
    MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY j'jN JU_ AUG SEP OCT NOV
Fig.  5.13      Simulated cumulative evaporative  losses.
                               96

-------
       MIXED LAYER DEPTHS
 2XCO-? CLIMATE SCENARIO (GISS - mode!)
        MIXED L*ttK DEPTHS
SIMULATION PERIOD 1955 - 1979 (=AST>
                                                           oeep oiigotropnic
                                                               cutrophsc
                                                        X—X medium ofigotrophicH j»
                                                        -f—-1- rr>«diurn eutropn'tc
,.RAPR MAYJUN JUL AUGSEP OCTNCV     APR MAYJUNJULAUGSEPOC7NOV
        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)                            ,    N
                       Ln = 	372	                    (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), zm is maximum lake depth (m),
zg is the height of the center of volume of lake, A0 is lake surface area (m2),

pQ 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 1°C.   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 1°C 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

-------
  100
   80-
   60-
<;
                                  eutropWc PAST (1955-1979)
                               Q-Q shallow od'gotrophis PAST (1955-197
                                  shallow eutrophic 2XC02CISS
                               X-K sfwuow ougotrophic 2XCO2OSS
                                             25      30
Fig.  5.15      Simulated lake numbers  as a
                status.
                                    Dn  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)







DHKP
•?1 <»i







krn2
SMALL
(0.2)

MEDIUM
(1.7)

LARGE
(10.0)

SMALL
(0.2)

MEDIUM
(1.7)

LARGE
(10.0)

SMALL
(0.2)

wriMUM


, * »


TROPHIC
STATUS

Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
( )li(.'otrophic
Kulrophic
* ' . , • . - v .

• - .. .


BSS
day
118
134
0
132
0
0
133
0
0
100
100
101
105
106
106
106
106
124
101
101
101
104
!H1

•1


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
? ' * '.'
™ t
»

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
1HH
!'-9


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


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
ef ! !
*
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
'*'">
G1SS
ESS
day
271
246
0
255
138
0
255
137
0
288
287
276
274
270
251
250
250
247
301
302
302
290
29U
290
275
275
•}T\
1 "--'" ' "- --" " 1"— g -— "• - -' ": "" ,n i —,,M. -—-„..,....,„ ^.M,.-.III..LL||I%,L
-2xCO2
LSS
day
204
162
0
180
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
xs~ 	 •-;;•'— ±';!sgas::;
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
tS
11


15


_6
-3
8
12
14
10
9
17
37
-11
-11
-10
-5
-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

-------
 •SS      Beginning seasonal  stratification, i.e. first  Julian day when difference
         between  surface and deep water temperature is greater than  1"C.
                               m
•SS      End  seasonal  stratification,  i.e.  last  Julian   day  when  difference
         between  surface and deep temperature is less than 1°C.

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:

     (l)  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.4°C  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  epiliranetic 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.7°C;  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 th-;
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


     -As  a  result of the  research described here,  a  better  understanding  of
  *  freshwater  inland  lakes  respond to variable atmospheric  conditions has
  •_  rained.

      Chapter  2  describes  how a specific lake water  temperature model was
  icrjJized to simulate the seasonal  (spring to fall) temperature  stratification
  er   a   wide  range  of  lake  morphometries   and  meteorological  conditions.
  c±l coefficients related  to  hypolimnetic  eddy  diffusivity,  light attenuation,
  ni  sheltering  and  convective heat  transfer were generalized using theoretical
      empirical   model   extensions.      The  proposed  regional  lake   water
'.'•rrr-erature model simulates  the onset  of  stratification,  mixed  layer  depth,
ir,d i=rater 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
hypclimnion, 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
hypclimnetic 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  of 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-435.

Aldama,  A.A,  D.R.F.   Harleman,  and  E.E.  Adams  (1989).    Hypoiimnetic
      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,  AD—SB-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
      thermocline 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,  ELS.  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, EYl, 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,  F.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 ana  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.  (1987J.  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.

Eorsch,  G.  M.   and H.   G.  Stefan  (1988).    Convective  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.
                m
 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., Lemmin,  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.p.,  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.
       S91-S92.

 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 it*
       annual heat  budget for  the sediments in two Wisconsin lakes. Limnolo-
       and Oceanography, Vol. 14(1), pp.  115-135.

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

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

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

McLaughlin, D. (1985). - A  distributed  parameter  state  space  approach
       evaluating the  accuracy of the groundwater predictions,  Ph.D.  t:--
       Dep. of Civ.  Eng., Princeton  Univ., Princeton,  N.  J., 1985.
                                     108

-------
   ~-i-   FLO.,  W.S.  Combs,   Jr.,  P.O.  Smith,  and  A.S.  Knoll  (1979).
      :--xtenuation  of light and  daily integral rates  of photosynthesis  attained
      -T   planktonic  algae,  Limnology  and  Oceanography,  Vol.  24(6),  pp.
      1 33S-1Q50.

   -1er1  J.D.,  J.L.  Goddier,  H.A.  Regier,  B.J.  Shuter,  and  W.J.  Christie
      _~9S7).    An assessment  of  the  effects  of  climate  warming on  Great
      -•akes  Basin  fishes,  Journal  of Great Lakes  Research, Vol.  13(3), pp.
      MQ-352.

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

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

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

           U.P.,  Schindler,   P.W.,   Wirz,  U.E.  and   Imboden,   D.M.  (1983).
       Chemical  and  geochemical studies of  Lake Biel  II.  A chemical approach
       to lake mixing. Schweiz. Z. Hydrol, Vol. 45(1),  pp. 45-61.

Os rood,  R.A. (1984).   A  1984  study of the water quality  of 43 metropolitan
       area lakes, Metropolitan Council,  St. Paul, MN, Publ. No. 10-84-172, p.
       40.

Osgood,  R.A. (1989).   An  evaluation of  the  effects  of watershed treatment
       systems  on  summertime  phosphorus concentration  in metropolitan area
       lakes,  Metropolitan  Council, St. Paul, MN,  Publ. No.  590-89-4)62, p. 85.

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

Patankar, S. V. (1988). Numerical  heat transfer and fluid  flow,  McGraw-Hill,
       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  hypolimnetic    nitrogen  - and   phosphorus
       transformations in  Lake Rotoiti,  New  Zealand: A geothennally influenced
       lake,  Limnology and Oceanography, Vol. 31(4), pp. 812-831.

 Protopapas,  A.  L.  (1988).  Stochastic  hydrologic  analysis of  soil-crop-climate
       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.
             m
Quay,  P.D.,  W.S.  Broecker,  R.H.  Hesslein,  and  D.W.  ScMndler  (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 Fails  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,  D.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  developmec:
       of climate scenarios for  impact assessment:  Phase  I  final report,  USEPA
       Atmospheric    Research    and    Exposure    Assessment   Laboraton
       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 uni-
       stable atmospheric conditions, Boundary-Layer Meteorology, Vol.  44, ;:
       171-180.

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

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

Shapiro, J.  and H. Pfannkuch  (1973).   The  Minneapolis  chain of  iak
       study  of urban  drainage and its effects 1971-1973,  Limnological
       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

-------
  ..r   •'•: H..  J.  Imberger,  and  K.N.  Rayner  (1986).   Modeling  the diurnal
        .i-i  layer,  Limnology and Oceanography, 31, pp. 533—556.
                                     s
   ,  ••-  K.E., and  Armstrong, D.E.  (1983).  Lake  mixing and  its relationship
          epilimnetic  phosphorus in  Shagawa  lake,  Minnesota,   Can.  J.  Fish.
       \^c:. Sci., Vol.  41(1), pp. 57-69.

    »..-  H.G.  and D.E.  Ford (1975).   Temperature  dynamics in dimictic  lakes,
       Journal of the Hydraulics Division,  ASCE, Vol. lOl(HYl),  pp. 97-114.

    i-   E.G.,  M.J.  Hanson,  D.E.  Ford,  and   S.   Dhamotharan  (1980a).
       Stratification  and  water   quality   predictions  in  shallow  lakes   and
       reservoirs,  Proc.  Second Int'l Symp. on Stratified  Flows,  International
       Association for Hydraulic Research, pp. 1033-1043.

    .::.  H.G.,  J.   Gulliver,   M.G.  Hahn,  and  A.Y.  Fu  (1980b).    Water
       temperature  dynamics   in  experimental   field  channels:  Analysis   and
       modeling,   University   of   Minnesota,   St.   Anthony  FaUs  Hydraulic
       Laboratory, Report No.  193.

  :rfin,  E.G.,  S.   Dhamotharan,  and  F.R.  Schiebe  (1982).    Temperature/
       sediment model for a shallow lake,  Journal of Environmental Engineering
       Division, ASCE,  Vol. 108(EE4), pp.  750-765.

 -..fan, H.G.,  J.J.  Cardoni,  F.R. Schiebe,  and C.M.  Cooper  (1983).   Model
       of light penetration in  a  turbid lake,  Water Resources  Research,  Vol.
       19(1), pp. 109-120.

 Suub,  P.  T.  and T.  M.   Powell  (1987).    The  exchange  coefficients  for
      latent  and sensible heat flux over  lakes:   dependence upon atmospheric
       stability, Boundary Layer Meteorology, 40, 349-361.

 Sweers, H.E.  (1976).   A  monogram to  estimate the  heat-exchange  coefficient
       at the air-water  interface  as  a function  of wind speed and temperature;
       a critical  survey of some  literature,  Journal  of Hydrology, Vol. 30,  pp.
       375-401.

Townley, L. R., and  J.  L.  Wilson (1985).  Computationally efficient algorithms
       for   parameter estimation  and  uncertainty   propagation  in  numerical
       models  of groundwater  flow, Water Resour. Res., 21(12), 1851-1860.

US  Army  Corps  of   Engineers  (1986).   CE-QUAL-Rl:  A  numerical  one-
       dimensional  model of reservoir water quality user's manual,  Waterways
       Experiment   Station,  Corps of Engineers, Vicksburg, Mississippi, Report
       E-82-1, 1986.

Waggoner,  P. E.  (1990).   Climate Change and  U.S.  Water Resources, Wiley
       series in climate  and the biosphere, John Wiley  & Sons. Inc., pp.  496.

Wanner,  H.  and  U. Siegenthaler (1988).   Long  and  short  term variability
       of   climate,    Lecture   Notes   in  Earth   Sciences,  Springer-Verlag
       Berlin/Heidelberg, pp. 175.
                                    Ill

-------
Wdander,  P.  (1968). Theoretical  forms for the vertical exchange coefficients in
       a  stratified fluid  with  application  to  lakes  and  seas,  Acta  Regiae
       Societatis  Scientiarum  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).    Ground-water  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.   (Kz) was calculated  over  time intervals of 1  to  15  days and varied
from  10"3  to 10"1  cm2s"1.    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 =  c(N2)  .   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  (1m) used in this  study  were found to  give  the
best Kz estimation.
A.I   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 hydrothermal  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  Earleman  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 0.01* C.   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

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

     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:

                                                   -i.s                (n i\
                                                   +  S                (A.I)


The  flux-gradient method for  the  computation of Kz reduces  this  equation to
the form


            K                                     "
                =  ( § j"1 { If     T(^C +  wT|"  - /  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).
                                   z           rr      rr
            v     f <9T 1—1 /  d   /* rp//-\j/.     sed     sol(z) 1         /A  q"
            Kz = 1 si j    t  -at J  T^dc  - WP --- pep )         (
Vertical kinematic thermal eddy diffusivity can be explicitly  expressed as:
z
ft [ JT(0 dC "
0
:
+ w T
c
z
- j S dz
i 0
                Kz = 	2	$	S	            ^4
                                     5T
where  T(z,t)  =  measured  water  temperature  distribution,  z  =  upwarz
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  a:
which  Kz is computed,  and without  significant inflow  or outflow to or fro-
the lake,  net vertical velocity w is customarily small enough to  be  neglecte_
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-

-------
radiation  absorbed in the water  and heat  flux through  the water column to
>r 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

— ~  7T  §   '  P  —  density of water and  g  =  acceleration of gravity.   If
      (7Z  p
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  Sts 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  <3Ts/#t  =  Kzs(52Ts/&2)  is  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  za  =  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
                                         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.088°C.   The  square of
the correlation coefficient (R2) was  0.98.

     The  unsteady  sediment  heat  conduction  model  was  calibrated  fcr

thermal diffusivity.  A sediment thermal  diffusivity of 0.0022 cm2sec"1 applied
uniformly  over  depth.  (Gu  and  Stefan   1990)  simulated  measured  sedimer;
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 th~-
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  differer.: •:
was  less than 0.1  °C.   Since  sediment thermal  properties do  not  change  wr.i
season,  the calibrated  values should hold for year—round  simulation.

     Using the  calibrated model and measured  deep  water  temperatures  frc~
the 1989  data as the  boundary condition  at the water  sediment interface,  t::
sediment  temperatures  for  the  period of lake  eddy diffusivity  studies (May
to August 9,  1989) could  be simulated.   Heat flux  from  the" sediments wa~
calculated using the simulated sediment temperatures in equation  (2).  Liker.

and  Johnson (1969) used values of  ps=l.2 g cm'3 and  cps=0.8 cal g'^C"1 ::'

soft  bottom  sediments,  giving (ps  CpS) =  0.96  cal cm"3 "C"1.   A  water :
solids  ratio of 3:1  corresponded  to  the  calibrated  thermal diffusivity  (Carsli *

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



                                    A-6

-------
      6
O
o
 (D
 ^_
 3
-4—'
 O

 0)
 Q_

 E
 (D
 0.5m Above
-& Sediment
  Surface
-0.5m Below

- 1 .Dm Below

M .5m Below
             i     |     ,     |    ,     j     ,     j_ _,__   ^  . _j    ^    ,     ^     r_.

      April  2    April  16    April 30   May  14    May 28   June  11   June 25    July  9
       Fig. A.2      Temperatures  recorded  in  lake  sediments  and  overlying   water,
                     1990.

-------
>
                        6.5-
                        6.0-
                   o
                   ^  5.5-
                   U
                   D
                   Ld
                   (1
                        4.5-
                        4.0-
0,5 m
— calculated
	 measured
	measured
	 measured
                          April  2
       April  10      April 18      April 26       May 4       May 12
                                       MmiMr.l and measured temperatures in lake sediments, 1990.

-------
                                    accuracy
equation  (2) is  estimated  to  be 2 kcal  nr2day"1.
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 thermocline.  '  The  rise in  water temperature  at
the  4  m depth  in July was due  to epilimnetic  warming  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  surface
water  temperature  from May 7 to 17  coincides  with high  sclar 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,;  at
increments  Az  and At in depth and time,  respectively, Equation  (A. 4) must
be  discretized numerically to yield the eddy diffusion  coefficient estimate:  in
the  form:
                v  _

                                       231

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,   Esed   = water/sediment
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

-------
h-«
o
                 o
                 o
                  a;
O
v,

CL

F
0)
                                                                                     -7 meters
                                                                                     -8
                                                                                       9 motors
                                                                                     MO molec
                                       i#wm#^
      Moy /     Moy 7 I     June -1    June IH   July ,!    July )G
                                                                            July 30   AIKJUS!  I '>
                        Fig, A.4     Hypolimnetic lake  water  temperatures  recorded at 2 min interval
                                    in  Ryan Lake, May 7 to  August 9, 1989.

-------
800
           solar radiation
           1  m depth water temperoture
           air  temperature
  May 7     May 21    June 4    June 18    July 2    July 16    July JO   August 13
Fig. A. 5      Meteorological  conditions  (solar  radiation,  wind
               tomperatuieH)  dining tlic  period of  investigation.
                                                                        and   air

-------
 o.

Q
      May  7    Jun  4    Jul 2    Jul 30   Aug  27   Sep  7A   Oc(  22
       Fig.  A.6     Seasonal lake  temperature structure. Isotherms (°C) shown in this
                   figure arc derived  from measurements at 20 minute and 1 m
                   «l« i.lh inlet viil'l

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

=  o(N2)^ .    As expected,  an inverse relationship  between Kz  and  N2 was
observed.  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

-------
Q
fN
_J

<

O
x
Z)
_J
u_

I—


                                                                            LJ
                                                                            I—
         May 7    May 21   June 4   June 18   July 2    July  16   July 30  August
        Fig. A. 7     Heat  flux  through  the  sediment-water   interface  and  solar

                    shortwave radiation received at the 4m  depth, May 7 to August
                    ! i

-------
I-*
Cn
                               10°-q
                         I
                         U
                         OJ
                        (N
                         F
                                                      RYAN LAKE
                                          TURBULENT DIFFUSION COEFFICIENT TIME SERIES
                                  30     May 20      Jun 9
)ul 19      Aug 8
                          Fig. A.8     Vertical turbulent diffusion coefficient time series  in Ryan Lake.

-------
       10
  o
  Q)
  C/5
  O
  O
  O
  CO
 0>
 GO
 N
 p
 o


 isi"
                   I Jf^K^Z^OXICT'CN2)-0-3* :
                 R3=0.32 . T=1 day •
                 SEDIMENT HEAT NOTCo
                 '   ~ ~-"L	""'""" ""	"   ~"--'-
                                    ^=0.55 , T=1y       .
                                    SEDIMENT HEAT CONSlDEHffl  tf
                         K =1.15X10"'(N2)
                 R'=0.67 , T=5 days*
                 SEDIMENT HEAT NOT CONSIDERED
                                    R3=0.79 . T=5 days
                                    SEDIMENT HEAT CONSIDERED
R =0-87 , T=10 doys
SEDIMENT HEAT NOT CONSIDERED
                                                     R'=0.92 . T-10 days
                                                     SEDIMENT HDa CONSIDERED
                         K =1.87X10~'(N2)
R3=0.90 . T=15 doys        *
SEDIMENT HEAT NOT CONSIDERED
                                                                  ooys
                                                     SEDIMENT HEAT CONSIDERED
                 io-6  io-5
                        N2 (sec-2)
Fig. A.9      Calculated vertical  eddy diffusion  coefficients for time intervals c;
               1,  5, 10, and  15 days,  with  and  without sediment  heat flux.
                                    A-16

-------
    i;:le  A.I   Regression  coefficients  for Kz = a(N2)
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
-Q.43±0.018
Hsed = 0
a 7
2.20*10-5
1.15*10-4
1.47*10-4
1.87*10-4
-0.36±0.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
 in estimate of an instantaneous  value.  Variability of the eddy  diffusivity was
 :r.e highest for  a  sampling  interval of one  day.   Different regression lines
 rjuld   be  fitted  to  the  one day  results without   changing  the  regression
 :oe£ficient  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.
 Table A.I.
Regression coefficients  with standard errors are summarized  in
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 day'1 is about  30%  of the average
net heat flux value of 7 kcal

  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

                              e
                              8
0.01 (°C)


0.05 (°C)

0.01 (m)

 2.0 (kcalnr'day*1)
                                    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:
                              K2 = K z + ckz                         (A.7)

                              ATt = AT;  +  £t                         (A.8)

                              Az = Az  + eAZ                         (A.9)
                                S = 5 +  6S   -                       (A.ll)


with the statistical properties

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


                         E(e2) =  
-------
                             N   N
     VAR(KZ)=
                  ATz   f AZ  ,  ^IA   - '   N
                                 Az
                                    4
At
       T,AZj  - JL 1
          1 fi    N    2
       * -^ £   4  eAz ATi I                                     (A.15)
        AtAz3 i = 1   L  ^


where  i or  j designate the depth under consideration, Az =  depth increment,
ATi  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  2oiz  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

-------



A
-p
H^
Q.
Ld
Q


J
4-
5-
6-
7-
li
9-
m-
~"™~ rnson
... Jun 9, 1989 ."."." ^
T
%.
'•>:<^.
\ ''••-
\ "'""

                                                                                 3-

                                                                                 4

                                                                                 5

                                                                             1   6
                                                                             i
                                                                             H-
                                                                             Q.   7
                                                                             UJ
                                                                             o
                                                                                 8

                                                                                 9
                               0.000  0.010  0.020  0.030  0.040  0.050  0.060

                                                     a
                                                    Jun 29,  1989
                                                                           -2a
                                               Kz(cmasec-')
                                        10
                                         0.000  0.010  0.020  0.030  0.040  0.050  0.060

                                                         K2(cm2sec-')
K3
O
3-

4-


5-
                            I
                            Q.   7 .
                            u
                            O    .

                                8-


                                9-
                               10-
                               0.000  0
                                              9, 1989
                                         5-

                                     i  e-i
                                     I
                                     a  7-
                                     Ld    '
                                     o
                                         8-
010  0.020  0030  0.040  0.050  0.060

         ?(crn3sec~')
                                                       19. 1989
mean
+ 2o
-2ff
                                                10
                                                0.000  0.010  0.020  O.OJO  0.040  0,050  0.060

                                                                 K^cm'sec-')
                                   Fig.  A. 10     Mean   values  plus  or  minus  two   standard  deviations   of  eddy
                                                 diffusion  coefficient as a function of  depth.

-------
 a
 CO

 £
 o
•b
CN
     0.070
     0.016-
0.012-
     0.008-
     0.004-
     0.000-
                        Starting date:  Jul  19, 1989
                                                  1  '  !  '  1  f
                                                  •  A m depth

                                                  •  5 m

                                                  -  6 m

                                                  -  7 m
                                                    	8 m

                                                    	 9 m
                                   —,._„____,
                           5   6   7   8   D   10  1 1   12131415
      1   2   3
                          TIME  INCREMENT  (DAYS)
         Fig.  A.11    Estimated  eddy diffusion errors (2cTj,- ) for  different  sampling

                    intervals.

-------
to
                  3 no
                  2.60
              2!
              ^   2.20
              Cu
              (r
              cj

              ?   1.80
              I __

              Q_
              UJ

              Q
1.40
                   1.00
                 8g
         00060° d o o 6
                                      4   5   6   7.8
                                        1  '  i  '  I  '  I
                                        9   10  11  12
!3  14  15
                                         TIME  INCREMENT (DAYS)


                         Fig. A.12    Estimated eddy diffusion errors 2crK (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  cm2/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 = o(N2)^
produced the  best fit when  a =  0.00014 and 7  =  — 0.43 for periods of  15

days,  where Kz is in cm2s'1 and  N  in  s'1 .  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  m^day1.  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  (1.1)  is  discretized using  the  control volume
method.   For intermediate  control volumes (i  = 2, m-1).

           A  .    k        k+i
                      i
                                   V0.5   k
                                   —   - Ki+0.5
          A        k        k+1     (Az)2  k    (Az)2H
        < "if Ki+0.5 ) T1+1  = — T, +  _—              (B.1,



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  K 1+0-5

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 = - LJ=±-                      (B.2)
                                     Ki-l+ Ki
                                    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   8K.      _k+i     dK-n,  -k+l -i    k
         i—0.51      1—0.5 ri—0.5  rp      T'   4-
                                         i .    I  J. . -  "T~
         A.  L  gi*      i-i     ^k      i   J   i-i
           1        1—1              1—1
       A.     r  dK. ...  .k+l     5K.  ...  .k+l -i    k
         I—0.5 I     1-0.5 rp           1-0.5  rr,      T /    I
               	—	  _[_ _    — 	_	  j_ _      J. .   T
         A.   L  ffTk      1~1     5Tk      l    J   l
          .+05r  5K.+Q5  .k+1    ^Kj+QS  "k+1 1   (Az)2        1    k
          	—  	r-r——  T.	a-r——  T.     +	h &Had  T'   +
           A  L  frrk      1+1    ^T^      l    J    At      l   I   l
                5K..n-  -k+1     SK... .k+l  n    k
                   1+0.5 rp	1—0.5 rp       rp/
         A
         *X .     wj-.,            ^/ -*. .  -
           i        i+l               i+l
where flj =  l if i = l else Si  = 0


     For  the  surface control  volume matrix  B(k) will  slightly differ.   First,
terms which are multiplied by ~Ti are  the  first  entry  (bu).  Secondly,  eddy
diffusivity K 1-0.5 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  aHbr    SHc       3He i  Az    k
               Had  =   —«- + --j- + —^   	 T;           (C.I)
                      L  5Tk     OTj       5Tk J pwcp  l


     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
Ki|o-5 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.
           b                                              dp(T]   g
K =  a  (N2)    where P is Brunt -Vaisala frequency P =( - ) -3-
                                                           dz     p

                    dK...          K2         5K.  5N.
                    - L    = 2 - 1 -- 1 — i               (C.2)
                      5T.        (K.  ,  +  K.)    SN. 5T.
                        i        v  i— 1      v      11
              M
     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(k):
where air perturbation  temperature vector F (k) is:



                                                  TI
                                                     ai


                                          f,
dew  point temperature  perturbation vector F (k) is:
                               5Ha    5HC    Az     k
                      Hcl =  (  —- + _ )  	 T'a.                 (C.3)
                                                  k
                                                 T'd.                  (C.4)
wind speed perturbation (mxm) matrix F (k) is:
        A-_ncf   5K -_nK  -k+!    5K •  n=  -k+i
         1	0.5 I       1—0.5 rp	1—0.5 rp
          A.  L   5\\rsk  ,     i-1    5\\%k
           1
        A.  n_r   5K ...  .k+l    5K .  . _  .k+i
         1 -0.5 I       1—0.5 rp	1-0.5  rp
          A.  L   5Wsk      i~1    5Wsk      »
            1         1                 1
                                    A-26

-------
           A.
              -r  dK ..„_  .k+l     3K  .  ..  .k+l -,          v      k
             •5       J+P-5  T            1+0.5 rp       ,  CTT   I  \TCT I   _L
             —  	*  .  f   1.  .	«—r	  1.      + o.iicS   VVS.   +
          A.
                                           -
  where
                                                Az     k

                      HC3 =  ( —   + - r )  - Ws'                (C.5)
     First  and  last  control  volume  have  interface  diffusivities  K.   _  and


K.   ,  equal to zero,  and the additional term HC3 is present  only in the  first
  l"~"Tj. tJ

control  volume.




     Solar  radiation perturbation vector F (k) is



                                 93. sn
                                                 H's.                    (C.6)
                                    A-27

-------
                                APPENDIX D

                            Cross-term 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
                     loo          o

  S(n) =        °         atd^   °          °
                 000          0

                 000          0
Pta(n,k)
Ptdta(n.k)
0
0
Ptatd(n.k)
Ptd(n,k)
0
0
0
0
0
0
                     100         0

                 0         0"td(k)   0         0

  S«=        0000

                 0         0        0.0
where  ata is air temperature standard deviation,  
-------
      If meteorological  perturbations  are not cross-correlated covariance matrix
Muc(k,k)  is

                     0-ta2(k) 0 *       0         0

                     0      (7td2(k)   0         0

        Muc(k,k) = 0      0         aws2(k)   0


                     00         0         asr2(k)
                                      A-29

-------
               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
Secchi depth reading
Number of days for output
Dates for output
Heat budget and mixing coefficient
Initial stage & Outflow channel
Initial conditions

MBOT NM NPRINT INFLOW DY MONTH ISTART MYEAR
Z(1)...Z(MBOT)
T(1)...T(MBOT)
Field data section

NF NPRFLE
NFLD
DEPTH(1)...DEPTH(NF)
FDATA( l).-.FDATA(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. 12.  15. 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 12.3  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.
14 1
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." 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
7 1
1
0. 12. 13. 14. 15. 20. 23.
12.7  12.7 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
DBL
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
Emissivity 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

-------
12345*789
SLARGE
 * PROGRAM REGIONAL

C
C     THIS PROGRAM IS MODIFIED VERSION OF THE MINLAKE
C     MODEL (RILEY & STEFAN, 1987, UNIVERSITY OF MINNESOTA
C     SAFHL. PROJECT REPORT * 263). 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,TAlR(31).TDEWt31).RAD(31),CR(31).W]ND(31).
   + PR{31).DRCTpl)
   COMMON/RESULT.T2('W).CHLATCrr(«).PA2(-W),PTSUM(«)),BODa40).
   +PC2(3.40),XNC<3,«),T2MMON/VOLXMAXD:«40).Z<40).A(40).V(40),TV(40}.ATOP(41).DBL
   COMMONSUESDZt"60).SZ<60),LAY(40).AVGI(4,60),SVOL(<>0)
   COMMON/CHOiCE-MODEl_NrrRO.IPRNT(6),NDAYS.NPRNT(30V.NCLASS.PLT{30)
   COMMON*' ATER.3ETA.EMlSSJCKLXKiHK.MAXWCOEF.WSTR
   COMMON/CHAN EL'WCHANL.ELCB^LPHA.flW.WLAKE
   COMMON/STE?S.D211_DZUUMBOT.NM.NPR!NT..VDAY.MOYrH.;LAY.DY
   COMMON,'STAT,S-JMXY(K)',.SUMX(10),SUMY(10)JCSO.:i01.
   + YSCX10).RSQ(10).RMS{10).R£LM(10).MTHREL(10),MDAYREUtO'..2REU10-,.
   + ZREL.M{!0).RS(W-.REmoi.MTHR.MS{]0).MDArR.MS!10),ZRS{50!.ZR.MS(!(r,
   COMMON/INFLOW QIN(5).TINr5).PAiN(5i,BODiN(5},3O!N(5).aN:S ;.
   *CDINrS)_XNH!N(5:,XNOIN:S),CHLAiN(3J)
   COMMON,SOURC£RM(3.40),PROD('IOVXMR(3.'10),PRODSUM(4us
   COMMON/FIELD IFLAG(10).FLDATA(ia50).DEPTH;SO).NFLD;;0;
   COMMON/FILE D:N.MET.FLO.TAPES.TAPEL1REC
   COMMON.TITL' TTTLE
   DIMENSION B{!05^TATVAR(SOUCX(4)
   CHARACTER-t>4 TTTLE
   CHARACTER' 14 DIN.METPLO.TAPE&TAPEl
   CHARACTER'! TS(W).FF.CXl,CXiCX3.CXOCS
   EQUIVALENCE (STATVAR(1)3UMXY(1)X(TAPE&TS(1))
   EQUIVALENCE(CXLlCX(l}),(C5C2.ICX(:)).(CX3.!CX{3n.(CX4.ICXN))
   DATA IPAN.PCOEFF/0,0.«/
   DATA ICX/16 »ia 14*5B. 16*32. 16*4A/
   FF=CHAR(12)
  »0 FORMAT(1X4A1)
   YSCHL-m
   HSCSI=0.03
   CONST».5
   CHLMAX-ao
   DO 995 1=1.4
  995  IPRNT(I)-0
  9*i  !FLAG(I)=0
    DO 997 1-180
  99?  STATVAR(I)=aO
    WRJTEC.lOOl)
    READCVCA)") DIN
    WRITEC.1CIOO)
    READ(','(A)T MET
    WRJTE('.10(C)   "
    READ ('.'(AJT TAPES
    DO405!=L16
      11=14-1 + 1
      IF(T8dl).NE' 0 THEN
        T8(H •(•!)='.•
        T8«I+4).T
        GOTO 406
      ENDIF
 405   CONTINUE
 406 OPEN (7.F1LE-DINJ
    OPEN (8,FILE«TAPE8)
    OPEN (9.FILE=MET)
 C
 C THESE ARE OUTPUT DATA FILES
 C
    OPEN(21.F!LE=-£?!L.PRNT
    OPENC-FILE-WPOLPRNT
    OPEN(2S,F!LE = TEM PER.PRN1)
 C
    READ(7.'(A)') TTTLE
                        A-35

-------
 505  FORMATf REQUEST CHANGE OF TITLE (YV) ^.S
 506  FORMATC ENTER NEW TITLE :")
     WRTTEC.505)
     READ(V(A)-) XS
    IF(XS.EQ."T .OR. XSEQ.VT THEN
     READ(V(A)') TTTLE
    ENDtF
    WRTTE(8,W) FF
    WRJTE(&1900)
    \VR1TE(8,'(A)T TITLE
  9» FORMAT(LXA1)
  1001 FORMATC' ENTER INPUT DATA FILE NAME :MOX/)
  1002 FORMATC ENTER FILE SAME FOR TABULAR OUTPUT .  V)
  1000 FORMATC ENTER METEOROLOGICAL DATA FILE NAME .• v
  1003 FORMATC ENTER INFLOW DATA FILE NAME :-.«.V}
 c
    CALLSTART(STS.FT.ISTART..'NFLOW.MYEAR.!RUN.!LSEKI)
 C
 C**"* Call LAKE routine to change or add input quantities "*

    ZMAX=ST-DBL
    CALL SETUP
    1-1
  11 IF(DZ(!}.GT.DZUL_AND.MBQT.LT.«>: THEN
     CALL SPUT(U'-AY)
     GO TO 11
    EN'DIF
    i-l+1
    IF(LGT.MBOT.OR-I.GT.40) GOTO 12
    GOTO 11
  12CALLSETZ(MBOT)
    CALL VOLUME(MBOT)
    CALL AREA
    CALL ABOUND
    CALL TVOL(MBOT)
 C..DETERM1NAT1ON OF INITIAL MLXED LATCR DEPTH
    ILAY-1
    DO7!-tMBOT-l
    !F(T2(I),GT.T2(Itl) + CONST) GO TO 8
  7 ILAY-ILAY+1
  8 DM!X=Z(ILAY)+DZ{lLAYpOJ
    DO 9 1»1.10
    RELM(I)-a
  9  RMS(I)=ao
    NDAYS-1
    MP-0
   KDAYS=0
   IPRNT(5)=1PRNT(5)-1
   MDAY-0
   CALLPTABLE
   1PRNT(5)=IPR,NT(S) + 1
C— Stan limubuo^ for each month
   YEAR=R£AL(MYEAR)
   IF(AMOD{YE\R.4.0).EO.aO) EDAY-Vs^
   ETSUM=a
   HTSUM-a
c
   DO100J-UNM
   CALL MTHDATA(MONTRKDA^•S_V^i"EAR)
   EDAY-365.
   YEAR-REAUMYEAR)
   IF(AMODO:'EAR.40).EO.a) EDAY«3».
C-START SIMULATION FOR EACH DAY
   DO 200 MDAY-START.KDAYS
   NFLOW-. INFLOW
   CALL LAKE(0_fl,O.S)
   IF(MDAY.EO.KDAYS.OR.MP/KPRIXT-NPR!NT,F,O.MP> IPRNT(1 >
   IF(MONTH*100+MDAYiO.KPRNTOOAYS))!PRNT(l.«l
C
   P=.PR(MDAY)'0-0254
   MP=MP+1
   TM1X-T2(1)
C-CALCULATION OF KINETIC ENERGY FROM WIND STRESS
   CALL WNENtTAUVCWINTXMDAY))
   RKE-TAU*VC'ATOP(n
-------
c
C-BETERMtNATION OF WIND MIXING ORDER
C-HEAT IS ABSORBED FIRST. THEN WATER COLUMN IS MIXED BY THE WIND
C
   CALL HEBUG(ILAY.TMXQNET.HS.HA.HBR.HE.HC
   +TAIR(MDAY),TDEW(MDAY),CR(MDAY),RAn(MDAY).W]ND(MDAY).
   + IFAN.PCOEFF.SEKI)
   CALL CONMIX(ILAY.TMDLMBOT)
C..CALCULATION OF EVAP. IN TERMS OF VOLUME
C..CALCULATES LATENT HEAT OF VAPORIZATION ALV
   HEV.HED-ATOP(l)
   CALL WINMtX(RK£.TMK.ILAY,MBOT)
   DMDC= Z(ILAY) +OJ«DZflLAY)
C
 35  CALL LAKEfCLCLftO)
c
    IF(IPRNT(l).EQ.l) THEN
C-Ourcui ubuUr data on Upc&OAT
    CALLPTABLE
CL-Output meteorolopcal data on upeS-DAT
    CALLPMETE;HS.RAD;MDAr).HA.WlND(MDAY").HBR,
   + P.HE.TAJR(MDAY),HCTDEW(MDAY).HED.HEV,QNET.DMIX.ZE'JPH.SEXI)
   ENDiF
C.. Output to plot Be (apc&PLT)
   IF(MDAY + MONTHMOO.EO.NPRNT(NDAYS))THEN
C ___ Access and output Held data
     CALL FDATA;NF)
     WRTTEC.3020;
     R£AD{-,99) XS
 9»  FORMAT(Al)
 3020  FORMATC/,".' DO YO'J WANT TO \1EW GRAPHICAL RESULTS  Y'NV.
   +     3X\)
C,..Caii plotting routines
     !F{XS.EO.'Y' ^OR. XS-EO.y) THEN
      CALL SUBPLOT(NF.MYEAR)
     ENDiF
     ENDIF
   DY=DY+t
   IPRNT(1)=0
 200 CONTINUE
   ISTART=1
 100 CONTINUE
C-Coccpute and list statistics:
C~   l)retative and absolute m^dmum deviations bem-een
C_    model and Heid data and day of occurrence
C_   2)  slope of regression of Reid data on simulation results
C—   3)  regression coefficient
C~   4)  standard error of the regression
   WRJTE<8,3000)
   DO 101 1=1,10
    rF(XSQ(!).GT. 0.) B(1)=SUMXY(I)/XSQ(I)
    IF(iFLAG(i).GTj) THEN
      RSQfl)-L-SUMXY(l)'
      SUMXY(I) = SORT(SUMX Y(I) )
     ENDIF
 101 CONTINUE
    WRJTE(3>3001)
    WRITE(&3013)(RELM(J)J= 1.10)
    WRlTE<&30MXMTHREL(J).MDAYREL(Jp-tlO)
    WRJTE(a3015)(!lMS{J)J-tlO)
    WR!TE
-------
 3016 FORMATflX'DATE OF MAX. ABS. ERR. M4X10(4X. 12,7.n))
 3017 FORMAT(1X.•SLOPE; MODEL TO DATA REGRESS1ON',3X10(2,X.F7.2))
 3C18 FORMAT(1X'REGRESSION COEFFICIENT (R"2)',7X.10(2.X.F"7.2))
 3019 rt>RMAT(lX/WLAKE WATER TEMPERATURE STATISTICS')
 3024 FORMAT(LX.-STANDARD ERROR OF ESnMATEMOX1002X.F7.3))
   END
C
   SUBROUTINE FDATA(NF)
C*"***
C**"*" Subroutine to read field data from the input data
C***** and compute statistics and deviations benvcen Held
C*"** data and simulation
C*****
   COMMON/FILE/DIN.MET.FLO.TAPE&TAPEUREC
   COMMON/RESULT/T2(40).CHUVrOT{40).PA2('tO).FTSUM(40).BOD2(40X
   + DSO2(*).C2;40),CD2(«),XNO2(40).XNH2(40).CHLA2a'>0),
   + PC2(3,40S.XNC2(3.'W),T20(40).S!2<40)
   COMMON,STAT,SUMXY(10).SUMX(10)5UMY(10),XSO(10X
   +YSO(1C).RSQ(!0).RMS<10).REL-MC10).MTHREL<]0).MDAVREU10;,ZRELC10},
   H-ZR£LM(10XRS(10).REL(10),MTHRMS{10),MDAMlMS(10).ZRS{10i.ZRMS(10)
   COMMON,VOL,'ZMAX.DZ<40}.24'40)J\f«XV(40).TV{40)J\TOPi«).DBL
   COMMON/STEPS/DZLL.DZUUMBOT.NV.NPR1NT.MDAY.MONTH.ILAY.DY
   COMMOKCHOICE'MODEL,KrmOjPRKT(6),NDAYS.NPRNT(30),NCLASS.PLOT{»)
   COMMON./FIELD/IFLAG(10).FlDATA(ia50).DEFTH(50).NTU>{10)
   DIMENSION COMP(4aiO)
   EQUIVALENCE (T2(1).COMP(U))
   CHARACTER'!* DIN.MET.FLO.TAPEa.TAPE!
   DO 303 lm 1,10
     RS(1)=0
 303  REUO-0
   DO304J-I.20
     DO 304 1-J.10
 304  FLDATA(LJ)=aO
   READ(7,')NF.NPRFLE
   NDAYS-NDAYS+1
   IF(NF.GT.O) THEN
     READ (7.')(NFU>a).i = l.NPRFLE)
     READ(7.')(DEPTH(I).I« 1.NF)
     DO 305 I-1NPRFLE
     READ{7.')(FLDATA(NFLD{I)J)J-1.NF)
 305   CONTINUE
     LL-1
C~.Loa»e stmuliiJon values corresponding to sampled
C—constituenu »nd depth of fveld data
     DO 310KS=1,NF
     L=LL
     DO 315 LL-L.MBOT
       ZX»ZZ(LL-l)V(Z
-------
    WRlTE(a3015) (RS(1).I-UO)
    WRITE(8,3016) (ZRS(I).] = UO)
C—S&re data on plot file (upeiPLT)
    IFnPRNT(5).GT.O) THEN
      WRfTE(l,REC-IR£C) REAL(NF)
      IREC-IREC + 1
      WRITE<1.REC=1R£C) REAUNPRFLE)
      IREC-IREC»1
      DOSOOl-tNF
       WRITE<1.R£C=IREC) DEPTH(I)
 500     1REC«IREC+1
      DO 501 1-LNPRFLE
       WRJTE(l.REC»iREC) REAUNFIXKO)
 501     IREC=IREC+1
      DO502I2-LNPRFLE
      DO502I-1.NF
       WRTTB'l.REC.iREC) FLDATA(NFU)(!2).I)
 502     1REC.IREC + 1
      IF( N'FLD(l).NE.l) THEN
        x«ao
      ELSE
C_.Mfct£
-------
      IF(NrTRO-GT.O) THEN
      ELSE
      ENDIF
      %TUTE(S.3004)
      !F(NITRO.GT.O) THEN
      ELSE
      ENDIF
    ENDIF  .
 1200 FORMAT(///5X"ZOOPLANKTON PARA.METERSVJX2K'-
   + -COSC (#.'/M3)1.rC.'PREDATION{*/MJ)>J.\.X3RA2ING(MG'M3)1.3X.
   + ' DAYDEFTH(M)V.5>XFZO,9XF7.1.9XF&4.1LXF4.1.0
 1201 FORMATfFlU)
 2990 FOR,MAT(//SX1NmAL CONDmONSVSXlSC-OO
 3000 FORMAT(,'.5XWFORMAT!ON ON LAKE QUALrTYV5X27f-V)
 3001 FORMAT(5X"Z(M)  T(C) SS(PPM) TDS(PPM)  CHLA(PPM)  PC(PPM)',
   + ' P(PPM)  PT(PPM)  BOD(MG/L) DO(MG'L) DZ(M)-.
   + ' AREA(M2)   VOUM3}")
 30K:FORMAT(5X3!F5.:.:X).lXF5.1,4X4{2XF6.4),6XF&i6X
   + F5.2.4XF4.2.2(2XE10.4))
 3003 FORMAT(5X'Z,'M)  T(C) SS(PPM) TDS(PPM)  PA(PPM)  PT(P?M)'.
   + ' BOEXMG.1-) DO(MG,T_) NO^MG,!.) NH4(.MG/L) DZ(MV,
   + ' AREA(M2)   VOUMJ)1)
 3004FORMAT(SX3(F5.1ZX),lXF5.1,2<«.Ft3).4XF&lSXF5.2,4X
 3005 FORMAT(SX3(F5.iiX),lXFS.t2(«.Fi3).«.Fil5XF5.12SXF4.1
   +2(1XE10.4))
 3006 FORMATC.VSX31OLOGICAL PARAMETERS'. SX21('-VJX-Z(My.9X
   +'CHLOROPHYLL'.6XTOTAL1.6X1NERNAL PM IX "STERNAL NVWX.
                  \2X))
 300? FORMAT(5XF5.12X10fF5.3,2X))
 3008 FORMAT(;,'.SX' TEMPERATURE PROF1LESV.5X1Z(M)-,«.T(C)\
   +6X'AREA (M2)'.SXVOL (M3)1)
 30W FORMAT(5X2(F5.2.2X)A'2XEia4))
   RETURN
   END
   SUBROUTINE SUBPLOT(NF.MYEAR)
^••••**
C"*** Produce profile of slate variables and
C"*" Ma dau (if available)
C*****
   CHARACTER"! XS.CHtCH2.CH3.CH4
   CHARACTER'S: TTTLEtl)
   CHARACTER'4 MNTH(U)
   COMMON/SOURCE/R,MC3,40).PROD(40).XMR(3.40)J'RODSUM(«3)
   CO.M.MON,STEPS-DZLL,DZUl,MBOT.SM,NPRiNT.SroAY.MONTK.ILAY.DY
   COMMON/RESULT/ T2( 40 ),CHLATOT,40).PA;<40),PTSUM(40>.BOD:/ 40),
   + DS02(40).C2(40).CD2{40;j(NO2(40)_XNH2(40),CHLA2(3.40).
   + PC2{3,40)JCNC2(3.40),T3>: 40)512(40)
   COMMON/nELD'lFLAGilO).FLDATA(10.50).DEPTH(50).Nn_D(10)
   COMMON/VOLZMAXDZ(40).Z( 40)^(40). V(40',.T\'!401,ATOP(4!).DBL
   DIMENSION  ZD(23),FDfr").\TMf43).ZV(43i.VAR(4«,10)
   DIMENSION  FCT(10).LENM(l)(.
   EQUrVALENCE(CHl.lC\(l)),(CHllCX(2)),(CH3.ICX;3)).(CH4.!CXi-r,}
   DATA !CX'l<»*I3,1A#58.k.*31]&*4A>'
   DATA !SCAL'!,3'-!.3'l.3'-l/
   DATA FCT/u%Kioo,4'i.:riooo'
   DATA MNTHTJAN '.'FEB VMAR VAPR '.'MAY '.'J'JN VJUL ".'AUG '.
   ->• SEP •.•OCT '."NOV '.'DEC V
   DATA TITLE .TEMPERATURE (C)7
   DATA LEN n&,27.2a25.23.4'24.25/
 1000 •ft'RrTEC.109) CHtCHlCH3.CH4
C_itst and ieiecl desired state variable for plotting
   DO 1001 = 1,1
 100  WRJTE(MOl) !.Tm-EiO
 101  FORMATflXir - 'A32)
   WRTTEC.99)
 99 FORMATC/IXXZHOOSE (1) DESIRED PLOT (ENTER  O TO QUIT) :
   READ(V.ERR-lOOl) 1C1
   X»(X
C>..cban|;e deptb to negative vatues for pkxune with depth
   DO110!=LMBOT
   ZV(1)»-Z(I)
   VTM(I)=VARaiCl)
   IF(XLT,VAR(L!C1)) THEN
     X-VAR(UCl)
     ISO-1C1
   ENDIF
                        A-40

-------
 110 CONTINUE
   J2-°
C...if field data i» available, locale field data corresponding
C_to selected state vanabies
   IF(NF.GT.O) THEN
    DO 111 I-INF
    IF(FLDATA(ICU).GT.O) THEN
     12-12+1
   .  FD(I2)=FLDATA(!CLI)
     ZD(I21--DEPTH(!)
     IF(X.LT.FD(12)) THEN
      X-FD(I2)
      IND--I
     ENDIF
    ENDIF
 111  CONTINUE
   ENDIF
   NPEN-1
   NPEN1-1
   !PORT»99
   MODL=99
   ZV(MBOT+1)» THEN
     D-MBOT+1
     CALL SCALE(ZV..,m_25,'. '.tt^)
    CALL NUMBER(SW,W9_25.XMYEAR_0,-1)
    CALL NEWPEN(NPENl)
    CALL PLOT (L.7.,-3)
 C~p!o< simulated proTiies as a line
    CALL UNE(VTM.ZV,MBOT.l,a,l)
 C-piot field dau with a symboJ
    IF(llGT.O) CALL UNE(FD.ZD,!2.1,"M)
    CALL ?LOT{Oi,0-.999)
    WR1TE(*.UO)
   130 FORMAT(1X'SEND TO HARDCOPY DEVICE ? (Y/N)  '.\)
     IFCXS-EQ.T- .OR. XS.EO.y) THEN


                        A-41

-------
     READ-;'.*) IPORT.MODL
  140  FORMAT(/1X'ENTER PLOT8S IOPORT AND MODEL : \\)
     WWTEC.143)
     READ;',') NPEN.NPEN1
  143  FORMAT(/1X 'ENTER UNE WEIGHT (AXIS.DATA) :  -.\)
     WRJTEC.141)
     READ;'.") FCTOR
  141  FORMATC/1X 'ESTER REDUCTION FACTOR { >10 ) : '.\)
     GOTO WO
    ENDIF
    GOTO 1000
 1001 RETURN
    END
    SUBROUTINE ABOUND
£*•••••
C**"* Compete! the surface area of each layer (ATOP)
C***"* using ibc depth area retaubnsbip in LAKE.
C""*    ATOP(l) « surface area of the late
C""*    ATOP(MBOT+1) = HO,
c*****
    COMMON.STEPS.DZU.DZUL.MBOT.N.M.NPRlST.MDAY.MONTHJijkV.DY
    COMMO.VVOL'Z.MAXDZ;;'W).Z(40).A(40),V{«5).T%'(''0),ATOP(4S>.DBL
    DUM-tt
    DO 100I = 1.MBOT
    ZDUM = ZMAX-DUM
    CALL LAKE(ZEiUM.ADUM.O,l)
 100 ATOP(!}=ADUM
    ATOP(MBOT+1)-0.
    RETURN
    END
    SUBROUTINE AREA
C******
C*"** Compute the area through the tniddk of each layer
C*"**" us*ng ibc depth-area relationship in LAKE
£-*•**»
    COMMON,STEP&OZLLDZUL,MBOT.NM.NPRlST.MDAY.MONTH.r_A'i'.3Y
    COMMOSA'Ob7-MAXDZ(40).Z(40)^(40).V(40i.TV(4»).ATOP(41 J.D3L
    DO1001-1.MBOT
    ZDUM-ZMAX-DUM-DZ(!)A
    CALL LAKE(ZDUMJ>kDUM.O,l)
 100 A(I)-ADUM
    RETURN
    END
    SUBROUTINE COEF(MODEI_MBOT.NCLASS)
C*****
C***** Compere some coefficient* used in the constant
C****** volume and finite difference solutions
C**"**
     COMMON/COEFT/ DUM^40).DUM3(40)
     COMMON/VOL' ZMAXDZCW).2440).A!40),V(40).TV(40).ATOP(41-s.DB.
     COMMON/'FLOW/HMK(41),OE('tO),F\'CHLA(5].PE(5.4I}
    DO !«H*1MBOT-1
     D'JM1<= l,'CA(I)'DZ
-------
   DIMENSION RHOT(40)
   DO100I-LMBOT
 100 RHorci)-RHOVCO.MB
                       A-43

-------
      1F(MODEL,E0.3) THEN
       DO57K-1.NCLASS
 57     PC2(KJl)-(rX^K,n)'V01
      'ENOIF
     PA2(l[)=(PA2CIl)'V(H)+PA2ADUM(40)^ID(40)
   COMMON.'STEPS/DZli.DZUUMBOT.NhtNPRINT.MDAV.MONTri.ILAY.DY
   COM.MON/FLOW/HMK(41).QE(-W)J:VCHLA(5).P£(5,41)
   COMMONAVATER/BETA.EMISSJCKl_XKZHKMAXWCOEF.WSrR
   DIMENSION Q{40)
C..CALCULAT1ON OF THE HEAT ABSORPTION FROM METEOROLOGICAL PARAMETERS
C_!N A COLUMN OF WATER
C-.CALCULATION OF HEAT FLUXES INTO THE WATER BODY
   CALL FLXIN(HS.HA.TAIR.RAD.CR.C2)
   CALLFLXOUT(TS.HBR.H&HCTAIR.TDEW.W1ND.IPAN.PCOEFF.WCOEF)
   HOOUT-HBR+HE+HC
C..CALCULATION OF EXTINCTION COEFF. (ETA) AS A FUNCTION OF SUSPENDED
C...SEDIMENT CONCENTRATION
C..CALCULATION OF HEAT ABSORBED IN EACH LAYER
     HQ-a-BETA/HS
   EX=EXP(-ETA'DZ{1))
                       A-44

-------
C-CONVERSION FACTOR OF 1000 USED FOR DENSITY-HEAT CAPACITY OF WATER
   HQpHQ'EX
C-CALCULATE THE SOURCE TERM Q FOR EACH LAYER
   DO 10U2..MBOT
    ETA-LWa-SEKJ)
    EX-EXP(-ETA'DZ(I))
    HQ=HO'EX
 10 CONTINUE
   CALL SETAMKfWIND,HKMAXlLAY.,MBOT)
C-SET-UP COEFFICIENTS FOR TRI-D1AGONAL MATRIX
   DO 100 I-2.MBOT-1
    Dl- 2.'(A
   BK(1)=L-CK(1)
   I « MBOT
   CALL SOLVE(T1MBOT)
   TS =T2(1)
   CALL CONMIX(ILTS,MBOT)
   DO"»I=LIL
   T2{I)=TS
 90 CONTINUE
C
C KEEP 4oC WATER TEMPERATURE
C
   !Ffnu).LT.4.) THEN
   DO 94671 -1..MBOT
   T2(I)=4
 »4«7 CONTINUE
   ENDIF
C
   ON=HS+HA-HQOUT
   RETURN
   END
   SUBROUTINE FLXINtHSNJUN.TORAD.CCC)
C-CALCULAT1ON OF THE TOTAL RADIATION FLUX INTO A BODY OF WATER IN
C-FROM NET SOLAR RADIATION (HSR) AND NET ATMOSPHERIC RADIATION (HSN)
C..IN KCAL-M'M
C-iDso JACKSON FORMULA USED FOR ATM. RADIATION
C-.CONVERT AIR TEMPERATURE IN DEGREE C TO DEGREE ABSOLUTE
   TCA-TC + 273.
   TCA*TCATCA
   + 'CTCATCA)t(L»ai7'CC'CC)
C.CALCULATION OF NET SOLAR RADIATION AND CONVERSION TO KCAL.'M-M'DAY
C-CALCULATION OF REFLECTED SOLAR RADIATION HSR
C.HSRW— DUE TO CLEAR WATER USING KOBERGS CURVES
C-HSRS — DUE TO SUSPENDED SEDIMENTS AT THE WATER SURFACE
   HSR- (a087-^76E-5
-------
    HE=FCN'(ESA-EA)
    GO TO 30
 20 HE-2.54TD*PCOEFF'5
-------
20  CONTINUE
   D0104K-I.KDAYS
  «READ(»,')TAJR(K).TDEW(K).WND(K).DRCTCKVRAD(K1.CRCK;.PR(K')
   PR(K)»PR(Kyioa
104  CONTINUE
   DO 100K-LKDAYS
    TAIR(K) = (TAIR(K)-311-O.SS56
    TDEW(K)-1
c
C MAKE CORRECTIONS FOR CLIMATE MODEL OUTPUT
C IF NEEDED. THESE ARE MONTHLY ADJUSTMENTS.
C
c
c
c
c
c
c
c
c
c
c
*c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
DO12S7K=l.KDAYS
1F(MONTH.EO.3) THEN
TAJR(K)=TA!R(K)+4.SO
TDEW(K)=TDEW(K)»-7.23
WIND(K) = W!ND(K)'aS2
RAD(K} = RAO(K)*1.0*
?R(K) = PR(Kri.02
ENDIF
IrfMONTH.EO.4) THEN
TAIRfKJsTAIRCK") -M 97
TDEW(K)=TOEW(K)"<.-»
WIND = WiND{K)*a85
RAD{K) = RAD(K)*l.O
PR{ K',» PR; K'' 1.17
ENDIF
:F(MONTH.EQ.5) THEN
TA1R(K)=TA:R(K) + U4
TDEW(K', =TDE\VYK'i - 1.%
WIND(K)-W!ND(K)''i.57
RAD; K) = RAD(K)* l.(W
PR(K) = PR(K"0.8
-------
    RETURN
    END
    SUBROUTINE PMETEfHSLRAD.HA.WIND.HBR.P.HE.TAIR.HC.TDEW.HED.HEV.
   + dNET.DMlXZEUPKSECCHI)
C*****
C""' Output able of meteorological and beat flux values
C"" •
    COMMON/FILE/ DIN.MET.FLO.TAPE&TAPEl.IREC
    CHARACTER- 16 D1N.MET.FLO.TAPERTAPE1
    COMMON/STEPS TlZU.DZULMBOT.NM.SPRJNT.MDAY.MONTaiLAY.Dy
C
    WRITE(&2000)HS,RAD.HA.'W!ND.HBR.P.HE.TA1R.HCTDEV,'.HED,HEV,
   + ONET.DM IX ZEUPRSECCHI
C
2000 FORMAT(/.5X26HMETEOROLOGICAL !NFORMATlONw'.5X26 , SXFS.110H DEGREES C
   +/.7X30HEV. HEAT TRANSFER « 3XF7.4.10H M,DAY OR .E10.4.7H M"5-D
   +2HAYJ.TX1WNET HEAT FLUX IN - . 1X.F9.Z9H KCAL^'M.
   +/.7X' MIXED LA'i'ER DEPTH (M)-'.2X.FS.l
   +38X"EUPHOT1C DEPTH (M) ='.^C.F5.2.
   +/8XSECCH1 DEPTH (M) =',7XFS.2)
    RETURN
    END
    SUBROUTINE SETUP
C*****
C*"*** Determine the initial thickness, volume, and area of layer*
C*""* and the toui voiume of above each bycr from the depths given
C"**** in :be input data fiie.
C"****
   COMMONA'OL/ZMAXDZ('IO).Z(40)j\(40).V(40).TVi;40).ATOP(41).DBL
   COMMON,'STEPS:DZLL.DZULMBOT,NM.NPR1NT.MDAY.MONTH.1LAY.DY
   DZ
-------
    SUBROUTINE SOLVE(VARZMBOT)
C— *•
£••••• Tri-diagonal matrix solving routine
C*""
    COMMON/SOLV/ AK040),BK(40),CK(40).DK(40)
    DIMENSION VAR^-WXTXC-W)
    DO6Oi=2..MBOT
     Tr-AK(I)/BK(M)
     BK(I)-BK(I)-CK(H)*TT
 «0   DK(I)=DK([)-DK(l-l)Tr
C"— ••••• BACK SUBSTITUTION*"*"*— •"•
    TX(MBOT)=DK(MBOT)/BK(MBOT)
    D070I-I.MBOT-1
     3=MBOT-I
 70  TX(J)=(DK DZUL) into two or sort
C***** layers of equal thickness. Ail state vambles are !he riffle m
C*"" each new layer. Volume is adjusted later.
C .....
   COMMON/CHOICE/MODEUNITRO.IPRNT(6).NDAYS.Nr>R,ST;30i.NCLASS.PL
   COMMON/RESULT/ T2(40),CHLATOT(40).PA2(-loi.rTS'JV, «n.BOD;!(40),
   COMMON/STEPS/DZLL,DZUUMBOT.NM.NPRINT.MDAY..MOVrH.:LAY.DY
   DO 100 K=!.MBOT
   II-MBOT+I-K
   CD2(I! + !)=CD2(!I)
   DO50KJ=LNCLASS
 50
   DSO2(H+1)=DSO2(II)
   IF(MODEL.EQ.3) THEN
    DO51KI-U
 51   PC2(KL!I + 1)=PC2(KI.H)
    ff(NITRO.EQ.l) THEN
     DO 52 KJ = U
 52   XNC2(K].1I + 1
     XNO2(II-t-l)=XNO2(!l)
    ENDIF
   ENDIF
 100 CONTINUE
   MBOT-MBOT-t-1
   iF(LW.GE-l) LW=LW+1
   RETURN
   END
   SUBROUTINE START(ST,S.FT,ISTART.INFLOW.MYEAR.1RUN.LEN8,SEKI)
C"" •
C*"""* Routine to read the input data file for initial
C***" conditions and input coefficients
C""'
   CX)MMON/RESULT/T2(40),CHLATOT(40),PA2(40).PTSUM{«)).BO02(40).
   »DSO:(40).C2(«I).CD2(*)).XNO2(«)).XNH2(-10),CHLA2(3.*),
   + PCI(3,'«)).X.NO(3.«)).T20(40),S!2(-«))
   COMMON,TEMP/PARIO(
-------
   NITRO=0
   DO 200 I- 1.30
 200 NfRNT(I)-0
C* ........ INPUT MODEL OPTIONS AND INITIAL CONDITIONS ..... *
   REAEK"?.') SEKt
    WRITEr.lOOS}
    READ(V(A)0 X
     IF(XEQ.'Y" .OR. XEQ.y) THEN
     ]PRNT(2)-1
     TAPE1 "TAPES
     T1(LEN8+2)='O'
     Tl(LEN8+3)«ir
     OPEN(ZFILE=TAPELSTATUS»1NEW)
   ENDIF
   !PRNT(4)=0
 1000 FORMATf PLOT FILE TO BE CREATED (Y/N) ?   V)
 1001 FORMAT(Al)
 1004 FORMATC ENTER UP TO 10 DEPTHS TO BE SAVED '.'.
   +'  END WITH A CHARACTER (*,*.#_JC) :  V)
 100* FORMATC OUTFLOW FILE TO BE CREATED (Y/N) ?  V)
   READ(7,') NDAYS
   IF(NDAYS.GT.OJ R£AD(7,')(NPRNT(1).l - 1.NDAYS)
   READ{7,')  DZLL,DZUL,BETA.EM!SS,WCOEF,WSTR
   READfV)  WCHANL.WLAKE.DBLSrS.FT
   R£AD{7,-)  ELCBALPHA.RW
•  READC7-)  MBOT.NM,NPR1NT.!NFLOW.DY.MONTR1START.MYEAR
   READ(7,-)  (Z(!).I»1.MBOT)
   READf?,")  (T2(I).) = 1.MBOT)
 1500 FORMAT( .'.  5X
   +/.5X41HMIN1MUM THICKNESS OF EACH LAYER (DZLL) - J5.2.9H METERfS)
   +.JX41HMAXIMUM THICKNESS OF EACH LAYER CDZUL) * .F5.2.9H METER(S)
   +/.5X40HSURFACE ABSORPTION COEFFICIENT (BETA) - .F5Z
   +/jX3oHEMissrvrrv OF WATER CEMISSJ -  .FS.I
   +WX'EXTlNCnON COEFF. OF WATER (XK1) - '.F5.13.VM"-r.
   +/.5X 'EXTINCTION COEFF. OF CHLA  (XK2) - '.FS.iSX.'L'MG M'.
   +/JXMAX HYPOLIMNETIC DIFFUSrvrTY (HKMAX) - '.F7.4,' M"2,O'.
   +/.SX'WIND FUNCTION COEFFIQENT (WCOEF) - '.F6.3,
   +.'.5X"W1ND SHELTERING COEFFIOENT fWSTR)- \F6JJ.
   +.3XJ4HW1DTH OF INLET CHANNEL (WCHANL) = .F6.19H V.ETERJS))
 1501 FORMATfSXT-ONGmJDlNAL LENGTH  OF LAKE (WLAKEJ -'.FiaiTH METERS.
   +/JX30HDEEPEST BED ELEVATION (DBL! - .F&ilTH METERS ABOVE MSL
   +,'.5XMHINmAL LAKE STAGE (ST) » .FS.117H METERS ABO\"E MSU
   »,'JX1«HBEO SLOPE (S) = .F10.8,
   +;.SX2»HROUGHNESS COEFFIQENT (FT) =  .F44,
   + ^X48HELE\'ATION OF BOTTOM OF OUTFLOW CHANNEL (ELCBi = .F6.1
   + 17H METERS ABOVE MSL.
   +/.5X40HS1DE-SLOPE OF OUTFLOW CHANNEL (ALPHA1 - .F4.2.SH DEGREES,
   +.JXJ1HBOTTOM WIDTH OF CHANNEL (BW) »  .FtlTri METERS,
   +,V.SX34HlNrnAL NUMBER OF LA's'ERS (MBOT)  - .12.
   +..5X37HNUMBER OF MONTHS OF SIMULATION (NM) - .11
   + .3XS4HDAY OF MONTH  OF THE FIRST DAY OF SIMULATION (ISTART) « .::
   +/.5X53HINTERVAL AT WHICH RESULTS V.7LL BE PRINTED ,'SPR!ST> = .13.
   + ?HDAY(S))
C* ......... INPUT PARAMETERS FOR BIOLOGICAL ROUTINES •""
  2  FORMAT(//1X 'ENTER CHANGES: INTEGER.NEW VALUE'..
   + 5X '(ENTER  ANY CHARACTER FOR SO CHANGES):  ;~>
 1950 FORMAT(/.iX 'BIOLOGICAL COErFICIESTSV.iXUC-'i, ,,10X
   + 'CARBON-CHLOROPHYLL RATIO',5XFi.O,.10X'MAX NUTRIENT
   + -SATURATED PHOTOSYNTHET1C RATE'JXFS.3.'  ,T)AYV.!OX
   + 'MINIMUM CELL QUOTA FOR PHOSPHORUS'J.XF5.3..10X
   + 'MORTALITY RATE'_3XF6J,' /DAT/)
 2000 FORMAT;//. JX'BODKM'AX'SBIO'.iX'BRR'.SXTVBODV. 3X
   RETURN
   END
   SUBROUTINE STATS(FLDATAJOUFLAG.DEPTH.l)
C*****
C***** Compute statistics and statistical quantities with
C*"** Y designating field data and X designating model results.
C—"
   COMMON/STAT«UMXY(10).SUMX(10)5UMY(!0)_\SO(10).
   +YSO(10),RSQ(10).RMS(10).RELM(iO).MTHREU10).MDAYREL(10).ZREL(10).
   »ZRELM(10).RS(10).REL(10).MTHRMS(10}.MaA'S'RMS(10j.ZRS(10).ZR.MS{:0)
   COMMONfiTEPS,DZLL.DZUL,MBCT.NM.NPRlNT.MDAY.MONTH.ILAY.DY
   SUMXV(I)«SUMXY{1)+FLDATA*XX
   SUMX(I)"SUMX(1>+XX
   SUMY(I)«SUMY(I) + FLDATA
   XSO(!)=XSQ(I)+.XX
-------
   X3-ABS(XX)
   IF(X2.GT.ABS(REL(!))) THEN
     REUD-X2
     ZftELlI)- DEPTH
   ENDIF
   IF(XlGTJVBS(RELM(i))) THEN
     RELM{I)-X2
     KTIHRELOS-MONTH
     MDAYREL(!)=MDAY
     ZRELO)- DEPTH
   ENDIF
   IF(30.GT.ABS(RSa))) THEN
     RS(I)-XX
     ZRS0)» DEPTH
   ENDIF
   ffCX3.GTABS(RMS(I») THEN
     RMS(t)=XX
     MTHRMS(I ) - MONTH
     MDAYRMS(1)=MDAY
     ZRMS(I)» DEPTH
   ENDIF
   IFLAG=IFLAG + 1
   RETURN
   END
   SUBROUTINE THICKNS{.MBOT)

C**"* Compute tbicicness of each layer from the deptb
C ..... area curve in LAKE.
c .....
   COMMON7vX)L71MAXDZ(4O).Z(40)A(40).V(40),TV(40').ATOP(41).DBL
   AZ=0.
   AVOL=0.
   DO 100N1.MBOT
   II=MBOT+1-I
   AVOL=AVOL+V(II)
   CALL LAKE(ZDUMJ\VOL,0,4)
   DZCII)=ZDUM-AZ
   AZ=AZ+DZ(H)
 100 CONTINUE
   RETURN
   END
   SUBROUTINE TVOL(MBOT)
C"—
C****' Detennioe tbe volume of water above a layer

   COMMONA'OL/ZMAXDZ(40),Z(40XA(40),V(40),TV(40)^TOP(4I).DBL
   SUM=0.
   DO100I-1.MBOT
   SUM=SUMW(I)
 100 CONTINUE
    RETURN
   END
   SUBROUTINE WINENfT.V.WIND)
C_
C_CALCULATION OF THE SHEAR VELOCITY AND THE WIND SHEAR STRESS
C-CONVERSJOS OF WIND SPEED FROM M.P.H. TO M/S
C-DENSITY OF WATER AND AIR ASSUMED TO BE 1000 AND 1.177 KG/M3
C_
   W.WiND-0,447
   CALL LAKE
-------
   + PC2(3,><0).XNC2(3.'!0),T20(40).SI2<.W)
   IFCIL.GE-.MBOT) GO TO 35
   siwfi-o.
   su.M2=a
   DO 10 !• Lit
   1LAY-!
   RV«RHOfnd),C2(l).CD2(!)3'V(I)
   SUM1=SUM]+RV
 10  SUM2=SUM2+R\"ZfRHO(TST£P,C2(l
   »-RHO(TS,C2(!),CD2(i)))
C..CRTTERIA FOR ENTRAINMENT
   1F( E.LT.PE) GO TO 40
C-.E-NTRAINMENT OF L-^VYER 1+1
   1 = 1 + 1
   TS = fTS'TV(i- 1 ) +TSTEP' V(I>iTV(!)
   IFCLGEMBOT) GO TO 40
   RV-RHO(T2(I),C2(I).CD2(!5rV(I)
   S'JM1«SUM1+RV
   SUM2=« SUM2 + R\"Z(1)
   GO TO 20
 35  I-IL
 40  IL-!
   DO50K-UL
 50 T2(K)=TS
   RETURN
   END
   SUBROUTINE VOLUME(MBOT)
C .....
C"*** Compute the volume of each iaycr based on the dcptb-volume
C""* rclauonship found in LAKE
C .....
   COMMON/VOUZMAXDZ('W),Z('W)j\f40).V('10).T\'('lO)J\TOP(«l).DBL
   Azz»a
   CALL LAKE(ZMAXVDUM.0.3)
   VZ^VDUM
   DO100I=1MBOT-1
     Z2=ZMAX-AZZ
     CALL LAKE(Z2,VDUM.0.3)
     V{I)= VZ-VTSUM
 100   VZ=VDUM
   V(MBOT)=VZ
   RETURN
   END
                      A-52

-------
                LAKE SPECIFIC SUBROUTINE
Area computation section
Depth-area functions of the form AREA=f(ZMAX-Z), written as DUM=f(ZD).
AREA(nr), 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)
   COMMON'MTHD,TAIR(31).TDEW(31}.RAD(31).CR(31).WIND{31).
   * PR{3'.;.DRCT(31)
   COMMON RESULT T2(40).CHLATOT(40),PA2{40)
   COMMON,"CHO!CE'MODEL.NrrRO,IPRNT(6).NDAYS,NPRNT(30).NCLASS,PLT(30)
   COMMON;WATEiyBETA,EMiSS.XKl,XiC!.HKMAX.WCOEF,WSTU
   COMMOVCHANEL.'WCHANLELCB.ALPHABW.WLAKE
   COMMONSnEPS,TJZLJ_DZUI-MBOT,NM.NTR]NT.MDAY,MO.NTH,!LAY,DY
   COMMON.STAT,SUMXY(!0>5U.\1X(105.SUMY(10).XSQ(10),
   + YSO(W).RSQ(10),R-MS(10;,REL.M;W).MTHREL;!0;.MDAYREU105.ZS-EL(!0).
   -t-ZR£UM(]0=,RS{10),REL(10).NfrriRMS(10;.MDAYP^MS!10).ZRS(10i.ZR,MS(10)
   COMMONMSFLOW/O!N(5;.'nN{5).PAIN{5;.BOD!S(5).DO!N(5).C!N;5).
   + CD!Nf5;.XNH!N!5).X.NOlN(Si.CHLAIN(3.S)
   COMMON.KELD'IFLAG(;Oi.FLDATA(10.50).DEPTH(50).NFLDC.Oj.SD
   COMMON'TILE'D!N.MET.RJD.TAPES,TAPE1,IREC
   COMMON.TTTL' TITLE
   CHARACTER- 16 D1N.MET.FLO.TAPES.TAPE1
   GOTO(:OO.aO. 300. 400.500, «0,7t«,SOO.- '.i2)
C ..... WRITE EPiLIMNION & HY?O'_!MN'1ON TES'PFR..\T^RES ""
   DO lir.l !«:.MBOT
   IF (1.EOJ) THEN
   WRITE ::ir!) MONTH.MDAY.Z^I).T2(!;.
   ENDIF
   IF (I.EO.22) THEN
   WRITE '21jri) MONTH.MDAY.Z.;!).T2(I)
   ENDIF
 87] FORMAT(1XI4,!.X,K2X,F9.2.3.X.F9.3)
inn CONT.NUE
   RETURN
 600 RETURN
 700 RETURN
 800 RETURN
 <*» RETURN
 1000 RETURN
 1100 RETURN
 1200 CONTINUE
 1300 CONT.NUE
   RETURN-
   END
                                                   A-54

-------