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