University of Minnesota
St. Anthony Falls Hydraulic Laboratory
Project Report No. 329
Water Temperature Characteristics of Lakes
Subjected to Climate Change
by
Midhat Hondzo
and
Heinz G. Stefan
Prepared for
ENVIRONMENTAL RESEARCH LABORATORY
U.S. ENVIRONMENTAL PROTECTION AGENCY
Duluth, Minnesota
August 1992
Minneapolis, Minnesota
-------
Abstract
A deterministic, one dimensional, unsteady lake water temperature
model was modified and validated to simulate the seasonal (spring to fall)
temperature stratification structure over a wide range of lake morphometries,
trophic and meteorological conditions. Model coefficients related to
hypolimnetic eddy diffusivity, light attenuation, wind sheltering, and
convective heat transfer were generalized using theoretical and empirical
extensions.
Propagation of uncertainty in the lake temperature model was studied
using a vector state-space method. The output uncertainty was defined as
the result of deviations of meteorological variables from their mean values.
Surface water temperatures were affected by uncertain meteorological forcing.
Air temperature and dew point temperature fluctuations had significant effects
on lake temperature uncertainty. The method presents a useful alternative
for studying long-term averages and variability of the water temperature
structure in lakes due to variable meteorological forcing.
The lake water temperature model was linked to a daily meteorological
data base to simulate daily water temperature in several specific lakes as well
as 27 lake classes characteristic for the north central US. Case studies of
lake water temperature and stratification response to variable climate were
made in a particularly warm year (1988) and a more normal one (1971). A
regional analysis was conducted for 27 lake classes over a period of
twenty-five years (1955-1979). Output from a global climate model (GISS)
was used to modify the meteorological data base to account for a doubling of
atmospheric COz- The simulations predict that after climate change: 1)
epilimnetic water temperatures will be higher but will increase less than air
temperature, 2) hypolimnetic temperatures in seasonally stratified dimictic
lakes will be largely unchanged and in some cases lower than at present, 3)
evaporative water loss will be increased by as much as 300 mm for the open
water season, 4) onset of stratification will occur earlier and overturn will
occur later in the season, and 5) overall lake stability will become greater in
spring and summer.
The University of Minnesota is committed to the policy that all persons shall have equal
access to its programs, facilities, and employment without regard to race religion, color,
sex, national origin, handicap, age, or veteran status.
-------
Preface
This study addresses the question of how lake water temperatures
respond to climate and climate changes. The study is conducted by model
simulation. The chapters of this study are a collection of papers or
manuscripts previously published or submitted for publication in professional
journals. Each chapter has its own abstract and conclusions. Each chapter
-of this study deals with a subquestion of the problem.
11
-------
Table of Contents
Page No.
Abstract i
Preface ii
List of Figures v
List of Tables ix
Acknowledgements xi
1. INTRODUCTION AND LITERATURE REVIEW 1
1.1 Introduction 1
1.2 Previous temperature prediction model 2
1.2.1 Model formulation 2
1.2.2 Model coefficients 5
1.3 Overview of study 6
2. REGIONAL LAKE WATER TEMPERATURE
SIMULATION MODEL DEVELOPMENT . 8
2.1 Introduction ' 8
2.2 Model generalization 9
2.2.1 Hypolimnetic diffusivity closure 9
2.2.2 Attenuation coefficient 13
2.2.3 Wind sheltering coefficient 13
2.2.4 Wind function coefficient 16
2.3 Water temperature model validation after
generalization of hypolimnetic eddy diffusivity 19
2.4 Numerical uncertainty of model after
hypolimnetic closure 29
2.5 Accuracy of the regional model after
implementation of all changes 30
2.6 Conclusions 36
3. PROPAGATION OF UNCERTAINTY DUE TO
VARIABLE METEOROLOGICAL FORCING
IN LAKE TEMPERATURE, MODELS 37
3.1 Introduction 37
3.2 Numerical model 38
3.3 First and second moment development 39
3.4 Lake Calhoun - application 43
3.4.1 First moment analysis 46
3.4.2 Second moment analysis 46
3.5 Conclusions 54
111
-------
4. CASE STUDIES OF LAKE TEMPERATURE AND
STRATIFICATION RESPONSE TO WARMER CLIMATE 55
4.1 Introduction 55
4.2 Method of lake temperature modeling 56
4.3 Model validation 59
4.4 Results and discussion 59
4.4.1 Thermal energy budget 59
4.4.2 Equilibrium temperatures 64
4.4.3 Vertical mixing/Onset of stratification 69
4.4.4 Water temperatures 69
4.5 Conclusions 70
5. WATER TEMPERATURE CHARACTERISTICS OF MINNESOTA
LAKES SUBJECTED TO CLIMATE CHANGE 74
5.1 Introduction 74
5.2 Method of lake temperature modeling 76
5.3 Climate conditions simulated 76
5.4 Regional lake characteristics 80
5.5 Simulated lake water temperature regimes for
historical and future weather . . . 83
5.5.1 Water temperatures 83
5.5.2 Thermal energy flexes 89
5.5.3 Vertical mixing/stratification/stability 95
5.6 Conclusions 102
6. SUMMARY 103
7. REFERENCES 105
APPENDIX A. Vertical diffusion in small stratified lake:
Data and error analysis A.I
A.I Introduction
A.2 Study site
A.3 Vertical eddy diffusivity
A.4 Sediment heat storage
A.5 Water temperature observation
A.6 Vertical eddy diffusivity estimates A-17
A.7 Error analysis
A.8 Conclusions
APPENDIX B. Temperature equation discretization A-24
APPENDDC C. Temperature equation linearization A-25
APPENDIX D. Cross-term evaluations A.28
APPENDIX E. Regional lake water temperature simulation model A-30
E.I Lake input data file
E.2 Example input data file
E.3 Meteorological data file
E.4 Example meteorological data file E.4
E.5 Program listing
iv
-------
List of Figures
Figure 1.1 Schematic diagram of source and sink terms in the heat budget
model.
Figure 2.1 Hypolimnetic eddy diffusivity dependence on lake surface area.
Figure 2.2 Hypolimnetic eddy diffusivity forcing parameter (a) dependence
on lake surface area.
Figure 2.3 Maximum hypolimnetic eddy diffusivity (at N2=7.5*10'5 sec"2)
dependence on lake surface area.
Figure 2.4 Relationship between total attenuation coefficient and Secchi
disk depth.
Figure 2.5 Wind sheltering coefficient dependence on lake surface area.
Figure 2.6 Wind function coefficient dependence on lake surface area.
Figure 2.7 Lake "wind speed measurements.
Figure 2.8 Ecoregions and spatial distribution of selected lakes.
Figure 2.9 Cumulative distributions (%) of lake parameters in
Minnesota.Lakes selected for model validation are shown by
symbols.
Figure 2.10 Lake Calhoun water temperature profiles.
Figure 2.11 Square Lake water temperature profiles.
Figure 2.12 Waconia Lake water temperature profiles.
Figure 2.13 Thrush Lake water temperature profiles.
Figure 2.14 Williams Lake water temperature profiles.
Figure 2.15 Standard deviations of, estimated lake water temperature
uncertainties.
Figure 2.16 Standard deviations of estimated lake water temperature
uncertainties.
Figure 2.17 Standard deviations of estimated lake water temperature
uncertainties.
-------
Figure 2.18 Simulated temperature (isotherm) structure in Thrush Lake.
Top shows results from validated model and bottom shows
results from regional model.
Figure 3.1 Schematic illustration of the lake temperature perturbation
system.
Figure 3.2 Meteorological variables at Minneapolis/St. Paul. Daily means
and standard deviations for the period 1955-1979.
Figure 3.3 Estimated long-term average epilimnion and hypolimnion
temperatures.
Figure 3.4 Long-term average isotherms in Lake Calhoun.
Figure 3.5 Standard deviations of estimated epilimnion temperature
uncertainties. Contributions by several meteorological variables
and totals are shown.
x
Figure 3.6 Standard deviations of estimated deep water temperature
uncertainties. Contributions by several meteorological variables
and totals are shown.
Figure 3.7 Long-term average temperature profiles plus or minus one
standard deviation in Lake Calhoun.
Figure 3.8 Epilimnion temperature long-term average plus or minus one
standard deviation.
Figure 4.1 Lake Elmo water temperature profiles.
Figure 4.2 Lake Holland water temperature profiles.
Figure 4.3 Lake Calhoun water temperature profiles in 1971.
Figure 4.4 Cumulative evaporative losses (simulated).
Figure 4.5 Mean monthly equilibrium temperatures (simulated).
Figure 4.6 Mixed layer depths (simulated).
Figure 4.7 Simulated epilimnion temperatures.
Figure 4.8 Simulated hypolimnion temperatures.
Figure 5.1 Regional boundaries and geographic distribution of lakes in
MLFD database.
Figure 5.2 Geographical location of the closest GISS grid points for
Minneapolis/St. Paul and Duluth.
VI
-------
Figure 5.3 Climate parameters a Minneapolis/St. Paul in the past and
under a 2xCC>2 (GISS) climate scenario.
Figure 5.4 Geographic distribution of lakes according to key parameters:
Secchi depth, maximum depth, and surface area.
Figure 5.5 Cumulative distributions (%) of key parameters in Minnesota
lakes (from MLFD database).
Figure 5.6 Horizontal area vs. depth relationship for lakes. Area
and depth are normalized.
Figure 5.7 Simulated weekly epilimnion temperatures.
Figure 5.8 Simulated weekly hypolimnion temperatures.
Figure 5.9 Examples giving range of epilimnetic and hypolimnetic
temperatures over a 25 year period (95% confidence interval).
Figure 5.10 Simulated temperature (isotherm) structure in (a) three medium
deep (13 m maximum depth) lakes of large (10 km2), medium
(1.7 km2) and small (0.2 km2) surface area, (b) three medium
size lakes (1.7 km2 surface area) with maximum depths of 4 m,
13 m, and 24 m. Isotherm bands are in increments of 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.
Figure A.7 Heat flux through the sediment-water interface and solar
shortwave radiation received at the 4m depth, May 7 to
August 19, 1989.
Figure A.8 Vertical turbulent diffusion coefficient time series in Ryan Lake.
Figure A.9 Calculated vertical eddy diffusion coefficients for time intervals
of 1, 5, 10, and 15 days, with and without sediment heat flux.
Figure A.10 Mean values plus or minus two standard deviations of eddy
diffusion coefficient as a function of depth.
Figure A.11 Estimated eddy diffusion errors (2ffv ) for different sampling
•J^-Z
intervals.
Figure A.12 Estimated eddy diffusion errors 20w (cmV1) space-time
tradeoff, 7 m depth July 19.
Vlll
-------
List of Tables
Table 1.1 List of calibration coefficients with ranges used in previous
simulations.
Table 2.1 Quantitative measure of the success of simulations-Validated
model.
Table 2.2 Coefficients for calibration of water temperature model.
Table 2.3 Coefficients for uncertainty analysis.
Table 2.4 Quantitative measure of the success of the simulations-Regional
model.
Table 3.1 Correlation coefficients of daily meteorological variables for
Minneapolis/St. Paul, 1955-1979.
-\
Table 4.1 Lake data.
Table 4.2 Mean monthly meteorological data.
Table 4.3 Differences (°C) in simulated mean daily epilimnetic and
hypolimnetic temperatures for different starting dates of the
model (April 1 reference).
Table 4.4 Monthly averages of daily heat balance components (1000 kcal
m-2 dayi).
Table 4.5 Cumulative heat balance components (1000 kcal nr2).
Table 4.6 Net cumulative heat input (content) per meter of average
depth (1000 kcal m-i).
Table 5.1 Weather parameters changes projected by the 2xCC>2 climate
model output for Minneapolis/St. Paul.
Table 5.2 Lake classification.
Table 5.3 Morphometric regression coefficients in the area vs. depth
relationship.
Table 5.4 Maximum temperatures of southern Minnesota lakes.
Table 5.5 Seasonal stratification characteristics of southern Minnesota
lakes.
IX
-------
Table A.I Regression coefficients for K2 = o(N2)7.
Table A.2 Errors of the eddy diffusivity estimation.
-------
Acknowledgements
We are grateful to the U.S. Environmental Protection Agency, Office of
Program Planning and Evaluation, Washington, D. C. and Environmental
Research 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.
-------
1. Introduction and Literature Review
1.1 Introduction
The concentrations of some gases (CC>2, H20, N20, CE^) have been
increasing in the atmosphere (Bolin and Doos, 1986; NRG, 1982; 1983;
Houghton et al., 1990). These commonly called "greenhouse gases" are
absorbing and reradiating energy at both long and short wavelengths. As a
consequence, greenhouse gases are able to affect global climate possibly
resulting in global mean warming of the earth's terrestrial and aquatic surface
and the lower atmosphere (Bolin and Doos, 1986; NRG, 1982; 1983; Wanner
and Siegenthaler, 1988; Waggoner, 1990).
Special attention has been paid to the increase of carbon dioxide
because it is estimated that about half of the temperature change is due to
the increase of atmospheric CC-2 alone. Mathematical models of global
climate change lead to the conclusion that the increase in mean global
equilibrium surface temperature for a doubling of CO2 is most likely to be in
the range of 1.5 to 5.5° C (Bolin and Doos, 1986, Waggoner, 1990). One of
the uncertainties is due to the transfer of increased heat into the oceans
(NRG, 1982; 1983, Waggoner, 1990). Surely due to their high heat capacity,
oceans will act as a sink for heat and delay the warming.
The question which we want to address in this report is how freshwater
lake temperatures respond to atmospheric conditions. Changes in lake water
temperatures and temperature stratification dynamics may have a profound
effect on lake ecosystems (Meisner et al., 1987; Coutant, 1990; Magnuson et
al. 1990; Chang et al., 1992). Dissolved oxygen, nutrient cycling, biological
productivity, and fisheries may be severely affected through temperature
changes.
Considerable effort has gone into global climatological modeling with the
objective to specify future climatic conditions in a world with high greenhouse
gases., Some models use statistical analysis of past climatological data in
order to provide analogies for future climatological changes. Unfortunately,
all causes of past climate changes are not fully understood (Bolin and Doos,
1986; Waggoner, 1990), and predictions of future climates are difficult,
especially on a regional basis. Nevertheless simulated climate conditions are
and will be used in numerous effect 'studies. Another approach to finding
both climatic trends and their effects is to examine long-term records. In
few lakes, e.g. in the experimental lake area (ELA) in Ontario, Canada,
weekly or biweekly vertical profiles of water quality and biological parameters
have been collected over periods of 20 or more years and these records reveal
e.g. rising average surface water temperatures, shorter ice cover periods, etc.
(Schindler et al., 1990). A data record of more than 100 years for Lake
Mendota was analyzed by Robertson (1989) and Magnuson et al. (1990).
-------
To make generalizations to lakes of different geometries and latitudes,
and to extrapolate to possible future climates, numerical simulation models
(McConnick, 1990; Robertson and Ragotzkie, 1990) are useful. Herein the
use of such a model is demonstrated by application to morphometrically
different lakes with sparse data sets. The lakes are located near 45° northern
latitude and 93° western longitude in the northcentral United States.
1.2 Previous Temperature Prediction Model
A one-dimensional lake water quality model, which has been successfully
applied to simulate hydrothermal processes in different lakes and for a variety
of meteorological conditions (Stefan and Ford, 1975; Stefan et al., 1980a; Ford
and Stefan, 1980) was used in this study. The model was previously
expanded to include suspended sediment (Stefan et al., 1982V light
attenuation (Stefan et al., 1983), phytoplankton growth and nutrient dynamics
(Riley and Stefan, 1987). Only the hydrothermal part of the model was
applied in this study.
, /
1.2.1 Model formulation
In the model the lake is described by a system of horizontal layers,
each of which is well mixed. Vertical transport of heat is described by a
diffusion equation in which the vertical diffusion coefficient Kz(z) is
incorporated in a conservation equation of the form:
A fr = I
-------
CO
Ha
atmospheric
radiation
Hsn
He
He
Hbr
incoming reflected evaporation convection back
solar fraction radiation
radiation
t r t
surface absorption
absorption with depth
fraction transmited
to next layer
absorption with depth
fraction transmited
to next layer
Fig. 1.1 Schematic diagram of source and sink terms in the heat budget model.
-------
where Hs is incoming solar radiation (kcal m^day"1), r is the reflection
coefficient computed as a function of the angle of incidence and the
concentration of suspended sediment in the surface layer (Dhamotharan, 1979;
Stefan et al., 1982). /? is the surface absorption factor (Dake and Harleman,
1969). The attenuation of solar radiation with depth follows Beer's law:
H(i) = H(i-l) exp(-M Az) (1.4)
gn
where Hsn(i-l) is solar radiation at the top of a horizontal layer of water
__ n ^*
(kcal m day ), Hsn(i) is solar radiation at the bottom of a layer, A is
thickness of a layer (m), fj, is the extinction coefficient (m~ )
fj. = \i + u -SS + /i.Chla (1.5)
r *w rsa ;ch v '
where /^ is the extinction coefficient of lake water (m~ ), ^ss is the specific
extinction coefficient due to suspended sediment (1 m^mg"1); SS is suspended
inorganic sediment concentration (mg I"1); \L& is the extinction coefficient due
to chlorophyll (m2 g"1Chla)(Bannister> 1974), Chla is chlorophyll-a
concentration (g m'3).
H = a e T4 (1.6)
a a a v '
where a is Stefan-Boltzmann constant, Ta is absolute temperature (°K), ea is
atmospheric emissivity (Idso and Jackson, 1969). Back radiation Hbr follows
the same formulation (6), but the emissivity is fixed at 0.975, and
atmospheric temperature is replaced by water surface temperature Ts .
Aerodynamic bulk formulae were used to calculate surface wind shear r,
latent heat flux H , and the sensible heat flux H across the water surface
c' c
(Keijman, 1974; Ford and Stefan, 1980; Strub and Powell, 1987; Sadhyram et
al., 1988):
Uf (1-7)
H = p c 8 ' u' = p c C u, 6. = p c f(U )(T - T ) (1.8)
c ra p 'a p s * * 'a p v a'v s a' v '
H = p L q ' u' = p-L C^u. = p L f(U )(q - q ) (1.9)
e ra v n ra v C1* * a v x a/vT a' v '
where T is the surface wind stress, pa is the density of the air, u' and u;'
are. turbulent fluctuations of velocity in horizontal and vertical direction; the
overbar represents a time average; u, is a velocity scale, Ua is the wind
speed above the water surface, Cd is the momentum or drag coefficient (Wu,
1969),
0' is turbulent fluctuation in temperature, 6^ is a temperature scale, Cs and
Cf are heat transfer and vapor transfer coefficients, respectively, and together
-------
with u. are expressed as a function of wind speed, f(Ua), (Ford, 1976), Ts is
water surface temperature, Ta air temperature above the water surface, Lv is
'.atent heat of vaporization, q' is the specific humidity fluctuation, q^ is the
specific humidity scale, qa is the specific humidity above the water surface, qs
is the specific humidity at saturation pressure at the water surface
xemperature.
Turbulent kinetic energy supplied by wind shear and available for
possible entrainment at the interface was estimated (Ford and Stefan, 1980)
bv
TKE = Wstr
dA
(1.10)
where As is lake surface area (m2), Ut is shear velocity in the water (m
day-i), and Wstr is the wind sheltering coefficient.
The model distributes the surface heat input in the water column using
turbulent diffusion (Eq. 1) in response to wind and natural convection (Ford
and Stefan, 1980). The numerical model is applied in daily timesteps using
mean daily values for the meteorological variables. Initial conditions, model
set-up parameters, and daily meteorological variables average air temperature
(Ta), dew point temperature (Td), precipitation (P), wind speed (Ua), and
solar radiation (Hs) have to be provided to use the model.
1.2.2
Model coefficients
Model calibration coefficients needed for simulations of lake water
temperatures are given in Table 1.1. These coefficients are kept at their
initially specified value throughout the entire period of the simulation.
Table 1.1 List of calibration coefficients with ranges used in previous
simulations.
Coefficient
Symbol Units
Range of values
Previous Literature
Simulations Values
Radiation extinction
by water
Radiation extinction
by chlorophyll
Wind sheltering
Wind function
coefficient
Maximum hypolimnetic
eddy diffusivity Kzn,ax
(mi)
(m2 g-i Chla)
0.4-0.65
0.02-2.0
W
str
8.65-16.0 0.2-31.5
0.1-0.9 0.1-1.0
20.0-30.0 20.0-30.0
(m2 day1)
0.1-2.0
0.086-8.64
-------
Radiation extinction coefficients by water (^w) and chlorophyll
specify the rate of attenuation of short-wave radiation energy as it penetrates
through the water column. Both coefficients vary as a function of the
wavelength. Usually these coefficients are reported by a single mean spectral
value for a given lake. Smith and Baker (1981) measured a range of
0.02-2.0 (m*1) for /iw as a function of the wavelength. Values of fa in the
range of 0.68 ± 0.35 (m-i) have been reported by Megard et al. (1979).
Chlorophyll extinction coefficient is species dependent. Values in the range of
the 0.2-31.4 (m2 g'1 Chla) with a mean spectral value of 16.0 have been
reported by Bannister (1979; 1974) for IL<± while Megard et al. (1979)
reported values of 22 ± 5 (m2 g'1 Chla) for the photosynthetically active
radiation.
The wind sheltering coefficient (Wstr) determines the fraction of
turbulent kinetic energy from the wind applied at the lake surface and
available for mixing. The coefficient can range from 0.1 to 1.0 depending on
the size of the lake and the terrain surrounding the lake. The coefficient
defines the "active" portion of the lake surface area on which wind shear
stress contributes to the turbulent kinetic energy.
The wind function coefficient is defined for the neutral boundary layer
above the lake surface. This condition occurs for the case of negligible
atmospheric stratification. The wind speed function used is linear with the
wind speed
f(Ua) = cUa (l.H)
where c is defined as a wind function coefficient. The atmosphere above
natural water bodies is often nearly neutrally stable. A significant amount of
experimental and theoretical research has been done in regard to wind
function coefficient estimation (e.g. Dake, 1972, Ford, 1976, Stefan et al.,
1980b, Adams et al., 1990). Different ranges of coefficients were reported
depending on measurement location of the windspeed Ua relative to the lake
surface. Herein the wind function coefficient is taken to be in the range
20-30 if wind speed is in mi h'1, vapor pressure in mbar, and heat flux in
kcal m^day1.
Maximum hypolimnetic eddy diffusivity is the threshold value for the
turbulent diffusion under negligible stratification. In modeling this condition
is assumed to be satisfied by small stability frequency e.g. N2 = 7.5*10'5
sec"2. Maximum hypolimnetic eddy diffusivity ranges from 8.64 m2 day'1 for
large lakes (Lewis, 1983) to 0.086 m2 day1 for small lakes (Appendix A).
1.3 Overview of Study
The goal of this study is to develop an understanding of how freshwater
lake temperatures respond to atmospheric conditions.
Chapter 2 presents the regional lake water temperature model
development and validation. The lake water temperature model, which was
originally developed for particular lakes and particular years has been
generalized to a wide range of lake classes and meteorological conditions.
-------
Chapter 3 presents a first order analysis of uncertainty propagation in
lake temperature models. The source of the uncertainty is variable
meteorological forcing which enters the lake temperature equations through
the source term and boundary conditions. The analysis presents a useful
alternative for the study of long-term averages and variability of temperature
structures in lakes due to variable meteorological forcing.
Chapter 4 presents a lake water temperature model application to a
particularly warm year (1988) and a normal year (1971) for comparison. The
comparison is made for morphometrically different lakes located in the north
central US. The analysis was a first step in quantifying potential thermal
changes in inland lakes due to climate change.
Chapter 5 presents a lake water temperature model application to a
representative range of lakes in Minnesota for past climate and a future
climate scenario associated with doubling of atmospheric CO2- Emphasis was
on long term behavior and a wide range of lake morphometries and trophic
levels. The base weather period (or reference) was for the years from
1955-1979. For future climate scenario the daily weather parameters were
perturbed by the 2XC02 GISS climate model output. The simulation results
showed how water temperatures in different freshwater lakes responded to
changed atmospheric conditions in a region.
Chapter 6 summarizes the results of the study.
-------
2. Regional Lake Water Temperature Simulation
Model Development
A lake water temperature model was developed to simulate the seasonal
(spring to fall) temperature stratification structure over a wide range of lake
morphometries, trophic and meteorological conditions. Model coefficients
related to hypolimnetic eddy diffusivity, light attenuation, wind sheltering,
and convective heat transfer, were generalized using theoretical and empirical
model extensions. The new relationships differentiate lakes on a regional
rather than individual basis. First order uncertainty analysis showed
moderate sensitivity of simulated lake water temperatures to model
coefficients. The proposed regional numerical model which can be used
without calibration has an average 1.1* C root mean square error, and 93% of
measured lake water temperatures variability is explained by the numerical
simulations, over wide ranges of lake morphometries, trophic levels, and
meteorological conditions.
2.1 Introduction
Changes in meteorological variables in the future "greenhouse"
atmospheric conditions are usually specified through the global climate change
models output on a regional rather than a local scale. Usually water
temperature data are only available for a few lakes, not necessarily for
"typical" lakes in order to calibrate lake water temperature model and to
validate predictions. Some coefficients such as eddy diffusion coefficients or
turbulence closure coefficients used in lake water temperature models are not
universal due to their dependence on stratification, and the longer than
subdaily time step of the simulations (Aldama et al.,, 1989).
The purpose of this chapter is to describe how a lake temperature
model, which was described in Chapter 1, and which was initially developed
for particular lakes and particular meteorological years, could be generalized
to a wide range of lake classes and meteorological conditions. To do this,
new functional relationships had to be introduced for the calibration
coefficients which differentiate lakes on a regional rather than an individual
basis. The generalized model can than be applied to lakes for which no
measurements exist. Fortunately -it can be demonstrated that the regional
model makes prediction almost with the same order of accuracy as the
validated previous calibrated to particular lakes. Therefore regional and long
term lake temperature structure modeling rather than short time behavior of
particular lakes can be accomplished with same confidence.
-------
2.2 Model Generalization
In order to apply the lake water temperature numerical model to lakes
for which there are no measurements, the model has to be generalized. This
was accomplished by introducing functional relationships for the model
coefficients which are valid for lakes on a regional rather than individual
basis.
2.2.1 Hypolimnetic diffusivity closure • . .
Although the hypolimnion is isolated from the surface (epilimnetic) layer
by the thermocline and its associated density gradient, strong and sporadical
local mixing events have been observed in the hypolimnion (Jassby and
Powell, 1975; Imberger, 1985; Imberger and Patterson, 1989). Heat flux
between water and lake sediments was found to be important in eddy
diffusivity estimation for inland shallow (10 m maximum depth) lakes,
representative for north central United States (Appendix A). Hypolimnetic
eddy diffusivity dependence on stratification strength measured by buoyancy
frequency has been pointed out consistently (Quay et al., 1980; Gargett, 1984;
Gargett and Holloway 1984; Colman and Amstrong, 1987; Appendix A).
Stability frequency is related to hypolimnetic eddy diffusivity by :
K2 = a (N2)7 (2.1)
where stability frequency N2=-(dp/cte)(g//?), in which p is density of water,
and g is acceleration of gravity, 7 is determined by the mode of turbulence
production (narrow or broad band internal waves, local shear etc.), and a is
determined by the general level of turbulence. For most inland lakes,
coefficient 7 ranges from 0.4 to 0.6 (Jassby and Powell, 1975; Quay et al.,
1980; Gerhard et al., 1990; Ellis and Stefan, 1991, Appendix A).
Hypolimnetic eddy diffusivity estimations in five northern Minnesota
lakes follow Eq. 2.1 as shown in Fig. 2.1. Lakes were selected from the
regional prospective i.e. with different surface areas and maximum depths.
Dimensionless analysis (Ward, 1977) suggests that lake surface area, can
provide the horizontal scale for the vertical eddy diffusivity estimation. The
vertical scale (lake depth) is implicitly built into the stability frequency. The
a coefficient in Eq. 2.1 can be interpreted as a measure of turbulence level
and is plotted as a function of lake surface area in Fig. 2.2. A general
relationship applicable to lakes on a regional scale was therefore summarized
as:
Kz = 8.17 x 10-4,(Area)0'56 (N2)-°'« (2.2)
where Area is lake surface area (km2), N2 is in sec'2, and Kz is in cm2 sec'1.
Maximum vertical hypolimnetic eddy diffusivity Kzmax was also
correlated with lake surface area because turbulent mixing in non-stratified
lakes depends strongly on kinetic wind energy supplied, which in turn depends
on lake surface area. Maximum hypolimnetic eddy diffusivity versus lake
surface area for eight different lakes is plotted in Fig. 2.3. Data are from
Jassby and Powell (1975), Ward (1977), Lewis (1983), and from this study.
-------
10°-
Ann Lake (0.47 km2)
o
(U
E
o
o
d>
c/i
(N
E
o
10-
10-J^
10-M
10
• %• : •
• *
.
Big Carnelion Lake ( 1 .89 km")
Square Lake (0.85 km2)
Big Marine Lake (74 krna)
-»
10-2
10
-3
10'2
N2(sec-2) N2(sec-2)
Fig. 2.1 Hypolimnetic eddy diffusivity dependence on lake surface area.
-------
= 8.17X10-4(Area)°-56
Ryan Lake (0.06 km2)
Ann Lake (0.47 km2)
Square Lake (0.85 km2)
Big Cornelian Lake ( 1 .89 km")
Big Marine Lake (7.4 km)
Lake Mcllwanie (26.3 km)
10-*-
Area (km2)
Fig. 2.2 Hypolimnetic eddy diffusivity forcing parameter (a) dependence on
lake surface area.
-------
U
(D
(f)
CM
U
N
10°-=
IU '
-2_
10
-J_
o
D
i i i iii| 1—i—i i 11ii| 1—i—i i 111ii r
Ryan Lake (0.06 km2)
Castle Lake (0.2 km2)
Ann Lake (0.47 km2)
Square Lake (0.85 km2)
Big Cornelian Lake (1.89km2)
Big Marine Lake (7.4 km2)
Lake Mcilwanie (26.3 km2)
Lake Valencia (350 km2)
Ryan Lake (winter)
K7mox = 0.048(Area)°-56
10-2
10
-'
Area
101
102
Fig. 2.3 Maximum hy poll nine lie eddy diffusivity (N2=7.5*10'5 sec"2)
dependence on lake surface area.
-------
Attenuation coefficient
The specific radiation attenuation coefficients for water and chlorophyll
• •-TO replaced by the total attenuation coefficient. This was done following
•. '..(• parsimonious principle i.e. the fewer coefficients in the model, the less
.-.ncLTtain the model estimate. In addition uncertainty analysis showed that
hlorophyll-a made a minor contribution to lake water temperature
uncertainty. A relationship between total attenuation coefficient \i (m'1) and
y-cchi depth zsd (m) was obtained from measurements in 50 lakes in
Minnesota (Osgood, 1990) and is plotted in Fig. 2.4.
It = 1.84 (Zsd)-i (2.3)
I'he form of this relationship has been found to be valid in inland waters in
St-neral (Idso and Gilbert, 1974) and in the ocean (Poole and Atkins, 1929).
L' 2.3 Wind sheltering coefficient
The wind sheltering coefficient is a function of lake surface area (fetch).
The turbulent kinetic energy computation (Eq. 1.10) uses a wind speed and
direction taken, from off-site weather station at 10 m elevation and adjusts
that wind speed for fetch over the lake in the direction of the wind. As
wind speed typically increases with fetch, the calculated downstream wind
speed is an estimate of the maximum wind speed on the lake surface.
Typically fetch on a lake is reduced by wind sheltering the upwind side of
the lake where the wind makes a transition from a landbound turbulent
velocity profile to the open water. This was explained by Ford and Stefan
(1980). The reduction in fetch or surface area sheltered from direct wind
access by trees or buildings along the shoreline will be more significant for
small lake than a large one because a) a relatively larger portion of the total
lake surface area will be wind sheltered b) the downwind maximum wind
speed does not grow linearly with fetch and will on a large lake be near the
real wind speed over a large portion of the lake surface area, and c) wind
gusts will be less effective over a small lake surface than a large one because
of spatial averaging. Also lake morphometry, i.e. distribution of area with
depth will be a factor in the translation of wind energy into mixing. A
maximum wind speed at the downwind end of a large lake will also be more
representative for a large lake than a small one, especially if the lake
morphometry is taken into consideration.
For all these reasons a very strong dependence of the wind sheltering
coefficient (Wstr) on lake surface area can be expected. A functional
relationship was obtained by plotting the wind sheltering coefficient obtained
by calibration in several previous numerical model simulations (Fig. 2.5).
Biweekly temperature profile measurements in ten lakes and throughout the
summer season were used to optimize the wind sheltering coefficients plotted
in Fig. 2.5. ' The empirical relationship is
Wstr = 1-0 - exp(-0.3*Area) (2.4)
where Area is the lake surface area in km2.
13
-------
84*(secchi depth)-1
.,,,,,,.,
0.3 0.5 0.8 1.0 1.2 1.5
(Secchi Depth m )"'
Fig. 2.4 Relationsliip between total attenuation coefficient and Sechi disk
depth.
-------
c
OL)
O
0>
O
O
O>
C
102-
10'-
10°-:
CD
±L 10--:
0>
00
? 10-H
io-3-
10-'
Wstr= 1,0-EXP(-0.3*Area)
10-2
10-
rnrj
10°
Surface Area (km2)
10'
E Williams (0.35)
X Riley (1.2)
^ Waconia (10.0)
O Square (0.8)
A Thrush (0.07)
A Greenwood (7.70)
ffl Fish (1.16)
D Elmo (1.23)
• Cedar (3.30)
O Calhoun (1.71)
102
Fig. 2.5 Wind sheltering coefficient dependence on lake surface area.
-------
The result in Fig. 2.5 seems to indicate that the modeling of wind
mixing in lakes, especially small ones, depends more on a correct amount of
energy supplied than on a energy dissipation. This is a new insight which
appears to result from this study.
2.2.4 Wind function coefficient
Wind function coefficients (c) defined in Eq. 1.11 enters into the heat
transfer relationships (1.8) and (1.9), and depends also on lake surface area
(fetch) as was found by Harbeck (1962), Sweers (1976), and summarized by
Adams et al. (1990). Harbeck (1962) analyzed data from several lakes of
different sizes and pointed out that evaporation rates in small and large lakes
-might be the same. The fetch dependence is introduced mainly due to the
wind speed increase over the water. As air flows from land to a smoother
water surface, at a constant height above the water (e.g. 10 m), its velocity
increases with fetch. In this numerical model off-lake wind speeds measured
at permanent weather stations are are used, but they are adjusted for lake
fetch (Ford and Stefan, 1980). Nevertheless, some residual wind function
coefficient dependence on lake fetch is shown in Table 2.1. A functional
relationship was obtained by plotting the wind function coefficient from
several previous numerical model simulations against lake surface area (Fig.
2.6). The estimated relationship is
c = 24+ln(Area) (2.5)
where Area is again in km2. This relationship shows only a week dependence
of c on lake surface area, and can be viewed as a minor adjustment. The
need for this adjustment can be explained by examining the wind boundary
layer development over the surface of small and large lakes (see Fig. 2.7).
Wind speed increases with fetch (distance from the leeward shore) but
non-linearly. In our model wind speed is taken from an off-lake weather
station and a maximum wind speed at the downwind end of the lake as
shown in Fig. 2.7 (top) is computed for the use in the heat transfer
equations (1.8) and (1.9). This calculated wind speed is an overestimate of
the area! average wind speed over the lake surface. Because of the
non-linearity of wind speed with distance the overestimate is more severe for
small lakes than for large lakes. Therefore the wind function coefficient has
to be smaller for smaller lakes in order to compensate for the wind velocity
overestimate. If on the other hand, wind speeds are measured on the lake
(middle of the lake) as shown in Fig. 2.7 (bottom) the situation is reversed.
In that case the wind measurements on a small lake are severely
underestimated relative to the area! average than on a big lake. For this
reason the wind function coefficient has to decrease with fetch (surface area)
to compensate for this non-representativeness of the wind speeds measured in
midlake. This decreasing trend of wind function coefficient with lake surface
area was found and reported by Harbeck (1962), Sweers (1976) and
summarized by Adams et al. (1990). In addition Adams called upon
increases in relative humidity with fetch over cooling ponds to justify the
decrease in wind function coefficient with fetch.
It is concluded from all of the above that the wind function coefficient
can increase or decrease with lake surface area depending on the location
where wind speed is measured.
16
-------
40-
35-
G>
'u
o 25^
20-
15-
10-2
Wcoef = 24 + In(Area)
.O O «
X
IO-1 10° 10'
Surface Area (km2)
\ ~—
X Williams (0.35)
X Riley (1.20)
4 Woconia (10.0)
O Square (0.8)
A Thrush (0.07)
A Greenwood (7.70)
ffl Fish (1.16)
D Elmo (1.23)
• Cedar (3.3)
O Calhoun (1.71)
Fig. 2.6 Wind function coefficient dependence on lake surface area.
-------
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
diffusivity closure (Eq. 2.2). The number of calibration coefficients (Table
1.1) was thereby reduced from five to four. The modified numerical model
than had to be validated with water temperature measurements in several
selected lakes over a period of several years. Representative lakes in
Minnesota were selected through an analysis of the state's extensive data
bases. Differences between waterbodies in adjacent ecoregions were found too
small to justify further subdivisions on this basis. The state was divided into
a northern part, roughly coinciding with three ecoregions, and a southern
part, roughly coinciding with three other ecoregions (Fig. 2.8) which also
extended into Wisconsin, Iowa, and South Dakota. "Representative" lake
meant either having values of lake surface area, maximum depth, and secchi
depth near the median as identified in a state report by the Minnesota
Pollution Agency (Heiskary et al, 1988) or being near the far ends of the
respective frequency distributions for ecoregions. Selected representative lakes
with their position on the cumulative frequency distribution curves for
northern and southern Minnesota are given in Fig. 2.9. Lakes covered the
entire range of maximum depths (shallow-medium-deep), surface area
(smaU-medium-large), and trophic status (eutrophic-mesotrophic-oligotrophic).
Geographical distribution of these lakes in Minnesota is given in Fig. 2.8.
To validate 'the model numerical simulations were started with
isothermal conditions (4 *C) on March 1 and continued in daily timesteps
until November 30. Ice goes out of Minnesota lakes sometime between the
end of March and beginning of May. Dates of spring overturn vary with
latitude and year. To allow for these variable conditions, a 4"C isothermal
condition was maintained in the lake water temperature simulations until
simulated water temperatures began to rise above 4*C. This method"
permitted the model to find its own date of spring overturn (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
p P
V > Tsi V i
(2-6)
19
-------
NorVhern Minnesota Wetlands
^""•s^^^^ i
N 0 RTH "\_XNa- GREENWOOD LAKE
*&
A
THRUSH LAKE
rrgSv,
-------
!00
10-'
100
10-' 10° 10'
AREA (km2)
TOO-
MAXIMUM DEPTH (m)
80-
£ 60-
i 40-
20-
0-r
O
o
D
O
A
D
O
A
A
A
•
m
D
O
468
SECCHI DEPTH (m)
10
Waconia (10.0)
Greenwood (7.2)
Cedar (3.3)
Williams (0.35)
Calhoun (1.7)
Elmo (1.2)
Fish (1.1)
Square (0.8)
Thrush (0.07)
Elmo (42.0)
Greenwood (34.0)
Calhoun (24.0)
Square (20.0)
Waconia (11.0)
Williams (10.0)
Thrush (14.0).
Fish (8.0)
Cedar (4.7)
Square (7.0)
Greenwood (6.5)
Elmo (3.0)
Calhoun (2.5)
Williams (2.1)
Thrush (7.3)
Waconia (1.9)
Cedar (1.2)
Fish (1.0)
Fig. 2.9
Cumulative distributions (%) of lake parameters in Minnesota.
Lakes selected for model validation are shown by symbols.
21
-------
(b) temperature root mean square errors
P P
0.5 r 2rfVil
Ejf=
I*
0.5
(2.7)
and (c) r2 i.e. portion of the temperature measurements explained by the
simulations (Riley and Stefan, 1987). In the above equations subscripts i, s,
and m refer to the counting index, simulated, and measured temperature
respectively. Vj is the water volume of a layer in the stratified lake. The
above parameters are estimated by summing over lake depths. Overall
seasonal average parameters are reported in Table 2.1. Examples of
simulated and measured vertical lake water temperature profiles are given in
Figs. 2.10, 2.11, 2.12, 2.13, and 2.14. The model simulates onset of
stratification, mixed layer depth and water temperatures well.
Table 2.1 Quantitative measure of the success of the
simulations-Validated model.
Lake
Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Year
1971
1984
1988
1987
1985
1985
1986
1986
1984
Tm
(oC)
14.37
20.64
13.94
24.40
14.37
20.14
11.80
11.97
17.26
Ts
(oC)
14.52
20.86
14.09
24.13
14.52
20.12
11.97
11.91
16.37
(oC)
0.86
0.94
1.77
0.80
0.86
0.78
0.89
0.90
1.08
E2
(oC)
0.79
0.99
1.80
0.82
0.79
0.73
0.79
0.91
1.07
.*
0.97
0.93
0.92
0.90
0.97
0.92
0.93
0.96
0.97
Number of
field data.
136
20
214
32
136
43
46
114
110
Average
'T--
Ts-
E,-
E2-
16.54 16.49 0.97
0.96
0.94
95
measured volume weighted average temperature
simulated volume weighted average temperature
temperature root mean square error
volume weighted temperature root mean square error
portion of the measured water temperature variability explained by
the simulations
22
-------
simulated
measured
s/20/7i
e/oa/7i
6/39/71
7/15/71
&
o .J
1
e/02/71
S/2V1
I
LJ
o
16-
tO/ll/71
TO/25/71
10 IS 20 25
TEMPERATURE (*C)
10 IS 20 25
TEMPERATURE (-C)
Fig. 2.10 Lake Calhoun water temperature profiles.
23
-------
I 10
a
15-
• measured
— simulated
V15/85
5/O9/SS
10-
a
Ul
o
15-
20-
6/11/85
7/12/85
.5-
a.
UJ
o
15-
20-
S/07/SS
9/CS/85
S 10 IS 20 25
TEMPERATURE (-C)
5 10 IS 20 25 20
TEMPERATURE (*C)
Fig. 2.11 Square Lake water temperature profiles.
24
-------
to
Cn
X
CL
LJ
O
I
CL
UJ
Q
u-
1-
2-
3-
4-
5-
6-
7-
8-
9-
10-
n_
1-
2-
3-
4-
5-
6-
7-
fl-
9-
10-
11-
• measured
i i _i
simulotea
5/17/85
'
1
i
•
•
T
.
•
• - -
-
.
-
-
-
8/20/85
1 * * 1 " *
•
•
• -
•
•
,
•
•
•
\ *
•
7/22/85 •
•
•
•
•
•
l
l
l
i
i
(
/
9/18/85
/*
f .
1
10 15 20 25
TEMPERATURE (°C)
10 15 20 25 , 30
TFMPERAfURE (°C)
Fig. 2.12 Waconia Lake water temperature profiles.
-------
1 6
Q. 8
UJ
a
10
12
1*
• measured
simulated
6/17/86
/
10-
12
7/15/86
O- 6-
I
CL
t. i
o
10-
12-
9/'6/56
10/'«/B6
10 15 20
TIMPERATURL CC)
5 10 15 2:
TEMPERATURE cc;
Fig. 2.13 Thrush Lake water temperature profiles.
26
-------
e *-
8-
10
S/17/84
2-
I «:
•*~r
x
8-
10
6/lt/St
7/03/8<
J-
8-
8/OS/S*
8/26/84
£ *•
x
CL
UJ
O
S/IJ/8*
11/02/8*
10 15 20
TEMPERATURE (*C)
25 0 5
10 IS 20
TEMPERATURE (*C)
25
Fig. 2.14 Williams Lake water temperature profiles.
27
-------
Volume weighted and unweighted root mean square error was less than
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 (Ej and £3) indicate the vertical position
of the maximum simulation error. If E2 is greater than Ej, than the
difference between measured and simulated lake water temperatures are
greater in the surface layers because £2 values are volume weighted and EI
values are not.
Table 2.2 Coefficients for calibration of water temperature model
Lake
Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Year
1971/
1984
1988
1987
1985
1985
1986
1986
1984
Max.
depth
Hmax
(m)
24.0
4.70
41.8
8.20
21.0
11.0
34.0
14.0
10.0
Surface
area
As
(km2)
1.71
3.30
1.23
1.16
0.85
10.0
7.70
0.07
0.35
Wind
funct.
coeff.
c
H
24
24
26
26
24
27
29
24
22
Wind
shelt.
coeff.
Wstr
(-)
0.40
0.60
0.50
0.50
0.10
0.90
0.80
0.01
0.20
Attenuation
coefficient
Mw Hch
(m-i)(m2g-iChl-a)
0.65 8.65
0.65 8.65
0.65 8.65
1.00 8.65
0.50 8.65
0.65 8.65
0.65 8.65
0.65 8.65
0.65 8.65
Chl-a
(g m-3)
4-371
6-1302
3-83
18-484
1-47
11-348
1-35
2-l«
3-79
Field data given by:
i
2
3>7>8
4'5
6
9
Shapiro and Pfannkuch, 1973
Osgood, 1984
Osgood 1989
Minnesota Pollution Control Agency, 1988
Wright et al., 1988
Winter, 1980
The average root mean square error for all lakes was 10C, and 94% (r2
= 0.94) of water temperature measurements variability was explained by the
numerical model (Table 2.2).
Model coefficients used in the simulations are given in Table 2.2.
These coefficients give minimum values of root mean square error, and
highest value of r2 between measurements and simulated lake water
temperatures.
In the following sections the modified model with the hypolimnetic eddy
diffusivity closure as described in this section will be referred to as the
validated model.
28
-------
2 4 Numerical Uncertainty of Model After Hypoliranetic Closure
Uncertainty in the lake water temperature simulations was considered in
•.'•rrns of all model coefficients except maximum hypolimnetic eddy diffusivity
15 specified in Table 1.1. To first-order the uncertainty in lake water
u-mperature depends on the uncertainty in the model coefficients, and on the
sensitivity of the lake water temperatures to changes in the coefficients:
PT = E{[T - TJ[T - T]tr) s
E{[T(u) + ^ (u - u) - T][T(u) + ^ (u - u) - T]tr) =
dn dn
— (u - n)}[~ (u - u)]tr}=
dn dn
cm on
OM on
where P is the (m x m) covariance matrix of the simulated lake water
temperatures, m is the total number of discritized lake control volumes, E{.}
is the mathematical expectation, T is the mean lake water temperature, (tr
is the transpose, u is the vector of the n coefficients, Pu is the (n x n
frr
covariance matrix of system coefficients, -— is the (m x n) sensitivity matrix
of partial derivatives of the lake water temperatures with respect to the
coefficients. Sensitivity matrix is estimated using the influence coefficient
method (Willis and Yeh, 1987).
Data for the system coefficients covariance matrix are given in Table
2.3. These values were chosen to be in the range of theoretical and
simulated values (Tables 1.1 and 2.2), and to have coefficients of variation
(standard deviation/mean) equal to 0.3. This value is chosen because
first-order uncertainty analysis could be questionable when the coefficient of
variation of a nonlinear function increases above 0.3.
Table 2.3 Coefficients for uncertainty analysis
Coefficient Lake Calhoun Williams Lake Cedar Lake
mean st. dev. mean st. dev. mean st: dev.
/__ .o
/iw lm lj
/A:h (m2
\vstr
c
1
g-tChla)
0.65
8.65
0.60
24.0
0.20
2.65
0.18
7.20
0.65
8.65
0.20
20.0
0.20
2.65
0.06
6.00
0.65
8.65
0.60
24.0
0.20
2.65
0.18
7.20
29
-------
Three lakes are selected for the lake water temperature uncertainty
estimation. Lake Calhoun is a eutrophic, deep (24 m maximum depth) lake,
Williams Lake is oligotrophic, and has maximum depth close to the median
depth of 3002 lakes in Minnesota (Fig. 2.9), and Cedar Lake is a highly
eutrophic shallow (4.7 m maximum depth) lake.
Standard deviations of smoothed simulated epilimnion and volume
weighted average hypolimnion temperatures are given in Figures 2.15, 2.16,
and 2.17. Although high variability in model coefficients was imposed,
maximum standard deviation in epilimnion temperatures was less than 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 t£> uncertainty in epilimnion lake water temperatures than the
wind function coefficient. Volume weighted hypolimnion temperatures
displayed higher uncertainty than epilimnion temperatures. For Williams
Lake and Lake Calhoun, all three coefficients i.e. water attenuation, wind
sheltering and wind function coefficient significantly contributed to the lake
water temperature uncertainty. Schindler (1988) pointed out that in
oligotrophic lakes dissolved organic • carbon is one of the major light
attenuating factors.
2.5 Accuracy of the Regional Model After Implementation of all Changes
The Number of calibration coefficients was reduced from four to zero.
Functional relationships substituted into the model in Equations 2.2, 2.3, 2.4,
and 2.5. The model output was compared with water temperature
measurements in nine selected representative lakes. Simulations started with
isothermal conditions (4eC) on March 1 and progressed in daily time steps
until November 30. Quantitative measure of the success of the simulations
and differences between the regional model and the validated model of section
2.3 and 2.5 are given in Table 2.4. The average weighted and unweighted
root mean square error was 1.1 °C (16.5 °C average measured lake water
temperature). Ninty three percent of measured lake water temperature
variability was explained by the numerical simulations (r2=0.93). The
regional model has in average 0.15° C higher temperature root mean square
error.
One example of the daily simulated isotherms for the regional and
validated model (section 2.3) is given in Fig. 2.18. Both models simulate
onset of stratification, mixed layer depth and water temperatures in a
virtually identical way.
30
-------
LAKE CALHOUN
- 1
wind sheltering
wind function
water attenuation
chlorophyll—a attenuation
EPILIMNION
2.0-
.-J
0.0-
wind sheltering
wind function
water attenuation
chlorophyll—a attenuation
HYPOL1MNION
MAR APR MAY JUN' JUL AUG SEP OCT NOV
Fig. 2.15 Standard deviations of estimated lake water temperature
• uncertainties.
31
-------
WILLIAMS LAKE
2.5-
2.0-
o:
UJ
Q_
I.O-i
wind sheltering
wind function
water attenuation
chlorophyll-a attenuation
0.0
2.0-
wind sheltering
wind function
water attenuation
chlorophyll —a attenuation
EPILIMNION
HYPOLIMNION
0.0 .-Ti --IT
MAR APR MAY JUN JUL AUG SEP OCT NOV
Fig. 2.16 Standard deviations of estimated lake water temperature
uncertainties.
32
-------
CEDAR LAKE
2.5-
^ 2.0H
1.5-
n
UJ
1.0-
0.5-
0.0-
cm
UJ
CL
2
UJ
t—
2.0-
1.5-
1.0-
Ld
< 0.5H
0.0-
wind sheltering
wind function
water attenuation
chlorophyll—a attenuation
wind sheltering
wind function
water attenuation
chlorophyll—a attenuation
EPILIMNION
| ...... T ,. , I j i . . -j I : I j—ir- , —j i : I ,
HYPOLIMNION
MAR APR MAY JUN - JUL AUG SEP OCT NOV
Fig. 2.17 Standard deviations of estimated lake water temperature
uncertainties.
33
-------
SIMULATION YEAR - 1986
14-r
MAR APR MAY JUN JUL AUG SEP OCT NOV
Fig. 2.18 Simulated temperature (isotherm) structure in Thrush Lake. Top
shows results from validated model and bottom shows results
from regional model.
34
-------
Table 2.4 Quantitative measure of the success of the simulations - Regional model
CO
tn
Lake
Calhoun
Cedar
Elmo
Fish
Square
Waconia
Greenwood
Thrush
Williams
Average
Year
1971
1984
1988
1987
1985
1985
1986
1986
1984
Regional model
Tm
CC)
14.37
20.64
13.94
24.40
14.37
, 20.14
11.80
11.91
17.26
16.54
T8
(•C)
14.44
20.68
14.31
23.90
14.90
20.09
12.61
12.54
16.57
16.67
•E,
«CC)
1.02
1.07
1.83
0.87
1.24
0.68
1.24
0.95
1.26
1.13
Ej
(•C)
0.89
1.15
1.93
0.89
1.03
0.68
0.99
0.97
1.25
1.10
H
0.96
0.91
0.90
0.89
0.95
0.94
0.92
0.95
0.95
0.93
regional model -
AT8
(•C)
-0.08
-0.18
0.22
-0.23
0.38
-0.03
0.64
0.63
0.20
0.17
\
AEt
0.16
0.13
0.06
0.07
0.38
-0.10
0.35
0.05
0.18
, 0.14
Differences
- validated model
AE2
(•c)
0.10
0.16
0.13
0.07
0.24
-0.05
0.20
0.06
0.18
0.12
w
-1
-2
-2
-1
-2
2
-1
-1
-1
-1
-------
2.6 Conclusions
A lake specific water temperature model was generalized for the
application to a wide range of lake classes and meteorological conditions.
Functional relationships which differentiate lakes on a regional rather than on
an individual basis were developed.
Hypolimnetic eddy diffusivity was estimated as a function of lake
surface area, and stability frequency. Equation 2.2 extends Ward's (1977)
analysis to a wider range of lake geometries. Although the proposed
relationship is a significant simplification of the turbulent diffusion processes
taking place in the hypolimnion, it was found to be useful in the seasonal
lake water temperature modeling.
Total attenuation coefficient was estimated as a function of Secchi depth
(Fig. 2.4). Secchi depth is chosen because it can be measured easily and
values are commonly available.
Wind sheltering and wind function coefficient increase with surface area
(fetch) of the lake (Figs. 2.5 and 2.6). The wind function coefficient increase
is very likely an additional adjustment of the wind velocity coming from land
over the lake surface.
Uncertainty analysis revealed moderate sensitivity of simulated lake
water temperatures to the variability of individual model coefficients. This
could be due to the high thermal inertia of the water especially for the
seasonal lake water temperature modeling. Nevertheless epilimnion
temperatures showed 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
••-ctor state-space method. The output uncertainty is defined as the result of
.rviations of the meteorological variables from their mean values. The
ir.ilysis is applied to systems with correlated and uncorrelated meteorological
> .\nables. Surface water temperatures are strongly affected by uncertain
~:<:teorological forcing. Air temperatures and dew point temperature
>.:ctuations have significant effect on lake temperature uncertainty. Ignoring
rrelation in meteorological variables underestimates uncertainties in lake
.'•rnperature estimates. Long-term average water temperature structure in
akos can be estimated by computer model simulation for just one year when
••-suits from a statistical analysis of meteorological variables are used as
nput. The analysis presents a useful alternative for the study of long-term
ivcrages and variability of water temperature structures in lakes due to
variable meteorological forcing.
31 Introduction
It was shown in Chapters 1 and 2 that vertical water temperature
profiles in lakes are related to meteorological variables by heat transport
•yquations which apply basic conservation principles. Atmospheric conditions
.ire the driving force for heat transfer through a lake water surface. Surface
water temperatures of lakes are primarily related to the meteorological forcing
and secondarily to lake morphometries (Ford and Stefan, 1980a).
Observed meteorological variables used in lake water temperature
modeling (Harleman and Hurley, 1976; Ford and Stefan, 1980b) such as solar
radiation, air and dew point temperature, and wind speed are usually single
realizations of the weather process for a particular year. For lake
management purposes and decision analysis we are interested in mean
temperatures as well as expected ranges of water temperature variation due to
the weather variability over a longer period of time. Deterministic lake
water temperature models cannot provide such information from a single
model simulation for a particular year. The stochastic alternative is to
consider meteorological variables as random variables with estimated statistical
properties in terms of first and second moments, and correlation structure.
First and second moment of lake temperatures can then be predicted from a
single mode simulation.
Lake water temperature models are nonlinear dynamic systems.
Approximation techniques, for obtaining the second moment of a dynamic
system output from the moments of its input have been employed in the area
of groundwater hydrology (Dettinger and Wilson, 1981; Mclaughlin, 1985;
Townley and Wilson, 1985; Protopapas and Bras, 1990; McKinney, 1990).
37
-------
Generally, three techniques are available i.e. (1) Monte Carlo, (2) derived
distribution, and (3) perturbation approach techniques. Monte Carlo methods
have been proposed in lake water quality modeling of phytoplankton,
herbivores, nitrate, and available phosphorus (Scavia et al., 1981; US Army
Corps of Engineers, 1986; Canale and Effler, 1989). Although simple,
limitations of this approach have been related to the large number of
simulations. In addition, the prescribed probability distribution for the
coefficients could change in time-varying systems. The derived distribution
approach is not applicable because of the complex relationship between inputs
(meteorological variables] and outputs (lake temperatures). The perturbation
approach utilizes generally two methods: time domain (state-space) methods
of the Taylor series expansion type, and spectral (frequency domain) methods.
As pointed out by Protopapas and Bras (1990), state space methods are
advantageous for the time variable boundary conditions.
In this Chapter we employed the perturbation vector state-space
approach to propagate uncertainty of meteorological input variables into a
lake temperature model. This study follows the work of Protopapas (1988)
who used the state-space approach for uncertainty propagation of
meteorological inputs in a soil/plant model.
The question we want to address in this Chapter is how to predict the
lake temperature uncertainty due to the variability of the meteorological
forcing in time. This analysis quantifies contribution of each meteorological
variable to temperature uncertainty separately. Secondly, we will demonstrate
that a long-term average thermal v structure in a lake can be obtained without
running a water temperature model for several years of meteorological data.
3.2 Numerical Model
In this study a one-dimensional lake water temperature model, which
has been previously described in Chapter 1, was used. Lake temperature is
represented by a nonlinear partial differential equation (1.1). Nonlinearity
comes through the boundary conditions and hypolimnetic eddy diffusivity.
Analytical solution of this equation is possible only under certain
approximations (Dake and Harleman, 1969). Equation (1.1) is discretized
numerically (Appendix B) using an implicit control volume method. This
leads to a system of equations in the form
Ac(K(k),G) T(k+l) = T(k) + H(k) (3.1)
where Ac is a system (mxm) tridiagonal matrix, m is the number of
discretized control volumes, T(k-fl) is a (mxl) vector with lake temperatures
at time step k+1, Kfk) is a (mxl) vector with lake eddy diffusivity
parameters; note that K(k) = f(T(k), Ws(k)), Ws is a wind speed, H(k) is a
(mxl) vector function with source term parameters, and G is a (mxl) vector
with lake geometry parameters. Boundary conditions are treated through the
source term. The control volume approach, satisfies the heat balance in the
computational domain regardless of the number of discretized control volumes
(Patankar, 1988).
38
-------
The numerical model is applied in daily time steps using mean daily
-•iJut's for the meteorological variables. The required meteorological variables
i.-r solar radiation, air temperature, dew point temperature, wind speed and
•.::'Xtion. Initial conditions, model setup parameters, have to be provided to
.u- the model.
Taylor series expansion is commonly used for the linearization of
runctional relations around nominal values. The function and its first
Derivative must be defined at the nominal point. Expanding equation (3.1)
- a Taylor's series around the nominal value and keeping first order terms,
ovts a linear perturbation temperature equation.
Ac(k) T'(k+l) = B(k) T'(k) + F(k) C'(k) (3.2)
Nominal (mean) values and first order derivatives evaluated at these values
x-'i.- denoted by circumflex. Perturbations of the water temperatures T'(k),
-.r.d meteorological variables C'(k) are denoted by primes are defined as:
/ .
T'(k) = T (k) - T (k), C'(k) = C(k) - C(k) (3.3)
The tridiagonal matrix Ac(k) is equivalent to the matrix Ac(k) of the
A A
icterministic temperature model. Matrices B(k) and F(k) require evaluation
:' the first order derivatives of all terms in equation (3.1) which contains
..\ke temperature and meteorological variables at time step k respectively.
A A A
Details about entries in matrices Ac(k), B(k), and F(k) are given in
Appendix C. Terms with the same perturbation parameter are collected
before entries into the matrices. Equation (3.2) can be rewritten as
T'(k+l) = 0(k) T'(k) + #k) C'(k) (3.4)
where ^(k) is transition matrix 0(k)=Ac~1(k) B(k), and ^k)=Ac~1(k)F(k).
The first term in equation (3.4) describes unforced dynamics of the
system while the second term describes the variation of the meteorological
forcing function. A schematic illustration of the lake temperature
perturbation system is given in Fig. 3.1. Air temperature (Ta), dew point
temperature (Td), solar radiation (Hs), and wind speed (Ws) are forcing
meteorological functions. A transition matrix ^(k) connects the state of the
system between time steps.
3.3 First and Second Moment Development
Taking the expected value of equation (3.1) yields the first order
estimate of the mean of water temperature
E [T(k+l)] = T(k+l) = Ac-i(k) (T(k) + H(k)) (3.5)
Notice that the first order estimate of the mean water temperature is exactly
the value obtained through the deterministic approach.
39
-------
Ws(k)
Fig. 3.1 Schematic illustration of the lake temperature perturbation
system.
40
-------
A recursive, solution of equation (34) is
T'(k) = #k,0)T'(0) + #k,n+l) #n) C'(n); T'(0)=0
n=0
T'(k) = S #k,n+l) V
-------
ST/T,(k) = ^(k,n+l)^n)Muc(k)k)V
-------
Lake temperature covariance propagation is calculated in the following
(l) Set isothermal (4°C ) initial steady state conditions for lake water
*™.;.<:ratures in spring, initialize covariance matrix of meteorological
,«-:v.ubations; (2) read meteorological variables, mean and perturbation values,
•: the next time step; (31 using mean values for meteorological variables,
mpute first moment of lake temperature profile for the next time step
• A
' r-i«ion (3.5), store matrix Ac(k); (4) compute matrices B(k), F(k), i.e. first
-' in derivatives with respect to the perturbed meteorological variables and
T.:mated lake temperatures (Appendix C); (5) calculate matrix N(k) for the
-dated case (Equation 3.12); (6) compute transition matrix 0(k), and
>«ii. (7) calculate additional term P(k) (Equation 3.13) for the correlated
ii-c. (8) propagate and store temperature covariance matrix £, , for the
vtt time step (correlated case Equation 3.14, uncorrelated Equation 3.10);
-• store transition matrix <£(k), and ^(k), if correlated case, store in addition
•« k), and S(k); (10) go to step 2 if last day of simulation is not reached.
: * Lake Calhoun — Application
The test lake, Lake Calhoun, is a temperature zone dimictic lake. The
lie is eutrophic with maximum depth of 24 m, and surface area of 1.7 km2.
Vu-teorological data used are from the Minneapolis St. Paul International
i.-rport, located 10 km from the lake.
Daily meteorological data time series (1955-1979), averaged over 25
• cars, are given in Fig. 3.2. Long term means of solar radiation, dew point
•.•-mperature, and air temperature display typical seasonal cycles. Means are
increasing till the end of summer and decreasing towards fall. Perturbations
inandard deviations) for meteorological forcing variables are also obtained by
direct data processing. They describe weather variability over 25 years for a
particular day. Standard deviations were higher in spring and fall than in
lummer (Fig. 3.2).
The time series for each meteorological variable is reduced to a residual
S'-ries by removing periodic means and standard deviations as pointed out by
l-jchardson (1981). The dependence among the meteorological variables was
described by calculating cross correlation coefficients of the residual time
•cries. The serial correlation coefficients for time lags up to 3 days are given
:n Table 3.1. The serial correlation coefficients for the one day lag were
significant for air temperature (0.69) and dew point temperature (0.66). A
significant cross correlation coefficient (0.8) was calculated for zero time lag
ithe same day) between dew point temperature and air temperature. Other
meteorological variables were uncorrelated for the same day.
The first order estimate of the mean and covariance temperatures is
constrained to parameter perturbations within only the linear region about the
model trajectories. Linear approximation could be questionable when the
coefficient of variation for the parameter of a highly nonlinear function
increases above 0.3 (Gardner et al., 1981). Average coefficients of variation
for input meteorological variables are: air temperature 0.13, dew point
temperature 0.17, wind speed 0.33, and solar radiation 0.37. Although the
43
-------
fN
I
E
o
O
h-
<
O
a
3
o
GO
UJ
a:
LJ
a
2
O
a
UJ
Q
700
600-
500
•too
300-
200-
100-
n
35
30-
25-
20
15-
10
5-
0-
T5-
-10
I —I 1 I I ' * 'I '
— overage
••• standard deviation
i »•"• :;i i- -•
M'R MAV JUN JUL AUG SIP OC1
ovmoge
standard devialion
' i '—•—'~r • • • i • • • i •-• i i • • • i • • • i •
APR MAY JUN JUL . AUG SEP OCT
(
AI'H MAY JUN JUL AUC SL'P OCI
..".:'',V«.;i. .1. > '!.-.
< '( V ^ '' (' : ^'•!"\ft,-.l.i..;: ;V /Vi-V.A'-'^^A:"'^ '
-------
Table 3.1 Correlation coefficients of daily meteorolo/;icnl viiriahlrtt f..i
Minneapolis-St. Paul, 1955-1979.
Solar
Radiation
Air
Temperature
Dew Point
Temperature
Wind
Speed
Solar
Radiation
Time Lag (days)
01 2
1.00 0.39 0.14
0.18 x 0.17 0.11
-0.25 -0.14 -0.06
-0.16 -0.05 -0.04
Air
Temperature
Time Lag (days)
0 1 2
1.00 0.69 0.38
0.80 0.54 0.26
0.11 0.08 0.07
Dew Point
Temperature x
Time Lag (days)
0 1 2
0.80 0.58 0.29
1.00 0.66 0.33
0.09 0.07 0.06
Wind
Speed
Time Lag (days)
0 1 2
1.00 0.38 0.18
Cn
-------
solar radiation had the highest variability note that it is linearly related
through the source term to the water temperature equation (Equations 1.1,
1.2, 1.3, and 1.4)
3.4.1 First moment analysis
The nonlinear lake temperature model was used for the first moment
temperature estimation. Model setup parameters which are basically related
to lake geometry have been estimated by comparing model simulations with
measurements (Chapter 2). The standard error between measurements and
simulations was about 1.0 °C. The error is mostly associated with small
differences between measured and predicted thermocline depths.
Long term average temperature structure in Lake Calhoun was obtained
using two different methods. In the first method, the lake temperature model
computed the vertical temperature structure in the lake for each of twenty
five years (1955-1979), separately using daily values for meteorological data.
The results of these twenty five years of simulated lake temperatures were
statistically Analyzed in terms of mean (Teav) and standard deviation (oeav)
for the particular day. In the second method, twenty five years of daily
meteorological data were first statistically analyzed to provide daily means
and standard deviations. This averaged meteorological year was used in a
single simulation run to obtain the average daily water temperature (Tav)
throughout a season.
Epilimnetic and hypolimnetic lake temperatures obtained by these two
methods are compared in Fig. 3.3. Epilimnion temperature is defined as the
temperature of the upper isothermal (mixed) layer. Hypolimnetic temperature
is a volume weighted average temperature below the upper isothermal layer
down to the lake bottom (Equation 2.5). Nearly identical temperature
distributions were obtained by the two methods. Maximum difference was
less than 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
-------
35-
0-
O
o
3CH
'eav 1955-1979
T,
eav 1955-1979
av
0 r
EPILIMNION
KYPOLIMNION
APR MAY JUN JUL AUG SEP OCT
Fig. 3.3 Estimated long-term average epilimnion and hypolimnion
temperatures.
47
-------
24 ^
APR MAY JUN JUL AUG SEP OCT
Fig. 3.4 Long-term average isotherms in Lake Calhoun.
48
-------
Standard deviations of simulated epilimnion and hypolimnion
'.rmperatures are given in Figs. 3.5 and 3.6, respectively. Contributions by
-".-rturbations of individual meteorological variables perturbations as well as
•.he total contribution of all perturbation variables were calculated with
or related and uncorrelated input variables at a daily timestep. Air and dew
>--:nt temperature contributed the most to the temperature uncertainty, while
viiar radiation and wind speed had smaller effects. Furthermore, the overall
uncertainty in water temperature was found to be larger in the case of the
Correlated daily process than in the uncorrelated one. Uncertainty in lake
•ater temperature varies with time, since sources of uncertainty vary with
iirae. These sources are, the sensitivity of lake temperatures to
meteorological variables as well as the amount of the uncertainty concerning
these variables. At the beginning of the simulated period uncertainty was set
;.o zero since initial conditions were considered perfectly known. Isothermal
initial conditions of 4"C (after ice thaw) April 1 are appropriate for the 45°
.latitude. Although isothermal water at 4'C may not exactly exist on April
i, thermal inertia of the water makes summer predictions insensitive to initial
nnditions if a starting date at or before "ice-out" is chosen (Chapter 2).
I'hree periods can/be distinguished in Fig. 3.5 : a steep rise in temperature
•.mcertainty in spring, more or less constant uncertainty after onset of
* i ratification in summer, and decreasing uncertainty in fall when lake
•.cmperature is driven towards isothermal conditions. Temperature uncertainty
:s decreasing in fall when observed meteorological variables and estimated lake
*ater temperatures are both decreasing. First order derivatives with respect
•u the lake temperatures and meteorological variables are evaluated at these
»bserved and estimated values respectively. Thus, they have less weight in
uncertainty propagation.
Uncertainty propagation for deep hypolimnetic temperature (1m above
lake sediments) is given in Fig. 3.6. In spring and fall, during well-mixed
conditions (overturn periods), standard deviations of 0.4 °C (correlated case)
and 0.3 °C (uncorrelated case) are calculated. During stable stratification,
uncertainty was not significant. This is a result of the fact that Lake
Calhoun has no significant continous point inflows (tributaries). .Summer
temperature in the hypolimnion was determined by mixing events in spring,
and remained almost constant throughout the fall overturn (Ford and Stefan,
1980a). In lakes with point inflows this would not be the case, due to
plunging flow phenomena.
Vertical profiles of the first moment lake temperatures, plus or minus
one standard deviation interval, are shown in Fig. 3.7. Spring (April) and
fall (October) indicated periods when uncertainty propagates throughout the
entire lake depth. These are the periods of weak stratification or well mixed
conditions. Uncertainty was decreasing with depth. After the onset of
stratification estimated uncertainty was much more significant for the
epilimnetic layer than for the hypolimnetic layer. For the same period of
time, deep water had insignificant lake temperature uncertainty.
The first moment epilimnion temperature estimates plus or minus one
standard deviation obtained by two different approaches for the 1955-1979
period are compared in Fig. 3.8. In the first approach the deterministic
water temperature model" was run for 25 years using daily meteorological
49
-------
STANDARD DEVIATION OF SIMULATED EPILIMNION TEMPERATURES
6.0-
O 5.0-
o
LU
rr
Z)
4.0-
rr
.UJ 3.0-
o_
UJ
"~ -2.0-
rr
UJ
S 1.0H
0.0-
all inputs
air temperature ,
dew point temperature
solar radiation
wind speed
UNCORRELATED DAILY PROCESS -
O 5.0-
o
UJ
cr
LJ
Q_
LJ
— all inputs
— air temperature
•-- dew point temperature
— solar radiation
— wind speed
CORRELATED DAILY PROCESS
0.0
APR MAY
JUN JUL
AUG
.SEP
OCT
Fig. 3.5 Standard deviations of estimated epilimnion temperature
uncertainties. Contributions by several meteorological
variables and totals are shown.
50
-------
STANDARD DEVIATION OF SIMULATED DEEP HYPOLIMNION TEMPERATURE
O.c
LU
Ct
ID
<
C£
UJ
CL
cz
UJ
all inputs
air temperature
dew point temperature
solar radiation
UNCORRELATED DAILY PROCESS
wind speed
O.-H
•0.2-
O.C
O
o
all inputs
CORRELATED DAILY PROCESS
air temperature
dew point temperature
solar radiation
wind speed
APR
OCT
Fig. 3.5 Standard deviations of estimated, deep water temperature
uncertainties. Contributions by several meteorological
variablesand totals are shown.
51
-------
4-
8-
12-
1S-
20-
24
4/10
6/08
4.
a-
12-
,'
16-
20-
24
S/M
I
t^
Q.
U)
o
4-
s-
12-
16-
10-
24'
B/24
12-
16-
20-
34
10/28
5 10 15
CC)
10 IS 23 25 X
TEMPERATURE fC)
Fig. 3.7 Long-term average temperature profiles plus or minus one
standard deviation in Lake Calhoun.
52
-------
Cn
CO
40-
^ 35-1
O
o
"~" 30-
Ld
cr
P 25-
r
<
o:
UJ 20-
CL
10-
5-
0
Teov (1955-1979)-f(~)^"eav
EPILIMNION
I I T I ^| | i i| I T^ 1^I I ^ i 1 j^ t 1 I | I F I T T
APR MAY JUN JUL AUG SEP OCT
Fig. 3.8 Epilimnion^temperature long-term average plus or minus one
standard deviation.
-------
data. Long term average (Tgav) and standard deviations (aeav) were
estimated from the simulated lake water temperatures over the 25 year
period. In the second approach the long term average (Tav) temperature
structure in the lake was estimated using the method described in Section
5.1. Water temperature variability (<7av) was estimated using the proposed
perturbation model. Results shown are for correlated meteorological
perturbation variables. The maximum difference was less than 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
instrumentally recorded air temperatures (Kerr, 1989). Uncertainty analysis
of the effects of variable meteorological forcing on lake temperature models
indicates that air temperature has the most significant effect on lake
temperature uncertainty (Henderson-Sellers, 1988; Chapter 3). 1971 was
normal in the sense that mean air temperature from May to September was
only 0.2* C below the normal from 1941 to 1970. The effects of the 1988
(warmer) and the 1971 (normal) summer climate on temperatures and
stratification in the three lakes are reported herein.
55
-------
4.2 Method of Lake Temperature Modeling
The test lakes, Lake Calhoun, Lake Holland, and Lake Elmo, are three
temperature zone dimictic lakes. Water temperature data were collected in
Lake Calhoun in 1971 (Shapiro and Pfannkuch, 1973) and used to validate
the model for normal weather conditions. For warmer conditions (1988) the
model was validated with data from Lake Elmo and Lake Holland (Osgood,
1989). The terrain in which the lakes and weather stations are located is
flat and quite uniform with respect to land use (residential and park land).
Morphometric characteristics, Secchi-depths and chlorophyll-a measurements
for all three lakes in 1984 (Osgood, 1984) are given in Table 4.1. Lake Elmo
surface area is equal to the median value of 970 statistically analyzed lakes
in the North Central Harwood Forests ecoregion in Minnesota (Heiskary and
Wilson, 1988). Lake Calhoun and Holland Lake have a larger and smaller
surface area than the median, respectively. All three lakes were classified as
eutrophic. Secchi depths and chlorophyll-a were close to the median values
of the lakes in the ecoregion.
. Table 4.1 Lake data
Lake
Calhoun
Elmo
Holland
Mean
depth
[m]
10
13.4
4.6
Max
depth
[m]
24.0
41.8
18.8
Surface
area
[km*]
1.71
1.23
0.14
Volume
[10«m3]
17.1
16.5
0.65
Secchi
Depth
M
2.5
2.8
2.2
Typical
Summer, Chla
[gm-3]
20
8
28
Meteorological data used are from the Minneapolis-St. Paul Internationa]
Airport located 5 to 18 miles from the studied lakes. The meteorological
data file contains measured daily values of average air temperature (Ta), dew
point temperature (T.'
duration have been continuously recorded on Lake Mendota (Wisconsin, C
latitude) since 1855. The mean date of ice thaw was April 5 with 11 cU*-
standard deviation (Robertson, 1989). Model sensitivity to the date of :>
initial isothermal conditions is summarized in Table 4.3. Epiiiranr.;
temperatures are very well simulated throughout the entire summer per.--*
56
-------
Table 4.2 Mean Monthly Meteorological Data
Max Min Aver. Diff. from
Normal*
Td P W HS
[•C] [mm] [ms'1] [cal cm^d'1]
U'R
MAY
;r.\
-TL
U.'G
• Y. P
CT
14.9
19.4
27.4
26.5
27.5
22.9
15.7
1.7
6.5
16.5
14.3
14.2
11.3
5.8
8.3
13
21
20
20
17
/io
.0
.9
.4
.9
.1
.8
Year
0.6
-1.3
2.0
-2.3
-0.8
1.2
0.5
1971
-1
2
15
12
12
10
6
.7
.8
.0
.8
.2
.0
.7
0.9
2.6
3.0
3.2
1.4
2.3
4.6
5
4
4
4
3
4
4
.1
.6
.0
.1
.8
.1
.5
411
482 .
450
563
479
338
192
MEAN(MAY to SEPT)
, 24.7 12.6 18.7
(MAY to SEPT)
3.15 3.45
3.26
-0.2
1.60
10.6
2.5 4.1
4.20 0.63 0.26
MEAN(MAY - SEPT)
28.2 15.0 21.6 2.7 11.8 1.7 4.8
a (MAY - SEPT)
3.42 3.23 3.29 1.20 3.03 1.16 0.23
462
72.7
APR 15.1
MAY 25.7
JUN 30.5
JUL 32.3
AUG 29.3
SEP 22.8
OCT 12.5
2.0
11.4
16.6
18.8
17.3
10.9
0.8
8.5
18.5
23.5
25.6
23.3
16.9
6.7
Year
0.8
4.2
3.6
2.9
1.6
1.0
-3.6
1988
-2.4
7.1
12.3
14.5
15.3
9.8
-0.6
1.3
1.4
0.2
0.9
3.5
2.4
0.7
4.6
5.1
4.8
4.4
4.8
4.8
4.7
469
584
654
610
497
331
284
535
114
*Normal is the 30-year average from 1941 to 1970
a — standard deviation
57
-------
Or
OO
Table 4.3 Differences (°C) in simulated mean daily epilimnetic and hypolimnetic temperatures for different
starting dates of the model (April 1 reference)
MAY
JUN
JUL
Epilimuiou
AUG SEP OCT SEASON
MAY
JUN\
Hypolimnion
JUL AUG SEP
OCT
SEASON
Year 1971
MAR 1 -0.05
MAR 12-0.09
MAR 22-0.09
APR 10 0.05
APR 20 0.91
-0.07
-0.08
0.07
0.01 v
-0.03
0.00
0.00
0.01
0.00
0.00
0.00
0.00
-0.01
0.00
0.00
0.00
0.00
-0.05
0.00
0.02
0.00
-0.01
-0.03
0.00
0.08
-0.02
-0.03
-0.02
0.01
0.16
-0.19
-0.29
-0.18
0.08
1.77
-0.24
-0.34
-0.18
0.03
1.59
-0.28
-0.39
-0.18
0.01
1.46
-0.30
-0.41
-0.19
0.00
1.35
-0.32
-0.43
-0.19
0.00
1.26
-0.33
-0.45
-0.19
0.01
1.20
-0.28
-0.39
-0.18
0.02
1.44
Year 1988
MAR 1 0.14
MAR 12 0.07
MAR 22 0.07
APR 10 0.70
APR 20 1.60
0.00
0.00
0.00
0.03
0.08
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
-0.02
-0.03
0.01
0.00
0.00
-0.04
-0.07
0.03
0.01
0.01
-0.06
-0.15
0.03
0.01
0.01
0.10
0.24
0.48
0.01
0.29
1.42
3.14
0.49
0.13
0.30
1.46
2.93
0.51
0.10
0.32
1.44
2.79
0.52
0.10
0.33
1.43
2.77
0.52
0.10
0.33
1.42
2.67
0.44
0.08
0.28
1.40
2.59
0.49
0.10
0.31
1.43
2.82
-------
regardless of the starting date of the model. Surface water temperatures
''catch up" in time. Hypolimnetic summer water temperatures are good as
long as the model is started before seasonal stratification sets in. Better
results are obtained if temperature is not allowed to drop below 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
insensitive to initial conditions if a starting date at or before "ice-out" is
chosen.
4.3 Model Validation
The model was validated with water temperatures measured in Lake
Elmo and Holland Lake in 1988. Eight examples of measured and calculated
water temperature profiles for these lakes are given in Figs. 4.1 and 4.2.
Actually 16 profiles were measured in each lake. Simulations started with
isothermal conditions (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 (Hbr) cool the water. Conductive heat transfer (Hc) can either heat
or cool the water. These five mechanisms, mainly responsible for the net
heat energy input to the water, changed from month to month and from year
(1971) to year (1988). Solar radiation (Hsn) and atmospheric radiation (Ha)
are only given once because they are the same for all three lakes.
Cumulative heat balance terms for both simulated periods are given in Table
4.5.
Under warmer conditions (1988) more solar radiation reached the lake
surfaces. The cumulative difference at the end of the simulation period was
5000 kcal m'2. The additional available solar radiation increased the
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 h
c.
Ul
o
40 h
Fig. 4.1 Lake Elmo water temperature profiles.
60
-------
TEMPERATURE (°C)
10 -
20 -
o MEASURED
— CALOAATED
20 -
Fig. 4.2 Lake Holland water temperature profile.
61
-------
TEMPERATURE (°C)
o (-
10 -
20 -
OJ
Q
O MEASURED
— CALCULATED
10 -
20 -
Fig. 4.3 Lake Calhoun water temperature profiles in 197.1.
62
-------
Table 4.4 Monthly averages of daily heat balance components [1000 kcal
H8n
Lake Calhoun
Ha Hbr He
Hc Hn
Lake Elmo
Hbr He HC
Hn
Lake Holland
Hbr He HC
Hn
Year 1971
APR
MAY
JUN
JUL
AUG
SEP
OCT
MEAN
3.89
4.60
i
4.27
5.37
4.54
3.18
1.83
(MAY to
4.39
5.76 -6.92
6.27 -7.64
7.44 -8.58
7.11 -8.85
7.20 -S.71
6.87 -8.36
\
6.18 -7.69
SEPT)
6.98 -8.43
-1.04
-1.86
-2.00
-3.30
-2.68
-2.14
-1.39
2.39
0.34
-0.08
0.07
-0.46
-0.20
-0.25
-0.42
-0.18
2.03
1.28
1.20
-0.13
0.15
0.70
-1.49
0.36
-6.88 -1.04 0.44
-7.49 -1.65 0.13
-8.44 -1.74 0.25
-8.82 -3.47 -0.46
-8.67 -2.76 -0.17
-8.35 -2.26 -0.25
-7.71 -1.55 -0.49
-8.35 -2.38 -0.10
X 2.17
1.85
1.77
-0.27
0.13
-0.81
-1.75
0.53
-7.03
-7.76
-8.63
-8.83
-fl.71
-8.30
-7.49
-8.45
-1.37
-2.22
-2.24
-3.38
-2.79
-2.10
-1.09
-2.55
0.17
-0.25
0.01
-0.45
-0.21
-0.19
-0.15
-0.22
1.41
0.64
0.84
-0.19
0.04
-0.54
-0.73
0.16
o>
co
Year 1988
APR
MAY
JUN
JUL
AUG .
SEP
OCT
MEAN
4.46
5.57
6.27
5.83
4.72
3.11
2.69
(MAY to
5.10
5.76 -7.08
6.84 -8.06
7.52 -8.99
7.84 . -9.17
7.56 -9.00
6.79 -8.23
5.54 -7.50
SEPT)
7.31 -8.69
-1.29
-2.49
-^J.40
-4.06
-3.74
-2.32
-1.97
-3.40
0.16
0.32
-0.18
0.03
-0.21
-0.27
-0.79
-0.06
2.01
2.18
0.22
0.47
-0.67
-0.92
-2.03
• •
0.26
-6.96 -1.15 0.39
-7.84 -2.00 0.71
-8.87 -4.29 -0.04
-9.11 -4.14. 0.11
-8.99 -3.98 -0.21
-«.24 -2.55 -0.31
-7.39 -1.86 -0.47
-8.61 -3.39 0.05
2.49
3.28
0.59
0.52
-0.89
-1.20
-1.44
0.46
-7.25
-S.21
-9.01
-9.14
-8.95
-8.16
-7.30
-8.69
-1.71
-3.07
-4.63
-4.15
-3.70
-2.21
-1.65
-3.55
-0.12 .
0.10
-0.20
0.05
-0.15
-0.18
-0.48
-0.08
1.13
1.23
-0.06
0.43
-0.52
-0.65
-1.21
0.09
-------
air
was
Atmospheric long wave radiation and back radiation from the water
surface are proportional to the fourth power of absolute temperatures. Both
were higher under wanner conditions. Higher back radiation was an
indication of higher surface water temperatures under increased air
temperatures and solar radiation.
Cumulative evaporative losses resulting from the average 2.9'C
temperature increase are plotted in Fig. 4.4. Cumulative evaporative loss
higher by about 180,000 kcal m'2 for the 1988 season compared to 1971.
This translates into an additional water loss of about 0.3 m in 1988
compared to 1971. This loss occurred in each of the three lakes despite their
differences in size and depth. Increased evaporation not only represents an
additional water loss but also contributes to increased natural convection due
to surface cooling.
Conductive heat transfer through the lake surface made only a small
contribution to the heat budget. The cumulative conductive heat input was
not significantly different during the two years, but the onset of cooling by
convection was delayed until August in 1988.
/
Net heat fluxes on a monthly time scale are shown in Table 4.4, and
on a cumulative basis in Table 4.5. Cumulative net heat flux (Hn) from the
atmosphere to the water increased from April to June in 1971 and from April
to July in 1988, and then began to decrease indicating that the lakes received
heat for a longer period in 1988 than in 1971. The net cumulative heat
input is also a measure of heat content relative to April 1. The maxima of
the net cumulative heat input were only slightly different in 1971 and 1988
(see Table 4.5), but very different among the three lakes because of the effect
of depth especially surface mixed layer depth. Normalized values with respect
to depth are given in Table 4.6. The trend is from higher to lower values as
the depth increases. This reflects the thickness of the surface mixed layer
depth relative to the total lake depth.
4.4.2 Equilibrium temperatures
Equilibrium temperature is defined as that water temperature at which
the net rate of heat exchange through the water surface is zero and
continually changes in response to the meteorological conditions. Mean
monthly equilibrium temperatures for Lake Calhoun are shown in Fig. 4.5.
These values were obtained by a separate calculation setting the net heat
transfer rate Hn equal to zero. Calculations were carried out for the entire
year (12 months) to see how the dates of the 0°C crossings and hence the
date of ice formation might shift from year to year. Under warmer
conditions equilibrium temperature was higher from March to August. From
August up to the ice formation in November no difference between the colder
and the warmer year was noticed probably because the fall of 1988 was
cooler than in 1971 (see Table 4.2). The 0*C crossings in Fig. 6 occurred at
about the same time in 1988 and 1971 indicating that dates of ice formation
and ice thaw were not significantly affected by the heat in July and August.
There could be a larger change in ice thaw and freeze-over dates if air
temperatures were changed year—around, not only in summer as in this case
study.
64
-------
Table 4.5 Cumulative heat balance components [1000 kcal m"2]
Lake Calhoun
Hgn Ha Han
He
Hc
Hn
Lake Elmo
He Hn
Lake Holland
He Hn
Year 1971
APR
MAY
JUN
JUL
AUG
SEP
OCT
117
259
387
554
694
790
846
173
367
590
810
1034
1240
1431 »
-208
-444
-702
-976
-1246
-1497
-1735
-31
-89
-149
-251
-334
-398
-441
10
8
10
-4
-11
-18
-31
61
101
137
132
137
116
70
-31
1
-82
-135
-242
-328
-396
-444
65
122
175 ,
167
171
147
93
-41
\
-110
-177
-282
-369
-431
-465
42
62
87
81
82
66
44
01
Year 1988
APR
MAY
JUN
JUL
AUG
SEP
OCT
134
306
495
675
822
915
998
173
385
610
863
1088
1291
1463
-212
-402
-732
-1016
-1295
-1542
-1775
-39
-110
-248
-374
-490
-559
-620
5
15
9
10
4
-4
-29
60
128
134
149
128
101
38
-35
-97
-225
-354
-477
-553
-611
75
176
194
210
183
147
102
-51
-M7
-286
-414
-529
-595
-647
34
72
70
84
68
48
11
-------
SIMULATION PERIOD 4/01-10/31
800000-
700000-
5 600000-
500000-
100000
0
^ 700000-I
J5 600000-1
*
g 500000 -|
S J
5 400OOO-1
o
$
S 300000-
LU
>
3 200000-
3
§
0 100000-
D
^ 700000-
"2
5 600000-
2.
to
cl 500000-
—J
LJ
>
^ 400000-
cr
O
CL
§ 30C-000-
LJ
3 200000-
3
2
O
0 100000-
0
1971 losses (kcol/m2)
1988 losses (kcol/m2)
1971 losses (m)
1988 losses (m)
LAKE CALHOUN
""[" "* | ' I _
1971 losses (kcol/mz)
1988 losses (kcol/m2)
1971 losses (m)
1988 losses (m)
~\ ' T ' f ~
1971 losses (kcol/m^)
1988 losses (kcol/m2)
1971 losses (m)
1988 losses (m)
LAKE HOLLAND
1.2
p.1
j-i.o
-0.9 ~
E-0.8 «
E-o.s >
r" I
rO.3 3
-0.2
=-0.1
1-1.1
r°-8 $
t o
r-0.6 0
-0.4
-0.3
-.,
E-o.i
o.o
APR
MAY
JUN JUL AUG SEP OCT
Fig. 4.4 Cumulative evaporative losses (simulated).
66
-------
30-
o>
O
o
bJ
0£
i
Ul
Q.
2
UJ
20-
15-
10-
5-
-5-
-10-
1988-89
1971-72
FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC JAN
Fig. 4.5 Mean monthly equilibrium temperatures (simulated).
-------
Table 4.6 Net cumulative heat input (content) per meter of average depth
(1000 kcal m'1)
Lake Holland
: (4-6 m)
Lake Calhoun
(10 m)
Lake Elmo
(13.4 m)
Year 1971
APR
MAY
JAN
JUL
AUG
SEP
OCT
11
16
22
20
21
17
11
6
10
14
13
14
12
7
5
10
14
13
14
12
7
Year 1988
APR
MAY
JUN
JUL
AUG
SEP
OCT
8
18
18
21
17 •
12
3
6
13
13
15
13
10
4
6
14
16
17
15
12
8
68
-------
4.4.3 Vertical mixing/onset of stratification
Surface mixed layer depths are shown in Fig. 4.6. The mixed layer
depth is defined as the thickness of the upper isothermal layer. Large mixed
layer depths at the beginning and at the end of the simulated period indicate
spring and fall overturns. After ice-out in spring, mixing depths were high,
i.e. temperature was uniformly distributed throughout the entire lake. That
was also the justification for selecting April as the initial time for numerical
simulations.
In summer mixed layers were deeper in two of the three lakes under
warmer (1988) conditions. Increased net heat flux to the lake caused a
slightly earlier onset of stratification. The simulated onset of stratification is
first observed in the smallest lake (Holland Lake). Lake Calhoun and Lake
Elmo started to stratify later and showed similar mixing events on a daily
timescale. Vertical mixing is caused bv wind and natural convection.
Surface mixed layers were deeper in Lake Elmo, mainly because more wind
energy was available for entrainment at the thennocline due to the longer
fetch (greater surface area) of the lake. Natural convection is mainly driven
by net surface (evaporative, conductive) heat loss. Under warmer conditions
evaporative loss was much higher, ana 1 to 2 m deeper mixed layers were
probably produced in this way.
Fall overturn occurred earlier after the wanner 1988 summer because
lower fall temperatures produced stronger cooling and surface water
instabilities, i.e. thermals and convective negatively buoyant (cold) currents
earlier (Horsch et al.,1988). In the presence of convective cooling, less
turbulent kinetic energy, supplied by the wind, is needed for the deepening of
the thermodine.
4.4.4 Water temperatures
Daily epilimnetic temperatures at a depth of 1.5 m are shown in Fig.
4.7. Although lakes have different morphometries, similar temperature
patterns were observed. This is in agreement with field measurements made
by Ford and Stefan (1980) in 1974 and 1975. In both 1971 and 1988 the
surface temperatures of the three lakes exhibited similarities and parallel
trends which are predominantly related to weather phenomena and only
secondarily to lake morphometry (Ford and Stefan, 1980). From April
through August epilimnion water temperatures were higher in 1988 (average
water temperature increase « 3°C compared to 3'C in air temperature
change) and lower after the lake staried cooling.
Daily hypolimnetic temperatures are shown in Fig. 4.8. Values are at
depths well below the thennocline, and water temperatures are nearly
isothermal below that depth. Lake Calhoun and Lake Elmo received
additional heat during the spring turnover periods (Fig. 3.6) when the climate
was warmer (1988). Average hypolimnion temperature was 0.6 and 1.4* C
higher in 1988 in Lake Calhoun and Lake Elmo, respectively, than in 1971.
Lake Holland experienced an opposite trend. The lake started to stratify
earlier too, but due to the increased stability and small lake surface area,
wind mixing throughout the entire lake in spring under warmer conditions
69
-------
was weaker. Average hypolimnion temperature was 1.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. -
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
wanner summer. Evaporative water losses increased by about 300 mm out of
800 mm or about 40 percent.
(4) Dates of ice formation in fall seemed only weakly affected by the
hot midsummer weather. Dates shifted by a few days. This may not be
typical because of the cool fall.
(5) Simulated conductive heat transfer had a negligible effect on heat
budget changes.
(6) Higher atmospheric radiation due to higher air temperature was
compensated by higher backradiation from the water.
(7) Simulated surface mixed- layer depths increased slightly (by 1 to 2
m) in the warmer summer, probably due to stronger convective mixing.
(8) Simulated stratification onset occurred slightly earlier in the warmer
year.
70
-------
SIMULATION PERIOD 4/01-10/31
45
40-
35-j
30-j
I 25-|
r
£j 20-
15~
3
5-;
]
o4-
j
LAKE CALHOUN
1971
1988
40-
35 -i
& 20-(
15-
LAKE ELMO
1971
tAKE HOLLAND
40-
1971
196=
jo-:
2 2*-3
& 20-J
;PR K/.-V JUN . JUL AUG SEP CCT
Fig. 4.6 Mixed layer depths (simulated).
71
-------
SIMULATION PERIOPD 4/01-10/31
g
UJ
*-^i*' A
' ^"' A/V\
j^\ MI/ ^/V\/\AA A . •
•v/- • \ :
• , i • i • i • i i i • t • i
"• 1971 LAKE ELMO
--- 1 QOQ
(j7OO -
— DIFFERENCE
i* /*-* '.' ' *"* *» -
,/'»•*.* .•.;*"* V'"\ _....--"*t
/ * """"" * " \ —
'X ;"v' ~"*'\
- ' ./ "l<1*..
/,>. ;*'. ; '' V-^. ..,
* ." V '""" *»
^ . f »"*•*" A —
"' ••"""" / u^ yA ~
Pv Afv y\/^ \. .A
W" ^ \ -
-
• • i • i • i • t • t i
~- 1971 LAKE HOLLAND "
1988
DIFFERENCE „ - , /-, >. __ ,,
/'^^V'-'-'V' '"">.'.
A-; /U ^ ''•\"/'\-
** ' ' • ' *** tir*
,;/l/^ /^ : 'v.'
APR WAY JUN JUL AUG SEP
Fig. 4.7 Simulated epilimnion temperatures.
72
-------
SIMULATION PERIOPO 4/01-10/3!
12-
£• *-
2
....._._ 1971
1988
DIFFERENCE (1988-1971)
LAKE CALHOUN \
i _...
-2-
- 1971
1988
DIFFERENCE
2H
LAKE ELMO
-2'
1971
1988
DIFFERENCE
6-4
LAKE HOLLAND \
-2
APR MAY JUN JUL AUG SEP OCT
Fig. 4.8 Simulated hypolimnion temperatures.
73
-------
5. Water Temperature Characteristics of Minnesota
Lakes Subjected to Climate Change
A deterministic, validated, one-dimensional, unsteady-state lake water
quality model was linked to a daily weather data base to simulate dailv
water temperature profiles in lakes over a period of twenty-five (1955-79)
years. 27 classes of lakes which are characteristic for the north-central US
were investigated. Output from a global climate model (GISS) was used to
modify the weather data base to account for a doubling of atmospheric CO^-
The simulations predict that after climate change epilimnetic temperatures
will be higher but increase less than air temperature, hypolimnetic
temperatures in seasonally stratified dimictic lakes will be largely unchanged
or even lower than/at present, evaporative water loss will be increased by as
much as 300 mm for the season, onset of stratification will occur earlier and
overturn later in the season, and overall lake stability will become greater in
spring and summer.
5.1 Introduction v
This Chapter deals with the question of how climate change may affect
thermal aquatic habitat in lakes. A regional perspective is taken, and the
scope is to estimate temperature changes in lakes of different morphometric
and trophic characteristics in a region. Southern Minnesota is chosen as an
example because an extensive lake database is available (ERLD/MNDNR,
1990). The geographic boundaries of Southern Minnesota are defined in
Figure 5.1.
Herein a dynamic and validated regional lake water temperature model
(Chapter 2) will be applied to a representative range of lakes in a region for
past climate and one future climate scenario. Rather than analyzing
particular years and lakes, emphasis is on long term behavior and a wide
range of lake morphometries and trophic levels. In this study the base
period (or comparable reference) was from 1955 - 1979. For the same period
of time weather parameters were perturbed by the 2XC02 GISS (Goddard
Institute for Space Studies) climate model output. The regional impact of
these climates on different lake classes in southern Minnesota is reported
herein. The simulated water temperatures, past and future, will be presented,
interpreted and related to the lake characteristics and climate characteristics.
The results will show how water temperatures in different freshwater lakes
respond to changed atmospheric conditions in a region.
Lake levels will be largely controlled by the water budget including
evaporation and runoff. The response of watershed (surface) runoff to climate
change is the subject of other investigations not included herein. Lake depths
74
-------
Fig. 5.1 Regional boundaries and-geographic distribution of lakes in MLFD
database.
75
-------
will therefore be treated herein as either invariant or be lowered to account
for increased evaporative water losses, where applicable. Changes in the
watershed may affect nutrient loadings and hence primary productivity and
transparency of the water. Such secondary effects, also were not investigated,
but a sensitivity analysis indicates that water temperature predictions for the
types of lakes studied herein are usually not sensitive to transparency
(Chapter 2).
5.2 Method of Lake Temperature Modeling
The numerical model is applied in daily timesteps using mean daily
values for the meteorological variables. The required weather parameters are
solar radiation, air temperature, dew point temperature, wind speed, wind
direction, and precipitation. Initial conditions, lake morphometry
(area-depth-volume), and Secchi depth have to be provided to use the model.
Simulations were made from spring overturn to fall overturn. Since the date
of spring overturn is unknown, the initial conditions were- set at 4"C on
March 1, and water, temperature was not allowed to drop below 4° C (well
mixed conditions). Although isothermal water at 4'C may not exactly exist
on March 1, the isothermal 4" condition continues until the model simulates
warmer temperatures and the onset of stratification. The summer predictions
are thus made quasi-independent of initial conditions and match
measurements well (Chapter 3). The model is one-dimensional in depth and
unsteady, i.e. it simulates water temperature distributions over depth in
response to time variable weather. Vertical water temperature simulations
are made over an entire season (March 1 to November 30) and in time steps
of one day. The calculated daily water temperature profiles are analyzed
statistically and presented graphically.
The regional water temperature simulation model was validated against
data from nine Minnesota lakes for several years (Chapter 2). The model
simulates onset of stratification, mixed layer depth, and water temperature
well. Root mean square error is 1.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 Airport
(93.13* longitude, 44.53* latitude) were used. The meteorological data file
assembled contains measured daily values of average air temperature, de»
point temperature, precipitation, wind s,peed, and solar radiation from 1955 u
1979 (March - November). The period from 1955 to 79 was chosen becauw
it is long enough to give a representative average of base conditions befcrr
climate warming. In the 1980s warmer than average air temperatures
observed (Jones et al., 1986; Kerr, 1989), and therefore this period
excluded. Sources of climate data were as follows:
76
-------
Climate scenarios were selected following EPA guidelines on global
dimate change effect studies (Robinson and Finkelstein, 1990). Climate
projections by four different models (GISS, GFDL, OSU, UMKO) for the
doubling of atmospheric C02 were provided by NOAA (1990). The monthly
dimate projections by the four models are different from each other and their
explicit effects on water temperature dynamics can be studied for each model
separately. In this study only the GISS projections for the grid point closest
to Minneapolis/St. Paul were used (Table 5.1), as suggested by EPA for
effect studies. The geographical location of this grid point is given in Figure
5.2. A comparison of the mean monthly weather parameter values (for
Minneapolis/St. Paul) projected by the four models shows that the GISS
projections are not substantially different from GFDL and OSU, except for
wind speeds in November. No adjustments were made to those wind speeds,
however, for a lack of a rational basis and because late fall winds do not
affect the summer water temperature dynamics. No interpolations between
grid points were made, following explicit EPA recommendations.
Table 5.1 Weather parameters changes projected by the 2XCO2
climate model output for Minneapolis/St. Paul.
MONTH
JAN
FEB
MAR
APR
MAY
JUN
JUL
AUG
SEP
OCT
NOV
DEC
AIR. TEMP
(diff.'Q*
6.20
5.50
5.20
5.05
2.63
3.71
2.15
3.79
7.02
3.73
6.14
5.85
SOL. RAD.
(Ratio)t
0.92
1.04
0.98
1.03
1.00
0.99
0.98
1.04
1.04
1.12
1.03
0.99
WIND S.
(Ratio)l
0.92
1.12
0.47
0.69
0.67
0.85
0.93
1.00
1.07
2.23
5.00
0.77
REAL. HUM.
(Ratio) t
1.16
1.01
1.13
1.00
1.09
1.01
0.93
1.02
0.90
0.95
1.00
0.98
PRECIP.
(Ratio)t
1.17
1.03
1.28
1.03
1.12
1.08
1.10
0.98
0.70
0.88
0.99
1.24
* Difference = 2XCOz GISS - PAST
t Ratio = 2XCO2 GISS/PAST
The uncertainty of the climate predictions is not the subject of this
paper. It is understood that relative humidity and wind speeds are not well
predicted at the local scale by global climate models. Fortunately, .
uncertainty analysis of the effects of variable meteorological forcing on lake
temperature models indicates that air temperature has the most significant
effect in lake temperature uncertainty (Henderson-Sellers, 1988; Chapter 3),
and that parameter is better predicted than others.
Seasonal distributions of the 25-year average of observed weather
parameters (which were used as model inputs) are shown in Figure 5.3. Past
climate and the 2XCOa GISS scenario were used as inputs to the water
temperature simulations.
77
-------
oo
• MINNEAPOLIS/SLPAUL
9GISS GRID POINT MINNEAPOLIS/St.PAUL
• DULUTH >
BBGISS GRID POINT DULUTH
north dokoto o
o
V)
0)
c - wscons
km
Fig. 5.2 Geographical location of the closest GISS grid points for
Minneapolis/St. Paul and Duluth.
-------
,« I
•• -
MAR APR MAY JUN JUL AUC SEP OCT NOV
MSTP (1955-1979)
2XC02 CLIMATE SCENARIO (CISS - model)
MAR APR MAY JUN JUL AUC SEP OCT NOV
30
o
20 8
-15
10
•5
ui
ui
a.
a
Fig. 5.3 Climate parameters a Minneapolis/St. Paul in the past and under
a 2xC02 (GISS) climate scenario.
-------
5.4 Regional Lake Characteristics
Regional classification of lakes was approached in a variety of ways.
The ecoregion approach was considered first, but found to give too detailed a
picture. The entire state was considered as a regional entity but rejected as
too large because of the diversity of climate. Dividing the state into a
northern and southern region was considered appropriate and not as arbitrary
as might seem because there is a significant gradient in geological,
topographic, hydrological dimatolpgical and ecological parameters across the
mid-section of the state (Baker et al., 1985, Heiskary et al., 1987). The
southern and northern regions are about equal in 'size (Fig. 5.1).
The Minnesota Lakes Fisheries Database, MFLD (ERLD/MNDNR,
1990), which contains lake survey data for 3002 Minnesota lakes, is for the
southern region. The MLFD database includes 22 physical variables and fish
species. Nine primary variables explain 80 percent of the variability between
lakes. These nine variables include surface area, volume, maximum depth,
alkalirJty, secchi depth, lake shape, shoreline complexity, percent littoral area,
and length of growing season. For regional classification of the lakes in this
study, the possible thermal structure (i.e. whether lakes are stratified or not)
and trophic status are of primary concern. Observations in the northern
hemisphere show that onset and maintenance of stratification in lakes is
dependent on surface area and maximum depth (Gorham and Boyce, 1989) as
well &s climatological forcing i.e. solar radiation and wind (Ford and Stefan,
1980). Lake trophic status contributes to solar radiation attenuation and
oxygea balance. Trophic status was assessed by using a Secchi depth scale
(Heisiary and Wilson, 1988) related to Carlson's Trophic State Index
(CarUon, 1977). Secchi depth information was available in the MLFD.
A statistical analysis of southern and northern Minnesota Lakes in the
MLFD in terms of surface area, maximum depth and Secchi depth was made.
The ceographic distribution of different classes of lakes in Minnesota is given
in Fipre 5.4. Cumulative frequency distributions shown in Figure 5.5 were
used **o subdivide all lakes into three ranges of surface area, maximum depth
and fecchi depth, as shown in Table 5.2. These represent 27 classes of lakes
in e^'h of the two regions of the state. "
Table 5.2 Lake classification
Key -Cumulative Description
Pa;an«ter Range Frequency Class Value
>m2) < 0.4 Lower 30% 0.2 Small
0.4-5 Central ^60% 1.7 Medium
> 5 Upper 10% 10 Large
Depth < 5 Lower 30% 4 Shallow
) 5-20 Central 60% 13 Medium
> 20 Upper 10% 24 Deep
cci- Depth < 1.8 Lower 20-50% 1.2 Eutrophic
(•z 1.8 - 4.5 Central 20-50% 2.5 Mesotrophic
>4.5 Upper 0-10% 4.5 Oligotrophic
80
-------
J\
I I I I I I I I I i
SECCM DEPTH > 4i METERS
AREA > 5.0 SO. KU. cinL»n
Fig. 5.4 Geographic distribution of lakes according to key parameters:
Secchi depth, maximum depth, and surface area.
81
-------
100
10-*
10"
10V
AREA (km2)
10'
100
MAXIMUM DEPTH (m)
100
2468
SECCHI DEPTH (m)
10
Fig. 5.5 Cumulative distributions (%) of key parameters in Minnesc1.;
lakes (from MLFD database).
82
-------
A representative value for surface area, maximum depth and Secchi
depth was chosen in each lake class as input to the model simulation. Those
values are shown under the heading "class" in Table 5.2.
Representative area-depth relationships for three different lake classes
(by surface area) were obtained from 35 lakes which covered the entire range
of distributions in a set of 122 lakes (Figure 5.6).
After areas are expressed as fractions of surface area and depths are
expressed as fractions of maximum depth, an equation of the form
Area = a • exp(b- Depth) -f c (5.1)
is fitted to the data and subsequently used in the simulation as a
representative area-depth relationship. Coefficients a, b; c, calculated by
regression analysis are given in Table 5.3. This procedure is equivalent to
self-similarity of depth-area relationships within a given class.
• .
Table 5.3 Morphometric regression coefficients
in the area vs. depth relationship.
Area
Small
Medium
Large
1.19
1.14
1.14
-1.76
-2.10
-2.91
-0.20
-0.15
-0.08
Lake basin shape was assumed circular for the purpose of wind fetch
calculation. The water temperature simulation results were shown to be
insensitive to these assumptions of morphometric self-similarity and basin
shape.
5.5 Simulated lake water temperature regimes for
historical and future weather
5.5.1 Water temperatures
Simulations of daily water temperature profiles from March 1 to
November 30 (275 days) in each year from 1955 to 1979 (25 years) were
made for each of the 27 lake classes. In addition to lake morphometric
input, i.e. surface area, maximum depth and depth-area relationship, these
simulations used actually recorded daily values of weather parameters, i.e.
solar radiation, air temperature, dew point temperature, wind speed, and
precipitation for each day simulated. A massive weather-database had to be
developed prior to the simulations. The calculated output of 185,625 vertical
water temperature profiles, each consisting of 24 water temperature values,
provided base line information on lake characteristics during a period of the
past when little climate change occurred.
83
-------
SMALL LAKES
G-O Corver
©-® Crystol
Deep
G-Q Forquhor
E-ffl Islond
McCorrons
McDonough
i-A Long-N
Lucy
O-O Porkers
«—* Reitz
Sucker
Twin-Middle
Wolsfeld
Aubum Eost
C (D Auburn West
Bold Eagle
Or-EJ Boss
£-83 Bovorio
Big Cornelian
George
Colhoun
Crystol
Elmo
Harriet
Eagle
LARGE LAKES
3-O Woconia
Z O Big Marine ~
Coon
3-Q Whits Bear ~
9-3 Greenwood :
0.1 0.2 0.3 0.4. 0.5 0,6 0.7 0.8 0.9 1.0
0.0
0.0
Fig. 5.6
Horizontal area vs. depth relationship for lakes.
and depth are normalized.
Area
84
-------
To simulate possible future water temperature regimes, the monthlv
corrections specified by the 2XC02 GISS model scenario were applied to the
weather data base and the simulations were repeated.
From these simulated water temperature data bases under historical and
future climates, each consisting of 4,455,000 water temperature values, the
following characteristics were extracted.
Epilimnetic water temperatures were defined as water temperatures at
1.0 m below the water surface regardless whether maximum-depth is 4 m, 13
m or 24 m, respectively. The seasonal course of epilimnetic temperatures,
averaged weekly over 25 years is shown in Figure 5.7 for both past climate
and the 2XC02 GISS climate scenario. The difference between the two is
also shown in Figure 5.7; the associated air temperature increments due to
climate change were presented in Table 5.1. The largest change in weekly
water temperature change in response to climate change, is on the order of 6
to 7"C, and occurs in spring (April), the minimum is on the order of 0 to
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. Extreme values in lakes of different
geometry vary by no more than 4*C on any given day. Maximum
differentials occur in spring and fall. From June through September, i.e.
during the period of seasonal water temperature stratification, surface water
temperatures in lakes of different morphometic characteristics (depth and
area) are very similar (within 1.0° C). In very large lakes (e.g. the North
American Great Lakes) the significantly greater water volumes and mixed
layer depths cause a substantial lag in heating and cooling leading to water
temperature differences larger than 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.
Hypolimnetic wanning during the summer is dependent on vertical turbulent
diffusion and therefore wind fetch and hence surface area. Dependence of
hypolimnetic temperatures on lake morphometry is very evident in Figure 5.8.
The seasonal pattern of hypolimnetic water temperatures was altered by
climate change most significantly in shallow lakes. All others showed typical
seasonal wanning patterns in response to vertical diffusion.
85
-------
00
0>
35
30
25
I tO
15
i
0
30
uj
a: 20-
CPIIIMHIOH ItMPERAIURCS
SIUULAI10H P'HIOU 1951 - 19''J (I'ASI)
EPIUUNION TCUPCRATUttCS
}XC02 CUUMC JCCKWIO (CISS - model)
CPICIMNION TtMPCRAIUNCS
DIFFERENCt (CISS mod«l - P»ST)
SUAll LAKES
MEDIUM LAKES
(XBCC LAltrS
SMALL LAKES
MEDIUM LAKES
-r-
LARCC IAKCS
I I
i ' • • i • • • I '
SMALL LAKES
• i ' ' • I ' ' • 1
MEDIUM LAKES
LARGE LAKES
-1
•4
MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JUH JUL AUG SEP OCT NOV MAR APR MAY JUN JUL AUG SEP OCT NOV
-I
Fig. 5.7 Simulated weekly epilimnion temperatures.
-------
HVPOL/MNtON
oo
35
30
25
10
IS
10
5
0
30
20-
15-
ID-
S'
0
30
20
15
10
5
0
1955 - 1979 (PASf)
MtPOtlUNION IEWPCRAIURCS
}XC02 CUUA1C SCCMKBIO (CISS - inodel)
KTPOUVINION TEUPCPATURCS
OVTERENCE (CISS rmxltl - PAST)
SUAU. UKCS
UCOIUU UKES
URGE LAKCS
SKUIL LAXCS
\
MEDIUM LAKES
URGE IAXES
*&
fj-v •neiijM wigciropnie
»-• thaltoo •Mlrophte
MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JUN JUL AUG SEP OCT NOV
Fig. 5.8 Simulated weekly hypolimnion temperatures.
-------
Table 5.4 Maximum temperatures of southern Minnesota lakes
OO
oo
PAST 1955-1979
Maximum
Depth
m
SHALLOW
(4.0)
MEDIUM
(13.0)
DEEP
(24.0)
day = Julian
Area
km2
SMALL
(0-2)
MEDIUM
(1.70)
LARGE
(10.0)
SMALL
(0.2)
MEDIUM
(1.70)
LARGE
(10.0)
SMALL
(0.2)
MEDIUM
(1.70)
LARGE
(10.0)
day when
Trophic
Level
eutrophic
mesotrophic
oligotrophc
eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophc
•> eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophc
eutrophic
meeotrophic
oligotrophic
eutrophic
mesotrophic
oligotrophic
Epilimnion
ec
27.5
27.4
27.3
27.4
27.4
27.3
27.4
27.4
27.3
20.0
26.5
26.6
26.4
26.4
20.5
26.5
26.5
26.6
26.4
26.3
26.1
26.2
26.2
26.1
26.1
26.1
26.1
day
203
203
203
203
203
203
203
203
203
203
206
203
204
207
207
203
206
207
206
204
207
206
206
206
206
206
207
Hypolimnion
°C day
24.9 206
26.8 204
27.0 203
26.2 204
27.0 203
27.1 203
26.5 203
26.9 203
27.1 203
11.9 278
12.8 277
17.5 261
18.7 254
19.9 252
23.0 233
24.0 220
24.6 218
25.5 211
7.3 308
7.4 308
7.8 308
11.6 294
11.8 293
12.6 291
18.2 261
18.4 263
19.4 259
GISS-2XCO2
Epilimnion
°C
29.4
29.4
29.3
29.4
29.5
29.4
29.5
29.6
29.5
28.7
28.7
28.7
28.5
28.6
28.7
28.6
28.7
28.7
28.5
28.3
28.10
28.2
28.1
28.1
28.1
28.1
28.2
day
217
217
217
217
217
217
217
217
217
217
218
218
218
218
218
223
218
218
217
218
220
218
223
223
223
223
218
Hypolimnion
•c
26.2
28.3
29.2
27.5
\ 29.1
29.4
28.3
29.1
29.4
12.6
13.0
17.5
18.2
20.3
25.3
26.0
26.6
27.3
10.3
10.4
10.6
12.8
12.9
13.3
18.4
18.7
20.1
day
229
218
217
205
181
217
181
181
217
289
284
276
274
271
248
233
224
218
305
305
305
291
291
291
276
276
273
DIFFERENCE
Epilimnion
•c
1.9
2.0
2.0
2.0
2.1
2.1
2.1
2.2
2.2
2.1
2.2
2.1
2.1
2.2
. 2.2
2.1
2.2
2.1
2.1
2.0
2.0
2.0
1.9
2.0
2.0
2.0
2.1
(GISS-PAST)
Hypolimnion
•c
1.3
1.5
2.2
1.3
2.1
2.3
1.8
2.2
2.3
0.7
0.2
0.0
-0.5
0.4
2.3
2.0
2.0
1.8
3.0
3.0
2.8
1.2
1.1
0.7
0.2
0.3
0.7
maximum temperature occur
-------
The highest hypolimnetic water temperatures (27.1'C) were calculated
for shallow oligotrophic lakes which are typically polymictic or well-mixed for
the entire simulation period. The lowest maximum hypolimnetic temperatures
(7.3° C) occurred in small and deep eutrophic lakes. Climate change raised
by 0° to 3*C the maximum hypolimnetic water temperature or lowered it by
as much as 3.5"C, depending on the particular stratification dynamics of a
lake.
In addition to long-term changes of water temperatures (Figures 5.7 and
5.8) variations from year to year are also of interest. Unfortunately weather
parameters for the GISS climate scenario were only given as long term
monthly averages. Therefore variability on an annual basis could not be
explored for the GISS scenario. On the other hand, annual weather
information was available for the 1955-79 period, and therefore could be used
to give the range of simulated daily water temperatures. Bands of water
temperatures within the 95% confidence interval are shown in Figure 5.9.
The spread is significant and on the order of ± 3 to 5° C around the mean.
This range is about twice as wide as that due to differences in lake
morphometry (Figures 5.7 and 5.8). This is in agreement with field
measurements by Ford and Stefan (1980) and has some bearing on habitat.
Examples of water temperature structures in typical lakes are given in Figure
5.10.
5.5.2 Thermal energy fluxes
The water temperatures discussed above are, of course, the result of net
heat energy input or losses through the water surface, and vertical
distributions of that heat within the lake. For better understanding of the
water temperatures, it is therefore appropriate to consider, at least, briefly
heat fluxes and stratification dynamics. Simulated net heat flux through the
water surface is plotted in Figure 5.12 for past and future (GISS) climate
conditions.
Five heat transfer processes are responsible for heat input into the
water: short wave solar radiation, long wave atmospheric radiation,
conductive heat transfer, evaporation, and back radiation. Short wave solar
radiation and atmospheric radiation increase the water temperature, while
evaporation and back radiation cool the water. Conductive heat transfer can
either heat or cool the water. All five fluxes together comprise net heat flux
at the water surface.
Individual daily heat fluxes vary dramatically with weather as is
illustrated in Figure 5.11. To keep-track of the extraordinary dynamics and
to explain them would take more space than available here, and may not be
particularly fruitful. As a summary, the cumulative net heat fluxes are
presented in Figure 5.12 for past and future (GISS) climate. The difference
between the two is also shown in Figure 5.12. Lakes with large surface areas
will receive more net heat input (up to 30%) than smaller ones, and in
extremely small lakes the difference is even negative, meaning less heat will
be transferred through the water surface and stored in the lake! All net heat
fluxes are per unit surface area of a lake, not total values.
89
-------
EPIJJUNION TEMPERATURES EUTROPHIC LAKES
PAST 1955-1979
HYPOUUMON TEUPERATURES EUTROPHlC LAKES
PAST IS55-1979
depth d«cp (24 m)
oreo small (0.2 km
depth deep (24 m)
oreo targe (10 km*)
MAR APS MAY JUN JUL AUG SEP OCT NOV MAR APR MAY JUN' JUL AUG SEP OCT NOV
Fig. 5.9 Examples giving range of epilimnetic and hypolimne:
temperatures over a 25 year period (95% confidence interval).
90
-------
to
Fig.S.lOa
Simulated temperature (iiotherml structure in three medium deep (13 m maximum depth) lakes of
large (10 km»), medium (1.7 km') and small (0.2 km') surface area. Isotherm bands are in increments
of 2'C. Simulated water temperatures are for past climate (1955-79) (lop) and the 2XCO} GISS
climate icenario (bottom).
-------
Fig.S.lOb Simulated temperature (isotherm) structure in three medium area (1.7 kmj) lakes of different
maximum depths: shallow (4 m), medium (13 m) and deep (24 ro). Isotherm bands are in increment
of 2'C Simulated water temperatures are for put climate (1056 70) (tup) and the 2XCO] G1SS
climate scenario (bottom).
-------
-------
CO
600000
400000
CUUULA1IVC NO MEAT flXIX
SIMULATION PRIOO 1955 - 1979 (PAST)
CUMULATIVE NET MEAT fUJX
2XC02 CLIMATE SCENARIO (GISS - model)
CUMULATIVE NET HEAT rwx
DlFFtREMCE (CISS model - PAST)
190000
100000
MAR APR MAY JUN JUL AUG SEP OCT NO^MAR APR MAY JUN JUL AUG SEP OCT NOV
I t ' t -100000
MAR APR MAY JUN JUL AUG SEP OCT NOV
Fig. 5.12 Simulated cumulative net heat flux.
-------
Back radiation and evaporation are the main processes by which lakes
lose heat in the summer. Evaporative losses were found to be significantly
increased after climate change (GISS). In all lakes, regardless of depth,
surface area and trophic status, the computed evaporation water losses were
uniformly 0.30 m ( ± 0.01 m) higher (Figure 5.13). In other words, lake
water budgets will be put under significant stress. This increased evaporation
also explains why the water temperature increases after climate change
remains at a relatively moderate 2*C, when air temperature increases by an
average seasonal simulation (March 1,- November 30) value of 4.4 "C.
Evaporative cooling is a key to the understanding of the temperature
responses to changed climate.
-5.5.3 Vertical mixing/Stratification/Stability
Vertical mixing and stratification affect lake water temperature
dynamics. A surface mixed layer depth is defined here as the thickness of
the isothermal layer from the water surface downward. Surface mixed layer
depths are calculated daily by the wind mixing algorithm in the model and
averaged over a week (Figure 5.14). Mixed layer depths at the beginning
and before the end of simulation are equal to the total lake depth and
indicate spring and fall overturns. The most shallow mixed layer depths were
calculated for small, deep, eutrophic lakes based on the classification in Table
5.2. Vertical mixing is caused by wind and natural convection. Surface
mixed layer depths were the shallowest for small lakes because of short fetch.
Smaller wind stresses and hence wind energy inputs axe usually associated
with smaller lake surface area (shorter fetch). In these lakes the smallest
amount of turbulent kinetic energy is available for entrainment of the
thermocline. Wind energy required for entrainment of layers at the
thermocline is inversely proportional to the stability (defined as a density
difference between adjacent layers) of the water column and depth of the
mixed layer. The lowest hypolimnetic temperatures and the highest
temperature (density) gradients were calculated for small, deep, eutrophic
lakes. That was the reason for the smallest mixed layer depths calculated for
these lakes.
For the same morphometric lake characteristics, oligotrophic lakes had
deeper surface mixed layers than eutrophic lakes because of higher penetration
depth of irradiance.
Climate change will impose higher positive net heat fluxes at the lake
surface earlier in the season than in the past. That causes an earlier onset
of stratification. This is in agreement with a conclusion derived by
Robertson (1989) from field data for Lake Mendota. In the period from the
onset of stratification until September, mixed layer depths were projected in
the average 1.2 m smaller than in the past. From the end of September,
mixed layer depths were deeper after climate change, mainly due to stronger
natural convection and higher winds caused by climate change. In spring1 and
summer evaporative losses were also increased by climate change but no
significant persistent cooling occurred because of net heat input from radiation
and convection. The earlier onset of stratification in spring and the mixed
layer depth increase in fall were also found by Schindler et al. (1990) in his
analysis of observations in the EL A. In the EL A mixed layer depths
increased due to transparency increase and increased winds due to reduced
forest cover resulting from increased incidence of forest fires.
95
-------
SIMULATION PERIOD 18SS - l»7» (PAST)
ZXCOj CUMATE SCEN*»0 (CISS - nyxl.l)
MAR APR MAY JUN JUL AUG SEP OCT NOV MAR APR MAY J'JN JUl AUG SEP OCT NOV
Fig. 5.13 Simulated cumulative evaporative losses.
96
-------
MIXED LAYER DEPTHS
2XC02 CLIMATE SCENARIO (GISS - model)
5-
10-
15-
20-
5-
£ 10-
ui 15-
Q
20-
5-
10-
15-
20-
25
MEDIUM LAKES
LARGE LAKES
MAR APR MAY JUN JUL AUG SEP OCT NOV
MIXED LAVER DEPTHS
SIMULATION PERIOD 1955 - 1979 (PAST)
MEDIUM LAKES
" i • • • • i •
X-X medium 0ltgoIropHic\
i eutropnlc
O~O shallow otlgotrophic
eu trophic
APR MAY JUN JUL AUG SE? OCT NOV
Fig. 5.14 Simulated weekly mixed layer depth.
97
-------
The stabilizing effect of the density stratification and the destabilizing
effect of the wind can be quantified using a Lake number (Imberger and
Patterson, 1989):
gSt(l -Zt/zm)
Ln = m (5.2)
Po U, AQ ' (1 - zg/zm)
where g is acceleration due to gravity (m s~2), zt is height from the lake
bottom to the center of the thermocline (m), zffl is maximum lake depth (mj,
zg is the height of the center of volume of lake, A0 is lake surface area (m2),
po is hypolimnion density (kg m"3), St is the stability of the lake (kg m;
Hutchinson, 1957), u^ is surface shear velocity (m s"1). Estimates for the
different elements in the Lake number are obtained from daily lake water
temperatures simulations, daily meteorological data, and lake geometry.
Larger Lake number values indicate stronger stratification and higher stability
i.e. forces introduced by the wind stress will have minor effect. Lake number
dependence on lake/area, depth, and trophic status, for different lake classes
is given in Figure 5.15. Stability is higher for oligotrophic lakes than
eutrophic lakes. Oligotrophic lakes had deeper thermoclines and required
greater wind force in order to overturn the density structure of the water
column. Climatic change caused higher lake numbers, i.e. more stable
stratification among the same lake classes.
Seasonal stratification is defined herein as the condition when
temperature difference between surface and deep water temperature is greater
than 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
-------
i-» .otrphic PAST<1955-1979)
CHD shollo. oSgotrophi* PAST (1955-197? B:»«ALL LAKES
ihollow eutrophic 2XC02CtSS
80- X-K stxiuow digotropnic 2XCO2»SS ,,
10 15 CD 25 30
Fig. 5.15 Simulated lake
status.
ce numbers as a function of lake depth and trophic
99
-------
Table 5.5 Seasonal stratification characteristics of southern Minnesota lakes
o
o
LAKE CHARACTERISTICS
MAXIMUMSURFACE
DEPTH AREA
m
SHALLOW
(4.0)
MEDIUM
(13.0)
DEEP
'?•! 0)
km2
SMALL
(0.2)
MEDIUM
(1.7)
LARGE
(10.0)
SMALL
(0.2)
MEDIUM
(1.7)
LARGE
(10.0)
SMALL
(0.2)
MriMUM
'
* ** **
TROPHIC
STATUS
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
Oligotrophic
Eutrophic
Mesotrophic
OliRotrophic
Kutrophic
"'i '. •••;!»: :<•
' <
f ....•.-.„ ^ ,
*
BSS
day
118
134
0
132
0
0
133
0
0
100
100
101
105
106
106
106
106
124
101
101
101
104
HJ1
* ( i *
".A
*
PAST 1955-1979
ESS
day
269
241
0
244
0
0
240
0
0
293
290
268
262
256
241
241
233
210
312
313
312
295
295
'!'!'.'
?*•«
VI i
LSS
day
152
108
0
113
0
0
108
0
0
194
191
168
158
151
136
136
128
87
212
213
212
192
192
1KH
l')9
1 i
i i
SSR
—
0.89
0.12
0.63
0.19
1.00
1.00
1.00
1.00
1.00
0.99
0.86
0.87
0.99
1.00
1.00
1.00
1.00
1.00
1. 00
0.99
\ r ^^
MAXSD MINSD BSS
m
1.7
1.9
1.6
1.8
5.0
4.5
7.0
4.9
4.9
6.1
3.5
4.5
5.5
9.0
10.0
13.0
10.0
10.0
12.0
8.0
<) <••
m
0.1
0.2
0.1
0.1
0.2
0.2
0.2
0.4
0.4
0.4
0.4
1.0
1.0
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0 .1
day
68
85
0
76
116
0
85
116
0
68
69
70
69
69
70
70
72
73
69
70
70
70
. 71
72
72
71
7?
GlSS-2xC02
ESS
day
271
246.
ox
255
138
0
255
137
0
288
287
276
274
270
251
250
250
247
301
302
302
290
290
290
275
275
273
LSS
day
204
162
0
160
23
0
171
22
0
221
219
207
206
202
182
181
179
175
233
233
233
221
220
219
204
205
202
SSR
—
0.98
0.54
0.84
0.22
0.54
0.14
1.00
1.00
1.00
1.00
1.00
1.00
0.99
0.97
0.90
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
GISS
MAXSD MINSD BSS
m
1.7
2.1
1.9
1.3
1.3
0.9
5.3
6.7
9.1
5.9
4.0
5.3
2.9
4.3
5.3
5.8
6.7
8.6
7.7
8.6
10.6
11.5
13.4
9.6
m
0.1
0.1
0.1
0.4
0.1
0.3
0.2
0.1
0.1
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
0.4
day
-50
-49
-56
•
-48
-32
-31
-31
-36
-37
-36
-36
-34
-51
-32
-31
-31
-34
-33
-33
-34
-35
-34
- PAST
ESS
day
2
5
11
15
-6
-3
8
12
14
10
9
17
37
-11
-11
-10
-6
-5
-2
11
11
13
LSS
day
52
54
67
63
27
28
39
48
51
46
45
51 •
88
21
20
21
29
28
31
45
46
47
SSR
—
0.09
0.42
0.22
0.35
0.00
0.00
0.00
0.00
0.00
0.02
0.13
0.10
-0.10
0.00
0.00
0.00
0.00
0.00
0.00
0.01
0.00
0.00
-------
HSS Beginning seasonal stratification, i.e. first Julian day when difference
between surface and deep water temperature is greater than 1"C.
•".SS End seasonal stratification, i.e. last Julian day when difference
between surface and deep temperature is less than 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:
(1) Simulated epilimnetic temperatures were predominantly related to
weather and secondarily to lake morphometry. Weekly average epilimnetic
temperatures were raised by climate change for all lake classes. The
seasonally averaged water temperature rise was 3'C, compared to 4.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 epilimnetic temperatures
rose only about 2/C with climate change.
(2) Hypolimnetic temperatures were predominantly related to lake
morphometry and mixing events in spring, and only secondarily to weather in
summer. The highest temperatures were calculated for large, shallow,
eutrophic lakes. After climate change, hypolimnetic water temperatures were
as follows: shallow lakes, warmer by an average 3.1'C; deep lakes, cooler by
an average 1.1° C; small-area, medium depth lakes, cooler by 1.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 thr
season. Length of the stratification period was increased by 40 to 60 days.
(7) Climate change caused greater lake stability in spring and summer
In fall lakes were driven faster towards isothermal conditions.
102
-------
6. Summary
JLs a result of the research described here, a better understanding of
v. * freshwater inland lakes respond to variable atmospheric conditions has
'XT' trained.
Chapter 2 describes how a specific lake water temperature model was
r-cer-j-Iised to simulate the seasonal (spring to fall) temperature stratification
TCT a wide range of lake morphometries and meteorological conditions.
Mcc=L coefficients related to hypohmnetic eddy diffusivity, light attenuation,
••ice sheltering and convective heat transfer were generalized using theoretical
xr.d empirical model extensions. The proposed regional lake, water
ictnp-exature model simulates the onset of stratification, mixed layer depth,
ir.d •Crater temperatures well.
Hypolimnetic eddy diffusivity was estimated as a function of lake
surface area and stability frequency. Although the proposed relationship is a
simplification of the turbulent diffusion processes taking place in the
hypcKmnion, it was found to be useful in seasonal lake water temperature
modeling. Heat exchange between water and lake sediments, a process
commonly neglected in previous work, was found to be important for the
analysis of vertical hypolimnetic eddy diffusivity (Appendix A). Estimates of
hype-limnetic eddy diffusivity made without sedimentary heat flux were up to
one third smaller than those made with the heat flux. Effects of errors in
temperature measurements and sediment heat flux estimates on the estimated
vertical eddy diffusivity were evaluated as well.
Chapter 3 describes a first order analysis of uncertainty propagation in
lake temperature modeling. The output uncertainty is defined as the result
of deviations of the meteorological variables from their mean values. The
analysis was applied to systems with correlated and uncorrelated
meteorological variables. Surface water temperatures are strongly affected by
uncertain meteorological forcing. Air temperature and dew point temperature
fluctuations have a significant effect on lake temperature uncertainty.
Long-term average water temperature structure in lakes can be estimated by
computer model simulations for just one year when the results from the
statistical analysis of meteorological , variables are used as input. This
analysis presents a useful alternative for the study of long-term averages and
the variability of temperature structures in lakes due to variable
meteorological forcing. In addition, the analysis revealed the separate
contribution of each meteorological variable to water temperature uncertainty.
The analysis described in Chapter 4 was a first step in quantifying
potential thermal changes in inland lakes due to climate change. Rather than
using global climate change predictions, this analysis looked at the changes in
heat balance and. temperature profiles in a particularly warm year (1988)
103
-------
compared to a "normal" year (1971). A comparison was made for three
morphometrically different lakes located in north central US. Simulated
water temperatures were daily values driven by daily weather parameters and
verified against several sets of measurements. The results show that in the
warmer year, epilimnetic water temperatures were higher; evaporative water
loss increased; and summer stratification occurred earlier in the season.
Rather than analyzing particular years and particular lakes, emphasis in
Chapter 5 is on long term behavior and a wide range of lake morphometries
and trophic levels. The regional lake water temperature model was linked to
a daily meteorological data base to simulate daily water temperature profiles
over a period of twenty-five (1955-1979) years. Twenty seven classes of
lakes which are characteristic of the north-central US were investigated.
Output from a global climate model (GISS) was used to modify the weather
data base to account for the doubling of atmospheric C02- The simulations
predict that after climate change epilimnetic temperatures will be higher but
increase less than air temperature; hypolimnetic temperatures in seasonally
stratified dimictic lakes will be largely unchanged or even lower than at
present; evaporative water loss will be increased by as much as 300 mm for
the season, onset 67 stratification will occur earlier and. overturn later in the
season; and overall lake stability will become greater in spring and summer.
104
-------
References
Adams, E.E., Scott, S.A. and Ho, E. K. (1987). "Vertical diffusion in a
stratified cooling lake, Journal of Hydraulic Engineering, ASCE, Vol.
113(3), pp. 293-307.
Adams, E.E., D.J. Cosier, and K.R. Helfrich (1990). Evaporation from
heated water bodies: Predicting combined forced plus free convection,
Water Resources Research, Vol. 26(3), pp. 425^135.
Aldama, A.A, D.R.F. Harleman, and E.E. Adams (1989). Hypolimnetic
mixing in a weakly stratified lake, Water Resources Research, Vol.
25(5), pp. 1014-1024.
Andricevic, R. (1990). Cost-effective network -design for groundwater flow
monitoring, Stochastic Hydrology and Hydraulics, Springer International,
Berlin, Germany, 4(1), pp. 27-41.
Baker, D.G., Kuehriast, E.L., and Zandlo, J.A (1985). Climate of Minnesota,
Part 15 - Normal temperatures (1951-1980) and their application,
Agricultural Experimental Station University of Minnesota, 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
thennocline formation in lakes and reservoirs," Tech. Rep. 249, Ralph
M. Parsons Lab. for Water Resour. and Hydrodyn;, Dep. of Civ. Eng.,
Mass. Inst. of Technol.
Bolin, B. and B.R. Doos (1986). The Greenhouse Effect, Climate Change
and Ecosystems, John Wiley and Sons, New York, p. 541.
Canale, R. " P., and S. W. Effler (1989). Stochastic phosphorus model for
Onondaga Lake, Water Research., 23(8), 1009-1016.
105
-------
Carlson, R.E.(1977). 'A trophic state index for lakes', Limnology and
Oceanography, Vol.22, pp. 361-369.
Carslaw, H.S. and Jaeger, J.C. (1959). Conduction of heat in solids, Oxford
University Press, N.Y.
Chang, L.H., S.F. Railsback, and R.T. Brown (1992). Use of reservoir water
quality model to simulate global climate change effects on fish habitat,
Climatic Change, Vol. 20, pp. 277-296.
Colman, J.A., and Armstrong, D.E. (1987). Vertical eddy diffusivity
determined with 222Rn in the benthic boundary layer of ice-covered
lakes, Limnology and Oceanography, 32(3), pp. 577-590.
Coutant, C.C. (1990). Temperature-oxygen habitat for freshwater and coastal
striped bass in a changing climate, Transactions of the American
Fisheries Society, Vol.. 119(2), pp. 240-253.
Dake, J.M.K. and D.R.F. Harleman (1969). Thermal stratification in
lakes: Analytical and laboratory studies, Water Resources Research,
Vol. 5(2), pp. 404^95.
Dake, J.M.K. (1972). Evaporative cooling of a body of water, Water
Resources Research, Vol. 8(4), pp. 1087-1091.
Dettinger, M. D., and J. L. Wilson (1981). First order analysis of uncertainty
in numerical models of groundwater flow, 1, Mathematical development,
, Water Resour. Res., 17(1), 149-161.
Dhamotharan, S., (1979). A mathematical model for temperature and
turbidity stratification dynamics in shallow reservoirs, Ph.D. thesis,
University of Minnesota, p. 319.
ERLD/MNDNR: (1990). Minnesota Department of Natural Resources Fisheries
Division Lake Data Base, expanded by US Environmental Protection
Agency Environmental Research Laboratory Duluth.
Ellis, C. and E.G. Stefan (1991). Water temperature dynamics and heat
transfer beneath the ice cover of a lake, Limnology and Oceanography,
Vol. 36(2), pp. 324-335.
'' . '
Ford, D.E. (1976). Water temperature dynamics of dimictic lakes: Analysis
and prediction using integral energy concepts, Ph.D. thesis, University of
Minnesota, p. 432. -.• '
Ford, D.E. and E.G. Stefan (1980a). Thermal predictions using integral
energy model, Jour. Eydraulics Division, ASCE, Vol. 106, EY1, pp.
39-55.
Ford, D. E., and E. G. Stefan (1980b). Stratification variability in three
morphometrically different lakes under identical meteorological forcing,
Water Resour. Bull., AWRA, 16(2), 243-247.
106
-------
Gardner, R.H., R.V., O'Neill, J.B., Mankin, J.H., Carney (1981). A
comparison of sensitivity analysis and error analysis based on a stream
ecosystem model, Ecol Modeling, 12, 173-190.
Gargett, A.E. (1984). Vertical eddy diffusivity in the ocean interior,
Journal of Marine Research,42(2), pp. 359-393.
Gargett, A.E. and Holloway, G. (1984). Dissipation and diffusion by internal
wave breaking, Journal of Marine Research, 42(1), pp. 15-27.
Gerhard, H., J. Imberger, and M. Schimmele (1990). Vertical mixing in
Uberlinger See, western part of Lake Constance, Aquatic Sciences, Vol.
53(3), pp. 256-268.
Gorham, E., and Boyce, P.M. (1989). Influence of lake surface area and
depth upon thermal stratification and the depth of the summer
thermocline, J. Great Lakes Res., Vol. 15, pp. 233-245
Gu, R. and H. G. Stefan (1990). Year-round simulation of cold climate lakes,
Cold Regions Science and Technology 18, pp. 147-160.
Harbeck, G.E. Jr. (1962). A practical field technique for measuring reservoir
evaporation utilizing mass-transfer theory, U.S. Geol. Surv. Prof. Pap.,
272E, pp. 101-105.
Harleman, D.R.F., and K.A. Hurley (1976). Simulation of the vertical thermal
structure of lakes under transient meteorological conditions, Dynamics of
Stratification and of Stratified Flows in Large Lakes: Proceedings of
Workshop of the Committee on Lake Dynamics (Winsdor, Ontario:
International Joint Commission, Regional Office), pp. 79-96.
Heiskary, S.A., Wilson, C.B., and Larsen, D.P. (1987). Analysis of regional
patterns in lake water quality: Using ecoregions for lake management in
Minnesota, Lake and Reservoir Management, Vol 3, pp. 337-344.
Heiskary, S.A. and C.B. Wilson (1988). Minnesota lake water quality
assessment report, Minnesota Pollution Control Agency, St. Paul, p. 49.
Henderson-Sellers, B. (1988). Sensitivity of thermal stratification models
to changing boundary conditions, Applied Mathematical Modelling, Vol
12, pp. 31-43.
Horsch, G. M. and H. G. Stefan (1988). Cbnvective circulation littoral
water due to surface cooling, Limnology and Oceanography, Vol. 33(5),
pp. 1068-1083.
Eoughton, J. T. et al. editors (1990). Climatic Change: The IPCC Scientific
Assessment. Cambridge University Press.
Idso, S.B. and R. D. Jackson (1969). Thermal radiation from the
atmosphere, Jour, of Geophysical Research, 74(23), pp. 5397-5403.
107
-------
Idso, S.B. and R.G. Gilbert (1974). On the universality of the Poole and
Atkins. Secchi disk-light extinction equation, J. Appl. Ecol, Vol. 11, pp.
339-401.
Imberger, J. (1985). The diurnal mixed layer, Limnology and Oceanography,
Vol 30(4), pp. 737-770.
Imberger, J., and Patterson, J.C. (1989). Physical Limnology, in Advances
in Applied Mechanics, edited by J.W. Hutchinson and T.Y. Wu,
Academic Press, Vol 27, pp. 303-475.
Imboden, D.M., Lemrnin, U., Joller, T., and Schurter, M. (1983). Mixing
processes in lakes: mechanisms and ecological relevance, Schweiz. J.
Hydrol., 45(1), pp. 11-44.
Jassby, A. and Powell, T. "(1975). Vertical patterns of eddy diffusion
during stratification in Castle lake, California, Limnology and
Oceanography, 20(4), pp. 530-543.
Jones, P.D., Wigley, T.M.L., and Wright, P.B. (1986). Global temperature
variations between 1861 and 1984, Nature, Vol 332(31), pp. 430-434.
Keijman, J. Q. (1974). The estimation of energy balance of a lake from
simple weather data, Boundary-Layer Meteorology, Vol. 7, 399-407.
*
Kerr, R.A. (1989). 1988 ties for warmest year, Science, Vol 243(17), pp.
891-892.
Lewis, W.M. (1982). Vertical eddy diffusion in a large tropical lake,
Limnology ana Oceanography, Vol. 27(1), pp. 161-163.
Lewis, W. M., Jr., (1983). Temperature, heat, and mixing in Lake Valencia..
Venezuela, Limnology and Oceanography, Vol. 28(2), pp.273-286.
Likens, G.E., and Johnson, N.M. (1969). Measurement and analysis of th
annual heat budget for the sediments in two Wisconsin lakes. Limnolo?,
and Oceanography, Vol. 14(1), pp. 115-135.
Magnuson, J.J, J.B. Meisner, and D.K. Hill (1990). Potential changes ::
thermal habitat of Great Lakes fish after global climate warmim
Transactions of the American Fisheries Society, Vol. 119(2), p:
254-264.
McCormick, M. J. (1990). Potential changes in thermal structure and cycle :.-•
Lake Michigan due to global warming, Transactions of the Amer*^
Fisheries Society, Vol. 19(2), pp. 183-194.
McKinney, D. C. (1990). Predicting groundwater contamination: uncertain-
analysis and network design, Ph.D. thesis, Cornell University, 256 pp
McLaughlin, D. (1985). - A distributed parameter state space approach v
evaluating the accuracy of the groundwater predictions, Ph.D. th
-------
a. R.O., W.S. Combs, Jr., P.D. Smith, and A.S. Knoll (19791
Attenuation of light and daily integral rates of photosynthesis attained
ty planktonic algae, Limnology and Oceanography, vol. /4^bj, pp.
1338-1050.
r. J.D., J.L. Goddier, H.A. Regier, B.J. Shuter, and W.J. Christie
'1987). An assessment of the effects of climate warming on Great
I^akes Basin fishes, Journal of Great Lakes Research, Vol. 13(dJ, pp.
340-352.
l Research Council (1982). Carbon Dioxide/Climate Review Panel.
Carbon Dioxide and Climate: A Second Assessment. National
Academy Press, Washington, B.C., p. 72.
l Research Council (1983). Changing climate: report of the
Carbon Dioxide Assessment Committee. National Academy Fress,
Washington, D. C., p. 496.
sCAA. (1990). National. Center for Atmospheric Research - data support
section/Personal communication Roy Jenne and Dennis Joseph, Boulder,
Colorado.
SrSder. U.P., Schindler, P.W., Wirz, U.E. and Imboden, D.M. (1983).
Chemical and geochemical studies of Lake Biel H. A chemical approach
to lake mixing. Schweiz. Z. Hydrol, Vol. 45(1), pp. 45-61.
Os-ood, R.A. (1984). A 1984 study of the water quality °(T431???°{^tai1
" area lakes, Metropolitan Council, St. Paul, MN, Publ. No. 10-84-172, p.
40.
Qs^ood, R.A (1989). An evaluation of the effects of watershed treatment
~ systems on summertime phosphorus concentration in metropohtan area
lakes, Metropolitan Council, St. Paul, MN, Publ. No. 590-89-062, p. 85.
Osgood, R.A. (1990). Personal communication, Metropolitan Council, St.
Paul, Minnesota.
Patankar, S. V. (1988). Numerical heat transfer and fluid flow, McGraw-BjH,
New York, N.Y., 197 pp..
Poole, H.H. and W.R.G. Atkins (1929). Photo-electric measurements of
submarine illumination throughout the year, J. Mar. Biol. Assoc. U.K.,
Vol. 16, pp. 297-324.
Priscu, J.C., Spigel, R.H., Gibs, M.M., and Downes, M.T (1986). A
numerical analysis of hypoh'mnetic nitrogen . and phosphorus
transformations in Lake Rotoiti, New Zealand: A geothermally mfluenced
lake, Limnology and Oceanography, Vol. 31(4), pp. 812-8.J1.
Protopapas, A. L. (1988). Stochastic hydrologic analysis of soil-crop-dimate
interactions, Ph.D. thesis, Parsons Lab., Dep. of Civ. Eng., Mass.
Inst.of Technol., Cambridge, 239 pp.
109
-------
Protopapas, A.L., and R.L. Bras (1990). Uncertainty propagation with
numerical models for flow and transport in the unsaturated zone, Water
Resour. Res., 26 (10), 2463-2474, 1990.
Quay, P.D., W.S. Broecker, R.H. Hesslein, and D.W. Schindler (1980).
Vertical diffusion rates determined by tritium tracer experiments in the
thermocline and hypolimnion of two lakes, Limnology and. Oceanography,
Vol. 25(2), pp. 201-218.
Richardson, C. W. (1981). Stochastic simulation of daily precipitation,
temperature, and solar radiation, Water Resour. Res., 17(1), 182-190,
1981.
Riley, M.J. and H.G. Stefan (1987). Dynamic lake water quality simulation
model "Minlake", University of Minnesota, St. Anthony Falls Hydraulic
Laboratory, Report No. 263, p. 140.
Robertson, D.M. (1989). The use of lake water temperature and ice cover as
a climatic indicators, Ph.D. thesis, University of Wisconsin-Madison, p.
330.
Robertson, P.M., and R.A. Ragotzkie (1990). Changes in the thermal
structure of moderate to large size lakes in response to changes in air
temperature, Aguatic Sciences, Vol. 52(4), pp. 362-379.
Robinson, P.J., and Finkelstein, P.L. (1990) Strategies for the development
of climate scenarios for impact assessment: Phase I final report, USEPA
Atmospheric Research and Exposure Assessment Laboratory
EPA/600/53-90/026.
Sadhyram, Y., P. Vethamony, A. Suryanarayana, G. N. Swamy, and J. 5
Sastry (1988). Heat energy exchange over a large water body und*-:
stable atmospheric conditions, Boundary-Layer Meteorology, Vol. 44, ~,
171-180.
Scavia D., W. F. Powers, R. P. Canale, and J. L. Moody (1981). Compare =
of first order error analysis and Monte Carlo simulation -.-•
time-dependent lake eutrophication models, Water Resour. Res., 17u
1051-1059.
Schindler, D.W., K.G. Beaty, E.J. Fee, D.R. Cruikshank, E.R. DeBruyn, ! -
Findley, G.A. Londsey, J.A. Sherer, M. Stainton, M.A. Turner (IS*
Effects of climatic warming on lakes of the Central Boreal Fc-r*».
Science, Vol. 250(16), pp. 967-970.
Shapiro, J. and H. Pfannkuch (1973). The Minneapolis chain of lakn
. study of urban drainage and its effects 1971-1973, Limnological n—
Center, University of Minnesota, Minneapolis, Report No. 9.
Smith, C.R. and K.S. Baker (1981). Optical properties of the
natural waters, Applied Optics, Vol. 20(2), pp. 177-184.
110
-------
H H.. J. Imberger, and K.N. Rayner (1986). Modeling the diurnal
•^iicd layer, Limnology and Oceanography, 31, pp. 533-556.
; R.E., and Armstrong, D.E. (1983). Lake mixing and its relationship
ui epilimnetic phosphorus in Shagawa lake, Minnesota, Can. J. Fish.
.\
-------
Wdander, P. (1968). Theoretical forms for the vertical exchange coefficients in
a stratified fluid with application to lakes and seas, Acta Regiae
Societatis Scientiamm et Litterarum Gothoburgensis, Geophysica 1, pp.
1-26.
Wu, J. (1969). Wind stress and surface roughness at air-sea interface,
Jour. Geophysical Research, 74, 2, 444-455.
Ward, P.R.B. (1977). Diffusion in Lake Hypolimnia, Proceedings 17th
Congress, International Association for Hydraulic Research,
Baden-Baden, IAHR, pp. 103-110.
Willis, R. and W.W-G. Yeh (1987). Groundwater Systems Planning and
Management, Prentice-Hall,' Englewood Cliffs, NJ 07632, p. 415.
Winter, T. (1980). Personal communication, U.S. Geological Survey, Denver
Federal Center, Lakewood, CO 80225.
Wright, D.} R. Lawrenz, W. Popp, and M. Danks (1988). Acid precipitation
mitigation program, Progress report Thrush Lake Minnesota, Minnesota
Department of Natural Resources Division of Fish and Wildlife
Ecological Services Section, p. 177.
112
-------
Appendices
-------
APPENDIX A
Vertical diffusion in a small
stratified lake: Data and error analysis
Water temperature profiles were measured at 2 minute intervals in a
stratified temperate lake with a surface area of 0.06 km2 and a maximum
depth of 10 m from May 7 to August 9, 1989. The data were used to
calculate the vertical eddy diffusion coefficient (Kz) in the hypolimnion. The
depth is representative of a large number of lakes in the north central United
States. (KE) was calculated over time intervals of 1 to 15 days and varied
from 10"3 to 10"1 cm2s~l. A numerical model was developed for heat
conduction in the sediments, and heat flux between water and sediments was
incorporated into the relationship from which eddy diffusivity was estimated.
Heat flux between water and lake sediments, a term commonly neglected, was
found to be important in the Kz estimation. Kz values were related to
stratification stability as measured by the Brunt-Vaisala frequency N using
Welander's expression of the form Kz = o(N2)T Values of a were on the
order of 10'4 and 7 varied from - 0.36 to - 0.45 when Kz was given in
cm2s'1 and N is in s"1. An error analysis was conducted and the effects of
different choices of sampling intervals in time and depth on the eddy
diffusivity estimates were evaluated. The longest time interval (15 days) and
the smallest depth increment (1 m) used in this study were found to give the
best Kz estimation.
A.1 Introduction
Density stratification due to vertical temperature gradients inhibits
vertical mixing in lakes and reservoirs, and mixing in turn affects the
distribution of phytoplankton, nutrients, and other water quality constituents.
Quantifying turbulent transport phenomena is one of the major challenges in
lake and reservoir hydro thermal and water quality analysis. Specification of
vertical turbulent (eddy) diffusion coefficients in one-dimensional water quality
models, which are often used for decision-making, is particularly difficult.
In the analysis of vertical turbulent mixing by one-dimensional wind
energy models, the depth of the surface mixed (epilimnetic) layer is calculated
by an integral entrainment model while the vertical transport in the
hypolimnion is taken into account by a diffusion equation (Stefan and Ford
1975; Bloss and Harleman 1979; Ford and Stefan 1980). Although the
hypolimnion is isolated from the epilimnetic layer by the thermocline and its
associated density gradient, strong and sporadical local mixing events have
been observed in the hypolimnion (Jassby and Powell 1975; Imberger 1985;
Imberger and Patterson 1989). Such mixing events can originate from
oscillating boundary layers induced by seiche motions on the bottom of lakes,
A-l
-------
internal wave interaction and breakdown, shear instabilities in the thermocline
(billows), epilimnetic turbulent kinetic energy leakage to the hypolimnion and
double diffusion processes. Scales for such events range from the
Kolmorgorov scale to the lake basin scale. Eddy diffusion dependence on
stratification strength as measured by buoyancy frequency has been pointed
out consistently (Colman and Armstrong 1987; Gargett 1984; Gargett and
Holloway 1984; Imboden et al. 1983; Jassby and Powell 1975; Quay et al.
1980; Welander 1968).
Direct measurements of vertical turbulent diffusion in lakes are not easy
because of the 3-D nature of the diffusion field, and the spatial and temporal
scales. To estimate diffusion values, one can rely on measurements of water
temperatures or concentrations of natural tracers present in the water. Water
temperature measurement is the most commonly used method because of its
simplicity; however, a careful assessment of all external and internal heat
sources is required.
The purpose of this study was to estimate vertical eddy diffusion from
water temperature measurements in a typical inland shallow lake. Sediment
heat exchange, commonly neglected along with error analysis, is also included
in the estimation. Lastly, criteria for measurement intervals in space and
time that minimize the error in eddy diffusivity estimation are proposed.
The latter uses principles which are also used in groundwater monitoring
network design (Andricevic 1990).
A.2 Study Site
Ryan Lake, located in Minneapolis, Minnesota, has a surface area of
0.06 km2, mean depth of 5 meters and maximum depth of 10.5 m (Fig. A.I).
The lake, located in a flat terrain, suburban residential area, is highly
eutrophic, and regularly experiences winterkill of fish. The maximum depth
of 10 m is equal to the median maximum depth of 779 statistically analyzed
lakes in Minnesota. The depth of Ryan Lake can be considered as typical
for the north central United States.
Lake water temperatures were measured every two minutes at 1 m
intervals from the lake surface to the 10 m depth. Every 20 minutes, the
previous ten measurements were averaged and recorded. The measurement
scheme was adopted to reduce the "high frequency" electronic and
measurement noise, while retaining the fluctuations expected at timescales of
hours and days. The probes are rubber coated thermisters with a time
constant of 10 seconds. They were calibrated in a water bath prior to
installation. Absolute accuracy (values' measured by two adjacent probes at
known temperatures) was ± 0.05° C (95% confidence interval), while relative
accuracy (the difference between successive measurements by the same probe)
was O.OrC. A Campbell Scientific CR10 datalogger installed on a small raft
recorded the water temperatures. Hypolimnetic data for the period from May
7 to August 9, 1989, were selected for analysis because this period was
characterized by stable seasonal lake stratification. In 1990 measurements
were extended to sediment temperatures using probes identical to those in the
water. The sediments were soft, organic material and poorly consolidated, as
A-2
-------
10m max * I I
9m —
MEASUREMENT
LOCATION
Fig. A.I Ryan Lake bathymetry.
A-3
-------
indicated by the ease with which the thermister probe support rod was
installed.
A.3 Vertical Eddy Diffusivity
Many studies have assumed lake basins to be closed systems and have
estimated an average vertical eddy diffusion coefficient over the whole basin
(Adams et al. 1987; Gargett 1984; Imboden et al. 1983; Jassby and Powell
1975; Lewis 1983; Nyffeler et al. 1983; Priscu et al. 1986; Quay et al. 1980).
One of the methods for such estimation is through the budgets of . scalar
quantities such as temperature (Gargett 1984).
The one-dimensional, unsteady heat transport equation applied along the
vertical axis of a water column is:
The flux-gradient method for the computation of Kz reduces this equation to
the form
K, • [ f? ]"' { It / T(Qd< + «T|S - / Sdz } (A.2)
It is assumed that there is no vertical advection (w = 0) of water
anywhere. There is, however, a conductive heat flux H , from the sediment
into the water at the lake bottom (z = 0) and a radiation (penetrative) heat
flux H .. A heat balance for the water column between z = 0 and z
sol
therefore leads to a replacement for equation (A.2).
(A-3)
3t
Vertical kinematic thermal eddy diffusivity can be explicitly expressed as:
z z
/•
£_
dt
K2 = -£ 2 2 (A.4
Jz
where T(z,t) = measured water temperature distribution, z = upwar-
coordinate starting at the lake bottom, t = time, w = vertical component c:
velocity, and S = internal source term. At the time scale of 1 to 15 days ^
which Kz is computed, and without significant inflow or outflow to or iio—
the lake, net vertical velocity w is customarily small enough to be neglectec
Short term effects during storms and turnovers will show up implicitly in th-
value of Kz. The source term "S" in Eq. A.4 accounts for solar short wav-
A-4
-------
radiation absorbed in the water and heat flux through the water column to
or from the sediments at the bottom. For shallow inland lakes, the source
term can be particularly significant.
As pointed out by Gargett (1984), the budget method has two
advantages. First, few additional assumptions are needed to estimate Kz
from Equation A.4. Second, time averaging is implicit in the estimate of Kz.
With the exception of the surface mixed layer, turbulence in lakes occurs in
patches and intermittently. Turbulent "bursts" involve small volumes of
water (tens of cubic meters) and last on the order of minutes (Imberger
1985). Therefore, time averaging for such systems appears to be essential to
capture the long-term behavior.
Eddy diffusion dependence on buoyancy frequency was pointed out by
Welander (1968) and others. Welander derived an expression relating Kz to
the square of the Brunt-Vaisalla frequency (N) as Kz = a (N2)^, where N2
& p ^— density of water and g = acceleration of gravity. If
turbulence is generated by the dissipation of energy from large-scale motions,
7 = -1.0; otherwise, if it is generated by shear flow, 7 = -0.5. Welander's
analysis was very informative, but it was based on several assumptions:
steady state, no boundary effects, and a linear dependence of density on
temperature. Such assumptions are only marginally valid for lakes. The
results to be presented herein will show that Welander's theory fits lake data
reasonably well.
A.4 Sediment Heat Storage
Few previous analyses include heat flux to or from the sediments in the
eddy diffusivity estimation for summer conditions. A notable exception is
Stauffer and Armstrong's (1983) study of Shagawa Lake's western basin
(maximum depth 14 m). In principle, sediment heat flux is related to the
water temperature gradient at the sediment/water interface (Nyffeler 1983);
however, unknown turbulent heat transfer coefficients relating the flux to the
gradient as well as exceedingly small temperature gradients in the
near-sediment water limit the usefulness of the relationship. Relying on
measurements and computer simulations, Priscu et al. (1986) assumed that
the heat flux from the sediments to the water was constant. This was
physically justified for the geothermally influenced lake which they studied.
In the more general situation, conductive heat flux through the sediments is
variable in depth as well as with time (Birge et al. 1927; Likens and Johnson
1969).
In this study, a numerical model was developed to simulate sediment
heating or cooling by the overlying water. A one-dimensional, unsteady heat
conduction equation was applied since conduction into and out of the
sediments is essentially a 1-D process. The unsteady heat conduction
equation for the sediments is a simplified "version of Eq. A.I. Vertical
velocity, w, is zero because there is no advection and S = 0 because there
are no internal heat sources or sinks in the sediments. Kz in Eq. A.I is
replaced by Kzs = sediment thermal diffusivity, and T is replaced by Ts =
A-5
-------
sediment temperature. So 5Ts/dt = Kzs(52Ts/&2) js the heat conduction
equation applied to the sediments. The partial differential equation was
discretized using a control volume method (Patankar 1988) and solved by a
tridiagonal matrix algorithm. The boundary conditions are: (1) measured
water temperatures at the water/sediment interface and (2) no flux (adiabatic
boundary) at zj = 6 m depth below the sediment surface. It could be
shown by unsteady heat transfer analysis of a semi-infinite slab that seasonal
heat storage did not penetrate significantly beyond 6 m depth in an annual
cycle. Heat flux (Hsed) through the sediment/water interface is calculated as
the rate of change in sediment heat storage given by integration of computed
sediment temperature profiles T(z,t): -
Zd
Hsed = Ps CPS "If / Ts(2'1) dz (A-5>
0
where ps — bulk sediment density, cps is sediment specific heat and (ps cps)
is bulk specific heat of the sediments per unit volume.
Time series of measured sediment and overlying water temperatures are
given in Fig. A.2 down to depths of 1.5 m into the sediment. No probes
were placed at any greater depths. These temperatures were recorded from
April 3 to July 9, 1990. In the early part of this season, temperature
gradients, and hence heat fluxes into the sediments, are at a maximum.
High fluctuations of water and sediment surface temperatures correspond to
the spring overturn. Water temperatures at 0.5 m above the sediments and
at the sediment surface plotted in Fig. A.2 were almost identical, indicating
the presence of significant turbulent mixing in the water boundary layer.
The differences between 20 minute readings at the two elevations had an
average of 0.00845" C and a standard deviation of 0.088°C. The square of
the correlation coefficient (R2) was 0.98.
The unsteady sediment heat conduction model was calibrated for
thermal diffusivity. A sediment thermal diffusivity of 0.0022 cm2sec"1 applied
uniformly over depth (Gu and Stefan 1990) simulated measured sedimen*.
temperatures well (Fig. A.3). The maximum difference between calculated
and measured temperatures was 0.2"C in a range of 1.7 "C. Accuracy of the
measurements was 0.05'C. The maximum discrepancy was observed at 0.5 ~
below the sediment surface in April. For the rest of the period the difference
was less than 0.1 °C. Since sediment thermal properties do not change wi:::
season, the calibrated values should hold for year-round simulation.
Using the calibrated model and measured deep water temperatures frcr:
the 1989 data as the boundary condition at the water sediment interface, ti=
sediment temperatures for the period of lake eddy diffusivity studies (May
to August 9, 1989) could be simulated. Heat flux from the" sediments w^j
calculated using the simulated sediment temperatures in equation (2). Likerj
and Johnson (1969) used values of ps=l-2 g cm"3 and cps=0.8 cal g'^C"1 i"
soft bottom sediments, giving (ps CnS) = 0.96 cal cm'3 'C'1. A water v
solids ratio of 3:1 corresponded to the calibrated thermal diffusivity (Carsu»
and Jaeger 1959), and hence values of 0.8 cal g'^C"1 for specific heat and 1 •'
A-6
-------
7
6-
O
o
d>
L_
Z3
-t~'
D
v_
CD
QL
E
-------
calculated
measured
measured
measured
April 2 April 10 April 18 April 26 May 4 May 12
\lruhu-il and measured temperatures in lake sediments, 1990.
-------
accuracy
equation (2) is estimated to be 2 kcal nrMay1.
A.5 Water temperature observations
Measured hypolimnetic water temperatures based on 20 minute averages
are plotted in Fig. A.4. Depths are below the water surface which fluctuated
by less than 0.1 m. There were no spatial variations of water temperatures
to speak of, other than in the vertical direction. For the entire period of
record, stratification was stable. Above 6 m depth, water temperatures were
influenced by the deepening thennocline. 'The rise in water temperature at
the 4 m depth in July was due to epilimnetic wanning associated with
increased solar radiation and deepening of the mixed layer.
Daily time series of on site incident solar radiation, air temperature,
wind speed, 'and epilimnetic water temperature are given in Fig. A.5.
Following standard weather bureau procedure, solar radiation is a daily total,
wind speed is an average of three-hourly readings, and air temperature is the
mean of a daily maximum and minimum reading. The strong rise of suriace
water temperature from May 7 to 17 coincides with high solar radiation and
low wind. Water temperature fluctuations from May 17 to June 23 are the
result of high fluctuations in solar radiation and air temperatures. From
June 12 to the end of the observation period, high solar radiation and air
temperatures increased water surface temperatures.
The entire stratification dynamics are put in perspective in Fig. A.6.
The isotherms were developed from the water temperature records such as
plotted in Fig. A.4. The window of data analyzed for vertical eddy
diffusivities is shown in Fig. A.6.
A.6 Vertical eddy diffusion estimates
With temperature T(z,t) given by discrete measured values Ti,j at
increments Az and At in depth and time, respectively, Equation (A.4) nust
be discretized numerically to yield the eddy diffusion coefficient estimate: in
the form:
where i and N are the bottom and topmost layer of the lake, respectively,
ATt = water temperature difference over a time interval At at a fixed depth
z, Hsoi = solar radiation heat flux at depth z, Hsed = water/sediaent
interface heat flux, and AT = temperature difference over a vertical distance
Az averaged over the time interval At, p = density of water and c? =
specific heat of water.
A-9
-------
»—'
o
May 7 Moy ?l June A June Ul July 2 July 16 July .50 Aucjust tJ
Fig. A.4 Hypolimnetic lake water temperatures recorded at 2 min interval
in Ryan Lake, May 7 to August 9, 1989.
-------
800-
1 m depth water temperature
oir temperature _^
wind speed
_j ,_i j
•8.0
-6.0
-4.0
-2.0
-0.0
in
2
May 7 Mo/ 21 June 4 June 16 July 2 July 16 July JO August 13
Fig. A.5 Meteorological conditions (solar radiation, vnnd epetd and air
tcmpcratnrcm) during lh« period of in ventilation.
-------
to
E
o.
-------
Heat flux through the water-sediment interface was calculated by
Equation (A.5). The average heat flux was 7.0 kcal m2 day1, and the actual
time series is shown in Fig. A.7. The heat flux was from the water into the
sediment, i.e. the sediment acted as a heat sink throughout the period of
investigation (May 7 - August 9).
Internal solar radiation absorption was calculated for each depth from
measured radiation at the lake surface and an attenuation coefficient (Eq.
1.4). Bi-weekly Secchi depths varied from 0.8 m to 1.25 m during the
period of analysis (May 7 to August 9, 1989) with a mean of 1.0 m.
Relationship between total attenuation coefficient and Sechi disk depth
translates a Secchi depth of 1.0 m into an attenuation coefficient of 1.8 m*1
which was used throughout the analysis (Fig. 2.4). Solar radiation adsorbed
at the 4 m depth amounted to about one third of the sediment heat flux
(see Fig. A.7) and was less at greater depth. In general, if an internal heat
source due to solar radiation exists but is neglected, the eddy diffusion
coefficient will be overestimated. Although not shown in Fig. A.7, solar
radiation and water to sediment heat flux had different signs in Equation
(A.6) because absorbed solar radiation is an input to the water and heat flux
to the sediments is a loss from the water during the period of study.
Vertical eddy diffusion coefficients calculated for sampling intervals of
five days and depth increments of 1 m (Fig. A.8) show decreasing values
with time as seasonal stratification progresses. High variability in space and
time is apparent. Vertical eddy diffusivity coefficient values ranged from
approximately 0.001 to 0.1 cmVi with an average near 0.01 cmV1. The
highest eddy diffusivity was found at the greatest depth (near the lake
bottom) while the 4 m depth had the lowest values. Eddy diffusivity near
the lake sediments is produced mainly by the interaction of internal waves
with the lake bottom resulting in internal breaking waves and by turbulence
induced by bottom shear during internal seiche motion. All of these are
contributing to an intensely mixed lake bottom boundary layer, as previously
shown by the temperature records in Fig. A.2. One result of this mixing is
the decrease in stratification intensity with greater depth in the hypolimnion.
There is also a positive feedback because shear-induced turbulence is
dampened by density stratification conditions.
Eddy diffusion coefficient estimates versus static density stability (N2)
for different sampling periods with and without consideration of the
sedimentary heat source term are given in Fig. A.9. A least squares linear
regression was used to estimate coefficients a and 7 in the relationship Kz
= Q(N2)7 . As expected, an inverse relationship between Kz and N2 was
observea. When the sediment heat flux term was omitted, eddy diffusivity
was underestimated. The error was up to one third of the estimated values.
It is noteworthy that a stronger dependence of Kz on N2 was observed when
the sedimentary heat source term was considered (Fig. A.9).
' A-13
-------
I-1
£>.
Q
fN
CJ
X
Z)
_J
U_
Id
n:
water-sediment heat flux
solar radiation heat flux 4m depth
water temperature above sediment y
water temperature above sediment ../""
F-5.4
r5.3
r5.2
1-5.1
r5,0
r4.9
I-4.8
j-4.7
j-4.6
r4.5
r4.4
r4.3
r4.2
-4.1
O
o
UJ
C£
ID
h-
<
cn
Ld
Q_
^
LJ
4.0
May 7 May 21 June 4 June 18 July 2 July 16 July 30 August 13
Fig. A.7 Heat flux through the sediment-water interface and solar
shortwave radiation received at the 4m depth, May 7 to August
-------
10°
RYAN LAKE
TURBULENT DIFFUSION COEFFICIENT TIME SERIES
i—i
en
I
u
O)
"E
u
ID"
IO'J
10'J
10-
Apr 30
7 M
6 M
••• 5 M
4 M DEPTH
i i [ i 'i i i i i i i i j i i r ' i r i i i i [ i i i i i i i i
May 20 Jun 9 Jun 29 Jul 19 Aug 8
Fig. A.8 Vertical turbulent diffusion coefficient time series in Ryan Lake,
-------
o
0>
00
CM
O
o
4>
CO
u
«
R'=0.32 . T=1 doy •
SEDIMENT HEAT NOT
K =1.15X10""(N2)-041
RJ=0.67 . 7=5 doys* ,
SEDIMENT HEAT NOT CONSIDERED
R=0.87 , T=10 doys
SEDIMENT HEAT NOT CONSIDERED
rrj i .-..
-------
A.I Regression coefficients for Kz =
Sampling
Interval
• (days)
1
5
10
• 15
Hsed * 0
a 7
9.14*10-5
9.77*10-5
1.12*10-<
1.37*10-<
-0.45*0.017
-0.45*0.022
-0.44*0.018
-0.43*0.018
Hsed = 0
a 7
2.20*10-5
1.15*10"*
1.47*10-4
1.87*10'*
-0.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
i.n estimate of an instantaneous value. Variability of the eddy diffusivity was
•.he highest for a Campling interval of one day. Different regression lines
could be fitted to the one day results without changing the regression
coefficient R2. By increasing the sampling interval, the effects of variable
meteorological conditions and mixing events were averaged out over longer
and longer periods. With a longer sampling interval, the linear regression fit
was better. Regression coefficients with standard errors are summarized in
Table A.I.
A. 7 Error Analysis
Uncertainties in the estimated Kz values are introduced by water
temperature measurement errors, model parameter values, and model
formulation. Magnitudes of errors for key variables in Equation (A.3) are
listed in Table A.2. The first three of these errors are measurement errors.
The last value is based on uncertainty in estimates of specific heat and soil
temperature measurements. 2.0 kcalm*2 day1 is about 30% of the average
net heat flux value of 7 kcal nr2dayl.
Table A.2 Errors
Variable
Symbol
Error
Temporal temperature
differential ATt (measured)
Depth temperature
differential ATZ (measured)
Depth increment Az (measured)
Source heat flux AH (calculated)
Az
c
s
0.01 (°C)
0.05 (-C.)
0.01 (m)
2.0 (kcalnr'day1)
A-17
-------
To assess the estimation error of eddy diffusivity, a small perturbation e
of the variables in Table A. 2 is introduced into equation (A. 6). Temporal
and depth temperature differentials, depth increments, and source heat fluxes
are considered random variables and represented by the mean value plus a
perturbation term. The perturbation terms have zero mean, and standard
deviations equal to one-half the values given in Table A. 2. Mean (denoted
by overbar) and perturbation .(denoted by e) were expressed as:
Kz = K7 + ekz ' (A.7)
ATt = ATT + et (A.8)
Az = Az + e (A.9)
/. ' S = 5 + eg - (A.ll)
with the statistical properties
E(e) = 0 (mean of e) (A. 12)
E(f2) =
-------
VAR(KZ)=
A2iA2j£
e2e2
£t£Az
4Az2
AT2
AT
4Az2
AT2
z
1 N
•i- S
At J = 1
AZ
TjAZj
4 AZ
Az
AZ<
AT
4Az
Az
N
E
At
16
N
AtAz3 ia,
AZ
(A.15)
where i or j designate the depth under consideration, Az = depth increment,
ATj or ATj = temporal temperature difference at the depth i or j,
respectively. N = total number of measuring points below the flux surface
under consideration, ATZ = temperature difference over the depth increment
Az averaged over the time interval At, S is source term. Other variables are
given in the main text (Table A.2).
Vertical profiles of the mean eddy diffusion coefficient plus or minus
two standard deviation intervals are given in Fig. A. 10. These profiles are
for five day sampling intervals and depth increments of 1 m. The error in
Kz estimation increased with depth mainly due to the smaller temperature
gradient near the lake bottom.
The dependence of calculated KZ values on the frequency and spacing of
the water temperature measurements is illustrated in Figs. A.11 and A.12. If
sampling intervals exceeded four days, the error in estimated Kz values did
not change significantly, regardless of depth. Sampling intervals of three days
or less increased the error.
The tradeoff between depth and time intervals with regard to the error
in Kz estimation is illustrated in Fig. A.12 by isolines of equal 1a\z values.
Three regions can be distinguished on the graph: (1) up to a sampling
interval of three days the error was a function of the time increment only.
The bigger errors correspond to the smaller sampling time increments. (2)
From 5 to 10 days sampling interval, errors were a function of both Az and
At. That was the region where trade-off between space and time was possible
in order to obtain the same estimation error. In general, errors were.
A-19
-------
3-
4-
5-
1 6-
x
a 7-
u
o
U
9J
10
Jun 9, 1989
-2(T
0.000 0.010 0^620 0.030 0.040 0.050 'o.060
4
5
1 6
x
h—
Q. 7.
UJ
O
8
9
10
K^cm'sec"1)
Jun 29, 1989
I '
mean
+ 2
-------
O.O70
0.016-
£'
b
(N
0.012-
0.008-
0.004-
0.000
i ' i ' I ' l ' l
m depth
5 m
Storting date: Jul 19, 19Q9 6
7 m
8 m
1 2 3 4 5 6 7 8 9101112131415
TIME INCREMENT (DAYS)
Fig. A.ll Estimated eddy diffusion errors (2crK ) for different sampling
intervals.
-------
£
to
3.00
-CN 2.60
2.20
Z
UJ
2
LJ
cr
o
1.80 -]
^-
o_
UJ
Q
.40
1.00
88888888 88
1 2 3 4 5 6 7 8 9 10 11.12 13 14 15
TIME INCREMENT (DAYS)
Fig. A.12 Estimated eddy diffusion errors 2aK (cmV1) space-time tradeoff,
7 m depth Jul 19.
-------
decreasing for smaller Az and larger At. (3) For sampling intervals larger
than 10 days the estimation error became primarily a function of Az, i.e. the
more measuring points used in a profile, the smaller the error in Kz.
A. 7 Conclusions
The vertical turbulent eddy diffusion coefficients in the hypolimnion of a
temperature stratified temperate lake with a depth typical of the north
central United States were evaluated from water temperatures measured at 2
minute intervals from May 7 to August 9, 1989. Kz values typically
increased by a factor of 10 between 4 m depth and 9 m depth. Eddy
diffusion coefficients Kz were computed for periods from 1 to 15 days and
varied from 0.001 to 0.1 cm^/s for the 1-day intervals and from 0.002 to 0.04
cm2/s for 15-day intervals. Kz also varied with stratification stability
measured by the Brunt-Vaisala frequency N. The relationship Kz =
produced the best fit when a = 0.00014 and 7 = - 0.43 for periods of 15
days, where Kz is in cmV1 and N in s~l . As the time step was shortened
to one day, the fit became poorer and 7 values changed slightly (Fig. A. 9).
The value 7 = - 0.43 fits Welander's (1969) model for shear driven
thermocline erosion. The a value is related to lake size (Fig. 2.2).
The water temperatures measured and recorded every 20 minutes at the
sediment /water interface and at 0.5 m above showed a mean difference of
only 0.008' C and nearly identical responses in time (Fig. A.2) over a period
of three months. This is indicative of a well-mixed boundary layer near the
lake bed.
Heat exchange between water and lake sediments was found to be
important to the analysis of vertical thermal diffusivity. A numerical model
was used to estimate sedimentary heat flux for incorporation into the eddy
diffusivity estimation. A mean value of sedimentary heat fluxes during the
summer period (May - August) was 7 kcal nr2dayl. Estimates of Kz made
without sedimentary heat flux were up to one third smaller than those made
with that heat flux. Heat input by solar radiation through the water surface
also influences the estimates of Kz, but this influence diminishes with depth.
Effects of errors in temperature measurements and sediment heat flux
estimations on calculated vertical eddy diffusion coefficients were evaluated.
Estimation errors were much larger near the lake bottom (in the
hypolimnion) than in the thermocline region (Fig. A. 10). The smallest
estimation errors of the eddy diffusivity were obtained for sampling intervals
of 15 days and depth increments of 1.0 m. At the 7 m depth, the error was
about 0.0011 cm2/s relative to a value on the order of 0.010 or 11 percent
(see Figs. A.8 and A. 11). The error doubled when the depth increment was
trippled to 3.0 m or when the sampling interval was reduced from 15 days to
5 days (Fig. A. 12). At the shorter sampling interval the error was, however,
insensitive to the depth increment. When the sampling interval was reduced
to 1.5 days, the estimation error increased to 0.008 cm2/s or nearly 80% of
the estimated value calculated at the 7 m depth. Thus the recommendation
is to sample at longer time intervals (several days) and at finer depth
resolution (order of 1 m).
A-23
-------
APPENDIX B
Temperature equation discretization
Temperature equation (l.ll is discretized using the control volume
method. For intermediate control volumes (i = 2, m-1).
A . k k-fl
~
Ai
where At is time step, Az is control volume width (constant) , k stands for
time, i stands for control volume location. Source term H, is described by
equation 1.4. System matrix of the deterministic model Ac(k) contains terms
on the left hand side of equation (B.I).
For the surface control volume (i=l) equation (B.I) will differ: eddy
diffusivity K 1-0.5 is zero, the first entry in the matrix is term multiplied by
T$+1, source terms are equations (1.2), (1.3), (1.6), (1.7), (1.8) and (1.9). For
the bottom control volume (i=m) insulation is imposed by setting Ki+o-s
equal zero. Diagonal entry in the matrix is term multiplied by Tik+1.
Eddy diffusivities at the control volume interfaces are defined as
harmonic mean .
2 K.K. .
Ki-o.5 = - 1-^~ (B-2)
A-24
-------
APPENDIX C
Temperature equation linearization
Matrix Ac(k) has the same entries as the system matrix Ac(k) given in
Appendix B. Matrix B(k) is tridiagonal, and contains derivatives with
respect to lake water temperature at time step k. For the intermediate
control volumes :
A- «*r ^--n^ -k+i 5K- nK -k+i
i—0.5 i—O.o rp _ i—Q.5 rp
" I ~" XY. ^- : i " i~ *• :
A.
A...!- dK. ... .k-t-l <9K. -k+l -, k
1—0.5 I 1-0.5 rp 1-0.5 rp T" -1.
1 . . 1 . . — » , i. . \ i . T
i i i
(A. ., „ _ r l/IV. . .k+1 UJV. , _ - .K
1 + 0.5 1+0.5 rp 1+0.5 rp
A I ffTk i+1 dTk '
A L 5T "*" 5T J At
i i i
A. -cr 5K. n, .k+1 5K..,.k+l
1 + 0.5] 1+0.5 ,p _ ^ 1—0.5 rr, I rp.
^ L ark i+l pr-k » I i+l
where ft = i if i = l else ft = 0
For the surface control volume matrix B(k) will slightly differ. First,
terms which are multiplied by 'Tj are the first entry (bn). Secondly, eddy
diffusivity K j-o-s is equal to zero. Thirdly, additional terms which take into
account boundary conditions have to be added (Equations 1.6, 1.8, and 1.9).
These equations have to be differentiated and evaluated with respect to the
water temperature in the first control volume.
k k k
•r dHbr 5HC dHe i Az k
Had = - + - + — -r - T (c.i)
For the bottom control volume matrix B(k) will also slightly differ.
Diagonal term (bmm) is the one which is multiplied by Tj. Eddy diffusivity
Kio-s is zero.
A-25
-------
Eddy diffusivity vector K(z,k) contains epilimnion and hypolimnion
diffusivities. Epilimnetic diffusivities are a function of wind speed, thus their
partial derivative with respect to lake temperature is zero entry. This is not
the case for the hypolimnetic eddy diffusivity.
i 0P(T) g
K = a (P)° where N2 is Brunt-Vaisala frequency N' =( ) -^
dz p
dK . .. K2 dK. 8N.
LM = 2 I i —-1 (C.2)
5T. (K._1 + K.)2 5N. 0T.
A
Matrix F(k) contains terms which require evaluation of first order
derivatives with respect to uncertain meteorological variables. Equations 1.3,
1.4, 1.6, 1.8, and 1.9 have to be differentiated and evaluated at the nominal
values of those variables. Entries in the matrix are grouped with respect to
the perturbed meteorological parameters. Matrix dimensions are m x m+3
[ F1(kj;
^
where air perturbation temperature vector F (k) is:
Az k
Hcl = (+ ' ~T/a (c-3)
dew point temperature perturbation vector F (k) is:
» •
5He 0ea Az k
Hc2 = ( TT 7^ ) —— T'di
wind speed perturbation (mxm) matrix F,(k) is:
..r 5K ._. .k+l 5K . „.. .k+i -, k
21| L±i T i-Q.5 T | w ,
L flw.k '-1 flw> Xi J Wsi-l
Ai-^.
Ai
Ai.o.5
A.
5K -_^^ -k+1 3K •_/>« -k+! 1 k
i-Q.5 rp i—Q.5 IT. W ' -4-
n» » J> • A ~LI ** "^. 1 • I * * S * T
0Wk 1~1 5W ' J 1
A-26
-------
&HCS ) ws;
-T,., - . ruj T; | + &HC, I ws; +
A.
1-4S -k+1 ^ k
A.
where
T
w /
» J si
Az k
HC3 = ( — + — ) - Ws (C.5)
First an4 last control volume have interface diffusivities K. and
K. equal to zero, and the additional term HC3 is present only in the first
l"~^J.5
control volume.
A
Solar radiation perturbation vector F4(k) is
dHsn Az k
HC4 = ( — j- ) - H's. (C.6)
9H wC
A-27
-------
APPENDIX D
Cross-*enn evaluations
Air and dew point temperature have significant correlation. The
covariance matrix between successive days for these two parameters has been
evaluated according to Protopapas (1988):
Cov [ C'(n) C'(k) ] = S(n) Mc S(k)T (D.I)
where
S(n) =
<7taln)
0
0
0
0
<7td(n
0
0
0
) °
0
0
0
0
0
0
Aa(n,k) Ptatd(n,k) 0 0
/>td(n,k) 0 0
LC —
00 00
00 00
S(k) =
0u(k)
0
0
0
0
tftd(k
0
0
0
) °
0
0
0
0
0
0
where Aa is air temperature standard deviation, k) is air temperature correlation between two day
Ptd(n>k) is dew point temperature correlation, ptatd(n,k) is air temperature
dew point temperature correlation.
Replacing index n by index k leads to the form of the covari;u:<-
matrix M(k,k) for the zero time lag. In addition, the diagonal terms an- ^
equal to one in the correlation matrix Mc. They are equal to the standt-
deviations of air temperature ffat(k), dew point temperature ortd(k), *--
speed o"ws(k), and solar radiation <7sr(k) in the matrices S(k), and S(n).
A-28
-------
If meteorological perturbations are not cross-correlated covariance matrix
Muc(k,k) is
-------
APPENDIX E
Regional lake water temperature simulation model
A-30
-------
LAKE INPUT DATA FILE
TITLE
SEKI
NDAYS
NPRNT
DZLL DZUL BETA EMISS WCOEF WSTR
WCHANL WLAKE DBL ST S FT
ELCB ALPHA BW
Initial conditions^
Secchi depth reading
Number of days for output
Dates for output
Heat budget and mixing coefficient
Initial stage & Outflow channel
MBOT NM NPRINT INFLOW DY MONTH ISTART MYEAR
Z(1)...Z(MBOT)
T(1)...T(MBOT)
Field data section
NFNPRFLE'
NFLD
DEPTH(1)...DEPTH(NF)
FD ATA( 1 )...FD ATA(NF)
Number of depths & parameters
Parameter code (1)
Depths
Temperatures
NDAYS
A-31
-------
EXAMPLE INPUT DATA. FILE
LAKE CALHOUN 1971 ** from APRIL through december **
2
9
520 608 629 719 802 824 913 1011 1028
0.15 1.00 .4 .97 24.5 0.4
100. 100. 200 224.0 .001 .035
205 18 100
24 8 1 0 91 4 1 1971
0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5
12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5
4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
4, 4. 4. 4.
11 1
1
0. 2. 4. 6. 7. 8. 10. 12T15. 20. 23.
13.1 12.8 12.5 12.5 12.5 10.4 8.6 7.9 7.4 7. 6.9
14 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12. 15. 20. 25.
20.4 20.3 20.1 19.8 15.1 13.6 113 11.7 11. 9.2 8.1 7.3 7.1 6.9
12 1
1
0. 1. 2. 3. 4. 5. 6. 8. 10. 12. 15. 20.
24.4 24.6 24.6 24.3 23.5 19. 15.1 11.5 10. 8.8 7.2 7.
141
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12. 15. 20. 24.
23.1 23.1 23.1 23. 22.9 22.8 17.6 12.4 11.4 10.2 9. 8.2 7.8 7.7
14 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12. 15. 20. 24.
21. 21. 20.8 20.8 20.8 20.6 19.4 14.7 11.3 9.8 9. 8.2 7.5 7.4
11 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 10. 12. 15. 20.
22.7 22.7 22.7 22.7 22.6 22.6 20.5 15.7 12. 10.2 9.2 8.1 7.7
15 1
1
0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 12. 15. 20. 25.
21.7 21.7 21.5 21.5 21.5 21.3 21.1 18.5 13.3 10.6 10. 9. 7:9 7.6 7.5
13 1
1
0. 1. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20. 22.
14.5 14.5 14.3 14.3 14.2 14.2 14.2 11.1 10.1 9.3 8.6 8.4 8.3
71 '
1
0. 12. 13. 14. 15. 20. 23.
117 117 10.5 9.0 8.4 7.9 7.8
A-32
-------
LAKE INPUT DATA VARIABLES
SEKI
NDAYS
NPRNT
DZLL
DZUL
BETA
EMISS
WCOEF
WSTR
WCHANL
WLAKE
DEL
ST
S
FT
ELCB
ALPHA
BW
MBOT
NM
NPRINT
INFLOW
DY
MONTH
ISTART
MYEAR
Z
T
NF
NPRFLE
NFLD
DEPTH
FDATA
Secchi depth reading (m)
Number of particular dates selected for output
Dates selected for output
Minimum layer thickness (m)
Maximum layer thickness (m)
Surface absorption factor for solar radiation
Eraissivity of water
Wind coefficient for convective heat loss
Wind sheltering coefficient
Width of the inlet channel (m)
Width of lake (m)
Elevation of the bed in the bottom layer (m)
Stage of the lake on the first day of simulation (m)
Bed slope at inlet channel
Roughness coefficient at inlet channel
Elevation of the bottom of the outlet channel (m)
Side slope of outlet channel
Bottom width of outlet channel (m)
Total number of layers in the initial conditions
Number of months to be simulated
Interval between days for tabular output
Number of inflow and outflow sources
Julian day of first day of simulation
First month of simulation
Day of the month that simulation starts
Year of simulation
Depths in the initial conditions (m)
Temperatures in the initial conditions (°C)
Number of depths in field data profile
Number of profiles in the field data
Parameter code to match field data profiles to state variables
1 = Temperature (°C)
Field data depths (m)
Field data temperatures (°C)
A-33
-------
METEOROLOGICAL DATA FILE
MONTH KDAYS MYEAR
TAIR TDEW WIND DRCT RAD CR PR
KDAYS
where
MONTH
KDAYS
MYEAR
TAIR
TDEW
WIND
DRCT
CR
PR
Month
Total number of days in the month
Year
Average air temperature (oF)
Dew point temperature (oF)
Average wind speed (mph)
Wind direction
Percent of sunshine
Precipitation (inches)
EXAMPLE METEOROLOGICAL DATA FILE
4
32
21.5
26
32.5
34
39.5
53
55
49
54
54
50
39
40
55
61
55
52.5
30
22
10
14
12
14
23
29
32
27
28
34
25
20
23
32
46
36
37
1971
22.4
19.3
13.9
7.6
3.9
6.7
9.3
13.5
13
15.3
13.5
6.8
9.7
4.3
13.9
12.9
6.6
11.7
320
320
310
280
310
30
80
50
290
80
340
280
280
50
100
40
190
140
203
324
514
568
558
489
533
398
553
527
432
373
389
471
531
347
565
76
9
8
93
100
100
100
100
94
100
100
62
- 78
73
79
99
79
97
7
13
6
0
0
0
0
0
0
3
0
0
0
1
0
0
2
0
11
KDAYS
A-34
-------
1U4W789
SLARGE
PROGRAM REGIONAL
C
C THIS PROGRAM IS MODIFIED VERSION OF THE MINLAKE
C MODEL (JULEY * STEFAN, 1987. UNIVERSITY OF MINNESOTA
C SAFHL, PROJECT REPORT * U3\ THE COMPUTER CODE
C IS ADDAPTED FOR THE REGIONAL LAKE WATER TEMPERATURE
C SIMULATION IN A LAKE ATTATCH USER SUBROUTINE DURING
C UNKING.
C
COMMON/MTHD/TAIR(31).TDEW(31).RAD(31),CR(31).WIND(31).
+ PR(31).DRCT(3I)
CX)MMON/RESULT.T2('W),CHLATOT(«)).PA2(«),PTSUM(W).BODy40).
COMMON/FLOW/HMK(41).OE(40),FVCHLA{5),PE(5.41)
COMMON/VObJZMAX.DZ(-«).Z<40)A(«l).V(W),TV(W)J\TOP(41)J)BL
COMMON/SURSDZ{«>).SZ<«>),LAY(«0).AVGt(<.60).SVOL(«0)
COMMON/CHOICE-MODEL.NrrRO.IPRNT(6XNDAYS_STRNT(30VNCLASS.PLT(30)
COMMON/WATER.BETA.EMISSJCKljaClHKMAX.WCOEF.'WSTR
COMVON/CHANE^'WCHANUELCaALPHA.BW.WLAKE
. COMMON/STEPS:DZU_DZUUMBOT.NM.NPRINT.MDAY.MONTKILAY.DY
OOMMON/STAT,SUMXY(10).SUMX(10)5UMY(10)J{SQ(10).
+YSOX10).RSOC10).RMS(10).RELM(10),MTHREL(10).MDAYREUlO).ZREmO).
* ZRELM( 10).RS( 10XREU lOJ.MTH RMS( 10).MDA YRMS( 10).ZRS( 10).2RMS( 10)
COMMONnNFi.OU,OIN(5).TlN(S),PAlN(S).BODIN(5).DOIN(5).CiN;S).
+CDIN(5)_XNH!S(5;.KNOIS(5).CHLA1N(3J)
COMMON«OURC£RM(3.JO).PROD(40).XMR(3.*)),PRODSUM(-!0)
COMMON/FIELD/ !FLAG(10).FLDATA(ia50).DEPTHC50),NFLD{10)
COMMON/FILE' DIN.MET.FLO.TAPE&TAPEUREC
COMMONmn.' TTTLE
DIMENSION B(10)5TATVAR(60VICX(-))
CHARACTER'04 17TLE
CHARACTER-16 DlN.METJLO.TAPE&TAPEl
CHARACrER'I TS
-------
505 FORMATC REQUEST CHANGE OF TITLE (Y/N) T.SXJ)
506 FORMATC ENTER NEW TITLE .-•)
WRTrec.505)
READC.XA)1) XS
IF(XSEQ.-r .OR. XS.EQ.V) THEN
WR1TEC.SOIS)
READC.-(A)-) TTTLE
ENDIF '
WRJTE(8,999) FF
WRrTE(8,1900)
999 FORMAT{LXA1)
1001 FORMATC ENTER INPUT DATA FILE NAME :'.10X/)
lOO: FORMATC ENTER RLE SAME FOR TABULAR OUTPUT: '-0
1000 FORMATC ENTER METEOROLOGICAL DATA FILE NAME : Vj
1003 FORMATf ENTER INFLOW DATA FILE NAME :'.9XJ>
1900 FORMATC//.8XO
C
CALLSTART(STAFT.lSTART.INFLO»'.MYEAR.IRUN.ILSEn)
C
C"'" Call LAKE routine u> change or >dJ input quanliba •••
ZMAX=ST-DBL
CALL SETUP
11 IF(DZ(I).GT.D2UI_AND.MBOT.LT.W! THEN
CALL SPUT(ULAY)
GO TO 11
ENDIF
l-I+l
IF(LGT.MBOT.OR.I.GT.40) GOTO 12
GOTOn
12 CALL SET2(MBOT)
CALL VOLUME(MBOT)
CALL AREA
CALL ABOUND
CALLTVOL(MBOT)
C_DETERMINATION OF INITIAL MLXED LA'i'ER DEPTH
1LAY-1
DO7I-1.MBOT-1
IF(T2(I).GT.T2(I + 1)+CONST) GO TO S
1 ILAY»ILAY»1
8 DMLX=Z(lLA^-)+DZ(ILAY)-aS
DO 9 1=1.10
R£LM(l)«a
9 RMS(I)=aO
MP.O
IPRNT(1)»1
KDAYS=0
IPRNT(5)=IPRXT(5)-1
MDAY=0
CALLPTABLE
1PRNT(5) = IPRNT(5) * 1
C~ Sun simubiion for each montb
EDA Y= 345.
YEAR«REAL(MYE.\R)
IF(AMOD(YEAR.4.0).EO.aO) EDAY-.Vid
ETSUM=a
HT5UM-a
C
DO100J»LNM
CALL MTHDATA(MONTH.KDAYSJrfYEAR)
EDAY-345.
YEAR=REAUMYEAR)
IF(AMOD(YEAR.4.0).EO.a) EDAY«3oi
C-START SIMULATION FOR EACH DAY
DO 200 MDAY-ISTART.KDAYS
NFLOWolNFLOW
CALL LAKE
-------
c
C-DETERMINATION OF WIND MIXING ORDER
C_HEAT IS ABSORBED FIRST, THEN WATER COLUMN IS MIXED BY THE WIND
C
CALLHEBUG(ILAY.TMIXQNET,HS.HA,HBR,KE.HC.
+TAIR(MDAY),TDEW(MDAY),CR(MDAY).RAD(MDAY).WIND(MDAY),
•HPAN.PCOEFF.SEKI)
CALL CONMlX(ILAY.TMrXMBOT)
C-.CALCULATION OF EVAP. IN TERMS OF VOLUME
C-CALCULATES LATENT HEAT OF VAPORIZATION ALV
HED»HE/((597Jl*5«31'n(l))'RHOCn(l),C2(l),CD2(l)))
HEV»HED*ATOP(l)
CALL WINM1X(RKE,TMDULAY.MBOT)
DMIX«Z(lLAY)».EQ.NPRNT(NDAYS))THEN
C — Access and output field data
CALL FDATA(NF)
WRTTEC.3020)
READ(*,99) XS
99 FORMAT(A1)
3020 FOR-MATC///: DO YOU WANT TO MEW GRAPHICAL RESULTS: ( Y/N)'.
+ 3»)
C_.Ca[l plouing routines
IF(XS.EQ/Y' .OR. XS.EQ.VT THEN
CALL SUBPLOT(NF.MYEAR)
ENDIF
ENDIF
IPRNT{1)=0
200 CONTINUE
ISTART=1
100 CONTINUE
C_CotBputc 2nd list statistics!
C_ l}relatjve am) absolute mnimmn deviacioRf ben>cen
C— model and Geld data and day of occurrence
C— 2) slope of regression of Held dau on simulation results
C- 3) regression coefticiem
C— 4) standard error of the regression
WRJTE(8,3000)
DO 101 1=1.10
ff(XSO(I).GT. 0.) B(I)=SUMXY(I)/XSQ(I)
IF(IFLAG(I).GT.2) THEN
SUMXY(l)=(YSQ(l)-SUMXY(l)'SUMXY(!)/XSQ(l))/(IFLAG(l)-2)
RSCKI)-L-SUMXY(I)'
+ IFLAGCO-OFLAGt^iydFLAG
SUMXY(I)=SQRT(SUMXY(I))
ENDIF
101 CONTINUE
WRTTE,7X'rT,«X,IBOD-.7X.-DO1.
+«X.TSS'.«XTDS-.6X'N03'.OC,'NH4'y)
3013 FORMAT(LX,'MAXIMUM RELATIVE ERROR ^MOXlOtSXF-t.O))
3014 FOR.MAT(1XTJATE OF MAX REL, ERR,M«.10(4XI2.r.I2))
3015 FORMATCIX^MAXIMUM ABSOLUTE ERRORM4X10(2XF7J))
A-37
-------
301* FORMAT(LVDAT£ OF MAX ABS. ERR.M«.10(XSQOO),
+YSQ(10).RSQ{10).RMS<10).R£LM(10)_MTHREmO).MDAYREUIO).ZREmO).
+ZRELM(10).RS(10XREU10).MTHR.MS{10XMDAYRMS(IO)JZRS(10).ZRMS(10)
COMMON/VOUZMAXDZ<40).Z<40XA(40).V(4«).TV(40)ATOP(41).DBL
COMMON/STEPSCZlJ.DZin.MSarNM.NPRirn'.MDAy.MONTHILAY.DY
COMMON/CHOICEMODEL.NJTRO.IPRNT(<).NDAYS,NPR,ST(>0),Naj^SS,PIjOT(30)
COMMON/FIELD/ IFLAOCIO).FLDATA(:a50).DEPTH(50).NFLD(10)
DIMENSION COMP(40.10)
' EQUIVALENCE (T2(1).COMP( 1.1))
CHARACTER'16 DIN.MET.F1X).TAPES.TAPE1
DO303IetlO
RS(1)=0
303 RELd)-0
DO304je!.20
DO304I = UO
304 FU)ATA(U)=0.0
READ(?.*)SF.NPRFLE
NDAYS-NDAYS*!
IFCNF.GT.O) THEN
READ (7.')(NFUX1).1=1.NPRFLE)
READ(7.-)(DETrH(l).I-= LNF)
DO30SI-J.NPRFLE
READ(7.-)(FLDATAf NFLD(I)J)J - 1.NF)
305 CONTINUE
U,»l
C— Locate simulation values corresponding to sampled
C»consu'tuenu and depth of field data
DO 310 KS= 1.NF
L=LL
DO 315 LL<«l_MBOT
IF(2XI.T.DEPTH(KK)) GOTO 313
IF(U-E0.1) THEN
DO 33) 12= l.S'PRFLE
I=NFLD(I2)
IF(FLDATA(LKK).GT.O)THEN
XX-COMP(LUI)
CALL STATS(FLDATA(LKKPOC.IFLAG(I).DEPTH(KK).I)
ENDIF
320 CONTINUE
ELSE
DO330I2=LNPRFLE
IF(F_DATAO.KK).CT.O) THEN
XX-COMP(LL-LI) + (DEPTH(KK)-Z(LL-l))/fZ
-------
WRrrE(8>3015) (RS(l).l-UO)
WRrTE(8,3016) (ZRS(I).I=U10)
C_Siore dau on plot file (Upe&PLT)
IFOPRNT(5).GT.O) THEN
WRTTE(1.REC=IREC) REAL(NF)
IREC-IREC+1
WRTTE(tREC=IREC) REAUNPRJFLE)
[REOIREC + 1
DO5001-LNF
WRTTEdREColREC) DEPTH(t)
500 IREC-IRE01
DOJ011-1.NPRFLE
WRTTEaREC-IREC) REAL(NFLD{I))
SOI IREC=IREC+1
DOS0212-1.NPRFLE
DO5021-UNF
WRrTE(l.REC-!REC) FLDATA(NFLD(t2).l)
502 IREC»IREC+1
ff( NFLD(l).NE.l) THEN
X-0.0
ELSE
C-Mixcd bycr depib in Held diu ukcn ai dT/dZ- 1.0
DO 503 J2=1NF
X=(FLDATA(U:-D-FLDATA(U2)V
+ (DEPTH(J2)-DEPTH(JM))
IF(XGT.LO) GOTO 504
503 CONTINUE
504 X=(DEPTH(J2)*DEPTH(J2-1))'0.5
IF(JZGENF)X-ZMA.X
END1F
WR1TE(1.REC=IREO X
IREC=iREC»l
ENDIF
ELSE
NF=0
IF(IPRNT(5).GT.O) THEN
WRrTE(l.REC=lREC) REAUNF)
KEC=IREC+1
ENDIF
ENDIF
29» FORMAT(1XI«XI5)
3009 FORMAT(///5X,'DATE: \if f ;aj/SX.
+SUM OF ERRORS BETWEEN MODEL AND FIELD DATA1)
+-LAKE TEMPERATURE STATISTICS'/)
3013 FORMATfUfMAXIMUM RELATIVE ERROR (%)MOX.10(5XF4.0))
3014 FORMATCIXUEPTH OF MAX REL. ERR. (M)'.9X10(5XF4.1))
3015 FORMATWMAXIMUM ABSOLUTE ERROR'.14X10(2XF7J))
30W FORMAT(1XT)EPTH OF MAXABS. ERR. (M)'.9X10(5XF4,1V//)
3020 FORMAT(/42XTEMP'.iX1CHLA>.7XTP-.7XTXn
RETURN
END
C
SUBROUTINE PTABLE
C"***
C"'" Send simulation raulu to the tabular
C"— output file (upeSDAT)
C*****
COMMON/RESULT/ T2<40).CHLATOT(«)).PA2<40).PTSi;M(40).BOD2<40),
+DSO2(40).C2(40),CD2(40)JCNO2{40)JCNH2(40).CHLA2(3.40),
+PC2(3,40)JCNC2(3.40).T20(40),SI2(40)
COMMON/VOL/ZMAXDZ440).Z(40) A(40).V(40).TV(40) j\TOP(41 )JDBL
COMMON/STEPS/'DZLL.DZUL.MBOT.NM.NPRINTJUDAY.MONTHILAY.DY
COMMON,'CHOlCE'MODEL.NrrRO,IPRNT(«).NDAYS,NPRNT(30)J
-------
IF(NrTRO.GT.O) THEN
ELSE
ENDIF
WRITE(a3006)
IFCNITRO.Crr.O) THEN
ELSE
EVD;F
ENDIF .
1200 FORMAT(///5X'ZOOPLANCTON PARAMETERSV-SX^v
+ 'CONG (*.3X'GRAZING(MG/M3)-.3X
* ' DAYDEPTH(M)V.9Xn.0.9XF7.1.9XF&4.1IXF4.1.'/)
1201 FORMAT(Fii.«} .. . ..
399 FORMAT(//5XlNmAL CONDmONSVSX18(Vy) .'
3000 FORMAT(//5X1NFORMATION ON LAKE QUAUTYVSX27r>V)
3001 FORMAT^X'ZfM) T(C) SS(PPM) TDS(PPM) CHLA(PPM) PC(PPM)',
+ ' P(PPM) PT(PPM) BODCMG/L) DO(MG/L) DZ(M)'.
+ ' AREA(M2) VOUWJIT
3002FORMAT(5X3(F5.2.2X).lXFS.J,4X4<2XFSXF<2.
3006 FORMATCWX-BIOLOGtCAL PARAMETERSV5X21C-VJX'Z(M)-.»X
+-CHLOROPHYLL'.«XTOTAL-.6X1NERNAL PM1X1NTERNAL SV16X
+-V.«XT.6X-3'JX1CHLA-.4XT.«X7-.«X'31.6XT.«XT.
3007 FORMAT(5XFJ.2.2X10(F5J.2X))
3008 FORMAT(//.5X' TEMPERATURE PROFILESV.5X'Z(M)-.
p****«
C**"** Prtxiucc profile of slate variables and
€*•••* field data (if niilable)
^••••»
CHARACTER'1 XS.CH1.CH2,CH3.CH4
CHARACTER'S; TITLEO)
CHARACTER'4 MNTH(15
COMMON/SOURCE/RM(3,40).PROD(40)JCMR(3.40)J'RODSUM(«)
COMMON/STEPS'DZLL.DZUUMBOT.NM.NPRINT^IDAY.MONTHILAY.DY
COMMON/RESULT/ T2(40).CHLATOT(40).PA2(40).FrSUM(40).BOD^-IO),
+DS02(40).C2(40).CD2(40-^XNO2(40).XNH2(40XCHLA2(3.«».
COMMON/F1ELD/1FLAG(10).FLDATA(10.50).DEPTH{50)>.'FLD(1I')
COMMON/\'OL.'ZMAXDZ(40).Z(«))A('«)).V(40Vn'(-<0)j\TOP(41).DBt
DIMENSION ZD(23).FD(23).\TM(43).ZV(43).VAR(«).IO)
DIMENSION FCT(10).LEN(10).1SCAL(10).ICX(4)
EQUIVALENCE (T2(1XVAR(1.1».(2D(1).PROD(1)).(FD(1).PRODSUM(I)).
EQUIVALENCE (CHUCX(1)MCH1ICK(2)).(CH3.ICX{3)).(CH4.IC.\(4))
DATA lCX/16*lB.16*5B.iw»3il6»-4A
DATA BCAL/UMJ-U'-I/
DATA Fcr/u*iooo,4'i.ri(x»
DATA MNTWJAN '.'FEB '.TvlAR '.'APR '.TklAY VJUS VJUL '.'AUG '.
+ 'SEP '.XXTT •,'NOV '.T)EC V
DATA TTTLE /TEMPERATURE (QV
DATA LEN /1^27.2&25.23.<-24.2S/
1000 WRTTEC.109) CHLCH2.CH3.CH4
C-IUl and idea desired suie variable for plotting
DO100I=1.1
100 U-RITE{M01) l.TrTLE
-------
110 CONTINUE
12-0
C_if Geld data b nibble, locale Held data corresponding
C-to selected state variables
IF(NF.CT.O) THEN
DO 111 1-tNF
IF(FLDATA(ICLI).GT.O) THEN
12-12+1
. FD02)=FU>ATA(ICU)
ZD(I2)--DEPTH(I)
IF(X,LT.FD(12)) THEN
X-FD(12)
IND— 1
END1F
ENDIF
111 CONTINUE
ENDIF
NPEN-1
NPEN1-1
IPORT-99
MODL=W
ZV(MBOT+1)-0.
ZD(I2+1)=0.
VTM(MBOT+lj=0.
FD(I2+1)=0.
FCTOR.-0.7I
C_begin plotting sequence
990 CAUL PLOTS(O.IPORT,MODL)
CA1JLSIMPLX
CALL FACTOR(FCTOR)
CALL NEWPEN(NPEN)
C-dflenninc maximum i &. y values either in simulated
C— atate variables or in fieid data
IF(IND.GT.O) THEN
13=MBOT+1
CALL SCALE(VTM,10,13. 1)
DX=VTM(MBOT+3)
ELSE
13-12+1
CALL SCALE(FD.10.,13.1)
DX=FD(U+3)
ENDIF
lF(ZV(MBOT)rr.ZD(I2)) THEN
D-MBOT+1
CALL SCALE(ZV.i,I3.-l)
YSC— ZV(MBOT+3)
ELSE
D=I2+1
CALL SCALECZD.4,D.-1)
YSC— ZEKI2+3)
ENDIF
ZV(MBOT+2)»YSC
ZD(I2+2)-YSC
VTM(MBOT+2)-DX
FD(I2+2)-DX
DAX1S=DX*FCT(IC1)
CALL STAXIS(i.I5.2,15.ISCAL(ICl})
C_draw x-axis
CALL AXlS(U,TTTLEaCl)a£N(ICl).-10,a.O,DAXJS)
CALL STAXB(.2,25_2_li,l)
C-dnw y-axis
CALL AX1S(U1,T)EPTH (M)MO.-^9a.YD.-YSC)
XMDAY=MDAY
XMYEAR-MYEAR.
C_prim title to diagram
CALL SYMBOL(4J J..25.MNTH(MONTH),0,4)
CALL NUMBER(999,999_2S.XMDAY.O.-1)
CALL SYMBOU9M,999_2S.-. '.0. 2)
CALL NUMBER(W,999_25.XMYEAR_0,-1)
CALL NEWPEN(NPENl)
CALL PLOT (U7,-3)
C—plot simulated profiles as a line
CALL UNE(VTM.ZV,MBOT.L
-------
WRTTELO ) : 'A)
GOTO 990
ENDff
GOTO 1000
1001 RETURN ' .
END
SUBROUTINE ABOUND
C—
c—
c—
C*"
c—
c—
COMMON.'STEPS.'DZLUDZUl-MBOT.NM.NPRINT.MDAY.MONTH.Il-AYDY
COMMOS.VOL/2MAXDZ<40).Z( 40)^(40). V(40).TV(40),ATOP(41).DBL
DUM-tt
DO1001»1.MBOT
/ ZDUM=ZMAX-DUM
DUM=DUM+DZ(I)
CALL LAKE(ZDUMw\DUM0.1)
100 ATOP(1)=ADUM
ATOP(MBOT+l)-a
RETURN
END
SUBROUTINE AREA
c*****
C""" Compute ihe area through the middle of each layer
C***** using the depth-area relationship in LAKE
C*****
COMMO.VSrEPS^>2U.D2UL.MBOT.NM.NPRJNT.MDAY.MONTHJlJ\Yjy
COMMONA'OlJ2MAX.DZ<40).Z(40)j\C40).V(«0).TV(JO).ATOP(41).D3L
DUM>a
DO1001.1.MBOT
ZDUM-ZMAX-DUM-DZ(!)/1
CALL LAKE(ZDUMADUM,0,1)
100 A(I)«ADUM
RETURN
END
SUBROUTINE COEF(MODELMBOT.NCLASS)
^••••«
C*""** Compute some coefficients used in the constant
C"**** volume and finite difference solutions
C*****
COMMON/COEFF/ DUM^40).DUM3(40)
COMMON/VOUZ.MAXDZ(40).Z(40)JV(40).V(40).TV(40).ATOP(41).DBL
COMMON/FLOW/ HMK(J!).OE(40).FVCHLAr5).PE<5.41)
DOlOOi.lMBOT-1
100 DUM^D^DUMI'ATOPd +
KK»2
IF(MODELGT.l) KK=1
DO 300 K=KK.NCLASS+1
DOaWi-iMBOT
X«FVCHLA(Kr(DZ(I-l)*D?(l))'J/HMK(I)
C - PE-tLO-.rABSfX))".^ -
AO.L0..1-ABS(X)
A1=AO-AO
PEfK.I)=AI'ArA(VX
PE(K.l)-aO
200 PE(KJHBOT+1)«0.0
RETURN
END
SUBROUTINE CONMLX(ILAY.TMIX.MBOT)
£•••••
C***** Remo\e deniity insubilities by mhang umublc
C""* \rrcn dcwnuranl and mcrjjngunth lower layen.
^••••*
COMMON RESULT/ 17(40).CHLATOT{40).PA2(40).PTSUM(40).BOD^40).
+ DSO««).C2(40).CD2(40).XNO2(40)JCNH2(40).CHLA2(3.40).
+PO(3,40).XNC2(3.40).TIO(40).SI^40)
COMMON VOL/ZMAXDZ(40).Z<40)JV(40).V(40).TV(40)J\TOP(41).DBL
A-42
-------
DIMENSION RHOT(40)
DO100I-UMBOT
:» RHOT(I)-RHO(T2(l).a.O.)
6 1FLAG-0
1-0
M-MBOT-l
1 1-1*1
IF(LEO.M) GO TO 5 .
IF(RHOT(l).LE,RHOT(ltl))GOTO 1
IFLAG-1
IB-I
TVDUM-T«I)'V(I)
VDUM-V(I)
TDUM=TVDUM/VDUM
RHODUM=RHO(TDUM.a.O.)
J=!-l
3 J=J+1
IF(RHODUM.LE.RHOT(J + 1)) GO TO :
1FLAG.1
VDUM»VDUM+V(J + 1)
TDUM-TVDUM.'VDU.V!
RHODUM = RHO(TDUM.O..a;
x IF(J.EO.M) OO TO I
GO TO 3
2 IE-J
!F(J.EO.M)1E"IE + 1
IF(iaEO.IE) GO TO 4
T2(K)=TDUM
y» RHOT(K) = RHODUM
4 IF(I.SE.M)GOTO1
5 IF(IFLAG.EQ.1)GOTO6
C-DETERMINE MLXED LAYER DEPTH...
DO700I = 1.MBOT-1
IF(CnWOU2.MAXDZ<40)i(40)J\(40).V(40).TVC40)w\TOP(41).DBL
COMMONi'RESUUT/T2('W),CHLA.TOT(40).PA2(40).fTSUM(40).BOD2(40),
+PO(3.40)JCNC(3.40).Ta)(40)^C(40)
COMMONOiOIC£/MODEL,NrrRO,]PRNT(6).NDAYS.NPRNT(30),NCUASS.PLOT(30)
ff(V(I)J-E.O.) THEN
!F(LEQ.M8OT) THEN
D-MBOT
MBOT»MBOT-t
IF(V(II-l).LE,aO) THEN
II-IM
GOTO:
ENDIF
RETURN
ENDIF
V(Itl)-V(I)»V(l + l)
DZZ«D2(0
KK=!
MBOT=MBOT-1 >
ELSE
U=I
!F(LEQ.MBOT) 11=1-1
KK»1! + 1
VC=V(IO+V(KK)
VCOMB=UVC
CD2(n)=(CD^SI)*V(l!)+CDZ(KK)'V(KK))*VCOMB
DOS5K=1.NCLASS
55 CHLA2(K.II)«(CHLAKK.n)-V(II)+CHLA2(K.KK)%V(KK))>VCOMB
A-43
-------
IF(MODELEOJ) THEN
DOS7K-1.NCLASS
57 PC2(O)»(PC2(K.BrV(7I
ENDIF
PA2(H)=CPA2(II)'V(ir)+PA2(KK)'\'(KK))'VCOMB
BOD^I1)={BODZ(II)*V(II)+BOD^KK)*V(KK))'VCOMB
DSO2(II)=(DS02(IU'V(II)+DSO2(KK)'V(KK)}'VCOMB
IF(NITRO.E0.1) THEN
DOS«K=1.NCLASS
56 XNC2(K.n)=CXNC2(K.II)>V(Il)+XNCKK.KK)tV(KK))'VCOMB
ENDIF
v(ii)=vc
DZ(1I)=DZJBOT
T2(K)-T2(K+1)
. 0).
* DSO2(«)).a(40).CDI(40)_XNO^«))_XNH2(«t.CHLA^3.<0).
COMMON/TEMP/PARIO(4).PCDUM(3.40)JCVHD(«VXNOD(40).
+ CHLADUM(3,40)JCNCD(3,40).PADUM(-W)51D(40)
COMMON/STEPaDZU.DZUl-MBOT.NM.NPRlNT.MDAY.MONTH.ILAY.DY
COMMON/FLOW/HMK(41XOE(«)JVCHLA(5).PE(5.41)
COMMON/WATER/BETA.EMISSJCK.tJCKZHKMAX.WCOEF.WSTR
COMMON/SOLV/ AK(*0).Bm40),CK<40).DtC(<0)
DIMENSION CK40)
C.-.CALCULATION OF THE HEAT ABSORPTION FROM METEOROLOGICAL PARAMETERS
C_IN A COLUMN OF WATER
C-.CALCULATION OF HEAT FLUXES INTO THE WATER BODY
CALL FLXIN(HS,HA,TAIR,RAD.CR.C2)
CALLFLXOUT(TS.HBR.HE,HC,TA1R.TDEW.W1ND.IPAN.PCOEFF.WCOEF)
HOOUT-HBR»HE+HC
C..CALCULATJON OF EXTINCTION COEFF. (ETA) AS A FUNCTION OF SUSPENDED
C-SED1MENT CONCENTRATION
ETA=l.W(l/SEia)
C.CALCULATION OF HEAT ABSORBED IN EACH LAYER
HO=(L-BETA)'HS
EX=EXP(-ETA'DZ(1))
A-44
-------
C-CONVERSION FACTOR OF 1000 USED FOR DENSfTY-HEAT CAPACITY OF WATER.
HQ-HO/EX
C-CALCULATE THE SOURCE TERM Q FOR EACH LAYER
D010U2.MBOT
EX-EXP(-ETA'DZ(I))
ODZ.:ir(DZ
-------
HE=FCN'(ESA-EA)
OOTQ30
20 HE=15-lTD'PCOEFF*59a'ia
C-.CALCULATES BACK RADIATION
30 HBR=*(U71E-«*0.»7TSK.TSKTSKTSK)
C_CALCULATES CONDUCTIVE LOSS USING BOWENS RATIO
HC=FCN-0.617W(TSK-TAK)
RETURN
END "
SUBROUTINE SETAMKCW.HKMAX11AY.MBOT)
C"
' Compute venial diffusion coefficient in each layer.
** Diffusion coefficient berween layers as tbe bannonic
1 mean of tbe diffusion coeffidenu in adjacent layer*.
»
DIMENSION AMK(41)
COMMON/FLOW/HMK(41).OE(40).FVCHLA(5).PE{5.41)
COMMON/VOU2.MAXDZ(40).Z<40)A(40).V(40).TV(40)ATOP(41).DBL
COMMON/RESULT/T2(40),CHLATOT(40).PA2(40).PTSUM(40).BOD2(40X
C"
C"
c-
C"
+PC2a«0)_XNC2(3.40).T20(40).SI2<40)
DIMENSION PSON(40)
AREA-(ATOP(iy(10"«0)
ALFA=ai7*O.OOOr(AREA)"0.5«
C—Venical diffusion coefficient in tbe mixed layer
C— computed as a function of tbe wind speed
DKM«2S.*WL3
IF(DKM.GT.HKMAX) THEN
SON =0.000075
DKM«.(ALFA/((SON)"0.43)r&«
ENDIF
C
D0100I=ULAY
100 CONTINUE
C
C 929 FORMATtlXHYPOUMNION1)
C_RMK TEMPORARILY USED TO STORE DENSITIES.
C_diffu>ion coefficient below tbe miied layer computed as
C_a function of the square of tbe Brunt-Vaula frequency (SON)
C_ SON = (G/RHO) ' d(RHO)/dZ
DO 200 1=ILAY.MBOT
200 HMK(I)-RHO(T2(I).C2(1).CDI(I))
C WRTrE(99.«)
DO 300 lcILAY + l.MBOT-1
SON = ABS(HMK(I-1)-HM1C(I + 1))/((Z
-------
20 CONTINUE
DO 1WK-LKDAYS
READ(9/)TA]R(K).TDSW(K).WIND(K).DRCT(K).RAD(K).CR(K),PR(K)
PRdC)-PR(KyiOa
104 CONTINUE
DO 100K»l,KDAYS
TAIR(K)=(TAIR(K)-31)'0.5S56
TDEW(K)-(TDEW(KV31)-aS55<>
100 CR(K)-(100.-CR(K))W.01
C
C MAKE CORRECTIONS FOR CLIMATE MODEL OUTPUT
C IF NEEDED. THESE ARE MONTHLY ADJUSTMENTS.
C
C DO1287K=1.KDAYS
C IF(MONTH.EQ.3) THEN
C TAIR(K)=TA!R(K)+4.SO
C TDEW(K)=TDEW(K)»7.25
C WIND(K)=WIND(K)-aS2
C RAD(K)=RAD(K)>L04
C PR(K)»PR(K)*1-02
C ENDIF
C !F(MONTH.EO.4) THEN
C TAIR(K)=TAiR(K1 + 497
C TDEW(K)=TDEW(Kl-i.40
C WlNDPR(Kri.!4
C ENDIF
C IF(MONTKEO.ll) THEN
C TAIR(K)=TAIR(K)+5.93
C TDEW(K)=TDEW(K)»5.51
C WIND(K)»WIND(K)'1.0S
C RAD(K) = RAD(K)'0.95
C PR(K) = PR(K)'L1«
C ENDIF
C 12S7 CONTINUE
1001 FORMAT(//.SX.70f')^.MX'PROGRAM ABORTED.V.10.X.
+ 'METEOROLOGICAL DATA FILE DOES NOTJ.15X.
+ MATCH YEAR OF SIMULAT.ONV.'.5X.70(~))
A-47
-------
RETURN
END
SUBROUTINE PMETE(HS.RAD.HA,WIND.HBR,P.HETA1R,HC.TDEW.HED>.HE.TA1RJ1CTDEW.HED.HEV,
+QNET.DMIXZEUPRSECCHI
C
MOO FORMAT(/.JX26HMETEOROLOGICAL INFORMATION./.5X26(1H-).
+/.TX17HNET SOLAR RAD. = .4XF9.2.13H KCAL/M'M/DAY.
•f Z5X13HSOLAR RAD. - , 9XF6.1.8H LANGLEY,
+/.TX16HNET ATM. RAD. - 4X.F9.113H KCAUM'MT>AY,
+2SXWHWIND VELOOTY - , 7XF5.17H M.P.H,
+/.7XUHBACKRAD. « . 9XF9.113H KCAL/M'M/DAY,
+2SXl«HPREaPITATION - . 8XF«J.13H METER(S).DAY.
+/.7X19HEVAPORATTVE LOSS • ,IXF.9Z13H KCALM'M/DAY,
+25X12HA1R TEMP. - .12XF5.Z10H DEGREES C
+/.7X20HCONDENSATION LOSS » . 1XF9.Z1JH KCAUM'MOAY.
. +25X16HDEW PT. TEMP. « . SXF5Z10H DEGREES C
+/.7X20HEV. HEAT TRANSFER « ^XF7.<.10H IWDAY OR. .Eltt^.TH M"3.T)
+2HAYJ.7X19HNET HEAT FLUX IN - . 2XF9.Z9H KCAL/M'M.
+/.7X1 MKED LAYER DEPTH (M)«'.ZXF5.i
+36X1EUPHOTIC DEPTH (M) «'.2XF5.J
+/8X'SECCHI DEPTH (M) ='.7XF5.2)
RETURN
END
SUBROUTINE SETUP
C"
C
C
C
C
1 Determine the initial thicbnos. volume. >nd Area of l»v«n
' and the total vc4ume of above each layer from the depths given
" in the input dau Hie.
COMMON/STEPS'DZLL.DZUL.MBOT.NM.NPRJNT.MDAY.MOSTH.ILAYI>Y
DZ(MBOT)«ZMAX-(Z(MBOD +Z(MBOT-1))' J
Z
-------
SUBROUTINE SOLVE(VARiMBOT)
C*****
C*"** Tri-dbgonal matrix solving routine
C*****
COMMON/SOLV/ AK(40).BK(40)1CK(40).DK(40)
DIMENSION VAR2(40).TX(40)
DO 60 I-2.MBOT
TT»AK(I)/BK(M)
BK(1)-BK(I)-CK(I-1)TT
«0 DK(I)=»DK(I)-DK(1-1)TT
C"— — "BACK SUBSTITUTION"——"""
TX(MBOT)=DK DZUL) into iwo or oore
C***** layers of equal thickness. All slate variables are (he ume in
C***** each new layer. Volume is adjusted later.
£•••••
COMMON/VOL/ZMAX.DZ(40).Z(*))>\(*)),V(40).TV(40}ATOP(J1).DBL
COMMON/CHOICE/MODEL.NITRO.IPRNT(«).NDAYS.NPR.NTC50).NCLASS.PLOr.:-0:
COMMON/RESULT/ T2(40).CHLATOT(40).PA2(40).rTSUM(«)l.BODajl<').
+ DSOy40).C2(40).CD2('«)),XNO2(40).XNH:(«).CHLAi-3.JOV
»PC2(3.40),XNC2(3,'10).T20(40),S12(40)
COMMON/STEPS/DZLL.DZULMBOT.NM.NPRINT.MDAy.MOSTH.!LAY.DY
DO100K=I.MBOT
n=MBOT+I-K
CD2.40).T20(40),SI2(-W)
COMMONrtEMP/PARJO(4),PCDUM(3,40),XNHD(T(30).VCLASS,PLOT1'30)
COMMON/CHANEUWCHANL.ELCBjUJ'HA.BW.WLAKE
COM.MON/WATER^ETA.EMISSJCKLXK2,HKMAXWCOEF>STR
COMMON/STEPS/DZLL,DZUI_.MBOT.NM.NPRJNT.MDAY.MONTH.ILAY.DY
COM.MONA'OUZMAXDZ<40).Z(-tO)^(40).V(40).TV(40)^TOP(-
-------
NITRO=0
DO 200 1-1.30
200 NPRNT(I)=0
C* INPUT MODEL OPTIONS AND INITIAL CONDITIONS ••"••
READ(7.') SEKI
WR]TE('.1006)
READ(V(A)-) X
IF(X.EQ.-r .OR. XEQ.y) THEN
1PRKT(2)=1
TAPE1=TAPE8
Tl(LEN8+2) = -O-
Tl(LEN8+3)-ir
OPEN'(1FILE=TAPEISTATUS=-NEW)
ENDIF
IPRNT(4)=0
1000 FORMATC PLOT FILE TO BE CREATED (Y/N) ? 'A
1001 FORMAT(A1)
1004 FORMATC ENTER UP TO 10 DEPTHS TO BE SAVED V.
+' END WITH A CHARACTER (».*.*__X>: V)
1006 FORMATC OUTFLOW FILE TO BE CREATED (Y/N) 1 V)
READ(7,') NDAYS
IF(NDAYSGT.O) READ(7,')(NPRNT(1).1-1,NDAYS)
READft*) DZLLDZULBETA.EMISS,WCOEF,WSTR
READ(7,') WCHANL.WLAKEDBLST.S.FT
. READ(7,') ELCB.ALPHA.BW
• READ(7,') MBOTJ4M.NPR1NT.INFLOW.DY.MONTH1START.MYEAR
READ(7,-) (Zfl).l»l.MBOT)
READ{7,*) (T2(I),I=1.MBOT)
1500 FORMAT(/, 5X
+/JX41HMIN1MUM THICKNESS OF EACH LAYER (DZLL) - J5Z»H METERCS)
+yXBRR'JX-FVBODVJX
+-(iyDAY)'.2(2X-(GM/M2)-)JX-(MS)ViX5(nj.2X))
RETURN
END
SUBROUTINE STATS(FLDATAJCXIFLAG.DEPTH,I)
C*****
C***** Compute lUiiuia >nd tuiaitical quaniiua with
C***** Y daipialinj field daa tad X daiputinf model raulu.
C""***
COMMON/STAT/SUMXY(10)5UMX(10)5UMY(10)_XSO(10).
+YSQ(10).RSQ(10).RMS(10).RELM(IO).MTHREL(10).MDAYREL(10).ZREL(10).
+ZRELM(10).RS(10).REL(10).MTHRMS(10),MDAYRMS(10).ZRS(10).ZRMS(10)
COMMON/STEPS/DZLL.DZUL.MBOT.NM.NPRINT.MDAY.MONTH,ILAY.DY
SUMXY(I)=SUMXY(I)+FLDATA'XX
SUMXO)-SUMX(I)*XX
SUMY(I)=SUMYa)+FLDATA
XSO(I)=XSO(I) »XX'XX
YSO(I) - YSQ(I) * FLO ATA' FLDATA
XX=XX-FLDATA
X2=XX'PLDATA'10a
A-60
-------
X3»ABS(XX)
!F(XlGT_ABS(REm))) THEN
R£L(I)«X2
ZREHO-DEPTH
ENDIF
ff(XlGTJ\BS(RELM(l))) THEN
RELM(I)-X2
MTHREm)=MONTH
MDAYREL(t)=MDAY
ZRELO)-DEPTH
ENDIF
ff(X3.GT.ABS(RS(l))) THEN
RS(I)-XX
ZRS(r)"DEPTH
ENDIF
tF(X3.GTj\BS(R-MS(I))) THEN
RMSMMONAfOL/2MAXDZ(«),Z(-W).A(40).V(40).TV{40)^TOP(41).DBL
AZ=tt
AVOL=a
DO100I=1.MBOT
II=MBOT+1-I
AVOL=AVOL+V(II)
CALL LAKECZDUM^VOL.0.4)
DZ(II)=ZDUM-AZ
A2=AZ+DZ(II)
100 CONTINUE
RETURN
END
SUBROUTINE TVOUMBOT)
C**"**
C***** Determine the volume of water above a layer
C*****
COMMON/VOUZMAX.DZ(40).Z(-iOyv(40).V(40).TV(40)^TOP(41).DBL
SUM=tt
D01001-1.M80T
SUM=SUM+V(I)
TVO)=SUM .
100 CONTINUE
RETURN
END
SUBROUTINE W1NEN(T.V,WIND)
C_
C-CA1.CULATIO.N OF THE SHEAR VELOCITY AND THE WIND SHEAR STRESS
C-CONVERSION OF WIND SPEED FROM M.P.R TO M/S
C-DENSITY OF WATER AND AIR ASSUMED TO BE 1000 AND L177 KC/M3
C_
CALL LAKE(FTCH.aao,:)
ZB-ALOG(FTCH)'a8- 1.0718
C-CALCUUAT10N OF WIND SHEAR STRESS
CZ=SQRT(W)'.0005
ff(W.GE35.) C2=.0026
C-ASSIGNMENT OF CHARACTERISTIC SURFACE VELOCITY
c-usiNG CALCULATION OF SHEAR VELOCITIES
V=.03J3>SORT(CZ)'W
RETURN
END
SUBROUTINE WINMLX(ETS,IL.MBOT)
C_
C.CAJXULATE THE AMOUNT OF ENTRAPMENT RESULTING FROM WIND MIXING.
C-USE THE DEPTH OF CENTER OF MASS OF MIXED LAYER TO DETERMINE THE
•C-POTENTUL ENERGY THAT MUST BE OVERCOME BY THE KJNECTIC ENERGY
C-OF THE WIND FOR ENTRAPMENT TO OCCUR.
C_
C»MMONA'OUZMAXDZ(«),Z<«).A(40).V(«).TV(«)XrOP(41).DBL
COMMON/RESULT/ T2(*)).CHLATOT('(0).PA2(«)).PTSUM('W).BOD2(40),
A-51
-------
+DS02(«),CI(40),CD2(40),XN02(40)JCNH2(40).CHLAI(3.«).
+ PC2(3,40)XNC2(3.40).T20(40)312<40)
!F(ILGEMBOT) GO TO 35
SUM2=0.
DOI01-1.1L
ILAY-l
SUM1=SUM1+RV
10 SUM2=SUM2+RVZ(I)
1 = ILAY , . .
20 DCM=SUMZ
-------
LAKE SPECIFIC SUBROUTINE
Area computation section
Depth-area functions of the form AREA=f(ZMAX-Z), written as DUM=f(ZD).
AREA(m:), Z(m).
Fetch computation section
The longest distance across the lake surface area in the wind direction (m).
Volume computation section
Volume-depth functions of the form VOLUME(below depth Z)=f(ZMAX-Z), written as
ZD=f(DUM). VOLUME (m3), Z(m).
A-53
-------
EXAMPLE LAKE SPECIFIC SUBROUTINE
SUBROUTINE LAKE(ZD.DUM.NFLOW.ID)
COMMOWMTHD/rAJR(3I).TDEW(3I).RAD(31).CR(3I).WtND(31).
+ PR(31).DRCT(31)
COMMON/RESULT/ T2(40).CHLATOT(40).PA2(40).PTSUM(40).BOD2(40),
COMMON/SOURCE/RM(3.40),PROD(40).XMR(3.40).PRODSUM(40)
COMMONVFL.OW/HMK(4I).pE(4b).FVCHLA(5).PE(5.4I)
CO.MMON;TEMP/PAR1O(4).PCDUM(3,40)XNHD(<0)_XNODC-10).
+ CHLADU.V{3.40),XNCD(3.40).PADUM(40).SID(40)
COMMONAX3L^MAXDZ(40).Z<40)^(40).V(40).TV(40).ATOP(41).DBL
COMMONJSUBJSDZ<60).SZ(«).LAY(40).AVGI(4.40).SVOL<«>)
COMMON,CHOlCE/MODEl-NrniO,IPR>!T(«).KDAYS.NPRNT(30).NCLASS.PLT(30)
COMMONAVATER^ETA.EMISSJ6-ASURF'(ZD+l.)"l.'<«209
RETURN
C* ......... FETCH COMPUTATION SECTION
200 ZD=»crt(O.l<»59>ASURF)
RETURN
C* ......... VOLUME COMPUTATION SECTION
300 CONTiNUE
1234 ASUR?=:.?1E«
CONS=0.00>J'ASURF
RETURN
C"— COMPUTE DEPTH FRO.V VOLUME .....
400 CONS=0.005«4'ASURF
DUM1=DUM'CONS
ZD=(DUM1"0. 406)6). 1.
RETURN
C ..... WRITE ON THE SCREEN. DAY. MONTH ••••
500 WR1TEC.S01) MONTH.MDAY
501 FORM.-kT(2X' monthM3.' dj> -.i2)
C ..... WRITE EPILIMNION i HYPOUMNION TEV?ERATL!RES ••
DO IIIH I=1.MBOT
IF (LEO.:) THEN
WRITE (2L8TI) MONTH.MDAY.Z(1).T2(I)
ENDIF
IF (LEO^2) THEN
WRITE (n.871) MONTH.MDAY.Z(I).T2(I)
ENDIF '
871 FORM.AT(1X.I4.1XI-1.2XF9.2.3X.F9.3)
11111 CONTiNUE
RETURN
«00 RETURN
700 RETURN
800 RETURN
900 RETURN
1000 RETURN
1100 RETURN
1200 CONTINUE
1300 CONTINUE
RETURN
END
A-54
------- |