EPA-660/3-75-005
MARCH 1975
Ecological Research Series
)deling of Phytoplankton
in Lake Ontario
1. Model Development and Verification
National Environmental Research Center
Office of Research and Development
U.S. Environmental Protection Agency
Corvallis, Oregon 97330
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development,
U.S. Environmental Protection Agency, have been grouped into
five series. These five broad categories were established to
facilitate further development and application of environmental
technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in
related fields. The five series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
This report has been assigned to the ECOLOGICAL RESEARCH STUDIES
series. This series describes research on the effects of pollution
on humans, plant and animal species, and materials. Problems
are assessed for their long- and short-term influences. Investigations
include formation, transport, and pathway studies to determine
the fate of pollutants and their effects. This work provides
the technical basis for setting standards to minimize undesirable
changes in living organisms in the aquatic, terrestrial and atmospheric
environments.
This report has been reviewed by the Office of Research and
Development, EPA, and approved for publication. Approval does
not signify that the contents necessarily reflect the views and
policies of the Environmental Protection Agency, nor does mention
of trade names or commercial products constitute endorsement or
recommendation for use.
-------
ERRATA
Please note the following corrections in EPA-660/3-75-005, MATHEMATICAL
MODELING OF PHYTOPLANKTON IN LAKE ONTARIO:
Page
31 Equation (18)
33 Equation (22)
Change to
a
l
'
. a
,J c
Dp P - Kl Nl
Pl .
1 ,J c
al Pj (1"1
-------
EPA-660/3-75-005
MARCH 1975
MATHEMATICAL MODELING OF PHYTOPLANKTON
IN LAKE ONTARIO
1. MODEL DEVELOPMENT AND VERIFICATION
By
Robert V. Thomann
Dominick M. DiToro
Richard P. Winfield
Donald J. O'Connor
Manhattan College, Bronx,
NY
Grant No. R800610
Program Element 1BA026
ROAP/Task No. 21AKP/15
Project Officer
William Richardson
Grosse lie Laboratory
National Environmental Research Center
Grosse lie, Michigan 48138
NATIONAL ENVIRONMENTAL RESEARCH CENTER
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
CORVALLIS, OREGON 97330
Fof Sale by the National Technical In'ormatton Service
L'.S. Department of t'ennnerjc. Sp';:ig!:e!d, VA 22!5I
-------
ABSTRACT
The basic mathematical structure for describing the
dynamics of phytoplankton in Lake Ontario is presented in this
report. Data on chlorophyll and principle nutrients are re-
viewed and summarized and the mathematical modeling strategy is
detailed.
The modeling strategy begins with the construction of a
horizontally completely mixed lake with vertical layers, LAKE 1.
This spatially simplified model is used to develop the inter-
actions and kinetic behavior of the various components of each
subsystem. A more detailed 3-dimensional model is then used to
describe open lake and near shore variations in phytoplankton
biomass. Ten biological and chemical variables are used in both
models and include four trophic levels above the phytoplankton,
chlorophyll a as a measure of phytoplankton biomass, two
phosphorus components and three nitrogen components. For LAKE 1,
using three vertical segments, a total of 30 simultaneous non-linear
differential equations must be solved while for the 3-,dimensional
model, using 67 segments a total of 670 equations are solved. In
both cases, the models are run for a period of one year or longer.
Under reasonable sets of model parameters as reported in the
literature, the Lake 1 model output compared favorably with observed
data on the dependent variables. Spring growth and peak chlorophyll
concentrations are related primarily to increasing light and tempera-
ture. The model indicates that growth ceases due to phosphorus
-------
limitation. Zooplankton grazing and subsequent recycling of
nutrients due to excretion and plankton endogenous respiration
result in a late summer increase in phytoplankton biomass. Both
nitrogen and phosphorus are then limiting resulting in a broad fall
peak in phytoplankton biomass. A decline then results from low
temperatures due to fall overturn.
The results to date indicate that the mathematical model
of phytoplankton in Lake Ontario as developed herein is a reasonable
first approximation to observed data. As such, the model can form
a basis for preliminary estimates of the effects of nutrient re-
duction programs on Lake Ontario.
This report was submitted in partial fulfillment of Grant
Number R800610 by the Environmental Engineering and Science Program,
Manhattan College, Bronx, New York, under the sponsorship of the
U.S. Environmental Protection Agency. Work was completed as of
October 1974.
111
-------
CONTENTS
Sections
I
II
III
IV
V
VI
VII
VIII
IX
Conclusions
Recommendations
Introduction
Lake Ontario - Background
Lake Ontario Model Structure
Lake 1 Model
Lake 2 Analytical Investigations
'Preliminary Lake 3 Model
Appendix
1
3
4
8
18
58
101
131
159
-------
LIST OF TABLES
1. Flow budget 9
2. Nutrient loadings 11
3. Nutrient load split 11
4, CCIW limnological data report summary
(Number of cruises with adequate spatial
defination) 14
5. Thermodynamic properties at 25 C 44
6. Thermodynamic properties at 25 C 44
7. Temperature dependence comparison 46
8. Basic physical data of the Lake 1 model .... 60
9. Principal parameter values used in
Lake 1 model output shown in figures
21-23 88
10. Log-,0 (G/d+u) for exponentially decreasing
growth rate H5
11. Monthly heat flux values 142
12. Regional to Lake 3 segmentation
equivalence 144
A-l Lake 1 physical parameters 160
A-2 Lake 1
Mass loadings and initial conditions 170
A-3 System parameters
v
-------
LIST OF FIGURES
Page
1. General basin nap of Lake Ontario 5
2. Typical CGIW cruise track, April 1969 13
3. Typical output from data reduction - main
lake average @ 1 meter 16
4. Eutrophication modeling strategy 21
5. System diagram - Lake 1 model 23
6. Relationship of pH and the ratio, Alk/C™ .... 38
7. Per cent error in assuming
(C02 - Acy) = (C02) 41
8. Computations for C02 - EUO system 49
9. Major physical features included in
Lake 1 model 59
10. Effect of vertical mixing 65
11. Effect of phytoplankton settling
velocity, K^ = lOyg/A 67
12. Effect of phytoplankton settling velocity
Kmp " ^g/* 69
13. Effect of phytoplankton settling velocity
on nutrients and nutrient limitation
Kmp = W^/l 70
14. Effect of Michaelis constant for phosphorus
on primary variables 73
15. Effect of Michaelis constant for phosphorus
on nutrient limitation of nitrogen and
phosphorus 74
16. Responses due to two photoplankton
growth formulations 76
17. Sixteen year computed output of
phytoplankton chlorophyll 78
18. Long term behavior of model components under
two phosphorus Michaelis constants 80
19. Preliminary comparison of model output and
observed data 82
20. Preliminary comparison of model output and
observed data, continued 83
VI
-------
21. Lake 1 model verification, 0-17 meters ...... 85
22. Lake 1 model verification (cont.),
0-17 meters ................... 86
23. Lake 1 model verification, 50 - 150 meters ..... 87
24. Dynamics of the phytoplankton growth rate ..... 90
25. Dynamics of the phytoplankton death rate ...... 92
26. Dynamics of the phytoplankton net production .... 93
27. Dynamics of nitrogen and phosphorus limitations . . 95
28. Biological and hypolimnetic recycling of
phosphorus ..................... 97
29. Lake Ontario Pbytoplankton Biomass,
horizontal averages 1973. IFYGL .......... 106
30. Central Lake Erie Phytoplankton Biomass,
four station log averages, 1970.
Project Hypo ....................
31. Onondaga Lake Phytoplankton Biomass
two stations, 1969
32. St. Margaret's Bay Phytoplankton Biomass,
station A, 1969 ....". ............. I10
33. Eigenvalue equation solutions for constant
growth layer ard exponentially decreasing
growth rate cases .................
34. Settling velocity-growth rate interaction.
The effect of dispersion. Exponentially
decreasing growth rate case ............
35. Population growth rate versus maximum growth
rate. The effect of Euphotic zone. The
exponentially decreasing growth rate case ..... 119
36. Population grovth rate versus settling
velocity - the effect of Euphotic zone depth.
The exponentially decreasing growth rate case . . . 121
37. Population growth rate versus dispersion
coefficient and the effect of Euphotic zone depth.
The exponentially decreasing growth rate case . . . 122
38. The effect of negative sinking velocity.
The exponentially decreasing growth rate case . . . 123
39. Application to Lakes and Oceans - the range
of biologically relevant dimensionless numbers . . . 125
40. Lake 3 segmentation ................ 132
VII
-------
41. Assumed summer circulation for Lake 3 model 136
42. Assumed winter circulation for Lake 3 model 138
43. Assumed horizontal and vertical dispersion
coefficients for Lake 3 model 139
Q
44. Regional division of Lake Ontario by Lee 145
45. Comparison of computed temperatures from
Lake 3 model with data from Lee", winter, fall . . . 146
46. Comparison of computed temperatures from
Lake 3 model with data from Lee , summer, spring . . 147
47. Principal IFYGL stations and the Lake 3 grid .... 150
48. Preliminary results of Lake 3,0-4 meters 151
49. Preliminary comparison of Lake 3 output to
observed data at two segments for June 153
50. Data output flow diagram used fro verification
and display purposes 154
51. Three dimensional plot of phytoplankton
chlorophyll calculated from Lake 3 model-June,
0-4 meters 156
A-l Variation of volume with depth 162
A-2 Variation of area with depth 162
A-3 Temporal variation of vertical dispersion
coefficient - interface 1-2 164
A-4 Temporal variation of water temperature 164
A-5 Temporal variation of solar radiation 165
A-6 Temporal variation of photoperiod 166
A-7 Determination of background water clarity - Ke . . . 168
Vlll
-------
ACKNOWLEDGMENTS
The cooperation of many IFYGL investigators who generously
supplied their data for use in this report is acknowledged.
Special thanks are due to EPA personnel at the Grosse lie
Laboratory including Messrs. William Richardson, Nelson Thomas,
and Dr. Tudor Davies. Finally, the excellent cooperation
and assistance of the Canada Centre for Inland Waters is also
gratefully acknowledged.
IX
-------
SECTION I
CONCLUSIONS
The dynamic behavior of lake-wide phytoplankton and associated
nutrients can be adequately modeled with a simplified model.
Such a model is horizontally well-mixed, but is vertically
stratified and therefore represents an average open lake condition
for Lake Ontario.
The phytoplankton modeling structure used in previous estuarine
and delta applications is therefore also appropriate for de-
scribing phytoplankton biomass in lakes.
Analytical studies of the one dimensional vertical phytoplankton
biomass equation indicate that when the biomass is increasing,
then asymptotically the log growth rate of the population will
approach a value independent of depths. Data for four bodies of
water including Lake Ontario support this theoretical prediction.
Detailed verification analyses using four years of data on Lake
Ontario indicate that the spring growth phase and peak phytoplankton
biomass are primarily controlled by increasing light and temperature
and phosphorus limitation. The mid-summer minimum in phytoplankton
is estimated to be due primarily to zooplankton grazing and nitrogen
limitation. The broad fall peak in phytoplankton is a complex inter-
action of nutrient regeneration (up to five times the external
nutrient inputs), subsequent nutrient limitation and then the fall
overturn. Both nitrogen and phosphorus are important nutrients in
this dynamic succession.
-------
It is estimated that a period of 10-20 years may be necessary
for the whole of Lake Ontario to reach a new level of dynamic
equilibrium in phytoplankton biomass after a nutrient reduction
program is instituted.
It is also concluded that a sufficient verification base for the
lake vride model has been established to permit the use of the
model for preliminary simulations of the effects of various levels
of nutrient reduction.
-------
SECTION II
RECOMMENDATIONS
Primary emphasis in this report is on the scientific foundation
and validity of the model, principally the LAKE 1 model. It is
recommended that simulations and projections of the effects of
nutrient reductions be carried out, first, in the Lake 1 model
and then explored at a greater level of detail in the 3-dimensional
Lake 3 model. These projections should be long term and attempt
to delineate the time required to reach a new dynamic equilibrium
in Lake Ontario.
A detailed examination of the behavior of the Lake 3 model is
warranted especially in view of the importance of near shore versus
open lake problems. The three dimensional nature of this model
requires a significantly greater effort for verification purposes
and for understanding the response of the Lake 3 model under
different environmental conditions of flow and near shore temperature
and waste discharges.
Following the analysis of Lake 3 model responses, it is recommended
that simulations be made of nutrient reduction programs. Such
simulations should be carried out to determine the expected phyto-
plankton biomass and the time to reach that level for the different
regions of Lake Ontario.
It is also recommended that the modeling structure of Lake 1 be ex-
panded to include phytoplankton species dependence and a more
flexible food web arrangement for the higher trophic levels. The
primary data source would be the data collected as part of the IFYGL
effort.
-------
SECTION III
INTRODUCTION
PURPOSE OF RESEARCH
The overall purpose of this research is to structure a mathe-
matical modeling framework of the major features of eutro-
phication in large lakes. Lake Ontario, the subject of in-
tensive field work as part of the International Field Year for
the Great Lakes (IFYGL) is used as the problem setting (Fig.l).
The overall objectives of the research include:
a. determination of important interactions in
lake eutrophication;
b. analysis of lake water quality and biological
responses to natural and man-made inputs;
c. formation of a basis for estimating the direction
of change to be expected under remedial
environmental control actions.
The problems of impairment of the quality of lake systems are
magnified for "large lakes" such as the Great Lakes. The size of
these lakes is such as to preclude any immediate improvement in
quality after control actions are taken. Further, it is much
more difficult to obtain reliable data on water quality, biological
structure and hydrodynamic circulation, again because of the
difficulty of sampling large lake systems. The general circulation
itself may not be adequately known in relation to climatological
and hydrological factors. In the biological area, measures of
phytoplankton populations are temporally and spatially dependent
-------
Figure 1. General basin map of Lake Ontario
-------
and reflect complex interactions with nutrients, such as
nitrogen and phosphorus,and upper trophic level predation.
This report reviews work accomplished to date on development
of a mathematical model of phytoplankton in Lake Ontario
that hopefully will form a reasonable basis for estimating
the effects of nutrient removal programs.
SCOPE OF RESEARCH
The scope of this research is therefore lake wide and
attention is directed primarily to the behavior of the phyto-
plankton for the lake as a whole. The smallest scale con-
sidered in the more advanced model (Lake 3) is on the order of
10-40 km. Detailed near-shore behavior on a scale less than
this is not a part of this research. Furthermore, this research
is concerned with phytoplankton dynamics as described on a time
scale of weeks to months. Therefore the seasonal progression of
the phytoplankton throughout a year and from year to year is the
time scale of interest. Short term fluctuations, in the order
of hours or days are not described herein.
Finally, the modeling structure developed as part of this re-
search is aimed at the phytoplankton biomass as a measure of
eutrophication and associated water quality. Attached algae,
such as cladophora, are not modeled. Particular emphasis is
placed on the interaction of the passive phytoplankton biomass
(as characterized by the concentration of chlorophyll) with
the nutrients, principally nitrogen and phosphorus and the
grazing by herbivorous zooplankton and higher order trophic
levels.
-------
The primary thrust of the research is to provide a more
rational and scientifically defensible structure for en-
vironmental decision making on the efficacy of nutrient
control and removal programs. As such, there are various
levels of detail that could be explored. The goal is to
include the "right" level of detail that explains the
observed data, provides a basis for prediction but is not
too inordinately complex.
This report summarizes work performed to date on model de-
velopment and strategy, data reduction and compilation,
sensitivity analyses and verification analyses. Sim-
ulations of conditions under varying levels of nutrient
reduction are not included in this report and will be part
of a separate document.
-------
SECTION TV
LAKE ONTARIO - BACKGROUND
With the Niagara River, the outflow of Lake Erie and hence the
cumulator of all the upper Great Lakes outflow, as its princi-
ple source of inflow and the St. Lawrence River its outflow,
Lake Ontario is the last in the chain of Laurentian Great Lakes.
Lake Ontario's narrow and deep rock basin was formed by
the action of glacial corrasion.
MORPHOMETRY
Lake Ontario is the smallest of the Great Lakes in terms of
2 2
surface area, 19,477km (7,520mi ), with a drainage area of
90,132km2 (34,880mi ).2 The mean depth of Lake Ontario, 90
meters (296 feet), is second only to Lake Superior for the Great
Lakes and can be considered one of its most important physical
23 33
features. ' Lake Ontario's volume is 1,669km (401mi ), with
3 4
a maximum depth of 244 meters (801 feet) ' . The lake's elevation
4
above sea level is 74.01 meters (242.8 feet) which makes its
depth of cryptodepression to be 170 meters (558 feet). Lake
Ontario's length is 307km (191 miles), with a maximum width of
4
87km (54 miles) . The lake's shoreline length is 1,380 km(858
miles) . The variation of volume and area with depth is illus-
trated in Appendix A (Figs. A-l and A-2).
HYDROLOGY
The Niagara River is the major source of inflow for Lake Ontario,
being the cumulative outflow of the other Great Lakes. The average
flow of the Niagara River is 195,000 cfs , and accounts for 84
percent of the flow discharged via the St. Lawrence River, as can
be seen in Table 1. The average annual precipitation on Lake
Ontario's water surface is 83.28cm (32.79 inches) and the average
annual evaporation is 71cm (28 inches) .
-------
TABLE 1 - FLOW BUDGET
Source
Major Tributaries
Niagara River
Twelve Mile Creek
Oswego River
Trent River Region
Black River
Genes se River
We Hand Canal
Other Tributaries
Other Sources
(Precipitation -
Evaporation and
Other)
Discharge
St. Lawrence River
Flow (cfs)
195,000
6,400
6,200
4,200
3,800
2,700
1,200
5,900
6,600
232,000
% of
Discharge
84.1
2.8
2.7
1.8
1.6
1.2
0.5
2.5
2.8
-------
The thermal bar development which precedes the development of the
full thermocline begins in late April or early May in Lake Ontario .
The offshore movement of the thermal bar is such that within 17 to
28 days the nearshore ring of stratified water covers over half the
area of the lake. The average depth of the thermocline is 17 meters
(55.8 feet), with its dissipation beginning in late September .
The hydraulic detention time of Lake Ontario (volume divided by flow)
is 8.1 years, and is significant in that it gives some indication of
the response time of the lake. (See Section VI)
NUTRIENT INPUTS
Nitrogen and phosphorus are the nutrients used in the present con-
figuration of the model, though silica will most likely be incorporated
in the future. The only nutrient loading information available at the
beginning of this project was for total nitrogen and phosphorus.
Loading information for the individual forms of the nutrients (ammonia,
orthophosphate, etc.) similar to those used in the model's nutrient
systems, would be of great utility.
Nutrient loads for Lake Ontario are shown in Table 2. These reflect
the total loads from the various sources (industrial, municipal and
tributary) and were obtained by the International Joint Commission
(IJC). These loads are broken down to individual waste discharges in
the IJC's excellent 1969 three volume publication.5 Table 2 in-
dicates that the Niagara River accounts for over one-half of Lake
Ontario's nutrient load. Table 3 gives the loading information as it
was incorporated into the Lake 1 model. It can be noted that the
loads are also shown as equivalent boundary concentrations, for an
10
-------
TABLE 2 - NUTRIENT LOADINGS
SOURCE
NITROGEN
PHOSPHORUS
(Niagara
Tributaries
Municipal
Industrial
Total
Ibs/Day
522,200
191,000
72,600
97,300
883,100
kg/Day
236,900
86,600
32,900
44,100
% ! Ibs/Day kg/Day
59 42,200 19,100
22 15,600 7,100
; 8 16,200 7,300
11 1,000 500
400,500 75,000 34,000
1 1 I ,_,„ _
%
56
21
22
1
TABLE 3 - NUTRIENT LOAD DISTRIBUTION
System
Nitrogen
Non-Living organic
N
Ammonia N
Nitrate N
Phosphorus
Non-Living organic
p
Inorganic P
Nutrient Load
Ibs/Day kg/Day
551,400
18,800
312,900
57,100
17,900
250,100
8,500
•; 141,900
25,900
; 8,100
j
i
'AS Boundary Concentration
mg/H
0
0
0
0
0
1
i
.443 !
.015
.250
.0457
.0143
11
-------
average inflow of 232,000 cfs. The inorganic forms of the
nutrients were set at the lake's winter concentrations with
the remaining load inputted in the non-living organic form.
AVAILABLE LAKE DATA
The bulk of the data used for this investigation was obtained
from three sources:
1. Limnological Data Reports, Lake Ontario,1966-1969
Canada Centre for Inland Waters (CCIW)
2. STORET, Environmental Protection Agency (EPA)
3. Report to the International Joint Commission on the
Pollution of Lake Ontario and the International
Section of the St. Lawrence River.5
These sources were supplemented with other data available in the
literature. The Limnological Data Reports (LDR) of CCIW comprise
the largest single source of Lake Ontario survey data available.
CCIW's cruises not only had adequately dense spatial grids (see
Fig. 2 for a typical grid) but also comprised good temporal
coverage for the years surveyed (see Table 4).
STORET, the Environmental Protection Agency's water quality storage
and retrieval system is the prime residence of all U.S. collected
water quality data. The pre-IFYGL STORET data base consisted mainly
of data from three 1965 synoptic cruises. The temporal coverage
by these cruises was not sufficient to be of direct use for verifi-
cation analysis. The data were therefore considered only as
supplemental to the basic data set. STORET1s utility is for storage
and analysis of IFYGL data and as a data base for tributary nutrient
concentration data. It was therefore decided to use the data con-
tained in the LDR as a verification data base for the preliminary
models and the STORET data as a verification data base for the
12
-------
LEGEND
• Limnological Station
ONTA*'0
44°
TORONTO
43<
10 0 10 20 30 40 50
10 0 10 20 30
79*
78<
77<
Figure 2. Typical CCIW cruise track, April 1969
-------
TABLE 4 -
CCIW LIMNOLOGICAL DATA REPORT SUMMARY
(Number of Cruises with adequate spatial definition)
1966
1967 1968 1969
Nitrogen:
Organic or
Kjeldahl-N
Ammonia
Nitrogen
Nitrate
Nitrogen
Phosphorus:
Total Phosphate
Reactive Phos-
phate
Temperature
Chlorophyll
3.
-
-
4
-
5
9
—
11
.
10
10
4
10
11
11
-
5
6
5
6
9
6
-
9
9
8
9
9
9
14
-------
more advanced models.
Due to the magnitude of data in the LDR, a data reduction
mechanism had to be implemented. The main criterion for the
data reduction was compatibility with model output so as to
facilitate comparison. Considerable effort was expended on
this task, with the principal product being main lake means
and standard deviation for the relevant variables at individual
sample depths. Main lake statistics were obtained, since this
is what the Lake 1 and 2 models were designed to represent.
Stations which had sounding depths of greater than 50 meters
were considered main lake stations. Data were also reduced for
depth intervals corresponding to segment depth for some years.
Main lake, near shore or entire lake reductions can be computed
using the present software. Also latitude-longitude blocks can
be specified to give more refined spatial definition to the data
retrievals. Therefore Lake 3 comparison data is readily obtained
with this option. Figure 3 is a sample output of the data re-
duction.
The Lake 1 model was designed to aid in the understanding of the
principal mechanisms which underlie the process of eutrophication.
A key criterion in evaluating the validity of a model is to compare
the models output to the data observed in the field. Analysis of
the field data of the water quality parameters under consideration
revealed temporal patterns that were observed year after year.
These patterns can be seen in the summary data plots used for the
evaluation of model output - (See Section VI). Therefore it is
apparent that data accumulation and reduction are key steps in the
modeling framework, for only when these are completed can the
validity of the model be determined.
15
-------
20 -
§10
Q
MEAN ± 1 STANDARD DEVIATION
\
\
J 30 F60M90A120M150J180J210A240S270°300N330D360
Figure 3. Typical output from data reduction - main
lake average @ 1 meter
16
-------
REFERENCES
1. Hutchinson, G., A Treatise on Limnology, Vol. 1.
John Wiley & Sons Inc., N.Y. 1957.
2. Beeton, A.M. and D.C. Chandler, "The St. Lawrence Great
Lakes," Limnology in North America (D.C. Frey, ed.).
University of Wisconsin Press, Milwaukee. 1963.
3. Canada Center for Inland Waters, Descriptive Limnology
Section. Personal communication (see Appendix A).
4. Canadian Hydrographic Service. Bathymetric Chart (881)
Lake Ontario. Department of Energy, Mines & Resources,
Ottawa. 1970.
5. International Lake Erie & International Lake Ontario -
St. Lawrence River Water Pollution Control Boards.
Report to the International Joint Commission on the
Pollution of Lake Ontario and the International Section
of the St. Lawrence River, Vol. 3. 1969.
6. Great Lakes Basin Commission. "Appendix 11, Levels & Flows,"
Great Lakes Framework Study, Draft No. 3, Ann Arbor. 1971,
7. Canada Centre for Inland Waters. Limnological Data Reports,
Lake Ontario. 1966-1969. Canadian Oceanographic Data
Centre, Burlington, Ontario.
17
-------
SECTION V
LAKE ONTARIO MODEL STRUCTURE
THEORETICAL BACKGROUND
The basic theory used for the Lake Ontario model makes use of
12
previous model structures (DiToro, et.al Hydroscience,Inc.,
Thomann, et.al., ) and essentially represents a mass balance
around each ecological or nutrient compartment and around physical
space.
The basic structure consists of three parts:
1. transport and dispersion sub-system
2. biological sub-system
3. chemical sub-system
The latter two sub-systems represent the kinetic behavior or the
extent and complexity of interactions between relevant variables
while the first system represents the physical processes of water
movement and mixing. The theory of each of the sub-systems is
presently developed to a different degree with the transport and
dispersion theory developed to a more advanced degree than the
biological and chemical systems.
The general equation which results from applying the principle of
the conservation of mass is
E 3s
TJ1 *\ & f? *\ o
3 ^7 ^ V C3 *7 ^
8y ly "3z"
(x,y,zft,s ,sv) + W. (x,y,z,t); k = l...m
^ K " K H = l..m
18
-------
where s, is the kth dependent \ariable (biological or chemical),
KL
u, v, w are the velocity components in the x, y and z directions
respectively, E , E and E are the dispersion coefficients in
each respective spatial direction, S, represents the kinetic
interactions between the k variables and W is the direct inputs
of the substance, s,. It should be noted that Eq. (l) as written,
does not include any interaction between the variables s^ ar.d the
transport and dispersion regime. Thus, the physical system is
separable, in the sense that it can be externally supplied and
once determined, it can be used for all variables. In vector form,
Eq. (1) is:
9s, (V. [E] (Vs.) - V . Us ) + S, + W
where [E] is a diagonal matrix of dispersion coefficients,
is a column vector, U is the velocity vector and
In Eq. (2), the first term in parenthesis on the right hand side
represents the dispersive and advective field in three dimensions,
as given by the general circulation. In this work, it is assumed
that the circulation is "known", either through observation or from
output from a hydrodynamic model. The inputted circulation can be
validated by utilizing Eq. (2) for a tracer such as dye or water
temperature to compare computed results to observed data. Water
temperature has been used in some of this work and the results of
that analysis are given in section VIII. Primary interest therefore
19
-------
centers about the reaction kinetics, their functional forms and
numerical values.
A finite difference approximation to Eq. (2) is necessary in
order to apply the equation to an actual water body. Thus, if
the lake is divided into a series of finite completely mixed
volumes, n in number and each with volume V., then it can be
shown that
d(s)k
[VI __ = [A] (s)k ± (S)k +
(3)
where [V] is an n x n diagonal matrix of volumes, (si is an n x 1
vector of the variable s, , [A] is an n x n matrix of advection and
dispersion terms; (S), is an n x 1 vector of kinetic interactions
~ JC
and (W), is an n x 1 vector of inputs of variable s, . For m
interacting variables and n physical locations, Eq. (3) indicates
that a total of mxn equations must be solved. In general, the
equations are both linear and non-linear.
When one recognizes that initially up to 10-15 interacting variables
can be specified and that 50-100 initial spatial segments may be
required, it is obvious that the computer program and data prep-
aration necessary to obtain solutions are of a substantial magni-
tude. Accordingly, prudent practice would dictate a modeling
strategy that allows one to proceed sequentially from simple models
to the more complex. Such a strategy is utilized in this work and
is illustrated in Fig. 4.
20
-------
SPATIAL DEFINITION
SYSTEM
GEOMETRY
PREPARATION
Flow
Depth, volume
Dispersion
Inputs
WATER TEMPERATURE
AND
CHLORIDE MO DEL
(67 SEGMENTS)
TEMPORAL AND KINETIC DEFINITION
Solar radiation
Temperature
Flow
Nutrients
Solar radiation
Temperature
Flow
Nutrients
LAKE 1 MODEL
(3 VERTICAL LAYERS)
LAKE 2 MODEL
(7 VERTICAL LAYERS)
Phytoplankton and
zooplankton dynamics
Temperature verification
phytoplankton,
chemistry and sediment
interactions
SYNTHESIS OF KINETIC AND SPATIAL MODELS
LAKE 3 MODEL
(67 SEGMENTS,
10-15 VARIABLES]
Coarse grid model
5000
I MODEL
i. ________ _j
(Future expansion)
Figure 4. Eutrophication modeling strategy
-------
As indicated, two parallel paths are being followed. The upper
path involves the gathering of data on the lake geomorphology
and transport and dispersion structure. The model used is a
single variable model of water temperature or chlorides both
representing tracer variables that can be used to validate
various sectors of the inputted circulation. A three dimensional
grid of 67 segments is employed.
The emphasis in the lower parallel path is on the sub-system
kinetics with spatial detail held to a minimum. The lake is
assumed to be horizontally completely mixed and is divided into a
series of vertical layers. Attention is specifically directed to
the interactions between phytoplankton, zooplankton and water
chemistry. Two models, Lake 1 and Lake 2 are explored in this
path.
The synthesis of the spatial definition of transoort and mixing with
the kinetic interactions is accomplished in Lake 3. This model, a
coarse grid model, is three-dimensional and incorporates the major
kinetic features of the Lake 1 and 2 models. with 67 segments and
10 variables, a total of 670 equations or "compartments" must be
solved. Future expansion calls for an approximate order of magnitude
increase in the number of compartments to 5000.
In this report, primary emphasis is placed on the Lake 1 model. As
such, a detailed review of the kinetics employed in the nodels is
given below for the biological and chemical sub-systems. The
systems diagram for the Lake 1 model is shown in Figure 5. The
Lake 2 model includes those variables for Lake 1 and also includes
i
a representation of the carbon cycle.
PHYTOPLANKTON CHLOROPHYLL
The basic kinetic interactions, for phytoplankton chlorophyll, ( a
measure of biomass) is given for a single volume, j, as
sP ' vi (GP - V: P: - vi" PJ + vi - Pi
22
-------
to
CO
UPPER TROPHIC
LEVEL *Z
CARBON
UPPER TROPHIC
LEVEL*!
CARBON
CARNIVOROUS
ZOOPLANKTON
CARBON
HERBIVOROUS
ZOOPLANKTON
CARBON
PHYTOPLANKTON
CHLOROPHYLL
BIOLOGICAL
SUB-MODEL
NITROGEN CYCLE
ORGANIC
NITROGEN
AMMONIA
NITROGEN
NITRATE
NITROGEN
PHOSPHOROUS CYCLE
ORGANIC
PHOSPHOROUS
AVAILABLE
PHOSPHOROUS
CHEMICAL- BIOCHEMICAL
SUB-MODEL
Figure 5. System diagram - Lake 1 model
-------
where P is the phytoplankton chlorophyll (yg/1) , G and D
growth and death rate (I/day) respectively and w.. and w. . .,
^ J c J 31 13 are the
sinking velocities of the phytoplankton between segments j and i.
Phytoplankton Growth Rate
The growth rate expression is similar to that developed previously
(1) and as used in this work is given by
GPJ= Gl(Gmax,T)- ^a/s^e/'^ ' <
-------
and k1 is the light extinction coefficient at zero chlorophyll
G
and f is the photoperiod. All pertinent variables are functions
of j , the segment location.
The nutrient effect makes use of product Michaelis - Menten
kinetics and is given by
N'P = N2 + N3
kmn
The range of values for K and K , the half-saturation constants
for nitrogen and phosphorus respectively, has been documented pre-
viously . More recent other work (5,6) tends to support the
original range; namely K in the range 5-50 yg N/l and K in the
range 1-10 yg/1. There is an obvious species dependence of this
constant as well as a possible dependence on the euthrophic state
of the lake. In the Lake 1 model, Kmn is used at 25 yg N/l and Kmp
of 1 and 10 ugP/1 is used as a range.
At the present time, silica limitation is not used in this model.
7 8
Goering et al and Paasche have estimated the half-saturation
value for silica at about .02 - .08 mg Si/1 for marine species.
Under peak growing conditions in Lake Ontario, surface silica
values are about 0.1-0.2 mg Si/1 which tends to indicate that
silica may not have a significant effect on the growth rate. Work
is in progress, however, to include this nutrient since at the
lower concentration and higher Michaelis levels, it could be
important.
The complete growth rate expression, G . is then given by the
product of (6), (7) and (8). Although there is no a priori reason
for using a product formulation, various experiments (see(l)) have
tended to support the approach. Another expression has been proposed
25
-------
on the grounds that Eq. (5) is too severe in reducing growth due
to the product term. The alternate construct is (9):
G1 = G, .
P X I/I + I/.
N' ^ *-' (9)
or in general
G1 = G. ,
P 1 - k (10)
where k is a normalization factor and y. represent the general
reduction factors. In general, from (5) and (10):
G1 > G
P P
the extent of the inequality depending on y.. Undoubtedly then,
Eq. (5) will result in increased growth at least at some times;
however, due to the coupling to the nutrient pool, such increased
growth will in turn be more severe in its impact on the nutrients.
Eq. do) therefore appears to be less severe than Eq. (5), but
only with respect to phytoplankton growth. In the final analysis,
the validity of any of these formulations rests on the ability to
reproduce observed data, barring any real theoretical reasoning.
This has been done at least for the product nutrient terms in early
work by Ketchum in 1934 (cited in (1)) and more recently by
DiToro for the flask experiments. The Lake 1 model, was modified
to test the sensitivity of the solutions to Eq. (5) or Eq. (10) and
to determine which expression seems to be more appropriate. The
results of the comparison are discussed more fully below.
It can also be noted that the gross growth rate of the phytoplankton
biomass, G .P. is equivalent to the daily average productivity in the
26
-------
jth segment. Thus,
v. = a G . P .H.
D cc PD D D
o
where vj is the average rate of carbon fixation (mgC/m -day) and
a is the carbon to chlorophyll ratio (range 20-100mgC/mg Chlor).
CC
Phytoplankton Death Rate
The expression for the death rate of the phytoplankton includes
two primary effects: endogeneous respiration and predation by
herbivorous zooplankton. The death rate is given by
Dpj = K2(1.08)T-20 f CgZ.
where ILis a constant (I/day) , C is the grazing rate of the her-
bivorous zooplankton (liters/day - mg zooplankton carbon) and Z is
the carbon concentration of the herbivorous zooplankton. As noted
below for the zooplankton observed in Lake Ontario, a filtering
rate of .03 - .06 liters/mgC-day-°C appears appropriate and at
20°C, represents about 2-6 ml filtered per individual per day.
Phytoplankton Sinking
In Eq. (4), the sinking of phytoplankton from one segment to another
segment can be a potentially significant sink or source of phyto-
plankton and associated nutrients. This phenomena is a complex
one incorporating vertical turbulence effects, the density
structure and the physiological state of different species of
phytoplankton. For example, the generation of gelatinous sheaths
by phytoplankton has been shown to be of importance in settling
(11) and apparently the settling velocity of nutrient rich cells is
somewhat less than cells that are nutritionally deficient (12).
27
-------
Published values of the sinking velocity of phytoplankton, mostly
in quiescent laboratory conditions range from 0.07 - 18 meters/day
In some instances, the settling velocity is zero or negative.
Burns and Pashley have investigated the settling velocity
profile of particulate organic carbon in Lake Ontario with an in-
genious in situ measuring device. Their results show a general
decrease of settling velocity with depth. and a marked seasonal
variation of settling velocity. For example, on July 19,1972
settling velocity in the thermocline (about 10 m) was estimated
at about - 0.3m/day, while at 50 meters depth, the settling
velocity was about + 1.3m/day. But in March 1973, the velocity
profile was all positive rangina from 0.9m/day at 20 meters to
G.lm/day at 140 meters. Overall, the values given by Burns and
Pashley varied between - 0.4 and + 2.0m /day. These results and
others tend to indicate that the settling velocitv should be
related in some way to the state of the phytoplankton biomass
(e.g. actively growing versus declining phase) and/or vertical
velocity variations. In this work however, the vertical settling
velocity has been treated as a constant. As such, this effect is
only approximately modeled and accordingly, the settling velocity
has been treated as a constant over a range from . 05m/day to
0.5m/day. The sensitivity of the results to the phytoplankton
settling velocity is discussed below.
Zooplankton Carbon
The measure of the herbivorous zooplankton population Z is
taken as the carbon content of the biomass and all sources of food
for Z are the phytoplankton chlorophyll. The general expression
for a single volume j is therefore:
- V. (Gfl) - D). Z (13)
z. 3 z z '3 3
28
-------
where G and D are the growth and death rates respectively
of Z(1) . The form of G, is
Z
Kg
where a is the zooplankton carbon to phytoplankton chlorophyll
zp -
efficiency, C is the grazing rate of the herbivorous zoo-
gzp
plankton and K is the half-saturation concentration of phyto-
g (1)
plankton at which the growth of Z v is half of maximum growth.
This latter effect reflects the limiting growth of zooplankton
at higher phytoplankton populations.
The mortality of the herbivorous zooplankton is given by two terms;
endogenous respiration and grazing by the next highest trophic
level.
Therefore :
D(1) = KtT) + CT) Z(2) (15)
z
where K is the endogenous respiration rate as a function of
Z
temperature (l/day-°C) , C 01 is the crrazing rate ( liter s/dav-mgC-
( 7 \ y^-1-
°C) and Z is the next highest trophic level, designated by
superscript (2) .
Expressions for all trophic levels above the herbivorous zooplankton
arranged in a linear food chain are similar to Eqs. (13) and (14).
S(2) =V. {GC2)_D(2) z(2)
Z 1 '• Z Z 1 j (16)
29
-------
Where r<2> a c z(1)
\~J — Q. *— ^ i " / 1 T \
z z g21 (17)
and
(91 (2) (1}
T\\ ' — Tf *• ' (TM 4- P fTM 7V '
L/ -~ rv V -L / ' ^ 11 V -L J ^
z z g ^J.
These expressions are an obvious oversimplification of a complex
food web and reflect only a "linear" type of food chain (See
Fig. 5). Indeed the effect of the zooplankton on the phytoplankton
as formulated above is to always decrease the population, yet it
has been shown (14) that zooplankton can also increase phytoplankton
population or have no effect. The mechanism for an increase in the
population is apparently through the passaae of aelatinous green
algae through the zooplankton gut with subsequent break-up of the
colonies and regrowth after excretion. The smaller species are
generally decreased in numbers by .zooplankton grazing, specifically
smaller diatoms and Asterionella formosa, which are characteristic
of Lake Ontario.
The zooplankton of Lake Ontario are dominated by the crustaceans
specifically of order copepoda and cladocera. Most hauls are
made over 0-50 meter der>th range. Principal cupepoid species
given by Patalas15 are cyclops bicuspidatus and Tropocyclops
prasinus and of the cladocerans: oaphnia retrocurva and Bosmina
longirostris. In Patalas1 work in 1967, these four species
accounted for 91% of the total zooplankton and averaged about
60/liter from June - October.
Others (16, 17, 18) generally show values in a similar range in
some cases up to a maximum of almost 200/liter. For use in biomass
computations, conversions must be made to equivalent carbon. Such
30
-------
conversion depends on age and accompanying weight of the animal
T 8
and can therefore vary over a wide range. Watson et al in
their conversions used weights of species to arrive at total
biomass and equivalently averaged about 2.8yg dry weight/ind.
The results using an average of 40% carbon by dry weight show
peak values general in August at lakewide averages of between
0.07 - 0.12mg carbon/liter.
At the average dry weight of about l-5yg/ individual grazing
1 19
rates of about 0.8 - 2mi/ind - day are appropriate ' . In-
deed the rate for Bosmina long., a dominant species in Lake
Ontario has been given as l-3ml/ind-day . Grazing rates ex-
pressed as liters per mg zooplankton carbon per day have been
assumed to be linearly related to temperature and have been used
o
in the range of .03-. 06 (1/mgC-day- C) corresponding to the
reported grazing ranges at 20°C.
NITROGEN SYSTEM
The basic equations for the nitrogen sub-system follow previous
work (2,3) and include non-living organic nitrogen, N, , ammonia
nitrogen, N^ and nitrate nitrogen, N., . The equation for the
source - sink term for organic nitrogen is:
1
1 S D Z .
z
C
. DP " - «
1
—
C
-L (T)
)C Z (18)
where a, is nitrogen to chlorophyll ratio, and a is the carbon
to chlorophyll ratio and K, is the overall decay of N- .
31
-------
The first and third terms represent organic nitrogen released
through endogenous respiration by the zooplankton and phyto-
plankton respectively and the second term represents the or-
ganic nitrogen of the grazed but unassirailated phytoplankton.
The last term represents the decomposition, settling and other
effects that contribute to the decay of organic nitrogen.
For ammonia,
SN2. = K12(T)N2 - K22 (T) N2 - al aGp P
where K is the rate of production of ammonia from organic
12
nitrogen, K_2 is the rate of oxidation of ammonia to nitrate,
and a is a preference factor given by
(20)
or is alternately set equal to one-half for an "indifferent"
uptake of nitrogen. The last term represents the uptake of
ammonia for phytoplankton growth. For nitrate,
V " = K« N2 * 'l a"°" °P P
3 ' 3
where K00 = K0- the production of nitrate by ammonia oxidation.
23 &•*• i
PHOSPHORUS SYSTEM
The equations for the phosphorus cycle used in the models are
as previously developed (2,3) and include non-living organic phos
phorus and phosphorus available for phytoplankton growth. For or
ganic phosphorus, the equation for the source and sink for
is:
32
-------
where a is the phosphorus to chlorophyll ratio and K.
represents a general decay of organic forms.
For the inorganic phosphorus, assumed as phosphorus availa-
ble to the phytoplankton, the source-sink equation for p2 is:
S _ . =K10 (T) p - a G P - K90 (T) P9
P2/J 12 1 P P 22 2 (23)
where K,„ represents the rate of decomposition of organic phos-
phorus to forms available for phytoplankton utilization and K2«
represents an overall loss of inorganic phosphorus.
33
-------
CHEMICAL MODELS - DISSOLVED CARBON DIOXIDE
The inclusion of chemical reactions into models of Great Lakes
water quality, by contrast to the biochemical and biological re-
actions considered above, poses a new set of problems by virtue
of the rapid rates and the highly non-linear forms of the govern-
ing reactions. Nevertheless, it appears that the problems can be
handled nicely within the context of conservation of mass equations.
The investigation of purely chemical equilibrium models has been
underway for quite some time. An altogether satisfactory pre-
sentation of these results is available. Application of certain
of these models to Great Lakes waters chemistry has been undertaken
(21 22)
by Kramer in a series of important papers. ' From this work
it is clear that the biological reactions affect the chemical systems
primarily via the assimilation of dissolved carbon dioxide during
phytoplankton growth. From the point of view of the eutrophication
problem,at issue are the purely chemical reactions which phosphorus
undergoes; for example precipitation or dissolution of phosphate
minerals. These reactions may be of practical importance in terms
of the overall dynamics of the phosphorus in the lakes. The other
chemically active phytoplankton nutrient, silica, may also be affected
by purely chemical reactions involving the alumincH=silicate
minerals. Thus there are possible practical consequences of these
purely chemical reactions so that their inclusion, at least on a
preliminary basis, is warranted.
The presentation given below is based on the carbon dioxide -
bicarbonate - carbonate equilibria primarily because of its importance
as the major buffer system of the lakes, and because it can serve
34
-------
as a prototype for all other such reactions. In fact the procedure
is the same for the inclusion of all such reactions within the con-
text of conservation of mass equations.
Equilibrium Analysis
The carbonate equilibria is treated in the conventional way* '
The major species considered are aqueous or dissolved [CO-] (the
concentration of hydrated CO- being small), bicarbonate, [HCO ~],
and carbonate, [C0.,=] , ion, together with the hydrogen, [H+] , and
hydroxyl, [OH~], ions. At first glance one is tempted to write mass
balance equations for each of these species. However this is compli-
cated by the fact that they undergo reversible ionization reactions,
e.g. CO- + H_O = H+ + HCO.," which occur verv rapidly relative to
(23)
the time scales of mixing and gas transfer. In fact it is clear
that the individual species are each extremely reactive so that a
direct mass balance formulation results in equations which are non-
linear and numerically quite badly behaved.
The very fact that the ionization reactions are rapid leads to a
much more tractable formulation in terms of quantities which are con-
servative relative to the ionization reactions. A formal mathe-
matical discussion is given subsequently. However, it is easy to
see the principle involved. Consider the total inorganic carbon con-
centration.
CT = [C02] -I- [HC03~ ] + [C03=] (24)
It is clear that any reversible ionization reaction will not affect
the CT of the system (we assume no carbonate precipitates are forming
or dissolving.) Thus a mass balance equation for C can ignore the
35
-------
ionization reactions since they are neither sources nor sinks of
CT-
The electroneutrality requirement can be used to obtain a second
conservative quantity. Let {M) and {L} be the concentration in
equivalents of the metal ions and liqands, e.g. N , K+, C ++,....
a a
and Cl~, SO.= ..., which are present and assumed not to react
appreciably x^ith the carbonate or hydrogen-hydroxide species. Then
it is clear that a charge balance equation:
{M} - {L} = [HC03~] + 2[C03=] + [OH~] - [H+] (25)
is unaffected by ionization reactions so that the Quantity called
alkalinity
[Alk] = {M} - {L} (26)
is also conservative in the sense that a mass balance equation for
Alk need not consider ionization sources or sinks. It can be shown
there are no other independent conservative auantities for this
system, the other common conservative quantities 20 are linear
combinations of CT and A.
If C_ and Alk are known, the concentration of all the species being
considered is easily obtained. In fact, using the equilibrium
(20,24)
equations
[H+] [OH~]
[C03= ] = K2[HC03~]
36
-------
(24)
it is true thatv '
Alk = a(H) C + y(H) (28)
where 1 + 2K /[H+]
(29)
+ [H+]/KI + K2/[H+]
and
Y \ •* / — -"-v^.,/ L 11 j L *4. j (30)
It does not seem to have been noticed that if the alkalinity con-
centration is large relative to y (H) = [OH~] - [H+] , which is the
common situation then, eq. (28) becomes
= a(H) (31)
so that only the ratio of alkalinity to total inorganic carbon is
important in determining the resulting pH whenever Alk ^ Alk 4-
[H+] - [OH~] • Using tabulated values of K, and K~ for pure water
Figure 6 illustrates the pH, Alk/CT relationship as a function of
temperature. These figures are most useful since they represent
the carbonate equilibria completely as a function of the ratio of the
cricital conservative quantities. Taken together with the conventional
log concentration plot of the species versus pH, they provide a
complete picture of the behavior of the equilibria.
37
-------
5 15 25
1020 30° C
1.20
[ALK]/[CT]
Figure 6. Relationship of pH and the ratio,
38
-------
The equations which govern the distribution of Alk and CT are
available as mass balances, independent of the ionization re-
actions. For example, in a one dimensional vertical model the
equations would be
3 [Alk] ~ 1 _1(EA 8 [Alk]. s
3t A 3z(bA ~g^ > = SAlk
3t A 3z" ?z~ CT
with a zero mass flux boundary condition at z=o for alkalinity and
an air-water exchange condition for C :
T = KL (C02(g) - CO
representating the exchange of carbon dioxide. The difficulty in
solving these equations is due to this non-linear boundary condition
since C02 (aq) is a non-linear function of the dependent variables
of the problem: [Alk] and [C_], so that a numerical procedure is
required.
The internal sources and sinks of alkalinity and inorganic carbon
are due to biological reactions. For example, the nitrification re-
action (if bacterial synthesis is ignored) can be represented as:
NH^ + 1 1/2 02 -»• NO + H20 + 2H+
so that as ammonia is oxidized to nitrate, alkalinity is consumed
(i.e. H is produced) and a term representing this sink of
alkalinity would be included in the expression of SA,,.
39
-------
Similarly as phytoplankton grow they utilize inorganic carbon and
this sink would be included in S_ . These mass balance
CT
equations would then be solved in conjunction with the equations
for the other species of interest.
Linear Approximation
It is clear that mass balance equations which ignore the ionization
reactions are only correct if the dependent variable is conservative
relative to these reactions . Consider the conservative quantity
[C ]-[Alk]=[CO- - Acidity]. Clearly it is conservative since it is
the difference of conservative quantities. From eqs. (24) , (25)
and (26) , it is the case that
[C02-Acy] = [C02] + [H+] - [C03=] - [OH~] (36)
But this equation implies that over a certain range of pH and C ,
the major component of CO, -acidity is [C00] itself. Figure 7 pre-
^ £
sents the % error in the approximation [CO2-Acy] = [CO™] . The
chemical reasoning which leads to the same conclusion begins with the
observation that in the above pH range, the alkalinity is primarily
bicarbonate ion. If a source introduces only carbon dioxide, this
does not change the alkalinity; therefore the concentration of
bicarbonate remains constant. Hence the introduced CO- does not
ionize to HCO-~ and H+ but rather remains as CO2. Hence ionization
reactions can be ignored. Thus from eqs. (32) and (33) the mass
balance equation for CO2~Acidity is
9 [CO -Acy] 1 3 8 [CO -Acy] _ S _ s ,„.
•* _ - __ (EA _ I _ ) ~ T '
at A 9z az
40
-------
a:
o
UJ
O
CE
20
10
0
-10
-20
-30
Cr=10-5M/l
10
-'
5.0
6.0 ,, 7.0
pH
8.0
C02(approx) -C02
error: 100
C02
Error in approximation:
[C02(aq)] = [C02 -acidity]
Figure 7. Per cent error in assuming
[C02 - AcyJ * [C021
41
-------
and the reaeration source of [CCu-Acy] can be approximated as
KL ([COj] . - [C02 - Acy]) so that within the errors in Fig.
7 the equation for [CO^-Acidity] is linear. And since the
alkalinity equation is also linear, simple analytical solutions
are available and superposition can be applied to buildup a
solution when many sources and sinks are active. This is in con-
trast to the non-linear C_ equation which must be solved with
all sources and sinks at once, a situation which is quite inconvenient.
Unfortunately this linear region does not include the pH range
typically encountered in the Great Lakes so that this simplifi-
cation is not available. However, for more acid situations this
simplification provides a useful and direct avenue for the solution
of the inorganic carbon - alkalinity mass balance equation.
Temperature Dependence of the Equilibrium Constants
The equilibrium constants of the dissolved carbon system are rather
strong functions of temperature so that this effect must be taken
into account. The implimentation of this dependency poses no
serious problem since standard curve fitting techniques can be
applied to the experimental data. However, such an approach lacks
a certain elegance and generality. It turns out, however, that for
the temperature ranges encountered in natural waters there is an
entirely satisfactory solution available based on the thermodynamic
constants of the reactions.
It is well known that the equilibrium constant, K, of a reversible
reaction is related to the change in Gibbs free energy of the re-
action at standard state, AG° via the equation:
Ar*°
In K = - r (38)
RT
42
-------
where
AG° = £ AG| - £ AG|
products reactants
AG£ = the Gibbs free energy of
formation
R = universal gas constant
T = temperature °K
and further that the free energy change is related to the change
in enthalpy, AH° and entropy, AS° as follows:
AG° = AH; - T AS; (39)
Hence the equilibrium constant is given by:
in K = -^ + ^
RT R (40)
This equation is strictly true at any temperature; the difficulty
is that AH° and AS0 are also functions of temperature. The con-
ventional temperature at which they are available is 298.15 °K or
25°C. Table 5 presents the relevant thermodyiiamic properties for
the carbonate system and Table 6 gives the properties for the re-
actions of interest.
Two approximations are investigated:
(1) Assume that AH° and AS0 are constants over the
temperature range of interest (20)
(2) Apply a correction which attempts to evaluate
the change in the specific heat as a function of
temperature
43
-------
TABLE 5 -
THERMODYNAMIC PROPERTIES AT 25°C
Compound
C02(g)
H2C03(aq)*
HC03
co3=
H+
H2°
OH~
AHf
(kcal/mole)
-94.051
-167.22
-165.39
-161.84
0
-68.315
-54.970
AGJ
(kcal/mole)
-94.254
-148.94
-140.26
-126.17
0
-56.687
-37.594
ASJ
(cal/deg mole)
51.06
44.8
21.8
-13.6
0
16.71
-2.57
Source: Wagman, 1968 ^25^
* H2C03(aq) =-H2C03 + C02(aq)
* [C09(aq)l
TABLE 6
THERMODYNAMIC PROPERTIES AT 25°C
Reaction
1. H20 = H+ + OH~
2. C02 (g) + H20=H2C03
3. H2C03 = H+ + HC03
4. HCO3 = H+ -1- C03
"r
(kcal/mole)
13.345
-4.854
1.830
3.55
(cal/deg mole)
-19.28
-22.97
-23.0
-35.4
44
-------
The equation for the second approximation is
r + T(T) ASr] (41)
where
T(T) = Tr- JL [1_e(T-Tr)/0]
(42)
and Tr = 298.15 °K
0 = 219.
w = 1.00322
These approximations are compared to the experimental equilibrium
constants in Table 7. As can be seen the approximations are
quite good with the non-linear relationship giving somewhat better
results at the lower temperatures. Thus to within a few percent
these approximations are quite adequate.
Simple Equilibrium Models - Exploratory Calculations;
In order to explore the computational feasibility of including
chemical reactions in water quality models a series of computations
have been performed utilizing the Rand Chemical Composition Program
It is expected that this computer program will form the basis for
the chemical submodel calculations to be included in Great Lakes
water quality models.
The computations are for the system C0« - H~0 for which the partial
pressure of CO- (g) is fixed and varying amounts of acidity and
alkalinity are added. The results are computed for the temperature
o
range 0 to 30°C using the linear approximation for the temperature
45
-------
TABLE 7
TEMPERATURE DEPENDENCE COMPARISON
REACTION -
HC03-
ENTHALPY OF REACTION
ENTROPY OF KEACTION=
H+ + C03=
= 3550.0000 CAL/MOLt
-35.4000 CAL/DEG-MULE
TEMP CC) EXPERIMENTAL METHOD 1PER CENT
(20)
PK
PK
ERROR
METHOD 2 PER CENT
PK ERROR
0.0
5.0
10.0
15.0
20.0
25.0
30.0
REACT
10.6250
10.5570
10.4900
10.4300
10.3770
10.3290
10.2900
ION - H2C03*
10.5770
10.5200
10.4767
10.4292
10.3832
10.3369
I0.29b9
-.45
-.29
-.13
-.01
.06
.10
.06
10.6181
10.5523
10.4917
10.4360
10.3852
10.3389
10.2970
-.07
-.04
.02
.06
.08
.10
.07
= HC03- * H*
ENTHALPY OF REACTION = 1030.0000
ENTROPY OF REACTIONS -23.0000
TEMP CC)
0.0
5.0
10.0
15.0
20.0
25.0
30.0
EXPERIMENTAL
PK
6.5790
6.5170
6.4640
6.4190
6.3810
6.3520
6.3270
CAL/MOLE
CAL/DEG-MOLE
(20)
METHOD 1 PER CENT
PK ERROR
6,490ti
6.4645
6.4391
6.4146
6.3909
6.36bl
6.3459
46
-1.34
-.81
-.38
-.07
.16
.25
.30
C26)
METHOD 2 PEH CENT
PK ERROR
6.5175
6.4816
6.4489
6.4191
6.3922
6.3681
6.3466
-.93
-.54
-.23
.00
.is
.25
.31
-------
TABLE 7
TEMPERATURE DEPENDENCE COMPARISON (continued)
REACTION -
H20 = H2C03*
TEMP
ENTHALPY OF REACTION =-4854.0000 CAL/MOLt
ENTROPY OF REACTIONS -22.9700 CAL/DEG-MULE
EXPERIMENTAL
PK
C20)
METHOD 1 PER CENT
PK
EKROR
METHOD 2(26PER CENT
PK ERROR
0.0
s.o
10.0
15.0
20.0
25.0
30.0
REACT
1.1100
1.1*00
1.2700
1.3200
1.4100
1.4700
1.5300
ION - H20 =
1.1364
1.2062
l.,2735
1.3385
1.4U13
1.4620
1.5207
H* + OH-
ENTHALPY OF REACTION =13345.0000
ENTROPY OF REACTIONS -19.2300
TEMP <°C)
0.0
5.0
10.0
15.0
20.0
25.0
30.0
EXPERIMENTAL
PK
14.9435
14.7338
14.5346
14.3463
14.1669
13.9965
13.8330
2.37
1.36
.28
1.40
-.62
-.54
-.61
1.1630
1.22J3
1.2832
1.3430
1.4026
1.4620
1.5214
4.78
2.80
1*04
1.74
-.53
-.54
-.56
CAL/MOLE
CAL/DEG-MOLE
. (20)
M€THOO I PER CENT
PK ERROR
14.8911
14.6991
14.6140
14.3352
14.1626
13.9Vb«
13.6344
-.35
-.24
-.14
-.08
-.03
-.01
.01
METHOD Z26
PK
14.9134
14.7135
14.5221
14.3390
14.1636
13.9958
13.6350
!)PER CENT
ERROR
-.20
-.14
-.09
-.05
-.02
-.01
.01
47
-------
dependence of AG° for the three reactions of interest. Fig. 8 (a)
gives the total inorganic carbon that is in equilibrium with an
atmosphere for which the partial pressure of CO_ (g) is 10~3-5
atmospheres. As alkalinity is added the C^ increases and it is
JL
approximately true that [C ] - [Alkj with the pH in the range
7.5 to 8.5. The temperature effects change the pH of the
resulting solutions by approximately 0.25 pH units as shown in
Fig. 8(b).
For varying partial pressures of CG>2 (g) and Alk=0, the results
are shown in Fig. (8)c and (8) d. A substantial temperature effect
occurs with the CT in solution varying by a factor of two over the
range of interest, whereas the pH variation is smaller, due to its
logarithmic units.
Computation time for the numerical solution of one set of conditions
is on the order of 0.1 sec. of cpu time on a CDC 6600. Thus for
daily evaluations of the chemistry for a year simulation for ten
layers would consume on the order of 300 seconds, which is not
excessive so that the computation is feasible.
Theoretical Considerations
A formal derivation of equations (32) and (33) requires that
ionization reactions be considered explicitly as sources and sinks
of the components. Assume that the following three reactions are
the mechanisms by which hydration and ionization take place:
C02 + H20 = H+ + HC03~ ; r^ = k^O^tCC^] - [H+] [HC03~])
HC03~ = H+ + C03= ; r2 = kv2 (K2 [HC03~] - [H+] [CO3=] ) (43)
= H+ + OH~ ; r = k ( -[H+] [OtT]
48
-------
2.8
2.4
_2.0
iU
ol.2
0.8
0.4
0.0
8.0
7.0
6.0
L5.0
4.0
3.0
2.0
0.0
CT vs TEMPERATURE AT
VARY ING ALKALINITY
2.0
1.5
1.0
0.5
0.0
*
0.14
0.12
_0.10
1 0.08
-------
and that the law of mass action correctly describes the rates of
reaction as indicated. For notational convenience let D denote the
mass transport differential operator, i.e. the lefthand side of
eqs. (32) and (33) are denoted by D[Alk] and D[C ]. The con-
servation of mass equations for each species is obtained by
equating mass transport to sources and sinks:
D[CO ] = -r + S
Z L C°2 (44)
D[H+1 = rl + r2 + r3 + SR (45)
D[HC03-] = rx - r2 + S (46)
D[C03=] = r2 +5 (47)
D[OH"1 = r3 +SOH (48)
The S's are whatever direct sources of the species exist, and the r's
are the rates associated with the assumed reactions. It is, of course,
possible in principle to solve these simultaneous nonlinear differential
equations numerically. The values chosen for the k 's would be large
enough so that the solution is unaffected by their magnitude. This
procedure has been tried and found to be quite unstable numerically,
the cause being the "stiffness" of the equations due to the large
reaction rates which makes convergence difficult to achieve. Also the
k 's are, in a sense, artifacts of the formulation since they are
present only to introduce the ionizations correctly in a formulation
that requires such reactions to be included as kinetic expressions.
Therefore it would be convenient to obtain equations which are devoid
of the reaction rates entirely . This can be done by defining new
50
-------
variables which are appropriate sums of the species in such a
way that the reaction terms cancel out. Thus if C = [CO^] +
[HC03 ] + [C03=] then by adding eqs. (44), (46) and (47) the
result is
DCT = SC02 + SHC03 + SC03 (49)
The other variable is, of course,
Alk = [HC03~] + 2 [C03=] + [OH~] - [H+] (50)
and since D is a linear operator,
DtAlk] = D[HC03~] + 2D [CO3=] + D [OH] - D [H+] (51)
and using eqs. (45), (46), (47), and (48) the result is eg.(32) as
before, with the r's cancelling out in a satisfying way.
So far no approximations have been made; the equations for Alk and
C are exact, regardless of the magnitude of the k 's. They are
(2R^ v
the invariants of the reactions °' in the sense that their values
are unchanged by the kinetics. The approximation, and it is an
excellent one for this application, is to compute the concentrations
of the species from the equilibrium expressions, eq. (27), rather
than the differential equations (44-48)
There is a slight theoretical flaw in the above derivation since the
explicit forms are assumed for the reactions considered and a re-
action of the form.
C02 + OH~ = HCO3~
was not explicitly considered. It is a simple matter to verify that
51
-------
the inclusion of this reaction does not affect the validity of
the resulting equations for Alk and C but no guarantee has been
presented that such a reaction does not exist. This can be done
C28) although it requires some results from linear vector space
theory. The crux of the demonstration lies in showing that the
subspace spanned by the reactions implied by the equilibrium
equations is orthogonal to the subspace spanned by the invariant
quantities, Alk and C , and that the dimension of the orthogonal
subspace is two which implies that no other independent invariant
quantities exist, and that the dimension of the reaction subspace
is three, indicating that all the possible reactions have been taken
into account for the five quantities of interest
The theoretical framework and a theorem which indicates that the
(29)
invariant quantities are always easy to find has been structured.
The results are essentially a generalization of the previous deri-
vations .
52
-------
SECTION V
REFERENCES
L. DiToro, D.M., D.J. O'Connor and R.V. Thomann. "A
dynamic model of the phytoplankton population in the
Sacramento - San Joaquin Delta," Advances in Chemistry
No. 106, pp.131-180. American Chemical Society, Washington,
D.C. 1971.
2. Hydroscience, Inc. Limnological Systems Analysis of the
Great Lakes. Prepared for Great Lakes Basin Commission.
Westwood, N.J. 474 pp. 1933.
3. Thomann, R.V., D.M. DiToro and D.J. O'Connor. "Preliminary
model of Potomac Estuary phytoplankton," Journal of
Environmental Engineering Division, ASCE 100(SA3) -.699-715.
1974.
4. Eppley, R.W. "Temperature and phytoplankton growth in the
sea." Fishery Bulletin 70(4):1063-85. 1972.
5. Caperon J. and J. Meyer. "Nitrogen-limited growth of marine
phytoplankton - II. Uptake kinetics and their role in
nutrient limited growth of phytoplankton," Deep Sea
Research 19:619-32. 1972.
6. Hendrey, G.R. and E. Welch. The Effects of Nutrient
Availability and Light Intensity on the Growth Kinetics of
Natural Phytoplankton Communities. Presented at 36th Meeting
of the American Society of Limnology and Oceanography,
Salt Lake City, Utah. June 1973.
53
-------
7. Goering, J.J., D.N. Nelson and J.A. Carter. "Silicic acid
uptake by natural populations of marine phytoplankton, "Deep
Sea Research 20(9):777-89. 1973.
8. Paasche, E. "Silicon and the ecology of marine plankton
diatoms, I. Thalassiosira pseudomona (Cyclotella Nana) grown
in a chemostat with silicate as limiting nutrient," Marine
Biology 19:117-126. 1973.
9. Bloomfield, J.A., R.A. Park, D. Scavia and C.S. Zahorcak.
"Aquatic modeling in the eastern deciduous forest biome, U.S.,"
in Modeling the Eutrophication Process (E.J. Middlebrooks,
H. Falkenborg and T.E. Maloney, editors), pp.139-158. Workshop
proceedings of the International Biological Program. Utah
Water Res. Lab., Logan. 1973.
10. DiToro, D.M. An evaluation of phytoplakton kinetic behavior
in flask experiments. In preparation. 1974.
11. Hutchinson, G.E. A Treatise on Limnology. Vol. II: Introduction
to Lake Biology and the Limnoplankton. J. Wiley & Sons, N.Y.
601 pp. 1967.
12. Yentsch, C.S. "Marine plankton," in Physiology and Biochemistry
of Algae (R.A. Lwein, editor), pp.771-797. Academic Press,
N.Y. 1962.
13. Burns, N.M. and A.E. Pashley. "In situ measurement of the
settling velocity profile of particulate organic carbon in Lake
Ontario," Jour. Fish Res. Bd. of Canada 31(3):291-297. 1974.
54
-------
14. Porter, K.G. "Selective grazing and differential digestion
of algae by zooplankton," Nature 244(5412):179-180. 1973.
15. Patalas, K. "Composition and horizontal distribution of
crustacean plankton in lake Ontario," Jour. Fisheries
Res. Bd. of Canada 26(8):2135-2164. 1969.
16. Glooschenko, W., J.E. Moore and R.A. Vollenweider. "The
seasonal cycle of pheo-pigments in Lake Ontario with particular
emphasis on the role of zooplankton grazing," Limn. & Ocean.
17(4)-.597-605. 1972.
17. McNaught, B.C. and M. Buzzard. "Zooplankton production in
Lake Ontario as influenced by environmental perturbations,"
in First Annual Reports of the EPA/IFYGL Projects, NERC-ORD,
EPA-660/3-73-021, pp.28-70. U.S. Environmental Protection
Agency, National Environmental Research Center, Corvallis,
Oregon. 1973.
18. Watson, N.H.F. and G.F. Carpenter. "Seasonal abundance of
crustacean zooplankton and net plankton biomass of Lakes Huron,
Erie and Ontario," Jour. Fish. Res. Bd. of Canada 31(3):309-317.
1974.
19. Kibby, N.V. "Effects of temperature on the feeding behavior of
Daphnia Rosea," Limn. & Ocean 16(3):580-581. 1971.
55
-------
20. Stumm, W. and J.J. Morgan. Aquatic Chemistry.
Wiley-Interscience, New York. 1970.
21. Kramer, J.R. "Theoretical model for the chemical
composition of fresh water with application to the
Great Lakes," in Pub. No. 11, pp.149-160. Great
Lakes Research Division. 1964.
22. Kramer, J.R. "Equilibrium models and composition of
the Great Lakes," in Equilibrium Concepts in Natural
Water Systems, Adv. in Chemistry Series No. 67, pp.243-254.
American Chemical Society, Washington, B.C. 1967.
23. Kern, D.M. "The hydration of carbon dioxide," J. Chem
37(l):14-22. 1960.
24. Trussell, R.R. and J.F. Thomas. "A discussion of the
chemical character of water mixtures," J.AWWA 63(1):49-51.
1971.
25. Wagman, D.C., W.H. Evans, V.B. Parker, I. Halow, S.M. Bailey
and R.H. Schumm. "Selected values of chemical thermo-
dynamic properties," Technical Note 270-3. National
Bureau of Standards. Jan. 1968.
26. Helgeson, H.C. "Thermodynamics of complex dissociation in
aqueous solution at elevated temperatures," J. Phys.
Chem. 71:3121-3136. 1967.
27. Shapley, M. and L. Cutler. Rand's Chemical Composition
Program: A Manual. Report R-495-PR, Rand Corp.,
Santa Monica, California. June 1970.
56
-------
REFERENCES (Continued)
28. Shapiro, N. Analysis by Migration in the Presence of
Chemical Reaction. Report, P. 2596, Rand Corp.,
Santa Monica, California. June 1962.
29. DiToro, D.M. "Combining Chemical Equilibrium and Water
Quality Models." Manhattan College, Bronx, N.Y.
Oct. 1974.
57
-------
SECTION
LAKE 1 MODEL
BASIC STRUCTURE AND EQUATIONS
Because of the complexity of the overall modeling structure, a
simplified version of the lake has been developed to explore
the kinetic behavior in greater detail. This is indicated in the
modeling strategy of Fig. 4. The simplified model, designated
Lake 1, assumes that Lake Ontario is well mixed horizontally and
gradients are allowed to develop only in the vertical dimension.
Such a simplification obviously does not permit near-shore vs.
open-lake comparisons or the effect of the Rochester or Toronto
discharges on the local lake environment. However, an inspection
of the data on nutrients and chlorophyll does not appear to in-
dicate substantial horizontal gradients although variations do
exist in certain areas. (Lake Ontario can be contrasted in this
regard to Lake Huron with marked horizontal gradients from
Saginaw Bay or Lake Erie with important horizontal gradients
from the Western Basin to the Eastern Basin.)
Because of the complexity of the interactive systems and most
importantly, the computational time involved in obtaining solutions,
the horizontally well-mixed assumption appears to offer a reasonable
starting point for understanding the dynamic behavior of Lake
Ontario. Fig. 9 is a schematic of the Lake 1 model and shows the
three vertical segments that comprise the basis geometry. The
principal physical features included are (a) horizontal transport,
(b) vertical dispersion between the epilimnion and hypolimnion and
58
-------
NUTRIENT INPUTS
NIAGARA RIVER*)
TRIBUTARIES I
MUNICIPAL
INDUSTRIAL WASTES I
i
ENVIRONMENTAL INPUTS
"SOLAR RADIATION
WATER TEMPERATURE
LIGHT EXTINCTION
SYSTEM PARAMETERS
EPILIMNION
TRANSPORT
J L
f ^SETTLING | V f '
HYPOLIMNION
BENTHOS
Figure 9. Major physical features included in Lake 1 model
59
-------
TABLE 8
BASIC PHYSICAL DATA OF THE LAKE 1 MODEL
Surface
Segment Volume Depth Arfa ?low
Segment No. Interface (nTxlCT) % (m) (nr) cfs tiT/sec %
1
2
3 (sediment)
1-2
2-3
297,000
1,373,000
-
19
81
17
73.3
0.15*
1.64xl010
0.89xl010
43,500
188,500
1232
5323
19
81
Note: Vertical dispersion coefficient between segments No. 1 and 2
varied from 0+90 cm2/sec (0+778 m^/day)
* Segment #3 depth is arbitrary
-------
Cc) vertical settling of the phytoplankton. The Lake shown in
Fig. 9 has a defined epiliitvnion depth of 17 meters and a
hypolimnetic depth of 73.3 meters. The Lake 1 model is well-mixed
in the winter and spring, stratifies during summer and then
goes through a fall over turn. The effects are simulated by
control of the vertical mixing. Table 8 presents the physical data
of the Lake 1 model.
With the vertical segmentation as given in Fig. 9 and Table 8 and
using the general Eq. (3), the equations for Lake 1 are:
Segment No. 1 - Epilimnion:
V
Bl
EA
vlslk
Segment No. 2 - Hypolimnion:
ds_.
V2 ar* =
V2 S2k
Segment No. 3 - Sediment:
3k
(54)
where s,
and s,
are the concentrations of the kth system
, _, ,
K fiH. K. ,
inputted into the epilimnion and hypolimnion respectively,
an
are the flows into segments 1 and 2, V. is the volume
61
-------
of segment i, E,_ is the vertical turbulent exchange, A,2 is
the cross-sectional area between the epilimnion and hypolimnion
and Az is the average depth of the two segments. As shown in
Fig. 5 the Lake 1 model consists of 10 systems (k=l ...10).
The kinetics are as given previously in Section V. A review
of the Lake 1 computer program is given in Appendix A together
with a listing of all system parameters, inputs and sample
outputs.
MODEL CALIBRATION AND VERIFICATION
Several ways have been used to test the validity of the basic
model structure as given in Lake 1. First, a series of runs were
made to test the general behavior and sensitivitv of the model.
This is followed by inspection of the model output to determine
the general order of magnitude of the response of each of the
variables. The comparisons are made between the observed data
and the computed output for data collected during years prior to
IFYGL. These comparisons are qualitative in the sense that an
assessment of the degree to which the computed cutout agrees
with observed data is left to the analyst. It should also be
noted that other checks can be used in the model verification;
specifically variables that can be constructed from the model
output and compared subsequently to field measurements. The
total verification then includes both the primary variables
as given in the systems diagram (Fig.5 ) and secondary variables
constructed for the behavior of the model. The variables examined
therefore for calibration and verification are:
1) Phytoplankton chlorophyll "a" - iug/£
2) Primary productivity - mg carbon/m2 - day
3) Total zooplankton carbon mgC/£ ~
4) Bottom deposition rate - mg carbon/m - year
62
-------
5) Total phosphorus - mg/JJ.
6) Available phosphorus (qrthophosphate) mg/£
7) Total Kjelhdahl nitrogen - mg/£.
8] Ammonia nitrogen - mg/£
9) Nitrate nitrogen - rag/1
The approach was to examine past data (principally from the
years 1967-1971), run the Lake 1 model under different
parameter settings consistent with reported values and
compare to the observed data. The results of this sensitivity
and calibration analysis are given below. It should be noted
however, that even with variable parameter settings (such as
zooplankton grazing rate) the results are generally always
consistent with the observed data. The results of the comparison
to observed data are presented after the results of the pre-
liminary sensitivity and calibration runs.
Preliminary Runs and Sensitivity Analysis
A finite difference scheme is used to solve the equations using
explicit time-space differencing. For the Lake 1 model, a time
step of 0.5 days is used. For a one-year simulation, the central
processing unit (CPU) time required for execution is about 7
seconds on a CDC 6600 computer. Total CPU time required however
is 30 seconds with additional overhead converted to equivalent
CPU time being 115 seconds. The CPU time excluding overhead is
equal to about 1.4 milliseconds per compartment step. A
number of preliminary runs were made using the Lake 1 model
structure to l)test program elements 2) study the behavior of the
system and its sensitivity to various system parameters and in-
puts and 3) prepare a preliminary verification of the Lake 1 model
using data collected prior to IFYGL. These early runs therefore
63
-------
represent a type of "tuning" of the model using past data as a
basis for additional verification of the IFYGL data. The Lake 1
model has been used to examine several areas including, for ex-
ample variable levels of spring and fall vertical mixing, and
settling velocities for phytoplankton in addition to zooplankton
grazing rate, and effect of half-saturation constant for phos-
phorus, Some of the output from these areas is discussed below.
Variable Vertical Mixing
Vertical mixing plays an important role in the dynamics of phyto-
plankton population in lakes. Several paths were followed in order
to obtain an estimate of the vertical dispersion in Lake Ontario.
A survey of the literature on vertical mixing indicated that
generally the vertical dispersion coefficient is several orders
of magnitude lower than the horizontal coefficient. Oceanic
vertical dispersion coefficients are in the range of 10-100 cm2/sec
and generally decrease with depth. Recently, during IFYGL, Murthy
1 2 ?
et al estimated values from O.lcm /sec to 22cm /sec depending on
effects such as wind conditions and density stratification. In
the Lake I work, the vertical dispersion was therefore treated
as a parameter and sensitivity analyses made under an expected
range from complete stratification (no mixing) throughout the
•\
year to variable mixing up to 90cm /sec. More extensive an-
alysis of vertical mixing was carried out in the Lake 2 and 3
models using temperature as a tracer variable. The details of
these analyses are given in Section VI.
Figure 10 shows the effect of vertical dispersion on chlorophyll
a. The runs include a sinking velocity of 0.05 m/day. A
bimodal distribution of phytoplankton occurs in all cases due to
the interaction of the two higher zooplankton levels used in
64
-------
o
MJ J ASOND
^ 3.5
E
D U^ 0 UN C
• • • .
0> CSJ CSJ r- 1 f
.01 'NOISifldSIQ
__, l.U
e
>
(2)
» r
% ,
\ /
\ i
(i)
^^^^^M _^^^^^^^*
\\ /
1 1 1 1 M 1 1 1 /I 1
JFMAMJJASO
—
—
1
N D
- 100
- 50
Figure 10. Effect of vertical mixing
65
-------
these runs. As populations build up in the spring, the her-
bivorous zooplankton increase. This together with nutrient
depletion decreases the phytoplankton population and at about
the same time, the carnivorous zooplankton prey on the next
lowest trophic level. The phytoplankton can then increase
again in the late summer.
As shown in Fig. 10, the main effect of vertical mixing is in
the spring where populations increase more rapidly when no
mixing is allowed and reach a peak earlier. With mixing, the
first peak is delayed and reduced in magnitude. The results
in the summer and early fall are generally comparable under the
different dispersion regimes although the timing of maximum and
minimum is changed.
Variable Phytoplankton Settling Velocity
A number of analyses using the Lake 1 model have been made to
examine the behavior of the phytoplankton population under
different settling velocities. The reported values for phyto-
plankton sinking in quiescent waters (see Section III) are
apparently too high to support adequate growth in the Lake 1
model. This is due in part, to the crude vertical definition
and the interaction between the sinking of phytoplankton and
their physiological state. The settling velocity has been
treated as a parameter ranging from 0 to 0.5 meters/day.
Figure 11 summarizes the results from several analyses using
different sinking velocities. For these runs, zooplankton graz-
ing was set at 0.06 1/mg carbon-day-°C, no vertical mixing was
used and the Michaelis constants were 25yg/& and 10 ug/S- for
nitrogen and phosphorus respectively. At a velocity of 0.5 m/
day phytoplankton populations never exceed about 3yg chlor/1
66
-------
^ 25
en
10
cc
31 O 5
n —I J
o
0
0.20
0.15
0.10
o
z
o
0.05
0
0.15
1 o-10
8 0.05
0
PHYTOPLANKTON
SETTLING VELOCITY
0.5m/day
0.05 m/day
(1)
F ' M1 A"M'J'J'A'S'O'N'D
HERBIVOROUS
ZOOPLANKTON
CARNIVOROUS
ZOOPLANKTON
FOR
ABOVE
ERBIVOROUS
ZOOPLANKTON
ARNIVOROUS
/ ZOOPLANKTON
FOR
~ (2)
ABOVE
7" F ' M ' A'MM'J'A'S'O'N
Figure 11. Effect of phytoplankton settling velocity,
67
-------
and total zooplankton carbon (Figure 8) never exceeds about
0.08 mg Carbon/1. Both values are less than observed. The
reason for the low levels is that under a settling velocity
of 0.5 m/day (which for the 17 meter depth of the epilimnion
represents a "decay" coefficient of .03/day), the phytoplankton
are not retained in the upper layer long enough to undergo net
growth. As a consequence, zooplankton levels are also low and
the nutrient concentrations remain high and are not reflective
of observed nutrient depletion. Using a velocity of 0.05m/day,
the behavior of the phytoplankton biomass is quite different
as shown in Figure 11. Now, the lower trophic level has a chance
to grow and a reasonable predator-prey relationship begins to
develop.
The half-saturation constant for phosphorus, K , has an important
effect however. An order of magnitude reduction in this co-
efficient to lyg/2, does permit growth of the phytoplankton at
the higher settling velocity of 0.5m/day. Figures 11 and 12 show
this effect and indicate the trade-off between a population adapted
to low phosphorus concentration and settling velocity. These
results also tend to indicate the possible interaction of settling
velocity and the active growth phase of the population.
The sinking velocity of phytoplankton also has an important effect
on the nutrient uptake as shown in Fig. 13, which is a plot of the
nitrate and available phosphorus concentrations calculated under
the conditions of K = lO^ig/Jl. In addition, Fig. 13 shows the
mp
nitrogen and phosphorus limitation terms, i.e., the ratio of
total inorganic nitrogen to total inorganic nitrogen plus the
Michaelis for nitrogen and similarly for available phosphorus
(see Eq.8). It can be seen that nutrient uptake and growth
68
-------
o
20
15
10
5
0
SETTL ING VELOCITY 0.05 ml day
SETTLING VELOCITY 0.5 m /day
M
CO
e
z_ 0.20
f\ CT>
8 E. °-15
o.io
0.05
0
S2~ 0.04
0.03
^£0.02
H
0.01
0
o
OS
O
o.
uo_
J 30 F 60M90 A120M150J 180J210 A240S 270°300N 330 D360
;_ 0.20 -
'o
I4 0.15 -
M J J A S 0 N D
CD
5
<
jrFTMIA'M1J1J
ON D
Figure 120 Effect of phytoplankton settling velocity,
mp "~
69
-------
§ 0.3
£-
z z 02
gf
£ 0.1
(1) SETTLING VELOCITY 0.05 m/day
{2) SETTLING VELOCITY 0.50m/day
J ' F ' M ' A MJ J ASON D
o
»—
<
§ % 1.0
-• E
oz 0.5
o
oc
J ' F ' M' A ' M' J ' J ' A ' S ' 0 ' N ' D
i 1—:—i—:—i 1
JASOND
Q_
+
a
o.
to O
. 1.0
0.5
NO LI Ml TAT ION
JFMAMJJASOND
Figure 13. Effect of phytoplankton settling velocity on nutrients
and nutrient limitation, Kmp= 10yg/£
70
-------
limitation is minimal for the case of sinking velocity = 0.5m/day.
This is a result of the minimal phytoplankton growth. At the
lower sinking rate, however, nitrate uptake is increased but as
indicated in the upper curves of Fig. 13 nitrogen does not sig-
nificantly limit growth for the case of Kmp= lOyg/S,.
The lower curves representing the phosphorus dynamics however
do exhibit a limiting effect. If attention is directed to the
phosphorus limitation term, it is seen that for both settling
velocity cases, levels of phosphorus are such that a limitation
of .50 - 0.60 prevails during the early and later parts of the
year. However, substantial differences occur in the spring and
late summer. At day 135, a minimum value of about 0.2 is cal-
culated indicating that phosphorus is acting as a significant
limiting factor in the phytoplankton growth. This helps explain
the decrease in phytoplankton biomass beginning at day 120,
(see Fig.11). Two effects are occurring: a) herbivorous zoo-
plankton are growing rapidly and b) phosphorus levels are being
depleted to below the half-saturation constant thereby acting
to reduce the growth rate. It is interesting to note then that
a biomodal distribution in phytoplankton can be obtained without
a species differentiation. The latter is often offered as the
explanation for the observed two peaks in the phytoplankton.
In order to accomplish this, however, at least two trophic
levels must be included above the phytoplankton. This permits
a higher order predation, e.g. carnivorous zooplankton which
reduces the lower zooplankton level. The reduction permits
phytoplankton to grow again in late summer. This is explored in
greater detail below.
71
-------
Effect of Michaelis Constant
As indicated in Eq. (5), the growth of the phytoplankton is con-
sidered to be limited by a Michaelis Menton product term of total
inorganic nitrogen and available phosphorus given by Eq. ( 8 )•
The dynamics of growth therefore are affected by the levels of
K and K , the half-saturation constants for nitrogen and
mn mp
phosphorus, respectively. Several sensitivity runs have been made
therefore specifically examining two levels of K ,1 and 10 yg/&
both of which encompassed the range of published values. Phospho-
rus is emphasized because of the waste removal programs presently
underway in the Great Lakes Basin.
Figs- 14 an<^ 15 show the results of two computations at an order
of magnitude difference in K . For both runs, settling velocity
was ,05m/day, K (nitrogen half-saturation constant) is 25ug/&
and zooplankton are included as before. As seen in Fig.14 at the
lyg/£ level, phytoplankton biomass grows earlier more rapidly and
reaches a peak value of 14yg/£ chlorophyll/1 at the end of April.
At the 10yg/fc value, growth is later, slower and reaches a maxi-
mum value of about 9yg/Jl almost a month later at the end of May.
Values in the late summer are comparable. The nutrient plots shown
in Fig. 14 indicate that in May, phosphorus is limiting under the
lyg/£ case while in August-September, nitrogen appears to be
limiting. This is shown more clearly in Fig. 15 - which is a plot
of the nitrogen and phosphorus limitation terms under both conditions
It can be noted that at lyg/£ the growth is limited by phosphorus
for only a short period in May and the system tends to be nitrogen
limited throughout the rest of the growing season. On the other
72
-------
1"'1
I - 1 - 1 - 1 - 1 - 1
JFA/TAMJJASOND
E- 0.10-
JFNIAMJJASOND
i 1 1 1 1 1 I 1
JFMAMJJASOND
i 1 1 r 1 1 1 1 1
JFMAMJJASOND
Phosphorus Michaelis Constant
10ug/l
Figure 14. Effect of Michaelis constant for phosphorus
on primary variables
73
-------
2 L°
O
£ z 0-8
1 I 0.6
O
o
0.4
0.2
0
NO LIMITATION
J1 J 'A'S' O'N
to
Ou
o
Q_
1.0
0.8
0.6
0.4
0.2
0
NO LIMITATION
J ' F ' M ' A ' M ' J ' J ' A l S ' 0
TIME OF YEAR
PHOSPHORUS MICHAELIS CONSTANT 1 vq/\
— ioug/1
Figure 15. Effect of Michaelis constant for phosphorus
on nutrient limitation of nitrogen and phosphorus
74
-------
hand, at a phosphorus Michaelis constant of lOpg/Jt the system
is more phosphorus limited throughout the year although nitrogen
is still an important limiting interest in late summer.
The results indicate the importance of reasonably accurate estimates
of the half-saturation constant. This will be particularly important
during any simulation of nutrient reduction programs. The results
also indicate that reliance on a single nutrient as "the limiting
nutrient" does not recognize the interaction between nutrients and
may be quite misleading in estimating the effects of a control pro-
gram.
Effect of Alternate Growth Formulation
As discussed in Section III, Eqs. (5) and (9) represent different
expressions for the growth term in the phytoplankton equation. In
order to examine the behavior of the response to each form, the
Lake 1 growth kinetics were modified to incorporate Eg. (9). The
results for three output variables are shown in Fig. 16. Both runs
displayed therein are for identical conditions except for the
alternate growth formulation as indicated in the figure.
Eq. (9) results in immediate and rapid growth of the phytoplankton
with subsequent total depletion of the phosphorus. The dynamic be-
havior of the results using the formulation of Eq. (9) do not even
approximately agree with observed data. This is especially true
for the high levels of phytoplankton in the early and late months
of the year and the zero concentration of phosphorus throughout
the better part of the year. The dynamic shape is also not con-
sistent with the observed data (see Fig. 3 ). From this one
comparison, then, Eq. (5) the original product formulation is still
to be preferred.
75
-------
O QC
ii
0.0
15-
10-
5-
0
J F M A M J
J A S 0 N D
LU
O
0.30
O.ZO-
0.10-
0
JFMAMJJASOND
co S —
J FMAMJ JASOND
N'-P' ---G =
Figure 16. Responses due to two photoplankton growth
formulations
76
-------
Time to Dynamic Equilibrium
With the relatively long detention time of about 8 years for
Lake Ontario, the time to reach a dynamic steady state (i.e.
where the solution becomes periodic) is an important management
consideration. For a single non-conservative variable in a
completely mixed lake, the solution for a mass input of W is
s = JL_ [ i-e-(k+1/t0)t + e-(k+i/to)t (55)
Q+KV o
where s is the initial condition of s and t is the detention
o o
time given by V/Q. For K=0, the conservative case, the time
to reach 95% of a new steady value is about 3 detention times,
which for the whole of Lake Ontario is about 24 years. Further-
more, the initial conditions will also take approximately that
long to "die down" and not influence the solution. Any solution
therefore that is calculated for only one year (as per the previous
figures) is dominated by the initial conditions and the effect of
the mass input is negligible. For verification purposes, however,
a single year run is sufficient to explore the general behavior
of the system. For any projections of the effects of reduced in-
puts, longer runs are needed. Therefore, the solutions must be
run for a period of time until the initial conditions for all
variables are repeated year by year. With the complex interactions
of the Lake 1 model, this time to reach a dynamic equilibrium is
not immediately obvious. Accordingly, several long-term runs were
prepared to illustrate this effect.
Fig. 17 shows a time history plot of 16 years of solution equivalent
to two detention times. The conditions are comparable to previous
runs but use a Michaelis phosphorus constant of lOyg/Z. As shown,
77
-------
00
a.
o
CXL
O
O
z
o
Q_
O
10
5
0
10
5
0
3
10
11
12
13
14
15
TIME, years
8
16
Figure 17, Sixteen year computed output of phytoplankton chlorophyll
-------
peak phytoplanJcton levels in the spring gradually decrease from
greater than 9yg/£ in the first year to about 7yg/£ by the
eighth year after which time, a dynamic equilibrium has been
reached. Thus, a time of about one detention time for Lake
Ontario was required for the phytoplankton to reach a steady
state. However, due to internal cycling between the various
nutrient sub-systems shorter or longer times may be required
for some of the nutrient forms. This is shown in Fig. IQ
which indicates system response under two phosphorus Michaelis
constants. The upper figure indicates a parallel decrease in
the phytoplankton maximum under the two conditions and although
the two runs began with a peak difference of almost 3yg/l, the
difference is about half of that at equilibrium.
The lower two plots of Fig. 18 are quite interesting and illus-
trate further the difference in the dynamic behavior under the
two phosphorus Michaelis levels. Under the lower level of lyg/j,
equilibrium is reached almost immediately in both phosphorus and
total inorganic nitrogen. Under an order of magnitude increase
in the level to 10yg/& dynamic steady state is not reached for
the inorganic nitrogen until almost two detention times or 16
years. With the lower level of lygp/£ phytoplankton growth is
increased due to the increased ability to utilize low phosphorus
concentrations. As a consequence, however, nitrogen utilization
is increased and as shown in the middle plot, the system becomes
nitrogen limited. The times of the year of the nitrogen and
phosphorus limitation are different however as shown in Fig. 15
In essence, then these runs illustrate the need to carry out any
simulations of a proposed nutrient reduction program out to about
8-16 years to determine the expected phytoplankton levels under a
79
-------
0. O
ii 4
xo
i 2
- \ PHOSPHORUS MICHAELIS CONSTANT
\
V 1M9/I
1 3 5 7 9 11 13 15 17
z 0.10
5^
I 0.02
2
^
10 M9/I
NITROGEN MICHAELIS CONSTANT .025 mg/l
i i
1 3 5 7 9 11 13 15 17
3^ 0.010
|?O.OOS
IS 0.006
ii o-004
Eg 0.002
0.000
i—i—ffi~1~i'"l
1 3 5 7 9 11 13 15 17
TIME, years
Figure 18. Long term behavior of model components under
two phosphorus Michaelis constants
80
-------
different waste load input. If however, the loads are changing
from year to year, as, for example, in a linear increasing or
fashion, then even longer times may be required. Indeed, under
changing loads, a dynamic equilibrium may never be theoretically
possible. For many practical purposes, the system will probably
reach steady state in a 10-20 year time horizon after significant
changes have been made.
Preliminary Comparison
Figs.19 and 20 show one of the early comparisons of observed data
and model output. The Michaelis constant for phosphorus was
10vtg/£ and settling velocity was 0.05m/day. In the preliminary
comparisons, a variety of similar plots was generated. Figs. 19
and 20 therefore represent a type of a first level of calibration.
In general, the overall comparison is quite good and indicates
that even with first, approximate analysis, results are obtained
from the model that are in general agreement with the observed
data. This tends to support the contention that the basic model
structure does represent the principal features of the phytoplankton
biomass.
Several details of Figs.19 and 20 however remained to be explored in
greater depth before a "final" verification could be obtained. For
example, in the preliminary runs, the fall peak in phytoplankton is
only approximately reproduced and in addition the total zooplankton
carbon is generally too high. Finally, the dynamic behavior of the
nutrients such as phosphorus and ammonia, while generally satis-
factory, could be improved. Using therefore, the sensitivity runs
discussed previously and the results of the preliminary comparisons,
further verification analyses were conducted to construct a coherent
81
-------
^ 20
en
J 15
& 10
£3
<* 5
o 0
CCIW 1967 Data Mean ± 1 Standard Deviation
_ Range of observed
data
30 60 90 120 150 180 210 240 270 300 330 360
z- 0.20
o
QQ
5 0.15
o
2^0.10
3 0.05
Q_
O
o n
M u
Zooplankton No. 1
-r i
M
MJJ
N
O
CO
a:
0.20
0.15
o _
^ "g'O.lO
Zooplankton No. 2
D-
o
o
ISI
5 0.05
0
0.8
&
5 0.4
o
t: 0.2
30 60 90 120 150 180 210 240 270 300 330 360
CCIW1967 Data Mean ± 1 Standard Deviation
I I
FMAMJJ
SOND
Figure 19. Preliminary comparison of model output and observed data.
82
-------
. 0.20
u
1 0.15
0.05
0.00
CCIW 1967 Data Mean ± 1 Standard Deviation
30 60 90 120 150 180 210 240 270 300 330 360
o
o
P
0.40
0.30
<
a:
~5i0.20
E
0.10
to
=j
or
O
a.
to
a.
a
CO
0.00
0.40
0.03
|>0.02
0.01
0.00
CCIW Data Mean ± 1 Standard Deviation
J FMAMJJASOND
CCIW Data Mean ± Standard Deviation
.1 .1 I
30 60 90 120 150 180 210 240 270 300 330 360
o _
0.020
0.015
g
g 0.010
s
g 0.005
^ o.ooo
i i i
J L
J F M A M J J
S 0 N D
Figure 20. Preliminary comparison of model output and observed
data, continued
83
-------
picture of. the dynamics of the lake. The goal then in this
final verification phase of Lake 1 model was to produce a
set of parameters and accompanying hypotheses that represented
as good a comparison as possible to available data.
RESULTS OF "FINAL" VERIFICATION ANALYSIS
Data were available for the years 1967-70 collected by the
Canadian Centre for Inland Waters for this final verification
phase. These data were supplemented by other observations in
the literature. Values for the various coefficients, parameters
and exogenous variables were obtained from published sources as
discussed previously and all values used are within reported
ranges in the literature. More than 80 runs have been made of
Lake 1 to examine the effects of various phenomena such as
settling velocity, vertical mixing, zooplankton predation and
half-saturation constants for nitrogen and phosphorus. Out of
these runs has emerged a consistent set of parameter values which
satisfactorily explains the observed data on a variety of
different variables. A plausible explanation for the dynamic be-
havior of the phytoplankton is therefore possible and is disucssed
below.
Figs.2land 22 show the verification of the Lake 1 model upper layer
(0-17 meters) while Fig.23 shows the verification for the hypo-
limnion (data in range, 50-150 meters). Table 9 shows the
principal coefficients used for the runs. As shown, the verifi-
cation for plankton in the epilimnion is quite good and satis-
factorily duplicates the spring peak/subsequent mid-summer die off
and fall bloom. Five other variables in addition to the phyto-
plankton chlorophyll are satisfactorily verified in the com-
parison. It should be noted here that it is generally not
84
-------
- 20
01
3d
% =
1= Q.
>- O
31 OL
O- O
15
10
1970 RANGE OF MEAN ±1 STANDARD DEVIATION
1%Q \ ////,. LAKE-WIDE MEAN.
o
0.16 1-
0.12
^L
01
E 0.08
Q_
8
rsi
0.04
0
NOTE: ZOOPLANKTON CARBON
IS FOR 0-50 METER HAUL
IFI\AAMI IASOND
30 60 90 120 150 180 210 240 270 300 330 360
o
o
«»>
0.6
0.4
0.2
0
j ' F'M'A
M'J'J'A'S'O'N
Figure 21. Lake 1 model verification, 0-17 meters
85
-------
LAKE-WIDE MEAN. 0-17 METERS
• 1969
D 1968
o
s~ s^jWss//^*^ smr"^
M J J
RANGE OF MEAN ±1 STANDARD DEVIATION
J 30 F 60 M 90 A 120M 150J 180 J 210A24QS 270 °300*330°360
JFMAMJ J A S 0 N D
Figure 22. Lake I model verification (cont.), 0-17 meters
86
-------
o
o;
o
O
Q£
0.6
0.2
0
0.15
0.10
0.05
LAKE-WIDE MEAN, 50-150 meters
• 1970
J ' F ' M
A ' M ' J ' J ' A ' S ' 0 T N ' D
RANGE OF MEAN ± 1 STANDARD DEVIATION
I'F'M'A'M'J'J'A'SOND
30 60 90 120 150 180 210 240 270 300 330 360
I I I 1 I
J A S 0 N D
Figure 23. Lake 1 model verification, 50 - 150 meters
87
-------
TABLE 9. PRINCIPAL PARAMETER VALUES USED IN LAKE 1 MODEL OUTPUT
SHOWN IN FIGURES 21-23
00
00
Notation
mn
Kmp
K
g
12
Lzp
r
"z
K]
K
K
22
op
w
Name
Half-saturation constant-nitrogen
HaIf-saturation constant-phosphorus
Grazing rate for zooplankton
Half-saturation constant-phytoplankton
Endogenous respiration rate-phytoplankton
Zooplankton conversion efficiency
Zooplankton endogenous respiration rate
Decomposition rate of organic nitrogen
Ammonia to nitrate nitrification rate
Decomposition rate of organic phosphorus
Nitrogen-chlorophyll ratio
Phosphorus-chlorophyll ratio
Carbon-chlorophyll ratio
Settling velocity of phytoplankton
Value
Unit
25
2
.06
10
0.1
0.6
.001
.00175
.002
.007
10
1
50
0.1
ygN/l
ygp/l
l/mg C-Day-C°
yg Chlor/1
days'1 (@20°C)
(days - °C)~1
(day - °C)~1
(days - ^C)'1
(day - "C)"1
ygN/yg chlor
ygP/yg chlor
ygC/yg chlor
m/day
-------
possible to obtain the verification shown in Figs. 21-23 by
arbitrary specification of the coefficients, so that the
results shown do not simply represent a "curve fitting" ex-
ercise. Rather, the results, which are based on the theory
discussed earlier, are representative of the observed data
in a much more meaningful way than just through arbitrary
statistical coefficients. Since the model permits compu-
tation of the components of the dynamic behavior, consider-
able insight can be gained by examining the influence of
various phenomena on the growth and death rates of the
phytoplankton.
Fig.24 shows the dynamics of the kinetic growth rate of the
phytoplankton (G in Eq.C4)). As shown, maximum growth rate
reaches a peak in mid-September in the Lake 1 model at about
1.8/day and then decreases rapidly due to the fall overturn.
The reduction in the maximum growth rate due to light is sub-
stantial and shifts peak growth to early August. A further
reduction in growth rate is due to the nutrient limitation so
that the resultant growth rate, G is lowest at the end of May
and peaks at about 0.25/day in August. The peak value repre-
sents an overall reduction of 90% from saturated growth con-
ditions. Primary productivity values calculated from the
growth rate and phytoplankton biomass Eq-Cll) give peak values
of about 600 mg carbon/m -day at the end of August and values
of about 500 mg C/m2 day at the height of the spring bloom.
The former value is in general agreement with Glooshenko
while the latter is somewhat lower than the Glooshenko's
2
estimated spring values of 500-1100 mg C/m -day.
89
-------
Dynamics of Phytoplankton Growth Rate
J 30 F60 M90 A120 M150J 180 J 210 A 240 S 270° 300 N 330 D 360
2.0
1.8
£ 1.6
1"
Sl.O
o
g 0.8
5 0.6
S 0.4
0.2 -
0
MAX I MUM GROWTH
RATEATEPILIWNION
WATER TEMPERATURE
REDUCTION DUE TO
LIGHT LIMITATION
- RESULTANT PHYTOPLANKTON
GROWTH RATE
REDUCTION DUE TO
NUTRIENT LIMITATION
2.5
I
2.0
1.5
o.
o
z
_1
OQ
O
O
1.0
CO
0.5
Figure 24. Dynamics of the phytoplankton growth rate
90
-------
Fig.25 shows the dynamics of the kinetic death rate of the
phytoplankton (Dp in Eq. (4) ). Three effects are to be
noted: a) phytoplankton settling, b) water temperature
effect on endogeneous respiration and c) effect of zooplankton
grazing. Peak resultant death rate is 0.25/day at the begin-
ning of August and is primarily due to zooplankton predation
and quantitatively explains the mid-summer decline in phyto-
plankton biomass.
The dynamics of phytoplankton net production, i.e. chloro-
phyll/liter-day, dP/dt, are shown in Fig. 26. Net production
here includes the kinetic interactions of growth in Fig.21,
death and predation in Fig.25 and the effect of vertical mix-
ing and lake outflow. The results shown in Fig.23 summarize the
basic hypotheses of phytoplankton dynamics in Lake Ontario.
The spring growth phase to approximately mid-May is due primarily
to increasing light and temperature. There are little zoo-
plankton as yet, so growth continues until nutrient limitation
becomes significant in early June. (This is principally
phosphorus limitation as discussed below). The spring growth
and spring peak are therefore simply described by a nutrient
interaction effect with light and temperature dominating the
growth and death terms. Following the spring peak however,
the situation becomes considerably more complex and indeed a
significant effort was devoted to unraveling the complex inter-
actions which lead to a mid-summer decline and subsequent
broad full peak.
It is hypothesized as computed in Fig.26, that during mid-summer,
and following the nutrient limitation effect, zooplankton grazing
91
-------
J 30 F60 M90 A120M150 J 180 J 210 A 24QS 270° 300 N 330° 360
fr
•a
S °'25
« 0.20
§ 0.15
| 0.10
z
S O-05
i °
CX-
RESULTANT PHYTOPIANKTON DEATH RATE
EFFECT OF WATER
TEMPERATURE
h EFFECT OF
PHYTOPIANKTON
I- SETTLING RATE
(0.1 meter/day)
ZOOPLANKTON
GRAZING EFFECT
J F
^^AXH^J< f, f. * » *• **...«s. »_g
A'M' J '
0 N D
Figure 25. Dynamics of the phytoplankton death rate
92
-------
10
1 5
I I
J 30 F60 M90 A120M150 J 180 J 210 A 24CS 270° 300 N 330 D 360
ACTIVE
GROWTH INCREASING LIGHT
AND TEMPERATURE
NUTRIENT REGENERATION
NUTRIENT LIMITATION
NUTRIENT LIMITATION
FALL OVERTURN
ZOOPLANKTON GRAZING
DECLINING
GROWTH
Figure 26. Dynamics of the phytoplankton net production
93
-------
becomes increasingly important and accounts for the minimum
values of biomass in July. However, at increased zooplankton
grazing, biological recycling of nutrients becomes more sig-
nificant. Nutrients are released back to the epilimnion
through excretion as well as phytoplankton endogenous
respiration. In late July then, growth exceeds death due to
nutrient regeneration and an active growth phase begins again.
Nutrients however, are already at low levels so growth is not
pronounced and proceeds slowly (as shown by the net production
approaching zero in September but not yet entering a negative
growth phase). The fall overturn of the lake then reduces
the growth of phytoplankton biomass primarily because of mixing
with the colder hypolimnion waters.
Figs. 27 and 28 explore these dynamic effects at another level
of detail. The model permits calculation of nutrient limitation
effects as well as the flux of nutrients due to zooplankton ex-
cretion and release of nutrients due to endogenous respiration
by the plankton. Fig. 27 shows the nutrient effects. As seen,
the spring bloom is halted essentially by a reduction in phosph-
orus that is quite pronounced and growth is reduced by a factor
of about 0.35 due to the low phosphorus levels. At the same
time, nitrogen is not limiting and is at a level twice that of
the phosphorus. Nitrogen however does become limiting during
July and again in September at which time phosphorus also exerts
=> limiting effect on growth. Fig.27 shows therefore that the
spring peak is primarily phosphorus controlled. The dynamics
after the spring peak are however controlled by both nitrogen
and phosphorus (as well as zooplankton predation).
94
-------
J FMAMJJASOND
. 1.0
z
o
t 0.8
t: z
1 +c 0.6
o
o
CK:
0.4
0.2
0
NO LIMITATION
APPROXIAMTEPERIODCF
NITROGEN LIMITATION
J'F'M'A'M'J'J'A'S'O'N'D
30 60 90 120 150 180 210 240 270 300 330 360
o
i—
•<
1.0
0.8
C* ^ r. .
o-0.4
a.
g 0.2
Q_
0
NO LIMITATION
1 '
J ' F
M ' A '
APPROXIMATE PERIOD OF
PHOSPHORUS LIMITATION
"TTT
J ' J ' A ' S ' 0 ' N ' D
Figure 27. Dynamics of nitrogen and phosphorus limitations
95
-------
Fig.28 shows the estimated sources of phosphorus recycled to
the epilimnion by both biological effects and vertical mixing
effects from the hypolimnion. As shown, during early spring
the hypolimnion contributes almost twice as much nutrients to
the epilimnion as from external inputs and by early summer in
June-July, biological recycling reaches almost five times the
mass rate of external phosphorus sources due to waste residuals,
Niagara River and other contributors. Zooplankton excretion
recycles almost twice the input rate and in late summer-early
fall, biological recycling reaches another peak but then drops
rapidly due to fall overturn. The flux of nutrients from the
more nutrient rich hypolimnion during October is not sufficient
to offset the drop in water temperature in the epilimnion.
Biological recycling therefore drops rapidly in October and
the phytoplankton biomass declines (see Fig.21). The results
shown in Fig.28 indicate quantitatively the significance of
biological recycling and vertical mixing relative to the level
of external nutrient inputs. A reduction in input phosphorus
load therefore may not necessarily result in a concomitant re-
duction in biomass. Further, the time to reach a new equi-
librium biomass level will extend beyond a single year due, in
part, to the recycling effects discussed previously, and long
detention times in Lake Ontario.
In summary, the Lake 1 modeling effort, to date, has indicated
that a considerable degree of quantitative insight into the be-
havior of phytoplankton biomass can be obtained with a spatially
simplified model that consists of an epilimnion, hypolimnion and
sediment but is horizontally well mixed. The model, using
acceptable ranges for plankton growth parameters, verifies observed
data on phytoplankton, zooplankton and nutrient levels and as
such can form a first basis for estimating the lake-wide effects
96
-------
Biological Recycling
o 8.
UJ o
a: O
o;
O
O
a.
1.6 -
1.2
1.0
0.8
0.6
0.4
0.2
0
TOTAL EXTERNAL INPUTS
(NIAGARA RIVER,
MUNICIPAL, AND
TRIBUTARIES
5.0
4.0
3.0
2.0
1.0
ZJ
o.
ce
O
CL.
-
O
O
I
Q-
«>-l
O
Q_
Figure 28. Biological and hypolimnetic recycling
of phosphorus
97
-------
of a nutrient reduction program.
The results indicate that the spring growth phase and peak
phytoplankton biomass are primarily.controlled by increasing
light and temperature and phosphorus limitation. The mid-
summer minimum in phytoplankton is estimated to be due
principally to zooplankton grazing and nitrogen limitation.
The broad fall peak in phytoplankton is a complex interaction
of nutrient regeneration (up to five times the external nutrient
inputs)r subsequent nutrient limitation and then fall overturn.
Both nitrogen and phosphorus are important nutrients in this
dynamic succession.
From a water quality viewpoint, peak values of biomass in the
spring are important, and are generally limited by phosphorus
concentrations which reduce the nutrient growth rate by 35%.
Nitrogen also has an additional effect of about a further 10%
drop in growth rate. In the later summer months, it appears
that the system is nitrogen limited as well as phosphorus
limited. There may be important implications of this hypothesis;
to wit, nutrient removal programs aimed solely at phosphorus
may not achieve a sought-after objective especially in the late
summer recreation months. This is especially relevant when one
recognizes that generally this period of the year is dominated
by the green and blue-green algae, which have greater potential
for water use interference than do the diatoms.
A period of 10 - 20 years may be necessary for Lake Ontario to
reach a new level of dynamic equilibrium in phytoplankton biomass
after a nutrient reduction program is instituted. The long de-
tention times, recycling effects and large store of nutrients in
the hypolimnion are the principal reasons for this long response
98
-------
Therefore, the Lake Ontario presently being observed represents
conditions of about 10-20 years ago and since loads into the
upper basin Ce.g. Lake Erie) have been increasing steadily over
the years, Lake Ontario is obviously in a "slowly" varying state
from year to year. Thus, the effects of a waste removal program
may not be apparent for approximately one-two decades and if
some residual loads are allowed to continue to grow, significant
changes in Lake Ontario quality may "never" be apparent. This
also implies that a detailed program of sampling of effluents
and tributaries and concurrent lake wide sampling are correlated
only through an approximate ten year lag. The results also in-
dicate that surveillance of large lake systems such as Lake Ontario
may best be accomplished in a time scale of perhaps once every
five years rather than year to year.
All of these observations on response times apply to the lake as
a whole. The more detailed Lake 3 model which includes horizontal
detail will obviously reflect more localized situations which may
(or may not) reach dynamic equilibrium in shorter times.
99
-------
SECTION VI
REFERENCES
1. Murthy, C.R. , G. Kullenberg, H. Westerberg, and K.C.
Miners. "Large scale diffusion studies," in IFYGL
Bull. No. 10, pp.22-49. NOAA, Rockville, Md. May 1974.
2. Glooshenko, W.A., J,E. Moore, M. Munawar and R.A.
Vollenweider. "Primary production in Lakes Ontario and
Erie: A comparative study," Journal of Fisheries Research
Board of Canada. 31(3):253-263. 1974.
100
-------
SECTION VII
LAKE 2 ANALYTICAL INVESTIGATIONS
INTRODUCTION
As illustrated in Fig. 4 a strategy of essentially parallel
investigations has been pursued in the development of Lake Ontario
phytoplankton models. The Lake 2 investigation has resulted in a
theoretical analysis, as opposed to numerical investigations, of
the vertical interactions in phytoplankton populations. Two
specific areas are addressed: the behavior of the asymptotic
population growth rate and the effect of vertical settling velocity.
The form of the equations for a quantitative theory of the dis-
tribution of phytoplankton biomass have been known for some time
and their applications to various Great Lakes problem settings is
ongoing^2f3'4^. An analyses of the conditions required for temporal
steady state accompanied their introduction; however no time
varying analysis of their general properties has been made.
The condition under which a population can maintain itself against
transport loss are related to the analysis to be presented insofar
as an eigenvalue problem results, which gives the bounds on the
parameter groups that result either in population increase or de-
crease. What has not been recognized is that if the population does
increase it asymptotically approaches a condition which is a sub-
stantial simplification of the general solution of the biomass con-
servation equation. In particular, it can be shown that there exists
an asymptotic population growth rate, y, which is independent of
position, and an asymptotic population distribution, which completely
101
-------
characterize the solution to the biomass equation. The asymptotic
solution depends on the kinetic and transport parameters and to
find these dependencies it is necessary to make some simplifying
assumptions. However, the general existence of the asymptotic
solution follows directly from the assumption of temporally con-
stant parameters. It is shown subsequently that in fact this pre-
diction is born out by observed phytoplankton biomass distributions.
The specific spatial setting for the analysis presented sub-
sequently is in the vertical dimension. The motivation for the in-
vestigation is a desire to calculate the effects of vertical settling
velocity on phytoplankton population development; in particular,
what are the important dimensionless parameter groups. Depending on
the question to be addressed a series of dimensionless plots are
presented which give quantitative answers to the importance of the
various parameters. In this way an economy is realized so that the
complexity of even these idealized situations is kept to a manage-
able level.
A final product of this analysis is a demonstration that mathematical
theory and biological fact can be brought into accord if the com-
plexities of the latter are aggregated in such a way that the larger
scale relationships are susceptible to a necessarily simplified
mathematical treatment. The use of biomass as the primary dependent
variable is a biological example; the use of a dispersion coefficient
is a physical example. The result is a tractable mathematical problem
which yields an array of information that can then be checked against
observation. In this case, the prediction is the existence of an
asymptotic population growth rate which is independent of position.
102
-------
THEORETICAL CONSIDERATIONS
The conservation of phytoplankton biomass in the presence of
vertical transport processes requires that the rate of change
of the biomass with respect to time is the result of the balance
between the kinetic source and sink of biomass and the vertical
transport, the effect of which is to remove biomass from the
euphotic zone. The equation is quite well known:
(E H
where E = vertical dispersion coefficient
w = settling velocity/ positive downward
G = growth rate
D = death rate
P(z,t) = phytoplankton biomass
In the absence of vertical transport, the population of each depth
would either grow exponentially if G(z) > D(z) or decay exponentially
if G(z) < D(z). In fact, for E = w = 0 it is clear that the solution
is of the form:
P(z,t) ^ exp [G(z) - D(z)]t (57)
so that G-D is the net growth or decay rate of the population in a
kinetic reactor with conditions comparable to those at depth z.
I r -J \
The magnitude of G is quite well know v ' from numerous short term
growth experiments. For the case of excess nutrients and constant
103
-------
optimal light intensity, G is in the range of 1.5-2.5 day
at 20°C. The decay rate, D, for a population isolated from a
light source is on the order of 0.05 to 0.2 day~l at 20°C .
Thus, in the absence of vertical transport and with excess
nutrients present, there should exist a depth at which the
population is growing at a rate on the order of 1.0 day so
that in ten days the population at this depth increases by
approximately 22,000 times. At large depths, however, the
population should, in ten days, decrease to approximately 37%
of its initial concentration. Such explosive growth does
occasionally occur (patch blooms) but the more normal course of
events is a much slower growth of the entire population until
either nutrient exhaustion, zooplankton predation, or self-shading
reduces G-D to a point at which population growth ceases. The
question of interest is how does the population growth rate depend
on the kinetic parameters, G and D, the transport parameters, E
and w, and the variation of growth rate with depth due to light ex-
tinction.
If conditions are such that the population does increase with time,
and for a period the transport and kinetic parameters are constant
in time, then there are mathematical reasons,given subsequently,to
expect that after a relatively short time, the population distribu-
tion will be of the form
P(z,t) - P(z)eyt C58)
that is, the entire population will be growing at an asymptotic
population log growth rate, y, independent of depth. The
asymptotic population growth rate and the asymptotic population
distribution as a function of depth, P (z) , aside from a
constant multiple, are independent of the initial biomass
104
-------
distribution P(z,0) at the time of the start of the population
growth (arbitrarily taken as t=0). They depend only on the
kinetic and transport parameters which characterize the period
of population growth. This is a substantial simplification of
what, in general, are complex and difficult solutions to the
biomass equation, and from this simplification follows sub-
stantial insights into the behavior of phytoplankton populations
during their growth phase.
ASYMPTOTIC POPULATION GROWTH RATES
To establish that, in fact, phytoplankton biomass behaves as pre-
dicted by the mathematical analysis, a series of plots of log
P(z,t) vs. time at various depths are presented in Figs. 29
through 32.
During the winter and spring the horizontally averaged phytoplankton
biomass in Lake Ontario increases quite slowly but the increase is
approximately exponential as shown in Fig. 29. The data collected
during the 1973 IFYGL cruises shows this population growth and in
particular, it is approximately the same for all depths to 50
meters. Since the population growth rate is quite small, and the
time for the development of the asymptotic distribution is on the
order of 1/y, in this case approximately 100 days, it is not
surprising that the lower layers have not yet responded in the
predicted manner.
During the intensive Project Hypo study of Central Lake Erie (8)
the phytoplankton population increased approximately five fold in
35 days. As shown in Fig. 30 the asymptotic population growth
rate was y^0.043 day" and again the population increased at this
105
-------
= 0.013 (day"1)
I FEB'MAR'APRI FEB'MAR'APRI FEB'MAR'APRIFEB 'MAR'APR!
o
cr\
5.0
2.0 •
| FEB'MAR'APRl FEB'MAR'APRl FEB ' MAR'APRl FEB'MAR'APRl
Figure 290 Lake Ontario Phytoplankton Biomass, horizontal averages 1973.
IFYGL
-------
M = 0.043 (day"1)
>-
Q_
210 230 250 230 250
190 210 190 210
TIME, days
230 250 230 250
190 210
Figure 30. Central Lake Erie Phytoplankton Biomass, four station log averages. 1970.
Project Hypo
-------
exponential rate at 1,12 and 17-20 meters, with the population
in the deeper layers (21-25 m) acquiring this growth rate near
the end of the sampling period. The time to achieve the
asymptotic distribution in this case is on the order of 25 days.
2
Onondaga Lake is a small (12km surface area, 20 m maximum
I Q \
depth), highly eutrophic lake in New York Statev . The
duration of the growth period for the phytoplankton population
is 150 days during which the population increases one thousand
fold to concentrations of 10 cells/ml corresponding to a
population growth rate of y=0.048 day . As shown in Fig. 31
the population growth rate is approximately the same at the four
depths sampled even though there is almost an order of magnitude
less biomass in the deeper layers.
St. Margaret's Bay is a coastal inlet on the southern shore of
Nova Scotia . Similar behavior is observed during the period of
exponential growth. The population growth rates is 0.16 day" ,
substantially larger than those of the previous examples. The
deviations from the predicted constant population growth rates in
the lower layers cannot be attributed to the time it takes to
achieve the asymptotic distribution since at this large y, 10 days
should be quite sufficient. It appears that self-shading is
causing the growth rate in the lower layers to decrease as the
population increases in the top layers, thus violating the assumption
that the parameters are constant in time.
Taken as a whole, the data presented strongly suggests that during
the exponential growth of a population, if the kinetic and trans-
port parameters are constant, an asymptotic population growth rate
occurs which characterizes the growth of the entire population,
108
-------
-ll
= 0.048 (day'1)
CO
8 104
Q_
O
io
IO2
0 100 200 300 300 300 300
0 100 200 200 200
0 100 100
0
TIME, days
Figure 31«, Onondaga Lake Phytoplankton Biomass, two stations.
1969.
109
-------
M =0.16 (day'1)
100
TIME, days
Figure 32. St. Margaret's Bay Phytoplankton Biomass, station A, 1969.
-------
independent of position.
SOLUTIONS FOR SPECIFIC CASES
The existence of an asymptotic population growth rate, independ-
ent of depth for temporally constant parameters, is a consequence
of the properties of the general mass balance equation. In order
to investigate the relationship between this growth rate and the
parameters of the equation, it is necessary to make certain
assumptions concerning the variations in depth of the kinetic and
transport parameters. In particular, assume (not very realistically
for some cases) that all the parameters, save the growth rate, are
constant. For the growth rate, two forms lend themselves to ana-
lytical solutions:
a two layer approximation:
G(z) = G 0 < z < L
= 0 L < z
and an exponential approximation:
G(Z) = G e'Z/L (60)
The solution proceeds directly using the asymptotic solution (eq. 58 ).
Substituting it into the conservation of mass equation (56 ) yields:
_E d_P + wdP = [G(z) - (D+u)]P (61)
This equation is solved together with the boundary condation at the
water surface:
dP <62)
-E Sif- + wP = 0 at z=0
dz
which requires that no biomass cross the air-water interface, and
111
-------
the boundary condition at infinite depths, for which we reauire
that the biomass remain finite. The result is an eigenvalue
problem to find y and P(z) for the above ordinary differential
equation with the homogenous boundary conditions.
For the exponential case the differential equation (56)becomes:
-E £-?• + w P. = [Ge"Z/L - (D+y)]P (63)
which has a solution
P(z) = e- J (2X e'Z/2L) (64)
/2 2
where v = / IT + 4A /<
K = G/(D+y)
TT = WL/E
X = L /G/E
and J (Z) is a Bessel function of the first kind. At large
depths, the solution becomes
A. 2
P (z) = exp [(1 - /I + _) TTZ, (65)
which approaches zero exponentially. The eigenvalue equation is
obtained from the boundary condition at z = 0 (eq. 62) :
(2A) (66)
A multiplicity of v's are solutions for the equation; the one of
interest is that which is largest and, more important, positive
since it is the asymptotic population log growth rate. For
112
-------
certain sets of TT and X no positive y exists. For such
situations no population increase is possible.
For the two layer case, the solution is obtained in two parts,
a sinusoidal solution in the euphotic zone, 0 < z < L, and a
decreasing exponential thereafter. The eigenvalue equation is:
-1 ( - } -1C / l+4n)
IT = tan v 3- + tan v — (67)
3
3/2
where 3 = /4n(<-!)- 1
n = (D+y) E/w2
and K = G/ (D+U)
The solutions of these eigenvalue equations is shown in Fig. 33
in terms of the dimensionless ratio of growth rate of the sum of
death rate and population growth rate G/D+y; a penetration depth
L/ GYE which characterizes the size of the euphotic zone in terms
of the growth rate and the dispersion coefficient; and a piclet
number wL/E, which is the ratio of advective transport to dis-
persive transport. The values for the exponential case are given
in Table 10. It is interesting to note that the shapes and magnitudes
of the results are comparable, indicating that the details of the
vertical distribution of the growth rate are not critical to the
conclusions that can be drawn. As an example of their utility con-
sider the effect of settling velocity for the range of G/D+y en-
countered in practice, G/D+y < 10, it is clear that a peclet
number, wL/E, of greater than 0.1 is required to affect the growth
rate of the population.
113
-------
Two-Layer Constant Growth Rate
Exponentially Decreasing Growth Rate
1.8 -
o
o
,—-*
s
- 0.6 -
-1.0 -0.5 0.0 0.5
1.8 -°\
1.4 h
1.0 -
0.6 h
0.2 -
1.0 -1.0 -0.5 0.0 0.5
log10(LVG7E)
1.0
Figure 33. Eigenvalue equation solutions for constant growth layer and exponentially
decreasing growth rate cases.
-------
TABLE 10
Log1Q (G/D+p) FOR EXPONENTIALLY DECREASING GROWTH RATE
Log
TT *u -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1
0.50 2.2281
0.20 1.7007 1.0358 0.7731
0.10 1.4158 1.0455 O.P251 0.6681
0.05 1.9903 1.4323 1.1192 0.9129 0.7536 0.6263
0.02 2.3052 1.7515 1.4364 1.2021 1.0115 0.8519 0.7172 0.6038
0.01 2.2545 1.8509 1.5721 1.3456 1.1512 0.9816 0.8337 0.7059 0.5966
0.00 2.0170 1.S22C, 1.6325 1.4484 1.2724 1.1068 0.9542 0.8166 0.6951 0.5897
Log10
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 p. a Q.9
5'°° 1.2401 0.5484 0.3512 0.2479
3t°° 0.7034 0.4394 0.3106 0.2322 0.1795
2«°0 1.4192 0.6510 0.4352 0.3186 0.2440 0.1920 0.1540
1'°° 1.4778 0.7376 0.5114 0.3842 0.3002 0.2401 0.1950 0.1602 0.1327
°'50 0.9140 0,6411 0.4880 0.3P56 0.3113 0.254R 0.210R 0.1754 0.1469 0.1235
°'20 0.6083 0.4908 O.M020 0.332R 0.2771 0.2322 0.1953 0.1P49 0.1395 0.1183
°'10 0.5486 0.4545 0.3789 0.3175 0.2670 0.2253 0.1906 0.1615 0.1372 0.1167
°.°° 0.4994 0.4228 0.3580 0.3034 0.2575 0.21R7 0.1860 0.15P3 0.1.349 0.1150
-------
It is this kind of order of magnitude analysis for which
the analytical expressions are most suited. For more realistic
investigations of detailed properties of the population, resort
to direct numerical methods is usually necessary.
THE SINGLE LAYER APPROXIMATION
It is common practice in one layer models (12) of the euphotic
zone to approximate the effect of the settling velocity on the
population by comparing w/L to the growth rate G. This approxi-
mation comes from the resulting differential equation for the
euphotic zone layer:
II- = (G - D - w/L)P + E' (P,-P) (68)
at -1-
for E1 the bulk mixing coefficient and P-j^ the hypolimnion con-
centration. It is clear that in this approximation, changes in
the w/L can be exactly compensated for by changes in G. To check
the validity of this approximation against the continuous solution
consider Fig.34, a plot of w/L vs. G, both normalized by D+y. The
exponentially decreasing growth rate case is used. It is clear
that no positive asymptotic growth rate exists if w/L>G so to
this extent the approximation conveys the correct impression.
In order to interpret this figure in more detail, assume that D+y
and L are fixed, and the plot is w/L versus G for varying dis-
persion coefficient. For large vertical dispersion, L/" D+y/E is
small and the solution is independent of w/L, or, put another way,
w/L can change without demanding a compensating change in G to
keep v constant. What is occuring is that the settling biomass'
is being rapidly remixed into the euphotic zone so that the loss
116
-------
2.0
1.6 -
1.2 -
0.8 -
0.4 -
0.0 t
0.4 -
0.8 -
0.2 0.6 1.0 .1.4 1.8
Figure 34. Settling velocity-growth rate interaction.
The effect of dispersion. Exponentially
decreasing growth rate case.
117
-------
via this mechanism is small relative to the loss by dispersive
mixing. On the other hand, for a small vertical dispersion the
curves bend to such an extent that their slope is greater than 1,
indicating that for this situation it requires larger than pro-
portional increases in G to compensate for increases in w/L. The
meaning of this unexpected result will become clear subsequently.
In order to investigate these interrelationships in a systematic
way a series of plots is presented to indicate how each parameter
affects y while holding others constant. The exponentially de-
creasing growth rate is employed for these calculations.
THE EFFECT OF INCREASING GROWTH RATE
For various euphotic zone depths (assuming the transport parameters
w and E are constant) the effect on y of increasing G is shown in
Fig.35. It is immediately clear that no solutions are possible for
which y>G-D, which is obviously true since such a situation would
require that mixing and settling have no effect on the population
and that the fact that G(z) decreases in depth is unimportant.
Perhaps a more unexpected result is that the slopes of the curves
are greater than one, indicating that a change in G causes a larger
than proportional change in D+y. As G is further increased, how-
ever, the population growth rate asymptotically approaches G-D as
expected, and since there are non zero values of G for which D+y
is zero, it is reasonable that the slopes are greater than one to
enable D+y to, in effect, catch up with G. The condition G= D+y
occurs at smaller G for deeper euphotic zone, as expected.
118
-------
3.
+
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Figure 35. Population growth rate versus maximum growth
rate. The effect of Euphotic zone. The ex-
ponentially decreasing growth rate case.
119
-------
THE EFFECT OF INCREASING SETTLING VELOCITY
For various euphotic zone depths (assuming G and E are constant)
the effect on U of increasing w is shown in Fig. 36. As w is in-
creased the population growth rate is not affected until a critical
range is reached at which time the effect begins. Further in-
creases cause a sharp drop off in P until population growth is no
longer possible.
THE EFFECT OF INCREASING DISPERSION
The effect on y of increasing the dispersion coefficient, E, for
various euphotic zone depths (assuming G and w are constant) is
shown in Fig. 37. A most surprising result is that for the range
2
of the parameter, EG/w below one and LG/w < 10, increasing the
dispersion increases the asymptotic growth rate. At first glance
this appears to be incorrect since the mixing causes the popu-
lation in the euphotic zone to be transported out of the zone, on
balance, and therefore, would tend to decrease population growth.
However, biomass is also being advected out of the zone via the
settling velocity. With a small dispersion coefficient there is
very little biomass at the surface since there is no way for the
population, which grows as it sinks, to return to the surface
layers. Hence, increasing the mixing actually benefits the
population by allowing more biomass to accumulate in the surface
layers where the growth rates are largest. It is this phenomena
which causes the slopes greater than one in Fig. 34 for small dis-
persion coefficients.
THE EFFECT OF BOUYANCY
It is well known that certain species of phytoplankton exhibit
negative sinking velocity and an in situ measurement has con-
firmed an example so that instead of a liability their-
120
-------
0.0
+
o
-1.0 -
-2.0
-2.0
-1.0
logi0(w/vfGE)
0.0
Figure 36. Population growth rate versus settling velocity
the effect of Euphotic zone depth. The exponentially
decreasing growth rate case.
121
-------
0.0
o
-1.0 -
s
8*
-2.0
0.0
1.0
logio
-------
1-0
LO
o
s
g>
1.0
0.8 -
0.6
0.4 -
0.2 -
wL/E 0
lO "T" "~~ "~~ ~T ~~~ ."""" "T" ~"~
l__l I I—I
' I 1_J 1
-1.0
-0.5
0.0
0.5
1.0
Figure 38. The effect of negative sinking velocity. The
exponentially decreasing growth rate case.
-------
bouyancy is an asset which should increase their population
growth rate. Fig. 38 presents the results. For negative
peclet numbers with magnitude less than 0.2, not much increase
in y is calculated. However, as [w|L/E increases further, the
population growth rate approaches the maximum possible, i.e..
y -»• G-D, independent of Lv G/E as indicated by the horizontal
lines in the figure. The effect of the upward advection is to
pile the biomass at the surface where it can grow at the kinetic
rate, G-D, independent of the dispersive forces tending to mix it
downward. This effect occurs for the negative peclet numbers
greater than 1.0, provided L/ G/E exceeds 0.3.
PRELIMINARY APPLICATIONS
As a first step in exploring the utility of the preceding an-
alysis, a series of representative values for the dimensionless
parameters are shown in Fig. 39. The majority come from the coastal
ocean with txvo lakes included. It is interesting to note that
the values of G/D+y are comparable whereas the value of L/ G/E
span an order of magnitude. The major difficulty with increasing
the number of points on this graph is in estimating the vertical
dispersion coefficient and settling velocity. Values for G, D and
L for a given temperature, light intensity and extinction coefficient
can be approximated; y is obtained from population observation,
vertical dispersion can be obtained from vertical temperature models
in which case the analysis, to within the accuracy of the assumption
of constant parameters, can provide vertical settling velocities.
It is interesting to note that these data, when plotted on the figure
which investigates the effect of dispersion, are in the range of the
124
-------
1.0
"- 0.8
+ 0.6
°" 0.4
0.2
0.0
WL/E
a SARGASSO SEA
a FLORIDA STRAITS
BMONTAUK POINT
AGEORGES BANK
OONONDAGA LAKE
x LAKE ONTARIO
J _ I
-1.0 -0.5
0.0
0.5
log10(LVG/E)
1.0
+
o
0.0
-1.0
0.0
10
LG/w
0.5 1.0
logio(EG/w2)
1.5
2.0
Figure 39. Application to Lakes and Oceans the range of
biologically relevant dimensionless numbers.
125
-------
maxima for the oceanic situations and in the decreasing regions
for the shallower situations. Whether this indicates some
adaptive mechanisms for the characteristics of the marine
phytoplankton is an interesting speculation but it is probably
beyond the precision of the analysis to persue this point.
SUMMARY OF MATHEMATICAL FACTS
The properties of eq.C56) are best illucidated by examining a
spatially discretized version.
Let P.(t) = P [ (i + h) Az,t] for a finite set of points starting
at Az/2 and spaced Az apart thereafter. The spatially discrete
version of the differential equation follows by using the ap-
proximations :
- 2 P. + P.
C69)
AZ2
P. P.
Az
for w > 0
(70)
Pi+l ~ Pi for w < 0
There results a set of ordinary differential equations which can be
written in matrix form as:
at
126
-------
Two properties of A are important:
(i) A is essentially nonnegative, i.e.
3ij ' ° ± * j 021
(ii) A is irreducible, i.e. for any i,j there exists
a finite sequence of indicies k (0) = i, k(l),...,
k(r) = j such that a^/j^D k/n) 7* ° for h=l,...,r.
Property (i) follows directly from the fact that the off diagonal
2
elements of A are made up of sums of E/^z and w/Az which are posi-
tive. Property (ii) essentially requires that mass from any seg-
ment i can eventually be transported to any segment j , which is
clear from the geometry of any reasonable problem.
The theorem of interest is:
The matrix A satisfying (i) and (ii) has a unique
(up to a multiplicative constant) strictly positive
eigenvector 4>.., with real simple eigenvalue V,, i.e.
A*l = yi*l (73)
and V, is the largest eigenvalue in the sense that
yi * Re [V (74,
for all the other eigenvalues, y., of A.
It follows directly from this fact that the asymptotic solution of
the differential equations (71 ) has the form
P (t) = k eylt <|»1 + o(eyt) (75)
127
-------
where
y, > y >max Re [y .]
and k is a constant.
This fact can be proven directly . However, to see what is
occurring assume that all eigenvalues y . are distinct so that an
~* (15)
eigenvalue expansion of the solution exists of the form.
P(t) = 2.. k.. ev^ $.. (76)
where k. are constants related to the initial condition P(z,0).
If y, is positive and Re[y.] negative for all other j, then clearly
the positive eigenvalue part of the solution dominates after a time
on the order of 1/y-,. If both y, and y- are positive (and y~ real
for simplicity) the solution in segment i approaches:
P.(t) «
(77)
kl*li
but y^"-^1-! so that the exponent is negative and it becomes small
relative to 1 after a time on the order of I/ (y-^-yg) and again
the maximum eigenvalue solution dominates.
128
-------
SECTION VII
REFERENCES
1. Riley, G.A., H. Stommel and D.F. Bumpus. "Quantitative
ecology of the plankton of the western North Atlantic,"
Bull. Bingham Oceanog. Coll. 12(3):1-169. 1949.
2. DiToro, D.M., D.J. O'Connor, R.V. Thomann and J.L. Mancini.
"Preliminary phytoplankton, zooplankton nutrient model of
western Lake Erie," in Systems Analysis and Simulation in
Ecology, Vol. 3. Academic Press, New York. 1973. In press.
3. Canale, R.P., S. Nachiappan, D.J. Hineman and H.E. Allen.
"A dynamic model for phytoplakton production in Grand
Traverse Bay," in Proc. 16th Conf. Great Lakes Res.,
pp.21-23. International Association, Great Lakes Research.
1973.
4. Thomann, R.V., D.M. DiToro, D.J. O'Connor and R P. Winfield.
"Mathematical modeling of eutrophication of large lakes,"
in Annual Report - Year 1. Manhattan College, Bronx, New
York. 1973.
5. Kierstead, H. and L.B. Slobodkin. "The size of water
masses containing plankton blooms," J. Mar. Res. 12:141-147.
1953.
6. DiToro, D.M. , D.J. O'Connor and R.V. Thomann. "A dynamic
model of the phytoplankton population in the Sacramento-
San Joaquin Delta," Advances in Chemistry, No. 106, pp.131-180.
American Chemical Society, Washington, D.C. 1971.
7. Eppley, R.W. "Temperatures and phytoplankton growth in the
sea," Fishery Bull. 70(4):1063-1085. 1972.
129
-------
REFERENCES
8. Burns, N.M. and C. Ross. Project Hypo., CCIW Paper No. 6
and EPA TS-05-71-208-24. U.S. Environmental Protection
Agency, Washington, B.C. February 1972.
9. O'Brien and Gere Engineers Incorporated. Onondage Lake
Study. U.S. EPA Publication 11060, SAE 4/71, U.S.
Government Printing Office. 1971.
10. Platt, T. and B. Irwin. Primary Productivity Measure-
ments in St. Margaret's Bay. Tech. Report No. 203. Fish.
Res. Bd. of Canada. 1973.
11. Murphy, G.M. Ordinary, Differential Equations and Their
Solutions. D. Van Nostrand Co., Princeton, N.J. 1960.
12. Steele, J.H. "Plant production on fladen ground," Jour.
Mar. Bio. Assoc., U.K. 35:1-33. 1956.
13. Burns, N.M. and A.E. Pashley. "In situ measurements of
the settling velocity profile of particulate organic
carbons in Lake Ontario," Jour. Fish. Res. Bd., Canada
31(3):291-297. 1974.
14. Birkhoff, G. and R.S. Varga. "Reactor criticability and
nonnegative matrices," Jour. SIAM 6(4):354-377. 1958.
15. Bellman, R. Introduction to Matrix Analysis. McGraw-Hill,
New York. 1960.
130
-------
SECTION VIII
PRELIMINARY LAKE 3 MODEL
.As shown in Fig. 40,the Lake 3 model represents a synthesis of
spatial geometry and transport structure with the biological and
chemical kinetic interactions of the Lake 1 and 2 models. Lake 3
is therefore a three-dimensional time variable model which does
permit some analysis of phenomena in the near-shore region as
compared to the open lake region. The structure of Lake 3
represents a further step beyond the simple Lake 1 model, the ex-
tent of the step incorporating such factors as computational time
for a one year run and availability of survey data.
The work reported here is of a preliminary nature since considerably
more effort is required to assure the veracity of the computed
results both in terms of the numerical output and the degree to
which the model verifies observed data. The full detailing and
documentation of Lake 3 will be the subject of a future report.
The basic segmentation used for the Lake 3 model is shown in Fig.40.
As shown, 67 segments are used distributed over five vertical layers
The upper two layers, 0-4 meters and 4-17 meters are considered to
represent the epilimnion during periods of vertical stratification.
A "ring" of segments extending 10 km out from the coast is used to
represent the near-shore environment. All cross-sectional areas
for interfaces in three dimensions between the segments were plani-
metered and volumes were computed and summed to within 10% of total
lake volume.
CIRCULATION AND DISPERSION
The emphasis in Lake 3 is on the biological and chemical kinetics
131
-------
to
Ml
>150 Meters
Figure 40 0 Lake 3 segmentation
-------
rather than the hydrodynamic circulation although the inter-
action between the kinetics and transport is important. The
hydrodynamic circulation is therefore considered to be ex-
ternally supplied through inferences from observations and/or
hydrodynamic model output. An extensive literature exists on
the circulation of Lake Ontario and relatively sophisticated
123
mathematical models have been developed ' ' . A review of the
complete development is given by Simons . It is not intended
in this report to provide.a detailed review of the circulation
in the Lake and the various analytical and numerical model
schemes that have been used to estimate the circulation.
Rather, the purpose here is to simply report on the circulation
that was used in some of the preliminary Lake 3 phytoplankton
runs and the approximate sensitivity of phytoplankton dynamics
to changes in the assumed circulation patterns.
This disucssion will therefore outline the procedure used for
establishing circulation patterns and flow exchanges in the 67
segment - 5 layer, Lake 3 model. Two flow balances are pre-
sented, one corresponding to "summer circulation" and the other
to "winter circulation".
It is known that Lake Ontario has two distinct net surface cir-
culation patterns. During summer and fall the prevailing winds,
the dominant factor governing water circulation, are from the
west and southwest and one distinct circulation pattern is set
up ("summer circulation"). During winter and spring, the pre-
vailing winds are from the west -and northwest and a different
circulation pattern is established ("winter circulation"). The
summer and winter circulations in Lake Ontario are associated
133
-------
with the stratified and isothermal temperature states of the
water, respectively.
4 5
The International Joint Commission and the E.P.A. have
published summer circulation patterns for Lake Ontario. Mean
monthly resultant velocity vectors that were observed and used
by EPA for establishing their circulation patterns were obtained
from them for use in assigning velocity values to the general
circulation vectors. Wide variance in magnitude and direction
over each sampling station was encountered however, even to the
point where resultant velocity direction did not agree with the
general circulation pattern. Typical stations were also chosen
where values were assumed to be reflective of the average
observed magnitude and direction. This procedure, however,
yielded a poor flow balance and the scheme was discarded in favor
of using the general values assigned to the summer circulation
pattern by IJC^ as a guideline.
The water budget information contained in the IJC report was used
for establishing tributary inflows to the various shoreline seg-
ments. It was decided to split the major inflow to Lake Ontario
from the Niagra River (195,000 cfs) over the first and second
vertical layers (0-4, and 4-17 meters) while restricting all other
tributary inflow (Oswego River, Genesee River, 12 mile Creek,etc.)
to the first layer (0-4 meters). The Niagra inflow was split
volumentrically between segments 5 and 31 each receiving 4/17 and
13/17 of the flow, respectively.
Having established tributary inflow to the segments, segment flow
balances were established and intersegment velocities were back
calculated to check for agreement with the general values assigned
134
-------
by IJC. The final circulation pattern, as given by inter-
segment velocity values can be seen in Fig.41. It should
be noted that since the Niagra River inflow was split be-
tween the first two layers only, there is no inflow and
hence can be no outflow from the third layer. This is a good
approximation since the average depth of the northeastern most
segment (26 and 52 in top two layers), the outflow segment,
is only 20-25 meters while the third layers extend from 17 to
50 meters. In any event the bathymetry dictates that there
would not be much outflow in the third layer. No flows were
assigned to segments in the fourth layer since with only three
adjoining segments any inflow to one of the segments from an-
other would have to be countered with the same back flow for a
flow balance. The same applies to the fifth layer.
For the winter, the EPA general circulation pattern correspond-
ing to winds out of the west was used for establishing a flow
balance. Bonham-Carter's wind driven numerical model outputs
a fairly similar circulation for a west wind. The model, how-
ever, indicates easterly flow along the north shore, while the
EPA circulation pattern indicates upwelling along the north
shore. Since both agree the general transport is to the east,
there would have to be return flow or upwelling and since the EPA
circulation pattern indicates areas of upwelling it was decided
to use their general circulation pattern as a guideline for
establishing a winter flow balance. Later discussion will in-
dicate how Bonham-Carter's modelled circulation was set up to
this segmentation scheme. The balanced velocity vectors are
135
-------
4
N-
O N T A
R|0
44°
TORONTO
43<
Lake Erie
10 0 10 2t) 30 40 50
10 0 10 20 30
_L
79° 78° 77°
Figure 41. Assumed summer circulation for Lake 3 model,
-------
shown in Figure 42 and represent the assumed circulation for
the winter condition. It should be noted that discrepancies
exist in the literature in regard to Lake Ontario circulation,
and that many of the numerical circulation models are awaiting
IFYGL current measurements for verification. However, the
Lake Ontario circulation must be calculated from the smaller
(approximately 5km) grid on to the larger grid of the phyto-
plankton model. As additional information becomes available,
adjustments to the circulation structure can be made.
Time and space varying horizontal and vertical dispersion co-
efficients are incorporated to simulate some of the more pre-
dominant features such as the thermal bar and vertical strati-
fication. Horizontal coefficients are constant for the winter
and are set to zero for the stratification period to simulate
the thermal bar effect resulting in restriction of the mixing
of near shore waters with main lake waters. Vertical dispersion
is also time varying in Lake 3 in a fashion similar to Lake 1 and
is constant in the winter, lower and depth varying in the summer,
and evaluated on the basis of a straight line function during the
transition period.
Fig,. 43 illustrates the time varying horizontal and vertical
exchange coefficients to be used in this model.
137
-------
c
oo
-N-
O N T A
RlO
44C
TORONTO
43°
HAMILTON
Lake Erie
\
-
10 0 20 30 40jO KM
m_j_a_a-gM.
79° 78° 77°
Figure 42. Assumed winter circulation for Lake 3 model
-------
9000
"E6000
u
^3000
I I I I L
HORIZONTAL EXCHANGE
NEAR SHORE-NEAR SHORE
MAI NUKE-MA IN LAKE
i I I i 1—
^,9000
-------
VALIDITY OF ASSUMED CIRCULATION
Since the preceding circulation represents only the crudest
estimate of flow transport in Lake Ontario, analyses must be
made to determine the validity (or lack thereof) of the
assumed circulation structure. The strategy therefore at
this stage of the Lake 3 model is to input a circulation
pattern into the kinetic structure that is at least reason-
able in its broadest outlines. Since the segments of the
Lake 3 model are large (10-4Okm), much of the fine scale
detail of the circulation will not be reproduced. Indeed,
one of the functions of this effort at this time is to
determine the degree to which a detailed circulation is re-
quired for the phytoplankton dynamics.
Initial Lake 3 model runs therefore were aimed at verification
of transport regimes using water temperature as an indicator.
The analysis, while crude, served the purpose of providing some
estimate of the validity of the transport regime that was
assumed. An annual temperature function is inputted to the top
layer of the lake as a time variable forcing function, and the
hypothesized flow and dispersion regimes transports this heat
throughout the entire lake, computing the time variable temper-
ature response in each of the segments. This response is then
compared to observed regional thermal properties of Lake Ontario.
The temperature input used is based on the work of Rodgers and
Sato for IFYGL, who have computed values of monthly heat flux
terms for 1965 thru 1968. Though discrepancies between their
results and observed data are reported,the computation of a
total heat storage term, Q , proves to be ideal for coupling
to the Lake 3 model and permits inputing temperature as a lake
140
-------
surface forcing function.
Using Lake Ontario meteorological data and empirical relation-
ships, these heat storage values were computed using the
formula:
Qt= Qs -Qr -Qb -Qh -Qe (78)
Q is the heat storage; Q , the incident solar radiation;
t S
Q, , the net back radiation; Q, , the heat transfer to atmosphere
2
and Q is the evaporation energy where all units are cal/cm -day.
With the monthly heat flux terms, values of Q. were calculated
for the five year period studied. These values are listed in
Table 11 and the five year average values were used in the Lake 3
runs.
Once Q. , the total heat flux term or total heat storage in the
lake, is known, then
Ts (t) = Qt/H (79)
where T (t) is the temperature forcing function for the surface
s 3
layer (°C/day or cal/cm -day) and H is the depth (cm) of the first
segment.
This computed temperature forcing function was inputted to the
top layer segments (0 - 4 ) of the Lake 3 model, and was then
advected and dispersed throughout the model by the hypothesized
flow and dispersion regime. The computed time variable temper-
ature response in each of the segments was then compared with a
regional analysis of Lake Ontario thermal properties.
The degree of data manipulation involved to manipulate a transport
regime in this 67 segment model, with its 187 dispersion inter-
faces and its 142 flow interfaces, is extensive so that a trans-
port regime based on verification against generalized temperature
141
-------
TABLE 11 -
MONTHLY HEAT FLUX VALUES
MONTH
JANUARY
FEBRUARY
MARCH
APRIL
MAY
JUNE
JULY
AUGUST
CT?"D*nt?M't3'C|HD
br*F lilrlDCiK
DC 1 OJBri R
NUV.bMB.LR
TMP/^T^TUTOPn
D£iC£»Mi3r*K
* (cal)
2
(cm day)
1965
V
-329.1
-2l6.6
4.6
275.2
533.2
463.0
270.7
67.0
A 1 1
*> "3 C O
— ^ -3D . 0
o Ark £
zUU . b
O9£ 1
ZZO . 1
1966
Qt
-392.0
-139.5
80.3
207.0
413.7
610.9
307.0
69.5
cc o
DO . ^
C. A A
~ b4 . 4
ICO >1
IDZ . 4
I/I Q 1
j^y . j.
1967
Qt
-214.2
-278.1
48.4
318.8
420.6
618.6
303.0
51.2
99 "7
Z . /
1 ft O T
— J.UO . /
O t^/l Q
/_>4 . O
OQC Q
/yo . o
1968
Qt
-329.2
-244.9
88.4
361.1
353.7
471.7
293.3
198.5
i r\o i
-LU^C . J.
yl c o
— 4o . y
OftC "7
OT n o
1969 5 year
Q. average
-275.5 -308.0
-153.6 -206.54
11.3 46.6
313.1 295.04
448.1 380.8
475.4 464.8
293.5
96.55
— 97 77R
_ ___ "1 1 ^
-------
properties rather than one verified against one year of data
with its specific meteorological conditions, was considered to
be more feasible.
Based on these considerations, the verification medium chosen
g
was the temperature characterization of Lee . Composite tempera-
ture data for the years 1960 - 1969 were used and distinct
horizontal variations of water temperature in Lake Ontario were
noted. Lee distinguishes seven distinct regions of the lake
where thermal properties seemed to be seasonally consistent on
a year to year basis. This regional division of the lake, seen
in Figure 44 is similar to the underlying basis of surface layer
segmentation for Lake 3, namely a near shore ring around the
lake and large main lake segments. This similarity coupled with
the fact that Lee's data are based on a generalized analysis of
ten years of data prompted the use of this data set for tempera-
ture verification of transport regimes in the Lake 3 model.
Several different regimes were investigated with emphasis on the
vertical and horizontal dispersion characteristics. The degree of
verification of the computed output as compared with Lee's data
is presented in Figs.45,46. Table 12 lists the equivalence between
Lee's regions and the Lake 3 segments, and the values which are
overplotted against Lee's data are the means of the appropriate
Lake 3 segments at times corresponding to the median cruise dates
reported by Lee.
Figure 45 also shows the effect of several sensitivity runs. The
effects of using vertical dispersions only with no horizontal
dispersion or flow, and the effect of reversing the direction of
143
-------
TABLE 12 -
REGIONAL TO LAKE 3 SEGMENTATION EQUIVALENCE
Lee (region)
1
2
3
4
5
6
7
SEGMENT OF LAKE 3
1
26
1,2,4,6
10,14,18
21,22,23,
25
3,7,8,11,
12
15,16,19,
20,24
5,9,13,17
2
52
27,28,30,
32
36,40,44
47,48,49,
51
29,33,34,
37,38
41,42,45,
46,50
31,35,39,
43
3
60
53,55,56
56,58
60,62
54,57
57,61
55,59
4
63,64
64,65
5
66
66,67
144
-------
I
Ul
79(
78<
IT
o
Figure 44. Regional division of Lake Ontario by Lee°.
-------
TEMPERATURE, deg C
10 15 20
201
40 H*
60 P
80 r Rl
100 1 !LLJ
u
20
40
60
80
100
J
-
•
1 1 T 1
i
R2
u
20
40
60
80
inn
r ' '
-
'I R3
0
20
40
60
80
inn
4i ' ' '
f
-
-1
Q.
UJ r
0 n1
U
20
40
60
80
inn
) 5 10 15 21
: 'IJ '
Rl
5 10 15 20
V
20
40
60
80
inn
i .
>
• l\
)
R5
u
20
40
60
80
inn
0<
20
40
60
80
inn
n
20
40
60
80
inn
1UU
Oc
20
40
60
80
inn
J
)
-
3
-
)
-
R2
5 10 15 2(
ii ' ' '
I
6 R6
5 10 15 2
1 . I I
1
R2
5 10 15 2(
i \ i i
1
/i R6
u
20
40
60
80
20
40
60
80
100
20
40
60
80
100
0
20
40
60
80
inn
I R3
) 5 10 15 2(
B ' ' '
• i
R7
3 5 10 15 2
R3
) 5 10 15 21
1 ll'.l ' '
" '/
" R7
WINTER
5 10 15 20
u
20
40
60
80
inn
I r.
L !'
1
L R4
FALL
LEE HORIZONTAL/VERTICAL DISPERSION + FLOW VERTICAL DISPERSION ONLY
ROW REVERSED
Figure 450 Comparison of computed temperatures from Lake 3 model with data
from Lee°, winter, fall.
-------
TEMPERATURE, deg C
on
40
60
80
100
C
20
40
60
E 80
. 100
x:
°- r
UJ V
O 0
20
40
60
80
100
20
40
60
80
inn
)
,.
.
-
-
-
-
)
_
.
.
-
3
•
1
5 10 15 2
1 ' .'''.••&
^f^^
Rl
5 10 15 2
.s^&P1^
:jf
I' R5
5 10 15 2
' '
7
Rl
5 10 15 2
ll 1 1
R5
0 (
20
40
60
80
ir\n
1UU
D o(
u
20
40
60
80
100
0 (
20
40
60
80
irtft
1UU
0 oc
u
20
40
60
80
inn
)
.
-
-
)
-
-
-
)
.
.
-
-
)
-
-
5 10 15 2
' .r~~~^Z--"'
/**7-^""
R2
5 10 15 2
I 1 .1 // /
I/
If R6
5 10 15 2
!>; '
•b
R2
5 10 15 2
it '
1'
i
! w
0 o{
u
20
40
60
80
inn
iw
0 o1
u
20
40
60
80
inn
0 o(
20
40
60
80
inn
1UU
0 oc
20
40
60
80
im
) 5 10 15 2
' J-^f^--"
. ^^•'••••"
- '
R3
3 5 10 15 2
i ' L^^ .-•
^ — ~^^
'. '
R7
) 5 10 15 2
1 f ' '
I
-
R3
5 10 15 2
1 R7
0
20
40
60
80
inn
0
0 (
n
20
40
60
80
inn
D
3 5 10 15 2(
' ' -^^^
. r£^^'^
-
R4
SUMMER
) 5 10 15 2(
//. i i
it :
'^
-
R4
SPRING
L£E HORIZONTAL/VERTICAL DISPERSION + FLOW VERTICAL DISPERSION ONLY
ROW REVERSED
Figure 46. Comparison of computed temperatures from Lake 3 model
with data from Lee°, summer, spring
-------
the flow field were investigated.
As can be seen, no marked effect on the temperature verifi-
cation results when the flow transport is turned off or when
flow is reversed. Though this appears to be odd, the reason
is clear when one considers the size of the segments in the
Lake 3 model. Though the model is large from a computational
point of view and it does, because of its spatial detail, give
near shore to main lake resolution, the surface areas of the
segments are so large that the interfaces for vertical trans-
port are much greater than those for horizontal and flow trans-
port. Hence vertical dispersion terms dominate.
The grid is not fine enough and effects such as near shore thermal
bar flow restriction and changing wind pattern effect on hori-
zontal transport cannot be investigated with any degree of
accuracy. A grid fine enough such that these effects could be
seen, does not appear to be computationally feasible at this
time.
The dispersion regime resulting from this analysis is presented
in Figure 43. Its underlying features consist of lowering vertical
dispersion in the summer to simulate stratification, lowering
near shore to mid-lake segments horizontal dispersion for a period
in the spring to simulate the termal bar, and simulating stratifi-
cation earlier in the near shore segments than in the main lake
segments. This is the dispersion regime which is used in Lake 3
for the full biological and chemical systems verification runs.
Though this analysis has many shortcomings, it did serve to give
some basis in theory for coupling of flow transport and dispersion
to the kinetics of the Lake 3 model and it appears that this
approach of treating temperature as a conservative tracer is a
reasonable method of analysis.
148
-------
PRELIMINARY PHYTOPLANKTON COMPUTATIONS
Several runs have been made of the Lake 3 model using the transport
and dispersion regime discussed previously and the system kinetics
given in Section v. Tne purpose of these early runs is to
determine the general direction of the output with special ref-
erence to spatial differences in phytoplankton biomass. Consider-
ably more work remains to be done on Lake 3 including detailed
verification analyses and simulation of effects of nutrient reduction.
The primary data source at this point in the analysis is the results
from the IFYGL effort. The cruise stations are mapped on to the
Lake 3 grid and only those stations for which a reasonably complete
set of data were available are included in the mapping. Fig.47 shovs
the location of the data stations and the Lake 3 grid. Some seg-
ments have only one station for comparison while other segments
(primarily near-shore stations) may have three stations for com-
parison. Data retrievals have been completed from the Storet system
and "segment-depth" averages (and other statistics) by month and
record length have been completed. This data set will then form
the basis for the detailed verification analysis.
Fig.48 shows early very preliminary results of phytoplankton biomass
at three segments across the lake. This run did not include any
vertical settling velocity. Segment #17 represents the Rochester
enbayment and as seen reaches peak values of almost 15 yg chloro-
phyll/Jl at the end of May while the more open lake area (segment
#16) lags by about one-half month. Peak values at the open lake
segment are about half of the Rochester segment. At the end of May
therefore, at day 150, a lateral gradient of chlorophyll of about
149
-------
LEGEND
® IFYGL Station
IFYGL Index Station
ONTA"°
44°
TORONTO
43°
10 0 10 20 30 40 50
79<
78*
77°
Figure 47. Principal IFYGL stations and the Lake 3 grid,
-------
en
t_>
25
I 20
o
"15
O
10
I 5
I 0
Segment 14
Segment 16
Segment 17
14
I
I
I
I
I
I
30 60 90 120 150 180 210 240 270 300 330 360
J FMAMJ JASOND
DAY OF YEAR
Figure 48. Preliminary results of Lake 3,0-4 meters
151
-------
10yg/i is computed. This results primarily from the thermal bar
effect simulated by variations in the horizontal dispersion (see
Fig.43). Gradients across the lake in the fall are not pro-
nounced and reflect the general well-mixed character of the lake
at that time.
Fig.49 explores the vertical variation for two segments during the
month of June and compares the results to data collected during
IFYGL. Maximum biomass occurs in the upper layer and then decreases
with depth. The mean values and range of the computed values
agree quite well with the observed chlorophyll data. In fact, the
agreement is so surprisingly good that further investigation is
warranted to determine why the agreement was obtained so rapidly.
This is especially important when it is recognized that the run
shown in Figs.48 and 49 is for zero phytoplankton settling
velocity.
Output Data Display
The Lake 3 model generates a substantial output totalling thousands
of numbers that represent the time variable response of each
variable at each segment. In addition, output on the growth and
death kinetics, nutrient recycle and nutrient fluxes are also
obtained. The size of Lake 3 makes it difficult to fully compre-
hend the output so some effort is being devoted to further analysis
and display of the output. In this way, it is hoped that further
understanding and insight can be obtained on the behavior of the
solutions that are generated.
Fig.50 shows the output flow diagram that is presently being utilized
Output from the model is written on tape in addition to hard copy
and digital overplots. The latter plots are for general screening
152
-------
UJ
0
D_
UJ
Q
0
10
20
30
40
50
60
70
80
LEGEND
CHLOROPHYLLa, meg/1
4 6 8 10 12
14 16
CHLOROPHYLLa, meg/1
4 6 8 10 12 14 16
i r
Segment 15
Day 150-180 (June)
CD
0
10
20
30
40
50
60
70
80
Computed for the time shown:
• Mean
/ Range
Observed IFYGL data over the
time shown: o
Segment 17
Day 150-180 (June)
Figure 49. Preliminary comparison of Lake 3 output to observed data
at two segments for June
-------
/TAPE
OUTPUT
Ul
Water Quality
Eutrophication Model
PRINTED
OUTPUT
LISTING
Additional Statistical
Analysis
Q 0
Plotting Software
(Contour-3D plots)
CATHODE RAY
TUBE DISPLAY
PAPER PLOTS
DIGITAL OVERPLOTS
(Observed data vs.
computed values)
PROJECTION
CAMERA
Figure 50, Data output flow diagram used for verification
and display purposes
-------
by the analyst. Plotting software is then employed on the tape
output from a model run. Paper contour plots and three dimensional
plotting routines are employed to provide further perspective on
the output. Fig.51 shows a three dimensional plot of output of
phytoplankton biomass in the 0-4 meter layer from a view looking
down the lake from east to west. Plots such as these reveal the
structure of the solution in a meaningful way and also permit
scanning of the output to determine any strange or anomalous
results that may indicate a flaw in the numerical computations.
Returning to Fig. 51 the three-dimensional plots are also displayed
on a cathode ray tube, microfilmed and prepared for a motion picture
of the complete output over time. A screening of this film then
reveals the dynamic and spatial behavior of the solution in a very
unique way. At this stage, a film has been prepared of a pre-
liminary run for a full year of the surface phytoplankton chloro-
phyll.
155
-------
LAKE ONTAR10-LAKE 3 MODEL
LAYER 1 (0-4 METERS)
JUNE
NIAGARA RIVER
OSWEGO RIVER*
ST. LAWRENCE RIVER
PHYTOPLANKTON CHLOROPHYLL a Maximum, 8.2jig/l
Figure 51. Three dimensional plot of phytoplankton chlorophyll
calculated from Lake 3 model - June, 0-4 meters.
156
-------
SECTION VIII
REFERENCES
Simons, T.J. Development of Numerical Models of Lake Ontario.
Proc. 14th Conf. Great Lakes Res. IAGLR, pp.654-669,
(1971).
Simons, T.J. Development of Numerical Models of Lake Ontario,
Part 2. Proc, 15th Conf. Great Lakes Res., IAGLR, pp.655-
672, (1972).
Simons, T.J. Development of Three-Dimensional Numerical
Models of the Great Lakes. Scientific Series No. 12,
Inland Waters Dir.,CCIW, Burlington, Ontario, 26 pp.,(1973).
Report to the International Joint Commission on the Pollution
of Lake Ontario and the International Section of the
St. Lawrence River - Vol. 3 - International Lake Erie Water
Pollution Board and the International Lake Ontario - St.
Lawrence River Water Pollution Board.(1969)
Casey, Donald J. Lake Ontario Environmental Summary 1965.
U.S. Environmental Protection Agency Region V. Lake Ontario
Basin Office, Rochester, N.Y.
Bonham-Carter, Graeme, John H. Thomas and David Lackner
A Numerical Model of Steady Wind-Driven Currents in Lake
Ontario and the Rochester Enbayment Based on Shallow-Lake
Theory. Report Number 7 International Field Year For the
Great Lakes - Rochester Enbayment Project - Department of
Geological Sciences, University of Rochester, (1973)
157
-------
Rodgersf G.K. and G.K. Sato. Energy Budget Study for
Lake Ontario. Canadian Meteorological Research
Report, Atmospheric Environment Service, Jan. (1971)
Lee, A.H. Regional Characterizations of the Thermal
Properties of Lake Ontario. Proc. 15th Conf. Great
Lakes Research, IAGLR, 625-634, (1972)
158
-------
APPENDIX A
LAKE 1
INPUT AND PROGRAM. LISTING
INPUT INFORMATION
Introduction
Lake 1 is a kinetic interactive model spatially defined by three
vertical, mixed layers representing the epilimnion, hypolimnion
and benthos of Lake Ontario as shown in Fiaure 9. The purpose of
this appendix is to summari2e the input information needed to run
the model, and to present where available, the background informa-
tion from which this information was drawn. Units for the in-
formation given are consistent with those required by the computa-
tional software used, which is named Water Analysis Simulation
Program (WASP). A listing of WASP and sample output are included.
MORPHOMETRY AND HYDRODYNAMIC REGIME
The information required to describe the morphometry of Lake
Ontario, for the Lake 1 model is presented in Table A-l, with the
associated hydrodynamic regime.
The bulk of the morphometric data was generated using data
supplied by the Canada Centre for Inland Waters (CCIW), Descriptive
Limnology Section . The bathymetric data, supplied on punched
cards, were depth averages for a square grid laid on a polyconic
2
projection of Lake Ontario, with geographical areas of 4km .
These surface areas were defined by sides which were actually curves
on the earth's surface. These data were processed to yield the
159
-------
TABLE A-l
LAKE 1 PHYSICAL PARAMETERS
Inter- Intersegment
Segment segment Representation Volume Depth Areas
[106CF] [Feetl
[FT2]
Dispersion
Inflow Outflow Coefficient
[cfs] fcfs] [miVday]
1
2
3
1-2
2-3
Epilimnion
Hypolimnion
Benthos
l.OSxlO7
4.85xl07
4.80xl04
55.8
240.5
0.5
1.76X1011
0.959X1011
43,500.
188,400.
0.0
43,500.
188,500.
0.0
see Fig. A-4
0.0
O
-------
necessary information. The average depth of Lake Ontario is
j
296 feet, (90 meters), The International Joint Commission
reports the average depth of the developed thermocline is
55.8 feet (17 meters). This depth is used in the representation
of the epilimnion. The remainder of the volume is placed in the
hypolimnion. The benthos, which is currently used only to track
accumulation of depositions, was given an arbitrary depth of
0.5 feet. Figures A-l and A-2 which show the variation of
volume and surface area with depth were generated with the CCIW
supplied data.
The hydrodynamic regime is inputted into WASP, as opposed to a
modelling configuration which would incorporate a hydrodynamic
sub-model internally. This differentiation is most significant for
more spatially defined models than Lake 1. The hydrodynamics must
therefore be obtained from a data base, the literature, or an ex-
ternal hydrodynamic sub-model. The mean annual outflow for Lake
2
Ontario is 232,000 cfs. Inflow was set equal to outflow and was
split proportionally to the segment depths, as shown in Table A-l.
The development of the thermocline was simulated by temporally
varying the vertical dispersion coefficient as shown in Fig. A-3,
The establishment of the thermocline prevents exchange between the
epilimnion and hypolimnion. This phenomena was incorporated by
setting the vertical dispersion between segments 1 and 2 to zero
during the period of stratification. For the remainder of the year
exchange takes place between segment 1 and segment 2. The vertical
2 2
dispersion coefficient was set at 0.0003mi /day(90 cm /sec) for
periods of non-stratification, which is in the range reported in
o
Limnological Systems Analysis .
161
-------
30
60
90
120
150
180
210
240
PERCENT VOLUME
10 20 30 40 50 60 70 80 90 100
I I i
- 200 400 600 800 1000 1200 1400 1600 KM3
100
200
300
500
600
700
50 100 150 200 250 300 350 400 Ml3
Figure A-l. Variation of volume with depth
fil
30
60
90
:- 120
150
180
210
240
PERCENT AREA
10 20 30 40 50 60 70 80 90 100
I
I
I
I
2
i
4
i
T
8 10 12 14 16
100
200
300
400
500
600
fc
°
700
1000km'
i i
12345
Figure A-2. Variation of area with depth.
162
6 lOOOmi2
-------
PHYSICAL EXOGENEOtlS VARIABLES
1. Water Temperature - the temporal variation of water tempera-
ture used in Lake 1 is shown in Fig. A-4. Since the segments
represent completely mixed volumes, representative temperatures
were chosen. Data averaged from CCIWs Limnological Data Reports
4
(1966,1967,1968,1969) were used as temperature input . Only
stations with depths over 50 meters were considered. These were
viewed as representative of main lake conditions.
2. Solar Radiation - data recorded at Brockport, New York, and
at Toronto-Scarborough, Ontario were available. The Toronto-
Scarborough data were used as this comprised the longest period of
record. The Toronto-Scarborough data were adjusted using an
over land radiation to over lake radiation ratio proposed by
Richards and Loewen . The over land radiation data, lake to land
ratio and adjusted temporal lake radiation function used is shown
in Fig. A-5.
3. Photo Period - that fraction of the day from sunrise to sunset
is variable during the year, and is a function of latitudinal
o
position. Sunset, sunrise data from the World Almanac were used
to obtain the functions shown in Fig. A-6. The Lake Ontario photo
period function was obtained by interpolating to the equivalent of
43° 40' N latitude.
4. Background Water Clarity - the variation of the extinction
coefficient with time is a function of the phytoplankton con-
centration in the water body. Given no phytoplankton, there
163
-------
'-o
0.0
100
40
20
30 60 90 120 M150 J 180 J 210 A 240 270°300 330° 360
Figure A-3. Temporal variation of vertical dispersion coefficient
interface 1-2
20
15
S
Q_
10
5 -
Depth 1966 1967 1968 1969
SEGMENT 1
Figure A-4. Temporal variation of water temperature
164
-------
700 -
600
500
fr
•o
>N
400
o
0=300
o
200
100
OVER-LAKE SOLAR
RADIATION
OVER-LAND SOLAR
RADIATION
LAND-TO-LAKE RAT 10
EXTRAPOLATED FOR JANUARY-MAY
T
1.4
1.2
1.0 g
i
0.8 ?
a
0.6 5
.30
60 M90 A120M150J 180 J 210 A 240 $ 270° 300 ** 330° 360
Figure A-5. Temporal variation of solar radiation
165
-------
-
0.8
0.7
0.6
0.5
o
o
0.4
Q.
O
I
O-
0.3
0.2
0.1
60° N
50° N
40° N
30° N
20° N
—•— Lake Ontario
0.0, r
\
18
15
12
UJ
Q-
J30F60M90A120M150J180J210A240S270°300N330D360
1974
Figure A-6. Temporal variation of photoperiod
166
-------
would still be extinction of light as it passes through the
water column attributable to other factors. This fraction of
the extinction coefficient, which is a function of background
water clarity, was determined by plotting chlorophyll and
extinction coefficient as shown in Fig. A-7. The extinction
coefficient was taken to be:
K = 1.9
e
secchi depth
9
after Beeton.
BIOLOGICAL AND CHEMICAL EXOGENEOUS VARIABLES
1. Mass Loadings and Boundary Concentrations - sources of
nutrient input into Lake Ontario are varied. Tributary runoff,
direct municipal and industrial discharges comprise the principal
sources of nutrients. Mass loadings (mass/ unit time) are
normally thought of as direct inputs to a water body. Boundary
concentration conditions (mass/unit volume) are usually thought
to be associated with tributary inflow or in dispersive situ-
ations with the 'downstream1 portion of the water body. Boundary
conditions, since they are always associated with a flow regime
(volume/ unit time), can also be viewed as mass loadings into the
water body. WASP has capabilities to input nutrients using
either mass loadings or boundary conditions, with the associated
flow regime. Lake 1 in its present configuration just utilizes the
boundary condition form of mass input. The loading information
used for the various nitrogen and phosphorous species used in the
nutrient systems (non-living organic, ammonia and nitrate,
167
-------
1.0
F*
E 0.8
•
*^>
lo.6
U_
EXTINCTIONCOEF
o p
k» *.
n n
D 1969
A 1968
•w
• 1967
A
a D
•v
D
•
ADD 0
A*
a A D
A • A
V
_ K; = 0.22 nr1
l l l l l l 1 1 1 1
0 12
34567
CHLOROPHYLL a,
89 10
Figure A-7. Determination of background water
clarity - Ke
168
-------
nitrogen, non-living organic and available phosphorus) was not
available. Only total nutrient information was given. Specie
loads were set at initial conditions for the inorganic nutrients
with the remainder placed in the respective non-living organic
system. Loading information and breakdown is shown in Table A-2.
2. Initial Conditions - the initial conditions for Lake 1 are
shown in Table A-2. The integration starts the first of January and
these values reflect the measured values for this part of the year.
CHEMICAL AND BIOLOGICAL ENDOGENEOUS VARIABLES
1. System Parameters - the various constants and parameters
which define the interaction between the various biological and
chemical systems are listed in Table A-3. Units are those re-
quired by the present configuration of Lake 1.
169
-------
TABLE A - 2
LAKE 1
MASS LOADINGS AND INITIAL CONDITIONS
System
Phytoplankton
Zooplankton - Herbivorous
Nitrogen
Non-living Organic
Ammonia - N
Nitrate - N
Phosphorus
Non-living Organic
Available - P
Zooplankton - Carnivorus
Zooplankton - Upper Trophic
Level 1
Zooplankton - Upper Trophic
Level 2
Boundary
Concentration
(mg/1)
1.0
.05
.443
.0150
.250
.0454
.0143
.01
by-passed system
by-passed system
Initial
Conditions
(mg/1)
2.0
.005
.090
.0150
.250
.005
.0143
.005
170
-------
TABLE A - 3
SYSTEM PARAMETERS
System - Parameter
Phytoplankton
Chlorophyll a
1. Saturated Growth Rate
(ej
2. Saturating Light
Intensity
3. Michaelis Constant
for Nitrogen
4. Michaelis Constant
for Phosphorus
5. Endogenous Respiration
Rate (620)
6. Carbon to Chloro-
phyll Ratio
7. Sink Velocity
8. Preference Structure
for Inorganic
Nitrogen
Symbol
KIT
K1C
IS
KMN
KMP
K2T
K2C
CCHL
SVEL
PNH3
Value
1.066
0.580
350.
0.025
0.002
1.08
0.10
0.05
0.10
0.5
Units
None
_i
day
ly/day
mg-N/1
mg-P/1
None
day'1
mg-c/ygChla
m/ day
None
171
-------
TABLE A - 3
SYSTEM PARAMETERS
(Continued)
System - Parameter
Non-living organic
Nitrogen
1. Organic N - NH3
Hydrolysis Rate
2. Nitrogen to Chloro-
phyll Ration
a
3. Sink Coefficient
Ammonia Nitrogen
1. NH-.-NO, Nitrification
3 3 Rate
Nitrate Nitrogen ~ none
Non-living organic
Phosphorus
1. Organic P-Inorganic P
conversion rate
2. Phosphorus to
chorophyll a
3. Sink Coefficient
Available Phosphorus
1. Sink Coefficient
Symbol
a?
NCHL
K33
K45T
K45C
K67T
K67C
PCHL
K 66
K 77
1
Value
0.00175
0
0.010
0.001
0.002
0
0.007
0
0.001
0.001
0
Units
(day'C)"1
day-J-
mgN/mg Chi
a
day'1
(day°C)~1
day"1
o -1
(day C)
mgP/yg Chi
a
day "I
day
172
-------
TABLE A - 3
SYSTEM PARAMETERS
(Continued)
System - Parameter
Herbivorous
Zooplankton-carbon
1. Grazing Rate
2. Michaelis Constant
for Phytoplankton
3. Conversion Efficiency
4. Endogenous Respiration
Rate
Carnivorous
Zoop lank ton- carbon
1. Grazing Rate
2. Conversion Efficiency
3. Endogenous Respiration
Rate
Upper Trophic
Level 1 - carbon
1. Grazing Rate
2. Conversion Efficiency
3. Endogenous Respiration
Rate
Symbol
CGC
CGT
KMPL
AZP
K4T
K4C
CGT8
CGC8
A8Z
K8T
K8
CGT9
CGC9
A98
K9T
K9
Value
0
0.06
10
0.60
0.001
0.0
0.06
0
0.60
0.001
0
0.01
0
0.60
0.00
0.0
Units
o
1/mgC/day C
1/mgC/day
ygChla/1
none
(day'C)-1
l/mgC/day°C
1/mgC/day
none
(day°C)~1
day"1
l/mgC/day°C
1/mgC/day
none
(day)-l
173
-------
TABLE A - 3
SYSTEM PARAMETERS
(Continued)
System - Parameter
Upper Trophic
Level 2 - carbon
1. Grazing Rate
2. Conversion
Efficiency
3 . Endogenous
Respiration Rate
Symbol
CGTO
CGCO
A109
KLOT
K10
Value
0.06
0
0.60
0.002
0
Units
o
1/mg C/day C
o
1/mgC/day C
None
(day°C)~1
day"1
174
-------
NOTE:
Source listings for Lake 1 model and sample input and output
are available from EPA, Grosse lie Laboratory, Grosse lie,
Michigan.
This program is currently operational on the Control Data Corp
6600 at Courant Institute of New York University.
A version will be available on Optimal Systems Inc., in the
near future.
175
-------
REFERENCES
1. Canada Centre for Inland Waters, Descriptive Limnology
Section, personal communication.
2. International Lake Erie Water Pollution Board and the
International Lake Ontario - St. Lawrence River Water
Pollution Board. Report to the International Joint
Commission on the Pollution of Lake Ontario and the
International Section of the St. Lawrence River,
Volume 3, 1969.
3. Hydroscience Inc. Limnological Systems Analysis of the
Great Lakes-Phase I. Hydroscience Inc., Westwood, N.J.
1973.
4. Canada Centre for Inland Waters. Limnological Data Reports.
Lake Ontario. Canadian Oceanographic Data Centre,
Burlington, Ontario.
5. Hubbard, John E. Department of Geology & Earth Science,
State University College, Brockport, N.Y., personal
communication.
6. Phillips, D.W. & J.A.W. McCulloch. The Climate of the
Great Lakes Basin. Climatological Study No. 20,
Environment Canada, Toronto. 1972.
7. Richards, T.L. and P. Lowen. A preliminary investigation
of solar radiation over the Great Lakes as compared to
adjacent land areas. Proc. Eight Conf. Great Lakes Res.,
Great Lakes Res. Division, Univ. of Michigan, pp. 278-282.
1965.
176
-------
REFERENCES
(Continued)
8. The World Almanac and Book of Facts. Newspaper Enterprise
Association, New York. 1974.
9. Beeton, A.M. Relationship between Secchi disc readings and
light penetration in Lake Huron, Am. Fish. Soc. Trans.
87:73-79. 1958.
177
-------
TECHNICAL REPORT DATA
I'li'asc read Inunctions on tin- reverse before completing)
m com NO. T;>.
JSPA-660/3-75-005 j _ ___
]"flftt ANUSUUIULE
Mathematical Modeling of Phytoplankton in Lake Ontario
1. Model Development and Verification
3. RECIPIENT'S ACCESSION-NO.
5. REPORT DATE
October 1974
6. PERFORMING ORGANIZATION CODE
AUTHOH(S)
Robert V. Thomann
Dominick M. DiToro
Richard P. Winfield
Donald J. O'Connor
8. PERFORMING ORGANIZATION REPORT NO.
EPA-660/3-75-005
9. PI FORMING ORG \NIZATION NAME AND ADDRESS
Manhattan College
Environmental Engineering and Science Program
Bronx, N.Y. 10471
10. PROGRAM ELEMENT NO.
1BA026
11. CONTRACT/GRANT NO.
R800610
SPONSORING Am NCY NAM£ AND ADDRESS ,
"U.S. EnvxronmentaT Trotectxon Agency
National Environmental Research Center
Grosse He Laboratory
9311 Groh Rd. , Grosse He, Michigan 48138
13. TYPE OF REPORT AND PERIOD COVERED
Final, Part 1
14. SPONSORING AGENCY CODE
IB.-fUPPl EMENTARY NOTES
The basic mathematical structure for describing the dynamics of phytoplankton in
Lake Ontario is presented in this report. Data on chlorophyll and principle nu-
trients are reviewed and summarized and the mathematical modeling strategy is de-
tailed. The modeling strategy begins with the construction of a horizontally
completely mixed lake with vertical layers, LAKE 1. This spatially simplified
model is used to develop the interactions and kinetic behavior of the various
components of each subsystem. A more detailed 3-dimensional model is then used to
describe open lake and near shore variations in phytoplankton biomass. Ten
biological and chemical variables are used in both models and include four trophic
levels above the phytoplankton, chlorophyll a as a measure of phytoplankton biomass
two phosphorus components and three nitrogen components. Under reasonable sets of
model parameters as reported in the literature, the Lake 1 model output compared
favorably with observed data on the dependent variables. Spring growth and peak
chlorophyll concentrations are related primarily to increasing light and temperatur
The model indicates that growth ceases due to phosphorus limitation. The results
to date indicate that the mathematical model of phytoplankton in Lake Ontario as
developed herein is a reasonable first approximation to observed data. As such,
the model can form a basis for preliminary estimates of the effects of nutrient re-
duction programs on Lake Ontario.
t7.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
Mathematical Models, Computer Models,
Eutrophication, Water Quality, Nutrients
Cycling Nutrients, Dispersion, Mass
Transfer, Simulation Analysis
b.lDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Group
Great Lakes, Lake
Ontario
11.
STATEMENT
19. SECURITY CLASS (This Report)
21. NO. OF PAGES
20. SECURITY CLASS (Thltpoge)
22. PRICE
•PA Form 2120-1 (»-73)
•S U.S GOVERNMENT FS:NTING OFFICE: 1975-698- 162/ni REGION 10
------- |